Cointegration, VECM, and ARDL Models

1. Definition

With \(r\) cointegrating relations, the system is driven by \(n-r\) common stochastic trends (Stock and Watson, 1988). Example: if short and long interest rates are both \(I(1)\) and the term spread is \(I(0)\), then \(\beta=(1,-1)'\), and there is one common trend: the level of interest rates.

2. Error correction and the Granger representation theorem

Start from a VAR(\(p\)) in levels, \(y_t=c+A_1y_{t-1}+\dots+A_py_{t-p}+u_t\). It can always be rewritten as

• If \(\operatorname{rank}(\Pi)=0\), there is no cointegration: a VAR in differences is appropriate.
• If \(\operatorname{rank}(\Pi)=n\), then \(y_t\) is stationary: a VAR in levels is appropriate.
• If \(0<\operatorname{rank}(\Pi)=r<n\), then \(\Pi=\alpha\beta'\), with \(\alpha\) and \(\beta\) both \(n\times r\).

Interpretation. \(\beta'y_{t-1}\) measures disequilibrium, and the loading (adjustment) coefficients \(\alpha\) say how fast each variable moves to restore equilibrium. If a row of \(\alpha\) is zero, that variable is weakly exogenous for \(\beta\): it does not respond to disequilibrium, and the other variables do the adjusting. Weak exogeneity of a variable means it drives the common trend.

Consequences for VARs in levels. A VAR in levels remains consistent when variables are cointegrated, and many coefficient tests are standard (Sims, Stock and Watson, 1990). Some tests, including Granger causality, can have non-standard distributions. The Toda–Yamamoto (1995) lag-augmentation approach restores standard Wald inference. Differencing a cointegrated system throws away the error-correction term: a misspecification.

3. Engle–Granger two-step procedure

Step 1. Estimate the levels regression \(y_{1t}=\mu+\gamma'y_{2t}+z_t\) by OLS and save \(\widehat z_t\).

Step 2. Test \(\widehat z_t\) for a unit root with an ADF regression without deterministics. Under the null of no cointegration, the test statistic has a distribution that depends on the number of regressors and the deterministics, because \(\widehat z_t\) is constructed to look as stationary as possible. Critical values come from Engle and Yoo (1987) or MacKinnon (2010); for two variables with a constant, the 5% value is about \(-3.34\), against \(-2.86\) for an ordinary ADF test.

Super-consistency and inference. If cointegrated, \(\widehat\gamma\) converges at rate \(T\). But its limit distribution is generally non-normal and depends on nuisance parameters (endogeneity of \(\Delta y_{2t}\) and serial correlation), so OLS \(t\)-statistics in the levels regression are invalid. Use:

• dynamic OLS (DOLS; Stock and Watson, 1993): add leads and lags of \(\Delta y_{2t}\), \(y_{1t}=\mu+\gamma'y_{2t}+\sum_{j=-q}^{q}\delta_j'\Delta y_{2,t-j}+e_t\), with HAC standard errors. Inference on \(\gamma\) is then asymptotically normal;
• fully modified OLS (FMOLS; Phillips and Hansen, 1990): non-parametric corrections for the same nuisance parameters.

Limitations. Results can depend on which variable is on the left. Only one cointegrating vector is estimated. In small samples, the static regression can have substantial bias.

4. Johansen's maximum-likelihood approach

Johansen (1988, 1991) estimates the VECM by Gaussian (reduced-rank) maximum likelihood. Concentrating out the short-run terms, regress \(\Delta y_t\) and \(y_{t-1}\) on the lagged differences and deterministics, and keep the residuals \(R_{0t}\) and \(R_{1t}\). With moment matrices \(S_{ij}=T^{-1}\sum_tR_{it}R_{jt}'\), solve the eigenvalue problem

giving ordered eigenvalues \(\widehat\lambda_1>\dots>\widehat\lambda_n\), the squared canonical correlations between \(\Delta y_t\) and \(y_{t-1}\). The MLE of \(\beta\) is given by the corresponding eigenvectors.

Rank tests.

Test sequentially: \(r=0\), then \(r\le1\), and so on, and stop at the first non-rejection. The distributions are non-standard and depend on the deterministic specification (no constant; restricted constant in the cointegrating space; unrestricted constant, giving linear trends in levels; restricted trend; unrestricted trend). The usual choice for trending macro data is an unrestricted constant or a restricted trend. The Johansen trace test is known to over-reject in small samples; Bartlett corrections or bootstrap tests help (Cavaliere, Rahbek and Taylor, 2012).

Identifying \(\beta\). Only the cointegrating space is identified: \(\alpha\beta'=(\alpha H^{-1\prime})(\beta H)'\) for any non-singular \(r\times r\) matrix \(H\). With \(r=1\), normalize one coefficient to 1. With \(r>1\), economically motivated restrictions (Johansen and Juselius, 1992) are needed, tested by LR tests that are \(\chi^2\)-distributed, because inference on \(\beta\) given the rank is asymptotically mixed normal.

5. ARDL models and the bounds test

The autoregressive distributed lag model

is a single-equation conditional ECM. The long-run coefficients are \(-\theta/\phi\), and \(\phi\in(-2,0)\) is the speed of adjustment.

Bounds test (Pesaran, Shin and Smith, 2001). Test \(H_0:\phi=\theta=0\) (no level relationship) with an \(F\) statistic, and the \(t\)-statistic on \(\phi\). The critical values depend on whether the regressors are \(I(0)\) or \(I(1)\). PSS tabulate lower bounds (all regressors \(I(0)\)) and upper bounds (all \(I(1)\)). An \(F\) above the upper bound indicates a level relationship regardless of the regressors' integration order; below the lower bound, none; in between, the test is inconclusive. Use finite-sample critical values (Narayan, 2005; Kripfganz and Schneider, 2020) when \(T\) is below about 100.

Assumptions that are often overlooked. The bounds test (i) requires that \(y_t\) is \(I(1)\) or \(I(0)\) and that the regressors are not \(I(2)\); (ii) assumes at most one level relationship, with \(x_t\) weakly exogenous (\(x\) does not error-correct toward \(y\)); and (iii) is a single-equation method. If there may be several cointegrating relations, or feedback from \(y\) to \(x\), use a Johansen VECM instead. "We use ARDL because our variables are a mix of I(0) and I(1)" is a valid reason only when these assumptions hold.

import numpy as np
from statsmodels.tsa.stattools import coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM

rng = np.random.default_rng(14)
T = 300
trend = np.cumsum(rng.normal(size=T))              # one common stochastic trend
x = trend + rng.normal(0, 0.5, T)
y = 2.0 + 1.5 * trend + rng.normal(0, 0.5, T)       # beta = (1, -1.5) up to noise in x

t_eg, p_eg, crit = coint(y, x)
print(f"Engle-Granger t = {t_eg:.2f}, p = {p_eg:.3f}, 5% cv = {crit[1]:.2f}")

j = coint_johansen(np.column_stack([y, x]), det_order=0, k_ar_diff=1)
print("trace stats:", j.lr1.round(1), " 5% cvs:", j.cvt[:, 1])
print("normalized beta:", np.real(j.evec[:, 0] / j.evec[0, 0]).round(3))

m = VECM(np.column_stack([y, x]), k_ar_diff=1, coint_rank=1, deterministic="ci").fit()
print("alpha (loadings):", m.alpha.ravel().round(3))

import statsmodels.api as sm
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM
d = sm.datasets.danish_data.load_pandas().data[["lrm", "lry", "ibo", "ide"]]
j = coint_johansen(d, det_order=0, k_ar_diff=1)
# TODO: print j.lr1 against j.cvt[:,1]; VECM(d, k_ar_diff=1, coint_rank=1, deterministic="ci").fit()

import numpy as np, pandas as pd, statsmodels.api as sm
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM
d = sm.datasets.danish_data.load_pandas().data[["lrm", "lry", "ibo", "ide"]]
j = coint_johansen(d, det_order=0, k_ar_diff=1)
print(pd.DataFrame({"H0: r <=": range(4), "trace": j.lr1.round(2), "5% cv": j.cvt[:, 1]}).to_string(index=False))
v = VECM(d, k_ar_diff=1, coint_rank=1, deterministic="ci").fit()
beta = v.beta[:, 0] / v.beta[0, 0]
print("\ncointegrating vector (normalized on lrm):", dict(zip(d.columns, beta.round(3))))
print("long run: lrm =", " + ".join(f"{-b:.3f}*{c}" for b, c in zip(beta[1:], d.columns[1:])), "+ const")
print("adjustment coefficients alpha:", dict(zip(d.columns, v.alpha[:, 0].round(3))))

H0: r <=  trace   5% cv
        0  48.80 47.8545
        1  17.29 29.7961
        2   7.14 15.4943
        3   0.56  3.8415

cointegrating vector (normalized on lrm): {'lrm': np.float64(1.0), 'lry': np.float64(-0.969), 'ibo': np.float64(5.403), 'ide': np.float64(-4.14)}
long run: lrm = 0.969*lry + -5.403*ibo + 4.140*ide + const
adjustment coefficients alpha: {'lrm': np.float64(-0.3), 'lry': np.float64(0.027), 'ibo': np.float64(0.004), 'ide': np.float64(0.02)}

import numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf
from statsmodels.tsa.stattools import coint, adfuller
m = sm.datasets.macrodata.load_pandas().data
d = pd.DataFrame({"c": np.log(m["realcons"]), "y": np.log(m["realdpi"])})
# TODO: coint(d.c, d.y); DOLS: add dy leads/lags; ECM with lagged residual

import numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf
from statsmodels.tsa.stattools import coint, adfuller
m = sm.datasets.macrodata.load_pandas().data
d = pd.DataFrame({"c": np.log(m["realcons"]), "y": np.log(m["realdpi"])})
for v in ("c", "y"):
    a = adfuller(d[v], regression="ct", autolag="AIC"); print(f"ADF(ct) {v}: t = {a[0]:.2f}, p = {a[1]:.3f}")
t, p, cv = coint(d["c"], d["y"])
print(f"Engle-Granger: t = {t:.2f}, p = {p:.4f}, 5% cv = {cv[1]:.2f}")

static = smf.ols("c ~ y", data=d).fit()
d["dy"] = d["y"].diff()
for k in (-2, -1, 1, 2):
    d[f"dy_{k}"] = d["dy"].shift(k)
d = d.rename(columns={"dy_-1": "dy_m1", "dy_-2": "dy_m2"})         # leads of dy (k = -1, -2)
dols = smf.ols("c ~ y + dy + dy_1 + dy_2 + dy_m1 + dy_m2", data=d.dropna()).fit(cov_type="HAC", cov_kwds={"maxlags": 4})
print(f"static OLS slope {static.params['y']:.4f};  DOLS slope {dols.params['y']:.4f} (HAC se {dols.bse['y']:.4f})")
w = dols.t_test("y = 1")
print(f"DOLS test of slope = 1: t = {float(np.squeeze(w.tvalue)):.2f}")

d["z"] = static.resid; d["dc"] = d["c"].diff(); d["z1"] = d["z"].shift(1)
ecm = smf.ols("dc ~ z1 + dy", data=d.dropna(subset=["dc", "z1", "dy"])).fit(cov_type="HAC", cov_kwds={"maxlags": 4})
g = ecm.params["z1"]
print(f"ECM: adjustment {g:.4f} (t = {ecm.tvalues['z1']:.2f}), short-run income effect {ecm.params['dy']:.3f}")
print(f"half-life of a deviation from equilibrium: {np.log(0.5)/np.log(1+g):.1f} quarters")

ADF(ct) c: t = -2.54, p = 0.309
ADF(ct) y: t = -1.62, p = 0.785
Engle-Granger: t = -2.54, p = 0.2620, 5% cv = -3.37
static OLS slope 1.0320;  DOLS slope 1.0359 (HAC se 0.0077)
DOLS test of slope = 1: t = 4.65
ECM: adjustment -0.0473 (t = -1.74), short-run income effect 0.360
half-life of a deviation from equilibrium: 14.3 quarters