---
title: "Homoscedastic steady-state BVAR (Villani, 2009)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Homoscedastic steady-state BVAR (Villani, 2009)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---





Here 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.


``` r
library(SteadyStateBVAR)
data("Villani2009")
yt <- Villani2009
```
The 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.


``` r
yt <- ts(yt[1:102, ], start = start(yt), frequency = frequency(yt))
plot.ts(yt)
```

![](figure/HOMO-1-1.png)

Also, let us create the bvar object which we will use throughout here.


``` r
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}
$$


``` r
bp <- which(time(yt) == 1992.75)
dum_var <- c(rep(1,bp), rep(0,nrow(yt)-bp))
```

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$.



``` r
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$.


``` r
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`.


``` r
#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.



``` r
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).


``` r
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    1
```
We 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:


``` r
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}
$$


``` r
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    0
```

We can now fit the model


``` r
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$.


``` r
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`.


``` r
(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).


``` r
rstan::plot(stanfit,
            pars=c("beta[6,6]", "Psi[5,1]"),
            plotfun="hist")
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
```

![](figure/HOMO-2-1.png)

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.


``` r
(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)
```

![](figure/HOMO-3-1.png)

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.



``` r
fcst <- forecast(bvar_obj,
                 pi = 0.68,
                 fcst_type = "mean",
                 growth_rate_idx = c(4,5),
                 plot_idx = c(4,5,6))
```

![](figure/HOMO-4-1.png)![](figure/HOMO-4-2.png)![](figure/HOMO-4-3.png)

For further inspection, we can print the point forecasts


``` r
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.951284
```


We 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$.


``` r
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.


``` r
cond_fcst <- conditional_forecast(bvar_obj,
                                  conditions,
                                  pi=0.68,
                                  fcst_type = "mean",
                                  plot_idx = c(5,6),
                                  growth_rate_idx = c(5))
```

![](figure/HOMO-5-1.png)![](figure/HOMO-5-2.png)

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.



``` r
irf <- IRF(bvar_obj,H=20,response=5,shock=6,type="median",method="OIRF",ci=0.95,growth_rate_idx=5)
```

![](figure/HOMO-6-1.png)

``` r
irf <- IRF(bvar_obj,H=20,response=4,shock=6,type="median",method="GIRF",ci=0.95,growth_rate_idx=4)
```

![](figure/HOMO-6-2.png)


## 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. 
