
<!-- README.md is generated from README.Rmd. Please edit that file -->

# SteadyStateBVAR

<!-- badges: start -->

[![R-CMD-check](https://github.com/markjwbecker/SteadyStateBVAR/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/markjwbecker/SteadyStateBVAR/actions/workflows/R-CMD-check.yaml)
[![CRAN
status](https://www.r-pkg.org/badges/version/SteadyStateBVAR)](https://CRAN.R-project.org/package=SteadyStateBVAR)
[![CRAN
downloads](https://cranlogs.r-pkg.org/badges/SteadyStateBVAR)](https://cran.r-project.org/package=SteadyStateBVAR)
[![CRAN downloads
total](https://cranlogs.r-pkg.org/badges/grand-total/SteadyStateBVAR)](https://cran.r-project.org/package=SteadyStateBVAR)
<!-- badges: end -->

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:

``` r
install.packages("SteadyStateBVAR")
```

You can also install the development version of SteadyStateBVAR from
[GitHub](https://github.com/) with:

``` r
# 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} (\theta_\beta,\Omega_\beta)$$

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 \ \text{and} \ i = j \\
0 & \text{otherwise}
\end{cases}\\
\kappa&=
\begin{cases}
\kappa^{level} & \text{if} \  \text{variable} \ i \ \text{is in level} \\
\kappa^{\Delta} & \text{if} \ \text{variable} \ i \ \text{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}(\theta_\Psi,\Omega_\Psi)$$

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:

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