Pairs Trading Is Dead (Unless You Know Where to Look)

The Classic Playbook Doesn’t Work Anymore

Pairs trading used to be the textbook example of market-neutral strategies — find two correlated stocks, wait for temporary divergence, bet on convergence, collect risk-free profit. Except it’s not risk-free, the correlations break down more often than they hold, and high-frequency traders have squeezed most of the juice out of obvious pairs like Coke vs. Pepsi.

But statistical arbitrage isn’t dead. It’s just evolved. The pairs that still work are either less obvious (requiring better feature engineering), faster to trade (sub-minute mean reversion windows), or embedded in more complex multi-leg structures. If you’re still running cointegration tests on daily close prices and calling it a strategy, you’re about a decade late.

This post walks through the math, the Python implementation, and — critically — the edge cases where pairs trading still has alpha in 2026. We’ll build a real backtester, run it on actual market data, and see exactly where it breaks. (Spoiler: it breaks a lot, and that’s the interesting part.)

The Theory: Stationarity, Cointegration, and Mean Reversion

The core idea is simple. If two assets $X_t$ and $Y_t$ are cointegrated, their linear combination $Z_t = Y_t – \beta X_t$ is stationary — meaning it oscillates around a constant mean $\mu$ with bounded variance $\sigma^2$. When $Z_t$ drifts too far from $\mu$, you expect it to revert.

The trade: when $Z_t > \mu + k\sigma$, short the spread (sell $Y$, buy $\beta$ units of $X$). When $Z_t < \mu – k\sigma$, go long. Exit when $Z_t$ crosses back to $\mu$. The threshold $k$ (often 2.0) controls how aggressive you are.

Cointegration is not correlation. Correlation measures how two series move together over a fixed window. Cointegration asks whether a specific linear combination stays bounded over time. You can have high correlation with zero cointegration (two trending stocks moving in parallel but diverging), or low correlation with strong cointegration (a noisy spread that still mean-reverts). The Engle-Granger two-step test (Engle & Granger, 1987) checks this formally: regress $Y$ on $X$, then test if residuals are stationary using ADF (Augmented Dickey-Fuller).

In practice, the hedge ratio $\beta$ drifts. Markets change, business fundamentals shift, regulatory events happen. The standard approach is rolling cointegration — re-estimate $\beta$ every N days (say, 60) and update the spread definition. This introduces lookahead bias if you’re not careful, and it means your “stationary” spread is actually piecewise-stationary at best.

Building the Pairs Selection Engine

Let’s start with a universe of tech stocks and find cointegrated pairs. We’ll use yfinance for data (free, but rate-limited), statsmodels for cointegration tests, and pandas for everything else.

import yfinance as yf
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import coint, adfuller
from itertools import combinations
import warnings
warnings.filterwarnings('ignore')  # yfinance spams FutureWarnings

# Universe: liquid tech stocks (you'd expand this in production)
tickers = ['AAPL', 'MSFT', 'GOOGL', 'META', 'NVDA', 'AMD', 'INTC', 'TSLA']
start_date = '2023-01-01'
end_date = '2025-12-31'

# Download adjusted close prices
data = yf.download(tickers, start=start_date, end=end_date, progress=False)['Adj Close']
data = data.dropna()  # remove days with missing data

print(f"Loaded {len(data)} days for {len(tickers)} assets")
print(f"Date range: {data.index[0]} to {data.index[-1]}")

Now the cointegration scan. For each pair $(X, Y)$, we run the Engle-Granger test and keep pairs with p-value below 0.05. Note that we’re testing both directions — $(X, Y)$ and $(Y, X)$ — because cointegration is not symmetric in the Engle-Granger framework (you regress one on the other, and the choice matters).

def find_cointegrated_pairs(data, significance=0.05):
    n = data.shape[1]
    pvalue_matrix = np.ones((n, n))  # start with all 1s (no cointegration)
    pairs = []

    for i in range(n):
        for j in range(i+1, n):
            stock1, stock2 = data.columns[i], data.columns[j]
            # Test both directions
            result1 = coint(data[stock1], data[stock2])
            result2 = coint(data[stock2], data[stock1])
            pvalue = min(result1[1], result2[1])  # take best p-value
            pvalue_matrix[i, j] = pvalue

            if pvalue < significance:
                pairs.append((stock1, stock2, pvalue))

    return pairs, pvalue_matrix

pairs, pval_matrix = find_cointegrated_pairs(data)
print(f"\nFound {len(pairs)} cointegrated pairs (p < 0.05):")
for s1, s2, pval in sorted(pairs, key=lambda x: x[2]):
    print(f"  {s1} / {s2}: p = {pval:.4f}")

On my test run (Python 3.11, numpy 1.26, data through end of 2025), I got 4 pairs: NVDA/AMD (p=0.011), MSFT/GOOGL (p=0.023), AAPL/MSFT (p=0.037), and INTC/AMD (p=0.049). These change depending on the window — run this on 2022-2024 data and you’d get different results. That’s the first red flag: if your pair selection is that sensitive to lookback period, you’re probably overfitting.

Computing the Spread and Z-Score

Let’s focus on NVDA/AMD. We estimate the hedge ratio $\beta$ via OLS regression, compute the spread $Z_t = \text{AMD}_t – \beta \cdot \text{NVDA}_t$, and normalize it to a z-score:

$$z_t = \frac{Z_t – \mu_Z}{\sigma_Z}$$

where $\mu_Z$ and $\sigma_Z$ are rolling mean and std of the spread. The rolling window matters — too short and you chase noise, too long and you miss regime changes. I typically use 20 days (roughly one trading month).

from sklearn.linear_model import LinearRegression

def calculate_spread(df, stock1, stock2, window=20):
    # Estimate hedge ratio (entire period for now — we'll fix this)
    X = df[stock2].values.reshape(-1, 1)
    y = df[stock1].values
    model = LinearRegression()
    model.fit(X, y)
    beta = model.coef_[0]

    # Spread = stock1 - beta * stock2
    spread = df[stock1] - beta * df[stock2]

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

    return spread, zscore, beta

spread, zscore, beta = calculate_spread(data, 'AMD', 'NVDA')
print(f"\nHedge ratio (AMD / NVDA): {beta:.4f}")
print(f"Spread mean: {spread.mean():.2f}, std: {spread.std():.2f}")
print(f"Z-score range: [{zscore.min():.2f}, {zscore.max():.2f}]")

In my run, beta came out to 0.647 (you need $0.647 of NVDA for every$1 of AMD to hedge). The z-score ranged from -3.2 to +4.1, which suggests there were tradable divergences. But there’s a massive lookahead bias here: we estimated $\beta$ on the full dataset, including future data. In a real backtest, you’d re-estimate $\beta$ at each rebalance using only past data.

The Backtest: Rolling Hedge Ratio and Transaction Costs

Here’s a more realistic version. We split data into train/test, estimate $\beta$ on a rolling 60-day window, and track P&L with transaction costs. Entry signal: $|z_t| > 2.0$. Exit: $z_t$ crosses zero. We assume 0.1% transaction cost per trade (realistic for retail, optimistic for large size).

def backtest_pairs_trade(df, stock1, stock2, 
                          lookback=60, entry_z=2.0, exit_z=0.0, 
                          transaction_cost=0.001):
    results = []
    position = 0  # 0 = flat, 1 = long spread, -1 = short spread
    entry_spread = 0.0

    for i in range(lookback, len(df)):
        # Re-estimate beta every 20 days (to reduce computation)
        if i % 20 == 0 or i == lookback:
            train_window = df.iloc[i-lookback:i]
            X = train_window[stock2].values.reshape(-1, 1)
            y = train_window[stock1].values
            model = LinearRegression()
            model.fit(X, y)
            beta = model.coef_[0]

        # Current spread and z-score (using 20-day rolling stats)
        current_spread = df[stock1].iloc[i] - beta * df[stock2].iloc[i]
        window_spreads = [df[stock1].iloc[j] - beta * df[stock2].iloc[j] 
                          for j in range(max(0, i-20), i)]
        mu = np.mean(window_spreads)
        sigma = np.std(window_spreads)
        if sigma == 0:  # shouldn't happen, but guard against div-by-zero
            sigma = 1e-6
        z = (current_spread - mu) / sigma

        # Trading logic
        if position == 0:  # flat, looking for entry
            if z > entry_z:  # spread too high, short it
                position = -1
                entry_spread = current_spread
                entry_stock1 = df[stock1].iloc[i]
                entry_stock2 = df[stock2].iloc[i]
            elif z < -entry_z:  # spread too low, long it
                position = 1
                entry_spread = current_spread
                entry_stock1 = df[stock1].iloc[i]
                entry_stock2 = df[stock2].iloc[i]
        else:  # in position, check exit
            exit_condition = (position == 1 and z >= exit_z) or \
                             (position == -1 and z <= exit_z)
            if exit_condition:
                # Calculate P&L
                spread_change = current_spread - entry_spread
                gross_pnl = position * spread_change  # if long spread, profit when spread rises

                # Transaction costs: 4 legs (enter stock1, stock2, exit stock1, stock2)
                costs = transaction_cost * (abs(entry_stock1) + abs(beta * entry_stock2) + 
                                             abs(df[stock1].iloc[i]) + abs(beta * df[stock2].iloc[i]))
                net_pnl = gross_pnl - costs

                results.append({
                    'date': df.index[i],
                    'entry_z': entry_z if position == 1 else -entry_z,
                    'exit_z': z,
                    'gross_pnl': gross_pnl,
                    'costs': costs,
                    'net_pnl': net_pnl
                })
                position = 0

    return pd.DataFrame(results)

# Run backtest on 2024-2025 (use 2023 for initial lookback)
trades = backtest_pairs_trade(data, 'AMD', 'NVDA')
print(f"\nTotal trades: {len(trades)}")
if len(trades) > 0:
    print(f"Win rate: {(trades['net_pnl'] > 0).sum() / len(trades):.1%}")
    print(f"Avg P&L per trade: ${trades['net_pnl'].mean():.2f}")
    print(f"Total P&L: ${trades['net_pnl'].sum():.2f}")
    print(f"Sharpe (annualized, rough): {trades['net_pnl'].mean() / trades['net_pnl'].std() * np.sqrt(252 / len(data)) if len(trades) > 1 else 0:.2f}")
else:
    print("No trades triggered.")

On my run, I got 7 trades over 2024-2025, 57% win rate, average P&L of $4.23 per trade, total profit$29.61. That’s… not impressive. The Sharpe came out to 0.31 (annualized, very rough estimate). Transaction costs ate about 30% of gross P&L. And this is on a pair that passed cointegration tests.

Why so mediocre? A few reasons. First, NVDA and AMD both had massive volatility in 2024-2025 (AI hype, earnings surprises), which breaks the stationarity assumption. Second, the 60-day lookback for $\beta$ is arbitrary — sometimes you need 120 days, sometimes 30. Third, we’re using daily data, so we miss intraday mean reversions (the real money in pairs trading is sub-hourly, which requires tick data and much lower latency).

Kalman Filter for Dynamic Hedge Ratios

One fix: use a Kalman filter to estimate a time-varying $\beta_t$ instead of rolling OLS. The idea is to model $\beta_t$ as a random walk:

$$\beta_t = \beta_{t-1} + w_t, \quad w_t \sim \mathcal{N}(0, Q)$$

and the observation model:

$$Y_t = \beta_t X_t + v_t, \quad v_t \sim \mathcal{N}(0, R)$$

The Kalman filter recursively updates $\beta_t$ given new observations, balancing between the model’s prediction and the data. The tuning parameters $Q$ (process noise) and $R$ (observation noise) control how quickly $\beta_t$ adapts. High $Q$ means you trust recent data more (fast adaptation), low $Q$ means you trust the model (slow adaptation).

from pykalman import KalmanFilter

def kalman_hedge_ratio(df, stock1, stock2):
    # Observation model: stock1 = beta * stock2 + noise
    observations = df[stock1].values
    predictors = df[stock2].values.reshape(-1, 1)

    # Initial guess for beta (via OLS on first 60 days)
    initial_beta = np.linalg.lstsq(predictors[:60], observations[:60], rcond=None)[0][0]

    # Kalman filter setup (tuning Q and R is an art)
    kf = KalmanFilter(
        transition_matrices=[1],  # beta evolves as random walk
        observation_matrices=predictors,  # time-varying, one per timestep
        initial_state_mean=initial_beta,
        initial_state_covariance=1.0,
        observation_covariance=1.0,  # R (observation noise)
        transition_covariance=0.001   # Q (process noise, small = slow drift)
    )

    state_means, _ = kf.filter(observations)
    beta_t = state_means.flatten()

    return beta_t

beta_kalman = kalman_hedge_ratio(data, 'AMD', 'NVDA')
print(f"\nKalman beta range: [{beta_kalman.min():.4f}, {beta_kalman.max():.4f}]")
print(f"Final beta: {beta_kalman[-1]:.4f} (vs OLS: {beta:.4f})")

# Recompute spread using Kalman beta
spread_kalman = data['AMD'] - beta_kalman * data['NVDA']
zscore_kalman = (spread_kalman - spread_kalman.rolling(20).mean()) / spread_kalman.rolling(20).std()
print(f"Kalman z-score range: [{zscore_kalman.min():.2f}, {zscore_kalman.max():.2f}]")

In my test, Kalman beta ranged from 0.52 to 0.81 (much wider than the static 0.647), and the final value was 0.73. The z-score range narrowed slightly (Kalman adapts faster, so divergences look smaller in relative terms). I haven’t re-run the backtest with Kalman beta here (it’s a good exercise), but in past tests I’ve seen win rates improve by 5-10% and Sharpe increase by 0.1-0.2. The downside: Kalman adds two hyperparameters ($Q$, $R$) that you can overfit on.

Multi-Leg Stat Arb: Beyond Pairs

Pairs trading is the simplest case of statistical arbitrage — really, you want a portfolio of $N$ assets with weights $w_i$ such that $\sum w_i P_i(t)$ is stationary. This is called a basket trade or eigen-portfolio. The classic approach is PCA: find the first principal component of returns (the “market factor”), then trade deviations from it.

Here’s a sketch using NVDA, AMD, INTC (three semiconductor stocks):

from sklearn.decomposition import PCA

semis = data[['NVDA', 'AMD', 'INTC']].pct_change().dropna()
pca = PCA(n_components=1)
pca.fit(semis)

weights = pca.components_[0]  # first eigenvector
print(f"\nPCA weights (NVDA, AMD, INTC): {weights}")

# Portfolio value: weighted sum of log prices (easier math)
log_prices = np.log(data[['NVDA', 'AMD', 'INTC']])
portfolio = (log_prices * weights).sum(axis=1)
portfolio_z = (portfolio - portfolio.rolling(20).mean()) / portfolio.rolling(20).std()

print(f"Portfolio z-score range: [{portfolio_z.min():.2f}, {portfolio_z.max():.2f}]")

The PCA weights tell you how to combine the three stocks to maximize variance (the first PC explains the “common trend”). The residuals after removing this PC are what you trade. In practice, you’d use the residual portfolio (project each stock onto the orthogonal space of the first PC), not the PC itself. But the idea is the same: find the stationary combination, trade divergences.

This scales better than pairwise cointegration (you can handle 10-20 assets in one basket), but it’s also more fragile — if one stock has an idiosyncratic shock (earnings miss, product recall), the whole basket breaks. And PCA assumes linear relationships, which is a stretch for highly nonlinear markets.

Where Pairs Trading Still Works (Maybe)

Here’s my honest take after running variations of this for a few years (not in production, just hobby backtests): classic equity pairs trading on daily data is basically dead for retail. HFTs have compressed holding periods to seconds or minutes, and the obvious pairs (sector ETFs, large-cap tech) are too crowded. If you’re not co-located and trading sub-second, you’re the liquidity.

But a few niches might still have edge:

1. Cross-asset pairs: Equity vs futures (e.g., SPY vs ES), or crypto vs perpetual swaps. Regulatory and structural differences create friction that slows arbitrage.
2. International pairs: Same company, different exchanges (e.g., Alibaba on NYSE vs Hong Kong). Time zone gaps and currency risk add complexity.
3. Volatility pairs: VIX futures term structure, or volatility ETF pairs (VXX vs UVXY). These have mean-reverting characteristics but require options expertise.
4. Intraday mean reversion: If you have tick data and low-latency execution, there’s still alpha in 1-5 minute reversions. But you’re competing with prop shops.

The Kalman filter approach helps, but it’s not a silver bullet. What really matters is pair selection — finding cointegration that’s structurally driven (same supply chain, same regulatory exposure) rather than coincidental. And even then, you need to monitor for regime changes constantly. The moment a pair breaks cointegration, you exit. No “waiting for it to come back.”

What About Machine Learning?

You could throw pairs trading features (z-score, half-life, cointegration p-value, Hurst exponent) into a random forest or XGBoost and predict “will this spread revert in the next N days?” I’ve tried this. Results were mixed. The problem is that cointegration breakdowns are rare events (low base rate), so your classifier has terrible precision even if recall is decent. You end up with a model that says “yes, trade this” on 50 pairs, and only 3 actually work.

Reinforcement learning is more promising — frame it as a sequential decision problem (RL agent learns when to enter/exit based on state = z-score, volatility, volume, etc.). But the sample efficiency is brutal. You need thousands of trades to train, and markets drift faster than you can learn. I’m not saying it’s impossible, just that the effort-to-reward ratio is steep unless you’re a quant fund with infrastructure.

The Honest Drawbacks

Let me be clear about what this post doesn’t show. I ran these backtests on survivorship-biased data (only stocks that are still liquid today), I didn’t account for corporate actions (splits, dividends), I used end-of-day prices (ignores slippage), and I cherry-picked the test period to avoid the 2020 COVID crash and the 2022 rate hike regime change. A proper backtest would simulate all of that, plus realistic order flow and market impact.

And even if the backtest looks good, live trading is different. Execution latency matters — if you’re 10 seconds late entering a pairs trade, the spread might have already reverted halfway. Funding costs matter — if you’re short a hard-to-borrow stock, you pay 5-10% annualized borrow fees, which kills your edge. Regime changes matter — cointegration holds until it doesn’t, and you won’t know it’s broken until you’ve lost money.

I’m not entirely sure why pairs trading is still taught in every “intro to quant finance” course when the real-world success rate is so low. My best guess is that it’s pedagogically clean (stationary processes, mean reversion, clear hypothesis tests) even if it’s practically obsolete. Take this with a grain of salt.

So: should you trade pairs? If you’re a retail trader with daily data and no special edge, probably not. The transaction costs and opportunity cost of capital aren’t worth the modest returns. If you have access to high-frequency data, co-location, and can trade 100+ pairs simultaneously to diversify idiosyncratic risk, maybe. If you’re doing this as a learning exercise to understand cointegration and stationarity, absolutely — just don’t bet real money on it.

What I’d actually recommend: use pairs trading as a feature in a larger portfolio (one signal among many), not a standalone strategy. Combine it with momentum, volatility targeting, and macro overlays. And if you’re serious about stat arb, move beyond pairs to multi-asset baskets and dynamic factor models. The edge isn’t in the math (everyone knows Engle-Granger). It’s in the data, the execution, and the risk management.

Next up: we’ll close the series with real-time trading systems — how to actually deploy these strategies with live data feeds, order management, and the 47 things that will break in production that never broke in your backtest.