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 financial 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.
You can install SteadyStateBVAR from CRAN with:
install.packages("SteadyStateBVAR")You can also install the development version of SteadyStateBVAR from GitHub with:
# install.packages("pak")
pak::pak("markjwbecker/SteadyStateBVAR")For a general theoretical introduction to steady-state BVAR models, run
vignette("SteadyStateBVAR-intro")For examples with real macroeconomic data, run
vignette("Homoscedastic-steady-state-BVAR")
vignette("RW-stochastic-volatility-steady-state-BVAR")
vignette("AR1-stochastic-volatility-steady-state-BVAR")Note that vignette("Homoscedastic-steady-state-BVAR")
shows a replication of the empirical example in Section 4.1 of Villani
(2009).
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)
breakpoint <- which(time(yt) == 1992.75)
dummy_variable <- c(rep(1,breakpoint), rep(0,nrow(yt)-breakpoint))
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, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_1: delta y_f
ppi( 1.50, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_1: pi_f
ppi( 4.50, 5.50, interval = 0.95 )$mean, #psi_1: i_f
ppi( 2.00, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_1: delta y
ppi( 1.70, 2.30, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_1: pi
ppi( 4.00, 4.50, interval = 0.95 )$mean, #psi_1: i
ppi( 3.85, 4.00, interval = 0.95 )$mean, #psi_1: q
ppi(-1.00, 1.00, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_2: delta y_f
ppi( 1.50, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_2: pi_f
ppi( 1.50, 2.50, interval = 0.95 )$mean, #psi_2: i_f
ppi(-1.00, 1.00, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_2: delta y
ppi( 4.30, 5.70, interval = 0.95, annualized_growthrate=TRUE, freq=4)$mean, #psi_2: pi
ppi( 3.00, 5.50, interval = 0.95 )$mean, #psi_2: i
ppi(-0.50, 0.50, interval = 0.95 )$mean #psi_2: q
)
Omega_Psi <-
diag(
c(
ppi( 2.00, 3.00, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_1: delta y_f
ppi( 1.50, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_1: pi_f
ppi( 4.50, 5.50, interval = 0.95 )$var, #psi_1: i_f
ppi( 2.00, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_1: delta y
ppi( 1.70, 2.30, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_1: pi
ppi( 4.00, 4.50, interval = 0.95 )$var, #psi_1: i
ppi( 3.85, 4.00, interval = 0.95 )$var, #psi_1: q
ppi(-1.00, 1.00, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_2: delta y_f
ppi( 1.50, 2.50, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_2: pi_f
ppi( 1.50, 2.50, interval = 0.95 )$var, #psi_2: i_f
ppi(-1.00, 1.00, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_2: delta y
ppi( 4.30, 5.70, interval = 0.95, annualized_growthrate=TRUE, freq=4)$var, #psi_2: pi
ppi( 3.00, 5.50, interval = 0.95 )$var, #psi_2: i
ppi(-0.50, 0.50, interval = 0.95 )$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
par(mfrow=c(3,3))
steady_state_priors_plot(bvar_obj, interval = 0.95, growth_rate_idx = c(1,2,4,5))
par(mfrow=c(1,1))
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
}
#block exogeneity for foreign variables
#i.e. foreign variables granger-cause domestic variables but not vice versa
bvar_obj <- restrict_beta(bvar_obj, restriction_matrix)
#fit the model
bvar_obj <- fit(bvar_obj,
H = 12,#forecast horizon
iter = 15000,
warmup = 5000,
chains = 4,
cores = 4)
#elementwise posterior summaries ("mean" or "median")
summary(bvar_obj , stat = "mean")
#you can look at the stanfit object directly
stanfit <- bvar_obj$fit$stan
print(stanfit)
#plot histogram of the posterior steady-state of inflation before and after 1992Q4
rstan::plot(stanfit, pars=c("mu[48,5]", "mu[49,5]"), plotfun="hist")
#note that inflation (pi) is specified as 100*diff(log(CPI)) in the model,
#so for annualized steady-state inflation after 1992Q4 we can do
posterior <- rstan::extract(stanfit)
hist(4*posterior$mu[,49,5], col="darkred", breaks=30)
abline(v=mean(4*posterior$mu[,49,5]), col="lightgreen")
#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", #mean as point forecast
growth_rate_idx = c(4,5), #convert QoQ forecasts to YoY
plot_idx = c(4,5,6),
ss = TRUE, #plot posterior steady-state
ss_type = "mean", #mean as point estimate for posterior steady-state
ss_ci = 0.99, #99% credible interval for steady-state
show_all = FALSE)
#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.95,
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=c(4,5), #gdp growth, inflation
impulse=c(6), #interest rate
type="median",
method="OIRF",
ci=0.95,
growth_rate_idx=c(4,5)) #YoY inflation instead of QoQClark, 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.