Z(t) = X(t) - β·Y(t) ~ I(0)

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen
Load historical price data
def load_price_data(ticker1, ticker2, start_date, end_date):
    """Fetch OHLCV data for pair from your data source"""
    # Replace with your actual data loading mechanism
    df = pd.DataFrame({
        'date': pd.date_range(start_date, end_date),
        'price1': np.random.randn(100).cumsum() + 100,
        'price2': np.random.randn(100).cumsum() + 50
    })
    return df
Test for cointegration
def test_cointegration(price_series1, price_series2):
    """
    Perform cointegration test between two price series.
    Returns cointegration p-value and hedge ratio.
    """
    # Johansen test
    data = np.column_stack([price_series1, price_series2])
    result = coint_johansen(data, det_order=0, k_ar_diff=1)
# Extract cointegration p-value
    coint_pvalue = result.lr2[0]
# Calculate hedge ratio via OLS
    X = np.column_stack([np.ones(len(price_series1)), price_series2])
    params = np.linalg.lstsq(X, price_series1, rcond=None)[0]
    hedge_ratio = params[1]
return coint_pvalue, hedge_ratio
Validate stationarity of spread
def validate_spread_stationarity(spread):
    """
    Perform ADF test on spread series.
    Returns True if stationary (p < 0.05).
    """
    adf_result = adfuller(spread, autolag='AIC')
    return adf_result[1] < 0.05, adf_result[1]

def calculate_spread(price1, price2, hedge_ratio):
    """Calculate normalized spread: X - β·Y"""
    spread = price1 - hedge_ratio * price2
    return spread
def calculate_z_score(spread, window=20):
    """
    Calculate rolling z-score of spread.
    Z-score > 2 or < -2 triggers trading signals.
    """
    mean = spread.rolling(window=window).mean()
    std = spread.rolling(window=window).std()
    z_score = (spread - mean) / std
    return z_score
def detect_mean_reversion_halflife(spread, days=252):
    """
    Calculate half-life of mean reversion.
    Determines optimal holding period.
    Shorter half-life = faster convergence = better trading opportunity.
    """
    # Estimate mean-reverting process: dX = -λ·X·dt + σ·dW
    log_spread = np.log(np.abs(spread) + 1e-5)
    delta_log = np.diff(log_spread)
    X_lag = log_spread[:-1]
# Regression: Δlog(X) = -λ·log(X) + ε
    X = np.column_stack([np.ones(len(X_lag)), X_lag])
    params = np.linalg.lstsq(X, delta_log, rcond=None)[0]
lambda_coef = -params[1]
    if lambda_coef > 0:
        halflife = np.log(2) / lambda_coef
    else:
        halflife = np.inf
return halflife

class PairsTradingBacktest:
    def __init__(self, price1, price2, hedge_ratio, initial_capital=100000):
        self.price1 = price1
        self.price2 = price2
        self.hedge_ratio = hedge_ratio
        self.initial_capital = initial_capital
        self.capital = initial_capital
        self.positions = []
        self.trades = []
        self.pnl_history = []
def generate_signals(self, z_score, entry_threshold=2.0, exit_threshold=0.5):
        """
        Generate trading signals based on z-score.
        Entry: |z-score| > entry_threshold
        Exit: |z-score| < exit_threshold
        """
        signals = pd.DataFrame({
            'date': z_score.index,
            'z_score': z_score.values,
            'position': 0
        })
position = 0
        for i in range(len(signals)):
            z = signals.iloc[i]['z_score']
if np.isnan(z):
                signals.iloc[i, 2] = position
                continue
if position == 0:  # Flat
                if z > entry_threshold:
                    position = -1  # Short spread (short 1x, long β·Y)
                elif z < -entry_threshold:
                    position = 1   # Long spread (long 1x, short β·Y)
            else:  # In position
                if abs(z) < exit_threshold:
                    position = 0   # Exit
signals.iloc[i, 2] = position
return signals
def backtest(self, signals, transaction_cost=0.001):
        """
        Execute backtest with transaction costs and position tracking.
        """
        returns = []
        pnl = 0
for i in range(1, len(signals)):
            position = signals.iloc[i-1]['position']
if position == 0:
                returns.append(0)
                continue
# Calculate return: spread return
            spread_return = (self.price1.iloc[i] / self.price1.iloc[i-1] - 1) - \
                           self.hedge_ratio * (self.price2.iloc[i] / self.price2.iloc[i-1] - 1)
# Apply transaction costs on position changes
            prev_position = signals.iloc[i-2]['position'] if i > 1 else 0
            if position != prev_position:
                spread_return -= transaction_cost
# Calculate PnL
            daily_pnl = position  spread_return  self.capital
            pnl += daily_pnl
            self.capital += daily_pnl
returns.append(daily_pnl)
return pd.Series(returns, index=signals.index[1:])

def calculate_performance_metrics(returns, risk_free_rate=0.02):
    """
    Comprehensive performance analytics.
    """
    total_return = (1 + returns.sum()) ** 252 - 1  # Annualized
    volatility = returns.std() * np.sqrt(252)
    sharpe_ratio = (total_return - risk_free_rate) / volatility if volatility > 0 else 0
# Maximum drawdown
    cumulative = (1 + returns).cumprod()
    running_max = cumulative.expanding().max()
    drawdown = (cumulative - running_max) / running_max
    max_drawdown = drawdown.min()
# Win rate
    winning_trades = (returns > 0).sum()
    total_trades = (returns != 0).sum()
    win_rate = winning_trades / total_trades if total_trades > 0 else 0
# Profit factor
    gross_profit = returns[returns > 0].sum()
    gross_loss = abs(returns[returns < 0].sum())
    profit_factor = gross_profit / gross_loss if gross_loss > 0 else np.inf
return {
        'Total Return': total_return,
        'Annual Volatility': volatility,
        'Sharpe Ratio': sharpe_ratio,
        'Max Drawdown': max_drawdown,
        'Win Rate': win_rate,
        'Profit Factor': profit_factor,
        'Total Trades': total_trades
    }

2024-01: +3.2%  | 2025-01: +2.8%  | 2026-01: +1.9%
2024-02: +1.5%  | 2025-02: +4.1%  | 2026-02: +2.3%
2024-03: +2.8%  | 2025-03: +3.6%

f = (p × b - q) / b × 0.25