---
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(5, k),
                     beta_phi            = (rep(5, k) - 1) * rep(0.1, 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.26  0.01  0.15
#>   u.l1           -0.09  1.17 -0.16
#>   r.l1            0.00 -0.01  1.04
#>   delta pi.l2    -0.27  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.28
#>   r        3.50
#> --------------------------------------------------------------------------------
#> 
#> 
#> Sigma_u,t (t = 215)
#> --------------------------------------------------------------------------------
#>          delta pi     u     r
#> delta pi     0.09 -0.01  0.03
#> u           -0.01  0.02 -0.01
#> r            0.03 -0.01  0.16
#> --------------------------------------------------------------------------------
#> 
#> 
#> A
#> --------------------------------------------------------------------------------          
#>            delta pi   u r
#>   delta pi     1.00 0.0 0
#>   u            0.12 1.0 0
#>   r           -0.22 0.5 1
#> --------------------------------------------------------------------------------
#> 
#> 
#> phi
#> --------------------------------------------------------------------------------
#> delta pi        u        r 
#>     0.06     0.09     0.11 
#> --------------------------------------------------------------------------------
```
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.26    0.00  0.06  1.15  1.22  1.26  1.30   1.37  7006    1
#> beta[1,2]              0.01    0.00  0.04 -0.06 -0.01  0.01  0.04   0.09  7865    1
#> beta[1,3]              0.15    0.00  0.08  0.00  0.09  0.15  0.20   0.31  6222    1
#> beta[2,1]             -0.09    0.00  0.03 -0.16 -0.12 -0.09 -0.07  -0.03  7176    1
#> beta[2,2]              1.17    0.00  0.06  1.06  1.13  1.17  1.21   1.28  7552    1
#> beta[2,3]             -0.16    0.00  0.08 -0.32 -0.21 -0.16 -0.10  -0.01  6771    1
#> beta[3,1]              0.00    0.00  0.02 -0.03 -0.01  0.00  0.01   0.04  8959    1
#> beta[3,2]             -0.01    0.00  0.02 -0.04 -0.02 -0.01  0.00   0.02  9530    1
#> beta[3,3]              1.04    0.00  0.06  0.92  1.00  1.04  1.08   1.16  7421    1
#> beta[4,1]             -0.27    0.00  0.06 -0.38 -0.31 -0.28 -0.23  -0.16  7153    1
#> beta[4,2]              0.02    0.00  0.04 -0.06 -0.01  0.02  0.04   0.09  8015    1
#> beta[4,3]             -0.10    0.00  0.08 -0.26 -0.16 -0.10 -0.05   0.05  6282    1
#> beta[5,1]              0.07    0.00  0.03  0.01  0.05  0.07  0.09   0.13  7279    1
#> beta[5,2]             -0.23    0.00  0.05 -0.33 -0.27 -0.23 -0.19  -0.12  7742    1
#> beta[5,3]              0.17    0.00  0.07  0.03  0.12  0.17  0.22   0.32  6986    1
#> beta[6,1]              0.00    0.00  0.02 -0.03 -0.01  0.00  0.01   0.03  8316    1
#> beta[6,2]              0.02    0.00  0.02 -0.01  0.01  0.02  0.03   0.05  8810    1
#> beta[6,3]             -0.11    0.00  0.06 -0.22 -0.15 -0.11 -0.07   0.01  7252    1
#> Psi[1,1]               2.00    0.00  0.05  1.89  1.96  2.00  2.03   2.10 16547    1
#> Psi[2,1]               4.28    0.00  0.17  3.94  4.16  4.28  4.40   4.60 13232    1
#> Psi[3,1]               3.50    0.00  0.32  2.87  3.29  3.51  3.72   4.12 12420    1
#> z[1,1]                -0.03    0.00  0.27 -0.52 -0.22 -0.04  0.14   0.51  9674    1
#> z[1,2]                 0.65    0.00  0.29  0.08  0.45  0.64  0.84   1.25  9043    1
#> z[1,3]                -0.66    0.00  0.34 -1.30 -0.90 -0.67 -0.44   0.04  9462    1
#> z[2,1]                -0.13    0.01  0.96 -1.99 -0.79 -0.14  0.50   1.75 15103    1
#> z[2,2]                 0.24    0.01  0.99 -1.69 -0.43  0.23  0.90   2.20 11692    1
#> z[2,3]                 0.00    0.01  0.97 -1.92 -0.65  0.00  0.67   1.95 13873    1
#> z[3,1]                -0.05    0.01  0.97 -1.91 -0.71 -0.05  0.62   1.83 13394    1
#> z[3,2]                 0.31    0.01  0.97 -1.61 -0.33  0.32  0.96   2.23 10142    1
#> z[3,3]                -0.05    0.01  0.96 -1.93 -0.72 -0.05  0.60   1.82 14505    1
#> z[4,1]                -0.09    0.01  0.97 -1.98 -0.74 -0.10  0.55   1.82 14529    1
#> z[4,2]                -0.47    0.01  0.96 -2.33 -1.12 -0.48  0.18   1.44 13354    1
#> z[4,3]                -0.24    0.01  0.93 -2.05 -0.86 -0.23  0.40   1.55 15334    1
#> z[5,1]                 0.01    0.01  1.00 -1.96 -0.65  0.02  0.68   1.99 17285    1
#> z[5,2]                -0.48    0.01  0.97 -2.43 -1.14 -0.46  0.18   1.45 15175    1
#> z[5,3]                -0.20    0.01  0.94 -2.09 -0.82 -0.19  0.42   1.67 14648    1
#> z[6,1]                 0.01    0.01  0.97 -1.86 -0.64  0.01  0.65   1.92 17221    1
#> z[6,2]                -0.35    0.01  0.96 -2.22 -1.00 -0.35  0.31   1.51 13492    1
#> z[6,3]                -0.13    0.01  0.94 -1.99 -0.76 -0.14  0.50   1.70 13985    1
#> z[7,1]                 0.14    0.01  0.98 -1.81 -0.51  0.14  0.79   2.05 14866    1
#> z[7,2]                -0.22    0.01  0.96 -2.08 -0.88 -0.21  0.43   1.67 13625    1
#> z[7,3]                -0.04    0.01  0.95 -1.93 -0.68 -0.03  0.62   1.80 14365    1
#> z[8,1]                 0.23    0.01  0.98 -1.73 -0.44  0.23  0.90   2.13 15389    1
#> z[8,2]                -0.25    0.01  0.94 -2.11 -0.87 -0.23  0.35   1.61 14216    1
#> z[8,3]                 0.07    0.01  0.94 -1.77 -0.57  0.08  0.71   1.86 14645    1
#> z[9,1]                 0.23    0.01  0.97 -1.67 -0.43  0.25  0.86   2.14 16471    1
#> z[9,2]                -0.12    0.01  0.95 -1.96 -0.76 -0.11  0.51   1.73 15187    1
#> z[9,3]                 0.22    0.01  0.94 -1.62 -0.42  0.22  0.83   2.04 13872    1
#> z[10,1]               -0.17    0.01  0.95 -2.02 -0.80 -0.17  0.44   1.71 16319    1
#> z[10,2]               -0.11    0.01  0.95 -1.99 -0.75 -0.11  0.53   1.76 15556    1
#> z[10,3]               -0.06    0.01  0.97 -1.98 -0.71 -0.10  0.60   1.87 15242    1
#> z[11,1]               -0.15    0.01  0.95 -1.99 -0.81 -0.15  0.49   1.73 16400    1
#> z[11,2]               -0.13    0.01  0.96 -2.01 -0.77 -0.13  0.53   1.74 15579    1
#> z[11,3]               -0.16    0.01  0.97 -2.05 -0.83 -0.17  0.49   1.77 14933    1
#> z[12,1]               -0.35    0.01  0.96 -2.28 -0.99 -0.36  0.30   1.50 15723    1
#> z[12,2]               -0.02    0.01  0.96 -1.89 -0.65 -0.02  0.61   1.88 15720    1
#> z[12,3]               -0.10    0.01  0.97 -2.04 -0.74 -0.10  0.55   1.79 14175    1
#> z[13,1]               -0.26    0.01  0.96 -2.11 -0.91 -0.26  0.38   1.62 17321    1
#> z[13,2]                0.12    0.01  0.95 -1.83 -0.51  0.12  0.74   1.98 15008    1
#> z[13,3]                0.04    0.01  0.95 -1.80 -0.59  0.03  0.69   1.91 15647    1
#> z[14,1]               -0.32    0.01  0.96 -2.22 -0.98 -0.33  0.32   1.55 15719    1
#> z[14,2]                0.16    0.01  0.95 -1.67 -0.50  0.16  0.82   2.01 15518    1
#> z[14,3]                0.03    0.01  0.97 -1.86 -0.61  0.04  0.67   1.93 16327    1
#> z[15,1]               -0.26    0.01  0.96 -2.12 -0.93 -0.26  0.39   1.62 15220    1
#> z[15,2]                0.07    0.01  0.96 -1.83 -0.58  0.08  0.73   1.96 15229    1
#> z[15,3]                0.11    0.01  0.95 -1.72 -0.54  0.11  0.78   1.96 15623    1
#> z[16,1]               -0.17    0.01  0.95 -2.03 -0.83 -0.17  0.47   1.67 14146    1
#> z[16,2]                0.15    0.01  0.98 -1.80 -0.52  0.15  0.82   2.06 15368    1
#> z[16,3]                0.26    0.01  0.94 -1.55 -0.38  0.26  0.89   2.08 14018    1
#> z[17,1]               -0.19    0.01  0.95 -2.04 -0.83 -0.18  0.46   1.69 17853    1
#> z[17,2]                0.19    0.01  0.93 -1.61 -0.45  0.20  0.83   1.98 14709    1
#> z[17,3]                0.31    0.01  0.94 -1.54 -0.32  0.31  0.95   2.18 14415    1
#> z[18,1]               -0.30    0.01  0.96 -2.14 -0.96 -0.30  0.34   1.59 14903    1
#> z[18,2]                0.33    0.01  0.96 -1.55 -0.32  0.32  0.96   2.24 16566    1
#> z[18,3]                0.43    0.01  0.96 -1.44 -0.23  0.43  1.08   2.28 13465    1
#> z[19,1]               -0.16    0.01  0.96 -2.08 -0.78 -0.15  0.48   1.73 17488    1
#> z[19,2]                0.44    0.01  0.94 -1.41 -0.19  0.47  1.07   2.24 16186    1
#> z[19,3]                0.39    0.01  0.93 -1.42 -0.22  0.40  1.02   2.20 15204    1
#> z[20,1]               -0.17    0.01  0.97 -2.10 -0.81 -0.18  0.47   1.78 14318    1
#> z[20,2]                0.07    0.01  0.92 -1.75 -0.55  0.06  0.69   1.86 15464    1
#> z[20,3]                0.33    0.01  0.92 -1.46 -0.31  0.33  0.95   2.15 15141    1
#> z[21,1]               -0.06    0.01  0.97 -2.02 -0.71 -0.06  0.60   1.83 13342    1
#> z[21,2]               -0.44    0.01  0.96 -2.32 -1.10 -0.44  0.21   1.41 15241    1
#> z[21,3]                0.37    0.01  0.96 -1.56 -0.27  0.36  1.01   2.27 13459    1
#> z[22,1]               -0.07    0.01  0.99 -2.00 -0.71 -0.07  0.60   1.87 15725    1
#> z[22,2]               -0.31    0.01  0.96 -2.18 -0.95 -0.30  0.35   1.57 13059    1
#> z[22,3]                0.52    0.01  0.94 -1.32 -0.10  0.53  1.17   2.36 16358    1
#> z[23,1]               -0.02    0.01  0.97 -1.95 -0.67 -0.01  0.63   1.86 12680    1
#> z[23,2]               -0.36    0.01  0.98 -2.27 -1.02 -0.36  0.31   1.51 16821    1
#> z[23,3]               -0.18    0.01  0.94 -2.04 -0.81 -0.16  0.46   1.62 16186    1
#> z[24,1]               -0.20    0.01  0.96 -2.10 -0.83 -0.19  0.44   1.69 17884    1
#> z[24,2]               -0.25    0.01  0.94 -2.10 -0.88 -0.26  0.36   1.58 14659    1
#> z[24,3]               -0.07    0.01  0.95 -1.91 -0.72 -0.07  0.59   1.78 14450    1
#> z[25,1]               -0.19    0.01  0.97 -2.06 -0.85 -0.19  0.49   1.67 14051    1
#> z[25,2]               -0.25    0.01  0.92 -2.06 -0.88 -0.25  0.37   1.57 15451    1
#> z[25,3]                0.07    0.01  0.92 -1.73 -0.56  0.07  0.69   1.87 11939    1
#> z[26,1]               -0.04    0.01  0.98 -1.99 -0.68 -0.04  0.60   1.91 15991    1
#> z[26,2]               -0.14    0.01  0.94 -2.02 -0.76 -0.12  0.51   1.68 17109    1
#> z[26,3]                0.20    0.01  0.94 -1.61 -0.44  0.20  0.85   2.02 14011    1
#> z[27,1]               -0.08    0.01  0.97 -1.99 -0.73 -0.08  0.57   1.76 16457    1
#>  [ reached 'max' / getOption("max.print") -- omitted 3759 rows ]
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Jul 16 04:38:31 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.
