Build a Pairs Trade with Cointegration

Two stocks in the same business, Coke and Pepsi, Shell and BP, are pushed around by the same macro forces, so their prices tend to wander together. A pairs trade bets that when they temporarily pull apart, the gap closes. The whole strategy lives or dies on one question: is that gap actually mean-reverting, or does it just look tight because both names drifted up together? Answering it correctly is the difference between statistical arbitrage and slowly bleeding out.

This guide builds the machinery in pure numpy/pandas, no statsmodels, so you can see every number. We reuse the shift-clean backtest discipline from backtest from scratch.

Why cointegration, not correlation

Correlation is about co-movement of returns over a short window; it says nothing about whether the price gap stays bounded. Two random walks that happen to drift upward together can post a 0.9 return correlation for years while the dollar spread between them grows without limit, trade that spread expecting reversion and you get run over.

Cointegration is the stronger, structural property you actually need: two individually non-stationary price series are cointegrated if some linear combination of them is stationary, it has a fixed mean and finite variance, so it must revert. That linear combination is your spread. Correlation is a statement about wiggle; cointegration is a statement about a leash. You want the leash.

A genuinely cointegrated pair

To test the code we build a synthetic pair where the answer is known. X is a pure random walk. Y shares that same stochastic trend (2·X) plus a stationary AR(1) disturbance, so Y − 2X mean-reverts by construction. Swap the synthetic block for yf.download([...])["Close"] on real data.

import numpy as np
import pandas as pd

rng = np.random.default_rng(7)
n = 1500
idx = pd.bdate_range("2015-01-01", periods=n)

x = pd.Series(50 + np.cumsum(rng.normal(0, 1.0, n)), index=idx, name="X")

phi = 0.94                                    # AR(1) persistence of the disturbance
u = np.zeros(n)
for t in range(1, n):
    u[t] = phi * u[t-1] + rng.normal(0, 1.0)  # stationary, mean-reverting noise
y = pd.Series(5.0 + 2.0 * x.values + u, index=idx, name="Y")

Step 1, the hedge ratio (Engle-Granger, part 1)

How many units of X hedge one unit of Y? Run an OLS regression of Y on X. The slope $\beta$ is the hedge ratio; the residual is the raw spread.

$Y_t = \alpha + \beta X_t + u_t$

X = np.column_stack([np.ones(n), x.values])          # design matrix: [1, X]
(alpha_hat, beta_hat), *_ = np.linalg.lstsq(X, y.values, rcond=None)
# recovers alpha_hat = 5.05, beta_hat = 2.006  (true 5.0 and 2.0)

OLS nails it: beta_hat = 2.006 against a true value of 2.0. A gotcha worth knowing: OLS is asymmetric, regressing X on Y gives a different slope than $1/\beta$. Practitioners either pick the more "dependent" leg deliberately or use total least squares (orthogonal regression). For a first build, plain OLS is fine.

Step 2, is the spread stationary? (Engle-Granger, part 2)

Now the actual test. Regress the change in the residual spread on its own lagged level, the Augmented Dickey-Fuller regression:

$\Delta s_t = \gamma\, s_{t-1} + \sum_{k=1}^{p}\delta_k \Delta s_{t-k} + \varepsilon_t$

If $\gamma < 0$ significantly, the spread pulls back toward its mean, it is stationary and the pair cointegrates. The test statistic is the t-stat on $\gamma$. Here is the whole thing by hand:

def adf_tstat(s, lags=1):
    s = np.asarray(s, dtype=float)
    ds = np.diff(s)
    rows, yv = [], []
    for t in range(lags, len(ds)):
        row = [1.0, s[t]]                        # [const, s_{t-1}]
        row += [ds[t-k] for k in range(1, lags+1)]  # lagged differences
        rows.append(row); yv.append(ds[t])
    Xm, yv = np.array(rows), np.array(yv)
    b, *_ = np.linalg.lstsq(Xm, yv, rcond=None)
    resid = yv - Xm @ b
    dof = len(yv) - Xm.shape[1]
    cov = (resid @ resid) / dof * np.linalg.inv(Xm.T @ Xm)
    return b[1] / np.sqrt(cov[1, 1])             # t-stat on gamma

The catch that trips everyone up: because we tested the estimated residual (not a raw series), the ordinary Dickey-Fuller critical values are too lenient. Use the Engle-Granger critical values instead, roughly −3.34 at 5% and −3.90 at 1% for a two-variable system. Our spread scores:

spread = y - (alpha_hat + beta_hat * x)
print(adf_tstat(spread, lags=1))                 # -6.47  -> reject unit root, cointegrated
print(adf_tstat(np.cumsum(rng.normal(0,1,n))))   # -2.13  -> a random walk does NOT reject

−6.47 blows past the 1% threshold: the spread is stationary. As a sanity check, a fresh random walk scores −2.13, nowhere near, correctly refusing to call noise stationary. If your control experiment doesn't fail, your test is broken.

Step 3, the spread and its z-score

Standardize the spread into units of its own recent standard deviation. Use a rolling window, not the full-sample mean, full-sample normalization is a subtle lookahead leak (you'd be using future data to center today's z).

$z_t = \frac{s_t - \mu_t^{(w)}}{\sigma_t^{(w)}}$

win = 60
mu, sd = spread.rolling(win).mean(), spread.rolling(win).std()
z = (spread - mu) / sd

Step 4, the half-life sets your clock

How long does reversion take? Fit an AR(1) / discretized Ornstein-Uhlenbeck model to the spread, regress its change on its lagged level, $\Delta s_t = a + b\, s_{t-1} + \varepsilon_t$, and convert the decay coefficient to a half-life:

$\text{half-life} = -\frac{\ln 2}{\ln(1+b)} \approx -\frac{\ln 2}{\theta}$

s = spread.values
ds, slag = np.diff(s), s[:-1]
b = np.linalg.lstsq(np.column_stack([np.ones(len(slag)), slag]), ds, rcond=None)[0][1]
half_life = -np.log(2) / np.log(1 + b)           # 11.5 days

We recover 11.5 days, matching the AR(1) we built in ($\phi=0.94 \Rightarrow \ln 2 / -\ln 0.94 \approx 11.2$). This number is not trivia, it should set your rolling window and expected holding period. A window far shorter than the half-life is pure noise; one far longer smears over the reversion. A 60-day window on an 11-day half-life is a reasonable few-multiples choice.

Step 5, trade it, with no lookahead

Enter when the spread is stretched ($|z|>2$), exit when it snaps back near the mean ($|z|<0.5$). Long the spread means long Y, short $\beta$ units of X. Crucially, decide from data through bar t, hold from t+1, the same shift(1) discipline as every honest backtest.

entry, exit_ = 2.0, 0.5
pos = pd.Series(np.nan, index=z.index)
pos[z >  entry] = -1.0                    # spread rich -> short the spread
pos[z < -entry] =  1.0                    # spread cheap -> long the spread
pos[z.abs() < exit_] = 0.0                # flatten near the mean
pos = pos.ffill().fillna(0.0).shift(1).fillna(0.0)

pnl  = pos * (spread.diff() / sd.shift(1))       # constant-risk spread P&L (z-units)
turn = pos.diff().abs().fillna(0.0)
net  = (pnl - turn * 0.02).fillna(0.0)           # ~2 z-cents per leg crossed: cost

sharpe = lambda r: (r.mean()/r.std())*np.sqrt(252)
print(f"gross {sharpe(pnl):.2f} | net {sharpe(net):.2f} | in-market {(pos!=0).mean():.0%}")
# gross 1.84 | net 1.74 | in-market 26%

A spread is dollar-neutral, so there's no natural percentage denominator, we hold constant spread-risk by dividing the daily spread change by trailing $\sigma$, which also makes the cost charge comparable across regimes. Net Sharpe 1.74, in a position only 26% of the time (flat and earning nothing the rest). That idleness is the point: you only bet when the spread is genuinely stretched.

Where this breaks

A stat-arb Sharpe of 1.74 on synthetic data should make you more suspicious, not less. The two ways real pairs die:

• Cointegration decays. Relationships are estimated on the past and are not guaranteed to hold. Re-test on a rolling window; if the ADF stat drifts above the threshold, the leash is fraying, stop trading the pair. Never fit the hedge ratio once and trust it forever.
• Structural breaks. A merger, a spin-off, a new plant, a regulatory shift, the economic link snaps and the spread walks off to a new level instead of reverting. Watch what it does to the test:

y_break = y.copy()
y_break.iloc[n//2:] += np.linspace(0, 40, n - n//2)   # post-break drift
Xd = np.column_stack([np.ones(n), x.values])
coef = np.linalg.lstsq(Xd, y_break.values, rcond=None)[0]
print(adf_tstat(y_break.values - Xd @ coef))          # -1.36  -> no longer cointegrated

The ADF stat collapses from −6.47 to −1.36: cointegration is gone. Your z-score, still assuming reversion, would keep adding to a losing position as the spread runs, the classic "picking up nickels in front of a steamroller" blow-up. This is why live pairs desks re-estimate constantly and cap losses with a hard stop, not a z-threshold.

The rest of the standard cautions apply too: real fills lag your signal and cross two bid-ask spreads; borrow on the short leg costs money and can be recalled; and if you screen thousands of pairs and trade the best-looking one, its in-sample stats are inflated by multiple testing. The machinery here, OLS hedge, ADF-by-hand, rolling z, OU half-life, shift-clean sizing, is exactly what Avellaneda & Lee and Chan formalize. Compare it against the full Statistical Arbitrage Pairs database entry, then go break it on real tickers.