Homoscedastic steady-state BVAR (Villani, 2009)
Source:vignettes/Homoscedastic-steady-state-BVAR.Rmd
Homoscedastic-steady-state-BVAR.RmdHere we estimate the original (homoscedastic - i.e. constant innovation covariance matrix \(\Sigma_{u}\)) steady-state BVAR model from Section 4.1 of Villani (2009).
First, let us attach the package and load the data.
library(SteadyStateBVAR)
data("Villani2009")
yt <- Villani2009The data set contains quarterly data for Sweden over the time period 1980Q1–2005Q4. The seven variables are: trade-weighted measures of foreign GDP growth \((\Delta y_f)\), CPI inflation \((\pi_f)\) and the 3-month interest rate \((i_f)\), the corresponding domestic variables (\(\Delta y\), \(\pi\) and \(i\)), and the level of the real exchange rate defined as \(q=s+p_f-p\), where \(p_f\) and \(p\) are the foreign and domestic CPI levels (in logs) and \(s\) is the (log of the) trade-weighted nominal exchange rate. As such, we have
\[ y_t= \begin{pmatrix} \Delta y_f \\ \pi_f \\ i_f \\ \Delta y \\ \pi \\ i \\ q \end{pmatrix} \]
Also, we will leave out the last two observations, so the user can compare the forecasts produced here to the last forecasts seen in Figures 1-3 in Villani (2009) to verify that this implementation works correctly.
Also, let us create the bvar object which we will use throughout here.
bvar_obj <- bvar(data = yt)To model the Swedish financial crisis at the beginning of the 90s and the subsequent shift in monetary policy to inflation targeting and flexible exchange rate, \(d_t\) (deterministic variables at time \(t\)) includes a constant term and a dummy for the pre-crisis period, i.e.
\[ d_{t}' = \begin{cases} \begin{pmatrix}1 & 1\end{pmatrix} & \text{if } t \le 1992Q4 \\ \begin{pmatrix}1 & 0\end{pmatrix} & \text{if } t > 1992Q4 \end{cases} \]
To formulate a prior on \(\Psi\), note that the specification of \(d_t\) implies the following parametrization of the steady state:
\[ \mu_t = \begin{cases} \psi_1 + \psi_2 & \text{if } t \le 1992Q4 \\ \psi_1 & \text{if } t > 1992Q4 \end{cases} \]
where \(\psi_i\) denotes the \(i\):th column of \(\Psi\). We are now ready to set up the model. Although it is not mentioned which lag length is used in Villani (2009), we assume \(p=4\).
bvar_obj <- setup(bvar_obj,
p=4,
deterministic = "constant_and_dummy",
dummy = dum_var)Now let us specify the priors. We first consider \(\beta\) (Minnesota prior). We choose the same values for the hyperparameters as in Villani (2009), i.e. an overall tightness of \(\lambda_1=0.2\), a cross-equation tightness of \(\lambda_2=0.5\), and a lag decay rate of \(\lambda_3=1\). We then specify the prior means for the first own lags of the variables. We follow Villani (2009), and as such for variables in growth rates, we set the prior mean to \(0\). For variables in levels, we set the prior mean to \(0.9\).
lambda_1 <- 0.2
lambda_2 <- 0.5
lambda_3 <- 1.0
#fol_pm = first own lag prior means
fol_pm=c(0, #delta y_f
0, #pi_f
0.9, #i_f
0, #delta y
0, #pi
0.9, #i
0.9 #q
)Now moving on to \(\Psi\), i.r. the steady-state priors, we set them according to the 95% prior probability intervals (normal distribution) in Table I in Villani (2009). We first note that for our data here, the growth rate variables (\(\Delta y_f, \pi_f, \Delta y, \pi\)) are specified in terms of quarterly rates of change/quarter-on-quarter growth, i.e. for a variable \(z\) which is on a quarterly frequency, the quarterly growth rate is \(100[ \ln (z_t) - \ln (z_{t-1})]\). The 95% prior probability intervals in Table I are specified in terms of annualized quarterly growth rates \(400 [\ln (z_t) - \ln (z_{t-1})]\).
The ppi() function is useful here. Simply input the
desired 95% prior probability interval (normal distribution) on the
annualized scale with annualized_growthrate=TRUE, and the
function returns the corresponding prior mean and variance on the
original scale (quarter-on-quarter growth). Of course, we could also
just annualize our data beforehand, and set
annualized_growthrate=FALSE.
#psi_1 = Psi col 1
#psi_2 = Psi col 2
theta_Psi <-
c(
ppi( 2.00, 3.00, annualized_growthrate=TRUE)$mean, #psi_1: delta y_f
ppi( 1.50, 2.50, annualized_growthrate=TRUE)$mean, #psi_1: pi_f
ppi( 4.50, 5.50, annualized_growthrate=FALSE)$mean, #psi_1: i_f
ppi( 2.00, 2.50, annualized_growthrate=TRUE)$mean, #psi_1: delta y
ppi( 1.70, 2.30, annualized_growthrate=TRUE)$mean, #psi_1: pi
ppi( 4.00, 4.50, annualized_growthrate=FALSE)$mean, #psi_1: i
ppi( 3.85, 4.00, annualized_growthrate=FALSE)$mean, #psi_1: q
ppi(-1.00, 1.00, annualized_growthrate=TRUE)$mean, #psi_2: delta y_f
ppi( 1.50, 2.50, annualized_growthrate=TRUE)$mean, #psi_2: pi_f
ppi( 1.50, 2.50, annualized_growthrate=FALSE)$mean, #psi_2: i_f
ppi(-1.00, 1.00, annualized_growthrate=TRUE)$mean, #psi_2: delta y
ppi( 4.30, 5.70, annualized_growthrate=TRUE)$mean, #psi_2: pi
ppi( 3.00, 5.50, annualized_growthrate=FALSE)$mean, #psi_2: i
ppi(-0.50, 0.50, annualized_growthrate=FALSE)$mean #psi_2: q
)
Omega_Psi <-
diag(
c(
ppi( 2.00, 3.00, annualized_growthrate=TRUE)$var, #psi_1: delta y_f
ppi( 1.50, 2.50, annualized_growthrate=TRUE)$var, #psi_1: pi_f
ppi( 4.50, 5.50, annualized_growthrate=FALSE)$var, #psi_1: i_f
ppi( 2.00, 2.50, annualized_growthrate=TRUE)$var, #psi_1: delta y
ppi( 1.70, 2.30, annualized_growthrate=TRUE)$var, #psi_1: pi
ppi( 4.00, 4.50, annualized_growthrate=FALSE)$var, #psi_1: i
ppi( 3.85, 4.00, annualized_growthrate=FALSE)$var, #psi_1: q
ppi(-1.00, 1.00, annualized_growthrate=TRUE)$var, #psi_2: delta y_f
ppi( 1.50, 2.50, annualized_growthrate=TRUE)$var, #psi_2: pi_f
ppi( 1.50, 2.50, annualized_growthrate=FALSE)$var, #psi_2: i_f
ppi(-1.00, 1.00, annualized_growthrate=TRUE)$var, #psi_2: delta y
ppi( 4.30, 5.70, annualized_growthrate=TRUE)$var, #psi_2: pi
ppi( 3.00, 5.50, annualized_growthrate=FALSE)$var, #psi_2: i
ppi(-0.50, 0.50, annualized_growthrate=FALSE)$var #psi_2: q
)
)Finally for \(\Sigma_u\) we will use
the noninformative Jeffreys prior \(\left|\Sigma_u \right|^{-(k+1)/2}\), as
done in Villani (2009). Now we simply pass everything to the
priors() function.
bvar_obj <- priors(bvar_obj,
lambda_1,
lambda_2,
lambda_3,
fol_pm,
theta_Psi,
Omega_Psi,
Jeffreys=TRUE)As in Villani (2009), we incorporate the assumption that Sweden is a small economy and therefore unlikely to affect the foreign economy by restricting the upper-right submatrix of \(\Pi_\ell\) for \(\ell =1,\dots,p\) or equivalently restricting the bottom-left submatrix of \(\Pi_\ell'\) to the zero matrix. This technique is called “block exogeneity” (Dieppe, Legrand, and van Roye, 2016). In essence we treat the foreign economy as exogenous to the domestic economy, although it is not exogenous in the strict sense (Karlsson, 2013).
p <- bvar_obj$setup$p
k <- bvar_obj$setup$k
kf <- 3 #first 3 variables are foreign in yt
restriction_matrix <- matrix(1, k*p, k)
for(i in 1:p){
rows <- ((i-1)*k + kf + 1) : (i*k)
cols <- 1:kf
restriction_matrix[rows, cols] <- 0
}
restriction_matrix
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 1 1 1 1 1 1 1
#> [2,] 1 1 1 1 1 1 1
#> [3,] 1 1 1 1 1 1 1
#> [4,] 0 0 0 1 1 1 1
#> [5,] 0 0 0 1 1 1 1
#> [6,] 0 0 0 1 1 1 1
#> [7,] 0 0 0 1 1 1 1
#> [8,] 1 1 1 1 1 1 1
#> [9,] 1 1 1 1 1 1 1
#> [10,] 1 1 1 1 1 1 1
#> [11,] 0 0 0 1 1 1 1
#> [12,] 0 0 0 1 1 1 1
#> [13,] 0 0 0 1 1 1 1
#> [14,] 0 0 0 1 1 1 1
#> [15,] 1 1 1 1 1 1 1
#> [16,] 1 1 1 1 1 1 1
#> [17,] 1 1 1 1 1 1 1
#> [18,] 0 0 0 1 1 1 1
#> [19,] 0 0 0 1 1 1 1
#> [20,] 0 0 0 1 1 1 1
#> [21,] 0 0 0 1 1 1 1
#> [22,] 1 1 1 1 1 1 1
#> [23,] 1 1 1 1 1 1 1
#> [24,] 1 1 1 1 1 1 1
#> [25,] 0 0 0 1 1 1 1
#> [26,] 0 0 0 1 1 1 1
#> [27,] 0 0 0 1 1 1 1
#> [28,] 0 0 0 1 1 1 1We can look at the restriction matrix for \(\beta\) to see which elements we restrict
to zero. Now we simply pass our \((kp \times
k)\) restriction matrix to the restrict_beta()
function:
bvar_obj <- restrict_beta(bvar_obj, restriction_matrix)Now, we need to supply our forecast horizon \(H\), and also a matrix containing the deterministic variables (\(d_t\)) for the future periods
\[ d_{\text{pred}}=\begin{bmatrix}d_{T+1}' \\ \vdots\\ d_{T+H}' \end{bmatrix} \]
Since the deterministic variables are i) a constant and ii) a dummy indicating whether \(t \leq 1992Q4\), we simply set
\[ d_{T+1}'=\ldots=d_{T+H}'=\begin{pmatrix} 1 & 0 \end{pmatrix} \]
H <- 12
(d_pred <- cbind(rep(1, 12), 0))
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 1 0
#> [3,] 1 0
#> [4,] 1 0
#> [5,] 1 0
#> [6,] 1 0
#> [7,] 1 0
#> [8,] 1 0
#> [9,] 1 0
#> [10,] 1 0
#> [11,] 1 0
#> [12,] 1 0We can now fit the model
bvar_obj <- fit(bvar_obj,
H = H,
d_pred = d_pred,
iter = 15000,
warmup = 5000,
chains = 4,
cores = 4)Let us look at the posterior means of \(\beta\), \(\Psi\), and \(\Sigma_u\).
summary(bvar_obj)
#> Posterior mean estimates
#> ------------------------
#>
#>
#> beta
#> --------------------------------------------------------------------------------
#> delta y_f pi_f i_f delta y pi i q
#> delta y_f.l1 0.18 0.03 -0.01 0.12 0.07 -0.12 0.00
#> pi_f.l1 -0.02 0.32 0.25 0.12 -0.07 0.01 0.00
#> i_f.l1 0.00 0.04 0.92 -0.04 0.06 0.05 0.00
#> delta y.l1 0.00 0.00 0.00 0.23 -0.09 -0.10 0.00
#> pi.l1 0.00 0.00 0.00 0.00 0.08 0.06 0.00
#> i.l1 0.00 0.00 0.00 0.00 0.02 0.76 0.00
#> q.l1 0.00 0.00 0.00 1.21 3.95 0.72 0.93
#> delta y_f.l2 0.03 -0.01 0.09 0.02 -0.02 0.10 0.00
#> pi_f.l2 0.01 0.02 0.04 0.00 -0.03 -0.15 0.00
#> i_f.l2 -0.02 -0.01 -0.01 0.00 0.04 0.07 0.00
#> delta y.l2 0.00 0.00 0.00 0.11 -0.01 0.15 0.00
#> pi.l2 0.00 0.00 0.00 0.01 -0.04 -0.05 0.00
#> i.l2 0.00 0.00 0.00 -0.01 0.01 0.04 0.00
#> q.l2 0.00 0.00 0.00 0.56 -0.39 0.32 -0.04
#> delta y_f.l3 0.01 -0.01 0.00 0.02 -0.01 0.00 0.00
#> pi_f.l3 -0.02 0.06 -0.01 0.00 0.08 0.02 0.00
#> i_f.l3 0.00 0.00 0.02 0.00 0.00 0.03 0.00
#> delta y.l3 0.00 0.00 0.00 0.06 0.01 -0.02 0.00
#> pi.l3 0.00 0.00 0.00 0.00 0.02 -0.02 0.00
#> i.l3 0.00 0.00 0.00 0.01 0.00 0.00 0.00
#> q.l3 0.00 0.00 0.00 -0.14 -0.02 -0.58 0.00
#> delta y_f.l4 0.03 -0.01 0.00 -0.01 0.03 0.02 0.00
#> pi_f.l4 0.00 0.16 -0.03 0.00 0.01 0.02 0.00
#> i_f.l4 0.00 0.00 -0.02 0.00 0.00 0.03 0.00
#> delta y.l4 0.00 0.00 0.00 -0.08 0.01 0.03 0.00
#> pi.l4 0.00 0.00 0.00 0.00 0.06 -0.01 0.00
#> i.l4 0.00 0.00 0.00 0.00 -0.01 0.00 0.00
#> q.l4 0.00 0.00 0.00 -0.15 -0.06 -0.18 -0.01
#> --------------------------------------------------------------------------------
#>
#>
#> Psi
#> --------------------------------------------------------------------------------
#> [,1] [,2]
#> delta y_f 0.58 0.08
#> pi_f 0.50 0.46
#> i_f 4.95 2.02
#> delta y 0.58 -0.04
#> pi 0.49 1.15
#> i 4.29 4.46
#> q 3.92 -0.10
#> --------------------------------------------------------------------------------
#>
#>
#> Sigma_u
#> --------------------------------------------------------------------------------
#> delta y_f pi_f i_f delta y pi i q
#> delta y_f 0.15 -0.01 0.01 0.07 -0.01 0.00 0.00
#> pi_f -0.01 0.09 0.05 0.01 0.13 0.04 0.00
#> i_f 0.01 0.05 0.52 0.01 0.18 0.11 0.00
#> delta y 0.07 0.01 0.01 0.19 -0.05 -0.01 0.00
#> pi -0.01 0.13 0.18 -0.05 0.60 0.11 0.00
#> i 0.00 0.04 0.11 -0.01 0.11 1.56 -0.01
#> q 0.00 0.00 0.00 0.00 0.00 -0.01 0.00
#> --------------------------------------------------------------------------------We can access the posterior means or medians with
bvar_obj$fit$posterior_means/bvar_obj$fit$posterior_medians
if needed.
Note that bvar_obj$fit$stan is an object of class
stanfit.
(stanfit <- bvar_obj$fit$stan)
#> Inference for Stan model: steady_state_bvar_homoscedastic_jeffreys_prior.
#> 4 chains, each with iter=15000; warmup=5000; thin=1;
#> post-warmup draws per chain=10000, total post-warmup draws=40000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> beta[1,1] 0.18 0.00 0.09 0.00 0.12 0.18 0.24 0.36 43372 1
#> beta[1,2] 0.03 0.00 0.05 -0.07 0.00 0.03 0.06 0.13 42966 1
#> beta[1,3] -0.01 0.00 0.13 -0.27 -0.10 -0.01 0.07 0.24 44014 1
#> beta[1,4] 0.12 0.00 0.08 -0.04 0.06 0.12 0.17 0.28 43888 1
#> beta[1,5] 0.07 0.00 0.14 -0.20 -0.02 0.07 0.17 0.34 40661 1
#> beta[1,6] -0.12 0.00 0.24 -0.60 -0.29 -0.12 0.04 0.36 41786 1
#> beta[1,7] 0.00 0.00 0.01 -0.01 -0.01 0.00 0.00 0.01 43618 1
#> beta[2,1] -0.02 0.00 0.09 -0.20 -0.08 -0.02 0.04 0.16 38569 1
#> beta[2,2] 0.32 0.00 0.08 0.16 0.26 0.32 0.37 0.48 33638 1
#> beta[2,3] 0.25 0.00 0.17 -0.08 0.14 0.25 0.36 0.58 40877 1
#> beta[2,4] 0.12 0.00 0.11 -0.09 0.05 0.12 0.19 0.33 37540 1
#> beta[2,5] -0.07 0.00 0.19 -0.45 -0.20 -0.07 0.06 0.30 33703 1
#> beta[2,6] 0.01 0.00 0.32 -0.61 -0.20 0.01 0.22 0.64 37478 1
#> beta[2,7] 0.00 0.00 0.01 -0.01 0.00 0.00 0.01 0.02 41142 1
#> beta[3,1] 0.00 0.00 0.03 -0.06 -0.03 0.00 0.02 0.05 29327 1
#> beta[3,2] 0.04 0.00 0.02 0.00 0.03 0.04 0.05 0.08 31281 1
#> beta[3,3] 0.92 0.00 0.08 0.78 0.87 0.92 0.97 1.07 29969 1
#> beta[3,4] -0.04 0.00 0.04 -0.11 -0.06 -0.04 -0.01 0.04 30014 1
#> beta[3,5] 0.06 0.00 0.06 -0.07 0.01 0.06 0.10 0.18 28920 1
#> beta[3,6] 0.05 0.00 0.11 -0.16 -0.03 0.05 0.12 0.26 31673 1
#> beta[3,7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 41230 1
#> beta[4,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39811 1
#> beta[4,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 40087 1
#> beta[4,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 37275 1
#> beta[4,4] 0.23 0.00 0.09 0.06 0.17 0.23 0.29 0.40 36298 1
#> beta[4,5] -0.09 0.00 0.12 -0.31 -0.17 -0.09 -0.01 0.14 37044 1
#> beta[4,6] -0.10 0.00 0.21 -0.51 -0.24 -0.10 0.04 0.31 40857 1
#> beta[4,7] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.01 45653 1
#> beta[5,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 38488 1
#> beta[5,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39530 1
#> beta[5,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39367 1
#> beta[5,4] 0.00 0.00 0.04 -0.07 -0.02 0.00 0.03 0.08 41946 1
#> beta[5,5] 0.08 0.00 0.09 -0.09 0.02 0.08 0.13 0.24 36317 1
#> beta[5,6] 0.06 0.00 0.12 -0.18 -0.02 0.06 0.14 0.30 43330 1
#> beta[5,7] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 47198 1
#> beta[6,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39640 1
#> beta[6,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 38353 1
#> beta[6,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 38796 1
#> beta[6,4] 0.00 0.00 0.02 -0.04 -0.01 0.00 0.01 0.04 37409 1
#> beta[6,5] 0.02 0.00 0.04 -0.05 0.00 0.02 0.04 0.09 33549 1
#> beta[6,6] 0.76 0.00 0.08 0.60 0.70 0.76 0.82 0.92 33042 1
#> beta[6,7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 56982 1
#> beta[7,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39969 1
#> beta[7,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39871 1
#> beta[7,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 40016 1
#> beta[7,4] 1.21 0.00 0.87 -0.49 0.63 1.21 1.80 2.91 31015 1
#> beta[7,5] 3.95 0.01 1.45 1.07 2.99 3.95 4.93 6.78 32554 1
#> beta[7,6] 0.72 0.01 2.59 -4.37 -1.04 0.72 2.49 5.82 34625 1
#> beta[7,7] 0.93 0.00 0.08 0.78 0.88 0.93 0.99 1.09 30372 1
#> beta[8,1] 0.03 0.00 0.07 -0.11 -0.02 0.03 0.08 0.17 44844 1
#> beta[8,2] -0.01 0.00 0.03 -0.07 -0.03 -0.01 0.01 0.05 49647 1
#> beta[8,3] 0.09 0.00 0.08 -0.07 0.04 0.09 0.15 0.25 47743 1
#> beta[8,4] 0.02 0.00 0.05 -0.08 -0.01 0.02 0.06 0.12 49438 1
#> beta[8,5] -0.02 0.00 0.09 -0.18 -0.07 -0.02 0.04 0.15 47327 1
#> beta[8,6] 0.10 0.00 0.15 -0.20 0.00 0.10 0.20 0.39 48771 1
#> beta[8,7] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.01 58151 1
#> beta[9,1] 0.01 0.00 0.06 -0.11 -0.03 0.01 0.05 0.13 47192 1
#> beta[9,2] 0.02 0.00 0.07 -0.11 -0.02 0.02 0.07 0.15 36392 1
#> beta[9,3] 0.04 0.00 0.12 -0.18 -0.04 0.04 0.12 0.27 48022 1
#> beta[9,4] 0.00 0.00 0.07 -0.14 -0.05 0.00 0.05 0.13 44487 1
#> beta[9,5] -0.03 0.00 0.12 -0.27 -0.11 -0.03 0.06 0.21 41947 1
#> beta[9,6] -0.15 0.00 0.21 -0.56 -0.29 -0.15 -0.01 0.27 46068 1
#> beta[9,7] 0.00 0.00 0.00 -0.01 0.00 0.00 0.01 0.01 48038 1
#> beta[10,1] -0.02 0.00 0.02 -0.06 -0.03 -0.02 0.00 0.03 38416 1
#> beta[10,2] -0.01 0.00 0.02 -0.04 -0.02 -0.01 0.00 0.02 37739 1
#> beta[10,3] -0.01 0.00 0.07 -0.15 -0.06 -0.01 0.04 0.14 31411 1
#> beta[10,4] 0.00 0.00 0.03 -0.05 -0.02 0.00 0.02 0.05 39838 1
#> beta[10,5] 0.04 0.00 0.05 -0.05 0.01 0.04 0.08 0.13 39727 1
#> beta[10,6] 0.07 0.00 0.08 -0.09 0.02 0.07 0.12 0.22 40960 1
#> beta[10,7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 56744 1
#> beta[11,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39380 1
#> beta[11,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 38479 1
#> beta[11,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 40028 1
#> beta[11,4] 0.11 0.00 0.07 -0.02 0.07 0.11 0.16 0.25 41530 1
#> beta[11,5] -0.01 0.00 0.08 -0.16 -0.06 -0.01 0.04 0.14 46577 1
#> beta[11,6] 0.15 0.00 0.14 -0.11 0.06 0.15 0.24 0.41 44622 1
#> beta[11,7] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.01 57329 1
#> beta[12,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39725 1
#> beta[12,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39034 1
#> beta[12,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 40530 1
#> beta[12,4] 0.01 0.00 0.03 -0.04 -0.01 0.01 0.03 0.06 45232 1
#> beta[12,5] -0.04 0.00 0.07 -0.17 -0.09 -0.04 0.00 0.09 42941 1
#> beta[12,6] -0.05 0.00 0.08 -0.20 -0.10 -0.05 0.00 0.10 48569 1
#> beta[12,7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 66125 1
#> beta[13,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 36843 1
#> beta[13,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39557 1
#> beta[13,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 40081 1
#> beta[13,4] -0.01 0.00 0.01 -0.04 -0.02 -0.01 0.00 0.02 43450 1
#> beta[13,5] 0.01 0.00 0.03 -0.04 0.00 0.01 0.03 0.06 41867 1
#> beta[13,6] 0.04 0.00 0.07 -0.10 -0.01 0.04 0.09 0.18 35626 1
#> beta[13,7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 70095 1
#> beta[14,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 39852 1
#> beta[14,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 41219 1
#> beta[14,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 38070 1
#> beta[14,4] 0.56 0.00 0.65 -0.72 0.12 0.56 1.00 1.83 39484 1
#> beta[14,5] -0.39 0.01 1.15 -2.62 -1.17 -0.40 0.38 1.86 37149 1
#> beta[14,6] 0.32 0.01 1.99 -3.56 -1.02 0.32 1.66 4.20 40158 1
#> beta[14,7] -0.04 0.00 0.08 -0.18 -0.09 -0.04 0.01 0.11 32133 1
#> beta[15,1] 0.01 0.00 0.05 -0.10 -0.03 0.01 0.05 0.12 47825 1
#> beta[15,2] -0.01 0.00 0.02 -0.05 -0.03 -0.01 0.00 0.03 50999 1
#> [ reached 'max' / getOption("max.print") -- omitted 440 rows ]
#>
#> Samples were drawn using NUTS(diag_e) at Mon Jul 13 23:45:47 2026.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).As such, we can do the usual rstan inference on our
fitted model. Let’s look at some examples for i) the first own lag of
the domestic interest rate (for which we set the prior mean to 0.9) and
ii) the post-crisis steady-state coefficient of inflation (multiply the
steady-state coefficient by 4 to obtain the annualized rate).
rstan::plot(stanfit,
pars=c("beta[6,6]", "Psi[5,1]"),
plotfun="hist")
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
We can also look at the model forecasts directly with
rstan. Remember that we left out the last two
observations/quarters, so let us look at our forecasts of the domestic
interest rate and compare them with the actual values.
(Villani2009[103:104,6]) #true values
#> [1] 1.478503 1.563795
rstan::plot(stanfit,
pars=c("y_pred[1,6]", "y_pred[2,6]"),
show_density = TRUE,
ci_level = 0.68,
fill_color = "blue")
#> ci_level: 0.68 (68% intervals)
#> outer_level: 0.95 (95% intervals)
So the model overshot a bit, but the true values are within the 68% prediction interval. Now let us plot the forecasts along with the historical data. We will choose a 68% prediction interval (“pi”) and the mean of the predictive distribution as the point forecast. For variables in quarter-on-quarter growth rates, we transform the historical data and predictions to yearly growth rates with ‘growth_rate_idx’ where we specify the index of the growth rate variables in \(y_t\). Note that this is not annualization, but we are now computing \(100 [ \ln (y_{i,t}) - \ln (y_{i,t-4})]\), i.e. the annual growth rate, by summing up to fourth differences.
fcst <- forecast(bvar_obj,
pi = 0.68,
fcst_type = "mean",
growth_rate_idx = c(4,5),
plot_idx = c(4,5,6))
For further inspection, we can print the point forecasts
print(fcst$forecast)
#> delta y_f pi_f i_f delta y pi i q
#> [1,] 0.6275889 0.5214174 2.826812 2.778850 0.9334433 2.036170 3.994601
#> [2,] 0.6464186 0.4509689 3.023410 3.058557 0.9778028 2.059633 3.990917
#> [3,] 0.6338417 0.4223705 3.163762 3.360539 1.7655304 2.173208 3.986382
#> [4,] 0.6394785 0.4944814 3.275603 3.439110 1.7623808 2.281314 3.980957
#> [5,] 0.6374716 0.4503922 3.405294 3.386102 1.7818525 2.392088 3.975869
#> [6,] 0.6352452 0.4311385 3.512768 3.303213 1.8419788 2.487643 3.971450
#> [7,] 0.6327143 0.4294042 3.606008 3.204599 1.8792027 2.583871 3.967013
#> [8,] 0.6285290 0.4409114 3.689660 3.120441 1.8732139 2.685173 3.963117
#> [9,] 0.6229851 0.4346041 3.769048 3.038515 1.8823192 2.766184 3.959660
#> [10,] 0.6229126 0.4286920 3.838016 2.971024 1.8860708 2.843557 3.956554
#> [11,] 0.6207810 0.4343662 3.902176 2.909566 1.8827908 2.919086 3.953725
#> [12,] 0.6224303 0.4330368 3.958150 2.856633 1.8622518 2.986326 3.951284We can also perform conditional forecasting by following Algorithm 3.3.1 in Dieppe, Legrand, and van Roye (2016). Note that for the structural shocks, identification is based on the Cholesky factorisation. Also, please note the limitations of this method, see the detailed discussion in Section 5.4 of Dieppe, Legrand, and van Roye (2016).
Now suppose we are interested in the forecasts of the domestic interest rate \(i\) conditional on a scenario where domestic inflation \(\pi\) gets really high (post COVID type scenario). Economic theory says the short interest rate should rise.
First we set up our conditions/scenarios, i.e., which variables, which horizons, and which values the variables will take during those horizons. Our conditions are that \(\pi\) will follow a specified path, at forecast horizons \(h=1,\dots,H=12\).
conditions <- data.frame(
var = rep(5,12),
horizon = rep(1:12),
value = c(1.0,1.5,2.0,1.8, #Note: QoQ scale for inflation here
1.5,1.2,1.0,1.0,
rep(0.5,4))
)We then do the conditional forecasting. We again select a 68% PI (prediction interval) and the mean of the predictive distribution as the point forecast.
cond_fcst <- conditional_forecast(bvar_obj,
conditions,
pi=0.68,
fcst_type = "mean",
plot_idx = c(5,6),
growth_rate_idx = c(5))
The short interest rate rises more dramatically compared to the unconditional case. Makes sense.
Now for some impulse response analysis. We can choose between the orthogonalized impulse response function (OIRF) and the generalized impulse response function (GIRF). Similar to forecasting, we can choose either the mean or the median (the default is the median), and we can also transform the IRFs for the quarter-on-quarter growth rate variables to the annual/yearly scale.
irf <- IRF(bvar_obj,H=20,response=5,shock=6,type="median",method="OIRF",ci=0.95,growth_rate_idx=5)
irf <- IRF(bvar_obj,H=20,response=4,shock=6,type="median",method="GIRF",ci=0.95,growth_rate_idx=4)
References
Dieppe, A., Legrand, R., and van Roye, B. (2016). The BEAR toolbox. Working Paper Series, No. 1934. European Central Bank.
Karlsson, S. (2013). Forecasting with Bayesian vector autoregression. In: Elliott, G. and Timmermann, A. (eds), Handbook of Economic Forecasting. Elsevier B.V., Vol. 2, Part B, pp. 791–897.
Villani, M. (2009). Steady-state priors for vector autoregressions. Journal of Applied Econometrics, 24(4), pp. 630-650.