Large BVAR with stochastic volatility

Vector autoregressions, Part 9a: where "big data" meets time-varying volatility

Two strands of this series have run separately. Part 3 built a large Bayesian VAR — many variables, disciplined by Minnesota shrinkage — but with a constant shock covariance. Parts 4–5 added stochastic volatility and drifting coefficients, but only on a tiny 3-variable system. Real applications need both at once: a system big enough to be informative and volatilities that move, because macro data went through the Great Moderation and 2008. Carriero-Clark-Marcellino (2019) show how to do it — a large VAR with stochastic volatility, made feasible by an equation-by-equation (triangular) algorithm. This notebook joins the two strands: SV on a 10-variable Bayesian VAR, from scratch in NumPyro, and reads off a by-product the small systems could not give — a common macro-volatility (“uncertainty”) factor.

The model

A 10-variable VAR($p{=}2$) in stationary macro series (annualised growth of output, consumption, investment, industrial production and employment; GDP-deflator and CPI inflation; the unemployment, funds and 10-year rates), with the triangular stochastic-volatility structure of Part 4 carried to all ten equations: $$y_t = c + B_1 y_{t-1} + B_2 y_{t-2} + \varepsilon_t,\qquad A\,\varepsilon_t = u_t,\quad u_{i,t}\sim\mathcal N\!\big(0,\,e^{h_{i,t}}\big),\quad h_{i,t}=h_{i,t-1}+\eta_{i,t}.$$ The constant unit-lower-triangular $A$ orthogonalises the shocks into ten independent components, each with its own random-walk log-volatility $h_{i,t}$; $|A|=1$ so the likelihood is the product of ten time-varying-variance Gaussians (no Jacobian). The VAR coefficients keep the Minnesota prior (bvar.minnesota_moments, $\lambda_1{=}0.2,\lambda_2{=}0.5,\lambda_3{=}1$, own_mean=0 for stationary growth rates); the log-volatilities are written non-centred, $h_{i,\cdot}=h_{i,0}+\sigma_i\operatorname{cumsum}(z_i)$.

Feasibility at scale. The triangular factorisation is exactly what makes a large SV-VAR tractable: conditional on $A$, the $m$ equations' volatilities separate into $m$ univariate stochastic-volatility problems (the Carriero-Clark-Marcellino insight). Here the whole model — coefficients, $A$, and the $10\times T$ volatility surface — is sampled jointly by NUTS in a couple of minutes; two chains give $\hat R \le 1.01$.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS
from numpyro.diagnostics import summary
from scipy.stats import multivariate_normal as mvn
import bvar                                                            # reused for the Minnesota prior (unchanged)

d = pd.read_csv("largesv_data.csv", index_col=0, parse_dates=True)
lab = list(d.columns); Y = d.values; m = Y.shape[1]; p = 2            # m = 10 variables
Yt, X = bvar.build_regressors(Y, p, intercept=True); Teff, k = X.shape
dates = d.index[p:]
li, lj = np.tril_indices(m, -1)
B0, V0, _ = bvar.minnesota_moments(Y, p, lam1=0.2, lam2=0.5, lam3=1.0, own_mean=0.0, intercept=True)
Bsd = np.sqrt(V0)
Bols = np.linalg.lstsq(X, Yt, rcond=None)[0]
h0_loc = np.log(np.diag(np.cov((Yt - X @ Bols).T)))

def model():
    B = numpyro.sample("B", dist.Normal(B0, Bsd))                     # (k,m) Minnesota coefficients
    a = numpyro.sample("a", dist.Normal(0., 1.).expand([len(li)]))    # constant contemporaneous A
    h0 = numpyro.sample("h0", dist.Normal(jnp.asarray(h0_loc), 2.0))
    sigh = numpyro.sample("sigh", dist.HalfNormal(0.3 * jnp.ones(m))) # volatility of volatility
    z = numpyro.sample("z", dist.Normal(0., 1.).expand([Teff, m]))    # non-centred RW innovations
    h = h0[None, :] + sigh[None, :] * jnp.cumsum(z, axis=0)           # 10 log-volatility random walks
    u = (Yt - X @ B) @ jnp.eye(m).at[li, lj].set(a).T
    numpyro.sample("obs", dist.Normal(0., jnp.exp(h / 2)), obs=u)

mc = MCMC(NUTS(model, target_accept_prob=0.9), num_warmup=1200, num_samples=1000,
          num_chains=2, chain_method="sequential", progress_bar=False)
mc.run(jax.random.PRNGKey(1))
ps = mc.get_samples()

# convergence: max split-Rhat over ALL sites, including the full T x m latent volatility surface z
sm = summary(mc.get_samples(group_by_chain=True))
max_rhat = float(max(np.nanmax(v["r_hat"]) for v in sm.values()))
print("Large BVAR-SV (m=%d): 1000 draws x2 | max R-hat over all states %.3f" % (m, max_rhat))

# posterior-median parameters -> reconstruct the ten reduced-form volatilities  Sigma_t = A^-1 diag(e^h) A^-T
Bm = np.median(ps["B"], 0); am = np.median(ps["a"], 0)
h0m = np.median(ps["h0"], 0); sighm = np.median(ps["sigh"], 0); zm = np.median(ps["z"], 0)
hm = h0m[None, :] + sighm[None, :] * np.cumsum(zm, 0)                 # (Teff, m) median log-variance
A = np.eye(m); A[li, lj] = am; Ainv = np.linalg.inv(A)
Sig_t = np.einsum('ij,tj,lj->til', Ainv, np.exp(hm), Ainv)           # reduced-form covariance path
rf_sd = np.sqrt(np.stack([np.diag(Sig_t[t]) for t in range(Teff)]))  # reduced-form stdev (Teff, m)

# in-sample Gaussian log-likelihood, same (median-B) residuals: constant Sigma vs time-varying Sigma_t
Em = Yt - X @ Bm; Sig_c = np.cov(Em.T)
ll_c = mvn(np.zeros(m), Sig_c, allow_singular=True).logpdf(Em).sum()
ll_sv = float(sum(mvn(np.zeros(m), Sig_t[t], allow_singular=True).logpdf(Em[t]) for t in range(Teff)))
print("in-sample log-lik: const %.0f  SV %.0f  (dLL=%.0f)" % (ll_c, ll_sv, ll_sv - ll_c))

# figure: the ten reduced-form volatilities (red) vs the single constant-Sigma estimate (navy dashed)
const_sd = np.sqrt(np.diag(Sig_c))
fig, ax = plt.subplots(2, 5, figsize=(15, 5), sharex=True)
for i, axi in enumerate(ax.ravel()):
    axi.plot(dates, rf_sd[:, i], color="firebrick", lw=1.1)
    axi.axhline(const_sd[i], color="navy", ls="--", lw=1)
    axi.set_title(lab[i], fontsize=9); axi.grid(alpha=.25)
fig.suptitle("Ten reduced-form shock volatilities (red) vs the constant-Σ VAR (navy dashed): "
             "collapse at the Great Moderation, spike around the 2008 crisis", y=1.02, fontsize=11)
fig.tight_layout(); plt.show()

Large BVAR-SV (m=10): 1000 draws x2 | max R-hat over all states 1.006
in-sample log-lik: const -3592  SV -3080  (dLL=513)

Reading the volatilities

Every one of the ten macro shocks has a volatility that collapses at the Great Moderation (~1984) and spikes around the 2008 crisis — most sharply in the real-activity shocks (output, investment, industrial production), more mutedly in inflation and the funds rate, whose own peaks are the 1970s and the Volcker disinflation. A constant-variance VAR (navy dashed) sits uselessly in the middle of each. The in-sample log-likelihood rises by ~510 moving from constant variance to SV — an overwhelming margin for ten volatility processes — so at this scale, as in the 3-variable case, stochastic volatility is not optional.

# The ten volatilities move together: 1st principal component of the standardised log reduced-form volatilities
# -> a common macro-volatility ("uncertainty") factor.
lv = np.log(rf_sd)
lv = (lv - lv.mean(0)) / lv.std(0)                                    # standardise each series
U, S, Vt = np.linalg.svd(lv - lv.mean(0), full_matrices=False)
evr = S**2 / (S**2).sum()
factor = U[:, 0] * S[0]
if np.corrcoef(factor, lv.mean(1))[0, 1] < 0: factor = -factor        # orient to rise with average volatility
print("common volatility factor (PC1 of the ten standardised log-vols) explains %.0f%% of their comovement" % (100 * evr[0]))

# NBER recession quarters, 1960-2019
rec = [("1960-04-01","1961-02-01"),("1969-12-01","1970-11-01"),("1973-11-01","1975-03-01"),
       ("1980-01-01","1980-07-01"),("1981-07-01","1982-11-01"),("1990-07-01","1991-03-01"),
       ("2001-03-01","2001-11-01"),("2007-12-01","2009-06-01")]
fig, ax = plt.subplots(figsize=(11, 3.4))
ax.plot(dates, factor, color="firebrick", lw=1.6)
for a0, a1 in rec: ax.axvspan(pd.Timestamp(a0), pd.Timestamp(a1), color="gray", alpha=.2)
ax.axhline(0, color="k", lw=.5)
ax.set_title("A common macro-volatility factor (PC1, %.0f%% of the comovement) — "
             "high in the 1970s–80s and every recession (shaded), low through the Great Moderation" % (100 * evr[0]))
ax.grid(alpha=.25); fig.tight_layout(); plt.show()

common volatility factor (PC1 of the ten standardised log-vols) explains 73% of their comovement

Results — a common macro-uncertainty factor

• Volatility is largely a single common force. The first principal component of the ten stochastic volatilities explains 73% of their comovement — the macro system does not have ten idiosyncratic volatilities, it has essentially one common volatility factor that rises in the turbulent 1970s–early-1980s and every recession, and falls through the Great Moderation. This is a model-based macroeconomic-uncertainty index (in the spirit of Jurado-Ludvigson-Ng 2015), a by-product only a large SV system can deliver — a 3-variable VAR cannot tell common from idiosyncratic volatility.
• The two strands, joined. This model is Part 3's large Minnesota VAR and Part 4's stochastic volatility in one object: shrinkage makes the many coefficients estimable, the triangular SV makes the time-varying covariance feasible, and together they give both disciplined dynamics and honest, moving uncertainty — the template of modern large-scale Bayesian macro forecasting (Carriero-Clark-Marcellino 2019).
• The R cross-check recovers the same volatility paths with stochvol applied to the VAR residuals.

References

• Carriero, A., Clark, T. E. & Marcellino, M. (2019). Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors. J. Econometrics 212, 137–154.
• Clark, T. E. (2011). Real-time density forecasts from Bayesian VARs with stochastic volatility. J. Business & Economic Statistics 29, 327–341.
• Jurado, K., Ludvigson, S. C. & Ng, S. (2015). Measuring uncertainty. American Economic Review 105, 1177–1216.

Data: largesv_data.csv — 10 stationary FRED-QD macro series, quarterly 1960–2019. Next: largesv_R.ipynb — the R cross-check (stochvol).