Structural Bayesian VAR

After having built the \(Y\) matrix of the Vector Autoregression, we can then build the structural part on the error component, adding sign restrictions as Adolfsen et al. (2024). The reduced-form VAR specification is: \[ Y_t=A_0 + \sum_{l=1}^L A_l Y_{t-l} + B \varepsilon_t ,\] with matrices \(A\) for coefficients and \(B\) to transform the reduced-form residuals into structural shocks. The lags used are 12, with data from January 2010 to December 2022.

shock_names <- c("Supply", "Activity", "Inventory", "Unrestricted")
var_order   <- c("dq", "dp", "dinv", "dip")   # ECB display order: quantity, price, inventories, IP
# This is only for display, Y can be in any order

# Ensure Y columns match the names you expect
stopifnot(all(var_order %in% colnames(Y)))

N <- ncol(Y)

sign_irf <- matrix(NA_real_, nrow = N, ncol = N, dimnames = list(colnames(Y), shock_names))

# restrictions 
# Supply shock: Gas quantity –, Gas price +, Gas inventories –
sign_irf["dq",   "Supply"] <- -1
sign_irf["dp",   "Supply"] <-  1
sign_irf["dinv", "Supply"] <- -1

# Economic activity shock: Gas quantity +, Gas price +, Gas inventories –, IP +
sign_irf["dq",   "Activity"] <-  1
sign_irf["dp",   "Activity"] <-  1
sign_irf["dinv", "Activity"] <- -1
sign_irf["dip",  "Activity"] <-  1

# Inventory shock: Gas quantity +, Gas price +, Gas inventories +
sign_irf["dq",   "Inventory"] <-  1
sign_irf["dp",   "Inventory"] <-  1
sign_irf["dinv", "Inventory"] <-  1

# Unrestricted shock: all NA

# Sign Matrix
sign_irf

##      Supply Activity Inventory Unrestricted
## dq       -1        1         1           NA
## dp        1        1         1           NA
## dinv     -1       -1         1           NA
## dip      NA        1        NA           NA

Here we created sign restrictions on matrix \(B\) that impose the three structural shocks considered:

• a gas supply shock, which reduces the supply of natural gas to the US market, increases the gas price and lowers gas inventories;
• an economic activity shock, which lifts demand for gas due to higher economic production;
• a shock to gas inventories, when gas prices are driven by precautionary demand by gas companies.

A fourth unrestricted shock is also included to capture other dynamics not covered in the previous shocks. This exactly mirrors the approach of Adolfsen et al. (2024), which is derived from the oil modeling literature.

Then we can specify and estimate the model using the package bSVARsigns:

set.seed(123)

spec <- specify_bsvarSIGN$new(
  data      = Y,
  p         = p_lags,
  sign_irf  = sign_irf,
  max_tries = 50000
)

# Posterior draws to save

post <- estimate(spec, S = 20000, thin = 1, show_progress = TRUE)

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

To estimate the model the package implements the standard Minnesota prior with hyperparameter selection as Giannone, Lenza, and Primiceri (2015). This follows the standard approach, which is also followed in the ECB study. The key idea is to treat \(\lambda=(\lambda_1, \lambda_2,\lambda_3)\) (hyperparameters that control tightness, cross-variable shrinkage, and lag-decay) as a random variable with its own prior \(p(\lambda)\) and estimate it jointly (hierarchical structure) from the data. So we can define the joint posterior as: \[p(\beta,\Sigma,\lambda|y) \propto p(y|\beta, \Sigma) \ p(\beta|\Sigma, \lambda) \ p(\Sigma) p(\lambda).\] Intuitively, when the data are informative the model allows more flexibility (less shrinkage, larger variance), while with noisier data it shrinks more aggressively to avoid overfitting.

Impulse Response Functions

After having estimated the model we can then plot Impulse Response Functions to better understand the effects of the shocks imposed. IRFs trace the dynamic effect of a one-time structural shock on the various variables of the VAR over future horizons (20 months in this study). If we define the model looking at its moving-average representation: \[y_t=\mu +\sum_{h=0}^\infty\phi_he_{t-h},\] where \(e_t=B\varepsilon_t\), the structural impulse response coefficients are equal to \(\Phi_h=\phi_hB\).

irf <- compute_impulse_responses(post, horizon = 20, standardise = FALSE)
IR  <- unclass(irf)
IR <- 100 * IR #to make them in percent as ECB

We can plot the IRFs functions with the built in function:

plot(irf*100, probability = 0.68)

We can also plot the results trying to replicate the same plot structure in the paper (shocks on rows, variables on columns):

H_plot <- 17 #17 for the graph, 20 for the IRFs
hgrid  <- 0:H_plot

plot_vars   <- c("dq","dp","dinv","dip")
plot_shocks <- c("Supply","Activity","Inventory")  # no unrestricted
col_titles  <- c("Gas quantity","Gas price","Gas inventories","Industrial production")

var_idx   <- match(plot_vars,   colnames(Y))
shock_idx <- match(plot_shocks, colnames(sign_irf))

# Standard quantiles 
qfun <- function(x) quantile(x, probs = c(0.05, 0.16, 0.50, 0.84, 0.95), na.rm = TRUE, type = 8)

op <- par(no.readonly = TRUE)
par(mfrow = c(length(plot_shocks), length(plot_vars)),
    mar = c(2.8, 2.8, 2.2, 1.0), oma = c(2, 2, 2, 0))

for (s in seq_along(plot_shocks)) {
  for (v in seq_along(plot_vars)) {

    i <- var_idx[v]
    j <- shock_idx[s]

    qs <- sapply(seq_along(hgrid), function(hh) qfun(IR[i, j, hh, ]))
    # qs is 5 x (H+1)
    lo90 <- qs[1,]; lo68 <- qs[2,]; med <- qs[3,]; hi68 <- qs[4,]; hi90 <- qs[5,]

    ylim <- range(c(lo90, hi90, 0), na.rm = TRUE)

    plot(hgrid, med, type = "n", xlab = "", ylab = "", ylim = ylim)
    abline(h = 0, lty = 3)

    polygon(c(hgrid, rev(hgrid)), c(lo90, rev(hi90)), #for the shaded parts
            border = NA, col = rgb(0,0,0,0.10))
    polygon(c(hgrid, rev(hgrid)), c(lo68, rev(hi68)),
            border = NA, col = rgb(0,0,0,0.18))

    lines(hgrid, med, lwd = 2)

    if (s == 1) title(main = col_titles[v], cex.main = 1.0)
    if (v == 1) mtext(plot_shocks[s], side = 2, line = 2.2, cex = 0.9)
  }
}

mtext("Months", side = 1, outer = TRUE, line = 0.8)
mtext("Impulse Responses (percent; median, 68% & 90% bands)", side = 3, outer = TRUE, line = 0.5)

These results can be compared to the results obtained by the ECB on the similar timeframe analyzed.