Pairs Trading: Cointegration Strategy in Python

What Is Pairs Trading?

Pairs trading is a market-neutral strategy. You simultaneously buy one stock and short another, betting that their price relationship will revert to the mean. Because you are long one and short the other, you are largely immune to market-wide moves. If the FTSE drops 5%, both legs lose value together, netting close to zero impact.

The classic example in London markets: Shell and BP. Both are major oil companies, driven by the same oil price, same regulation, same macro environment. When Shell temporarily outperforms BP by an unusual amount, you short Shell and buy BP, betting they will converge.

Testing for Cointegration

The Engle-Granger test checks if the spread between two time series is stationary (mean-reverting). A p-value below 0.05 suggests the pair is cointegrated. We use statsmodels for the statistical test and yfinance for market data.

import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import coint, adfuller
import yfinance as yf

# Download two correlated stocks (e.g. Shell and BP)
start = "2020-01-01"
end = "2023-12-31"

stock1 = yf.download("SHEL.L", start=start, end=end)["Close"]
stock2 = yf.download("BP.L", start=start, end=end)["Close"]

# Engle-Granger cointegration test
score, pvalue, _ = coint(stock1, stock2)
print(f"Cointegration test p-value: {pvalue:.4f}")
print(f"Cointegrated: {pvalue < 0.05}")

The Spread and Z-Score

The spread is Stock1 - β × Stock2, where β is the hedge ratio from linear regression. We normalise it into a z-score: how many standard deviations the spread is from its rolling mean. Z-scores above 2 or below -2 signal potential trades.

from sklearn.linear_model import LinearRegression

# Find hedge ratio via linear regression
model = LinearRegression()
model.fit(stock2.values.reshape(-1, 1), stock1.values)
hedge_ratio = model.coef_[0]
print(f"Hedge ratio: {hedge_ratio:.4f}")

# Calculate the spread
spread = stock1 - hedge_ratio * stock2

# Calculate z-score (rolling window)
window = 30
spread_mean = spread.rolling(window).mean()
spread_std = spread.rolling(window).std()
zscore = (spread - spread_mean) / spread_std

Full Backtest

Here is a complete pairs trading backtest. It tracks position, calculates daily P&L, and reports key metrics including Sharpe ratio and win rate.

def pairs_backtest(stock1, stock2, hedge_ratio, window=30,
                   entry_z=2.0, exit_z=0.5, stop_z=4.0):
    """
    Backtest a pairs trading strategy.
    Long spread when z < -entry_z, short when z > entry_z.
    Exit when |z| < exit_z. Stop loss at |z| > stop_z.
    """
    spread = stock1 - hedge_ratio * stock2
    spread_mean = spread.rolling(window).mean()
    spread_std = spread.rolling(window).std()
    zscore = (spread - spread_mean) / spread_std

    position = 0  # 1 = long spread, -1 = short spread
    positions = []
    pnl = []
    trades = 0

    for i in range(window, len(zscore)):
        z = zscore.iloc[i]
        spread_return = spread.iloc[i] - spread.iloc[i-1]

        # P&L from existing position
        pnl.append(position * spread_return)

        # Entry signals
        if position == 0:
            if z > entry_z:
                position = -1  # Short the spread
                trades += 1
            elif z < -entry_z:
                position = 1   # Long the spread
                trades += 1

        # Exit signals
        elif position == 1 and (z > -exit_z or z > stop_z):
            position = 0
        elif position == -1 and (z < exit_z or z < -stop_z):
            position = 0

        positions.append(position)

    pnl = pd.Series(pnl)
    total_return = pnl.sum()
    sharpe = pnl.mean() / pnl.std() * np.sqrt(252) if pnl.std() > 0 else 0

    print(f"Total P&L: {total_return:.2f}")
    print(f"Sharpe Ratio: {sharpe:.2f}")
    print(f"Total Trades: {trades}")
    print(f"Win Rate: {(pnl[pnl > 0].count() / pnl[pnl != 0].count() * 100):.1f}%")

    return pnl, positions

pnl, positions = pairs_backtest(stock1, stock2, hedge_ratio)

Scanning for Pairs

In practice you scan a universe of stocks and test all possible pairs for cointegration. With N stocks, there are N(N-1)/2 pairs to test. Beware of multiple testing: with 100 stocks you test 4,950 pairs, and some will appear cointegrated by chance. Apply Bonferroni correction or use a stricter threshold.

def find_cointegrated_pairs(symbols, start, end, pvalue_threshold=0.05):
    """Scan a universe of stocks for cointegrated pairs."""
    prices = {}
    for sym in symbols:
        try:
            prices[sym] = yf.download(sym, start=start, end=end)["Close"]
        except:
            continue

    pairs = []
    n = len(symbols)
    tested = list(prices.keys())

    for i in range(len(tested)):
        for j in range(i+1, len(tested)):
            s1, s2 = tested[i], tested[j]
            if len(prices[s1]) != len(prices[s2]):
                continue
            _, pvalue, _ = coint(prices[s1], prices[s2])
            if pvalue < pvalue_threshold:
                pairs.append((s1, s2, pvalue))
                print(f"  {s1} + {s2}: p={pvalue:.4f}")

    return sorted(pairs, key=lambda x: x[2])

# Scan FTSE 100 oil majors
oil_stocks = ["SHEL.L", "BP.L", "TotalEnergies.PA", "XOM", "CVX"]
pairs = find_cointegrated_pairs(oil_stocks, "2021-01-01", "2023-12-31")