BVAR with stochastic volatility — R cross-check¶

Vector autoregressions, Part 4b¶

Validates bvar_sv_python.ipynb with stochvol (Kastner 2016), the reference R sampler for stochastic volatility. The Python notebook estimates the whole SV-VAR jointly from scratch in NumPyro, with a random-walk log-volatility. Here we take the complementary two-step route: fit the reduced-form VAR by OLS (vars), then run stochvol::svsample on each equation's residuals — a stationary AR(1) log-volatility, an independent spec and an independent implementation. Note this makes the cross-check a deliberate hybrid: a frequentist (OLS) VAR mean paired with a Bayesian stochastic-volatility fit on the residuals — whereas the Python notebook estimates the coefficients and the volatility jointly and fully Bayesianly. (Only svsample is Bayesian here; the VAR coefficients are OLS point estimates, so this two-step path ignores coefficient uncertainty and its implied bands are a touch narrower than a fully joint posterior would give.) If both routes recover the same Great-Moderation collapse and the same crisis spikes, the volatility path is a feature of the data, not of one model or one paradigm.

Data: bvar_sv_data.csv — US quarterly annualised GDP growth, GDP-deflator inflation, and the Fed funds rate, 1959–2019 (FRED-QD).

In [1]:
.libPaths('C:/Users/user/R/win-library/4.6')
suppressMessages({library(vars); library(stochvol)})
d <- read.csv('bvar_sv_data.csv'); dt <- as.Date(d$date); Y <- d[, c('gdp_growth','inflation','tbill')]
v <- VAR(Y, p = 2, type = 'const'); E <- residuals(v)                 # reduced-form residuals
td <- dt[-(1:2)]                                                       # dates aligned to residuals
lab <- c('gdp_growth','inflation','tbill')

set.seed(1); vol <- list()
for (j in 1:3) {
  res <- svsample(E[, j] - mean(E[, j]), draws = 6000, burnin = 2000, quiet = TRUE)
  vol[[j]] <- exp(apply(res$latent[[1]], 2, median) / 2)              # posterior-median SD path, exp(h/2)
}
cat('stochvol AR(1)-SV on VAR residuals (posterior-median volatility):\n')
for (j in 1:3) cat(sprintf('  %-10s peak %.2f (%s) | trough %.2f (%s)  [Python RW-SV: see notebook]\n',
    lab[j], max(vol[[j]]), format(td[which.max(vol[[j]])]), min(vol[[j]]), format(td[which.min(vol[[j]])])))

options(repr.plot.width = 15, repr.plot.height = 3.8); par(mfrow = c(1, 3))
recs <- list(c('1973-11-01','1975-03-01'),c('1980-01-01','1980-07-01'),c('1981-07-01','1982-11-01'),
             c('1990-07-01','1991-03-01'),c('2001-03-01','2001-11-01'),c('2007-12-01','2009-06-01'))
for (j in 1:3) {
  plot(td, vol[[j]], type='l', col='firebrick', lwd=2, xlab='', ylab='shock std. dev.',
       main=paste('stochvol SV:', lab[j]))
  for (r in recs) rect(as.Date(r[1]), -1, as.Date(r[2]), 100, col=rgb(0,0,0,.08), border=NA)
  lines(td, vol[[j]], col='firebrick', lwd=2)
  abline(h = sd(E[, j]), col='navy', lty=2, lwd=1.3)
}
legend('topright', c('SV time-varying','constant sd'), col=c('firebrick','navy'), lwd=2, lty=c(1,2), bty='n', cex=.9)

Warning message:
"package 'vars' was built under R version 4.6.1"

Warning message:
"package 'strucchange' was built under R version 4.6.1"

stochvol AR(1)-SV on VAR residuals (posterior-median volatility):

  gdp_growth peak 5.00 (1978-06-01) | trough 1.36 (2017-09-01)  [Python RW-SV: see notebook]
  inflation  peak 1.65 (1974-09-01) | trough 0.41 (1994-09-01)  [Python RW-SV: see notebook]
  tbill      peak 2.82 (1980-12-01) | trough 0.11 (2018-03-01)  [Python RW-SV: see notebook]
No description has been provided for this image

Same volatility history, independent method¶

  • The Great Moderation is unmistakable in both. stochvol's AR(1) SV on the OLS residuals reproduces the Python random-walk SV almost exactly: every series' shock volatility is high and choppy through the 1970s–early-1980s and drops to a persistently low plateau after ~1984. Funds-rate volatility peaks in the Volcker period (1980), output-growth volatility a little earlier (its peak is 1978, the late-1970s stagflation), inflation volatility peaks with the 1974 oil shock, and the funds-rate volatility collapses toward zero in the post-2008 ZLB — the same landmarks the joint NumPyro model found.
  • Two different SV specifications agree. Python uses a random-walk log-volatility estimated jointly with the VAR; stochvol uses a mean-reverting AR(1) log-volatility on pre-estimated residuals. That they trace the same path is exactly the robustness one wants — the time-variation is in the data, and neither the persistence spec nor the one-step-vs-two-step estimation drives it.
  • Why it matters. A constant-variance VAR (navy dashed) splits the difference and is wrong in both regimes at once: it overstates risk during the Great Moderation and badly understates it in the 1970s and 2008 — which is why adding SV sharply raises the predictive likelihood in the Python notebook.

References¶

  • Kastner, G. (2016). Dealing with stochastic volatility in time series using the R package stochvol. J. Statistical Software 69(5).
  • Cogley, T. & Sargent, T. J. (2005). Drifts and volatilities: monetary policies and outcomes in the post-WWII US. Review of Economic Dynamics 8, 262–302.
  • Carriero, A., Clark, T. E. & Marcellino, M. (2019). Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors. J. Econometrics 212, 137–154.