SteadyStateBVAR

With this package, the user can estimate the steady-state BVAR model of Villani (2009). The steady-state BVAR is simply a BVAR rewritten in mean-adjusted form. The benefit of the mean-adjusted parametrization is that it allows the user to specify prior beliefs about the unconditional mean, or steady state of the VAR system. The model has proven very useful for forecasting of macroeconomic variables, and is routinely used in many central banks and other finanicial institutions (Gustafsson and Villani, 2025).

After estimation, the user can produce forecasts (unconditional and conditional) and impulse response functions (orthogonalized and generalized). The goal of SteadyStateBVAR is to use modern Bayesian tools (Stan) to: i) estimate the model as specified in the original paper, and ii) extend the model in different ways. Previously, extensions of the model seemed to be limited by what Mattias Villani had time to derive.

See for example Clark (2011), which extends the steady-state BVAR model to include Random Walk stochastic volatility: “In a methodological sense, this paper extends the estimator of Villani (2009) to include stochastic volatility.” Later on, he writes: “(Special thanks are due to Mattias Villani for providing the formulas for posterior means and variances of $\Pi$ and $\Psi$, which generalize the constant-variance formulas of Villani 2009.)”

Another example is Dieppe, Legrand, and van Roye (2016): “Villani (2009) only provides derivation in the case of the normal-diffuse prior distribution, so the incoming analysis will be restricted to this case.”

Times are different now, and with the help of Stan, we are essentially limited only by our imagination. At the time of writing, the package provides three versions of the steady-state BVAR model: i) the homoscedastic model, i.e. the original model in Villani (2009), ii) the Random Walk stochastic volatility model, i.e. the one in Clark (2011), and iii) an AR(1) stochastic volatility model.

Installation

You can install SteadyStateBVAR with:

install.packages("SteadyStateBVAR")

You can also install the development version of SteadyStateBVAR from GitHub with:

# install.packages("pak")
pak::pak("markjwbecker/SteadyStateBVAR")

Introduction

The steady-state BVAR($p$) model (Villani, 2009) is

$$y_t = \Psi d_t + \Pi_1(y_{t-1}-\Psi d_{t-1})+\dots+\Pi_p(y_{t-p}-\Psi d_{t-p})+u_t$$

where $y_t$ is a $k$-dimensional vector of endogenous variables (time series) at time $t$, $d_t$ is a $q$-dimensional vector of deterministic (exogenous) variables at time $t$, and the (reduced-form) innovations are $u_t \overset{\text{iid}}{\sim} \mathrm{N_k}(0,\Sigma_u)$. Here $\Pi_\ell$ for $\ell=1,\dots,p$ is a $(k \times k)$ autoregressive parameter matrix, and $\Psi$ is a $(k \times q)$ steady-state parameter matrix. Now

$$\mathrm{E}(y_t)=\mu_t=\Psi d_t$$

is the unconditional mean, or the steady state, of the process. Since long-horizon forecasts from stationary VARs converge to the unconditional mean (steady state), it is naturally very important from a forecasting perspective to obtain precise inference on $\Psi$. Note that the current version of this package only allows for $d_t$ to contain either a constant, a constant and a dummy variable, or a constant and a time trend.

We may stack the (transposed) $\Pi_i$ matrices in the $(kp \times k)$ matrix $\beta$

$$\beta= \begin{bmatrix} \Pi'_1 \\ \vdots \\ \Pi'_p \end{bmatrix}$$

We can then rewrite the model as a nonlinear regression (Karlsson, 2013)

$$y_t' =d_t'\Psi' + \left[w_t'-q_t'(I_p \otimes \Psi') \right]\beta +u_t'$$

where $w_t'=(y_{t-1}',\dots,y_{t-p}')$ is a $kp$-dimensional vector of lagged endogenous variables and $q_t'=(d_{t-1}',\dots,d_{t-p}')$ is a $qp$-dimensional vector of lagged deterministic (exogenous) variables, $I_p$ is the $(p \times p)$ identity matrix and $\otimes$ denotes the Kronecker product. This is how the likelihood is written in the Stan code. The goal is to estimate $\beta, \Psi$, and $\Sigma_u$, and as such priors are needed. Following Villani (2009), prior independence between $\beta, \Psi$ and $\Sigma_u$ is assumed. For $\beta$, the Minnesota prior is used

$$\mathrm{vec}(\beta) \sim \mathrm{N}_{kpk} \left[\theta_\beta,\Omega_\beta\right]$$

The prior means (the elements of $\theta_\beta$) are set to

$$\begin{aligned} \mathrm{E}\left(\Pi_{\ell}^{(i,j)}\right)&= \begin{cases} \kappa & \text{if } \ell = 1 \ \mathrm{and} \ i = j \\ 0 & \mathrm{otherwise} \end{cases}\\ \kappa&= \begin{cases} \kappa^{level} & \text{if } \mathrm{variable} \ i \ \mathrm{is in level} \\ \kappa^{\Delta} & \text{if }\mathrm{variable} \ i \ \mathrm{is in difference} \end{cases}\\ \end{aligned}$$

Here, the autoregressive coefficient $\Pi_{\ell}^{(i,j)}$ is element $\left(i,j\right)$ of $\Pi_{\ell}$ for $\ell=1,\dots,p$. As such, the Minnesota prior sets all prior means for the elements in $\beta$ to $0$, except for the elements that relate to the first own lags of the variables, which are set to $\kappa$. If variable $i$ is in level (e.g. nominal interest rate), then $\kappa=\kappa^{level}$, and typical choices for $\kappa^{level}$ are $1$ or $0.9$. Evaluating the equations at their prior means, equation $i$ becomes a random walk if $\kappa^{level}=1$ and a persistent stationary AR(1) process if $\kappa^{level}=0.9$. Since the steady state only exists if the process is stationary, $0.9$ is recommended for the steady-state BVAR. If variable $i$ is in difference (e.g. output growth), then $\kappa=\kappa^{\Delta}$, and the most common choice for $\kappa^{\Delta}$ is $0$, i.e. equation $i$ becomes (when evaluating it at its prior means) a random walk expressed in first differences. If a differenced variable still shows some degree of persistence (can be examined with an ACF plot), a suitable value for $\kappa^{\Delta}$ can be (for example) $0.6$ instead of $0$. Moving on to the prior variances, $\Omega_\beta$ is a diagonal matrix containing the prior variances for the elements in $\beta$. They are specified as

$$\mathrm{Var}\left(\Pi_{\ell}^{(i,j)}\right)= \begin{cases} \left(\frac{\lambda_1}{\ell^{\lambda_3}}\right)^2 & \text{if } i = j \\ \left(\frac{\lambda_1 \lambda_2\sigma_i}{\ell^{\lambda_3}\sigma_j}\right)^2& \text{if } i \neq j \end{cases}$$

Here $\lambda_1$, $\lambda_2$, and $\lambda_3$ are scalar hyperparameters known as the overall tightness, the cross-equation tightness and the lag decay rate. Furthermore, $\sigma_i^2$ is the $(i,i)$:th element of $\Sigma_u$, which we do not know and therefore replace with an estimate. In this package, it is replaced by the least squares residual variance from a univariate autoregression for variable $i$ with $p$ lags (including the constant and dummy/trend variable if applicable). Moving on to $\Psi$, which contains the steady-state parameters, the prior is

$$\mathrm{vec}(\Psi) \sim \mathrm{N}_{kq}\left[\theta_\Psi,\Omega_\Psi\right]$$

This is the core of the steady-state BVAR model. In $\theta_\Psi$, we specify our prior beliefs about the location of the steady state, and in $\Omega_\Psi$, which we assume to be a diagonal matrix, we specify our degree of certainty in those prior beliefs. Finally, the prior for $\Sigma_u$ is the usual non-informative Jeffreys prior

$$p(\Sigma_u) \propto\left|\Sigma_u \right|^{-(k+1)/2}$$

Alternatively, a proper inverse-Wishart prior can be used (Karlsson, 2013)

$$\Sigma_u \sim \mathrm{IW}(V,m)$$

where $V$ is the scale matrix and $m$ is the number of degrees of freedom. As mentioned, this package also allows for stochastic volatility (Random Walk or AR(1) specifications), where the covariance matrix varies over time, i.e. we have $\Sigma_{u,t}$ (see ?bvar for more details).

Example

To estimate the model in Section 4.1 of Villani (2009), produce (un)conditional forecasts and perform impulse response analysis, simply run the following code:

library(SteadyStateBVAR)
data("Villani2009")
yt <- Villani2009

#hold out last two observations to facilitate comparisons
#to Figures 1-3 in Villani (2009)
yt <- ts(yt[1:102, ], start = start(yt), frequency = frequency(yt))

bvar_obj <- bvar(data = yt)

#Use a dummy to model Sweden’s change in monetary policy in the 1990s
#(move to inflation targeting and flexible exchange rate)
bp <- which(time(yt) == 1992.75) #breakpoint
dummy_variable <- c(rep(1,bp), rep(0,nrow(yt)-bp))

bvar_obj <- setup(bvar_obj,
                  p=4,
                  deterministic = "constant_and_dummy",
                  dummy = dummy_variable)

lambda_1 <- 0.2 #overall tightness
lambda_2 <- 0.5 #cross-equation tightness
lambda_3 <- 1.0 #lag decay rate

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

#95% prior probability intervals (normal distribution)
#See Table I in Villani (2009)
#These are the "steady-state priors"
#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
  )
  )

bvar_obj <- priors(bvar_obj,
                   lambda_1,
                   lambda_2,
                   lambda_3,
                   fol_pm,
                   theta_Psi,
                   Omega_Psi,
                   Jeffreys=TRUE) #FALSE for uninformative inverse-Wishart

p <- bvar_obj$setup$p
k <- bvar_obj$setup$k
kf <- 3 #first three variables in yt are foreign

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
}
print(restriction_matrix)
#block exogeneity for foreign variables
bvar_obj <- restrict_beta(bvar_obj, restriction_matrix)

H <- 12 #forecast horizon
(d_pred <- cbind(rep(1, 12), 0)) #future d_t values


#fit the model
bvar_obj <- fit(bvar_obj,
                H = H,
                d_pred = d_pred,
                iter = 10000,
                warmup = 2500,
                chains = 2,
                cores = 2)

#posterior summaries
summary(bvar_obj , stat = "mean")

#you can look at the stanfit object directly
stan_fit <- bvar_obj$fit$stan
print(stan_fit)

#unconditional forecasts
#see last forecasts in Figures 1-3 in Villani (2009)
fcst <- forecast(bvar_obj,
                 pi = 0.95, #pi = prediction interval
                 fcst_type = "mean",
                 growth_rate_idx = c(4,5), #convert QoQ forecasts to YoY
                 plot_idx = c(4,5,6))

#conditional forecasts
#Toy scenario: inflation gets really high
#What will happen to domestic interest rate?
conditions <- data.frame(
              var        = rep(5,12),
              horizon    = rep(1:12),
              value      = c(1.0,1.5,2.0,1.8,
                             1.5,1.2,1.0,1.0,
                             rep(0.5,4)) #QoQ scale for inflation here
              )
              
cond_fcst <- conditional_forecast(bvar_obj,
                    conditions,
                    pi=0.68,
                    fcst_type = "mean",
                    plot_idx = c(5,6),
                    growth_rate_idx = c(5)) #convert QoQ forecasts to YoY

#impulse response analysis
irf <- IRF(bvar_obj,
           H=20,
           response=5,#inflation
           shock=6, #interest rate
           type="median",
           method="OIRF",
           ci=0.68,
           growth_rate_idx=5) #YoY inflation instead of QoQ

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.

Dieppe, A., Legrand, R., and van Roye, B. (2016). The BEAR toolbox. Working Paper Series, No. 1934. European Central Bank.

Gustafsson, O., and Villani, M. (2025). Variational inference for steady-state BVARs. arXiv preprint arXiv:2506.09271.

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.