Bayesian Vector Error Correction — I. The term structure¶

Vector autoregressions, Part 6a: cointegration¶

Every VAR in this series has treated its variables as either stationary (growth rates) or, in the large-BVAR notebook, as unit-root levels shrunk toward random walks. Cointegration is the case in between and the most economically loaded: several series each wander like a random walk (I(1)), yet specific linear combinations of them are stationary — pinned down by a long-run equilibrium the economy is always pulled back toward. The natural laboratory is the term structure of interest rates: the 3-month, 5-, 10- and 30-year Treasury yields each drift over decades, but the expectations hypothesis says they cannot drift apart without bound — the spreads between them are mean-reverting. If so, the yields are cointegrated, share a small number of common stochastic trends, and must be modelled in error-correction form.

This is the first of three cointegration notebooks (Part 6). It carries the full model and the from-scratch Bayesian sampler; the companions apply the same machinery to the Fisher effect and money demand — two long-run relations the data will treat very differently.

Data: US Treasury constant-maturity yields (3m, 5y, 10y, 30y), monthly 1994–2024 (FRED).

The model: VAR → VECM, and the reduced-rank restriction¶

Take the level VAR($p$) and difference it into error-correction form: $$\Delta y_t = c + \Pi\, y_{t-1} + \sum_{i=1}^{p-1}\Gamma_i\,\Delta y_{t-i} + \varepsilon_t,\qquad \varepsilon_t\sim\mathcal N(0,\Sigma),$$ with $\Pi = -(I-A_1-\dots-A_p)$ the long-run impact matrix. Every term but $\Pi y_{t-1}$ is stationary; that single levels term carries the cointegration. Its rank classifies the system: $\text{rank}\,\Pi=0$ → no cointegration (a VAR in differences); $\text{rank}\,\Pi=m$ → the levels are already stationary (a VAR in levels); and $0<r<m$ → $r$ cointegrating relations, with $$\boxed{\Pi=\alpha\beta'},\qquad \alpha,\beta\ \text{are } m\times r.$$ $\beta$ holds the cointegrating vectors — $\beta'y_t$ is stationary, the $r$ combinations that don't wander — and $\alpha$ the adjustment loadings, the speed at which each variable is pulled back when the system is off its equilibrium $\beta'y_{t-1}$.

Identification. $\alpha\beta'=(\alpha M)(M^{-1}\beta')'$ for any invertible $M$, so only the column space of $\beta$ is data-identified. We impose the linear normalization $\beta=[\,I_r\,;\,B\,]$ (fix the top $r\times r$ block to the identity, estimate the free $(m-r)\times r$ block $B$). This pins down a representative $\beta$ but — a point that matters below — the individual coefficients of $\beta$ are a convention, whereas the stationarity of $\beta'y_t$ is not.

The Bayesian sampler (from scratch — bvecm.py)¶

Given the rank $r$ and the normalization, the model is conditionally linear, so a three-block Gibbs sampler with conjugate draws follows — the same machinery as the other samplers in this series, in a new self-contained module (bvar.py untouched):

  1. $\Sigma \mid$ rest — inverse-Wishart on the VECM residuals.
  2. $[\alpha,\Gamma,c]\mid\beta,\Sigma$ — with $\beta$ fixed, $\beta'y_{t-1}$ is a regressor, so this is a multivariate-Normal (SUR) draw.
  3. $B\mid\alpha,\Gamma,\Sigma$ — the term $\alpha\beta'y_{t-1}=\alpha[I_r;B]'y_{t-1}$ is linear in $B$, giving a Normal conditional.

Priors: Normal on $B$ (diffuse, or centred on a theory value), Normal on $\alpha,\Gamma$, inverse-Wishart on $\Sigma$. The sampler is initialised at the Johansen reduced-rank estimate. We report the frequentist Johansen trace test alongside as the standard rank reference.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.vector_ar.vecm import coint_johansen
import bvecm

dfy = pd.read_csv("bvar_yields_data.csv", parse_dates=["date"], index_col="date")
cols = ["y3m", "y5y", "y10y", "y30y"]; Y = dfy[cols].values; dates = dfy.index; m = 4; p = 2; r = 2

# --- ADF integration-order pre-test (cointegration requires I(1) inputs) ---
from statsmodels.tsa.stattools import adfuller
print("ADF unit-root pre-test (H0: unit root / I(1)):")
for j, cc in enumerate(cols):
    a0 = adfuller(Y[:, j], maxlag=4, autolag="AIC"); a1 = adfuller(np.diff(Y[:, j]), maxlag=4, autolag="AIC")
    print("  %-5s level ADF %6.2f p=%.3f -> %s | diff ADF %6.2f p=%.3f" % (
        cc, a0[0], a0[1], "I(1)" if a0[1] > 0.05 else "I(0)", a1[0], a1[1]))
print("  => all four yields are I(1): unit root in levels, stationary in differences\n")

# --- Johansen trace test (frequentist rank reference) ---
jo = coint_johansen(Y, det_order=0, k_ar_diff=1)
print("JOHANSEN trace test (term structure):")
for i in range(m):
    print("  r<=%d : trace %6.1f  95%%crit %5.1f  %s" % (i, jo.lr1[i], jo.cvt[i, 1],
          "reject" if jo.lr1[i] > jo.cvt[i, 1] else "ACCEPT -> r=%d" % i))
print("  => cointegration rank r = 2\n")

# --- Bayesian VECM (from scratch): 3-block Gibbs at r=2 ---
fit = bvecm.gibbs_vecm(Y, p=p, r=r, ndraw=4000, burn=2000, seed=1)
beta = np.median(fit["beta"], 0); alpha = np.median(fit["alpha"], 0)
ec = bvecm.ec_term(fit, Y)                                        # (n, r) error-correction terms  beta'y_t
print("Bayesian beta (median):\n", np.round(beta, 2))
print("alpha (median, adjustment speeds):\n", np.round(alpha, 3))
print("EC-term std", np.round(ec.std(0), 2), " vs  yield-level std", np.round(Y.std(0), 2))

# figure: the yields wander (I(1)); the r error-correction terms are stationary  ->  vecm_ts_1.png
_dec = dates[p - 1:len(Y) - 1]                                    # dates aligned to Y_{t-1} in the EC term
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
_ystyle = [("y3m", "firebrick"), ("y5y", "navy"), ("y10y", "seagreen"), ("y30y", "darkorange")]
for _k, (_cc, _col) in enumerate(_ystyle):
    ax1.plot(dates, Y[:, _k], color=_col, lw=1.6, label=_cc)
ax1.set_title("Treasury yields: each wanders (I(1))")
ax1.set_ylabel("percent"); ax1.grid(alpha=.25); ax1.legend(ncol=2, fontsize=8)
_ecstyle = ["firebrick", "navy"]
for _j in range(ec.shape[1]):
    ax2.plot(_dec, ec[:, _j], color=_ecstyle[_j], lw=1.6, label=r"EC term %d ($\beta'y$)" % (_j + 1))
    ax2.axhline(ec[:, _j].mean(), color=_ecstyle[_j], lw=.8, ls=":", alpha=.6)
ax2.set_title("Error-correction terms: stationary, revert to a fixed level")
ax2.grid(alpha=.25); ax2.legend(fontsize=8)
fig.tight_layout(); plt.show()

ADF unit-root pre-test (H0: unit root / I(1)):
  y3m   level ADF  -2.05 p=0.265 -> I(1) | diff ADF  -6.78 p=0.000
  y5y   level ADF  -1.78 p=0.393 -> I(1) | diff ADF -16.95 p=0.000
  y10y  level ADF  -2.15 p=0.226 -> I(1) | diff ADF -10.69 p=0.000
  y30y  level ADF  -2.17 p=0.219 -> I(1) | diff ADF -10.78 p=0.000
  => all four yields are I(1): unit root in levels, stationary in differences

JOHANSEN trace test (term structure):
  r<=0 : trace   82.2  95%crit  47.9  reject
  r<=1 : trace   40.1  95%crit  29.8  reject
  r<=2 : trace   10.5  95%crit  15.5  ACCEPT -> r=2
  r<=3 : trace    3.6  95%crit   3.8  ACCEPT -> r=3
  => cointegration rank r = 2


Bayesian beta (median):
 [[ 1.    0.  ]
 [ 0.    1.  ]
 [-4.87 -2.38]
 [ 4.18  1.41]]
alpha (median, adjustment speeds):
 [[-0.087  0.3  ]
 [ 0.006  0.171]
 [ 0.006  0.184]
 [ 0.007  0.099]]
EC-term std [0.93 0.17]  vs  yield-level std [2.12 1.9  1.69 1.51]
No description has been provided for this image

Reading the cointegration¶

  • The yields wander; their combinations don't. Each yield has a level standard deviation of ~1.5–2.1 and drifts over the sample (left panel). The two error-correction terms $\beta'y_t$ have standard deviations of only 0.93 and 0.17 and visibly revert to a fixed level (right panel). That gap is cointegration: four I(1) yields, but only two independent stochastic trends, so $4-2=2$ stationary long-run relations hold the curve together — the Johansen trace test agrees at $r=2$.
  • Only the space is identified, not the coefficients. The normalised $\beta$ has odd-looking entries (e.g. $-4.87, 4.18$ on the long yields) rather than clean $(1,-1)$ spreads. That is expected: with $r=2$ the data identify the two-dimensional cointegrating space, not individual vectors — any rotation $\beta M$ spans the same space and yields the same (stationary) fit. What is economically real is that two stationary combinations exist, not the particular basis the normalisation prints.
  • The adjustment loadings $\alpha$ show who corrects. The short rate carries the largest loadings — the front end of the curve does most of the error-correcting back toward equilibrium, while the long end (30y) barely adjusts, consistent with long rates behaving as the near-random-walk "level" of the curve.
In [2]:
from statsmodels.tsa.vector_ar.vecm import VECM
from statsmodels.tsa.api import VAR
# Out-of-sample race: does keeping the cointegration help forecasts? Recursive, 2012-2024.
# VECM (rank 2) vs VAR-in-levels (ignores unit roots) vs VAR-in-differences (ignores the level relations).
H = 12; hshow = [1, 3, 6, 12]; start = int(np.where(dates >= pd.Timestamp("2012-01-01"))[0][0])
sse = {k: {h: [] for h in hshow} for k in ["VECM", "VAR-diff", "VAR-levels"]}
for t in range(start, len(Y) - H):
    tr = Y[:t]
    fV = VECM(tr, k_ar_diff=1, coint_rank=2, deterministic="co").fit().predict(steps=H)
    fL = VAR(tr).fit(p).forecast(tr[-p:], H)
    dtr = np.diff(tr, axis=0); fD = tr[-1] + np.cumsum(VAR(dtr).fit(1).forecast(dtr[-1:], H), 0)
    for h in hshow:
        a = Y[t + h - 1]
        sse["VECM"][h].append(np.mean((fV[h-1]-a)**2)); sse["VAR-levels"][h].append(np.mean((fL[h-1]-a)**2))
        sse["VAR-diff"][h].append(np.mean((fD[h-1]-a)**2))
rmse = {k: [np.sqrt(np.mean(sse[k][h])) for h in hshow] for k in sse}
print("OOS RMSE (avg over yields) by horizon (months):", hshow)
for k in ["VAR-levels", "VAR-diff", "VECM"]: print("  %-11s" % k, np.round(rmse[k], 3))

# figure: RMSE vs forecast horizon, one line per model  ->  vecm_ts_2.png
_models = [("VECM", "firebrick"), ("VAR-diff", "navy"), ("VAR-levels", "seagreen")]
fig, ax = plt.subplots(figsize=(8, 4))
for _name, _col in _models:
    ax.plot(hshow, rmse[_name], color=_col, lw=1.6, marker="o", ms=4, label=_name)
ax.set_xticks(hshow)
ax.set_xlabel("forecast horizon (months)"); ax.set_ylabel("RMSE")
ax.grid(alpha=.25); ax.legend(title="model")
ax.set_title("Out-of-sample RMSE by forecast horizon (term structure)")
fig.tight_layout(); plt.show()

OOS RMSE (avg over yields) by horizon (months): [1, 3, 6, 12]
  VAR-levels  [0.217 0.417 0.63  0.984]
  VAR-diff    [0.214 0.412 0.644 1.077]
  VECM        [0.216 0.415 0.632 1.009]
No description has been provided for this image

Results — cointegration is real, and it matters most for not over-differencing¶

  • The yield curve is cointegrated ($r=2$). Both the Johansen trace test and the Bayesian VECM find two stationary long-run relations among the four yields — the curve has only two common stochastic trends, and the spreads (in some basis) are mean-reverting, as the expectations hypothesis predicts.
  • Forecasting: keep the levels. At one year out the VAR-in-differences is clearly worst (RMSE 1.08) — differencing away the level relations discards exactly the information that anchors long-horizon forecasts (at very short horizons, one-to-three months, it is actually a shade best: the level relations bind only as the horizon lengthens). The VECM (1.01) and the levels-VAR (0.98) track together and both beat it: the honest lesson (Christoffersen-Diebold 1998) is that the payoff of cointegration is chiefly in not throwing the levels away, while imposing the rank restriction on top of a levels VAR is a smaller, sample-dependent refinement.
  • The Bayesian value-added. Beyond a point estimate, the Gibbs sampler delivers the full posterior of the adjustment speeds and the cointegrating space, and a coherent framework for rank and prior information — where Johansen offers a sequence of tests with non-standard critical values. The R cross-check reproduces the rank ($r=2$) and the stationary error-correction relations with urca.

The next two notebooks turn the same machine on relations the data treat very differently: the Fisher effect (Part 6b) and money demand (Part 6c).

References¶

  • Engle, R. F. & Granger, C. W. J. (1987). Co-integration and error correction. Econometrica 55, 251–276.
  • Johansen, S. (1991). Estimation and hypothesis testing of cointegration vectors in Gaussian VAR models. Econometrica 59, 1551–1580.
  • Campbell, J. Y. & Shiller, R. J. (1987). Cointegration and tests of present value models. J. Political Economy 95, 1062–1088.
  • Christoffersen, P. F. & Diebold, F. X. (1998). Cointegration and long-horizon forecasting. J. Business & Economic Statistics 16, 450–458.

Next: bvecm_termstructure_R.ipynb — the R cross-check (urca::ca.jo).