Time-varying-parameter VAR with stochastic volatility (Primiceri)

Vector autoregressions, Part 5

Part 4 let the shock variances move but held the VAR coefficients fixed — the economy's dynamics were assumed constant while only its turbulence changed. That is still a strong assumption. Did the way a monetary-policy shock propagates to output and inflation actually stay the same from the Great Inflation, through the Volcker disinflation, into the Great Moderation and the 2008 crisis? Or did the transmission mechanism itself drift — because the Fed's conduct changed, expectations re-anchored, and the structure of the economy evolved?

Primiceri's (2005) time-varying-parameter VAR with stochastic volatility (TVP-VAR-SV) answers this by letting everything drift: the autoregressive coefficients, the contemporaneous relations, and the volatilities all evolve as random walks. It is the canonical tool of empirical macro for studying structural change, and it is the natural capstone of this VAR series — it contains every earlier model as a special case (freeze just the coefficients → Part 4's BVAR-SV; freeze the coefficients and the volatilities → Part 3's constant-$\Sigma$ BVAR). We build it from scratch in NumPyro and read off its central object: the time-varying impulse response to a monetary shock.

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

The model

A VAR whose every parameter is dated. $$y_t = c_t + B_{1,t}\,y_{t-1} + \dots + B_{p,t}\,y_{t-p} + \varepsilon_t,\qquad \varepsilon_t = A_t^{-1}\,\Sigma_t^{1/2}\,z_t,\quad z_t\sim\mathcal N(0,I_m).$$ Collect the coefficients of equation-by-equation into $\theta_t=\operatorname{vec}(c_t,B_{1,t},\dots,B_{p,t})$. The covariance is factored exactly as in Part 4 but now both pieces move: a unit-lower-triangular $A_t$ (the contemporaneous relations) and a diagonal $\Sigma_t=\operatorname{diag}(e^{h_{1,t}},\dots,e^{h_{m,t}})$ (the stochastic volatilities).

Everything follows a driftless random walk (Primiceri's law of motion — parameters wander persistently rather than reverting): $$\theta_t=\theta_{t-1}+u_t,\quad u_t\sim\mathcal N(0,Q);\qquad a_t=a_{t-1}+\zeta_t,\quad \zeta_t\sim\mathcal N(0,S);\qquad h_t=h_{t-1}+\eta_t,\quad \eta_t\sim\mathcal N(0,W),$$ where $a_t$ holds the free below-diagonal elements of $A_t$.

The crucial discipline — tight, calibrated drift. A TVP-VAR has hundreds of drifting states and will happily fit noise as "structural change" if the random walks are allowed to move freely. Primiceri's solution, which we follow, is a pre-sample calibration: run OLS on the first $T_0=40$ quarters to get $\hat\theta_{\text{OLS}}$ and its covariance $V_{\text{OLS}}$, then

• initial states are centred there with a diffuse spread: $\theta_0\sim\mathcal N\!\big(\hat\theta_{\text{OLS}},\,(4\,\operatorname{sd}(\hat\theta_{\text{OLS}}))^2\big)$ (an initial standard deviation four OLS standard errors wide);
• drift sizes are a tiny fraction of the coefficient scale: the per-coefficient random-walk standard deviation is $k_Q\cdot\operatorname{sd}(\hat\theta_{\text{OLS}})$, with the global tightness $k_Q$ given a prior centred at Primiceri's benchmark $0.01$ — i.e. each quarter a coefficient may move only ~1% of its estimation uncertainty. Likewise $k_S$ for $A_t$ and a half-Normal(0.3) volatility-of-volatility for $h_t$ (as in Part 4).

Estimating $k_Q, k_S$ (rather than fixing them) lets the data decide how much time-variation there is, while the tight priors keep it from hallucinating drift.

Algorithm — from scratch in NumPyro. All $T\times(\dim\theta+\dim a+m)$ latent states plus the hyper-parameters are sampled jointly by NUTS. Every random walk is written non-centred, $\theta_t=\theta_0+k_Q\,\operatorname{sd}\odot\operatorname{cumsum}(e_t)$ with $e_t\sim\mathcal N(0,1)$, which defuses the drift-variance funnel that makes TVP models notoriously hard to sample. Because $A_t$ is unit-triangular, $|A_t|=1$ and the likelihood is again just the product of the $m$ time-varying-variance Gaussians in the orthogonalized shocks $u_t=A_t\varepsilon_t$. Two chains; convergence diagnostics for the hyper-parameters and initial states are reported below.

The deliverable: time-varying impulse responses

The point of the model is not the drifting coefficients themselves — individually they are uninterpretable — but the structural object they imply at each date. Freezing the parameters at their date-$t$ values gives a fully specified structural VAR, from which we compute the impulse response to a monetary-policy shock: a recursive (Cholesky) identification with the funds rate ordered last, so the policy shock is the part of the funds-rate innovation orthogonal to contemporaneous output and inflation. We normalise it to a +1 percentage-point funds-rate impact and trace the response of output growth and inflation over 20 quarters — computed separately at four dates (1975, 1981, 1996, 2008) so the change in transmission is visible directly.

A caveat carried over from Part 2: recursive identification here produces a price puzzle (measured inflation ticks up after a tightening) in the early-sample estimates. It is a known artefact of the identification, not of the time-variation; what is informative is how the responses' magnitude drifts across eras, and that the output response — correctly signed — changes markedly.

import numpy as np, pandas as pd
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS
from numpyro.diagnostics import summary
import matplotlib.pyplot as plt
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); Teff, k = X.shape
li, lj = np.tril_indices(m, -1); na = len(li)

# --- pre-sample (first T0) OLS calibration of initial states and drift scales (Primiceri) ---
T0 = 40; X0, Y0 = X[:T0], Yt[:T0]
B_ols = np.linalg.lstsq(X0, Y0, rcond=None)[0]; Sig0 = np.cov((Y0 - X0 @ B_ols).T)
sdB = np.sqrt(np.outer(np.diag(np.linalg.inv(X0.T @ X0)), np.diag(Sig0)))   # OLS coefficient std errors
Ach = np.linalg.cholesky(Sig0); A_ols = (np.linalg.inv(np.diag(np.diag(Ach))) @ Ach)[li, lj]
h0_loc = np.log(np.diag(Sig0))

def model():
    kQ = numpyro.sample("kQ", dist.HalfNormal(0.01))                       # coefficient-drift tightness (Primiceri k_Q)
    kS = numpyro.sample("kS", dist.HalfNormal(0.1))                        # contemporaneous-drift tightness
    B0 = numpyro.sample("B0", dist.Normal(B_ols, 4 * sdB))                 # theta_0 ~ N(OLS, 4 V_OLS)
    a0 = numpyro.sample("a0", dist.Normal(jnp.asarray(A_ols), 1.0))
    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
    eB = numpyro.sample("eB", dist.Normal(0, 1).expand([Teff, k, m]))       # non-centred RW innovations
    ea = numpyro.sample("ea", dist.Normal(0, 1).expand([Teff, na]))
    eh = numpyro.sample("eh", dist.Normal(0, 1).expand([Teff, m]))
    Bt = B0[None] + (kQ * sdB)[None] * jnp.cumsum(eB, 0)                    # drifting coefficients  (Teff,k,m)
    at = a0[None] + kS * jnp.cumsum(ea, 0)                                  # drifting contemporaneous (Teff,na)
    ht = h0[None] + sigh[None] * jnp.cumsum(eh, 0)                          # log-volatility RW (Teff,m)
    eps = Yt - jnp.einsum("tk,tkm->tm", X, Bt)
    At = jnp.broadcast_to(jnp.eye(m), (Teff, m, m)).at[:, li, lj].set(at)
    numpyro.sample("obs", dist.Normal(0., jnp.exp(ht / 2)), obs=jnp.einsum("tij,tj->ti", At, eps))

mc = MCMC(NUTS(model, target_accept_prob=0.9, max_tree_depth=9), num_warmup=1500, num_samples=800,
          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 ("kQ","kS","sigh","B0","a0","h0"))
print("TVP-VAR-SV sampled | max R-hat(hypers/initial states) %.2f | kQ=%.4f  kS=%.3f" %
      (rmax, np.mean(ps["kQ"]), np.mean(ps["kS"])))

# --- time-varying +1pp monetary-shock IRFs at four dates ---
# Reconstruct the drifting states for every posterior draw; freeze the VAR at date t and identify the
# monetary shock recursively (Cholesky, funds rate ordered last): the date-t structural impact matrix is
# P_t = A_t^{-1} diag(e^{h_t/2}) (lower-triangular with positive diagonal), and its last column is the
# response to a policy shock, normalised to a +1pp funds-rate impact. Propagate through the companion form.
off = 1; mp = m * p; H = 21
Bt = np.asarray(ps["B0"][:, None] + (ps["kQ"][:, None, None, None] * sdB[None, None]) * jnp.cumsum(ps["eB"], 1))
at = np.asarray(ps["a0"][:, None] + ps["kS"][:, None, None] * jnp.cumsum(ps["ea"], 1))
ht = np.asarray(ps["h0"][:, None] + ps["sigh"][:, None] * jnp.cumsum(ps["eh"], 1))
S = Bt.shape[0]

def mp_irf(t):                                                             # posterior-median IRF at date index t
    A = np.broadcast_to(np.eye(m), (S, m, m)).copy(); A[:, li, lj] = at[:, t]
    Pt = np.linalg.inv(A) * np.exp(ht[:, t] / 2)[:, None, :]               # A_t^{-1} diag(e^{h_t/2})
    v = Pt[:, :, -1] / Pt[:, -1, -1][:, None]                             # monetary column, +1pp funds impact
    Pi1 = Bt[:, t, off:off + m, :].transpose(0, 2, 1)                      # companion blocks (Pi_l = block^T)
    Pi2 = Bt[:, t, off + m:off + 2 * m, :].transpose(0, 2, 1)
    C = np.zeros((S, mp, mp)); C[:, :m, :m] = Pi1; C[:, :m, m:2 * m] = Pi2; C[:, m:, :m] = np.eye(m)
    st = np.concatenate([v, np.zeros((S, m * (p - 1)))], 1); resp = np.empty((S, H, m))
    for h in range(H):
        resp[:, h] = st[:, :m]; st = np.einsum("sij,sj->si", C, st)
    return np.median(resp, 0)                                             # (H, m)

sel = pd.to_datetime(["1975-03-01", "1981-09-01", "1996-03-01", "2008-12-01"])
tix = [d.index[p:].get_loc(x) for x in sel]
cols = ["#7b1113", "#c0392b", "#e08a2e", "#3b6ea5"]
fig, (axo, axi) = plt.subplots(1, 2, figsize=(13, 4.2))
for t, dt, c in zip(tix, sel, cols):
    irf = mp_irf(t)
    axo.plot(range(H), irf[:, 0], color=c, lw=1.8, label=dt.strftime("%Y"))
    axi.plot(range(H), irf[:, 1], color=c, lw=1.8, label=dt.strftime("%Y"))
for ax, ttl in ((axo, "Output-growth response"), (axi, "Inflation response")):
    ax.axhline(0, color="0.5", lw=.8); ax.grid(alpha=.25)
    ax.set_title(ttl); ax.set_xlabel("quarters after shock"); ax.set_ylabel("pp")
axo.legend(title="date", frameon=False, fontsize=9)
fig.suptitle("Time-varying response to a +1pp monetary-policy shock (posterior median)")
fig.tight_layout(); plt.show()

TVP-VAR-SV sampled | max R-hat(hypers/initial states) 1.00 | kQ=0.0290  kS=0.002

Reading the time-varying transmission

• The real effect of monetary policy was large in the 1970s–80s and much smaller by the Great Moderation. The peak output response to a +1pp tightening (left panel) is −0.72 in 1975 and −0.49 in 1981, versus just −0.15 in 1996 — before deepening again to −0.50 in the 2008 crisis. The monetary transmission to real activity weakened sharply as inflation became anchored and the economy more stable, then re-intensified in the crisis: the central Primiceri/Boivin-Giannoni finding, reproduced from scratch.
• The inflation response drifts too — and shows the price puzzle. Under recursive identification the measured inflation response is positive in the early eras (the price puzzle of Part 2) and shrinks in the later ones; the change across dates is the content, not the level.
• This is structural change a constant-coefficient VAR cannot see. Parts 3–4 would fit a single average response across all four dates; the TVP-VAR-SV shows those dates genuinely differ, which is the whole reason the model exists.

# The stochastic-volatility block still operates underneath the drifting coefficients: the reduced-form shock
# volatility of each variable, sqrt(diag(A_t^-1 diag(e^h_t) A_t^-T)), over time.
# Reconstruct A_t and h_t for every posterior draw, form the reduced-form covariance
# Sigma_rf,t = A_t^{-1} diag(e^{h_t}) A_t^{-T}, and take its per-variable standard deviation over time.
S = ps["h0"].shape[0]
at = ps["a0"][:, None] + ps["kS"][:, None, None] * jnp.cumsum(ps["ea"], 1)          # (S,Teff,na)
ht = ps["h0"][:, None] + ps["sigh"][:, None] * jnp.cumsum(ps["eh"], 1)              # (S,Teff,m)
At = jnp.broadcast_to(jnp.eye(m), (S, Teff, m, m)).at[:, :, li, lj].set(at)         # unit lower-triangular
Ainv = jnp.linalg.inv(At)
Sig_rf = jnp.einsum("stij,stj,stkj->stik", Ainv, jnp.exp(ht), Ainv)                 # (S,Teff,m,m)
vol = np.asarray(jnp.sqrt(jnp.diagonal(Sig_rf, axis1=2, axis2=3)))                  # (S,Teff,m) reduced-form sd
lo, med, hi = np.percentile(vol, [16, 50, 84], axis=0)                              # posterior median + 68% band

dates = d.index[p:]
fig, axes = plt.subplots(1, m, figsize=(15, 3.6))
for j, ax in enumerate(axes):
    ax.plot(dates, med[:, j], color="firebrick", lw=1.6)
    ax.fill_between(dates, lo[:, j], hi[:, j], color="firebrick", alpha=0.15)
    ax.set_title(lab[j]); ax.set_ylabel("reduced-form shock sd"); ax.grid(alpha=.25)
fig.suptitle("Time-varying reduced-form shock volatility (posterior median, 68% band)")
fig.tight_layout(); plt.show()

Results — the whole VAR series, in one model

• Time-variation is real and economically legible. Allowing the dynamics and the volatilities to drift, the model recovers a coherent history: a weakening real transmission of monetary policy across the four eras, alongside the Great-Moderation volatility collapse (second figure) that persists even now that the coefficients are free to move — so the volatility fall was not just an artefact of holding the dynamics fixed in Part 4.
• Estimated, not assumed, time-variation. The data lift the coefficient-drift tightness to $k_Q\approx0.03$ (printed above) — about three times Primiceri's $0.01$ benchmark, but still small in absolute terms: the coefficients move slowly and deliberately, not wildly. Tight priors plus non-centred sampling are what make estimation tractable; the reported $\hat R\approx1.00$ covers the hyper-parameters and initial states (the thousands of latent random-walk innovations are not individually $\hat R$/ESS-checked).
• The capstone of the series. This single model nests the others: freeze $\theta_t, a_t$ and it is Part 4's BVAR-SV; freeze $\theta_t, a_t$ and $h_t$ and it is Part 3's constant-$\Sigma$ BVAR; identify the shocks and it delivers Part 2's structural responses — now allowed to change through time. The R cross-check reaches the same drifting-transmission conclusion from a rolling-window VAR, the frequentist counterpart.
• The frontier from here. The obvious next steps are letting the contemporaneous structure carry the identification more richly (proxy/sign restrictions inside the TVP frame), scaling to a large TVP-VAR with shrinkage, or replacing the random-walk volatilities with the GARCH dynamics of the volatility strand of this project.

References

• Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy. Review of Economic Studies 72, 821–852.
• 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.
• Del Negro, M. & Primiceri, G. E. (2015). Time varying structural VARs with stochastic volatility: a corrigendum. Review of Economic Studies 82, 1342–1345.
• Boivin, J. & Giannoni, M. (2006). Has monetary policy become more effective? Review of Economics and Statistics 88, 445–462.

Next: bvar_tvpsv_R.ipynb — the R cross-check (rolling-window VAR impulse responses).