Cointegration: two random walks that cannot lose sight of each other

Why a regression between two non-stationary series almost always lies — and the one precise case where it tells the truth. Interactive simulators, the Engle–Granger test, error correction models, and Python code that runs directly in this page.

📅 August 2026 ⏱️ ~18 min read 🏷️ Time series • Unit roots • VECM

← Back to the blog Jump to the simulator 🇫🇷 Lire en français

Contents

Why this matters
The problem: spurious regression
Order of integration
Defining cointegration
Simulator: the leash and the dog
The Engle–Granger test
Error correction: the representation theorem
Beyond two series: Johansen
In practice (R and Python)
Limits and pitfalls
References

Why this matters

Two curves have been rising together for ten years. A commodity price and a sector index, electricity consumption and average temperature, two stocks on the same market. The temptation is irresistible: there must be a link, and you can trade on it. A regression confirms it, with a flattering \(R^2\) and a highly significant coefficient.

Except that this reasoning is wrong roughly three times out of four. Between two series that drift freely — two random walks with no relationship whatsoever — the t statistic crosses the significance threshold in more than 80% of cases instead of the advertised 5%. And it gets worse as the sample grows: the more data you have, the more solid the false relationship looks. This is the result that earned Granger a Nobel prize, and it is probably the most expensive statistical mistake ever made in finance.

Cointegration is the answer: the exact criterion separating real relationships from mirages. It says there exists a combination of the two series that does not drift — a gap that always returns to its mean, like a dog held on a leash by its walker. If that leash exists, you can bet on reversion to equilibrium. If it doesn't, you are betting on nothing.

On this page you will draw your own pairs of independent series and watch convincing "relationships" appear before your eyes. Then you will cut the leash on a genuinely linked pair, by dragging a slider, and watch the test verdict flip live. The code driving all of it is right here, editable and runnable.

Unit root Engle–Granger VECM Johansen Spurious regression

The problem: spurious regression

Take two series with nothing in common. Two pure random walks, produced by two independent generators:

Two independent random walks

\[X_{t1} = X_{t-1,1} + u_t, \qquad X_{t2} = X_{t-1,2} + v_t, \qquad \operatorname{Cov}(u_t, v_s) = 0 \;\;\forall t,s\]

Regress the first on the second. Least squares theory says the slope should be zero, the \(R^2\) nil, the t statistic small. It is not. Granger and Newbold (1974) showed it by simulation, Phillips (1986) proved it: under a unit root the \(t\) statistic diverges at rate \(\sqrt{n}\) instead of converging to a normal law. The longer the sample, the more significant the false relationship looks.

The simulator below draws two independent random walks and regresses one on the other. Run it a few times: convincing "relationships" will appear. Then launch the Monte Carlo experiment to measure the true rejection rate of a test whose nominal level is 5%.

The takeaway. The t test is not "slightly" too permissive: its actual rejection rate exceeds 70% instead of 5%, and it gets worse as \(n\) grows. Any regression between non-stationary series in levels is suspect until cointegration has been checked. The classic symptom in a report: a very high \(R^2\) paired with a Durbin–Watson statistic close to zero.

Order of integration

Vocabulary first. A series \(X_t\) is integrated of order \(d\), written \(X_t \sim I(d)\), if it must be differenced \(d\) times to become stationary:

Order of integration

\[(1-B)^d X_t \;\text{is stationary}, \qquad (1-B)^{d-1} X_t \;\text{is not}\]

White noise is \(I(0)\); a random walk is \(I(1)\), since its first difference is white noise. In the multivariate VARMA formalism this corresponds to the diagonal operator \(\lambda(B) = \operatorname{diag}\!\big((1-B)^{d_1}, \dots, (1-B)^{d_m}\big)\) of the VARIMA model \(\Phi(B)\lambda(B)(\mathbf{X}_t - \boldsymbol\mu) = \Theta(B)\mathbf{a}_t\), where each component may carry its own differencing order \(d_j\).

Two arithmetic rules give the intuition for everything that follows. If \(X_t \sim I(0)\) and \(Y_t \sim I(0)\), any linear combination is \(I(0)\). If \(X_t \sim I(1)\) and \(Y_t \sim I(1)\), then in general \(aX_t + bY_t \sim I(1)\): the non-stationarity does not cancel. Cointegration is precisely the exception to that rule.

Defining cointegration

Definition (Engle and Granger, 1987)

The components of the vector \(\mathbf{X}_t = (X_{t1}, \dots, X_{tm})\'\) are said to be cointegrated of order \((d, b)\), written \(\mathbf{X}_t \sim CI(d,b)\), if:

(i) every component is \(I(d)\); (ii) there exists a non-zero vector \(\boldsymbol\beta = (\beta_1,\dots,\beta_m)\'\) such that

\[z_t = \boldsymbol\beta\' \mathbf{X}_t \sim I(d-b), \qquad b > 0.\]

The case that matters in practice is \(CI(1,1)\): every series is \(I(1)\), but one particular combination is \(I(0)\), hence stationary. The vector \(\boldsymbol\beta\) is called the cointegrating vector, and \(z_t\) the equilibrium error: the momentary departure from a long-run relationship \(\boldsymbol\beta\' \mathbf{X}_t = 0\).

Note that \(\boldsymbol\beta\) is only defined up to a multiplicative constant: if \(\boldsymbol\beta\' \mathbf{X}_t\) is stationary, so is \(2\boldsymbol\beta\' \mathbf{X}_t\). The first coefficient is therefore normalised to 1 by convention, which in the bivariate case gives \(z_t = X_{t1} - \beta X_{t2}\).

The drunk and her dog. The image is due to Murray (1994). The walker follows a random walk: she has no destination. Neither does the dog. Taken separately, both trajectories are \(I(1)\) and unpredictable. But the leash forces their gap to stay bounded: \(z_t\) is stationary. Without a leash the two walks drift apart and the gap is itself \(I(1)\) — no cointegration. The whole statistical question is: does that leash exist, or am I looking at two walkers who happen to be heading the same way?

Simulator: the leash and the dog

The second simulator builds a pair whose leash you control. The series \(X_{t2}\) is a pure random walk — the common stochastic trend. The equilibrium error follows an AR(1) with parameter \(\phi = 1 - \alpha\), and \(X_{t1}\) is rebuilt as \(X_{t1} = \beta X_{t2} + z_t\):

The simulated model

\[X_{t2} = X_{t-1,2} + u_t, \qquad z_t = (1-\alpha)\,z_{t-1} + e_t, \qquad X_{t1} = \beta X_{t2} + z_t\]

The parameter \(\alpha\) is the speed of adjustment. When \(\alpha = 0\) the leash is cut: \(z_t\) becomes a random walk and the pair is no longer cointegrated, even though the two curves may still look very much alike. As \(\alpha\) grows, the gap is pulled back to zero faster and faster. Drag the \(\alpha\) slider towards 0 and watch the test statistic flip.

β 1.00 α (adjustment) 0.15 σ_e 1.0 n 240

Three things in the readout are worth watching. First, the least squares estimate of \(\beta\) is excellent even when the noise is enormous: this is superconsistency, \(\hat\beta\) converges at rate \(n\) rather than \(\sqrt{n}\). Second, the half-life \(\ln(0.5)/\ln(\hat\phi)\) gives the number of periods needed to absorb half a shock: it is the first quantity a pairs trader looks at. Third, the test statistic is compared with \(-3.34\) and not with \(-1.96\); the next section explains why.

The Engle–Granger test

The two-step procedure is straightforward. First estimate the long-run relationship by ordinary least squares, then test whether the residual is stationary:

Step 1 — long-run relationship

\[X_{t1} = \mu + \beta X_{t2} + z_t \quad\Longrightarrow\quad \hat z_t = X_{t1} - \hat\mu - \hat\beta X_{t2}\]

Step 2 — unit root test on the residual (ADF)

\[\Delta \hat z_t = \rho\, \hat z_{t-1} + \sum_{i=1}^{p} \gamma_i \Delta \hat z_{t-i} + \varepsilon_t, \qquad H_0 : \rho = 0 \;\;(\text{no cointegration})\]

Under \(H_0\), \(\hat z_t\) carries a unit root: the two series drift apart independently. Under \(H_1\), \(\rho < 0\) and the gap is pulled back towards zero. We reject \(H_0\) when the \(t\) statistic attached to \(\hat\rho\) is sufficiently negative.

The technical point that matters. The statistic is not Student\'s, and it is not even that of an ordinary ADF. Two reasons compound. First, under a unit root the limiting law is expressed through functionals of Brownian motion, not a normal — that is already true of Dickey–Fuller. Second, \(\hat z_t\) is not observed: it is a residual whose variance least squares has minimised, which makes it look more stationary than it is. The critical values must therefore be shifted to the left (MacKinnon, 1991).

Test	Statistic	5% critical value
Ordinary t test (invalid here)	\(t\)	−1.96
Augmented Dickey–Fuller, observed series, with constant	\(t_{\rho}\)	≈ −2.86
Engle–Granger, residual, 2 variables	\(t_{\rho}\)	≈ −3.34
Engle–Granger, residual, 3 variables	\(t_{\rho}\)	≈ −3.74
Engle–Granger, residual, 4 variables	\(t_{\rho}\)	≈ −4.10

Approximate asymptotic values from MacKinnon (1991) for the case with a constant and no deterministic trend. Using −1.96 multiplies the false discovery rate by five or six.

Checking all this by simulation

The code below really runs in your browser (Python compiled to WebAssembly). Edit the values and re-run: the first launch takes a few seconds while Python loads.

PYTHON · NUMPY

# Experiment 2: the two-step Engle-Granger test, applied to a
# cointegrated pair AND THEN to a pair that is not cointegrated.
import numpy as np

CV_5 = -3.34          # 5% critical value, 2 variables, with constant

def engle_granger(x1, x2, p=1):
    b, zhat, t_stat, r2 = mco(x1, x2)        # step 1: long-run relationship
    t_adf = adf_t(zhat, p=p, const=False)    # step 2: ADF on the residual
    return b, zhat, t_adf, r2

# --- case A: genuine cointegration (alpha = 0.15) ---
x1, x2, z = cointegrated_pair(n=250, beta=1.4, alpha=0.15, seed=11)
b, zhat, t_adf, r2 = engle_granger(x1, x2)
print("CASE A - cointegrated pair")
print(f"  estimated beta : {b[1]:6.3f}   (true value 1.400)   R2 = {r2:.3f}")
print(f"  t_ADF = {t_adf:6.2f}   vs {CV_5}  ->",
      "COINTEGRATION" if t_adf < CV_5 else "not rejected")

# --- case B: two independent random walks ---
rng = np.random.default_rng(159)
y1 = np.cumsum(rng.normal(0, 1, 250))
y2 = np.cumsum(rng.normal(0, 1, 250))
b2, zhat2, t_adf2, r2b = engle_granger(y1, y2)
print("\nCASE B - independent series")
print(f"  estimated beta : {b2[1]:6.3f}   R2 = {r2b:.3f}  (high, and yet meaningless)")
print(f"  t_ADF = {t_adf2:6.2f}   vs {CV_5}  ->",
      "COINTEGRATION" if t_adf2 < CV_5 else "not rejected")

plot_pair(x1, x2, "Case A: cointegrated pair")
plot_resid(zhat, zhat2, "Residuals: case A (stationary) vs case B (random walk)")

Error correction: the representation theorem

The deepest result in the theory is not the test but an equivalence. The Granger representation theorem states that two \(I(1)\) series are cointegrated if and only if they admit an error correction representation:

Error correction model (ECM), bivariate case

\[\Delta X_{t1} = \alpha_1 \underbrace{(X_{t-1,1} - \beta X_{t-1,2})}_{z_{t-1}\;=\;\text{disequilibrium}} + \sum_{i} \gamma_{1i}\Delta X_{t-i,1} + \sum_{i} \delta_{1i}\Delta X_{t-i,2} + a_{t1}\] \[\Delta X_{t2} = \alpha_2 \, z_{t-1} + \sum_{i} \gamma_{2i}\Delta X_{t-i,1} + \sum_{i} \delta_{2i}\Delta X_{t-i,2} + a_{t2}\]

Each equation is perfectly legitimate in the least squares sense: every variable is stationary, since we work in differences and \(z_{t-1}\) is \(I(0)\) by assumption. The coefficient \(\alpha_1\) measures the fraction of the disequilibrium corrected each period; it must be negative for the system to be stable. At least one \(\alpha_j\) must be non-zero, otherwise nothing pulls the system back to equilibrium.

In VAR notation the reformulation is illuminating. A VAR(\(p\)) in levels \(\mathbf{X}_t = \Phi_1\mathbf{X}_{t-1} + \dots + \Phi_p \mathbf{X}_{t-p} + \mathbf{a}_t\) can be rewritten exactly as:

VECM form of a VAR(p)

\[\Delta \mathbf{X}_t = \Pi \mathbf{X}_{t-1} + \sum_{i=1}^{p-1} \Gamma_i \Delta\mathbf{X}_{t-i} + \mathbf{a}_t, \qquad \Pi = -\Big(I - \sum_{k=1}^{p}\Phi_k\Big)\]

Everything then hinges on the rank of the \(m \times m\) matrix \(\Pi\):

rank(\(\Pi\)) = r	Interpretation	What to fit
\(r = 0\)	\(\Pi = 0\): no long-run relationship	VAR in first differences
\(0 < r < m\)	\(\Pi = \boldsymbol\alpha \boldsymbol\beta\'\), \(r\) cointegrating relations	VECM of rank \(r\)
\(r = m\)	\(\Pi\) has full rank: \(\mathbf{X}_t\) was stationary all along	VAR in levels

The factorisation \(\Pi = \boldsymbol\alpha\boldsymbol\beta\'\) cleanly separates the two natures: \(\boldsymbol\beta\) (of size \(m \times r\)) holds the long-run equilibrium relations, \(\boldsymbol\alpha\) (of size \(m \times r\)) the short-run adjustment speeds.

PYTHON · NUMPY

# Experiment 3: estimating the error correction model in 2 steps.
# Step 1 is the OLS long run; we now estimate the adjustment speed alpha.
import numpy as np

ALPHA_TRUE = 0.12
x1, x2, z = cointegrated_pair(n=400, beta=0.8, alpha=ALPHA_TRUE, seed=12)

# step 1: long-run relationship
b, zhat, t_stat, r2 = mco(x1, x2)
print(f"estimated beta = {b[1]:.4f}   (true 0.8000)")
print(f"t_ADF on the residual = {adf_t(zhat, const=False):.2f}   (critical value -3.34)")

# step 2: ECM for the first equation
#   Delta x1_t = a * zhat_{t-1} + g * Delta x2_t + c + noise
dx1 = np.diff(x1)
dx2 = np.diff(x2)
X   = np.column_stack([zhat[:-1], dx2, np.ones(len(dx1))])
coef, *_ = np.linalg.lstsq(X, dx1, rcond=None)
a1 = coef[0]

half_life = np.log(0.5) / np.log(1 + a1) if -1 < a1 < 0 else np.inf
print(f"estimated alpha_1 = {a1:+.4f}   (true {-ALPHA_TRUE:+.4f})")
print(f"half-life of the disequilibrium: {half_life:.1f} periods")

# superconsistency: the error on beta decays as 1/n, not 1/sqrt(n)
print("\nSuperconsistency of beta:")
for n in [100, 400, 1600, 6400]:
    e = [abs(mco(*cointegrated_pair(n=n, beta=0.8, alpha=ALPHA_TRUE, seed=s)[:2])[0][1] - 0.8)
         for s in range(40)]
    print(f"  n = {n:5d}   mean error |beta_hat - beta| = {np.mean(e):.5f}")

plot_ecm(zhat, dx1, "Disequilibrium z(t-1) and correction Delta x1(t)")

Beyond two series: Johansen

Engle–Granger has two weaknesses as soon as you go past two variables. First, the result depends on which variable you put on the left of the equals sign: regressing \(X_1\) on \(X_2\) or \(X_2\) on \(X_1\) does not give exactly the same verdict in finite samples. Second, with \(m\) series there can be up to \(m-1\) independent cointegrating relations, and a single regression reveals at best one combination of them.

The Johansen procedure (1988, 1991) attacks the problem directly as a question of rank. The VECM is estimated by maximum likelihood and the rank of \(\Pi\) is tested from the eigenvalues \(\hat\lambda_1 \geq \dots \geq \hat\lambda_m\) of a canonical correlation problem between \(\Delta\mathbf{X}_t\) and \(\mathbf{X}_{t-1}\), corrected for the lags:

Trace statistic

\[\mathrm{LR}_{\text{trace}}(r) = -n \sum_{j=r+1}^{m} \ln\!\big(1 - \hat\lambda_j\big), \qquad H_0 : \text{rank}(\Pi) \leq r\]

The test is run sequentially for \(r = 0\), then \(r = 1\), and so on, up to the first non-rejection — which gives the retained rank. The maximum eigenvalue statistic \(-n\ln(1-\hat\lambda_{r+1})\) provides a complementary test. Here too the limiting laws are non-standard and tabulated separately according to the deterministic terms included (constant inside the cointegrating relation, trend, and so on) — six classic cases, and picking the wrong one is a common source of error.

In practice (R and Python)

Nobody codes a Johansen test by hand in production. Here are the standard paths, with the detail that matters on each line.

# ─── R: tseries + urca + vars ───
library(tseries); library(urca); library(vars)

# 0. check that EVERY series really is I(1) first
adf.test(x1)              # expect no rejection in levels
adf.test(diff(x1))        # expect rejection in differences

# 1. Engle-Granger, two steps
eg  <- lm(x1 ~ x2)
adf.test(residuals(eg))   # WARNING: the p-value is invalid here
po.test(cbind(x1, x2))    # Phillips-Ouliaris: correct critical values

# 2. Johansen: cointegration rank
jo <- ca.jo(cbind(x1, x2), type = "trace", ecdet = "const", K = 2)
summary(jo)               # compare the trace stat with the critical values

# 3. fitted VECM, then conversion to a VAR in levels for forecasting
vecm <- cajorls(jo, r = 1)
var_levels <- vec2var(jo, r = 1)
predict(var_levels, n.ahead = 25)

# ─── Python: statsmodels ───
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM

adfuller(x1)                       # I(1)? in levels, then in differences
coint(x1, x2)                      # Engle-Granger, correct critical values

jres = coint_johansen(X, det_order=0, k_ar_diff=1)
print(jres.lr1, jres.cvt)          # trace vs 90/95/99% critical values

model = VECM(X, k_ar_diff=1, coint_rank=1, deterministic="ci").fit()
print(model.beta, model.alpha)     # long run, adjustment speeds
model.predict(steps=25)

The adf.test(residuals(eg)) trap. The function will return a p-value, and that p-value will be wrong — too small — because it is computed for an observed series, not for a least squares residual. Use po.test() in R or coint() in Python, which apply MacKinnon critical values adapted to the number of regressors.

Limits and pitfalls

Assumption or pitfall	Consequence → remedy
Series must share the same order \(I(1)\)	Mixing \(I(0)\) and \(I(1)\) makes the test uninterpretable → test each series separately first
Structural break in the relationship	The test fails to reject \(H_0\) although cointegration holds piecewise → Gregory–Hansen tests
Low power in small samples	With \(n < 100\) and slow adjustment almost nothing is detected → lengthen the sample rather than raising the frequency
Choice of the number of lags \(p\) in the ADF	Too few: autocorrelated residuals; too many: loss of power → AIC or the Schwert rule
Deterministic terms (constant, trend)	Each case has its own critical values → decide before looking at the data
Cointegration ≠ causality	An equilibrium relation does not say who adjusts to whom → read the \(\alpha_j\) and test weak exogeneity
Relationship stable over time	In finance, pairs "de-cointegrate" → re-estimate on a rolling window and monitor the half-life

The most expensive mistake. Differencing everything to "make it stationary" when the series are in fact cointegrated. You get a valid but amputated model: all the long-run information, the part that holds the adjustment mechanism, disappears in the differencing. Long horizon forecasts then become clearly worse than those of the corresponding VECM.

References

Engle, R. F. & Granger, C. W. J. (1987). Co-integration and error correction: representation, estimation, and testing, Econometrica 55(2) — the founding paper, Nobel 2003.
Granger, C. W. J. & Newbold, P. (1974). Spurious regressions in econometrics, Journal of Econometrics 2 — the original warning, by simulation.
Phillips, P. C. B. (1986). Understanding spurious regressions in econometrics, Journal of Econometrics 33 — the asymptotic theory behind it.
Johansen, S. (1991). Estimation and hypothesis testing of cointegration vectors in Gaussian vector autoregressive models, Econometrica 59(6).
MacKinnon, J. G. (1991). Critical values for cointegration tests, in Engle & Granger (eds.), Long-Run Economic Relationships, Oxford University Press.
Woodward, W. A., Gray, H. L. & Elliott, A. C. (2017). Applied Time Series Analysis with R, 2nd ed., CRC Press — chapter 10, the VAR/VARMA framework and tests of association between series.
Murray, M. P. (1994). A drunk and her dog: an illustration of cointegration and error correction, The American Statistician 48(1) — the leash metaphor.

The simulators on this page are written in plain JavaScript/Canvas with no external dependency; the Python cells run through Pyodide. The source is available on GitHub.