Bayesian VAR — finance application: the Treasury yield curve

Vector autoregressions, part 2: a yield-curve VAR with the Minnesota and steady-state priors

Part 1 applied the from-scratch Bayesian-VAR engine (bvar.py) to a macro monetary VAR. Part 2 turns it on a finance system — the US Treasury yield curve — where two things change the flavour: the series are near-unit-root levels (yields, not growth rates), and the natural "steady state" is a neutral yield curve. The same three models — classical VAR, Minnesota-prior BVAR, and the Villani steady-state BVAR — now answer finance questions: how does a short-rate (policy) shock propagate across the curve, and how do we forecast the curve while anchoring its long-run shape?

The system. Four Treasury constant-maturity yields, monthly, 1994–2024 (yfinance; saved to bvar_yields_data.csv), spanning the curve:

• 3-month (^IRX), 5-year (^FVX), 10-year (^TNX), 30-year (^TYX), all in percent.

The sample straddles three very different rate regimes — the 5–6% late-1990s, the near-zero ZIRP 2010s, and the 2022–24 hiking cycle — which is exactly why the sample average is a poor guide to the long run, and where an informative steady-state prior can help.

Models (recap) and the finance-specific setup

The three models are those derived in full in bvar_macro_python.ipynb; here is the compact statement plus what changes for yields.

VAR($p$) in levels. $y_t = c + \sum_{l=1}^p \Pi_l y_{t-l} + \varepsilon_t$, $\varepsilon_t\sim N(0,\Sigma)$, on the four yields (levels, not differences). Structural shocks via the Cholesky factor $\Sigma=PP'$ with a short$\to$long ordering [3m, 5y, 10y, 30y]: the short-rate (policy) shock is ordered first, so it may move the whole curve within the month.

Minnesota-prior BVAR (gibbs_minnesota, 2-block Gibbs). Same prior $\operatorname{vec}(B)\sim N(b_0,V_0)$ with std's $\lambda_1/l^{\lambda_3}$ (own) and $\lambda_1\lambda_2(i,j)/l^{\lambda_3}\cdot s_i/s_j$ (cross), where $\lambda_1$ is overall tightness, $\lambda_2(i,j)=0.8\,\widehat{\text{corr}}$ the cross-variable tightness, $\lambda_3$ the lag decay. What changes: because yields are highly persistent, the own-first-lag prior is centred near a random walk, own_mean = 0.9 (versus $0$ for the differenced macro data) — the benchmark each yield is shrunk toward is "close to last month," not "revert to a mean."

Villani steady-state BVAR (gibbs_steady_state, 3-block Gibbs). $y_t=\psi+\sum_l\Pi_l(y_{t-l}-\psi)+\varepsilon_t$ with an informative $\psi\sim N(\psi_0,\Lambda_0)$. The finance use: $\psi$ is the long-run (neutral) yield curve. We pin an upward-sloping neutral curve $\psi_0=[2.5,\,3.3,\,3.8,\,4.2]$ (a ~2.5% neutral short rate rising to ~4.2% at 30y), so forecasts revert to that rather than to the regime-blended sample average.

Estimation, stationarity filtering (companion eigenvalues inside the unit circle), and simulation forecasting are exactly as in Part 1.

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

df = pd.read_csv("bvar_yields_data.csv", parse_dates=["date"], index_col="date")
cols = ["y3m", "y5y", "y10y", "y30y"]; lab = ["3-month", "5-year", "10-year", "30-year"]
Y = df[cols].values; dates = df.index
print("Treasury yields monthly %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(figsize=(11, 4))
for i, l in enumerate(lab): ax.plot(dates, Y[:, i], lw=.9, label=l)
ax.legend(fontsize=8, ncol=4, loc="upper right"); ax.set_ylabel("yield (%)")
ax.set_title("US Treasury yield curve (monthly, 1994-2024)")
plt.tight_layout(); plt.show()

Treasury yields monthly 1994..2024  T=372  sample means=[2.38 3.35 3.83 4.35]

The term structure over time

Before modelling, look at the raw object the VAR is trying to capture. Two views: the curve's shape at notable dates — a normal upward curve, a flat one, an inverted one — and the 10-year minus 3-month slope through time, the classic recession bellwether (it inverts before downturns).

mats = np.array([0.25, 5, 10, 30])
fig, ax = plt.subplots(1, 2, figsize=(14, 4.4), gridspec_kw=dict(width_ratios=[1, 1.4]))
snaps = {"2000-12 (inverted, pre-dot-com)": "2000-12-31", "2004-06 (steep, easy Fed)": "2004-06-30",
         "2007-02 (flat, pre-GFC)": "2007-02-28", "2013-06 (ZIRP)": "2013-06-30",
         "2023-07 (deep inversion)": "2023-07-31", "2024-12 (latest)": "2024-12-31"}
for l, dt in snaps.items():
    row = df.loc[dt]; ax[0].plot(mats, [row[c] for c in cols], marker="o", lw=1.4, ms=4, label=l)
ax[0].set_xscale("log"); ax[0].set_xticks(mats); ax[0].set_xticklabels(["3m", "5y", "10y", "30y"])
ax[0].set_xlabel("maturity"); ax[0].set_ylabel("yield (%)"); ax[0].legend(fontsize=7)
ax[0].set_title("The yield curve at notable dates (normal / flat / inverted)")
slope = df["y10y"] - df["y3m"]
ax[1].plot(dates, slope, color="navy", lw=.9); ax[1].axhline(0, color="0.5", lw=.7)
ax[1].fill_between(dates, slope, 0, where=slope < 0, color="firebrick", alpha=.35, label="inverted (10y<3m)")
ax[1].set_ylabel("10y - 3m spread (%)"); ax[1].legend(fontsize=8, loc="upper right")
ax[1].set_title("Term-structure slope over time: inversions (red) precede the 2001, 2008, 2020 recessions")
print("inverted in %d of %d months; deepest %.2f%% in %s" % ((slope < 0).sum(), len(slope), slope.min(), slope.idxmin().strftime("%Y-%m")))
plt.tight_layout(); plt.show()

inverted in 46 of 372 months; deepest -1.61% in 2023-05

p = 2; res = VAR(pd.DataFrame(Y, columns=cols, index=dates)).fit(p)
print("Granger causality (does the variable help predict the rest?):")
for i, c in enumerate(cols):
    caused = [x for x in cols if x != c]
    print("  %-8s -> rest:  F-test p = %.3f" % (lab[i], res.test_causality(caused, [c], kind="f").pvalue))
fevd = res.fevd(24)
print("\nFEVD at 24 months -- share of each yield's variance from the short-rate (3m) shock:")
print("  ", dict(zip(lab, np.round(fevd.decomp[:, -1, 0], 2))))

irf = res.irf(24); orth = irf.orth_irfs; se = irf.stderr(orth=True); sh = 0     # 3m (short-rate) shock, ordered first
fig, ax = plt.subplots(1, 4, figsize=(15, 3.3))
for i, l in enumerate(lab):
    r = orth[:, i, sh]; s = se[:, i, sh]
    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" % l); ax[i].set_xlabel("months")
fig.suptitle("Classical VAR: yield-curve response to a short-rate (3-month) shock  [ordering short->long]", y=1.03)
plt.tight_layout(); plt.show()

Granger causality (does the variable help predict the rest?):
  3-month  -> rest:  F-test p = 0.002
  5-year   -> rest:  F-test p = 0.000
  10-year  -> rest:  F-test p = 0.068
  30-year  -> rest:  F-test p = 0.030

FEVD at 24 months -- share of each yield's variance from the short-rate (3m) shock:
   {'3-month': 0.45, '5-year': 0.21, '10-year': 0.16, '30-year': 0.11}

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 the short and medium tenors lead (3-month $p=0.002$, 5-year $p<0.001$), while the 10-year is only marginal ($p=0.068$): the long end behaves like a near-random walk driven by its own news, not something the front end foreshadows.

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%). The row above shows a short-rate (3-month) shock explains 45% of the 3-month's 24-month forecast variance but only 11% of the 30-year's — policy dominates the front end and its influence fades monotonically toward the long end, whose variance comes mostly from its own (expectations / term-premium) shocks.

The two diagnostics tell the same story from different angles: the Fed moves the short end of the curve; the long end has a mind of its own. (Contrast the macro notebook, where all three variables Granger-caused each other and shocks were more mutually shared.)

psi = [2.5, 3.3, 3.8, 4.2]                                        # pinned neutral curve
mn = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=p, own_mean=0.9, ndraw=4000, burn=1000, seed=1))
ss = bvar.filter_stationary(bvar.gibbs_steady_state(Y, p=p, own_mean=0.9, psi_prior_mean=psi, 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 %s)"
      % (100 * ss["stationary_frac"], np.round(ss["mu"].mean(0), 2), psi))
print("  sample-mean 3m yield %.2f%% (blends 90s ~5%% and ZIRP ~0%%) -> steady-state anchors it at %.2f%%"
      % (Y[:, 0].mean(), ss["mu"].mean(0)[0]))

H = 36; fmn = bvar.forecast(mn, Y, H, seed=2); fss = bvar.forecast(ss, Y, H, seed=2)
fdates = pd.period_range(dates[-1].to_period("M") + 1, periods=H, freq="M").to_timestamp()
fig, ax = plt.subplots(1, 4, figsize=(15, 3.4))
for i, l in enumerate(lab):
    ax[i].plot(dates[-60:], Y[-60:, i], color="0.4", lw=.8)
    for f, col, nm 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=nm); 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(l)
ax[0].legend(fontsize=8, loc="lower left")
fig.suptitle("BVAR yield-curve forecasts revert slowly toward the anchored neutral curve (dotted) -- here the two priors nearly agree,\n"
             "because the 1994-2024 sample level already resembles a neutral curve (unlike the macro-inflation case)", y=1.06)
plt.tight_layout(); plt.show()

Minnesota BVAR:    stationary draws 93%
Steady-state BVAR: stationary draws 93%;  psi posterior mean [2.48 3.31 3.75 4.24]   (prior [2.5, 3.3, 3.8, 4.2])
  sample-mean 3m yield 2.38% (blends 90s ~5% and ZIRP ~0%) -> steady-state anchors it at 2.48%

# out-of-sample: hold out the last 24 months, forecast, compare RMSE (yields are near-random-walk)
h = 24; Ytr, Yte = Y[:-h], Y[-h:]
rw = np.tile(Ytr[-1], (h, 1))
ols = VAR(pd.DataFrame(Ytr, columns=cols)).fit(p).forecast(Ytr[-p:], h)
mnT = bvar.filter_stationary(bvar.gibbs_minnesota(Ytr, p=p, own_mean=0.9, ndraw=3000, burn=1000, seed=3))
ssT = bvar.filter_stationary(bvar.gibbs_steady_state(Ytr, p=p, own_mean=0.9, psi_prior_mean=psi, 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 24 months:")
print("  %-18s %s   overall" % ("model", "  ".join("%7s" % l for l in lab)))
for name, F in [("random walk", rw), ("OLS-VAR", ols), ("BVAR-Minnesota", mnF), ("BVAR-steady-state", ssF)]:
    r = rmse(F); print("  %-18s %s   %.3f" % (name, "  ".join("%7.2f" % x for x in r), r.mean()))

Out-of-sample RMSE, last 24 months:
  model              3-month   5-year  10-year  30-year   overall
  random walk           0.81     0.37     0.44     0.47   0.521
  OLS-VAR               1.37     0.51     0.45     0.38   0.679
  BVAR-Minnesota        1.19     0.40     0.37     0.32   0.568
  BVAR-steady-state     1.13     0.38     0.35     0.32   0.545

Sensitivity — the 0.8*corr cross-shrinkage is inconsequential here

The Minnesota cross-variable tightness $\lambda_2(i,j)=0.8\,\widehat{\text{corr}}(y_i,y_j)$ is a mortality-motivated choice. For the yield curve it gives loose values ($\lambda_2\approx0.58$–$0.79$, since neighbouring tenors correlate ~0.9), close to a symmetric prior. We compare it against constant scalars on the 3m$\to$10y cross-lags and the forecast:

m = 4; Ytr, Yte = Y[:-24], Y[-24:]
print("lambda2 scheme    3m->10y 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, own_mean=0.9, lam2=l2, ndraw=4000, burn=1000, seed=1))["B"].mean(0)
    coefs = B[[1 + 0*m + 0, 1 + 1*m + 0], 2]                          # 3m (var 0) lags 1,2 in the 10y (eq 2)
    fT = bvar.filter_stationary(bvar.gibbs_minnesota(Ytr, p=p, own_mean=0.9, lam2=l2, ndraw=3000, burn=1000, seed=3))
    r = np.sqrt(((bvar.forecast(fT, Ytr, 24, 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    3m->10y lag coefs         sum|.|   OOS RMSE
  0.8*corr      [ 0.015 -0.001]   0.016   0.568
  scalar 0.2    [-0.005  0.019]   0.024   0.587
  scalar 0.5    [ 0.007  0.006]   0.013   0.573

• The three schemes give near-identical cross-coefficients (~0.02) and OOS RMSE (~0.57) — for the yield curve the cross-shrinkage choice does not move the needle. With own_mean=0.9 each yield is explained mostly by its own lag, and the tenors are so correlated that any reasonable $\lambda_2$ lands in the same place: the data swamp the prior.
• This is the mirror image of the macro case (bvar_macro_python.ipynb): there, 0.8*corr erased a genuine lagged channel because the variables were heterogeneous with near-zero contemporaneous correlation; here, on a smooth, highly-correlated curve — structurally the same object as the mortality age ladder the rule was built for — it transfers cleanly and harmlessly.

Results

• Policy pass-through decays up the curve. The FEVD shows a short-rate shock explains 45% of the 3-month's variance but only 11% of the 30-year's, and the IRF shows the whole curve rising after a short-rate shock but by less at the long end (~0.23 at 3m vs ~0.07 at 30y). This is the textbook monetary-transmission / curve-flattening picture: the front end tracks policy, the long end is anchored by expectations and the term premium. Granger causality agrees — the short/medium tenors lead, the long end is only weakly predictable.
• The steady-state prior barely moves here — and that is the lesson. The pinned neutral curve $[2.5,3.3,3.8,4.2]$ almost coincides with the sample means $[2.38,3.35,3.83,4.35]$, so the steady-state and Minnesota forecasts nearly overlap. Contrast Part 1, where the 1959–2009 sample inflation (3.98%) was far above a 2% target and the steady-state prior mattered a lot. The steady-state prior earns its keep exactly when the sample average is unrepresentative of the future — here, by luck of the sample, it isn't.
• Yields are near-random-walk, but BVAR still tames the VAR. Out of sample the random walk wins (RMSE 0.521) — no surprise for near-unit-root yields — but among the models, the steady-state BVAR (0.545) clearly beats the unconstrained OLS-VAR (0.679) and the Minnesota BVAR (0.568). With four variables the OLS-VAR overfits and drifts; Minnesota/steady-state shrinkage keeps it disciplined. This is the general BVAR message, sharper here than in Part 1 because the system is larger.

References

• 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.
• Diebold, F. X. & Li, C. (2006). Forecasting the term structure of government bond yields. J. Econometrics 130, 337–364.

Next: bvar_yields_R.ipynb — the R cross-check (vars + BVAR).