---
title: "Random Walk stochastic volatility steady-state BVAR (Clark, 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Random Walk stochastic volatility steady-state BVAR (Clark, 2011)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



Here we estimate the steady-state BVAR model with Random Walk stochastic volatility from Clark (2011), which is an extension of the original homoscedastic steady-state BVAR model (Villani, 2009). See `?bvar` for details.

We will estimate the model on a quarterly US data set from Koop and Korobilis (2010) on the inflation rate $\Delta \pi_t$ (the annual percentage change in a chain-weighted GDP price index), the unemployment rate $u_t$ (seasonally adjusted civilian unemployment rate, all civilian workers aged 16 years or older) and the interest rate $r_t$ (yield on the three-month Treasury bill rate). The sample is 1953Q1-2006Q3 and we have the data vector

$$
y_t = 
\begin{pmatrix} \Delta \pi_t \\
u_t \\
r_t
\end{pmatrix}
$$

First, let's load the package, then import and plot the data.


``` r
library(SteadyStateBVAR)
data("KoopKorobilis2010")
yt <- KoopKorobilis2010
plot.ts(yt)
```

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

Let's create the bvar object which we will use throughout here.


``` r
bvar_obj <- bvar(data = yt)
```

We choose 2 lags and only a constant as the deterministic variable.


``` r
bvar_obj <- setup(bvar_obj,
                  p=2,
                  deterministic = "constant")
```

We set the overall tightness to $\lambda_1 = 0.20$, cross-equation tightness to $\lambda_2 = 0.50$ and the lag decay rate to $\lambda_3 = 1.00$. For the prior means on the first own lags, we set them to $0.6$ for $\Delta \pi_t$ and $0.9$ for $u_t$ and $r_t$. Note that the prior mean on the first own lag of inflation is set to $0.6$ instead of $0$ to reflect some degree of persistence in the series (even though it is a growth rate variable).



``` r
lambda_1 <- 0.20
lambda_2 <- 0.50
lambda_3 <- 1.00

fol_pm=c(0.6, # delta pi
         0.9,  #u
         0.9)  #R
```

Now, for the steady-state coefficients we use some toy values (let us pretend that they are expert based).
Remember that we only have a constant now, so $q=1$ and therefore $\Psi$ only has one column $\psi_1=\Psi$. Since $d_t = 1 \ \forall \ t$, we have $\Psi d_t = \mu_t$ which simplifies to $\Psi = \mu$ and as such we can directly interpret $\Psi$ as the unconditional mean, i.e. the steady state.


``` r
theta_Psi <- 
  c(
  ppi(1.90, 2.10, interval=0.95)$mean,   #Psi: delta pi
  ppi(3.80, 4.50, interval=0.95)$mean,   #Psi: u
  ppi(2.60, 3.90, interval=0.95)$mean    #Psi: r
  )

Omega_Psi <- 
  diag(
  c(
  ppi(1.90, 2.10, interval=0.95)$var,    #Psi: delta pi
  ppi(3.80, 4.50, interval=0.95)$var,    #Psi: u
  ppi(2.60, 3.90, interval=0.95)$var     #Psi: r
  )
  )
```

Now we need to specify our stochastic volatility priors. See `?priors` for more information about the prior specification. I take some inspiration from Clark (2011) below.


``` r
k <- bvar_obj$setup$k
n_free_params_A <- bvar_obj$setup$n_free_params_A
sigma2 <- diag(bvar_obj$setup$Sigma_AR)

SV_priors_RW <- list(
                     theta_A             =  rep(0, n_free_params_A),
                     Omega_A             =  diag(10, n_free_params_A),
                     mu_log_lambda_1     =  log(sigma2),
                     sigma2_log_lambda_1 =  rep(4, k),
                     alpha_phi           =  rep(2.5, k),
                     beta_phi            =  rep(0.0875, k)
                    )
```

Here `sigma2` contains the residual variances from AR($p$) models (the same ones we used in the Minnesota prior).

Let's put everything into the `priors()` function.


``` r
bvar_obj <- priors(bvar_obj,
                   lambda_1 = lambda_1,
                   lambda_2 = lambda_2,
                   lambda_3 = lambda_3,
                   first_own_lag_prior_mean =fol_pm,
                   theta_Psi = theta_Psi,
                   Omega_Psi = Omega_Psi,
                   SV = TRUE,
                   SV_type = "RW",
                   SV_priors = SV_priors_RW)
```

Now we can fit the model. Note that we can use arguments from `rstan::sampling()` such as `control` where we can tweak `max_treedepth` and `adapt_delta`.


``` r
bvar_obj <- fit(bvar_obj,
                H = 40,
                d_pred = matrix(rep(1, 40)),
                iter = 4000,
                warmup = 1000,
                chains = 2,
                cores = 2,
                control = list(max_treedepth = 14, adapt_delta = 0.95))
```

Now let's see the posterior means


``` r
summary(bvar_obj, stat="mean", t = 215) #t = 215 for covariance matrix
#> Posterior mean estimates
#> ------------------------
#> 
#> 
#> beta
#> --------------------------------------------------------------------------------             
#>               delta pi     u     r
#>   delta pi.l1     1.27  0.02  0.15
#>   u.l1           -0.09  1.17 -0.16
#>   r.l1            0.00 -0.01  1.04
#>   delta pi.l2    -0.28  0.02 -0.10
#>   u.l2            0.07 -0.23  0.17
#>   r.l2            0.00  0.02 -0.11
#> --------------------------------------------------------------------------------
#> 
#> 
#> Psi
#> --------------------------------------------------------------------------------          
#>            [,1]
#>   delta pi 2.00
#>   u        4.29
#>   r        3.50
#> --------------------------------------------------------------------------------
#> 
#> 
#> Sigma_u,t (t = 215)
#> --------------------------------------------------------------------------------
#>          delta pi     u     r
#> delta pi     0.09 -0.01  0.02
#> u           -0.01  0.02 -0.01
#> r            0.02 -0.01  0.16
#> --------------------------------------------------------------------------------
#> 
#> 
#> A
#> --------------------------------------------------------------------------------          
#>            delta pi    u r
#>   delta pi     1.00 0.00 0
#>   u            0.12 1.00 0
#>   r           -0.22 0.51 1
#> --------------------------------------------------------------------------------
#> 
#> 
#> phi
#> --------------------------------------------------------------------------------
#> delta pi        u        r 
#>     0.04     0.07     0.10 
#> --------------------------------------------------------------------------------
```
You can always look at the `stanfit` object `bvar_obj$fit$stan` directly if you want. Note that
the `z`'s below are not parameters per se, they are simply used in a reparameterization trick to sample
the log volatilities more efficiently.


``` r
print(bvar_obj$fit$stan)
#> Inference for Stan model: steady_state_bvar_RW_stochastic_volatility.
#> 2 chains, each with iter=4000; warmup=1000; thin=1; 
#> post-warmup draws per chain=3000, total post-warmup draws=6000.
#> 
#>                        mean se_mean     sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
#> beta[1,1]              1.27    0.00   0.06  1.16  1.23  1.27  1.31  1.38 10297    1
#> beta[1,2]              0.02    0.00   0.04 -0.06 -0.01  0.02  0.04  0.09  8494    1
#> beta[1,3]              0.15    0.00   0.08 -0.01  0.10  0.15  0.20  0.30  7773    1
#> beta[2,1]             -0.09    0.00   0.03 -0.16 -0.12 -0.09 -0.07 -0.02  8907    1
#> beta[2,2]              1.17    0.00   0.06  1.06  1.13  1.17  1.21  1.28  8749    1
#> beta[2,3]             -0.16    0.00   0.08 -0.32 -0.21 -0.16 -0.11 -0.01  8377    1
#> beta[3,1]              0.00    0.00   0.02 -0.03 -0.01  0.00  0.01  0.04 10036    1
#> beta[3,2]             -0.01    0.00   0.02 -0.04 -0.02 -0.01  0.00  0.02 10459    1
#> beta[3,3]              1.04    0.00   0.06  0.92  1.00  1.04  1.08  1.15 10192    1
#> beta[4,1]             -0.28    0.00   0.06 -0.39 -0.32 -0.28 -0.24 -0.17 10144    1
#> beta[4,2]              0.02    0.00   0.04 -0.06 -0.01  0.02  0.04  0.09  8920    1
#> beta[4,3]             -0.10    0.00   0.08 -0.26 -0.16 -0.10 -0.05  0.05  7662    1
#> beta[5,1]              0.07    0.00   0.03  0.01  0.05  0.07  0.09  0.14  9346    1
#> beta[5,2]             -0.23    0.00   0.05 -0.34 -0.27 -0.23 -0.20 -0.12  8545    1
#> beta[5,3]              0.17    0.00   0.07  0.03  0.12  0.17  0.22  0.32  9076    1
#> beta[6,1]              0.00    0.00   0.01 -0.03 -0.01  0.00  0.01  0.03  9794    1
#> beta[6,2]              0.02    0.00   0.02 -0.01  0.01  0.02  0.03  0.05 11413    1
#> beta[6,3]             -0.11    0.00   0.06 -0.22 -0.15 -0.11 -0.07  0.01  9481    1
#> Psi[1,1]               2.00    0.00   0.05  1.90  1.96  2.00  2.03  2.10 15316    1
#> Psi[2,1]               4.29    0.00   0.18  3.93  4.17  4.30  4.41  4.65 15661    1
#> Psi[3,1]               3.50    0.00   0.32  2.87  3.28  3.50  3.72  4.11 13141    1
#> z[1,1]                -0.05    0.00   0.25 -0.50 -0.22 -0.06  0.11  0.47 10574    1
#> z[1,2]                 0.64    0.00   0.27  0.15  0.46  0.64  0.81  1.19 10220    1
#> z[1,3]                -0.65    0.00   0.34 -1.29 -0.88 -0.66 -0.43  0.03  9569    1
#> z[2,1]                -0.13    0.01   0.95 -1.97 -0.79 -0.13  0.52  1.75 15812    1
#> z[2,2]                 0.18    0.01   0.99 -1.76 -0.49  0.18  0.86  2.11 13629    1
#> z[2,3]                 0.00    0.01   1.00 -1.95 -0.69  0.00  0.68  1.96 11939    1
#> z[3,1]                -0.05    0.01   0.99 -1.98 -0.71 -0.04  0.62  1.85 15000    1
#> z[3,2]                 0.25    0.01   0.98 -1.70 -0.41  0.24  0.91  2.21 12218    1
#> z[3,3]                -0.04    0.01   1.00 -2.00 -0.71 -0.04  0.65  1.90 15046    1
#> z[4,1]                -0.09    0.01   0.98 -1.97 -0.75 -0.10  0.59  1.84 14747    1
#> z[4,2]                -0.43    0.01   0.92 -2.19 -1.06 -0.43  0.19  1.36 16494    1
#> z[4,3]                -0.26    0.01   0.95 -2.07 -0.89 -0.27  0.36  1.67 14793    1
#> z[5,1]                 0.00    0.01   0.97 -1.90 -0.66 -0.01  0.66  1.92 14608    1
#> z[5,2]                -0.46    0.01   0.95 -2.32 -1.11 -0.46  0.19  1.38 16332    1
#> z[5,3]                -0.21    0.01   1.00 -2.13 -0.89 -0.22  0.46  1.72 15558    1
#> z[6,1]                -0.01    0.01   0.99 -2.00 -0.69 -0.01  0.66  1.96 13342    1
#> z[6,2]                -0.34    0.01   0.97 -2.26 -1.00 -0.33  0.33  1.60 15816    1
#> z[6,3]                -0.13    0.01   0.97 -2.02 -0.80 -0.13  0.53  1.73 14987    1
#> z[7,1]                 0.10    0.01   0.95 -1.78 -0.54  0.10  0.76  1.94 20300    1
#> z[7,2]                -0.21    0.01   0.98 -2.13 -0.86 -0.22  0.46  1.74 14686    1
#> z[7,3]                -0.04    0.01   0.95 -1.92 -0.67 -0.03  0.61  1.85 13893    1
#> z[8,1]                 0.19    0.01   0.99 -1.79 -0.47  0.17  0.86  2.12 19128    1
#> z[8,2]                -0.25    0.01   0.97 -2.18 -0.90 -0.26  0.39  1.61 15872    1
#> z[8,3]                 0.06    0.01   0.92 -1.75 -0.57  0.05  0.69  1.88 12407    1
#> z[9,1]                 0.16    0.01   0.97 -1.74 -0.49  0.16  0.84  2.07 15461    1
#> z[9,2]                -0.15    0.01   0.96 -2.02 -0.79 -0.15  0.51  1.71 16838    1
#> z[9,3]                 0.19    0.01   0.97 -1.72 -0.48  0.19  0.85  2.03 15960    1
#> z[10,1]               -0.18    0.01   0.97 -2.08 -0.85 -0.18  0.48  1.66 15155    1
#> z[10,2]               -0.12    0.01   0.97 -2.01 -0.78 -0.11  0.55  1.82 15984    1
#> z[10,3]               -0.04    0.01   0.97 -1.98 -0.69 -0.03  0.59  1.88 18843    1
#> z[11,1]               -0.17    0.01   0.94 -2.04 -0.81 -0.17  0.47  1.71 13422    1
#> z[11,2]               -0.14    0.01   0.95 -2.00 -0.79 -0.14  0.50  1.74 14516    1
#> z[11,3]               -0.17    0.01   0.97 -2.03 -0.82 -0.18  0.49  1.71 15496    1
#> z[12,1]               -0.32    0.01   0.96 -2.21 -0.97 -0.32  0.31  1.53 13949    1
#> z[12,2]               -0.04    0.01   0.96 -1.91 -0.68 -0.04  0.60  1.86 15602    1
#> z[12,3]               -0.09    0.01   0.94 -1.92 -0.73 -0.09  0.55  1.74 14597    1
#> z[13,1]               -0.29    0.01   0.98 -2.22 -0.95 -0.29  0.39  1.64 16269    1
#> z[13,2]                0.09    0.01   0.97 -1.80 -0.55  0.08  0.73  2.00 14781    1
#> z[13,3]                0.04    0.01   0.96 -1.81 -0.61  0.04  0.68  1.92 14488    1
#> z[14,1]               -0.32    0.01   0.97 -2.16 -0.98 -0.32  0.33  1.57 15088    1
#> z[14,2]                0.13    0.01   0.95 -1.74 -0.51  0.14  0.78  1.96 15283    1
#> z[14,3]                0.03    0.01   0.96 -1.84 -0.62  0.02  0.65  1.95 16390    1
#> z[15,1]               -0.25    0.01   0.99 -2.19 -0.90 -0.27  0.42  1.70 14537    1
#> z[15,2]                0.06    0.01   0.96 -1.87 -0.58  0.07  0.71  1.90 14451    1
#> z[15,3]                0.11    0.01   0.97 -1.77 -0.55  0.12  0.78  1.97 14430    1
#> z[16,1]               -0.21    0.01   0.96 -2.06 -0.87 -0.22  0.44  1.68 15286    1
#> z[16,2]                0.12    0.01   0.98 -1.82 -0.53  0.13  0.78  2.06 15161    1
#> z[16,3]                0.26    0.01   0.96 -1.64 -0.39  0.26  0.90  2.15 13452    1
#> z[17,1]               -0.21    0.01   0.97 -2.10 -0.87 -0.22  0.45  1.71 14900    1
#> z[17,2]                0.18    0.01   0.95 -1.67 -0.46  0.19  0.83  2.04 13575    1
#> z[17,3]                0.30    0.01   0.96 -1.57 -0.35  0.30  0.94  2.17 12809    1
#> z[18,1]               -0.28    0.01   0.96 -2.18 -0.94 -0.28  0.37  1.57 16645    1
#> z[18,2]                0.26    0.01   0.95 -1.60 -0.39  0.26  0.90  2.13 16719    1
#> z[18,3]                0.43    0.01   0.95 -1.42 -0.22  0.44  1.08  2.27 12369    1
#> z[19,1]               -0.20    0.01   0.97 -2.09 -0.88 -0.22  0.47  1.71 15822    1
#> z[19,2]                0.38    0.01   0.95 -1.50 -0.25  0.37  1.01  2.24 15054    1
#> z[19,3]                0.37    0.01   0.97 -1.50 -0.29  0.37  1.02  2.27 13218    1
#> z[20,1]               -0.19    0.01   0.98 -2.08 -0.85 -0.19  0.48  1.71 16124    1
#> z[20,2]                0.04    0.01   0.93 -1.79 -0.58  0.04  0.67  1.79 14503    1
#> z[20,3]                0.31    0.01   0.99 -1.61 -0.37  0.32  1.00  2.26 17200    1
#> z[21,1]               -0.12    0.01   0.96 -2.04 -0.74 -0.11  0.51  1.79 14224    1
#> z[21,2]               -0.43    0.01   0.97 -2.30 -1.09 -0.44  0.20  1.51 15180    1
#> z[21,3]                0.38    0.01   0.94 -1.43 -0.27  0.38  1.02  2.21 14878    1
#> z[22,1]               -0.11    0.01   0.97 -1.97 -0.77 -0.12  0.55  1.79 14179    1
#> z[22,2]               -0.32    0.01   0.94 -2.17 -0.95 -0.34  0.30  1.54 14549    1
#> z[22,3]                0.50    0.01   0.94 -1.34 -0.14  0.52  1.14  2.35 16356    1
#> z[23,1]               -0.06    0.01   0.98 -2.01 -0.72 -0.05  0.59  1.83 15152    1
#> z[23,2]               -0.36    0.01   0.98 -2.28 -1.03 -0.36  0.31  1.57 13117    1
#> z[23,3]               -0.17    0.01   0.95 -2.07 -0.81 -0.17  0.46  1.67 16037    1
#> z[24,1]               -0.19    0.01   0.94 -2.02 -0.82 -0.19  0.44  1.67 15235    1
#> z[24,2]               -0.29    0.01   0.95 -2.15 -0.93 -0.29  0.36  1.52 13414    1
#> z[24,3]               -0.07    0.01   0.94 -1.92 -0.70 -0.08  0.57  1.78 15798    1
#> z[25,1]               -0.17    0.01   0.97 -2.10 -0.81 -0.18  0.47  1.76 18826    1
#> z[25,2]               -0.29    0.01   0.95 -2.15 -0.92 -0.30  0.34  1.57 13711    1
#> z[25,3]                0.06    0.01   0.95 -1.81 -0.59  0.07  0.71  1.88 13770    1
#> z[26,1]               -0.05    0.01   0.98 -2.02 -0.71 -0.04  0.61  1.87 16092    1
#> z[26,2]               -0.19    0.01   0.95 -2.04 -0.86 -0.17  0.45  1.66 13991    1
#> z[26,3]                0.20    0.01   0.96 -1.70 -0.45  0.21  0.84  2.07 14431    1
#> z[27,1]               -0.10    0.01   0.92 -1.94 -0.69 -0.09  0.52  1.76 14331    1
#>  [ reached 'max' / getOption("max.print") -- omitted 3759 rows ]
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jul 28 03:46:18 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).
```

We can forecast


``` r
forecast(bvar_obj, pi = 0.68, show_all = TRUE)
```

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

Let us plot the log volatility estimates and predictions


``` r
stochastic_volatility_plot(bvar_obj, ci = 0.95, vol = "log_lambda")
```

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

Let us plot the estimates and predictions of the implied innovation standard deviations


``` r
stochastic_volatility_plot(bvar_obj, vol = "sd")
```

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

We can also produce orthogonalized IRFs


``` r
IRF(bvar_obj, method = "OIRF", t=215, ci=0.68) #latest t
```

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


## References

Clark, T. E. (2011). Real-time density forecasts from Bayesian vector autoregressions
with stochastic volatility. *Journal of Business \& Economic Statistics*, 29(3), pp. 327–341.

Koop, G. and Korobilis, D. (2010). Bayesian multivariate time series methods for empirical macroeconomics.
*Foundations and Trends in Econometrics*, 3(4), pp. 267–358. 

Villani, M. (2009). Steady-state priors for vector autoregressions. *Journal of Applied Econometrics*, 24(4), pp. 630–650.
