"Two stocks that should move together but occasionally drift apart — that gap is not noise. It is signal."
In February 2024, ExxonMobil (XOM) and Chevron (CVX) — two energy giants whose returns had maintained a 0.91 correlation over the preceding five years — diverged by more than 4.5 standard deviations during a single trading session. The spread, which typically oscillated within a ±2σ band, gapped open and remained elevated for 72 hours before mean-reverting. Traders running a cointegration-based pairs strategy captured a 3.2% spread return in that window, while those relying on simple correlation models missed the entry entirely.
The difference is not luck. It is methodology.
Correlation measures whether two assets move in the same direction. Cointegration measures whether two assets share a common long-run equilibrium — even when they diverge in the short run. This distinction is critical for statistical arbitrage, because a correlation of 0.85 can coexist with a cointegration test failure. Two stocks can trend together chaotically while their linear combination still wanders off to infinity. Only cointegration identifies the stable long-run relationship that makes mean-reversion trading defensible.
This article builds a complete pairs trading system from scratch: from screening thousands of US equities for cointegration pairs, through the Engle-Granger two-step testing framework, to real-time spread monitoring with Z-Score signals and production-grade WebSocket deployment.
The Pain Point: Why Pairs Screening Fails at Scale
Before building the system, it is worth diagnosing why most retail implementations of pairs trading underperform.
Running a full cointegration test on all possible pairs from a universe of 3,000 stocks produces 4.5 million pair combinations. Each test involves OLS regression, residual computation, and an Augmented Dickey-Fuller (ADF) test — operations that, without optimization, can take hours on a standard laptop. The naive approach has three fatal flaws:
First, computational cost scales quadratically. A 3,000-stock universe at 252 trading days generates 756,000 data points per stock. A brute-force O(n²) screening pass exceeds reasonable computation budgets for individual traders.
Second, stationarity is not binary. The ADF test returns a p-value, not a yes/no answer. Traders who filter on p < 0.05 miss pairs with p = 0.07 that may still offer trading edges under specific regime conditions.
Third, cointegration is time-varying. A pair that passes the Engle-Granger test over 2019–2021 may fail during the 2022 bear market, when sector rotations introduce structural breaks. Static pair selection produces strategies that decay.
The system below addresses all three problems with a tiered screening pipeline, rolling window validation, and real-time drift detection.
Tier 1: Universe Construction and Pre-Filtering
The first stage narrows 3,000+ US equities to a computationally manageable candidate set using heuristic filters before running formal cointegration tests.
Pre-Filter Criteria
| Criterion | Threshold | Rationale |
|---|---|---|
| Average daily volume | > $5M | Ensures bid-ask spreads narrow enough for spread arbitrage |
| Market cap | > $1B | Reduces micro-cap noise and liquidity risk |
| Sector classification | Must match | Cross-sector pairs can be cointegrated but face regime risk from structural industry shocks |
| Price history | ≥ 252 trading days | Minimum sample for ADF test validity |
| Missing data rate | < 5% | Gaps introduce look-ahead bias if not handled carefully |
The sector constraint is deliberate. Pairs from different sectors (e.g., energy vs. technology) may be cointegrated during specific macro regimes but break down during sector rotations. A XOM-CVX energy pair and a JPM-BAC financial pair both pass sector screening — but they carry different regime risk profiles.
import pandas as pd
import numpy as np
import os
from itertools import combinations
def load_universe_data(api_key: str, symbols: list, start_date: str, end_date: str) -> pd.DataFrame:
"""
Load adjusted close prices for a list of symbols.
Requires daily kline data over the specified window.
Args:
api_key: TickDB API key
symbols: List of ticker symbols (e.g., ['AAPL.US', 'MSFT.US'])
start_date: Start date in YYYY-MM-DD format
end_date: End date in YYYY-MM-DD format
Returns:
DataFrame with datetime index and adjusted close columns per symbol
"""
import requests
base_url = "https://api.tickdb.ai/v1/market/kline"
headers = {"X-API-Key": api_key}
price_data = {}
for symbol in symbols:
params = {
"symbol": symbol,
"interval": "1d",
"start": pd.Timestamp(start_date).value // 10**6,
"end": pd.Timestamp(end_date).value // 10**6,
"limit": 500
}
try:
response = requests.get(
base_url,
headers=headers,
params=params,
timeout=(3.05, 10)
)
if response.status_code != 200:
print(f"Warning: {symbol} returned status {response.status_code}")
continue
data = response.json()
if data.get("code") != 0:
print(f"Warning: {symbol} API error {data.get('code')}: {data.get('message')}")
continue
records = data.get("data", {}).get("klines", [])
if not records:
continue
closes = {r["t"]: r["c"] for r in records}
price_data[symbol] = closes
except requests.exceptions.Timeout:
print(f"Timeout loading {symbol} — skipping")
except Exception as e:
print(f"Error loading {symbol}: {e}")
# Align all symbols to common timestamps
if not price_data:
raise ValueError("No data loaded — check API key and symbol list")
df = pd.DataFrame(price_data)
df.index = pd.to_datetime(df.index, unit="ms")
df.index.name = "date"
return df.sort_index()
def prefilter_universe(df: pd.DataFrame, min_avg_volume: float = 5e6,
max_missing_pct: float = 0.05) -> list:
"""
Apply heuristic filters to reduce pair-search space.
Args:
df: Price DataFrame from load_universe_data
min_avg_volume: Minimum average daily dollar volume (default $5M)
max_missing_pct: Maximum allowed fraction of missing data per symbol
Returns:
Filtered list of symbols that pass all criteria
"""
candidates = []
for col in df.columns:
series = df[col].dropna()
# Check missing data threshold
missing_pct = df[col].isna().mean()
if missing_pct > max_missing_pct:
continue
# Check minimum data length (252 trading days)
if len(series) < 252:
continue
# Check for sufficient price variation (constant prices are not tradeable)
if series.std() < series.mean() * 0.01:
continue
candidates.append(col)
print(f"Pre-filter pass: {len(candidates)}/{len(df.columns)} symbols")
return candidates
⚠️ Engineering note: The load_universe_data function fetches daily kline data sequentially. For a 3,000-symbol universe, this requires approximately 3,000 API calls. Implement batch fetching or async concurrent requests for production use. The current implementation is designed for clarity, not throughput. Rate limiting compliance is the caller's responsibility.
Tier 2: Cointegration Testing — The Engle-Granger Two-Step Framework
With candidates narrowed to a sector-matched subset, formal cointegration testing begins. The standard approach is the Engle-Granger two-step procedure:
Step 1: Regress Y on X and compute residuals (the spread).
Step 2: Apply the Augmented Dickey-Fuller test to the residuals. If the ADF statistic is significant (p < threshold), the pair is cointegrated.
Why Not Johansen?
The Johansen test handles multiple cointegrating relationships simultaneously and is superior for multivariate systems. For a two-stock pairs trading system, Engle-Granger is sufficient and more interpretable — the hedge ratio has a direct meaning (how many shares of X to short per share of Y).
The Cointegration Tester
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.regression.linear_model import OLS
import warnings
warnings.filterwarnings("ignore")
class CointegrationPair:
"""Represents a cointegrated pair with test statistics and hedge ratio."""
def __init__(self, symbol_a: str, symbol_b: str, hedge_ratio: float,
adf_stat: float, p_value: float, half_life: float):
self.symbol_a = symbol_a
self.symbol_b = symbol_b
self.hedge_ratio = hedge_ratio
self.adf_stat = adf_stat
self.p_value = p_value
self.half_life = half_life # Mean reversion half-life in days
def __repr__(self):
return (f"Pair({self.symbol_a}-{self.symbol_b}, "
f"ADF={self.adf_stat:.3f}, p={self.p_value:.4f}, "
f"HL={self.half_life:.1f}d)")
def compute_spread_and_hedge_ratio(series_a: pd.Series, series_b: pd.Series) -> tuple:
"""
Step 1 of Engle-Granger: OLS regression to compute spread.
Args:
series_a: Price series for asset A
series_b: Price series for asset B
Returns:
(hedge_ratio, spread_series)
"""
# Align series by index
aligned = pd.DataFrame({"a": series_a, "b": series_b}).dropna()
if len(aligned) < 100:
raise ValueError("Insufficient aligned observations for regression")
# OLS: A = β * B + ε
# The hedge ratio tells us how many units of B to short per unit of A
model = OLS(aligned["a"], aligned["b"]).fit()
hedge_ratio = model.params.iloc[0]
# Spread = A - β * B (this is the trading signal)
spread = aligned["a"] - hedge_ratio * aligned["b"]
return hedge_ratio, spread
def adf_test_spread(spread: pd.Series, max_lag: int = 12) -> dict:
"""
Step 2 of Engle-Granger: ADF test on spread residuals.
Args:
spread: Residual series from Step 1
max_lag: Maximum lag order for ADF test
Returns:
Dictionary with ADF statistic, p-value, and critical values
"""
result = adfuller(spread, maxlag=max_lag, autolag="AIC")
return {
"adf_statistic": result[0],
"p_value": result[1],
"used_lag": result[2],
"n_obs": result[3],
"critical_values": result[4]
}
def estimate_half_life(spread: pd.Series) -> float:
"""
Estimate mean-reversion half-life using Ornstein-Uhlenbeck decay.
The half-life tells us how many days the spread typically takes
to revert halfway to its mean — critical for entry/exit calibration.
Args:
spread: Spread series (residuals from cointegration regression)
Returns:
Half-life in days (float)
"""
X = spread.values.reshape(-1, 1)
lag = np.roll(X, 1)[1:]
y = np.diff(X, axis=0).flatten()
X = lag[1:]
# Fit: Δspread = λ * spread(-1) + ε
# Half-life = ln(2) / |λ|
coefs = np.linalg.lstsq(np.c_[np.ones(len(X)), X], y, rcond=None)[0]
lam = coefs[1]
if lam >= 0:
return np.inf # No mean reversion
half_life = -np.log(2) / lam
return max(half_life, 1) # Floor at 1 day
def test_cointegration_pair(series_a: pd.Series, series_b: pd.Series,
symbol_a: str, symbol_b: str,
p_threshold: float = 0.05) -> CointegrationPair | None:
"""
Full Engle-Granger two-step test for a single pair.
Args:
series_a: Price series for asset A
series_b: Price series for asset B
symbol_a: Ticker string for asset A
symbol_b: Ticker string for asset B
p_threshold: Maximum p-value to accept as cointegrated
Returns:
CointegrationPair object if cointegrated, None otherwise
"""
try:
hedge_ratio, spread = compute_spread_and_hedge_ratio(series_a, series_b)
except ValueError:
return None
adf_result = adf_test_spread(spread)
if adf_result["p_value"] > p_threshold:
return None
half_life = estimate_half_life(spread)
return CointegrationPair(
symbol_a=symbol_a,
symbol_b=symbol_b,
hedge_ratio=hedge_ratio,
adf_stat=adf_result["adf_statistic"],
p_value=adf_result["p_value"],
half_life=half_life
)
def screen_cointegration_pairs(df: pd.DataFrame, p_threshold: float = 0.05,
min_half_life: float = 2.0,
max_half_life: float = 60.0) -> list[CointegrationPair]:
"""
Screen all sector-matched pairs for cointegration.
This is the Tier 2 bottleneck — O(n²) in symbol count.
For production, parallelize across CPU cores or use GPU acceleration
for the regression step.
Args:
df: Price DataFrame (all symbols)
p_threshold: Maximum ADF p-value to accept
min_half_life: Minimum mean-reversion half-life (filter out too-fast pairs)
max_half_life: Maximum mean-reversion half-life (filter out stale pairs)
Returns:
List of CointegrationPair objects that pass all criteria
"""
symbols = list(df.columns)
valid_pairs = []
# Group symbols by sector for intra-sector pairing
# (Sector classification would come from a fundamental data source)
sector_groups = group_by_sector(symbols)
total_tested = 0
total_cointegrated = 0
for sector, sector_symbols in sector_groups.items():
for sym_a, sym_b in combinations(sector_symbols, 2):
total_tested += 1
series_a = df[sym_a].dropna()
series_b = df[sym_b].dropna()
pair = test_cointegration_pair(
series_a, series_b, sym_a, sym_b, p_threshold
)
if pair is None:
continue
# Filter by half-life bounds
if not (min_half_life <= pair.half_life <= max_half_life):
continue
total_cointegrated += 1
valid_pairs.append(pair)
print(f"Screened {total_tested} pairs — {total_cointegrated} cointegrated (p < {p_threshold})")
return sorted(valid_pairs, key=lambda p: p.p_value)
def group_by_sector(symbols: list) -> dict:
"""
Placeholder sector classification.
In production, integrate with a fundamental data provider
(e.g., Bloomberg sector classification, SEC filings).
"""
# Example mapping for common large-cap pairs
sector_map = {
"Technology": ["AAPL.US", "MSFT.US", "GOOGL.US", "META.US", "NVDA.US"],
"Finance": ["JPM.US", "BAC.US", "GS.US", "MS.US", "WFC.US"],
"Energy": ["XOM.US", "CVX.US", "COP.US", "SLB.US"],
"Healthcare": ["JNJ.US", "UNH.US", "PFE.US", "ABBV.US", "MRK.US"],
}
groups = {sector: [] for sector in sector_map}
groups["Other"] = []
for sym in symbols:
assigned = False
for sector, members in sector_map.items():
if sym in members:
groups[sector].append(sym)
assigned = True
break
if not assigned:
groups["Other"].append(sym)
return {k: v for k, v in groups.items() if len(v) >= 2}
Interpretation: What the Statistics Mean
| Metric | Formula | Trading interpretation |
|---|---|---|
| Hedge ratio (β) | OLS coefficient: A = β·B + ε | Short β shares of B per share of A |
| ADF statistic | Augmented Dickey-Fuller | More negative = stronger mean reversion |
| p-value | Probability of observing this ADF under null of unit root | p < 0.05: reject null, accept cointegration |
| Half-life | ln(2) / |λ| from OU process | Days for spread to revert 50% to mean |
A practical filter stack: accept pairs with p < 0.05, half-life between 5 and 30 days, and hedge ratio between 0.5 and 2.0 (pairs with extreme hedge ratios require disproportionately large capital allocation).
Tier 3: Rolling Window Validation — Detecting Structural Breaks
Static cointegration tests from a single historical window are insufficient for live trading. Markets are non-stationary; structural breaks invalidate past relationships. A robust system uses rolling window validation to monitor pair health continuously.
class RollingCointegrationMonitor:
"""
Monitors cointegration health over rolling windows.
Triggers alerts when the spread breaks down or half-life drifts.
"""
def __init__(self, pair: CointegrationPair, window_size: int = 252):
"""
Args:
pair: CointegrationPair object from screening
window_size: Trading days for rolling analysis window
"""
self.pair = pair
self.window_size = window_size
self.adf_history = []
self.half_life_history = []
self.break_alerts = []
def update(self, series_a: pd.Series, series_b: pd.Series, current_date: pd.Timestamp) -> dict:
"""
Update rolling statistics with new data point.
Args:
series_a: Full price series for asset A
series_b: Full price series for asset B
current_date: Current observation date
Returns:
Status dict with current metrics and any alerts
"""
# Rolling window
window_a = series_a[series_a.index <= current_date].tail(self.window_size)
window_b = series_b[series_b.index <= current_date].tail(self.window_size)
if len(window_a) < self.window_size * 0.8:
return {"status": "insufficient_data"}
try:
hedge_ratio, spread = compute_spread_and_hedge_ratio(window_a, window_b)
adf_result = adf_test_spread(spread)
half_life = estimate_half_life(spread)
except Exception:
return {"status": "calculation_error"}
self.adf_history.append({
"date": current_date,
"adf_stat": adf_result["adf_statistic"],
"p_value": adf_result["p_value"]
})
self.half_life_history.append({
"date": current_date,
"half_life": half_life
})
alerts = []
# Alert: ADF p-value crossed threshold (breakdown)
if adf_result["p_value"] > 0.10:
alerts.append({
"type": "breakdown",
"date": current_date,
"message": f"ADF p-value {adf_result['p_value']:.4f} exceeds 0.10 — pair may be broken"
})
self.break_alerts.append(alerts[-1])
# Alert: Half-life drifted significantly
if len(self.half_life_history) > 20:
recent_avg_hl = np.mean([h["half_life"] for h in self.half_life_history[-20:]])
if abs(half_life - recent_avg_hl) / recent_avg_hl > 0.5:
alerts.append({
"type": "half_life_shift",
"date": current_date,
"message": f"Half-life shifted to {half_life:.1f}d from recent avg {recent_avg_hl:.1f}d"
})
return {
"status": "ok",
"date": current_date,
"hedge_ratio": hedge_ratio,
"adf_stat": adf_result["adf_statistic"],
"p_value": adf_result["p_value"],
"half_life": half_life,
"alerts": alerts
}
⚠️ Engineering note: This monitor class maintains in-memory history. For production deployment tracking dozens of pairs simultaneously, persist state to a time-series database (InfluxDB, TimescaleDB) or use a sliding window buffer with periodic checkpointing.
The Spread Z-Score Trading Signal
With cointegration confirmed and hedge ratios established, the trading signal reduces to a standardized measure of spread deviation: the Z-Score.
Definition
Spread = A − β·B
Z-Score = (spread − rolling_mean) / rolling_std
The Z-Score tells us, in units of standard deviation, how far the current spread has deviated from its historical mean. A Z-Score of +2.5 means the spread is 2.5 standard deviations above its mean — a high-probability mean-reversion candidate.
Entry and Exit Rules
| Condition | Action |
|---|---|
| Z-Score > +2.0 | Short the spread: Short A, Long β·B |
| Z-Score < −2.0 | Long the spread: Long A, Short β·B |
| Z-Score crosses 0.0 | Close position (take profit or stop loss) |
| Z-Score > +3.0 | Add to short position (spread is extended) |
| ADF p-value > 0.10 | Exit all positions — pair is broken |
class SpreadZScore:
"""
Computes rolling Z-Score of a cointegrated spread.
Uses exponential moving average for mean and standard deviation
to weight recent observations more heavily.
"""
def __init__(self, lookback: int = 20, ema_span: int = 20):
"""
Args:
lookback: Historical window for mean/std estimation
ema_span: Exponential moving average span
"""
self.lookback = lookback
self.ema_span = ema_span
self.spread_history = []
def compute(self, spread: float) -> float:
"""
Update rolling Z-Score with a new spread observation.
Args:
spread: Current spread value
Returns:
Z-Score (float)
"""
self.spread_history.append(spread)
if len(self.spread_history) < self.lookback:
return np.nan
recent = pd.Series(self.spread_history[-self.lookback:])
# EMA-based mean and std for smoother, more responsive estimates
mean = recent.ewm(span=self.ema_span).mean().iloc[-1]
std = recent.ewm(span=self.ema_span).std().iloc[-1]
if std < 1e-8:
return np.nan
z_score = (spread - mean) / std
return z_score
def get_signal(self, z_score: float, entry_threshold: float = 2.0,
exit_threshold: float = 0.0, add_threshold: float = 3.0) -> str:
"""
Convert Z-Score to trading signal.
Args:
z_score: Current Z-Score
entry_threshold: Z-Score level for initial entry
exit_threshold: Z-Score level for exit
add_threshold: Z-Score level for adding to position
Returns:
Signal string: 'long', 'short', 'close', 'add_long', 'add_short', 'neutral'
"""
if np.isnan(z_score):
return "neutral"
if z_score > add_threshold:
return "add_short"
elif z_score < -add_threshold:
return "add_long"
elif z_score > entry_threshold:
return "short"
elif z_score < -entry_threshold:
return "long"
elif abs(z_score) < exit_threshold:
return "close"
else:
return "neutral"
Real-Time Spread Monitoring with WebSocket Streams
The strategy is only as good as the latency of the spread signal. A pairs trading system that updates hourly is vulnerable to adverse selection — faster market participants will have already traded the deviation before your signal fires. Real-time WebSocket feeds are essential for intraday pairs trading.
import json
import time
import random
import threading
import websocket
from collections import deque
class RealTimeSpreadMonitor:
"""
Real-time spread monitor using WebSocket price feeds.
Maintains rolling Z-Score and triggers trading signals.
⚠️ For production HFT workloads, consider aiohttp/asyncio with
a C++ or Rust co-processor for Z-Score computation.
This implementation prioritizes clarity for backtesting and
medium-frequency trading.
"""
PING_INTERVAL = 25 # seconds
MAX_RECONNECT_ATTEMPTS = 10
BASE_RECONNECT_DELAY = 1.0 # seconds
MAX_RECONNECT_DELAY = 60.0 # seconds
def __init__(self, pair: CointegrationPair, ws_url: str, api_key: str,
zscore_lookback: int = 20):
"""
Args:
pair: CointegrationPair with hedge ratio
ws_url: TickDB WebSocket endpoint
api_key: API authentication key
zscore_lookback: Number of observations for Z-Score rolling window
"""
self.pair = pair
self.ws_url = f"{ws_url}?api_key={api_key}"
self.api_key = api_key
self.zscore_lookback = zscore_lookback
self.spread_calculator = SpreadZScore(lookback=zscore_lookback)
self.latest_prices = {}
self.spread_deque = deque(maxlen=zscore_lookback * 2)
self.ws = None
self.reconnect_attempts = 0
self.is_running = False
self.last_signal = "neutral"
self._price_lock = threading.Lock()
self._alert_callbacks = []
def on_message(self, ws, message):
"""Handle incoming WebSocket messages."""
try:
data = json.loads(message)
if data.get("type") == "ping":
# Heartbeat response — maintain connection
ws.send(json.dumps({"cmd": "pong"}))
return
if data.get("channel") == "kline":
self._process_kline(data)
elif data.get("channel") == "trades":
self._process_trade(data)
except json.JSONDecodeError:
pass
except Exception as e:
print(f"Error processing message: {e}")
def _process_kline(self, data: dict):
"""Process 1-minute kline updates."""
try:
symbol = data.get("symbol")
kline = data.get("data", {})
close_price = float(kline.get("c"))
with self._price_lock:
self.latest_prices[symbol] = close_price
# Compute spread when both legs have prices
self._update_spread()
except (KeyError, ValueError) as e:
print(f"Kline parse error: {e}")
def _process_trade(self, data: dict):
"""Process tick-level trade updates for lowest latency."""
try:
symbol = data.get("symbol")
price = float(data.get("price", 0))
if price <= 0:
return
with self._price_lock:
self.latest_prices[symbol] = price
self._update_spread()
except (KeyError, ValueError) as e:
print(f"Trade parse error: {e}")
def _update_spread(self):
"""Compute spread and Z-Score from current prices."""
sym_a = self.pair.symbol_a
sym_b = self.pair.symbol_b
with self._price_lock:
if sym_a not in self.latest_prices or sym_b not in self.latest_prices:
return
price_a = self.latest_prices[sym_a]
price_b = self.latest_prices[sym_b]
# Spread = A - β * B
spread = price_a - self.pair.hedge_ratio * price_b
self.spread_deque.append(spread)
if len(self.spread_deque) < self.zscore_lookback:
return
# Compute rolling Z-Score
recent_spreads = list(self.spread_deque)
z_score = self.spread_calculator.compute(spread)
if np.isnan(z_score):
return
new_signal = self.spread_calculator.get_signal(z_score)
# Trigger alert callback on signal change
if new_signal != self.last_signal:
self.last_signal = new_signal
self._trigger_alert(new_signal, z_score, spread, price_a, price_b)
def _trigger_alert(self, signal: str, z_score: float, spread: float,
price_a: float, price_b: float):
"""Dispatch alerts to registered callbacks."""
alert = {
"signal": signal,
"z_score": z_score,
"spread": spread,
f"{self.pair.symbol_a}_price": price_a,
f"{self.pair.symbol_b}_price": price_b,
"hedge_ratio": self.pair.hedge_ratio,
"timestamp": time.time()
}
for callback in self._alert_callbacks:
try:
callback(alert)
except Exception as e:
print(f"Alert callback error: {e}")
def register_alert_callback(self, callback):
"""Register a function to be called on signal changes."""
self._alert_callbacks.append(callback)
def on_error(self, ws, error):
"""Handle WebSocket errors."""
print(f"WebSocket error: {error}")
def on_close(self, ws, close_status_code, close_msg):
"""Handle connection close with exponential backoff reconnection."""
print(f"Connection closed: {close_status_code} — {close_msg}")
if self.is_running and self.reconnect_attempts < self.MAX_RECONNECT_ATTEMPTS:
self._schedule_reconnect()
def on_open(self, ws):
"""Subscribe to price channels on connection open."""
print(f"Connected to {self.ws_url}")
self.reconnect_attempts = 0
# Subscribe to 1-minute klines for both legs of the pair
subscribe_msg = {
"cmd": "subscribe",
"params": {
"channels": [
{"channel": "kline", "symbol": self.pair.symbol_a, "interval": "1m"},
{"channel": "kline", "symbol": self.pair.symbol_b, "interval": "1m"}
]
}
}
ws.send(json.dumps(subscribe_msg))
print(f"Subscribed to {self.pair.symbol_a} and {self.pair.symbol_b}")
def _schedule_reconnect(self):
"""Exponential backoff with jitter to prevent thundering herd."""
delay = min(self.BASE_RECONNECT_DELAY * (2 ** self.reconnect_attempts),
self.MAX_RECONNECT_DELAY)
# Add jitter: random.uniform(0, delay * 0.1)
jitter = random.uniform(0, delay * 0.1)
reconnect_delay = delay + jitter
self.reconnect_attempts += 1
print(f"Reconnecting in {reconnect_delay:.1f}s (attempt {self.reconnect_attempts})")
threading.Timer(reconnect_delay, self._connect).start()
def _connect(self):
"""Establish WebSocket connection with heartbeat setup."""
self.ws = websocket.WebSocketApp(
self.ws_url,
on_message=self.on_message,
on_error=self.on_error,
on_close=self.on_close,
on_open=self.on_open
)
# Run with ping interval for keepalive
self.ws.run_forever(
ping_interval=self.PING_INTERVAL,
ping_timeout=10
)
def start(self):
"""Start the real-time monitor."""
self.is_running = True
self._connect()
def stop(self):
"""Stop the monitor and close the connection."""
self.is_running = False
if self.ws:
self.ws.close()
# Example alert callback: Slack webhook notification
def slack_alert_callback(alert: dict):
"""Forward trading signals to Slack."""
import os
import requests
webhook_url = os.environ.get("SLACK_WEBHOOK_URL")
if not webhook_url:
return
signal = alert["signal"]
z_score = alert["z_score"]
sym_a = alert.get(f"{signal.split('_')[0] if '_' in signal else 'pair'}_symbol_a", "A")
if signal in ("long", "short", "close", "add_long", "add_short"):
message = {
"text": f"⚠️ *Pairs Trading Signal*: {signal.upper()}",
"attachments": [{
"fields": [
{"title": "Z-Score", "value": f"{z_score:.2f}", "short": True},
{"title": "Spread", "value": f"{alert['spread']:.4f}", "short": True}
]
}]
}
try:
requests.post(webhook_url, json=message, timeout=5)
except requests.exceptions.RequestException:
pass
# Example usage
if __name__ == "__main__":
# Initialize with a pre-screened pair (from Tier 2 output)
example_pair = CointegrationPair(
symbol_a="XOM.US",
symbol_b="CVX.US",
hedge_ratio=0.87,
adf_stat=-3.82,
p_value=0.012,
half_life=12.3
)
api_key = os.environ.get("TICKDB_WS_API_KEY")
if not api_key:
raise ValueError("Set TICKDB_WS_API_KEY environment variable")
monitor = RealTimeSpreadMonitor(
pair=example_pair,
ws_url="wss://api.tickdb.ai/ws",
api_key=api_key,
zscore_lookback=20
)
# Register alert callbacks
monitor.register_alert_callback(slack_alert_callback)
monitor.register_alert_callback(
lambda a: print(f"[SIGNAL] {a['signal']} @ Z={a['z_score']:.2f}")
)
print("Starting real-time spread monitor for XOM-CVX...")
monitor.start()
⚠️ Engineering note: The WebSocket reconnect logic implements exponential backoff with a cap at 60 seconds. Rate-limit error codes from the API (code 3001) should trigger an immediate backoff regardless of reconnect state — check Retry-After headers and pause subscription attempts accordingly. This implementation uses threading.Timer for reconnection scheduling; production systems should use an async event loop (asyncio or trio) for cleaner concurrency management.
Backtest Framework and Performance Metrics
A pairs trading strategy requires rigorous backtesting before deployment. The code below provides a complete backtest engine with proper transaction cost modeling.
import numpy as np
from dataclasses import dataclass
from typing import Literal
@dataclass
class Trade:
entry_date: pd.Timestamp
exit_date: pd.Timestamp | None
direction: Literal["long", "short"]
entry_spread: float
exit_spread: float | None
pnl: float | None
hedge_ratio: float
class PairsBacktester:
"""
Backtest engine for cointegration-based pairs trading.
Accounts for:
- Round-trip commission (bid-ask + exchange fees)
- Slippage model (percentage of spread)
- Overnight funding costs (if holding through close)
- Spread crossing probability
"""
def __init__(self, commission_bps: float = 2.0, slippage_bps: float = 1.0,
spread_atr_multiplier: float = 1.5):
"""
Args:
commission_bps: Round-trip commission in basis points
slippage_bps: Slippage assumption in basis points
spread_atr_multiplier: ATR multiple for dynamic position sizing
"""
self.commission_bps = commission_bps
self.slippage_bps = slippage_bps
self.spread_atr_multiplier = spread_atr_multiplier
self.trades = []
self.current_position = None
def run(self, df: pd.DataFrame, pair: CointegrationPair,
entry_threshold: float = 2.0, exit_threshold: float = 0.0,
stop_threshold: float = 3.5) -> dict:
"""
Run backtest over historical price data.
Args:
df: DataFrame with price columns for both pair symbols
pair: CointegrationPair with hedge ratio
entry_threshold: Z-Score for entry
exit_threshold: Z-Score for exit
stop_threshold: Z-Score for stop-loss
Returns:
Dictionary with performance metrics
"""
zscore = SpreadZScore(lookback=20)
# Compute daily spreads
spreads = df[pair.symbol_a] - pair.hedge_ratio * df[pair.symbol_b]
spreads = spreads.dropna()
# Rolling Z-Score
z_scores = pd.Series(index=spreads.index, dtype=float)
for i, (date, spread) in enumerate(spreads.items()):
z_scores[date] = zscore.compute(spread)
# Walk-forward simulation
for date, z in z_scores.items():
if np.isnan(z):
continue
signal = zscore.get_signal(z, entry_threshold, exit_threshold)
if self.current_position is None:
if signal == "short":
self._open_position(date, z, "short", spreads[date], pair)
elif signal == "long":
self._open_position(date, z, "long", spreads[date], pair)
else:
if signal == "close" or abs(z) < exit_threshold:
self._close_position(date, z, spreads[date], pair)
elif signal == "stop_loss" or z > stop_threshold or z < -stop_threshold:
self._close_position(date, z, spreads[date], pair, forced=True)
# Close any open positions at end of series
if self.current_position:
final_spread = spreads.iloc[-1]
final_date = spreads.index[-1]
self._close_position(final_date, z_scores.iloc[-1], final_spread, pair)
return self._compute_metrics()
def _open_position(self, date: pd.Timestamp, z_score: float,
direction: Literal["long", "short"],
spread: float, pair: CointegrationPair):
"""Open a spread position with transaction costs."""
cost = abs(spread) * (self.commission_bps + self.slippage_bps) / 10000
self.current_position = Trade(
entry_date=date,
exit_date=None,
direction=direction,
entry_spread=spread,
exit_spread=None,
pnl=None,
hedge_ratio=pair.hedge_ratio
)
# Deduct entry cost
self.trades.append(self.current_position)
def _close_position(self, date: pd.Timestamp, z_score: float,
spread: float, pair: CointegrationPair, forced: bool = False):
"""Close position and record P&L."""
if not self.current_position:
return
pos = self.current_position
# Entry cost
entry_cost = abs(pos.entry_spread) * self.commission_bps / 10000
# Exit cost
exit_cost = abs(spread) * (self.commission_bps + self.slippage_bps) / 10000
# P&L: directional spread movement
if pos.direction == "short":
pnl = pos.entry_spread - spread # Short spread: profit when spread narrows
else:
pnl = spread - pos.entry_spread # Long spread: profit when spread widens
pos.exit_spread = spread
pos.exit_date = date
pos.pnl = pnl - entry_cost - exit_cost
if forced:
print(f"[{date.date()}] Forced exit: Z={z_score:.2f}, PnL={pos.pnl:.4f}")
self.current_position = None
def _compute_metrics(self) -> dict:
"""Compute performance statistics from closed trades."""
closed = [t for t in self.trades if t.pnl is not None]
if not closed:
return {"error": "No closed trades"}
pnls = [t.pnl for t in closed]
directions = [t.direction for t in closed]
return {
"total_trades": len(closed),
"winning_trades": sum(1 for p in pnls if p > 0),
"losing_trades": sum(1 for p in pnls if p <= 0),
"win_rate": sum(1 for p in pnls if p > 0) / len(closed),
"average_win": np.mean([p for p in pnls if p > 0]) if any(p > 0 for p in pnls) else 0,
"average_loss": np.mean([p for p in pnls if p <= 0]) if any(p <= 0 for p in pnls) else 0,
"profit_factor": abs(sum(p for p in pnls if p > 0) / sum(p for p in pnls if p < 0)) if any(p < 0 for p in pnls) else np.inf,
"total_pnl": sum(pnls),
"sharpe_ratio": np.mean(pnls) / np.std(pnls) * np.sqrt(252) if np.std(pnls) > 0 else 0,
"max_drawdown": self._max_drawdown(pnls),
"avg_holding_days": np.mean([(t.exit_date - t.entry_date).days for t in closed])
}
def _max_drawdown(self, pnls: list) -> float:
"""Compute maximum drawdown from cumulative P&L series."""
cumulative = np.cumsum(pnls)
running_max = np.maximum.accumulate(cumulative)
drawdowns = running_max - cumulative
return np.max(drawdowns)
Backtest limitations: Results above represent historical simulation and do not guarantee future performance. Key assumptions: 2 bps round-trip commission, 1 bps slippage, fixed hedge ratio (hedge ratio drift is not modeled), no liquidity exhaustion during rapid mean reversion. A minimum of 50 closed trades is recommended before drawing statistical conclusions.
Deployment Guide by User Segment
| Segment | Deployment approach | Data tier |
|---|---|---|
| Individual trader | Run backtest locally; deploy real-time monitor on a VPS ($10–20/month) | Free tier: 1-minute kline data, 1,000 API calls/day |
| Quant team | Backtest on cloud notebook (AWS SageMaker); deploy monitor as containerized microservice on ECS | Professional tier: daily OHLCV, real-time klines, WebSocket streams |
| Institutional | Full historical backtest across 10+ years of daily data; live paper-trading validation for 30 days before production | Enterprise tier: extended historical data, dedicated support, SLA guarantees |
For backtesting with 10+ years of cleaned, aligned US equity OHLCV data, contact TickDB's enterprise team at enterprise@tickdb.ai. The /v1/market/kline endpoint supports historical windows across the full 10+ year range required for multi-cycle statistical arbitrage validation.
Closing
The XOM-CVX divergence in February 2024 was not unique. Every quarter, energy pairs, financials, and technology sector cousins drift apart on earnings surprises, macro shocks, and sector rotations — then mean-revert. The traders who captured that signal were not guessing. They were running a system: cointegration screening, rolling window validation, Z-Score thresholds, and real-time spread monitoring.
Cointegration is not a static property. It is a living relationship that degrades when the underlying economic linkages change. A pairs trading system that does not monitor its own health — that does not alert when ADF p-values drift above 0.10 or half-lives shift by 50% — will bleed slowly until a catastrophic structural break wipes out months of accumulated edge.
The code in this article provides the complete pipeline: from screening 3,000 stocks down to validated pairs, through rolling window drift detection, to real-time WebSocket spread monitoring with Z-Score alerting.
The edge is not in the pair. It is in the system.
Next Steps
If you want to run this strategy yourself:
- Sign up at tickdb.ai (free, no credit card required)
- Generate an API key in the dashboard
- Set the
TICKDB_WS_API_KEYenvironment variable - Copy the screening code and real-time monitor from this article
If you need 10+ years of historical OHLCV data for cross-cycle pair validation:
Reach out to enterprise@tickdb.ai for institutional data plans covering US equities, HK equities, and crypto markets with the kline and trades endpoints.
If you are building a production pairs trading system:
Consider the tickdb-market-data SKILL available on ClawHub for your AI-assisted coding environment. It provides typed SDK wrappers around TickDB's WebSocket and REST APIs with built-in reconnection logic and error handling.
This article does not constitute investment advice. Statistical arbitrage strategies involve significant risk including but not limited to: model misspecification, structural breaks in cointegration relationships, execution costs, and adverse selection from faster market participants. Past performance does not guarantee future results. Always conduct thorough out-of-sample validation before committing capital.