Cointegration: Engle-Granger, Johansen and the half-life of mean reversion.

Correlation is not cointegration

A drunk leaving a bar walks aimlessly: each step is independent of the previous one and, over time, he can end up anywhere. If he has a dog off the leash, the dog wanders too. But if the dog is on a long leash, the two wander aimlessly and, even so, never get further apart than the leash. Murray (1994) used this image to explain cointegration: two series that are random walks on their own but whose distance is stationary.

In trading, the two walkers are two prices. The question that matters for a pair (lesson 2.3.3) is not whether they go up and down together, but whether the distance between them returns to a level. They are different questions:

• Correlation of returns measures whether the two prices move in the same direction on the same day. It is a short-term property. Two random walks with a common factor have a high correlation and can drift apart without limit.
• Cointegration says that there is a $\beta$ such that the spread $\ln y_t - \beta \ln x_t$ is stationary: it returns to its mean at a finite speed. It is a long-term property, and it does not require daily returns to be highly correlated.

Formally: two series integrated of order 1 (their differences are stationary, lesson 2.1.2) are cointegrated when there is a vector $(1, -\beta)$ whose combination is stationary. With prices one usually works with their logarithms, because what returns is a ratio and the spread is expressed as a percentage:

flowchart TB
  P["I(1) prices"] --> R{Correlated}
  P --> S{"I(0) spread"}
  R -- yes --> C1[Short-term hedge]
  R -- no --> N1[Little relationship]
  S -- yes --> C2[Mean-reverting pair]
  S -- no --> N2[They drift apart]

The two branches are independent: a pair can be in any of the four combinations. The figure at the top shows the difference with data: above, the spread of pair B (cointegrated) crosses zero constantly; below, the residual of control E (more correlated than B, not cointegrated) spends whole years on one side of zero.

The eight pairs of the course

The data for this lesson are eight synthetic pairs (pares2.parquet), NYSE sessions from 2014 to 2025, with the known answer in pares2_verdad.csv. All the relationships were planted between logarithms. Three pairs are cointegrated with different half-lives (A, B and C); D has a β that changes over time (lesson 2.1.6); E, F and G are controls: their returns are highly correlated, but the relationship between their logarithms carries a random walk, so they are not cointegrated; H is cointegrated until June 15, 2021 and stops being so from that day on. They are made-up data with the effect planted: they serve to check whether a method recovers what is known to be there, not as evidence that a real pair works.

import numpy as np
import pandas as pd

pares = pd.read_parquet("pares2.parquet")
verdad = pd.read_csv("pares2_verdad.csv", index_col="par")
precios = pares.pivot(index="fecha", columns="activo",
                      values="precio")
lp = np.log(precios)      # the relationships are between logarithms
print(precios.shape, precios.index[0].date(),
      precios.index[-1].date())
verdad[["tipo", "formulacion", "beta", "vida_media_sesiones",
        "sigma_spread"]]

(3018, 16) 2014-01-02 2025-12-31
                         tipo formulacion  beta  vida_media_sesiones  sigma_spread
par
A                 cointegrado         log   1.5                  5.0         0.020
B                 cointegrado         log   0.8                 15.0         0.030
C                 cointegrado         log   1.1                 40.0         0.045
D               beta_variable         log   NaN                 10.0         0.025
E    correlado_no_cointegrado         log   1.0                  NaN           NaN
F    correlado_no_cointegrado         log   0.9                  NaN           NaN
G    correlado_no_cointegrado         log   1.2                  NaN           NaN
H                     ruptura         log   0.9                 10.0         0.025

Each pair has an asset x and an asset y (PZAX and PZAY for A). The formulacion column says the relationship is between logarithms (log), and sigma_spread is the standard deviation of the log spread: 0.02 is a typical 2% distance between the two prices. The first thing is to measure what an analyst in a hurry usually measures, correlation, next to the Engle-Granger test explained in the following section:

from statsmodels.tsa.stattools import coint

filas = []
for par, v in verdad.iterrows():
    y, x = lp[v.activo_y], lp[v.activo_x]
    t, p, _ = coint(y, x)
    filas.append({"par": par, "tipo": v.tipo,
                  "corr_rend": y.diff().corr(x.diff()),
                  "corr_logprecios": y.corr(x),
                  "t_eg": t, "p_eg": p})
tabla = pd.DataFrame(filas).set_index("par")
tabla.round(3)

                         tipo  corr_rend  corr_logprecios    t_eg   p_eg
par
A                 cointegrado      0.912            1.000 -14.971  0.000
B                 cointegrado      0.792            0.990  -8.979  0.000
C                 cointegrado      0.875            0.992  -5.586  0.000
D               beta_variable      0.871            0.994  -3.472  0.035
E    correlado_no_cointegrado      0.859            0.978  -2.599  0.237
F    correlado_no_cointegrado      0.856            0.976  -2.036  0.510
G    correlado_no_cointegrado      0.856            0.931  -2.371  0.339
H                     ruptura      0.635            0.941  -2.647  0.219

The three controls have a correlation of returns of 0.86, more than pair B (0.79), and the correlation of their logarithms is above 0.93. If the criterion were correlation, E would be a better candidate than B. The test says the opposite: B rejects the hypothesis of “no cointegration” comfortably (t = −8.98, far below the 5% critical value of −3.34) and E, F and G do not reject it (p between 0.24 and 0.51). The correlation of prices is an even worse guide: two trending series have correlations of levels close to 1 even if they have nothing to do with each other (spurious regression).

D and H are the trick cases. H does not reject in the full sample (p = 0.22): the relationship existed and broke, and the jump of the break drags down even its correlation of returns (0.64). D does reject (p = 0.035), and it is a misleading pass: its true β moves and the full-sample regression estimates an average one; the residual of that fixed β reverts, but much more slowly than the true spread (it shows in the half-life, below). Both are dealt with in this lesson and in 2.1.6.

Engle-Granger in two steps

The Engle and Granger (1987) test is the most direct way to check the definition:

1. Estimate $\ln y_t = \alpha + \beta \ln x_t + u_t$ by least squares.
2. Test whether the residual $\hat u_t$ has a unit root with an augmented Dickey-Fuller regression, $\Delta \hat u_t = \gamma \hat u_{t-1} + \sum_{i=1}^{k} c_i \Delta \hat u_{t-i} + e_t$. The null hypothesis is $\gamma = 0$: the residual is a random walk and there is no cointegration.

flowchart TB
  A["ln y, ln x"] --> B[OLS: y on x]
  B --> C[Residual u]
  C --> D[ADF of u]
  D --> E{t below critical}
  E -- yes --> F[Cointegrated]
  E -- no --> G[No evidence]

The detail that gets forgotten is the critical value. In the first step, OLS chooses the β that minimizes the variance of the residual, that is, the one that makes it look most stationary. That is why the ADF of the residual cannot be compared with the ADF tables for an observed series: more negative critical values are needed, MacKinnon's for cointegration residuals, which are the ones coint uses. Doing it by hand shows it:

import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

def dos_pasos(y, x):
    mco = sm.OLS(y, sm.add_constant(x)).fit()
    adf = adfuller(mco.resid, result_object=True)
    eg = coint(y, x)
    return {"alfa": mco.params.iloc[0], "beta": mco.params.iloc[1],
            "t": adf.statistic, "p_adf": adf.pvalue,
            "crit5_adf": adf.critical_values["5%"],
            "p_coint": eg[1], "crit5_coint": eg[2][1]}

pd.DataFrame({par: dos_pasos(lp[f"PZ{par}Y"], lp[f"PZ{par}X"])
              for par in ["B", "E"]}).round(3)

                 B      E
alfa         0.683  0.913
beta         0.805  0.821
t           -8.978 -2.598
p_adf        0.000  0.093
crit5_adf   -2.862 -2.862
p_coint      0.000  0.237
crit5_coint -3.338 -3.338

In B, OLS recovers the planted relationship ($\alpha$ = 0.7 and $\beta$ = 0.8; it gets 0.68 and 0.805) and the residual is clearly stationary with any table. In E the statistic is the same in both columns (−2.60), but the plain ADF gives p = 0.09, a rejection at 10%, and coint, p = 0.24: the 5% critical value goes from −2.86 to −3.34. In a screen of many pairs, using the wrong table multiplies the false positives.

Engle-Granger has two known weaknesses. The first is that it is not symmetric: regressing y on x or x on y gives different statistics and p-values, and in a borderline case each direction can lead to a different conclusion. The second is that with more than two series it only finds one relationship, even if there are several. Johansen solves both.

Johansen: how many relationships there are among several series

Johansen (1991) starts from a vector error correction model (VECM) for a vector $Y_t$ of $n$ series, here the logarithms of the prices:

If the series are not cointegrated, $\Pi = 0$ and the model is a VAR in differences. If there are $r$ cointegrating relationships, $\Pi = \alpha \beta'$ has rank $r$: the $r$ columns of $\beta$ are the stationary combinations and $\alpha$ says how strongly each series returns to equilibrium. The trace test asks in sequence “rank ≤ 0?”, “rank ≤ 1?”… and the estimated rank is the first one that is not rejected.

With the four series of A and B, the known answer is $r = 2$:

import warnings
from statsmodels.tsa.vector_ar.vecm import coint_johansen

# With numpy 2.5, the eigenvectors come out as complex numbers with a
# zero imaginary part and statsmodels warns (ComplexWarning) when
# casting them to real: the warning does not change any number
warnings.filterwarnings("ignore",
                        category=np.exceptions.ComplexWarning)

def johansen(datos, k_ar_diff=1):
    j = coint_johansen(datos, det_order=0, k_ar_diff=k_ar_diff)
    n = datos.shape[1]
    tabla = pd.DataFrame(
        {"traza": j.lr1, "crit95_traza": j.cvt[:, 1],
         "max_autov": j.lr2, "crit95_max": j.cvm[:, 1]},
        index=[f"rank <= {r}" for r in range(n)])
    rango = int((j.lr1 > j.cvt[:, 1]).cumprod().sum())
    return tabla, rango, np.real(j.evec)

cuatro = lp[["PZAY", "PZAX", "PZBY", "PZBX"]]
tabla_j, rango, vectores = johansen(cuatro)
print("Estimated rank:", rango)
tabla_j.round(1)

Estimated rank: 2
           traza  crit95_traza  max_autov  crit95_max
rank <= 0  296.2          47.9      209.6        27.6
rank <= 1   86.6          29.8       80.0        21.1
rank <= 2    6.6          15.5        5.6        14.3
rank <= 3    0.9           3.8        0.9         3.8

The trace rejects “rank ≤ 0” (296.2 versus 47.9) and “rank ≤ 1” (86.6 versus 29.8), and does not reject “rank ≤ 2” (6.6 versus 15.5): two relationships, the planted ones. The cointegrating vectors are not unique (any combination of two stationary vectors is stationary too), but once normalized they are recognizable:

v1 = vectores[:, 0] / vectores[0, 0]   # normalized by PZAY
v2 = vectores[:, 1] / vectores[2, 1]   # normalized by PZBY
pd.DataFrame({"vector_1": v1, "vector_2": v2},
             index=cuatro.columns).round(3)

        vector_1  vector_2
activo
PZAY       1.000    -0.040
PZAX      -1.501     0.065
PZBY       0.070     1.000
PZBX      -0.059    -0.812

The first vector is almost exactly the relationship of A, $\ln y - 1.5 \ln x$ (−1.501), with small weights on B; the second, almost that of B, $\ln y - 0.8 \ln x$ (−0.812), with small weights on A. Here they come out clean because the two relationships have very different reversion strengths, but that is not guaranteed: any combination of two stationary vectors is stationary too, and what Johansen identifies uniquely is the number of relationships and the space they span, not a particular basis.

If a series that has nothing to do with a cointegrated pair is added to it, Johansen should find a single relationship and give zero weight to the unrelated series. And with series that are not cointegrated with each other it should give rank 0, except by chance (at 5%, in one case in twenty). This last check depends on the formulation:

from itertools import combinations

tabla_3, rango_3, vec_3 = johansen(lp[["PZAY", "PZAX", "PZEX"]])
print("A + PZEX: rank", rango_3, "vector",
      (vec_3[:, 0] / vec_3[0, 0]).round(3))

controles = ["PZEY", "PZEX", "PZFY", "PZFX", "PZGY", "PZGX"]
trios = list(combinations(controles, 3))
falsos = {"levels": 0, "logs": 0}
for trio in trios:
    datos = precios[list(trio)]
    falsos["levels"] += johansen(datos)[1] > 0
    falsos["logs"] += johansen(lp[list(trio)])[1] > 0
print(len(trios), "control trios; with rank > 0:", falsos)
crec = precios[controles].iloc[-1] / precios[controles].iloc[0]
print("Final / initial price:", crec.round(1).to_dict())

A + PZEX: rank 1 vector [ 1.    -1.51   0.008]
20 control trios; with rank > 0: {'levels': 4, 'logs': 0}
Final / initial price: {'PZEY': 2.4, 'PZEX': 3.1, 'PZFY': 3.3, 'PZFX': 3.4, 'PZGY': 5.1, 'PZGX': 3.2}

With the unrelated series, rank 1 and an almost zero weight (0.008) for PZEX: correct. With the controls, in logs none of the 20 trios gives a rank above zero, as expected (chance would give one); in levels, 4. The reason is in the last line: the controls multiply by between 2.4 and 5.1 in twelve years. The test assumes errors with constant variance, and in levels the daily change of a price is proportional to its level: a price that has multiplied by 5 has changes five times larger at the end than at the start, and the distribution of the trace is no longer that of the tables. In logs the daily change is a percentage and its size does not depend on the level.

Half-life: the spread as an Ornstein-Uhlenbeck process

Knowing that a spread reverts is not enough to trade it: you need to know how fast. The standard model is the Ornstein-Uhlenbeck process, a random walk with a force that pulls toward the mean:

$\theta$ is the speed of reversion: the expected distance to the mean decays like $e^{-\theta t}$. Sampled once per session, it is an AR(1), $S_t = c + \varphi S_{t-1} + \varepsilon_t$ with $\varphi = e^{-\theta}$, and the half-life (the time in which the expected distance to the mean is halved) is

It is estimated with the regression of the change in the spread on its previous level, $\Delta S_t = c + b\, S_{t-1} + e_t$, with $\varphi = 1 + b$. In the cointegrated pairs it is applied to the Engle-Granger residual:

def vida_media(spread):
    """Half-life in sessions of the AR(1) dS_t = c + b * S_(t-1)."""
    s = pd.Series(spread).dropna()
    X = sm.add_constant(s.shift(1).iloc[1:])
    b = sm.OLS(s.diff().iloc[1:], X).fit().params.iloc[1]
    return -np.log(2) / np.log(1 + b)

def residuo_eg(par):
    y, x = lp[f"PZ{par}Y"], lp[f"PZ{par}X"]
    return sm.OLS(y, sm.add_constant(x)).fit().resid

vidas = pd.DataFrame({
    "estimada": {p: vida_media(residuo_eg(p)) for p in "ABCDEFG"},
    "verdadera": verdad["vida_media_sesiones"]}).dropna(
        subset=["estimada"])
vidas.round(1)

   estimada  verdadera
A       4.7        5.0
B      13.0       15.0
C      35.4       40.0
D      38.4       10.0
E     162.1        NaN
F     259.4        NaN
G     270.0        NaN

In A, B and C the estimate is close to the truth but always below it: 4.7 versus 5, 13.0 versus 15 and 35.4 versus 40. In E, F and G the function returns half-lives of 162 to 270 sessions, even though there are none: an AR(1) always gives a $\varphi$ below 1 in a finite sample, so the half-life of a random walk comes out finite. A half-life only makes sense after a cointegration test that supports it. And D shows the opposite: the test passed, but the residual of its fixed β has a half-life of 38.4 sessions when that of the true spread is 10. What reverts in D is something other than what was estimated: with a β that moves, the residual of an average β mixes the spread with the slow drift of β (lesson 2.1.6).

Is it bad luck or bias? It can be answered by simulation: generate many Ornstein-Uhlenbeck processes with the true half-life and the length of this data, and estimate each one the same way.

from scipy.signal import lfilter

rng = np.random.default_rng(2028)
n, caminos = len(lp), 400
filas = {}
for h in [5, 15, 40]:
    phi = 0.5 ** (1 / h)
    e = rng.standard_normal((caminos, n)) * np.sqrt(1 - phi**2)
    s = lfilter([1.0], [1.0, -phi], e, axis=1)
    ds, s0 = np.diff(s, axis=1), s[:, :-1]
    s0c = s0 - s0.mean(axis=1, keepdims=True)
    b = (ds * s0c).sum(axis=1) / (s0c**2).sum(axis=1)
    est = -np.log(2) / np.log(1 + b)
    filas[h] = np.percentile(est, [5, 50, 95])
pd.DataFrame(filas, index=["p5", "mediana", "p95"]).T.round(1)

      p5  mediana   p95
5    4.4      4.9   5.6
15  11.8     14.7  17.8
40  27.3     37.2  51.2

The median of the simulations is below the truth, and more so the slower the reversion: it is the finite-sample bias of the autoregressive coefficient, which pushes $\hat\varphi$ down. The three estimates from the data fall inside the 90% band, so they are what one would expect. The band is what matters most: with twelve years of data, the half-life of 40 sessions could have been estimated at 27 or at 51. A pair with a half-life of 40 sessions has barely 75 half-lives in the sample, few independent observations of the reversion.

The half-life is used later for two decisions of the system in lesson 2.3.3: how long to hold a position before calling it failed and which window to use to estimate the mean and deviation of the spread. A badly estimated half-life spoils both.

When cointegration breaks: pair H

Pair H is cointegrated (half-life 10, spread with a standard deviation of 0.025 in logs) until June 15, 2021. That day the spread jumps five standard deviations and from then on it is a random walk. It is the most dangerous real-world situation: a relationship that worked for years, that any test confirms with past data, and that ceases to exist.

y, x = lp["PZHY"], lp["PZHX"]
antes = lp.index < "2021-06-15"
for nombre, m in [("full sample", slice(None)),
                  ("before the break", antes),
                  ("after the break", ~antes)]:
    t, p, _ = coint(y[m], x[m])
    print(f"{nombre:22} sessions {len(y[m]):5}  t {t:6.2f}  "
          f"p {p:.4f}")

full sample            sessions  3018  t  -2.65  p 0.2191
before the break       sessions  1875  t  -5.91  p 0.0000
after the break        sessions  1143  t  -3.66  p 0.0205

Before the break the cointegration is unambiguous (t = −5.91; p < 0.0001). The full sample no longer detects it (p = 0.22). And the later stretch, which is a random walk by construction, gives t = −3.66 and p = 0.021: a false positive. An analyst who looked only at the last four years would conclude that the relationship is still alive. A 5% test rejects one time in twenty when there is nothing; if you look at a single stretch, chance can strike.

One way to see the whole problem is to test in rolling one-year windows, as a system that re-estimates the relationship every month would:

def ventanas(y, x, n=252, paso=21):
    return np.array([coint(y.iloc[i:i + n], x.iloc[i:i + n],
                           maxlag=10)[1]
                     for i in range(0, len(y) - n + 1, paso)])

casos = {"A (half-life 5)": "A", "C (half-life 40)": "C",
         "E (control)": "E", "F (control)": "F",
         "G (control)": "G"}
filas = {}
for nombre, par in casos.items():
    p = ventanas(lp[f"PZ{par}Y"], lp[f"PZ{par}X"])
    filas[nombre] = (len(p), (p < 0.05).mean())
p = ventanas(y[~antes], x[~antes])
filas["H after the break"] = (len(p), (p < 0.05).mean())
pd.DataFrame(filas, index=["ventanas", "rechazo_5pct"]).T.round(3)

                   ventanas  rechazo_5pct
A (half-life 5)       132.0         0.977
C (half-life 40)      132.0         0.091
E (control)           132.0         0.038
F (control)           132.0         0.098
G (control)           132.0         0.061
H after the break      43.0         0.070

Three readings. In A, with a half-life of 5 sessions, one year of data is enough and the test rejects in 98% of the windows. In C, truly cointegrated but with a half-life of 40, one year contains barely six half-lives and the test rejects in 9%, almost as if there were nothing: it is a problem of power, not of the absence of a relationship. In the controls the rejection rate ranges from 3.8% to 9.8%, around the 5% expected from chance (the windows overlap, so those rates are very noisy), and in H after the break it is 7% (3 of 43 windows). The short window does not tell a slow pair from a non-existent one well; the long one mixes stretches with and without a relationship. There is no right window: you have to know what is lost with each one. Detecting the break itself is the subject of lesson 2.1.7.

Screening pairs: selection bias

In practice nobody tests one pair: hundreds are screened and the ones that pass are kept. With the sixteen assets in the data there are 120 possible pairings, and only those of A, B and C (and H for part of the sample) have a relationship. The 112 pairings of assets from different pairs are not cointegrated by construction:

filas = []
for a, b in combinations(precios.columns, 2):
    filas.append({
        "mismo_par": a[:3] == b[:3],
        "p_logs": coint(lp[a], lp[b], maxlag=10)[1],
        "p_niveles": coint(precios[a], precios[b], maxlag=10)[1]})
criba = pd.DataFrame(filas)
ajenas = criba[~criba["mismo_par"]]
print("Pairings from different pairs:", len(ajenas))
(ajenas[["p_logs", "p_niveles"]] < 0.05).sum()

Pairings from different pairs: 112
p_logs       4
p_niveles    2
dtype: int64

In logs, 4 of the 112 unrelated pairings pass the test at 5%, and in levels, 2: both numbers are close to the 5.6 that chance would give, and there is no lesson about the formulation between them (the pairings share a market factor and are not independent tests, so the number of rejections varies a lot from one sample to another). What matters is that they exist: any of those four pairings, presented on its own, would look like a finding. The problem of picking the best of many tests has a name and corrections (lesson 2.1.8); the simplest defense is to require an economic reason before looking at the p-value and to validate outside the sample in which the pair was chosen (lesson 2.5.2).

Common mistakes with cointegration

• Estimating β on the full sample and trading on it. The β and the mean of the spread in a backtest have to be estimated with data prior to each trade; with the full sample, the spread reverts “by construction” (lesson 2.5.1).
• Screening without correcting. Of 112 unrelated pairings, 4 pass in logs. With hundreds of candidates, the “pairs” that pass are largely chance.
• Confusing correlation and cointegration. The controls E, F and G have more correlation than B and do not revert.
• Trusting a pass with a fixed β. D passes Engle-Granger with an average β, but the residual that reverts has a half-life four times longer than the real spread.
• Choosing the wrong formulation. In levels, Johansen gives false positives and Engle-Granger misses pair A.
• Ignoring power. A slow pair looks non-cointegrated in short windows; a fast pair looks cointegrated in almost any window.
• Assuming it lasts. The relationship of H worked for seven years and broke in a day. Monitoring the spread and the stability of β (lessons 2.1.6 and 2.1.7) is part of the system, not an extra.
• Unadjusted prices. An unadjusted split or dividend puts a jump in one of the two prices and breaks a cointegration that exists (lesson 1.4.3).

Exercises

1. Repeat the Engle-Granger test for pairs A, B, C and for the controls E, F and G with prices in levels instead of in logarithms. Does any conclusion change?

for par in "ABCEFG":
    t, p, _ = coint(precios[f"PZ{par}Y"], precios[f"PZ{par}X"])
    print(par, f"t {t:6.2f}  p {p:.3f}")

A t  -3.30  p 0.055
B t  -8.56  p 0.000
C t  -5.44  p 0.000
E t  -2.32  p 0.366
F t  -1.85  p 0.607
G t  -2.04  p 0.508

2. Estimate the half-life of pair B with the first three years (2014-2016) and with the last three (2023-2025), using in each stretch its own β. How much does it vary? Does it fall inside the simulated band for a half-life of 15?

yb, xb = lp["PZBY"], lp["PZBX"]
for desde, hasta in [("2014", "2016"), ("2023", "2025")]:
    ys, xs = yb.loc[desde:hasta], xb.loc[desde:hasta]
    u = sm.OLS(ys, sm.add_constant(xs)).fit().resid
    print(desde, hasta, f"sessions {len(u)}",
          f"half-life {vida_media(u):.1f}")

2014 2016 sessions 756 half-life 9.7
2023 2025 sessions 752 half-life 13.5

3. A colleague screens 300 pairings of stocks from the same sector in levels, finds 24 with p < 0.05 in Engle-Granger and proposes trading all of them. With what you have seen in this lesson, what would you tell them?