\[ \]
\[ \]
\[\begin{align*} y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \boldsymbol\mu_0 + \epsilon_t\\ \epsilon_t|Y_{t-1} &\sim iidN\left(\mathbf{0}_N,\mathbf\Sigma\right) \end{align*}\] for \(t=1,\dots,T\)
Let the number of variable \(N=2\) and the lag order \(p=2\).
Then, the model equation is:
\[\begin{align*} \begin{bmatrix} y_{1.t} \\ y_{2.t} \end{bmatrix} &= \begin{bmatrix} \mathbf{A}_{1.11} & \mathbf{A}_{1.12} \\ \mathbf{A}_{1.21} & \mathbf{A}_{1.22} \end{bmatrix} \begin{bmatrix} y_{1.t-1} \\ y_{2.t-1} \end{bmatrix} + \begin{bmatrix} \mathbf{A}_{2.11} & \mathbf{A}_{2.12} \\ \mathbf{A}_{2.21} & \mathbf{A}_{2.22} \end{bmatrix} \begin{bmatrix} y_{1.t-2} \\ y_{2.t-2} \end{bmatrix} + \begin{bmatrix} \boldsymbol\mu_{0.1} \\ \boldsymbol\mu_{0.2} \end{bmatrix} + \begin{bmatrix} \epsilon_{1.t} \\ \epsilon_{2.t} \end{bmatrix}\\[2ex] \begin{bmatrix} \epsilon_{1.t} \\ \epsilon_{2.t} \end{bmatrix} &\Big|Y_{t-1} \sim iid N_2\left( \begin{bmatrix} 0\\ 0\end{bmatrix}, \begin{bmatrix}\boldsymbol\sigma_1^2 & \boldsymbol\sigma_{12} \\ \boldsymbol\sigma_{12} & \boldsymbol\sigma_2^2\end{bmatrix} \right) \end{align*}\]A \(K\times N\) matrix \(\mathbf{A}\) is said to follow a matrix-variate normal distribution: \[ \mathbf{A} \sim MN_{K\times N}\left( M, Q, P \right), \] where
if \(\text{vec}(\mathbf{A})\) is multivariate normal: \[ \text{vec}(\mathbf{A}) \sim N_{KN}\left( \text{vec}(M), Q\otimes P \right) \]
An \(N\times N\) square symmetric and positive definite matrix \(\mathbf\Sigma\) follows an inverse Wishart distribution: \[ \mathbf\Sigma \sim IW_{N}\left( S, \nu \right) \] where
then the joint distribution of \((\mathbf{A},\mathbf\Sigma)\) is normal-inverse Wishart \[ p(\mathbf{A},\mathbf\Sigma) = NIW_{K\times N}\left( M,P,S,\nu\right) \]
The model assumptions state: \[\begin{align*} \epsilon_t|Y_{t-1} &\sim iidN_N\left(\mathbf{0}_N,\mathbf\Sigma\right) \end{align*}\]
Collect error term vectors in a \(T\times N\) matrix: \[\underset{(T\times N)}{E}= \begin{bmatrix}\epsilon_1 & \epsilon_2 & \dots & \epsilon_{T}\end{bmatrix}'\]
Error term matrix is matrix-variate distributed: \[\begin{align*} E|X &\sim MN_{T\times N}\left(\mathbf{0}_{T\times N},\mathbf\Sigma, I_T\right) \end{align*}\]
The inverse Wishart density function is proportional to: \[\begin{align*} \text{det}(\mathbf\Sigma)^{-\frac{\nu+N+1}{2}}\exp\left\{ -\frac{1}{2}\text{tr}\left[ \mathbf\Sigma^{-1} S \right] \right\} \end{align*}\]
Consider a case where:
\[ \underset{(K\times N)}{\mathbf{A}}=\begin{bmatrix} \mathbf{A}_1'\\ \vdots \\ \mathbf{A}_p' \\ \boldsymbol\mu_0' \end{bmatrix} \quad \underset{(T\times N)}{Y}= \begin{bmatrix}y_1' \\ y_2'\\ \vdots \\ y_T'\end{bmatrix} \quad \underset{(K\times1)}{x_t}=\begin{bmatrix} y_{t-1}\\ \vdots \\ y_{t-p}\\ 1 \end{bmatrix}\quad \underset{(T\times K)}{X}= \begin{bmatrix}x_1' \\ x_2' \\ \vdots \\ x_{T}'\end{bmatrix} \quad \underset{(T\times N)}{E}= \begin{bmatrix}\epsilon_1' \\ \epsilon_2' \\ \vdots \\ \epsilon_{T}'\end{bmatrix} \] where \(K=pN+1\)
Define the MLE: \(\widehat{A}=(X'X)^{-1}X'Y\)
Perform simple transformation of the likelihood
\[\begin{align*} L\left(\mathbf{A},\mathbf{\Sigma}|Y,X\right)&\propto\text{det}(\mathbf{\Sigma})^{-\frac{T}{2}}\exp\left\{-\frac{1}{2}\text{tr}\left[\mathbf{\Sigma}^{-1}(Y-X\mathbf{A})'(Y-X\mathbf{A})\right]\right\}\\ &=\text{det}(\mathbf{\Sigma})^{-\frac{T}{2}}\\ &\quad\times\exp\left\{ -\frac{1}{2}\text{tr}\left[\mathbf{\Sigma}^{-1}(\mathbf{A}-\widehat{A})'X'X(\mathbf{A}-\widehat{A}) \right] \right\}\\ &\quad\times \exp\left\{ -\frac{1}{2}\text{tr}\left[ \mathbf{\Sigma}^{-1}(Y-X\widehat{A})'(Y-X\widehat{A}) \right] \right\} \end{align*}\]Under the likelihood, \((\mathbf{A},\mathbf{\Sigma})\) are normal-inverse Wishart distributed:
\[\begin{align*} L\left( \mathbf{A},\mathbf{\Sigma}|Y,X \right) &= NIW_{K\times N}\left(\widehat{A},(X'X)^{-1},(Y-X\widehat{A})'(Y-X\widehat{A}), T-N-K-1 \right) \end{align*}\]A natural-conjugate prior leads to joint posterior distribution for \((\mathbf{A},\mathbf{\Sigma})\) of the same form \[\begin{align*} p\left( \mathbf{A}, \mathbf{\Sigma} \right) &= p\left( \mathbf{A}| \mathbf{\Sigma} \right)p\left( \mathbf{\Sigma} \right)\\ \mathbf{A}|\mathbf{\Sigma} &\sim MN_{K\times N}\left( \underline{A},\mathbf{\Sigma},\underline{V} \right)\\ \mathbf{\Sigma} &\sim IW_N\left( \underline{S}, \underline{\nu} \right) \end{align*}\]
\[ \]
Hint: Proceed by:
Posterior mean of matrix \(\mathbf{A}\) is: \[\begin{align*} \overline{A} &= \overline{V}\left( X'Y + \underline{V}^{-1}\underline{A} \right)\\[2ex] &= \overline{V}\left( X'X\widehat{A} + \underline{V}^{-1}\underline{A} \right)\\[2ex] &= \overline{V} X'X\widehat{A} + \overline{V}\underline{V}^{-1}\underline{A} \end{align*}\] a linear combination of the MLE \(\widehat{A}\) and the prior mean \(\underline{A}\)
Note that: \[ \overline{V} X'X + \overline{V}\underline{V}^{-1} = \overline{V} ( X'X + \underline{V}^{-1}) = I_K \]
According to Bayes Rule, the kernel of the posterior is normalised by the Marginal Data Density \(p(data)\):
\[ p\left( \mathbf{A}, \mathbf{\Sigma}| data \right) = \frac{L(\mathbf{A},\mathbf{\Sigma}| data)p\left( \mathbf{A}, \mathbf{\Sigma} \right)}{p(data)} \]
For Bayesian VARs the posterior is known \[ p\left( \mathbf{A}, \mathbf{\Sigma}| data \right) = MNIW\left(\overline{A},\overline{V}, \overline{S}, \overline{\nu} \right) \]
and so is the analytical formula for the MDD: \[p(data)\]
This can be used to our advantage!
Sims, Litterman, Doan (1984) proposed an interpretable way of setting the hyper-parameters on the NIW prior \(\underline{A}\), \(\underline{V}\), \(\underline{S}\), and \(\underline{\nu}\) for macroeconomic data.
\[ \] The prior reflects the following stylised facts about macro time series:
Set
\[\begin{align*} \underline{S} &= \begin{bmatrix} \psi_1 &0 &\dots & 0 \\ 0 & \psi_2 &\dots & 0\\ \vdots &\vdots&\ddots& \vdots\\ 0&0&\dots&\psi_N \end{bmatrix}\\[2ex] \underline{\nu}&= N+2 \end{align*}\]\(\psi =(\psi_1, \dots, \psi_N)\) have to be chosen (or estimated)
For \(\quad l = 1+\text{floor}((i-1)/N) \quad\text{and }\quad k = i - (l-1)N\), set: \[\begin{align*} \underline{A} &= \begin{bmatrix} I_N \\ \mathbf{0}_{((p-1)N +1)\times N}\end{bmatrix}& \underline{V}_{ij} &= \left\{\begin{array} (\lambda ^ 2 / (\psi_k l^2) &\text{ for }i=j,\text{ and } i\neq pN+1 \\ \lambda^2 &\text{ for }i=j,\text{ and } i= pN+1 \\ 0&\text{ for } i\neq j \end{array}\right. \end{align*}\]
\(\lambda^2\) has to be chosen (or estimated)
Consider a simple case of a model for \(N=2\) observations and with \(p=1\) lag.
Use Bayes Rule to derive the joint prior of \((\mathbf{A},\mathbf\Sigma)\) given \(Y^*\) and \(X^*\).
It is given by the MNIW distribution:
\[\begin{align*} p\left( \mathbf{A}, \mathbf{\Sigma}|Y^*,X^* \right) &= MNIV_{K\times N}\left( \underline{A}^*,\underline{V}^*, \underline{S}^*, \underline{\nu}^* \right) \end{align*}\] \[\begin{align*} \underline{V}^*&= \left( X^{*\prime}X^* + \underline{V}^{-1}\right)^{-1} & \underline{A}^*&= \underline{V}^*\left( X^{*\prime}Y^* + \underline{V}^{-1}\underline{A} \right)\\ \underline{\nu}^*&= T_d+\underline{\nu} & \underline{S}^*&= \underline{S}+Y^{*\prime}Y^* + \underline{A}'\underline{V}^{-1}\underline{A} - \underline{A}^{*\prime}\underline{V}^{*-1}\underline{A}^* \end{align*}\]
Let a \(p\times N\) matrix \(Y_0\) denote the initial observations, that is, the first \(p\) observations of the available time series.
Let an \(N\)-vector \(\bar{Y}_0\) denote its columns’ means.
Generate additional \(N\) rows by \[ Y^+ = \text{diag}\left(\frac{\bar{Y}_0}{\mu}\right) \quad\text{ and }\quad X^+ = \begin{bmatrix}Y^+ & \dots & Y^+ & \mathbf{0}_N \end{bmatrix} \]
Generate an additional row by \[ Y^{++} = \frac{\bar{Y}_0'}{\delta} \quad\text{ and }\quad X^{++} = \begin{bmatrix} Y^{++} & \dots & Y^{++} & \frac{1}{\delta} \end{bmatrix} \]
\[ Y^* = \begin{bmatrix}Y^+ \\ Y^{++} \end{bmatrix}\quad\text{ and }\quad X^* = \begin{bmatrix}X^+ \\ X^{++} \end{bmatrix} \]
Suppose that:
Write out the matrices \(Y^*\) and \(X^*\) of dimensions \(2\times 3\) and \(3\times 3\) respectively.
Hyper-parameters \(\psi\), \(\lambda\), \(\mu\) and \(\delta\) can be fixed to values chosen by the econometrician.
A better idea is to assume priors for these hyper-parameters and estimate them as in Giannone, Lenza, Primiceri (2015).
Extend the existing prior to: \[\begin{align*} p\left( \mathbf{A}, \mathbf{\Sigma}|Y^*,X^*,\psi,\lambda,\mu,\delta \right) &= MNIV_{K\times N}\left( \underline{A}^*,\underline{V}^*, \underline{S}^*, \underline{\nu}^* \right) \end{align*}\] And specify: \[\begin{align*} \psi_n &\sim IG\left(0.02^2, 0.02^2\right)\\ \lambda &\sim G\left(0.2,2\right)\\ \mu &\sim G\left(1,2\right)\\ \delta &\sim G\left(1,2\right) \end{align*}\]
Giannone, Lenza, Primiceri (2015) propose the following estimation procedure:
\[ p(\psi,\lambda,\mu,\delta|data) \propto p(\psi,\lambda,\mu,\delta)p(data|\psi,\lambda,\mu,\delta)\] - Use an \((N+3)\)-variate Student-t distribution as the candidate generating density
Use the MNIW posterior derived for the implied prior:
\[\begin{align*} p\left( \mathbf{A}, \mathbf{\Sigma}|Y,X, Y^*,X^* \right) &= MNIW_{K\times N}\left( \overline{A}^*,\overline{V}^*,\overline{S}^*, \overline{\nu}^* \right)\\[2ex] \overline{V}^*&= \left( X'X + \underline{V}^{*-1}\right)^{-1} \\ \overline{A}^*&= \overline{V}^*\left( X'Y + \underline{V}^{*-1}\underline{A}^* \right)\\ \overline{\nu}^*&= T+\underline{\nu}^*\\ \overline{S}^*&= \underline{S}^*+Y'Y + \underline{A}^{*\prime}\underline{V}^{*-1}\underline{A}^* - \overline{A}^{*\prime}\overline{V}^{*-1}\overline{A}^* \end{align*}\]Package bsvarSIGNs by Wang & Woźniak (2024) implements this algorithm.
\(\left.\right.\)
… is to use the available data to provide a statistical characterisation of the unknown future values of quantities of interest.
\(\left.\right.\)
The full statistical characterisation of the unknown future values of random variables is given by their predictive density.
\(\left.\right.\)
Simplified outcomes in a form of statistics summarising the predictive densities are usually used in decision-making processes.
\(\left.\right.\)
Summary statistics are also communicated to general audiences.
\[\begin{align*} y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \boldsymbol\mu_0 + \epsilon_t\\[2ex] \epsilon_t|Y_{t-1} &\sim iidN_N\left(\mathbf{0}_N,\mathbf\Sigma\right)\\ & \end{align*}\]
… is implied by the model formulation: \[\begin{align*} y_{t+h}|Y_{t+h-1},\mathbf{A},\mathbf\Sigma &\sim N_N\left(\mathbf{A}_1 y_{t+h-1} + \dots + \mathbf{A}_p y_{t+h-p} + \boldsymbol\mu_0,\mathbf\Sigma\right) \end{align*}\]
\(\left.\right.\)
Bayesian forecasting takes into account the uncertainty w.r.t. parameter estimation by integrating it out from the predictive density.
\[\begin{align*} &\\ p(y_{T+1}|Y,X) &= \int p(y_{T+1}|Y_{T},\mathbf{A},\mathbf\Sigma) p(\mathbf{A},\mathbf\Sigma|Y,X) d(\mathbf{A},\mathbf\Sigma)\\ & \end{align*}\]
\(\left.\right.\)
… and obtain \(S\) draws \(\left\{ \mathbf{A}^{(s)},\mathbf\Sigma^{(s)} \right\}_{s=1}^{S}\)
\(\left.\right.\)
In order to obtain draws from \(p(y_{T+1}|Y,X)\), for each of the \(S\) draws of \((\mathbf{A},\mathbf\Sigma)\) sample the corresponding draw of \(y_{T+1}\):
Sample \(y_{T+1}^{(s)}\) from \[ N_N\left(\mathbf{A}_1^{(s)} y_{T} + \dots + \mathbf{A}_p^{(s)} y_{T-p+1} + \boldsymbol\mu_0^{(s)},\mathbf\Sigma^{(s)}\right) \] and obtain \(\left\{y_{T+1}^{(s)}\right\}_{s=1}^{S}\)
\(\left.\right.\)
This procedure can be generalised to any forecasting horizon.
This is an illustration for \(h=2\).
\[\begin{align*} &\\ p(y_{T+2},y_{T+1}|Y,X) &= \int p(y_{T+2},y_{T+1}|Y_{T},\mathbf{A},\mathbf\Sigma) p(\mathbf{A},\mathbf\Sigma|Y,X) d(\mathbf{A},\mathbf\Sigma)\\[1ex] &= \int p(y_{T+2}|y_{T+1},Y_{T},\mathbf{A},\mathbf\Sigma)p(y_{T+1}|Y_{T},\mathbf{A},\mathbf\Sigma) p(\mathbf{A},\mathbf\Sigma|Y,X) d(\mathbf{A},\mathbf\Sigma)\\ & \end{align*}\]
\(\left.\right.\)
… and obtain \(S\) draws \(\left\{ \mathbf{A}^{(s)},\mathbf\Sigma^{(s)} \right\}_{s=1}^{S}\)
For each of the \(S\) draws, sample \(y_{T+1}^{(s)}\) from \[ N_N\left(\mathbf{A}_1^{(s)} y_{T} + \dots + \mathbf{A}_p^{(s)} y_{T-p+1} + \boldsymbol\mu_0^{(s)},\mathbf\Sigma^{(s)}\right) \]
For each of the \(S\) draws, sample \(y_{T+2}^{(s)}\) from \[ N_N\left(\mathbf{A}_1^{(s)} y_{T+1}^{(s)} + \mathbf{A}_2 y_{T} + \dots + \mathbf{A}_p^{(s)} y_{T-p+2} + \boldsymbol\mu_0^{(s)},\mathbf\Sigma^{(s)}\right) \]
and obtain \(\left\{y_{T+2}^{(s)},y_{T+1}^{(s)}\right\}_{s=1}^{S}\)
# set seed
set.seed(123)
# MCMC settings
burn = 1000
S = 10000
# specify the model
spec = specify_bsvarSIGN$new(
optimism[,c(1,2,5)],
p = 4,
)
# estimate all hyper-parameters
spec$prior$estimate_hyper(
S = S + burn,
burn_in = burn,
mu = TRUE,
delta = TRUE,
lambda = TRUE,
psi = TRUE
)**************************************************|
Adaptive Metropolis MCMC: hyper parameters |
**************************************************|
forecasts_prod = fore$forecasts[1,,]
h = dim(forecasts_prod)[1]
limits.1 = range(forecasts_prod)
point.f = apply(forecasts_prod,1,mean)
interval.f = apply(forecasts_prod,1,HDInterval::hdi,credMass=0.90)
x = seq(from=limits.1[1], to=limits.1[2], length.out=100)
z = matrix(NA,h,99)
for (i in 1:h){
z[i,] = hist(forecasts_prod[i,], breaks=x, plot=FALSE)$density
}
x = hist(forecasts_prod[i,], breaks=x, plot=FALSE)$mids
yy = 1:h
z = t(z)
library(plot3D)
theta = 150
phi = 15.5
f4 = persp3D(x=x, y=yy, z=z, phi=phi, theta=theta, xlab="\nproductivity[t+h|t]", ylab="h", zlab="\npredictive densities of productivity", shade=NA, border=NA, ticktype="detailed", nticks=3,cex.lab=1, col=NA,plot=FALSE)
perspbox (x=x, y=yy, z=z, bty="f", col.axis="black", phi=phi, theta=theta, xlab="\nproductivity[t+h|t]", ylab="h", zlab="\npredictive densities of productivity", ticktype="detailed", nticks=3,cex.lab=1, col = NULL, plot = TRUE)
polygon3D(x=c(interval.f[1,],interval.f[2,h:1]), y=c(1:h,h:1), z=rep(0,2*h), col = "#F500BD", NAcol = "white", border = NA, add = TRUE, plot = TRUE)
for (i in 1:h){
f4.l = trans3d(x=x, y=yy[i], z=z[,i], pmat=f4)
lines(f4.l, lwd=0.9, col="black")
}
f4.l1 = trans3d(x=point.f, y=yy, z=0, pmat=f4)
lines(f4.l1, lwd=2, col="black")
mean_A = apply(post$posterior$A, 1:2, mean)
rownames(mean_A) = colnames(optimism[,c(1,2,5)])
colnames(mean_A) = c(paste0("A",1:4 %x% rep(1,3)),"mu0")
knitr::kable(mean_A, caption = "Posterior estimates for autoregressive parameters", digits = 2)| A1 | A1 | A1 | A2 | A2 | A2 | A3 | A3 | A3 | A4 | A4 | A4 | mu0 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| productivity | 0.96 | 0.00 | -0.09 | 0.09 | 0.00 | 0.02 | -0.08 | -0.01 | 0.11 | 0.03 | 0.01 | -0.04 | 0.01 |
| stock_prices | 0.09 | 1.10 | -0.48 | 0.02 | -0.11 | 0.48 | 0.18 | 0.02 | -0.36 | -0.28 | 0.00 | 0.36 | -0.02 |
| hours_worked | -0.06 | 0.02 | 1.35 | 0.01 | 0.01 | -0.25 | 0.08 | -0.02 | -0.11 | -0.03 | -0.01 | 0.01 | 0.00 |
mean_S = t(apply(post$posterior$Sigma, 1:2, mean))
mean_S = cbind(mean_S, cov2cor(mean_S))
rownames(mean_S) = colnames(optimism[,c(1,2,5)])
colnames(mean_S) = c(rep("cov",3),rep("cor",3))
knitr::kable(mean_S, caption = "Posterior estimates for covariance", digits = 5)| cov | cov | cov | cor | cor | cor | |
|---|---|---|---|---|---|---|
| productivity | 7e-05 | -0.00005 | 0e+00 | 1.00000 | -0.07692 | 0.01586 |
| stock_prices | -5e-05 | 0.00662 | 7e-05 | -0.07692 | 1.00000 | 0.12352 |
| hours_worked | 0e+00 | 0.00007 | 4e-05 | 0.01586 | 0.12352 | 1.00000 |
mean_h = rbind(apply(run$hyper, 2, mean), apply(run$hyper, 2, sd))
rownames(mean_h) = c("E[hyper|data]", "sd[hyper|data]")
knitr::kable(mean_h, caption = "Posterior estimates for hyper-parameters", digits = 3)| lambda | soc | sur | |
|---|---|---|---|
| E[hyper|data] | 1.955 | 0.314 | 0.947 |
| sd[hyper|data] | 0.302 | 0.180 | 0.495 |
\[ \]