This second vignette uses the standard sign-restriction example in structural VARs, following Uhlig (2005). It is the reference exercise for studying the effects of monetary policy shocks on U.S. output, identified by signs only, and it imposes the signs on several horizons. The data are available through the R package bsvarSIGNs, which must be installed: if it is missing, the code chunks below are not run.
The monetary dataset of bsvarSIGNs is sourced from the
replication package of Antolín-Díaz and Rubio-Ramírez (2018). It
comprises six monthly time series spanning January 1965 to November 2007
(Uhlig stops in 2003): real GDP, GDP deflator, a commodity price index,
total reserves and non-borrowed reserves, all in logarithms, and the
federal funds rate in percentage points. The logarithmic values are
multiplied by 100 here, so that the responses read as percentages.
data("monetary", package = "bsvarSIGNs")
Y <- as.matrix(monetary)
Y[, 1:5] <- 100 * Y[, 1:5]
var_names <- c("real GDP", "GDP deflator", "commodity prices",
"total reserves", "non-borrowed reserves", "fed funds")
colnames(Y) <- var_names
dim(Y)
#> [1] 515 6
round(head(Y, 3), 2)
#> real GDP GDP deflator commodity prices total reserves
#> Jan 1965 825.72 292.33 2018.74 248.09
#> Feb 1965 825.87 292.14 2017.96 247.51
#> Mar 1965 826.70 292.25 2015.81 247.52
#> non-borrowed reserves fed funds
#> Jan 1965 245.56 3.90
#> Feb 1965 244.05 3.98
#> Mar 1965 243.99 4.04The model under consideration is a VAR(12), consistent with the approach in Uhlig (2005). The hyperparameters \(\kappa_1\) and \(\kappa_2\) are selected by Empirical Bayes. The total number of draws is fixed at 1000; the chunk prints the selected hyperparameters, the convergence code and the elapsed time.
A single shock is identified (n_id = 1, S
is a row vector). According to Uhlig, after a contractionary monetary
policy shock prices (deflator, commodity prices) and non-borrowed
reserves do not rise, as the Federal Reserve drains liquidity, and the
federal funds rate does not fall. The entry for real GDP is zero,
i.e. unrestricted: the effect on output is the question of the exercise
and is not imposed. Total reserves are also left unrestricted, since
their sign is not determined a priori. The horizon is set to 5, as Uhlig
imposes the sign restrictions on impact and for the following five
months; thus h ranges from 0 to 5. The chunk returns the row vector
S, the indices of the discarded draws (none here) and the
summary of the number of rotations tried per draw.
S <- matrix(c(0, -1, -1, 0, -1, 1), nrow = 1,
dimnames = list("monetary", var_names))
S
#> real GDP GDP deflator commodity prices total reserves
#> monetary 0 -1 -1 0
#> non-borrowed reserves fed funds
#> monetary -1 1
t_sign <- system.time(
accepted <- sign_restrict(mcmc, S, horizons = 5L, max_tries = 20000L, verbose = FALSE)
)
length(accepted)
#> [1] 1000
attr(accepted, "failed_draws")
#> integer(0)
summary(attr(accepted, "n_tries"))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.00 5.00 11.00 15.85 21.00 118.00
t_sign[["elapsed"]]
#> [1] 0.382The maximum impulse response horizon is set to 24 months. The impulse
response functions (IRFs) include six shocks, although only the first
shock is identified, so only its slice is kept. plot_irfs()
draws the response of each variable to the monetary policy shock, with
the posterior median and the 16-84 band.
H <- 24L
irfs <- collect_irfs(accepted, n = 6L, p = p, H = H)
irf_mp <- irfs[, , 1, , drop = FALSE] # the identified shock only
plot_irfs(irf_mp, var_names = var_names, shock_names = "monetary policy",
probs = c(0.16, 0.84))The second chunk provides the output trajectory, including the median and the 16th to 84th percentile equal-tailed credible band at h ∈ {0, 6, 12, 24} months. The findings of Uhlig are corroborated: there is no clear evidence that a monetary tightening reduces real GDP, as the credible bands encompass zero at every horizon and the median is even positive in the first months.
path <- apply(irf_mp[1, , 1, ], 1, stats::quantile, probs = c(0.16, 0.5, 0.84))
colnames(path) <- paste0("h", 0:H)
round(path[, c(1, 7, 13, 25)], 3) # real GDP after a monetary tightening
#> h0 h6 h12 h24
#> 16% -0.045 -0.084 -0.105 -0.154
#> 50% 0.122 0.081 0.048 0.005
#> 84% 0.285 0.230 0.202 0.162This section replicates the same procedure using bsvarSIGNs (Wang and
Woźniak), which is characterised by the Minnesota prior of Giannone,
Lenza and Primiceri with sum-of-coefficients and single-unit-root dummy
observations. In bsvarSIGNs the hyperparameters can be estimated by a
separate adaptive Metropolis step (estimate_hyper()), which
is not called here: they are kept at the package defaults. The sampler
is implemented in C++. The sign_irf argument of bsvarSIGNs
is a three-dimensional array structured as variable by shock by horizon,
where NA indicates an unrestricted response. In acpbvar,
the matrix S is structured as shock by variable, with 0
indicating an unrestricted response. The number of draws remains the
same, and only a single core is used.
sign_irf <- matrix(NA, 6, 6)
sign_irf[, 1] <- c(NA, -1, -1, NA, -1, 1) # variable x shock, horizon 0
sign_irf <- array(sign_irf, dim = c(6, 6, 6)) # horizons 0..5
set.seed(1)
t_bs <- system.time({
spec_bs <- bsvarSIGNs::specify_bsvarSIGN$new(Y, p = p, sign_irf = sign_irf, mc.cores = 1)
post_bs <- bsvars::estimate(spec_bs, S = 1000, show_progress = FALSE)
irf_bs <- bsvars::compute_impulse_responses(post_bs, horizon = H)
})
dim(irf_bs) # variable x shock x horizon x draw
#> [1] 6 6 25 1000
t_bs[["elapsed"]]
#> [1] 3.979The plot presents the median responses to the monetary shock, displaying acpbvar in black and bsvarSIGNs in red. For values of h between 0 and 5, where the signs bind, and in the following months the two medians overlap on all six series; they diverge at longer horizons on real GDP and reserves, where the data are weak and the influence of the prior distributions increases: the dummy observations of the Minnesota prior, at the default hyperparameters, sustain the persistence of the responses, while the asymmetric conjugate prior lets them decay. The differences lie within the credible bands of either package. The table below reports the elapsed seconds on the machine that rendered the vignette.
med_ac <- apply(irf_mp[, , 1, ], c(1, 2), stats::median)
med_bs <- apply(irf_bs[, 1, , ], c(1, 2), stats::median)
op <- par(mfrow = c(2, 3), mar = c(2.5, 2.5, 2, 0.5), mgp = c(1.4, 0.4, 0))
for (j in 1:6) {
ylim <- range(med_ac[j, ], med_bs[j, ])
plot(0:H, med_ac[j, ], type = "l", lwd = 2, ylim = ylim, xlab = "", ylab = "",
main = var_names[j], cex.main = 0.9)
lines(0:H, med_bs[j, ], lwd = 2, lty = 2, col = "firebrick")
abline(h = 0, lty = 3, col = "grey50")
if (j == 1) legend("topright", c("acpbvar", "bsvarSIGNs"), lwd = 2,
lty = c(1, 2), col = c("black", "firebrick"), bty = "n", cex = 0.8)
}The identification is repeated using method = "joint",
where each draw is associated with a single rotation, and the draws
whose rotation does not satisfy the signs are discarded. Since the
per-draw search needed 16 rotations on average, about one draw in
sixteen survives. Achieving the same precision as the
per_draw method requires setting R to several
tens of thousands.