Bayesian VECM — I. Term structure — R cross-check¶

Vector autoregressions, Part 6a (R)¶

Validates bvecm_termstructure_python.ipynb with urca (Pfaff), the reference R implementation of Johansen's procedure. ca.jo runs the reduced-rank maximum-likelihood analysis: the trace test for the cointegration rank and the eigenvectors that are the cointegrating vectors. If it finds the same rank ($r=2$) and the same stationary error-correction relations as the from-scratch Bayesian VECM, the result holds across method and implementation.

Data: bvar_yields_data.csv — 3m / 5y / 10y / 30y Treasury yields, monthly 1994–2024 (FRED).

In [1]:
.libPaths('C:/Users/user/R/win-library/4.6')
suppressMessages(library(urca))
dfy <- read.csv('bvar_yields_data.csv'); Y <- dfy[, c('y3m','y5y','y10y','y30y')]
for (nm in colnames(Y)) {
  s <- summary(ur.df(Y[[nm]], type = 'drift', lags = 4, selectlags = 'AIC'))
  cat(sprintf('ADF %-6s: stat %6.2f  (5%% crit %.2f)  -> %s\n', nm, s@teststat[1], s@cval[1,2],
              ifelse(s@teststat[1] > s@cval[1,2], 'I(1)', 'I(0)')))
}
jo <- ca.jo(Y, type = 'trace', ecdet = 'none', K = 2, spec = 'transitory')
cat('Johansen trace test:\n')
print(round(cbind(trace = jo@teststat, jo@cval), 1))
cat('\ncointegrating vectors (first r=2, normalised):\n')
print(round(jo@V[, 1:2], 2))
cat('cross-check -> Python (Bayesian) also selects r=2; EC-term std 0.93 / 0.17 (stationary)\n')

ADF y3m   : stat  -2.07  (5% crit -2.87)  -> I(1)
ADF y5y   : stat  -2.04  (5% crit -2.87)  -> I(1)
ADF y10y  : stat  -2.20  (5% crit -2.87)  -> I(1)
ADF y30y  : stat  -2.23  (5% crit -2.87)  -> I(1)

Johansen trace test:

         trace 10pct 5pct 1pct
r <= 3 |   3.6   6.5  8.2 11.7
r <= 2 |  10.5  15.7 18.0 23.5
r <= 1 |  40.1  28.7 31.5 37.2
r = 0  |  82.2  45.2 48.3 55.4

cointegrating vectors (first r=2, normalised):

        y3m.l1  y5y.l1
y3m.l1    1.00    1.00
y5y.l1   -3.81  -81.97
y10y.l1   4.10  188.41
y30y.l1  -1.13 -111.85

cross-check -> Python (Bayesian) also selects r=2; EC-term std 0.93 / 0.17 (stationary)
In [2]:
# the two error-correction terms = the cointegrating combinations of the (levelled) yields -> stationary
options(repr.plot.width = 13, repr.plot.height = 4)
z1 <- as.matrix(Y) %*% jo@V[, 1]; z2 <- as.matrix(Y) %*% jo@V[, 2]
par(mfrow = c(1, 2))
matplot(as.Date(dfy$date), as.matrix(Y), type = 'l', lty = 1, lwd = 1.3, xlab = '', ylab = '%',
        main = 'Four yields — each I(1), wandering')
legend('topright', c('3m','5y','10y','30y'), col = 1:4, lty = 1, bty = 'n', cex = .8)
matplot(as.Date(dfy$date), cbind(z1 - mean(z1), z2 - mean(z2)), type = 'l', lty = 1, lwd = 1.3,
        col = c('firebrick','navy'), xlab = '', ylab = '', main = "Error-correction terms (beta'y) — stationary")
abline(h = 0, col = 'grey60'); legend('topright', c('EC 1','EC 2'), col = c('firebrick','navy'), lty = 1, bty = 'n', cex = .8)
No description has been provided for this image

Same rank, same stationary relations¶

  • urca confirms $r=2$. The trace statistic rejects $r=0$ and $r\le 1$ but not $r\le 2$ — exactly the Bayesian VECM's finding: four I(1) yields, two cointegrating relations, two common stochastic trends.
  • The error-correction terms are stationary. The two cointegrating combinations of the yields (right panel) revert to a fixed level while the yields themselves wander (left) — the maximum-likelihood counterpart of the Bayesian posterior EC terms.
  • As in the Python notebook, only the two-dimensional cointegrating space is identified; urca's normalised eigenvectors are a particular basis for it, not uniquely meaningful coefficients.

References¶

  • Pfaff, B. (2008). Analysis of Integrated and Cointegrated Time Series with R. Springer (urca).
  • Johansen, S. (1991). Econometrica 59, 1551–1580.