Bayesian VAR — macro application: the monetary-policy VAR

Vector autoregressions, part 1: classical VAR, Minnesota BVAR, and the Villani steady-state prior

This applies a from-scratch Bayesian VAR engine (bvar.py) to a small US monetary-policy VAR. The engine is a from-scratch Gibbs sampler for the Minnesota-prior VAR and the Villani mean-adjusted ("steady-state") VAR. We pair it with the classical frequentist VAR (statsmodels) for the identification/impulse-response machinery, and cross-check the Bayesian fit in R (bvar_macro_R.ipynb).

The system. Three quarterly US series, 1959Q2–2009Q3 (statsmodels.datasets.macrodata), the canonical small monetary VAR:

• GDP growth — annualised real output growth, $400\,\Delta\log(\text{realgdp})$
• Inflation — annualised CPI inflation, $400\,\Delta\log(\text{cpi})$
• T-bill — the 3-month Treasury-bill rate (the policy-rate proxy)

The three questions a VAR answers: how does a monetary shock propagate (impulse responses), how much of each variable's fluctuation does each shock explain (variance decomposition), and — the Bayesian payoff — how do we forecast while imposing sensible long-run beliefs (the steady-state prior)?

Model 1 — the VAR($p$) and its structural interpretation

The reduced-form VAR. For an $m$-vector $y_t$ (here $m=3$), a VAR of order $p$ is $$y_t = c + \Pi_1 y_{t-1} + \dots + \Pi_p y_{t-p} + \varepsilon_t,\qquad \varepsilon_t\sim N(0,\Sigma),$$ $m$ equations, each regressing a variable on $p$ lags of every variable. Stacked over $t$ it is a multivariate regression $Y = X B + E$ with $k=1+mp$ regressors per equation. Stationarity is governed by the companion matrix $F$ (stacking $\Pi_1,\dots,\Pi_p$): the process is covariance-stationary iff all eigenvalues of $F$ lie inside the unit circle.

Estimation (classical). Because every equation has the same regressors $X$, equation-by-equation OLS is the efficient (and ML, under normality) estimator: $\hat B=(X'X)^{-1}X'Y$, $\hat\Sigma=\hat E'\hat E/T$.

Granger causality. Variable $j$ Granger-causes $i$ if lags of $j$ help predict $i$ given everything else — an $F$-test that the relevant $\Pi$ blocks are zero.

Structural identification. The reduced-form innovations $\varepsilon_t$ are correlated, so they are not economically interpretable. We recover structural shocks $u_t$ with $\varepsilon_t = P u_t$, $u_t\sim N(0,I)$, via the Cholesky factor $\Sigma = PP'$ (lower-triangular $P$) — a recursive identification. The ordering [GDP, Inflation, T-bill] encodes the standard timing assumption: the policy rate reacts to output and inflation within the quarter, but affects them only with a lag (so the monetary shock is ordered last).

Impulse responses & FEVD. Inverting to the MA($\infty$) form $y_t=\sum_h\Theta_h\varepsilon_{t-h}$, the orthogonalized IRF $\Theta_h P$ traces each variable's response to a one-standard-deviation structural shock, and the forecast-error variance decomposition (FEVD) apportions each variable's $h$-step forecast variance across the structural shocks.

Model 2 — the Minnesota-prior Bayesian VAR

An unrestricted VAR has $m(1+mp)$ coefficients — here $3\cdot7=21$ at $p=2$, and it grows fast with lags and variables, so OLS overfits. The Minnesota prior (Litterman 1986) is a shrinkage device: it pulls each equation toward a parsimonious benchmark and lets the data override only where the evidence is strong. We use the independent Normal–inverse-Wishart version (a diagonal prior on the coefficients, so asymmetric shrinkage is possible), $\operatorname{vec}(B)\sim N(b_0,V_0)$, $\Sigma\sim IW(v_0,S_0)$.

Prior mean $b_0$. Every coefficient is centred at $0$ except each variable's own first lag, centred at own_mean — set to 0 here because the series are already differenced (growth rates), i.e. the prior benchmark is white noise around a mean, not a random walk.

Prior variance (the shrinkage, Meseguer eq. 16). For the coefficient on variable $j$ at lag $l$ in equation $i$, the prior standard deviation is $$\sqrt{V_0}=\begin{cases}\dfrac{\lambda_1}{l^{\lambda_3}} & j=i\ \text{(own lags)}\\[2mm]\dfrac{\lambda_1\,\lambda_2(i,j)}{l^{\lambda_3}}\dfrac{s_i}{s_j} & j\ne i\ \text{(cross lags)}\end{cases}$$ governed by three Minnesota hyperparameters $(\lambda_1,\lambda_2,\lambda_3)$ plus a scale factor:

• $\lambda_1$ — overall tightness. The own-first-lag prior std; the master knob. As $\lambda_1\to0$ the prior forces every equation to its benchmark (here white noise), while $\lambda_1\to\infty$ gives a flat prior, i.e. OLS. We use $\lambda_1=0.2$.
• $\lambda_2$ — cross-variable (relative) tightness, $0\le\lambda_2\le1$. It scales the shrinkage on other variables' lags relative to a variable's own lags: $\lambda_2=1$ trusts them equally, $\lambda_2<1$ shrinks cross-lags harder (the Minnesota belief that a series is explained mostly by its own past). We set it per pair, $\lambda_2(i,j)=0.8\,\widehat{\text{corr}}(y_i,y_j)$, so cross terms are loosened only between series that genuinely co-move.
• $\lambda_3$ — lag-decay rate. The exponent in $l^{\lambda_3}$: higher lags are shrunk toward zero ever harder, encoding that distant lags carry less information. We use $\lambda_3=1$ (linear decay).
• $s_i/s_j$ — scale correction. The ratio of univariate AR($p$) residual standard deviations, making the shrinkage unit-free across variables of very different volatility.

The intercept is left diffuse (a huge prior variance, essentially unrestricted).

Estimation — 2-block Gibbs sampler (gibbs_minnesota). The two full conditionals are conjugate, so we alternate:

1. Coefficients $\beta=\operatorname{vec}(B)\mid\Sigma,Y \sim N(\bar b,\bar V)$ with $$\bar V^{-1}=V_0^{-1}+(\Sigma^{-1}\otimes X'X),\qquad \bar b=\bar V\big(V_0^{-1}b_0+(\Sigma^{-1}\otimes X'X)\,\hat\beta\big),$$ drawn from its precision via a Cholesky solve (this is the SUR/GLS update, tying the equations together through $\Sigma^{-1}$).
2. Covariance $\Sigma\mid\beta,Y \sim IW\big(v_0+T,\ S_0+E'E\big)$ from the current residuals $E=Y-XB$.

Iterate, discard burn-in, keep the draws — a full posterior over $(B,\Sigma)$, hence over impulse responses, forecasts, and any function of them.

Model 3 — the Villani steady-state (mean-adjusted) BVAR

A reduced-form VAR's long-run mean is the tangled function $(\,I-\sum_l\Pi_l)^{-1}c$ — you cannot state a belief about it directly. Villani (2009) re-parameterises the VAR around its unconditional mean $\psi$: $$y_t = \psi + \sum_{l=1}^p \Pi_l\,(y_{t-l}-\psi) + \varepsilon_t,$$ so $\psi$ is an explicit parameter you can put an informative prior on. This is the model's headline feature for macro: you can pin the long-run steady state — e.g. a $2\%$ inflation target and trend growth — as $\psi\sim N(\psi_0,\Lambda_0)$, and the forecasts will revert to that, not to whatever the (crisis-laden) sample average happens to be.

Estimation — 3-block Gibbs sampler (gibbs_steady_state). Conditioning on $\psi$ makes the dynamics a Minnesota-prior regression on demeaned data; conditioning on $\Pi$ makes $\psi$ a linear-Gaussian problem. We cycle:

1. Dynamics $\Pi\mid\psi,\Sigma$: regress the demeaned $z_t=y_t-\psi$ on its lags with the same Minnesota prior $\Rightarrow$ Normal (as in Model 2, on $z$).
2. Steady state $\psi\mid\Pi,\Sigma$: the residual identity $u_t=\big(I-\sum_l\Pi_l\big)\psi+\varepsilon_t$ gives, with $D=I-\sum_l\Pi_l$, $$\psi\mid\cdot \sim N\big(\bar\psi,\bar\Lambda\big),\quad \bar\Lambda^{-1}=\Lambda_0^{-1}+T\,D'\Sigma^{-1}D,\quad \bar\psi=\bar\Lambda\big(\Lambda_0^{-1}\psi_0+D'\Sigma^{-1}\textstyle\sum_t u_t\big).$$
3. Covariance $\Sigma\mid\Pi,\psi \sim IW\big(v_0+T,\ S_0+E'E\big)$.

Forecasting. For each retained posterior draw we simulate $H$-step paths forward (Model 1's recursion with the drawn $\Pi/B$ and a fresh $\varepsilon\sim N(0,\Sigma)$), giving a predictive distribution; the fan charts below are its quantiles. We keep only covariance-stationary draws (companion eigenvalues inside the unit circle) so the simulated paths cannot explode — bvar.filter_stationary.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
import bvar

md = sm.datasets.macrodata.load_pandas().data
gdp  = 400 * np.diff(np.log(md["realgdp"].values))     # annualized real GDP growth, %
infl = 400 * np.diff(np.log(md["cpi"].values))         # annualized CPI inflation, %
tb   = md["tbilrate"].values[1:]                       # 3-month T-bill rate, %
cols = ["GDP growth", "Inflation", "T-bill"]; Y = np.column_stack([gdp, infl, tb])
dates = pd.period_range("1959Q2", periods=len(Y), freq="Q").to_timestamp()
dfY = pd.DataFrame(Y, columns=cols, index=dates)
print("US quarterly %d..%d  T=%d  sample means=%s" % (dates[0].year, dates[-1].year, len(Y), np.round(Y.mean(0), 2)))

fig, ax = plt.subplots(3, 1, figsize=(11, 6), sharex=True)
for i, c in enumerate(cols):
    ax[i].plot(dates, Y[:, i], color="C%d" % i, lw=.9); ax[i].set_ylabel(c); ax[i].axhline(0, color="0.8", lw=.5)
ax[0].set_title("US monetary VAR system (statsmodels macrodata, 1959Q2-2009Q3)")
plt.tight_layout(); plt.show()

US quarterly 1959..2009  T=202  sample means=[3.1  3.98 5.32]

p = 2
res = VAR(dfY).fit(p)
print("Granger causality (does the row variable help predict the other two?):")
for causing in cols:
    caused = [c for c in cols if c != causing]
    pv = res.test_causality(caused, [causing], kind="f").pvalue
    print("  %-10s -> {%s}:  F-test p = %.3f" % (causing, ", ".join(caused), pv))
fevd = res.fevd(20)
print("\nFEVD at 20 quarters (share of each variable's variance from each shock):")
for i, c in enumerate(cols):
    print("  %-10s <-  %s" % (c, dict(zip(cols, np.round(fevd.decomp[i, -1], 2)))))

irf = res.irf(20); orth = irf.orth_irfs; se = irf.stderr(orth=True); mon = 2   # T-bill (monetary) shock, ordered last
fig, ax = plt.subplots(1, 3, figsize=(14, 3.6))
for i, c in enumerate(cols):
    r = orth[:, i, mon]; s = se[:, i, mon]
    ax[i].plot(r, color="navy", lw=1.5); ax[i].fill_between(range(len(r)), r-1.64*s, r+1.64*s, color="navy", alpha=.15)
    ax[i].axhline(0, color="0.6", lw=.6); ax[i].set_title("Response of %s" % c); ax[i].set_xlabel("quarters")
fig.suptitle("Classical VAR: orthogonalized IRF to a monetary (T-bill) shock  [ordering: GDP, Inflation, T-bill]", y=1.02)
plt.tight_layout(); plt.show()

Granger causality (does the row variable help predict the other two?):
  GDP growth -> {Inflation, T-bill}:  F-test p = 0.018
  Inflation  -> {GDP growth, T-bill}:  F-test p = 0.005
  T-bill     -> {GDP growth, Inflation}:  F-test p = 0.003

FEVD at 20 quarters (share of each variable's variance from each shock):
  GDP growth <-  {'GDP growth': 0.9, 'Inflation': 0.06, 'T-bill': 0.04}
  Inflation  <-  {'GDP growth': 0.03, 'Inflation': 0.89, 'T-bill': 0.08}
  T-bill     <-  {'GDP growth': 0.22, 'Inflation': 0.25, 'T-bill': 0.53}

Reading the Granger and FEVD output

Granger causality asks whether a variable's past values help predict the others beyond what their own pasts already do — an $F$-test that the relevant lag coefficients are jointly zero (small $p$ = it helps). It is predictive precedence, not structural causation. Here all three variables Granger-cause the other two ($p=0.018,\,0.005,\,0.003$): output, inflation and the policy rate are mutually predictive — the tightly coupled system a monetary VAR is meant to capture.

FEVD — the forecast-error variance decomposition — splits each variable's $h$-step forecast-error variance across the orthogonal structural shocks (the shares sum to 100%). At 20 quarters, output and inflation are largely own-driven (90% and 89% of their variance from their own shocks), while the policy rate absorbs about a quarter of its variance from output (22%) and inflation (25%) shocks — a Taylor-rule-like reaction, the Fed's rate responding to the economy.

Together with the price-puzzle IRF above, these are a reminder that reduced-form diagnostics describe co-movement; the structural reading (which shock is which) depends entirely on the identification. (Contrast the yield-curve notebook, where the short end Granger-leads but the long end is nearly a random walk.)

mn = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=p, ndraw=4000, burn=1000, seed=1))
ss = bvar.filter_stationary(bvar.gibbs_steady_state(Y, p=p, psi_prior_mean=[2.5, 2.0, 4.0], psi_prior_sd=0.5,
                                                    ndraw=4000, burn=1000, seed=1))
print("Minnesota BVAR:    stationary draws %.0f%%" % (100 * mn["stationary_frac"]))
print("Steady-state BVAR: stationary draws %.0f%%;  psi posterior mean %s   (prior [2.5, 2.0, 4.0])"
      % (100 * ss["stationary_frac"], np.round(ss["mu"].mean(0), 2)))
print("  -> long-run INFLATION pulled from the %.2f%% sample mean down to %.2f%% (toward the 2%% target)"
      % (infl.mean(), ss["mu"].mean(0)[1]))

H = 20; fmn = bvar.forecast(mn, Y, H, seed=2); fss = bvar.forecast(ss, Y, H, seed=2)
fdates = pd.period_range(dates[-1].to_period("Q") + 1, periods=H, freq="Q").to_timestamp()
fig, ax = plt.subplots(1, 3, figsize=(14, 3.8))
for i, c in enumerate(cols):
    ax[i].plot(dates[-40:], Y[-40:, i], color="0.4", lw=.8)
    for f, col, lab in [(fmn, "steelblue", "Minnesota"), (fss, "firebrick", "steady-state")]:
        m_, lo, hi = np.percentile(f[:, :, i], [50, 16, 84], axis=0)
        ax[i].plot(fdates, m_, color=col, lw=1.3, label=lab); ax[i].fill_between(fdates, lo, hi, color=col, alpha=.15)
    ax[i].axhline(ss["mu"].mean(0)[i], color="firebrick", ls=":", lw=.8); ax[i].set_title(c)
ax[1].legend(fontsize=8, loc="upper right")
fig.suptitle("BVAR 5-year forecasts: steady-state reverts to the pinned long-run mean; Minnesota to the sample mean", y=1.02)
plt.tight_layout(); plt.show()

Minnesota BVAR:    stationary draws 97%
Steady-state BVAR: stationary draws 97%;  psi posterior mean [2.9  2.64 3.69]   (prior [2.5, 2.0, 4.0])
  -> long-run INFLATION pulled from the 3.98% sample mean down to 2.64% (toward the 2% target)

# out-of-sample: fit on all but the last 12 quarters (through the 2008 crisis), forecast, compare RMSE
h = 12; Ytr, Yte = Y[:-h], Y[-h:]
rw = np.tile(Ytr[-1], (h, 1))                                        # random walk
ols = VAR(pd.DataFrame(Ytr, columns=cols)).fit(p).forecast(Ytr[-p:], h)
mnT = bvar.filter_stationary(bvar.gibbs_minnesota(Ytr, p=p, ndraw=3000, burn=1000, seed=3))
ssT = bvar.filter_stationary(bvar.gibbs_steady_state(Ytr, p=p, psi_prior_mean=[2.5, 2.0, 4.0], psi_prior_sd=0.5, ndraw=3000, burn=1000, seed=3))
mnF = bvar.forecast(mnT, Ytr, h, seed=4).mean(0); ssF = bvar.forecast(ssT, Ytr, h, seed=4).mean(0)
rmse = lambda F: np.sqrt(((F - Yte) ** 2).mean(0))
print("Out-of-sample RMSE, last 12 quarters (2006Q4-2009Q3, incl. the 2008 crisis):")
print("  %-18s %10s %10s %8s   %s" % ("model", *cols, "overall"))
for name, F in [("random walk", rw), ("OLS-VAR", ols), ("BVAR-Minnesota", mnF), ("BVAR-steady-state", ssF)]:
    r = rmse(F); print("  %-18s %10.2f %10.2f %8.2f   %.3f" % (name, r[0], r[1], r[2], r.mean()))

Out-of-sample RMSE, last 12 quarters (2006Q4-2009Q3, incl. the 2008 crisis):
  model              GDP growth  Inflation   T-bill   overall
  random walk              3.29       5.79     3.28   4.120
  OLS-VAR                  5.14       4.54     2.93   4.203
  BVAR-Minnesota          4.79       4.55     3.10   4.149
  BVAR-steady-state       4.60       4.44     2.65   3.899

Sensitivity — does the 0.8*corr cross-shrinkage matter here?

The engine sets the Minnesota cross-variable tightness to $\lambda_2(i,j)=0.8\,\widehat{\text{corr}}(y_i,y_j)$ — a data-driven choice designed for the smooth age structure of mortality, not the textbook constant $\lambda_2\in[0.2,0.5]$ (Litterman 1986). Since $\widehat{\text{corr}}(\text{GDP},\text{T-bill})\approx0$, this pins the shrinkage on the T-bill$\to$GDP lags near zero. We compare it against constant scalars and check both the coefficient and the forecast:

m = 3; Ytr, Yte = Y[:-12], Y[-12:]
print("lambda2 scheme    T-bill->GDP lag coefs      sum|.|   OOS RMSE")
for name, l2 in [("0.8*corr", None), ("scalar 0.2", 0.2), ("scalar 0.5", 0.5)]:
    B = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=p, lam2=l2, ndraw=4000, burn=1000, seed=1))["B"].mean(0)
    coefs = B[[1 + 0*m + 2, 1 + 1*m + 2], 0]                          # T-bill (var 2) lags 1,2 in the GDP (eq 0)
    fT = bvar.filter_stationary(bvar.gibbs_minnesota(Ytr, p=p, lam2=l2, ndraw=3000, burn=1000, seed=3))
    r = np.sqrt(((bvar.forecast(fT, Ytr, 12, seed=4).mean(0) - Yte) ** 2).mean(0)).mean()
    print("  %-13s %s   %.3f   %.3f" % (name, np.round(coefs, 3), np.abs(coefs).sum(), r))

lambda2 scheme    T-bill->GDP lag coefs      sum|.|   OOS RMSE
  0.8*corr      [-0.013 -0.004]   0.017   4.149
  scalar 0.2    [-0.057 -0.057]   0.115   4.252
  scalar 0.5    [ 0.103 -0.187]   0.290   4.236

• 0.8*corr all but erases the T-bill$\to$GDP channel — summed $|\text{coef}|$ of 0.017 versus 0.29 under $\lambda_2=0.5$. It shrinks a lagged, dynamic relationship to zero on the basis of a near-zero contemporaneous correlation, even though the IRF above shows monetary policy clearly moves output with a lag. If you read the reduced-form coefficients structurally, this choice matters.
• But it does not hurt — it slightly helps — forecasting (OOS RMSE 4.15 vs ~4.25 for the scalars): heavier shrinkage is good insurance through the 2008 window.
• Verdict for these macro series: 0.8*corr is defensible for forecasting but risky for interpreting individual dynamic channels, because it conflates contemporaneous correlation with lagged predictive strength. A constant $\lambda_2$ — or, better, estimating the tightness hierarchically (Giannone-Lenza-Primiceri 2015; Chan 2021, as the R BVAR package does) — is the safer default when the variables are heterogeneous. (The yield-curve notebook shows the opposite case, where 0.8*corr transfers cleanly.)

Results

• The structural story is textbook. All three variables Granger-cause the others; the FEVD shows output and inflation are largely own-driven while the policy rate absorbs a quarter of its variance from output and inflation shocks (a Taylor-rule-like reaction). The IRF to a monetary tightening shows output turning negative after ~3 quarters and a persistent rate response — and the famous "price puzzle" (inflation rising, not falling, after a tightening): a well-known artefact of recursive identification, usually cured by adding a commodity-price index or using sign restrictions. A good reminder that IRFs are only as credible as the identification.
• The steady-state prior is the macro payoff. The Minnesota BVAR reverts, in the long run, to the sample-implied mean — which for 1959–2009 bakes in the high inflation of the 1970s–80s (~3.5–4%). The Villani steady-state BVAR instead reverts to the pinned $\psi$: it pulls long-run inflation down to 2.6%, and its forecasts head back toward $\sim2\%$ — expressing a modern inflation-targeting view the data alone would not give.
• Bayesian shrinkage earns its keep out of sample. Across the last 12 quarters — through the 2008 collapse, a brutal forecasting window — the steady-state BVAR has the lowest RMSE (3.90), beating the OLS-VAR (4.20), the Minnesota BVAR (4.15), and even the random walk (4.12). Anchoring the long run to sensible levels stabilises the medium-term forecast exactly when the unconstrained VAR chases the crisis.

References

• Sims, C. A. (1980). Macroeconomics and reality. Econometrica 48, 1–48.
• Litterman, R. B. (1986). Forecasting with Bayesian vector autoregressions. JBES 4, 25–38.
• Villani, M. (2009). Steady-state priors for vector autoregressions. J. Applied Econometrics 24, 630–650.
• Meseguer, J. (2010). [the bvar.py Minnesota + steady-state Gibbs implementation applied here].

Next: bvar_macro_R.ipynb — the R cross-check (vars + BVAR), then Phase 2 — a finance / yield-curve BVAR.