predictive analyses using R packages bsvars and bsvarSIGNs

\[ \]

structural vector autoregressions

predictive analyses using structural VARs

bsvars and bsvarSIGNs features

a small fiscal policy system analysis

materials

\[ \]

lecture slides as a website

GitHub repo to reproduce the slides and results

\[ \]

bsvars.org official project website

bsvars package on CRAN

bsvarSIGNs package on CRAN

structural vector autoregressions

structural vector autoregressions

• go-to models for the analysis of policy effects

• facilitate the analysis of dynamic causal effects of a well-isolated cause
• relatively simple to work with data and provide empirical evidence on the propagation of shocks through economies and markets
• provide data-driven stylised facts to be incorporated in theoretical models
• extensively used for: monetary and fiscal policy, financial markets, …
• extendible: featuring many variations in specification
  • non-normality
  • heteroskedasticity
  • time-varying parameters
  • Bayesian hierarchical modelling
• Proposed by Sims (1980)

structural vector autoregressions

the model.

\[\begin{align} \text{VAR equation: }&& y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \mathbf{A}_d x_{t} + \epsilon_t\\[1ex] \text{structural equation: }&& \mathbf{B}\epsilon_t &= u_t\\[1ex] \text{structural shocks: }&& u_t |Y_{t-1} &\sim N_N\left(\mathbf{0}_N,\text{diag}\left(\boldsymbol\sigma_t^2\right)\right) \end{align}\]

notation.

• \(y_t\) - \(N\)-vector of variables at time \(t\)
• \(\mathbf{A}_i\) and \(\mathbf{B}\) - \(N\times N\) matrices of autoregressive and structural coefficients
• \(\epsilon_t\) and \(u_t\) - \(N\)-vectors of error terms and structural shocks at time \(t\)
• \(\boldsymbol\sigma_t^2\) - \(N\)-vector structural shock variances

SVARs: hierarchical prior distributions

• normal-generalised normal prior for \(\mathbf{A}\) and \(\mathbf{B}\)
• 3-level local-global equation-specific shrinkage
• Minnesota prior for non-stationary time series
• improves model-fit and forecasting precision

• normal-inverse Wishart prior for \(\mathbf{A}\) and \(\mathbf{\Sigma} = (\mathbf{B}'\mathbf{B})^{-1}\)
• estimated Minnesota prior shrinkage levels
• Minnesota prior for non-stationary time series
• improves model-fit and forecasting precision

SVARs: volatility modelling

• homoskedasticity \(\boldsymbol\sigma_{n.t}^2 = 1\)

• stochastic volatility

• stationary Markov-switching heteroskedasticity

• non-parametric Markov-switching heteroskedasticity

• joint distribution of structural shocks

  • normal
  • finite mixtures of normal distributions
  • non-parametric mixtures of normal distributions
  • Student-t

• homoskedastic
• normal

SVARs: structural identification

• exclusion restrictions
• time-varying volatility
• non-normality

• sign restrictions
• exclusion restrictions
• narrative restrictions

predictive analyses using structural VARs

the objective of economic forecasting

\(\left.\right.\)

… is to use the available data to provide a statistical characterisation of the unknown future values of quantities of interest.

\(\left.\right.\)

The full statistical characterisation of the unknown future values of random variables is given by their predictive density.

\(\left.\right.\)

Simplified outcomes in a form of statistics summarising the predictive densities are usually used in decision-making processes.

\(\left.\right.\)

Summary statistics are also communicated to general audiences.

one-period ahead predictive density

SVAR model.

\[\begin{align*} y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \mathbf{A}_d x_{t} + \epsilon_t\\[2ex] \epsilon_t|Y_{t-1} &\sim iidN_N\left(\mathbf{0}_N,\mathbf\Sigma_t\right)\\[2ex] \mathbf\Sigma_t &= \mathbf{B}^{-1} \text{diag}\left(\boldsymbol\sigma^2_t\right)\mathbf{B}^{-1\prime} \end{align*}\]

one-period ahead conditional predictive density.

… is implied by the model formulation: \[\begin{align*} y_{t+h}|Y_{t+h-1},\mathbf{A},\mathbf\Sigma_{t+h} &\sim N_N\left(\mathbf{A}_1 y_{t+h-1} + \dots + \mathbf{A}_p y_{t+h-p} + \mathbf{A}_d x_{t+h},\mathbf\Sigma_{t+h}\right) \end{align*}\]

one-period ahead predictive density

\(\left.\right.\)

Bayesian forecasting takes into account the uncertainty w.r.t. parameter estimation by integrating it out from the predictive density.

\[\begin{align*} &\\ p(y_{T+1}|Y,X) &= \int p(y_{T+1}|Y_{T},\mathbf{A},\mathbf\Sigma_{T+1}) p(\mathbf{A},\mathbf\Sigma_{T+1}|Y,X) d(\mathbf{A},\mathbf\Sigma_{T+1})\\ & \end{align*}\]

• \(p(y_{T+1}|Y,X)\) - predictive density
• \(p(y_{T+1}|Y_{t},\mathbf{A},\mathbf\Sigma_{T+1})\) - one-period-ahead conditional predictive density
• \(p(\mathbf{A},\mathbf\Sigma_{T+1}|Y,X)\) - marginal posterior distribution
• includes volatility forecasting \(\boldsymbol\sigma^2_{T+1}\) to form \(\mathbf\Sigma_{T+1}\)

sampling from the predictive density

\(\left.\right.\)

step 1: sample from the posterior

… and obtain \(S\) draws \(\left\{ \mathbf{A}^{(s)},\mathbf\Sigma_{T+1}^{(s)} \right\}_{s=1}^{S}\)

\(\left.\right.\)

step 2: sample from the predictive density

In order to obtain draws from \(p(y_{T+1}|Y,X)\), for each of the \(S\) draws of \((\mathbf{A},\mathbf\Sigma_{T+1})\) sample the corresponding draw of \(y_{T+1}\):

Sample \(y_{T+1}^{(s)}\) from \[ N_N\left(\mathbf{A}_1^{(s)} y_{T} + \dots + \mathbf{A}_p^{(s)} y_{T-p+1} + \mathbf{A}_d^{(s)} x_{T+1},\mathbf\Sigma_{T+1}^{(s)}\right) \] and obtain \(\left\{y_{T+1}^{(s)}\right\}_{s=1}^{S}\)

impulse response functions

\(\left.\right.\)

Impulse response functions to orthogonal shocks computed for an empirically relevant SVAR model are considered the dynamic causal effects of the underlying shocks \(u_t\) on economic measurements \(y_{t+i}\) \(i\) periods ahead.

\[\begin{align*} \frac{\partial y_{n.t+i}}{\partial u_{j.t}}&=\theta_{nj.i} \end{align*}\]

• \(\theta_{nj.i}\) - response of \(n\)th variable to \(j\)th shock \(i\) periods after shock’s occurrence

  for \(i=0,1,\dots,h\) and \(n,j=1,\dots,N\)

impulse response functions

\(\left.\right.\)

…as a difference in conditional forecasts.

\[\begin{align} \mathbf\Theta_{h[\cdot i]} &= \mathbb{E}\left[\mathbf{y}_{T+h}| u_{i.T} = 1,\mathbf{Y},\mathbf{X}\right] - \mathbb{E}\left[\mathbf{y}_{T+h}| \mathbf{Y},\mathbf{X}\right] \end{align}\]

…as the shock propagation mechanism.

\[\begin{align} y_t &= \mathbf{B}^{-1}\mathbf{u}_t + \mathbf{A}_1 y_{t-1}\\ &= \mathbf{B}^{-1}\mathbf{u}_t + \mathbf{A}_1 \mathbf{B}^{-1}\mathbf{u}_{t-1} + \mathbf{A}_1^2 y_{t-2}\\ &\dots\\ &= \mathbf{\Theta}_0\mathbf{u}_t + \mathbf{\Theta}_1\mathbf{u}_{t-1} + \mathbf{\Theta}_2\mathbf{u}_{t-2} + \dots \end{align}\]

forecast error variance decompositions

\[ \]

… are a way to quantify the relative importance of the shocks in explaining the forecast error variance of the variables in the system.

\[ \]

forecast error covariance matrix.

\[\begin{align} \mathbf\Sigma_{T+h} &= \mathbf{\Theta}_h \text{diag}\left(\boldsymbol\sigma^2_{T+h}\right)\mathbf{\Theta}_h^{\prime} \end{align}\]

forecast error variance.

For the first variable, at the forecast horizon \(h\), for \(N=2\):

\[\begin{align} \sum_{i=0}^{h-1}\mathbf\Sigma_{T+i[11]} &= \sum_{i=0}^{h-1} \mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i} + \mathbf{\Theta}_{i[12]}^2 \sigma^2_{2.T+i} \end{align}\]

forecast error variance decompositions

\[ \]

… are a way to quantify the relative importance of the shocks in explaining the forecast error variance of the variables in the system.

\[ \]

forecast error variance decomposition.

First shock’s contribution to the forecast error variance of the first variable:

\[\begin{align} FEVD_{T+h[11]} &= \sum_{i=0}^{h-1}\frac{\mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i}}{ \mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i} + \mathbf{\Theta}_{i[12]}^2 \sigma^2_{2.T+i}} \end{align}\]

bsvars and bsvarSIGNs features

bsvars and bsvarSIGNs features

bsvars and bsvarSIGNs features

\[ \]

• Bayesian estimation of Structural VARs
• coherent code structure, workflows, and objects
• excellent computational speed
• frontier econometric & numerical techniques
• written in C++ using Rcpp and RcppArmadillo
• data analysis in R

bsvars and bsvarSIGNs features

• package and data loading

library(bsvars)
data(us_fiscal_lsuw)

• simple model setup

spec = specify_bsvar$new(us_fiscal_lsuw)

• simple estimation

burn = estimate(spec, S = 1000)
post = estimate(burn, S = 10000)

• package and data loading

library(bsvarSIGNs)
data(optimism)

• simple model setup

spec = specify_bsvarSIGN$new(optimism)

• simple estimation

post = estimate(spec, S = 10000)

bsvars and bsvarSIGNs features

• structural analyses

irfs = compute_impulse_responses(post , horizon = 12)
fevd = compute_variance_decompositions(post, horizon = 12)
hds  = compute_historical_decompositions(post)
ss   = compute_structural_shocks(post)
csds = compute_conditional_sd(post)
sddr = verify_identification(post)

• structural analyses

irfs = compute_impulse_responses(post , horizon = 12)
fevd = compute_variance_decompositions(post, horizon = 12)
hds  = compute_historical_decompositions(post)
ss   = compute_structural_shocks(post)
csds = compute_conditional_sd(post)

bsvars and bsvarSIGNs features

• predictive analyses

fvs  = compute_fitted_values(post)
fore = forecast(post, horizon = 12)

• plots and summaries

plot(irfs)
summary(irfs)

• predictive analyses

fvs  = compute_fitted_values(post)
fore = forecast(post, horizon = 12)

• plots and summaries

plot(irfs)
summary(irfs)

bsvars and bsvarSIGNs features

• workflow with the pipe

library(bsvars)
data(us_fiscal_lsuw)

us_fiscal_lsuw |> 
  specify_bsvar$new() |> 
  estimate(S = 1000) |> 
  estimate(S = 10000) -> post

post |> compute_impulse_responses(horizon = 12) |> plot()
post |> compute_variance_decompositions(horizon = 12) |> plot()
post |> compute_historical_decompositions() |> plot()
post |> compute_structural_shocks() |> plot()
post |> compute_conditional_sd() |> plot()
post |> forecast(horizon = 12) |> plot()
post |> verify_identification() |> summary()

• workflow with the pipe

library(bsvarSIGNs)
data(optimism)

optimism |> 
  specify_bsvarSIGN$new() |> 
  estimate(S = 10000) -> post

post |> compute_impulse_responses(horizon = 12) |> plot()
post |> compute_variance_decompositions(horizon = 12) |> plot()
post |> compute_historical_decompositions() |> plot()
post |> compute_structural_shocks() |> plot()
post |> compute_conditional_sd() |> plot()
post |> forecast(horizon = 12) |> plot()

bsvars and bsvarSIGNs features

• progress bar

• progress bar

a small fiscal policy system analysis

a small fiscal policy system analysis

\[\begin{align} y_t = \begin{bmatrix} ttr_t & gs_t & gdp_t \end{bmatrix}' \end{align}\]

• US quarterly data from Q1 1948 to Q2 2024
• all variables are in real, log, per person terms

library(bsvars)
data(us_fiscal_lsuw)    # upload dependent variables
plot(us_fiscal_lsuw, col = "#ff68b4", lwd = 4)

a small fiscal policy system analysis

system with exclusion restrictions.

\[\begin{align} \begin{bmatrix} B_{0.11}&0&0\\ B_{0.21}&B_{0.22}&0\\ B_{0.31}&B_{0.32}&B_{0.33} \end{bmatrix} \begin{bmatrix}ttr_t \\ gs_t \\ gdp_t \end{bmatrix} &= \dots + \begin{bmatrix} u_t^{ttr} \\ u_t^{gs} \\ u_t^{gdp}\end{bmatrix} \end{align}\]

system with sign restrictions.

\[\begin{align} \begin{bmatrix}ttr_t \\ gs_t \\ gdp_t \end{bmatrix} &= \dots + \begin{bmatrix} +&*&+\\ *&+&*\\ *&*&+ \end{bmatrix}\begin{bmatrix} u_t^{ttr} \\ u_t^{gs} \\ u_t^{gdp}\end{bmatrix} \end{align}\]

specify and estimate the SVAR

\[\begin{align} & \end{align}\]

data(us_fiscal_ex)      # upload exogenous variables

# specify the model and set seed
set.seed(123)
spec  = specify_bsvar$new(
  us_fiscal_lsuw, 
  p = 4, 
  exogenous = us_fiscal_ex
)

# run the burn-in
burn        = estimate(spec, S = 1000)

**************************************************|
bsvars: Bayesian Structural Vector Autoregressions|
**************************************************|
 Gibbs sampler for the SVAR model                 |
**************************************************|
 Progress of the MCMC simulation for 1000 draws
    Every draw is saved via MCMC thinning
 Press Esc to interrupt the computations
**************************************************|

# estimate the model
post      = estimate(burn, S = 10000)

**************************************************|
bsvars: Bayesian Structural Vector Autoregressions|
**************************************************|
 Gibbs sampler for the SVAR model                 |
**************************************************|
 Progress of the MCMC simulation for 10000 draws
    Every draw is saved via MCMC thinning
 Press Esc to interrupt the computations
**************************************************|

impulse responses

\[\begin{align} & \end{align}\]

post |> compute_impulse_responses(horizon = 40) |> plot(probability = 0.68, col = "#ff68b4")

forecast error variance decompositions

\[\begin{align} & \end{align}\]

post |> compute_variance_decompositions(horizon = 40) |> plot(col = grad)

conditional forecasts

\[\begin{align} & \end{align}\]

data(us_fiscal_ex_forecasts)      # upload exogenous variables future values
data(us_fiscal_cond_forecasts)    # upload a matrix with projected ttr 
post |> 
  forecast(horizon = 8, exogenous_forecast = us_fiscal_ex_forecasts, conditional_forecast = us_fiscal_cond_forecasts) |> 
  plot(col = "#ff68b4", data_in_plot = 0.2)

specify and estimate the SVAR

\[\begin{align} & \end{align}\]

library(bsvarSIGNs)

# specify the model and set seed
sign = matrix(NA, 3, 3)
diag(sign) = 1
sign[1, 3] = 1
sign

     [,1] [,2] [,3]
[1,]    1   NA    1
[2,]   NA    1   NA
[3,]   NA   NA    1

# specify the model
spec = specify_bsvarSIGN$new(
  us_fiscal_lsuw,
  p = 4,
  sign_irf = sign
)

# estimate the hyper-parameters
spec$prior$estimate_hyper(S = 10000)

**************************************************|
 Adaptive Metropolis MCMC: hyper parameters       |
**************************************************|

# estimate the model
posts = estimate(spec, S = 10000)

**************************************************|
 bsvarSIGNs: Bayesian Structural VAR with sign,   |
             zero and narrative restrictions      |
**************************************************|
 Progress of simulation for 10000 independent draws
 Press Esc to interrupt the computations
**************************************************|

impulse responses

\[\begin{align} & \end{align}\]

posts |> compute_impulse_responses(horizon = 40) |> plot(probability = 0.68, col = "#001D31")

forecast error variance decompositions

\[\begin{align} & \end{align}\]

posts |> compute_variance_decompositions(horizon = 40) |> plot(col = grad)

forecasts

\[\begin{align} & \end{align}\]

posts |>
  forecast(horizon = 8) |>
  plot(col = "#001D31", data_in_plot = 0.2)