Large Bayesian VAR — shrinkage that scales with the system¶
Vector autoregressions, Part 3¶
Parts 1–2 kept the VAR small (3–5 variables) so every coefficient could be read. But macroeconomists rarely believe the world is 3-dimensional: a monetary shock moves output, prices, employment, housing, credit, exchange rates and term spreads all at once, and a small VAR that omits them can mistake an omitted-variable correlation for a structural one (the price puzzle of Part 2 is partly this). The natural fix — add more variables — runs into over-parameterization. An $m$-variable VAR with $p$ lags has $m(1+mp)$ coefficients, growing with the square of $m$ (quadratically — so, despite the folklore, this is not literally the "curse of dimensionality," which is reserved for costs that grow exponentially in the dimension). Our 20-variable, 4-lag system below has 1,620 coefficients against ~240 quarterly observations. Least squares can still estimate that — it stays consistent and roughly unbiased — but it pays a ruinous sampling variance: each coefficient is mostly estimation noise, that noise compounds through the forecast iteration, and a fit that looks superb in-sample generalises terribly out of sample. The macro data sharpens the problem: the series are highly collinear and near-unit-root, so the regressor matrix is badly conditioned and the OLS estimates turn unstable. The failure of a big unrestricted VAR is a bias–variance problem — too many freely-estimated, collinear coefficients relative to the signal — not a dimensionality one.
Bańbura, Giannone & Reichlin (2010) showed the way through: keep all the variables, but shrink harder as the system grows. A large VAR under a suitably tight Minnesota prior is not only estimable, it forecasts as well as or better than a small VAR — and better than factor models — precisely because the prior disciplines the 1,620 coefficients while letting the data speak where it is informative. This notebook reproduces that result from scratch with the bvar.py Gibbs sampler, on the FRED-QD macroeconomic panel.
The model and the data¶
The VAR in levels. We estimate a reduced-form VAR($p$) with $p=4$ (one year of quarterly lags), $$y_t = c + B_1 y_{t-1} + \dots + B_p y_{t-p} + \varepsilon_t,\qquad \varepsilon_t\sim\mathcal N(0,\Sigma),$$ where $y_t$ collects $m$ macro series. Following BGR we do not difference to stationarity; we keep the variables in (log-)levels and let the prior handle the near-unit-roots. Real, nominal and price series enter as $100\log$; interest rates and the unemployment rate enter in levels.
The Minnesota prior (Litterman 1986), unit-root form. Each equation's coefficients are shrunk toward a random walk: the prior mean is $1$ on a variable's own first lag and $0$ on everything else (own_mean = 1), so absent data every series is an independent random walk — the natural no-change benchmark for persistent macro data. The prior standard deviation on the coefficient linking equation $i$ to lag $\ell$ of variable $j$ is
$$
\operatorname{sd}(B_\ell^{ij}) \;=\;
\begin{cases}
\dfrac{\lambda_1}{\ell^{\,\lambda_3}}, & i=j \quad(\text{own lags}),\\[2mm]
\dfrac{\lambda_1\,\lambda_2}{\ell^{\,\lambda_3}}\cdot\dfrac{s_i}{s_j}, & i\neq j \quad(\text{cross lags}),
\end{cases}
$$
with all four ingredients defined explicitly:
- $\lambda_1$ — overall tightness. The master knob: how far, in total, coefficients may stray from the random-walk prior. This is what BGR make a function of the system size $m$ (next section); small $\lambda_1$ = heavy shrinkage.
- $\lambda_2 = 0.5$ — cross-variable shrinkage. Other variables' lags are shrunk an extra factor of $\lambda_2$ harder than a variable's own lags (own dynamics are trusted more than cross dynamics). We use the traditional scalar $0.5$.
- $\lambda_3 = 1$ — lag decay. Shrinkage tightens linearly with lag $\ell$: distant lags are pulled to zero faster, encoding "recent lags matter more."
- $s_i/s_j$ — scale correction. The ratio of AR residual standard deviations, so the prior is invariant to the units of each series (a rate in % vs a log-level ×100).
Estimation. The posterior is sampled from scratch by the two-block Gibbs sampler in bvar.py (gibbs_minnesota): draw the coefficients given $\Sigma$ from their Normal conditional, draw $\Sigma$ given the coefficients from an inverse-Wishart, repeat. Nothing here is hand-tuned per variable — the prior formula above scales automatically to any $m$.
Three systems, one panel. From FRED-QD (McCracken-Ng, 1959Q1–2019Q4, pre-COVID) we build a nested ladder:
| system | $m$ | coefficients | variables |
|---|---|---|---|
| small | 3 | 39 | real GDP, GDP deflator, Fed funds |
| medium | 7 | 203 | + consumption, investment, hours, unemployment |
| large | 20 | 1,620 | + IP, payrolls, housing, M1/M2, 3m/1y/10y rates, CPI, PPI, oil, USD, avg. hours |
The three share the same three forecast targets — real GDP, the deflator, the funds rate — so we can ask the BGR question directly: does adding the other 17 variables help or hurt the forecasts of the core three?
Setting the one knob that matters: the BGR fit rule for $\lambda_1$¶
A tighter prior always improves out-of-sample stability but worsens in-sample fit; a looser prior does the reverse. The size of the system changes this trade-off — a 20-variable VAR at $\lambda_1=0.2$ fits the in-sample data far more flexibly than a 3-variable VAR at the same $\lambda_1$, simply because it has more regressors. BGR's insight is to neutralise that mechanical advantage: choose $\lambda_1$ for each system so that they all achieve the same relative in-sample fit on the core variables.
Concretely, define the in-sample overfit of a fitted model as the average, over the core variables, of its one-step in-sample mean-squared error scaled by the variable's variance, $$\text{Fit}(\lambda_1)=\frac1{|\text{core}|}\sum_{i\in\text{core}}\frac{\text{MSE}_i(\lambda_1)}{\operatorname{Var}(y_i)}\;\approx\;1-\bar R^2 .$$ Pick a reference fit from the small VAR at a conventional $\lambda_1=0.2$, then solve $\text{Fit}_{\text{medium}}(\lambda_1)=\text{Fit}_{\text{small}}(0.2)$ and likewise for the large system. Because bigger systems reach any given fit at a smaller $\lambda_1$, the rule automatically tightens the prior as variables are added — the mathematical content of "shrink harder when the system is larger."
import numpy as np, pandas as pd
import matplotlib.pyplot as plt
import bvar
d = pd.read_csv("bvar_largebvar_data.csv", index_col=0, parse_dates=True)
cols = list(d.columns); Y = d.values; p = 4
SMALL = ["GDPC1", "GDPCTPI", "FEDFUNDS"]
MEDIUM = SMALL + ["PCECC96", "GPDIC1", "HOANBS", "UNRATE"]
LARGE = cols # all 20
CORE = SMALL # forecast targets
idx = {"small": [cols.index(c) for c in SMALL],
"medium": [cols.index(c) for c in MEDIUM], "large": list(range(len(cols)))}
def insample_fit(Yk, lam1, core_local, ndraw=300):
# posterior-mean 1-step in-sample MSE of the core vars / their variance (~ 1 - R^2)
fit = bvar.gibbs_minnesota(Yk, p=p, lam1=lam1, lam2=0.5, lam3=1.0, own_mean=1.0,
ndraw=ndraw, burn=ndraw//3, seed=1)
Yt, X = bvar.build_regressors(Yk, p, intercept=True)
resid = Yt - X @ fit["B"].mean(0)
return np.mean([resid[:, j].var() / Yt[:, j].var() for j in core_local])
LGRID = [0.05, 0.1, 0.15, 0.2, 0.3, 0.5]
core_i = [cols.index(c) for c in CORE]
fitc = {}
for name in ["small", "medium", "large"]:
Yk = Y[:, idx[name]]; cl = [idx[name].index(i) for i in core_i]
fitc[name] = [insample_fit(Yk, lam, cl) for lam in LGRID]
ref = np.interp(0.2, LGRID, fitc["small"]) # target: small VAR at lambda=0.2
def match(name): # invert Fit(lambda)=ref
fc = np.array(fitc[name]); o = np.argsort(fc)
return float(np.interp(ref, fc[o], np.array(LGRID)[o]))
lam_sel = {"small": 0.2, "medium": match("medium"), "large": match("large")}
print("lambda_1 grid:", LGRID)
for name in ["small", "medium", "large"]:
print(" %-6s (n=%2d): %s" % (name, len(idx[name]), np.round(fitc[name], 3).tolist()))
print("reference fit (small@0.2): %.3f" % ref)
print("selected lambda_1:", {k: round(v, 3) for k, v in lam_sel.items()})
# in-sample overfit Fit(lambda_1) per system, with the reference level and the
# lambda_1 each system needs to match it (BGR: bigger system -> tighter prior).
sys_style = {"small": ("navy", "small (m=3)"),
"medium": ("seagreen", "medium (m=7)"),
"large": ("firebrick", "large (m=20)")}
fig, ax = plt.subplots(figsize=(7.2, 4.6))
for name, (col, lab) in sys_style.items():
ax.plot(LGRID, fitc[name], marker="o", ms=5, lw=1.8, color=col, label=lab)
ax.axvline(lam_sel[name], color=col, ls="--", lw=1.0, alpha=.7) # selected lambda_1
ax.axhline(ref, color="k", lw=.8, ls=":", label="reference fit (small@0.2)")
ax.set_xlabel(r"$\lambda_1$ (overall prior tightness)")
ax.set_ylabel(r"in-sample overfit $\approx\,1-\bar R^2$")
ax.set_title("BGR fit rule: matching in-sample fit forces a tighter prior as the system grows")
ax.grid(alpha=.25)
ax.legend(frameon=False, fontsize=9)
fig.tight_layout(); plt.show()
lambda_1 grid: [0.05, 0.1, 0.15, 0.2, 0.3, 0.5]
small (n= 3): [0.019, 0.018, 0.017, 0.016, 0.016, 0.015]
medium (n= 7): [0.019, 0.017, 0.016, 0.015, 0.014, 0.014]
large (n=20): [0.017, 0.015, 0.013, 0.012, 0.011, 0.009]
reference fit (small@0.2): 0.016
selected lambda_1: {'small': 0.2, 'medium': 0.135, 'large': 0.062}
Reading the fit rule¶
The overfit measure falls as $\lambda_1$ rises (a looser prior fits the sample better), and — at any given $\lambda_1$ — the large system sits well below the small one: with 1,620 coefficients it can shadow the in-sample data almost regardless of the prior. Holding the fit equal to the small VAR's reference level therefore forces the large system onto a much tighter $\lambda_1$ than the small one (the dashed verticals), with the medium system in between. This is the BGR prescription made quantitative: the prior tightness that keeps a large VAR honest is several times smaller than what a small VAR needs. The R cross-check confirms that a completely different criterion — maximising the marginal likelihood — lands in the same tight region for the 20-variable system.
from statsmodels.tsa.api import VAR
from statsmodels.tsa.ar_model import AutoReg
import matplotlib.pyplot as plt
# recursive (expanding-window) pseudo-out-of-sample, 2004Q1-2019Q4, forecasting the 3 core targets.
# competitors: AR(4) per variable; unrestricted OLS VAR (small & large); BVAR small/medium/large
# (each BVAR uses the lambda_1 chosen above). One Gibbs fit per origin serves all horizons via bvar.forecast.
def bvar_point(name, Yhist, H=4):
fit = bvar.gibbs_minnesota(Yhist, p=p, lam1=lam_sel[name], lam2=0.5, lam3=1.0, own_mean=1.0,
ndraw=250, burn=80, seed=1)
return bvar.forecast(fit, Yhist, H=H, seed=1).mean(0) # (H, m) posterior-mean path
# ... loop over origins, accumulate squared errors, form RMSE (see repo script) ...
rmse = pd.read_csv("_largebvar_rmse.csv") # precomputed on the full recursive run
rel = rmse.copy()
for m_ in ["AR4","VAR-OLS-small","VAR-OLS-large","BVAR-small","BVAR-medium","BVAR-large"]:
rel[m_] = rmse[m_] / rmse["AR4"]
print(rel.round(2).to_string(index=False))
# relative RMSE (vs AR(4)) by model, one panel per forecast horizon; dotted line = AR(4) = 1.0
models = ["VAR-OLS-small", "VAR-OLS-large", "BVAR-small", "BVAR-medium", "BVAR-large"]
mcol = {"VAR-OLS-small": "navy", "VAR-OLS-large": "grey",
"BVAR-small": "seagreen", "BVAR-medium": "goldenrod", "BVAR-large": "firebrick"}
core_lab = {"GDPC1": "real GDP", "GDPCTPI": "deflator", "FEDFUNDS": "funds rate"}
fig, axes = plt.subplots(1, 2, figsize=(13, 4.2), sharey=True)
for ax, h in zip(axes, [1, 4]):
sub = rel[rel["h"] == h]
xpos = np.arange(len(sub)); w = 0.15
for k, mdl in enumerate(models):
ax.bar(xpos + (k - (len(models) - 1) / 2) * w, sub[mdl].values, w,
color=mcol[mdl], label=mdl)
ax.axhline(1.0, color="k", lw=.8, ls=":") # AR(4) benchmark
ax.set_xticks(xpos); ax.set_xticklabels([core_lab[v] for v in sub["var"]])
ax.set_title("h = %d quarter%s ahead" % (h, "" if h == 1 else "s"))
ax.grid(alpha=.25, axis="y")
axes[0].set_ylabel("RMSE relative to AR(4)")
axes[1].legend(frameon=False, fontsize=8.5, loc="upper left")
fig.suptitle("Out-of-sample forecast RMSE relative to the AR(4) benchmark", y=1.02)
fig.tight_layout(); plt.show()
h var AR4 VAR-OLS-small VAR-OLS-large BVAR-small BVAR-medium BVAR-large 1 GDPC1 1.0 1.12 1.28 1.09 1.18 1.13 1 GDPCTPI 1.0 1.01 1.41 0.99 1.03 1.09 1 FEDFUNDS 1.0 1.14 2.52 0.89 1.06 1.13 4 GDPC1 1.0 1.11 1.53 1.11 1.20 1.04 4 GDPCTPI 1.0 1.12 1.92 1.13 1.31 1.42 4 FEDFUNDS 1.0 0.88 2.19 0.83 0.94 1.03
Seeing the forecasts¶
The RMSE table is decisive but abstract; a picture makes the bias–variance story concrete. Below, a 20-variable BVAR and the same system by unrestricted OLS are both trained through 2014Q4 and asked to project the next five years. The shrunk BVAR (red, with 68% / 90% predictive bands) stays anchored to the realised paths of output, prices and the funds rate. The failure of unrestricted OLS (navy dashed) is starkest for the funds rate: it projects a nonsensical −4.4% policy rate by 2019 — a negative interest rate the sample never came close to — while the realised value was 1.6% and the BVAR a sensible 3.7%. Same variables, same sample, same horizon; the only difference is whether the 1,620 coefficients were disciplined, and for the one variable with genuinely distinct dynamics the undisciplined version runs right off the rails.
One detail foreshadows Part 4: the BVAR's predictive bands here are essentially constant width, because this model assumes a constant shock covariance $\Sigma$. Real macro data is not homoskedastic — those bands should be wide in the 1970s and 2008 and narrow in the Great Moderation. That is exactly what stochastic volatility adds next.
from statsmodels.tsa.api import VAR
import matplotlib.pyplot as plt
# train the large system through 2014Q4, forecast the 3 core variables to 2019Q4, with predictive bands
split = np.where(d.index <= pd.Timestamp("2014-12-01"))[0][-1] + 1
Ytr = Y[:split]; H = len(d) - split
fit = bvar.gibbs_minnesota(Ytr, p=p, lam1=lam_sel["large"], lam2=0.5, lam3=1.0, own_mean=1.0,
ndraw=2000, burn=800, seed=1)
paths = bvar.forecast(fit, Ytr, H=H, seed=1) # (ndraw, H, m) predictive draws
ols = VAR(pd.DataFrame(Ytr, columns=cols)).fit(p).forecast(Ytr[-p:], H) # unrestricted OLS large VAR
for c in CORE:
j = cols.index(c)
print("%-9s 2019Q4: actual %.1f | BVAR %.1f | OLS %.1f" %
(c, Y[-1, j], np.median(paths[:, :, j], 0)[-1], ols[-1, j]))
# 5/16/50/84/95 percentile fan per core variable vs realised, OLS overlaid
q = np.percentile(paths, [5, 16, 50, 84, 95], axis=0) # (5, H, m)
dates = d.index
start = np.where(dates >= pd.Timestamp("2005-01-01"))[0][0] # history window for readability
cd = dates[split - 1:] # forecast axis incl. anchor point
lab = {"GDPC1": "real GDP (100·log)", "GDPCTPI": "GDP deflator (100·log)",
"FEDFUNDS": "Fed funds rate (%)"}
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, c in zip(axes, CORE):
j = cols.index(c)
anc = Y[split - 1, j]
cc = lambda a: np.concatenate([[anc], a]) # connect lines to the last in-sample point
ax.plot(dates[start:split], Y[start:split, j], color="0.35", lw=1.4) # in-sample history
ax.plot(cd, Y[split - 1:, j], color="0.10", lw=1.7, label="realised") # realised out-of-sample path
ax.fill_between(cd, cc(q[0, :, j]), cc(q[4, :, j]), color="firebrick", alpha=.15, lw=0, label="90% band")
ax.fill_between(cd, cc(q[1, :, j]), cc(q[3, :, j]), color="firebrick", alpha=.15, lw=0, label="68% band")
ax.plot(cd, cc(q[2, :, j]), color="firebrick", lw=1.9, label="BVAR median") # shrunk forecast
ax.plot(cd, cc(ols[:, j]), color="navy", lw=1.6, ls="--", label="OLS large VAR")
ax.axvline(dates[split - 1], color="k", lw=.6, ls=":") # forecast origin
ax.set_title(lab[c]); ax.grid(alpha=.25)
axes[0].legend(frameon=False, fontsize=8, loc="upper left")
fig.suptitle("20-variable system trained through 2014Q4, projected to 2019Q4: shrinkage vs unrestricted OLS", y=1.02)
fig.tight_layout(); plt.show()
GDPC1 2019Q4: actual 995.0 | BVAR 989.6 | OLS 985.9 GDPCTPI 2019Q4: actual 465.0 | BVAR 468.1 | OLS 461.0 FEDFUNDS 2019Q4: actual 1.6 | BVAR 3.7 | OLS -4.4
Results — shrinkage turns an unusable large VAR into a usable one¶
- Unrestricted OLS breaks down as the system grows. The 20-variable OLS VAR (grey) is the worst forecaster everywhere — 1.3–1.9× the AR(4) benchmark for output and prices and a catastrophic 2.5× for the funds rate at one quarter, 2.2× at one year. This is not the curse of dimensionality (the coefficient count grows only quadratically, and OLS stays estimable and consistent) but an estimation-variance / overfitting failure: 1,620 freely-estimated, highly-collinear coefficients are dominated by sampling noise that compounds through the forecast horizon. Tellingly, OLS blows up worst on the funds rate — the one target whose dynamics differ from the near-random-walk crowd, exactly where the ill-conditioned regressor matrix bites hardest. A big VAR estimated by least squares is simply not an option.
- The same 20 variables under shrinkage are competitive. Swapping OLS for the BGR-tightened Minnesota prior — nothing else changes — collapses the large system's relative RMSE from 1.3–2.5 down to 1.0–1.4, right back into the pack with the small models. The 1,620 coefficients are a liability or an asset depending entirely on whether they are disciplined; shrinkage, not the variable count, is what matters.
- Adding variables is free, not harmful — but here it is not a clear win either. The large BVAR (red) lands in the same ballpark as the 3-variable BVAR (navy): a touch worse at $h=1$, comparable at $h=4$ (it beats the small BVAR for one-year GDP, 1.04 vs 1.11, and trails it for prices). So the extra 17 series are carried at no cost to the core forecasts, but in this quarterly, levels, pre-COVID sample they do not deliver the decisive edge BGR report for monthly data and wider panels — an honest reminder that the large-VAR dividend is horizon-, variable- and design-specific.
- The small BVAR is the quiet winner for the funds rate — the most dynamically-rich target, and the one place cross-variable dynamics beat the AR(4) benchmark (BVAR-small 0.89 at $h=1$, 0.83 at $h=4$; the unrestricted small VAR also dips just under, ~0.88 at $h=4$). For the near-random-walk output and price levels, every method hugs the no-change benchmark, as expected.
- Independent confirmation (R). The
BVARpackage, choosing $\lambda_1$ by marginal likelihood rather than by our fit rule, tightens to the same region for the large system (≈0.10) and produces credible out-of-sample forecasts of the core variables — the shrink-with-size mechanism is not an artefact of one implementation.
The practical template behind modern large-BVAR forecasting and "big-data" macro (Giannone-Lenza-Primiceri 2015 turn the fit rule into a full hierarchical prior on $\lambda_1$) is exactly this: you do not have to choose between a small interpretable VAR and a big informative one — take the big one and shrink it correctly, and at worst you lose nothing on the variables you care about while gaining a coherent joint model of everything else.
References¶
- Bańbura, M., Giannone, D. & Reichlin, L. (2010). Large Bayesian vector auto regressions. J. Applied Econometrics 25, 71–92.
- Litterman, R. B. (1986). Forecasting with Bayesian vector autoregressions — five years of experience. J. Business & Economic Statistics 4, 25–38.
- Giannone, D., Lenza, M. & Primiceri, G. (2015). Prior selection for vector autoregressions. Review of Economics and Statistics 97, 436–451.
- McCracken, M. W. & Ng, S. (2020). FRED-QD: A quarterly database for macroeconomic research. NBER WP 26872.
Next: bvar_largebvar_R.ipynb — the R cross-check (BVAR package: marginal-likelihood $\lambda_1$ and out-of-sample forecasts).