BVAR with stochastic volatility — letting the shock variances move

Vector autoregressions, Part 4

Every VAR so far — through the large 20-variable system of Part 3 — assumed a constant shock covariance $\Sigma$: the economy is taken to be equally turbulent in 1974, in 1998, and in 2008. That is plainly false. US macro data went through a Great Moderation — a sharp, lasting fall in the volatility of output, inflation and interest rates around 1984 — bracketed by the turbulent 1970s and the 2008 crisis. A fixed-$\Sigma$ VAR splits the difference: its forecast intervals are simultaneously too wide in calm times and too narrow in storms, so its risk assessments are wrong exactly when they matter.

Stochastic volatility (SV) fixes this by letting the shock variances follow their own stochastic process. It is the single most valuable extension to a forecasting VAR: the literature (Clark 2011; D'Agostino-Gambetti-Giannone 2013) finds that time-varying volatility delivers most of the gain from "time-varying-parameter" models — more than letting the regression coefficients drift. This notebook adds SV to the monetary BVAR from scratch in NumPyro, keeping the coefficients constant so the volatility is the whole story; the next part lets the coefficients move too.

The model

Reduced form with a time-varying covariance. A VAR($p$) whose innovation covariance is now dated, $$y_t = c + B_1 y_{t-1} + \dots + B_p y_{t-p} + \varepsilon_t,\qquad \varepsilon_t \sim \mathcal N\!\big(0,\ \Sigma_t\big),$$ with the coefficients $c, B_1,\dots,B_p$ held constant and all the time-variation loaded into $\Sigma_t$.

Triangular stochastic-volatility factorization (Cogley-Sargent 2005; Primiceri 2005). Write $\Sigma_t$ through a constant unit-lower-triangular matrix $A$ and a time-varying diagonal: $$A\,\varepsilon_t = u_t,\qquad u_{i,t}\sim\mathcal N\!\big(0,\ e^{h_{i,t}}\big)\ \text{independent},\qquad\Longrightarrow\qquad \Sigma_t = A^{-1}\operatorname{diag}\!\big(e^{h_{1,t}},\dots,e^{h_{m,t}}\big)A^{-\top}.$$ $A$ orthogonalizes the shocks into $m$ independent structural components $u_{i,t}$; each component carries its own log-variance $h_{i,t}$ that drifts over time. Because $A$ has ones on the diagonal, $|A|=1$: the change of variables from $\varepsilon_t$ to $u_t$ has unit Jacobian, so the likelihood is simply the product of the $m$ univariate Gaussians in $u_{i,t}$ — no determinant term.

Log-volatilities as random walks. Each log-variance evolves as a driftless random walk, $$h_{i,t} = h_{i,t-1} + \eta_{i,t},\qquad \eta_{i,t}\sim\mathcal N(0,\sigma_i^2),$$ the most common SV law of motion: volatility is highly persistent and free to wander to a new level (as it did, permanently, at the Great Moderation) rather than reverting to a fixed mean.

Priors — every hyperparameter stated.

• Coefficients $B$: the same Minnesota prior as the rest of this series (via bvar.py), with $\lambda_1=0.2$ (overall tightness), $\lambda_2=0.5$ (cross-variable shrinkage), $\lambda_3=1$ (lag decay), and scale factors $s_i/s_j$. Here own_mean = 0: the data are growth rates, stationary, so the prior shrinks toward zero (white noise), not toward a random walk.
• Contemporaneous matrix $A$: each free below-diagonal element $\sim\mathcal N(0,1)$.
• Initial log-volatility $h_{i,0}\sim\mathcal N\!\big(\log\hat\sigma^2_{i,\text{OLS}},\ 2^2\big)$ — centred at the constant-variance estimate.
• Volatility of volatility $\sigma_i\sim\text{Half-Normal}(0.3)$ — how fast each log-variance can move.

Algorithm — from scratch in NumPyro. The whole model (VAR coefficients, $A$, and the $T\times m$ log-volatility surface) is sampled jointly by the No-U-Turn Sampler. The random walk is written in non-centred form $h_{i,\cdot}=h_{i,0}+\sigma_i\,\operatorname{cumsum}(z_i)$ with $z_{i,t}\sim\mathcal N(0,1)$, which removes the volatility-of-volatility funnel that otherwise cripples HMC on SV models. Two chains give $\hat R = 1.001$ and minimum effective sample size ~890 — clean convergence.

The data

A compact monetary VAR at quarterly frequency, 1959–2019 (pre-COVID), from the same FRED-QD source as Part 3, transformed to stationarity:

• gdp_growth — annualised real GDP growth ($400\,\Delta\log$),
• inflation — annualised GDP-deflator inflation ($400\,\Delta\log$),
• tbill — the Fed funds rate (level, %).

Three variables keep the $T\times 3$ volatility surface easy to read; the SV machinery is identical for any $m$. We use $p=2$ lags.

import numpy as np, pandas as pd
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
import matplotlib.pyplot as plt
from numpyro.infer import MCMC, NUTS
from numpyro.diagnostics import summary
import bvar

d = pd.read_csv("bvar_sv_data.csv", index_col=0, parse_dates=True)
lab = list(d.columns); Y = d.values.astype(float); m = Y.shape[1]; p = 2
Yt, X = bvar.build_regressors(Y, p, intercept=True)                      # (T,m), (T,k)
Teff, k = X.shape

# Minnesota prior on B (own_mean=0: stationary growth-rate VAR) -- reuse bvar.py, unchanged
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]; E = Yt - X @ Bols          # OLS residuals ->
Sig_const = np.cov(E.T); h0_loc = np.log(np.diag(Sig_const))            # constant-vol baseline + h0 centre
li, lj = np.tril_indices(m, -1)                                          # A's free (below-diagonal) elements

def model():
    B = numpyro.sample("B", dist.Normal(B0, Bsd))                        # (k,m) Minnesota coefficients
    eps = Yt - X @ B                                                     # reduced-form residuals
    a = numpyro.sample("a", dist.Normal(0., 1.).expand([len(li)]))       # contemporaneous A (unit lower-tri)
    A = jnp.eye(m).at[li, lj].set(a)
    u = eps @ A.T                                                        # orthogonal structural shocks
    h0 = numpyro.sample("h0", dist.Normal(jnp.asarray(h0_loc), 2.0))     # initial log-variance
    sig = numpyro.sample("sig", 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, :] + sig[None, :] * jnp.cumsum(z, axis=0)               # log-volatility random walk
    numpyro.sample("obs", dist.Normal(0., jnp.exp(h / 2)), obs=u)        # |A|=1 => no Jacobian term

mc = MCMC(NUTS(model, target_accept_prob=0.9), num_warmup=1500, num_samples=1500,
          num_chains=2, chain_method="sequential", progress_bar=False)
mc.run(jax.random.PRNGKey(1)); ps = mc.get_samples()
sm = summary(mc.get_samples(group_by_chain=True), prob=0.9)
rmax = max(float(np.nanmax(sm[q]["r_hat"])) for q in ("a","sig","h0"))
emin = min(float(np.nanmin(sm[q]["n_eff"])) for q in ("a","sig","h0"))
print("sampled 1500 draws x2 chains | max R-hat(a,sig,h0) %.3f | min ESS %.0f" % (rmax, emin))

# --- reconstruct the reduced-form time-varying volatility surface from the posterior draws ---
# reduced-form sd of variable i at t = sqrt(diag(A^-1 diag(e^h_t) A^-T))_i, evaluated per draw then summarised
idx = d.index[p:]                                                        # dates aligned to Yt (T-p rows)
a_dr = np.asarray(ps["a"]); ndr = a_dr.shape[0]
h_dr = (np.asarray(ps["h0"])[:, None, :]
        + np.asarray(ps["sig"])[:, None, :] * np.cumsum(np.asarray(ps["z"]), axis=1))   # (ndraw,Teff,m) log-var
sd_dr = np.empty((ndr, Teff, m))                                         # reduced-form sd path, per draw
Sig_bar = np.zeros((Teff, m, m))                                         # posterior-mean time-varying covariance
for j in range(ndr):
    Aj = np.eye(m); Aj[li, lj] = a_dr[j]; Ainv = np.linalg.inv(Aj)       # Sigma_t = Ainv diag(e^h) Ainv^T
    ev = np.exp(h_dr[j])
    sd_dr[j] = np.sqrt((Ainv ** 2) @ ev.T).T                            # sqrt of its diagonal
    Sig_bar += np.einsum("ik,tk,jk->tij", Ainv, ev, Ainv)
Sig_bar /= ndr
sd_sv = np.median(sd_dr, axis=0)                                         # posterior-median volatility path
sd_lo, sd_hi = np.percentile(sd_dr, [5, 95], axis=0)                     # 90% credible band
sd_const = np.sqrt(np.diag(Sig_const))                                   # single constant-vol baseline

# --- in-sample fit: constant-Sigma vs time-varying Sigma_t (reduced-form Gaussian log-likelihood) ---
eps_pm = Yt - X @ np.asarray(ps["B"]).mean(0)                            # posterior-mean reduced-form shocks
Scinv = np.linalg.inv(Sig_const); _, ldc = np.linalg.slogdet(Sig_const)
ll_const = float(np.sum(-0.5 * m * np.log(2 * np.pi) - 0.5 * ldc
                        - 0.5 * np.einsum("ti,ij,tj->t", E, Scinv, E)))
ll_sv = 0.0
for t in range(Teff):
    S = Sig_bar[t]; Sinv = np.linalg.inv(S); _, ld = np.linalg.slogdet(S)
    ll_sv += -0.5 * m * np.log(2 * np.pi) - 0.5 * ld - 0.5 * eps_pm[t] @ Sinv @ eps_pm[t]
ll_sv = float(ll_sv)
print("in-sample total log-lik:  constant-vol %d   SV %d   (Delta = %d in favour of SV)"
      % (round(ll_const), round(ll_sv), round(ll_sv - ll_const)))
for i in range(m):
    pk = int(sd_sv[:, i].argmax()); tr = int(sd_sv[:, i].argmin())
    print("  %-10s const sd %.2f | SV peak %.2f (%s) | SV trough %.2f (%s)"
          % (lab[i], sd_const[i], sd_sv[pk, i], idx[pk].date(), sd_sv[tr, i], idx[tr].date()))

# --- figure: reduced-form shock volatility, time-varying (SV) vs constant, one panel per variable ---
fig, ax = plt.subplots(m, 1, figsize=(9, 7), sharex=True)
for i in range(m):
    ax[i].fill_between(idx, sd_lo[:, i], sd_hi[:, i], color="firebrick", alpha=.20, lw=0,
                       label="SV 90% band" if i == 0 else None)
    ax[i].plot(idx, sd_sv[:, i], color="firebrick", lw=1.4,
               label="SV posterior median sd" if i == 0 else None)
    ax[i].axhline(sd_const[i], color="navy", ls="--", lw=1.2,
                  label="constant sd" if i == 0 else None)
    ax[i].set_ylabel(lab[i]); ax[i].grid(alpha=.25)
ax[0].legend(fontsize=8, loc="upper right")
ax[0].set_title("Reduced-form shock volatility: time-varying (SV) vs constant")
fig.tight_layout(); plt.show()

sampled 1500 draws x2 chains | max R-hat(a,sig,h0) 1.001 | min ESS 889

in-sample total log-lik:  constant-vol -1211   SV -914   (Delta = 297 in favour of SV)
  gdp_growth const sd 3.01 | SV peak 4.97 (1978-06-01) | SV trough 1.20 (2018-03-01)
  inflation  const sd 0.94 | SV peak 1.96 (1974-09-01) | SV trough 0.33 (1995-06-01)
  tbill      const sd 0.80 | SV peak 3.17 (1980-12-01) | SV trough 0.05 (2012-06-01)

Reading the volatility paths

Each panel is the posterior time-varying standard deviation of one variable's shock (red band), against the single constant estimate a fixed-$\Sigma$ VAR would use (navy dashed). The history is written plainly in them:

• The Great Moderation. All three volatilities are high and choppy through the 1970s and early 1980s, then drop to a low, stable plateau after ~1984 and stay there. Output-growth shock volatility falls from ~5 to barely above 1; the constant-$\Sigma$ line sits uselessly in between.
• Named episodes. Interest-rate volatility peaks in the 1979–82 Volcker disinflation, and output-growth volatility a little earlier — its posterior peak is in 1978, in the turbulent late-1970s stagflation; inflation volatility peaks with the 1974 oil shock; every series flares in 2008; and the funds-rate volatility collapses toward zero after 2008, when the policy rate was pinned at the zero lower bound and simply stopped moving.

# Diagnostic: divide each structural shock by its estimated volatility path. If the SV model is right,
# the standardized shocks u_it / exp(h_it/2) should be ~ iid N(0,1) -- a flat rolling std around 1 --
# whereas standardizing by a single constant sd leaves the 1970s-vs-Great-Moderation heteroskedasticity.

# structural shocks u = A eps and their SV volatility path exp(h/2), at the posterior median
a_med = np.median(np.asarray(ps["a"]), axis=0)
h_med = np.median(h_dr, axis=0)                                          # (Teff,m) posterior-median log-var
A_med = np.eye(m); A_med[li, lj] = a_med
u_pm = eps_pm @ A_med.T                                                  # posterior-mean structural shocks
z_sv = u_pm / np.exp(h_med / 2)                                          # standardized by the SV path
z_const = u_pm / u_pm.std(axis=0, ddof=1)                               # standardized by a single constant sd

W = 20                                                                   # 5-year (20-quarter) rolling window
roll = lambda x: pd.Series(x, index=idx).rolling(W, center=True, min_periods=W // 2).std()

# --- figure: rolling std of the standardized structural shocks -- SV whitens to ~1, constant does not ---
fig, ax = plt.subplots(m, 1, figsize=(9, 7), sharex=True)
for i in range(m):
    ax[i].plot(idx, roll(z_const[:, i]), color="navy", lw=1.3,
               label="standardized by constant sd" if i == 0 else None)
    ax[i].plot(idx, roll(z_sv[:, i]), color="firebrick", lw=1.4,
               label="standardized by SV path" if i == 0 else None)
    ax[i].axhline(1, color="k", lw=.5, ls=":")
    ax[i].set_ylabel(lab[i]); ax[i].grid(alpha=.25)
ax[0].legend(fontsize=8, loc="upper right")
ax[0].set_title("Whitening check: rolling std of standardized structural shocks (5-yr window)")
fig.tight_layout(); plt.show()

Results — volatility is not a nuisance, it is the signal

• Stochastic volatility massively improves fit. The in-sample log-likelihood rises by ~300 moving from constant variance to SV (printed above) — an overwhelming margin for three added volatility processes. The data are emphatic that $\Sigma$ is not constant.
• The standardization test passes. Dividing each shock by its SV volatility path whitens it to a flat rolling standard deviation near 1 (red, second figure); dividing by a constant standard deviation leaves the volatility swings glaring (navy). The SV path is capturing real, systematic heteroskedasticity, not overfitting noise.
• Calibrated forecast uncertainty. This is the practical payoff. A constant-$\Sigma$ VAR reports the same predictive interval in 1998 as in 2008; the SV-VAR's intervals widen in the turbulent 1970s and 2008 and narrow through the Great Moderation, so its probability statements are honest in every regime — the property that matters for risk management, fan charts and density forecasting.
• Coefficients held fixed. All of this came from time-varying volatility alone — the VAR coefficients never moved. That is the deliberate design: SV is where the forecasting gains concentrate, so it earns its own notebook before we let the dynamics drift in the time-varying-parameter VAR to come.

The R cross-check reaches the same volatility history from the opposite direction — a stationary AR(1) SV fit to the OLS residuals with stochvol — confirming the path is in the data, not the prior.

References

• 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.
• Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy. Review of Economic Studies 72, 821–852.
• Clark, T. E. (2011). Real-time density forecasts from Bayesian VARs with stochastic volatility. J. Business & Economic Statistics 29, 327–341.
• Kastner, G. (2016). Dealing with stochastic volatility in time series using stochvol. J. Statistical Software 69(5).

Next: bvar_sv_R.ipynb — the R cross-check (stochvol on the VAR residuals).