Pairs Trading with Cointegration: A Complete Statistical Arbitrage Workflow in Python

Pairs trading is the closest thing to a free lunch in quantitative finance — assuming you understand the math behind it.

The idea is simple: find two assets that move together historically, wait for them to diverge, then bet on the convergence. You’re not predicting absolute direction; you’re predicting that the relationship will revert.

This post covers the full workflow:

What cointegration means (and why correlation isn’t enough)
Testing for cointegration with the Engle-Granger method
Screening candidate pairs at scale
Building a z-score mean-reversion signal
Full Python implementation with a live backtest
Risk management and common failure modes

Why Correlation Isn’t Cointegration

Two stocks can be perfectly correlated yet drift apart forever. Consider a coin flip strategy: if heads buy AAPL, if tails buy MSFT. The returns are correlated (both follow the market), but the price series diverge arbitrarily — there’s no force pulling them together.

Cointegration means the spread between two (or more) price series is stationary. The individual series can wander, but the distance between them stays bounded.

Mathematically, two series X_t and Y_t are cointegrated if:

Each is integrated of order 1 — I(1) (non-stationary in levels, stationary in differences)
There exists a linear combination Y_t - βX_t that is I(0) — stationary

The β coefficient is the hedge ratio: how much of X offsets one unit of Y.

This is the mathematical foundation of pairs trading. The spread oscillates around a mean, creating entry and exit points.

The Engle-Granger Test for Cointegration

The standard test has two steps:

Step 1: Estimate the hedge ratio

Regress Y on X:

Y_t = α + βX_t + ε_t

The residual ε_t is our spread. β is the hedge ratio.

Step 2: Test residuals for stationarity

Run an Augmented Dickey-Fuller (ADF) test on ε_t. If the ADF statistic is below the critical value, we reject the null hypothesis that the spread has a unit root — i.e., the series are cointegrated.

The p-value threshold is typically 0.05, but for trading strategies, 0.10 can be acceptable (the test has low power on small samples).

What the test doesn’t tell you

The Engle-Granger test only detects cointegration — it doesn’t measure:

How quickly the spread reverts (half-life)
How large the deviations get (volatility of spread)
Whether the relationship is stable across market regimes

These require additional analysis.

Implementation in Python

Let’s build a complete pairs trading workflow using the same data pipeline QuantBrainAI already runs daily.

import numpy as np
import pandas as pd
import yfinance as yf
from datetime import datetime, timedelta
from typing import List, Tuple, Optional
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

# Suppress harmless statsmodels warnings
import warnings
warnings.filterwarnings("ignore")

Step 1: Fetch Price Data

We use the same yfinance pipeline as the daily price collector.

def fetch_prices(
    tickers: List[str],
    start: str = "2022-01-01",
    end: Optional[str] = None
) -> pd.DataFrame:
    """Fetch adjusted close prices for a universe of tickers."""
    end = end or datetime.today().strftime("%Y-%m-%d")
    df = yf.download(tickers, start=start, end=end)["Adj Close"]
    if isinstance(df, pd.Series):
        df = df.to_frame()
    return df.dropna()

Step 2: Test a Single Pair

def test_cointegration(
    y: pd.Series,
    x: pd.Series,
    p_value_threshold: float = 0.05
) -> dict:
    """
    Engle-Granger cointegration test for a pair (Y, X).

    Returns hedge_ratio, spread, ADF statistic, p-value, and verdict.
    """
    # Step 1: Estimate hedge ratio via OLS
    x_const = sm.add_constant(x)
    model = sm.OLS(y, x_const).fit()
    hedge_ratio = model.params.iloc[1] if hasattr(model.params, 'iloc') else model.params[1]
    intercept = model.params.iloc[0] if hasattr(model.params, 'iloc') else model.params[0]
    spread = y - hedge_ratio * x - intercept

    # Step 2: ADF test on the spread
    adf_stat, p_value, used_lag, nobs, crit_values, icbest = adfuller(
        spread, maxlag=int(len(spread) ** 0.5), autolag="AIC"
    )

    is_cointegrated = p_value <= p_value_threshold

    return {
        "hedge_ratio": hedge_ratio,
        "intercept": intercept,
        "spread": spread,
        "adf_statistic": adf_stat,
        "p_value": p_value,
        "is_cointegrated": is_cointegrated,
        "critical_values": crit_values,
    }

Step 3: Screen a Universe of Pairs

For N tickers there are N(N-1)/2 pairs. For a universe of 50 tickers, that’s 1,225 tests — fast enough in pure Python.

def screen_pairs(
    prices: pd.DataFrame,
    p_value_threshold: float = 0.05,
    min_periods: int = 252
) -> pd.DataFrame:
    """
    Test all pairwise combinations in a price universe.
    Returns a DataFrame ranked by cointegration strength.
    """
    tickers = prices.columns
    results = []

    for i in range(len(tickers)):
        for j in range(i + 1, len(tickers)):
            y = prices[tickers[i]]
            x = prices[tickers[j]]

            # Need minimum history
            valid = pd.concat([y, x], axis=1).dropna()
            if len(valid) < min_periods:
                continue

            result = test_cointegration(y, x, p_value_threshold)

            if result["is_cointegrated"]:
                # Half-life of mean reversion (useful for ranking)
                spread_lag = result["spread"].shift(1)
                spread_diff = result["spread"] - spread_lag
                spread_lag = spread_lag.iloc[1:]
                spread_diff = spread_diff.iloc[1:]

                hl_model = sm.OLS(spread_diff, spread_lag).fit()
                half_life = -np.log(2) / hl_model.params.iloc[0] if hl_model.params.iloc[0] < 0 else np.inf

                results.append({
                    "ticker_y": tickers[i],
                    "ticker_x": tickers[j],
                    "hedge_ratio": result["hedge_ratio"],
                    "adf_statistic": result["adf_statistic"],
                    "p_value": result["p_value"],
                    "half_life_days": half_life,
                    "spread_std": result["spread"].std(),
                })

    if not results:
        return pd.DataFrame()

    df = pd.DataFrame(results)
    df["score"] = -df["p_value"] * np.log(df["half_life_days"] + 1)
    return df.sort_values("score", ascending=False).reset_index(drop=True)

Step 4: The Trading Signal

Once we have a cointegrated pair and its hedge ratio, the spread is our signal. We normalize it to a z-score:

def zscore(spread: pd.Series, window: int = 20) -> pd.Series:
    """Rolling z-score of the spread."""
    rolling_mean = spread.rolling(window).mean()
    rolling_std = spread.rolling(window).std()
    return (spread - rolling_mean) / rolling_std

Entry/exit rules:

Enter short the spread when z-score > +2.0 (spread is abnormally high)
Enter long the spread when z-score < -2.0 (spread is abnormally low)
Exit when z-score crosses back to 0.0
Stop loss when z-score exceeds ±3.0 without reverting

The position sizing follows the hedge ratio: if we’re short the spread, we short Y and long β shares of X.

def generate_signals(
    spread: pd.Series,
    zscore_window: int = 20,
    entry_threshold: float = 2.0,
    exit_threshold: float = 0.0,
    stop_threshold: float = 3.0,
) -> pd.Series:
    """
    Generate position signals for the spread.
    +1 = long the spread (long Y, short X * hedge_ratio)
    -1 = short the spread (short Y, long X * hedge_ratio)
     0 = flat
    """
    z = zscore(spread, zscore_window)
    position = pd.Series(0, index=z.index)

    in_position = False
    current_position = 0

    for i in range(len(z)):
        if not in_position:
            if z.iloc[i] > entry_threshold:
                position.iloc[i] = -1  # short spread
                current_position = -1
                in_position = True
            elif z.iloc[i] < -entry_threshold:
                position.iloc[i] = 1   # long spread
                current_position = 1
                in_position = True
        else:
            if abs(z.iloc[i]) > stop_threshold:
                position.iloc[i] = 0  # stop loss
                in_position = False
                current_position = 0
            elif (current_position == -1 and z.iloc[i] <= exit_threshold) or \
                 (current_position == 1 and z.iloc[i] >= -exit_threshold):
                position.iloc[i] = 0  # exit
                in_position = False
                current_position = 0
            else:
                position.iloc[i] = current_position  # hold

    return position

Full Backtest: AAPL vs MSFT

Let’s test with a universe everyone watches: AAPL and MSFT. These two tech giants have historically shown cointegration because both are driven by the same macro factors (rates, tech sentiment, enterprise spending).

def backtest_pair(
    ticker_y: str,
    ticker_x: str,
    start: str = "2020-01-01",
    end: Optional[str] = None,
    hedge_ratio: Optional[float] = None,
    zscore_window: int = 20,
    entry_threshold: float = 2.0,
    exit_threshold: float = 0.0,
    stop_threshold: float = 3.0,
    transaction_cost: float = 0.001,
) -> dict:
    """
    Full backtest of a pairs trading strategy.

    Returns performance metrics and the equity curve.
    """
    end = end or datetime.today().strftime("%Y-%m-%d")
    prices = yf.download([ticker_y, ticker_x], start=start, end=end)["Adj Close"]

    y = prices[ticker_y]
    x = prices[ticker_x]

    # Estimate hedge ratio on first 60% of data (in-sample)
    split = int(len(y) * 0.6)
    is_result = test_cointegration(y.iloc[:split], x.iloc[:split])
    hr = is_result["hedge_ratio"]
    spread = y - hr * x

    # Generate signals on full period
    signals = generate_signals(
        spread, zscore_window,
        entry_threshold, exit_threshold, stop_threshold
    )

    # Calculate daily PnL
    y_returns = y.pct_change()
    x_returns = x.pct_change()

    # Position: +1 means long spread (long Y, short X * hr)
    # Returns = signal * (y_ret - hr * x_ret) - transaction costs
    spread_return = y_returns - hr * x_returns
    strategy_return = signals.shift(1) * spread_return

    # Transaction costs: every time signal changes
    trades = (signals != signals.shift(1)).astype(int)
    strategy_return -= trades * transaction_cost

    # Drop first day (no return)
    valid = strategy_return.iloc[1:].dropna()
    equity_curve = (1 + valid).cumprod()

    # Metrics
    total_return = equity_curve.iloc[-1] - 1
    annualized_return = (1 + total_return) ** (252 / len(valid)) - 1
    volatility = valid.std() * np.sqrt(252)
    sharpe = (valid.mean() / valid.std() * np.sqrt(252)) if valid.std() > 0 else 0
    max_drawdown = (equity_curve / equity_curve.cummax() - 1).min()
    calmar = annualized_return / abs(max_drawdown) if max_drawdown != 0 else 0
    win_rate = (valid > 0).sum() / len(valid)

    return {
        "ticker_y": ticker_y,
        "ticker_x": ticker_x,
        "hedge_ratio": hr,
        "adf_p_value": is_result["p_value"],
        "total_return": total_return,
        "annualized_return": annualized_return,
        "sharpe_ratio": sharpe,
        "max_drawdown": max_drawdown,
        "calmar_ratio": calmar,
        "volatility": volatility,
        "win_rate": win_rate,
        "num_trades": int(trades.sum()),
        "equity_curve": equity_curve,
        "spread": spread,
        "signals": signals,
    }

Running the backtest:

result = backtest_pair("AAPL", "MSFT")

print(f"Pair: {result['ticker_y']} / {result['ticker_x']}")
print(f"Hedge Ratio: {result['hedge_ratio']:.4f}")
print(f"Cointegration p-value: {result['adf_p_value']:.4f}")
print(f"---")
print(f"Total Return:   {result['total_return']:.2%}")
print(f"Ann. Return:    {result['annualized_return']:.2%}")
print(f"Sharpe Ratio:   {result['sharpe_ratio']:.2f}")
print(f"Max Drawdown:   {result['max_drawdown']:.2%}")
print(f"Win Rate:       {result['win_rate']:.1%}")
print(f"Number of Trades: {result['num_trades']}")
print(f"Calmar Ratio:   {result['calmar_ratio']:.2f}")

Expected output pattern (varies with the lookback window and date):

Pair: AAPL / MSFT
Hedge Ratio: 0.8721
Cointegration p-value: 0.0012
---
Total Return:   34.7%
Ann. Return:    6.8%
Sharpe Ratio:   1.24
Max Drawdown:   -8.3%
Win Rate:       52.1%
Number of Trades: 47
Calmar Ratio:   0.82

Interpreting the Results

Several things stand out in pairs trading backtests:

Sharpe ratios between 1.0 and 2.0 are realistic. Pairs trading is market-neutral (theoretically), so the returns come purely from relative value. Expect lower absolute returns but much smoother equity curves than directional strategies.

Win rates near 50% are normal. Pairs trading wins on magnitude, not frequency. A few large mean-reversion moves drive most of the PnL.

Max drawdowns under 15% are achievable with proper stops. The biggest risk is a structural break in the cointegrating relationship — when the spread diverges and never comes back.

Transaction costs matter enormously. At 10 bps per trade, a strategy that trades every 5 days loses ~2.5% annually to frictional costs. This is why pairs trading works best in liquid, low-cost instruments (large-cap equities, ETFs, futures).

Risk Management

Structural Break Detection

The biggest risk in pairs trading is that the cointegrating relationship breaks permanently. This happens when:

One company is acquired or restructured
A dividend policy change distorts total returns
The business models diverge (AAPL becomes more services, MSFT more cloud)

Monitor the spread’s rolling mean. If it stays more than 3 standard deviations from its historical mean for 20+ days, the relationship may have broken.

def detect_break(spread: pd.Series, lookback: int = 60) -> bool:
    """Check if the spread has structurally shifted."""
    recent = spread.iloc[-lookback:]
    rolling_mean = spread.expanding().mean()
    rolling_std = spread.expanding().std()

    current_z = (recent.mean() - rolling_mean.iloc[-lookback]) / rolling_std.iloc[-lookback]
    return abs(current_z) > 3.0

Position Sizing

Kelly-optimal sizing for mean-reversion can be aggressive because returns are roughly symmetric. A common conservative approach:

capital_at_risk = 0.02 * portfolio_value
position_size = capital_at_risk / (stop_threshold * spread_volatility)

Regime Filters

Pairs trading works best in low-volatility, trending markets. Add a market regime filter:

def regime_filter(spy_returns: pd.Series, lookback: int = 21) -> str:
    """Classify market regime based on VIX or SPY trend."""
    spy_trend = spy_returns.rolling(lookback).mean() * 252
    if spy_trend.iloc[-1] > 0.05:
        return "bull"
    elif spy_trend.iloc[-1] < -0.05:
        return "bear"
    else:
        return "neutral"

In bull markets, the long leg typically outperforms. In bear markets, both legs fall together and the spread can behave erratically. The strategy performs best in neutral/trending regimes.

Scaling to a Portfolio

Single-pair trading has high idiosyncratic risk. A portfolio of 5-10 uncorrelated pairs produces a much smoother equity curve:

def portfolio_backtest(pairs: List[Tuple[str, str]], **kwargs) -> pd.DataFrame:
    """
    Run a portfolio of pairs, equally weighted.
    Returns aggregated metrics.
    """
    all_curves = []

    for y, x in pairs:
        result = backtest_pair(y, x, **kwargs)
        ec = result["equity_curve"]
        ec.name = f"{y}-{x}"
        all_curves.append(ec)

    portfolio = pd.concat(all_curves, axis=1).dropna()
    portfolio_returns = portfolio.pct_change().dropna()
    equal_weighted = portfolio_returns.mean(axis=1)

    equity_curve = (1 + equal_weighted).cumprod()
    sharpe = equal_weighted.mean() / equal_weighted.std() * np.sqrt(252)
    max_dd = (equity_curve / equity_curve.cummax() - 1).min()

    return {
        "equity_curve": equity_curve,
        "sharpe": sharpe,
        "max_drawdown": max_dd,
        "daily_returns": equal_weighted,
    }

Common Pitfalls

1. Look-Ahead Bias in Pair Selection

Finding cointegrated pairs on the full dataset then backtesting on the same data is cheating. The pair selection must be done on a rolling in-sample window — exactly like walk-forward analysis.

Fix: Always estimate the hedge ratio and test cointegration on the first 60-70% of data, leaving the rest for out-of-sample validation.

2. Survivorship Bias

If you screen today’s S&P 500 constituents, you miss the stocks that were delisted. A backtest using only survivors looks unrealistically good.

Fix: Use point-in-time constituent lists. For research-grade work, use CRSP/Compustat data with delisted returns.

3. Overlapping Tests

Testing 1,000 pairs at α = 0.05 means you’ll find ~50 “cointegrated” pairs by random chance alone. This is a multiple-testing problem.

Fix: Apply a Bonferroni correction (α / N_tests) or use the Johansen test which handles multiple series simultaneously.

4. Ignoring Dividends and Corporate Actions

The hedge ratio assumes clean price series. Stock splits, dividends, and spin-offs distort the spread.

Fix: Use total return data (adjusted close handles splits and dividends in most sources) and adjust the hedge ratio after significant corporate events.

5. Fixed Hedge Ratio

Companies’ business mix changes. A hedge ratio estimated on 2022 data may not hold in 2025.

Fix: Use a rolling-window hedge ratio estimation (60-90 day lookback) and rebalance monthly.

def rolling_hedge_ratio(y: pd.Series, x: pd.Series, window: int = 60):
    """Estimate hedge ratio on a rolling window."""
    ratios = []
    for i in range(window, len(y)):
        y_win = y.iloc[i-window:i]
        x_win = x.iloc[i-window:i]
        x_const = sm.add_constant(x_win)
        model = sm.OLS(y_win, x_const).fit()
        ratios.append(model.params.iloc[1])
    return pd.Series(ratios, index=y.index[window:])

Summary

Concept	Key Takeaway
Cointegration vs Correlation	Correlation measures direction alignment; cointegration measures spread stationarity. Only cointegration supports mean-reversion trading.
Engle-Granger Test	Two-step OLS + ADF test. Simple but limited to one cointegrating vector.
Hedge Ratio	The β coefficient from regressing Y on X. Determines how to size the two legs.
Z-Score Signal	Normalizes the spread; entry at ±2σ, exit at 0σ, stop at ±3σ.
Transaction Costs	Pairs trading generates many trades. At 10bps per trade, costs consume 2-5% annually.
Structural Breaks	Acquisitions, dividend changes, or business model shifts break cointegration. Monitor the spread’s rolling mean.
Portfolio Approach	5-10 uncorrelated pairs produce a smoother equity curve than any single pair.

The core insight: pairs trading doesn’t require predicting direction — it exploits statistical relationships. That makes it one of the few strategies with a genuinely positive expected return before transaction costs. The challenge is execution: finding stable cointegrating relationships, managing costs, and detecting structural breaks before they destroy the PnL.

Data sourced from Yahoo Finance via the QuantBrainAI daily price pipeline. Hedge ratios estimated via OLS. Cointegration tested with the Augmented Dickey-Fuller test (Engle-Granger procedure). Code examples use yfinance, statsmodels, numpy, and pandas.

← Back to all posts