/09_model_based_features/14_panel_features.py

# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: tags,-all
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.19.3
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Panel Features: Pairwise and Cross-Sectional Transforms
#
# **Docker image**: `ml4t`
#
# This notebook demonstrates panel-level temporal features: pairwise
# relationships (cointegration, Kalman hedge ratios, O-U half-life) and
# cross-sectional transforms (ranking, relative features, universe
# aggregation).
#
# Panel features matter because temporal features computed in isolation
# (volatility, momentum, regime probabilities) gain signal when placed
# in cross-sectional context — a 25% conditional volatility means
# different things for a utility stock and a biotech.
#
# **Learning Objectives**:
# - Rank temporal features cross-sectionally for relative positioning
# - Compute relative temporal features (vs. sector, vs. market)
# - Test for cointegration using Engle-Granger and Johansen methods
# - Estimate dynamic hedge ratios with Kalman filter
# - Compute Ornstein-Uhlenbeck half-life for mean-reversion speed
#
# **Book Reference**: Chapter 9, Section 9.6 (Cross-Sectional and Panel Features)
#
# **Prerequisites**: `04_kalman_filter` for Kalman filter mechanics,
# `11_hmm_regimes` for regime features used in universe aggregation.

# %%
"""Panel Features — pairwise and cross-sectional transforms for multi-asset temporal features."""

import warnings

warnings.filterwarnings("ignore", category=FutureWarning)

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl

# Kalman filter for dynamic hedge ratio
from filterpy.kalman import KalmanFilter
from IPython.display import display
from ml4t.engineer.features.cross_asset import (
    beta_to_market,
    co_integration_score,
    correlation_regime_indicator,
    rolling_correlation,
)
from sklearn.linear_model import LinearRegression

# Cointegration tests
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen

from data import load_etfs
from utils.paths import get_case_study_dir

# %% tags=["parameters"]
# Production defaults — Papermill injects overrides for CI
MAX_SYMBOLS = 0  # 0 = all symbols
START_DATE = "2015-01-01"
END_DATE = "2024-12-31"

# %% [markdown]
# ## Load Data
#
# We use ETF pairs from the ETF Universe dataset to demonstrate pairs trading.
# Classic cointegrated pairs include:
# - **GLD/SLV** - Gold vs Silver (precious metals)
# - **XLE/USO** - Energy sector vs Crude Oil
# - **QQQ/SMH** - Nasdaq 100 vs Semiconductors (tech exposure)

# %%
# Load ETF universe data
etf_data = load_etfs()
if MAX_SYMBOLS > 0:
    symbols = etf_data["symbol"].unique().sort().to_list()[:MAX_SYMBOLS]
    etf_data = etf_data.filter(pl.col("symbol").is_in(symbols))

# Filter date range
etf_data = etf_data.filter(
    pl.col("timestamp").is_between(pl.lit(START_DATE).str.to_date(), pl.lit(END_DATE).str.to_date())
)

# Define pairs to test for cointegration
# Primary pair: XLE/USO (energy sector vs crude oil) with structural relationship
# Fallback pair: SPY/QQQ (available in all datasets for testing)
available_symbols = etf_data["symbol"].unique().to_list()
if "XLE" in available_symbols and "USO" in available_symbols:
    PAIR_Y = "XLE"  # Energy Sector ETF
    PAIR_X = "USO"  # Crude Oil ETF
else:
    # Fallback to symbols available in test fixtures
    print("Using fallback pair SPY/QQQ (XLE/USO not available)")
    PAIR_Y = "SPY"
    PAIR_X = "QQQ"

# Extract prices for the pair
# %%
pair_y = (
    etf_data.filter(pl.col("symbol") == PAIR_Y)
    .select(["timestamp", "close"])
    .rename({"close": "Y"})
    .sort("timestamp")
)

pair_x = (
    etf_data.filter(pl.col("symbol") == PAIR_X)
    .select(["timestamp", "close"])
    .rename({"close": "X"})
    .sort("timestamp")
)

# Join on date (inner join to align dates)
df_pl = pair_y.join(pair_x, on="timestamp", how="inner")

# Convert to pandas for statsmodels/sklearn
df = df_pl.to_pandas().set_index("timestamp")
df.index = pd.DatetimeIndex(df.index)

print(f"Pair: {PAIR_Y}/{PAIR_X}")
print(f"Data: {len(df)} observations from {df.index.min()} to {df.index.max()}")
print(f"{PAIR_Y} range: ${df['Y'].min():.2f} - ${df['Y'].max():.2f}")
print(f"{PAIR_X} range: ${df['X'].min():.2f} - ${df['X'].max():.2f}")

# %% [markdown]
# ## Cointegration Testing
#
# Two standard approaches:
# - **Engle-Granger**: Two-step test - regress Y on X, then test residuals for stationarity
# - **Johansen**: System-based test for multiple cointegrating relationships
#
# We apply both tests - if they agree, we have stronger evidence of cointegration.

# %%
# Engle-Granger cointegration test
eg_stat, eg_pval, _ = coint(df["Y"], df["X"], trend="c")[:3]

print("=== Engle-Granger Cointegration Test ===")
print(f"Test statistic: {eg_stat:.4f}")
print(f"P-value: {eg_pval:.4f}")
print(f"Conclusion: {'Cointegrated' if eg_pval < 0.05 else 'Not cointegrated'} at 5% level")

# %%
# Johansen cointegration test
joh_result = coint_johansen(df[["Y", "X"]], det_order=0, k_ar_diff=1)
trace_stats = joh_result.lr1
crit_values = joh_result.cvt[:, 1]  # 95% critical values

print("\n=== Johansen Cointegration Test ===")
print(f"Trace statistic (r=0): {trace_stats[0]:.4f} (critical: {crit_values[0]:.4f})")
print(f"Trace statistic (r≤1): {trace_stats[1]:.4f} (critical: {crit_values[1]:.4f})")

# Reject if trace stat > critical value
johansen_sig = trace_stats[0] > crit_values[0]
print(f"Conclusion: {'Cointegrated' if johansen_sig else 'Not cointegrated'} at 5% level")

# Cointegrating vector (normalized)
coint_vector = joh_result.evec[:, 0]
hedge_ratio_joh = -coint_vector[1] / coint_vector[0]
print(f"\nImplied hedge ratio from Johansen: {hedge_ratio_joh:.4f}")

# %% [markdown]
# When both tests agree on cointegration, we have stronger evidence of a
# long-run equilibrium relationship. When they disagree, treat the pair
# with caution — the relationship may be fragile or sensitive to the
# sample period. The Johansen-implied hedge ratio often differs from OLS
# because the two methods weight observations differently.

# %% [markdown]
# ## Spread Construction
#
# ### Static Hedge Ratio (OLS)
#
# The simplest approach: regress Y on X to find the hedge ratio $\beta$.
# Spread $= Y - \beta X$

# %%
# OLS hedge ratio
ols = LinearRegression()
ols.fit(df[["X"]], df["Y"])
hedge_ratio_ols = ols.coef_[0]
intercept_ols = ols.intercept_

print(f"OLS Hedge Ratio: {hedge_ratio_ols:.4f}")
print(f"OLS Intercept: {intercept_ols:.4f}")

# Compute spread
df["spread_ols"] = df["Y"] - hedge_ratio_ols * df["X"]

# Test spread stationarity
adf_stat, adf_pval, _, _, _, _ = adfuller(df["spread_ols"], autolag="AIC")
print(f"\nSpread ADF test: stat={adf_stat:.4f}, p-value={adf_pval:.4f}")
print(f"Spread is {'stationary' if adf_pval < 0.05 else 'non-stationary'}")

# %% [markdown]
# ### Dynamic Hedge Ratio (Kalman Filter)
#
# The hedge ratio may change over time. A Kalman filter provides:
# - Adaptive estimation as new data arrives
# - Smooth transitions in the hedge ratio
# - Natural handling of non-stationarity in the relationship


# %%
def kalman_hedge_ratio(y: np.ndarray, x: np.ndarray) -> tuple:
    """
    Estimate time-varying hedge ratio using Kalman filter.

    Returns:
        hedge_ratios: Array of time-varying hedge ratios
        intercepts: Array of time-varying intercepts
    """
    n = len(y)

    # State: [intercept, hedge_ratio]
    kf = KalmanFilter(dim_x=2, dim_z=1)

    # Transition matrix (random walk for states)
    kf.F = np.eye(2)

    # Measurement function: y_t = [1, x_t] @ [alpha, beta]
    # We'll update H at each step

    # Initial state
    kf.x = np.array([[0.0], [1.0]])  # Start with hedge ratio = 1

    # Covariance matrices
    kf.P = np.eye(2) * 1.0  # Initial state covariance
    kf.R = np.array([[1e-3]])  # Measurement noise (small)
    kf.Q = np.eye(2) * 1e-5  # Process noise (slow-changing states)

    hedge_ratios = np.zeros(n)
    intercepts = np.zeros(n)

    for t in range(n):
        # Update measurement matrix for current x value
        kf.H = np.array([[1.0, x[t]]])

        # Predict
        kf.predict()

        # Update with measurement y_t
        kf.update(np.array([[y[t]]]))

        # Store estimates
        intercepts[t] = kf.x[0, 0]
        hedge_ratios[t] = kf.x[1, 0]

    return hedge_ratios, intercepts


# Apply Kalman filter
hedge_ratios_kf, intercepts_kf = kalman_hedge_ratio(df["Y"].values, df["X"].values)

df["hedge_ratio_kf"] = hedge_ratios_kf
df["spread_kf"] = df["Y"] - df["hedge_ratio_kf"] * df["X"]

print(f"Kalman hedge ratio range: [{hedge_ratios_kf.min():.4f}, {hedge_ratios_kf.max():.4f}]")
print(f"Final hedge ratio: {hedge_ratios_kf[-1]:.4f}")

# %% [markdown]
# ## Mean-Reversion Half-Life
#
# The half-life measures how quickly the spread reverts to its mean.
# Estimated from an AR(1) model: spread_t = α + ρ * spread_{t-1} + ε
# Half-life = -ln(2) / ln(ρ)


# %%
def estimate_half_life(spread: pd.Series) -> float:
    """Estimate mean-reversion half-life from spread series."""
    spread = spread.dropna()
    spread_lag = spread.shift(1).dropna()
    spread_diff = spread.diff().dropna()

    # Align indices
    spread_lag = spread_lag.iloc[:-1] if len(spread_lag) > len(spread_diff) else spread_lag

    # Regress spread_diff on spread_lag
    X = spread_lag.values.reshape(-1, 1)
    y = spread_diff.values[: len(X)]

    reg = LinearRegression()
    reg.fit(X, y)
    phi = reg.coef_[0]

    # Half-life = -ln(2) / ln(1 + phi)
    # For small phi: half-life ≈ -ln(2) / phi
    if phi >= 0:
        return np.inf  # No mean reversion
    half_life = -np.log(2) / phi
    return max(1, half_life)


half_life_ols = estimate_half_life(df["spread_ols"])
half_life_kf = estimate_half_life(df["spread_kf"])

print(f"Half-life (OLS spread): {half_life_ols:.1f} days")
print(f"Half-life (Kalman spread): {half_life_kf:.1f} days")

# %% [markdown]
# The half-life estimate uses full-sample OLS on the spread — in a
# walk-forward context, only past data should be used. The Kalman spread
# typically yields shorter half-lives because the adaptive hedge ratio
# removes low-frequency drift from the spread, leaving a faster-reverting residual.

# %% [markdown]
# ## Trading Signals: Z-Score and Bollinger Bands
#
# Entry/exit rules based on spread deviation from mean:
# - **Long spread** (buy Y, sell X): Z-score < -2 (spread too low)
# - **Short spread** (sell Y, buy X): Z-score > +2 (spread too high)
# - **Exit**: Z-score crosses zero (mean reversion complete)

# %%
# Z-score and Bollinger bands from Kalman spread
spread = df["spread_kf"]

lookback = max(int(2 * half_life_kf), 20)  # At least 20 days
rolling_mean = spread.rolling(window=lookback).mean()
rolling_std = spread.rolling(window=lookback).std()

df["z_score"] = (spread - rolling_mean) / rolling_std
df["spread_mean"] = rolling_mean
df["upper_band"] = rolling_mean + 2 * rolling_std
df["lower_band"] = rolling_mean - 2 * rolling_std

print(f"Lookback window: {lookback} days (2x half-life)")

# %% [markdown]
# The z-score normalizes the spread by its rolling mean and standard
# deviation — it is the natural signal for mean-reverting spreads because
# cointegration implies the spread is stationary around a fixed level.
# Entry at $\pm 2\sigma$ captures significant deviations while filtering noise.

# %%
# Trading signals: enter at ±2σ, exit at mean crossing
df["signal"] = 0
df.loc[df["z_score"] < -2, "signal"] = 1  # Long spread
df.loc[df["z_score"] > 2, "signal"] = -1  # Short spread

# Position tracking (hold until mean reversion)
signals = df["signal"].values
z_scores = df["z_score"].values
positions = [0] * len(df)
in_trade = 0

for i in range(1, len(df)):
    if signals[i] != 0:
        in_trade = signals[i]
    elif in_trade != 0:
        if (in_trade > 0 and z_scores[i] > 0) or (in_trade < 0 and z_scores[i] < 0):
            in_trade = 0
    positions[i] = in_trade

df["position"] = positions

print(f"Number of trades: {(df['signal'] != 0).sum()}")
print(f"Days in position: {(df['position'] != 0).sum()}")

# %% [markdown]
# ## Visualize Trading System

# %%
fig, axes = plt.subplots(4, 1, figsize=(14, 14), sharex=True)

# Price series
ax = axes[0]
ax.plot(df.index, df["Y"], label=PAIR_Y, linewidth=0.8)
ax.plot(df.index, df["X"], label=PAIR_X, linewidth=0.8)
ax.set_ylabel("Price ($)")
ax.set_title(f"{PAIR_Y}/{PAIR_X} Price Series")
ax.legend()

# Hedge ratio evolution
ax = axes[1]
ax.axhline(hedge_ratio_ols, color="#d62728", linestyle="--", label=f"OLS: {hedge_ratio_ols:.3f}")
ax.plot(df.index, df["hedge_ratio_kf"], label="Kalman", linewidth=0.8)
ax.set_ylabel("Hedge Ratio")
ax.set_title("Dynamic vs Static Hedge Ratio")
ax.legend()

# Spread with Bollinger bands
ax = axes[2]
ax.plot(df.index, df["spread_kf"], label="Spread", linewidth=0.8)
ax.plot(df.index, df["spread_mean"], color="#d62728", linestyle="--", linewidth=0.8, label="Mean")
ax.fill_between(
    df.index, df["lower_band"], df["upper_band"], alpha=0.2, color="gray", label="±2σ bands"
)
ax.set_ylabel("Spread")
ax.set_title("Spread with Bollinger Bands")
ax.legend()

# Z-score with entry/exit zones
ax = axes[3]
ax.plot(df.index, df["z_score"], linewidth=0.8)
ax.axhline(2, color="#d62728", linestyle="--", alpha=0.5, label="Short entry")
ax.axhline(-2, color="#1f77b4", linestyle="--", alpha=0.5, label="Long entry")
ax.axhline(0, color="black", linestyle="-", alpha=0.3)
ax.fill_between(df.index, -0.5, 0.5, alpha=0.1, color="gray", label="Exit zone")
ax.set_ylabel("Z-Score")
ax.set_title("Z-Score with Entry/Exit Thresholds")
ax.legend()

plt.tight_layout()
plt.show()

# %% [markdown]
# ## Illustrative Backtest
#
# A compact backtest demonstrates how cointegration features translate into
# tradable signals. This ignores transaction costs and slippage — see
# Chapters 18–19 for proper strategy evaluation with realistic cost models.

# %%
# Spread returns (long Y, short X when position = 1)
df["ret_y"] = df["Y"].pct_change()
df["ret_x"] = df["X"].pct_change()
df["spread_ret"] = df["position"].shift(1) * (df["ret_y"] - df["ret_x"])
df["spread_ret"] = df["spread_ret"].fillna(0)
df["cumulative_ret"] = (1 + df["spread_ret"]).cumprod() - 1

total_return = df["cumulative_ret"].iloc[-1]
ann_return = (1 + total_return) ** (252 / len(df)) - 1
volatility = df["spread_ret"].std() * np.sqrt(252)
sharpe = ann_return / volatility if volatility > 0 else 0
rolling_max = (1 + df["cumulative_ret"]).cummax()
max_dd = ((1 + df["cumulative_ret"]) / rolling_max - 1).min()

backtest_summary = pd.DataFrame(
    {
        "Total Return": [f"{total_return:.1%}"],
        "Ann. Return": [f"{ann_return:.1%}"],
        "Ann. Vol": [f"{volatility:.1%}"],
        "Sharpe": [f"{sharpe:.2f}"],
        "Max DD": [f"{max_dd:.1%}"],
    }
)
display(backtest_summary)

# %% [markdown]
# ## Summary: Primary Pair Analysis

# %%
summary_df = pd.DataFrame(
    {
        "Metric": [
            "Observations",
            "EG p-value",
            "Johansen sig.",
            "OLS hedge ratio",
            "Kalman hedge ratio (final)",
            "Half-life (days)",
        ],
        "Value": [
            f"{len(df):,}",
            f"{eg_pval:.4f}",
            str(johansen_sig),
            f"{hedge_ratio_ols:.4f}",
            f"{hedge_ratios_kf[-1]:.4f}",
            f"{half_life_kf:.1f}",
        ],
    }
)
display(summary_df)

# %% [markdown]
# ## Cointegration Screening: Multiple ETF Pairs
#
# In practice, we screen multiple pairs to find cointegrated relationships.
# Not all pairs that seem related are actually cointegrated.

# %%
# Screen multiple ETF pairs for cointegration
CANDIDATE_PAIRS = [
    ("GLD", "SLV", "Gold vs Silver"),
    ("XLE", "USO", "Energy vs Crude Oil"),
    ("QQQ", "SMH", "Nasdaq 100 vs Semis"),
    ("SPY", "VTI", "S&P 500 vs Total Mkt"),
    ("TLT", "IEF", "Long vs Mid Treasuries"),
    ("EEM", "VWO", "EM (iShares vs Vanguard)"),
]

screening_results = []

# %%
for pair_y, pair_x, description in CANDIDATE_PAIRS:
    y_data = (
        etf_data.filter(pl.col("symbol") == pair_y).select(["timestamp", "close"]).sort("timestamp")
    )
    x_data = (
        etf_data.filter(pl.col("symbol") == pair_x).select(["timestamp", "close"]).sort("timestamp")
    )

    if len(y_data) < 252 or len(x_data) < 252:
        continue

    y_df = y_data.rename({"close": "Y"})
    x_df = x_data.rename({"close": "X"})
    pair_df = y_df.join(x_df, on="timestamp", how="inner").to_pandas().set_index("timestamp")

    if len(pair_df) < 252:
        continue

    # Engle-Granger test
    eg_stat, eg_pval_pair, _ = coint(pair_df["Y"], pair_df["X"], trend="c")[:3]

    # Johansen test
    joh = coint_johansen(pair_df[["Y", "X"]], det_order=0, k_ar_diff=1)
    joh_sig = joh.lr1[0] > joh.cvt[0, 1]

    # Hedge ratio and half-life
    ols_temp = LinearRegression()
    ols_temp.fit(pair_df[["X"]], pair_df["Y"])
    hedge = ols_temp.coef_[0]
    hl = estimate_half_life(pair_df["Y"] - hedge * pair_df["X"])

    screening_results.append(
        {
            "Pair": f"{pair_y}/{pair_x}",
            "Description": description,
            "Obs": len(pair_df),
            "EG p-value": eg_pval_pair,
            "Johansen": "Yes" if joh_sig else "No",
            "Cointegrated": eg_pval_pair < 0.05,
            "Hedge Ratio": hedge,
            "Half-life (d)": hl,
        }
    )

# %%
# Display screening results as a formatted table
screening_df = pd.DataFrame(screening_results)
display_df = screening_df[
    ["Pair", "Description", "EG p-value", "Johansen", "Hedge Ratio", "Half-life (d)"]
].copy()
display_df["EG p-value"] = display_df["EG p-value"].map(lambda x: f"{x:.4f}")
display_df["Hedge Ratio"] = display_df["Hedge Ratio"].map(lambda x: f"{x:.3f}")
display_df["Half-life (d)"] = display_df["Half-life (d)"].map(
    lambda x: f"{x:.1f}" if x < 1000 else "N/A"
)
display(display_df)

# %% [markdown]
# Economic relatedness does not guarantee cointegration — several
# apparently related pairs fail the Engle-Granger test, and the two
# tests sometimes disagree. Half-lives range from days to months,
# reflecting the speed at which different pair relationships mean-revert.
# Only pairs with short half-lives (days to weeks) are practical for
# trading at daily frequency.

# %% [markdown]
# ## Save Pairs Trading Features for Downstream Chapters
#
# Spread signals and hedge ratios are consumed by:
# - Chapter 18: Strategy simulation (pairs trading backtest)
# - Chapter 19: Portfolio construction (market-neutral portfolios)

# %%
# Save pairs trading data to output directory for ALL cointegrated pairs
MODEL_DIR = get_case_study_dir("etfs") / "models" / "time_series"
MODEL_DIR.mkdir(parents=True, exist_ok=True)

# %% [markdown]
# ### Generate Signals for Cointegrated Pairs
#
# Package spread signals and hedge ratios for all cointegrated pairs,
# saving for downstream strategy evaluation in Chapters 18–19.


# %%
def generate_pair_signals(pair_y: str, pair_x: str) -> pl.DataFrame | None:
    """Generate Kalman-filtered spread signals for a single pair."""
    y_data = (
        etf_data.filter(pl.col("symbol") == pair_y).select(["timestamp", "close"]).sort("timestamp")
    )
    x_data = (
        etf_data.filter(pl.col("symbol") == pair_x).select(["timestamp", "close"]).sort("timestamp")
    )

    if len(y_data) < 252 or len(x_data) < 252:
        return None

    y_df = y_data.rename({"close": "Y"})
    x_df = x_data.rename({"close": "X"})
    pair_df = y_df.join(x_df, on="timestamp", how="inner").to_pandas().set_index("timestamp")
    pair_df.index = pd.DatetimeIndex(pair_df.index)

    if len(pair_df) < 252:
        return None

    # Kalman filter hedge ratio and spread
    kf_hedges, _ = kalman_hedge_ratio(pair_df["Y"].values, pair_df["X"].values)
    pair_df["spread_kf"] = pair_df["Y"] - kf_hedges * pair_df["X"]

    # Z-score from half-life based lookback
    hl = estimate_half_life(pair_df["spread_kf"])
    lb = max(int(2 * hl), 20)
    rm = pair_df["spread_kf"].rolling(window=lb).mean()
    rs = pair_df["spread_kf"].rolling(window=lb).std()
    pair_df["z_score"] = (pair_df["spread_kf"] - rm) / rs

    pair_df["signal"] = 0
    pair_df.loc[pair_df["z_score"] < -2, "signal"] = 1
    pair_df.loc[pair_df["z_score"] > 2, "signal"] = -1

    return pl.DataFrame(
        {
            "timestamp": pair_df.index.values,
            "pair": f"{pair_y}/{pair_x}",
            "asset_y": pair_y,
            "asset_x": pair_x,
            "hedge_ratio": kf_hedges,
            "spread": pair_df["spread_kf"].values,
            "z_score": pair_df["z_score"].values,
            "signal": pair_df["signal"].values,
        }
    ).drop_nulls()


# %%
# Apply to all cointegrated pairs
all_signals = []
cointegrated_pairs = [r for r in screening_results if r["Cointegrated"]]

print(f"Processing {len(cointegrated_pairs)} cointegrated pairs...")

for r in cointegrated_pairs:
    pair_y, pair_x = r["Pair"].split("/")
    signals = generate_pair_signals(pair_y, pair_x)
    if signals is not None:
        all_signals.append(signals)

# Also include the primary pair if not already included
if f"{PAIR_Y}/{PAIR_X}" not in [r["Pair"] for r in cointegrated_pairs]:
    primary_signals = pl.DataFrame(
        {
            "timestamp": df.index.values,
            "pair": f"{PAIR_Y}/{PAIR_X}",
            "asset_y": PAIR_Y,
            "asset_x": PAIR_X,
            "hedge_ratio": df["hedge_ratio_kf"].values,
            "spread": df["spread_kf"].values,
            "z_score": df["z_score"].values,
            "signal": df["signal"].values,
        }
    ).drop_nulls()
    all_signals.append(primary_signals)
    print(f"  {PAIR_Y}/{PAIR_X}: {len(primary_signals)} obs (primary)")

# %%
# Combine and save
if all_signals:
    combined_signals = pl.concat(all_signals)
    output_path = MODEL_DIR / "pairs_trading_signals.parquet"
    combined_signals.write_parquet(output_path)

    print(f"\n[OK] Saved pairs trading signals to {output_path}")
    print(f"  Shape: {combined_signals.shape}")
    print(f"  Pairs: {combined_signals['pair'].unique().to_list()}")
    print(
        f"  Date range: {combined_signals['timestamp'].min()} to {combined_signals['timestamp'].max()}"
    )

    # Per-pair summary
    print("\n  Per-pair summary:")
    summary = combined_signals.group_by("pair").agg(pl.len().alias("n_rows")).sort("pair")
    for row in summary.iter_rows(named=True):
        print(f"    {row['pair']}: {row['n_rows']:,} rows")
else:
    print("\nWARNING: No pairs trading signals generated")

# %% [markdown]
# ## ml4t-engineer: Cross-Asset Features as Polars Expressions
#
# The manual pairwise analysis above uses statsmodels and numpy for each pair.
# `ml4t-engineer` provides cross-asset features as Polars expressions for
# multi-pair pipelines: rolling correlation, beta, dispersion, and
# cointegration scoring.

# %%
# Demonstrate cross-asset features on a pair
gld = etf_data.filter(pl.col("symbol") == "GLD").select(["timestamp", "close"]).sort("timestamp")
slv = etf_data.filter(pl.col("symbol") == "SLV").select(["timestamp", "close"]).sort("timestamp")

pair_df = (
    gld.rename({"close": "gld_close"})
    .join(slv.rename({"close": "slv_close"}), on="timestamp", how="inner")
    .with_columns(
        gld_ret=pl.col("gld_close").pct_change(),
        slv_ret=pl.col("slv_close").pct_change(),
    )
    .drop_nulls()
)

# Cross-asset features in one pipeline
pair_features = pair_df.with_columns(
    corr=rolling_correlation("gld_ret", "slv_ret", window=60),
    beta=beta_to_market("slv_ret", "gld_ret", window=60),
    coint_score=co_integration_score("gld_close", "slv_close", window=120),
)

# Correlation regime indicators (returns dict)
corr_regime = correlation_regime_indicator("corr")
pair_features = pair_features.with_columns(**corr_regime)

print("=== ml4t-engineer: Cross-Asset Features (GLD/SLV) ===")
for col in ["corr", "beta", "coint_score"]:
    vals = pair_features[col].drop_nulls()
    _m, _s = vals.mean(), vals.std()
    print(
        f"  {col:<14}: mean={_m:.4f}, std={_s:.4f}"
        if _m is not None
        else f"  {col:<14}: insufficient data"
    )

# %% [markdown]
# The library expressions compose naturally in `with_columns()`, making it
# straightforward to compute cross-asset features across many pairs in a
# single Polars pipeline. `correlation_regime_indicator()` adds binary
# regime flags based on rolling correlation levels — useful for
# conditional feature engineering.

# %% [markdown]
# ## Cross-Sectional Ranking of Temporal Features
#
# Raw temporal features (momentum, conditional volatility, regime
# probabilities) vary in scale across assets. Cross-sectional ranking
# converts them to universe-relative signals — essential for panel models
# that predict cross-sectional returns.

# %%
# Compute temporal features for a multi-asset ETF universe
RANK_SYMBOLS = ["SPY", "QQQ", "IWM", "EFA", "EEM", "TLT", "GLD", "XLE", "XLF", "XLV"]
rank_symbols = [s for s in RANK_SYMBOLS if s in available_symbols]

panel = (
    etf_data.filter(pl.col("symbol").is_in(rank_symbols))
    .select(["timestamp", "symbol", "close"])
    .sort(["symbol", "timestamp"])
)

# Temporal features: momentum (60d return) and volatility (60d rolling std of returns)
panel = (
    panel.with_columns(
        ret=pl.col("close").pct_change().over("symbol"),
    )
    .with_columns(
        momentum_60d=pl.col("ret").rolling_mean(60).over("symbol"),
        vol_60d=pl.col("ret").rolling_std(60).over("symbol"),
    )
    .drop_nulls()
)

print(f"Panel: {panel['symbol'].n_unique()} assets, {len(panel):,} obs")

# %%
# Cross-sectional rank, percentile, and z-score at each date
panel = panel.with_columns(
    rank_momentum=pl.col("momentum_60d").rank().over("timestamp"),
    rank_vol=pl.col("vol_60d").rank().over("timestamp"),
    pct_momentum=(pl.col("momentum_60d").rank().over("timestamp") - 1)
    / (pl.col("momentum_60d").count().over("timestamp") - 1),
    zscore_momentum=(
        (pl.col("momentum_60d") - pl.col("momentum_60d").mean().over("timestamp"))
        / pl.col("momentum_60d").std().over("timestamp")
    ),
)

# Show a snapshot for one date
snapshot_date = panel["timestamp"].unique().sort().to_list()[-252]  # ~1 year ago
snapshot = (
    panel.filter(pl.col("timestamp") == snapshot_date)
    .select(["symbol", "momentum_60d", "rank_momentum", "pct_momentum", "zscore_momentum"])
    .sort("rank_momentum")
)
print(f"\nCross-sectional snapshot ({snapshot_date}):")
display(snapshot.to_pandas())

# %%
# Visualize rank evolution over time for selected assets
fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

highlight = ["SPY", "GLD", "XLE"]
for sym in highlight:
    sym_data = panel.filter(pl.col("symbol") == sym).to_pandas()
    axes[0].plot(sym_data["timestamp"], sym_data["pct_momentum"], label=sym, linewidth=0.8)
    axes[1].plot(sym_data["timestamp"], sym_data["rank_vol"], label=sym, linewidth=0.8)

axes[0].set_ylabel("Momentum Percentile")
axes[0].set_title("Cross-Sectional Momentum Percentile Over Time")
axes[0].legend()
axes[0].axhline(0.5, color="gray", linestyle="--", alpha=0.3)

axes[1].set_ylabel("Volatility Rank")
axes[1].set_title("Cross-Sectional Volatility Rank Over Time")
axes[1].legend()

plt.tight_layout()
plt.show()

# %% [markdown]
# Ranks are bounded, stationary, and immune to outliers. A high
# momentum rank (top of universe) followed by rank decay signals
# momentum exhaustion — a feature that absolute momentum values
# cannot express because their scale drifts with market conditions.

# %% [markdown]
# ## Relative Temporal Features
#
# Many temporal features gain signal when expressed relative to a
# benchmark — isolating idiosyncratic dynamics from systematic exposure.
# We demonstrate market-relative features (vs SPY). Sector-relative
# features follow the same pattern with a sector ETF as benchmark
# (e.g., XLE momentum minus XLF momentum isolates energy-specific signal).

# %%
# Compute relative features: asset vs. market (SPY as proxy)
spy_features = (
    panel.filter(pl.col("symbol") == "SPY")
    .select(["timestamp", "momentum_60d", "vol_60d"])
    .rename({"momentum_60d": "mkt_momentum", "vol_60d": "mkt_vol"})
)

panel_rel = panel.join(spy_features, on="timestamp", how="inner")

panel_rel = panel_rel.with_columns(
    momentum_vs_market=pl.col("momentum_60d") - pl.col("mkt_momentum"),
    vol_vs_market=pl.col("vol_60d") / pl.col("mkt_vol"),
)

# Compare absolute vs relative for one asset
xle = panel_rel.filter(pl.col("symbol") == "XLE").to_pandas()

fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

axes[0].plot(xle["timestamp"], xle["momentum_60d"], label="Absolute momentum", linewidth=0.8)
axes[0].plot(
    xle["timestamp"], xle["momentum_vs_market"], label="Momentum vs. market", linewidth=0.8
)
axes[0].axhline(0, color="gray", linestyle="--", alpha=0.3)
axes[0].set_ylabel("60-Day Momentum")
axes[0].set_title("XLE: Absolute vs. Relative Momentum")
axes[0].legend()

# Absolute vol (~0.005 daily) and the vol-ratio (~1-3.5) live on different scales,
# so use twin y-axes to keep both visible.
ax_abs = axes[1]
ax_rel = ax_abs.twinx()
(l_abs,) = ax_abs.plot(
    xle["timestamp"], xle["vol_60d"], label="Absolute volatility", linewidth=0.8, color="#1f77b4"
)
(l_rel,) = ax_rel.plot(
    xle["timestamp"],
    xle["vol_vs_market"],
    label="Vol / Market vol",
    linewidth=0.8,
    color="#d62728",
)
ax_rel.axhline(1.0, color="gray", linestyle="--", alpha=0.3)
ax_abs.set_ylabel("Absolute Volatility (60-day)", color="#1f77b4")
ax_rel.set_ylabel("Vol / Market Vol", color="#d62728")
ax_abs.tick_params(axis="y", labelcolor="#1f77b4")
ax_rel.tick_params(axis="y", labelcolor="#d62728")
ax_abs.set_title("XLE: Absolute vs. Relative Volatility")
ax_abs.legend(handles=[l_abs, l_rel], loc="upper left")

plt.tight_layout()
plt.show()

# %% [markdown]
# Relative momentum removes the market component — XLE's energy-sector
# bets become visible only after subtracting the broad equity trend.
# Relative volatility (vol / market vol) distinguishes periods where an
# asset is genuinely more volatile from periods where the entire market
# is stressed. Values above 1.0 indicate the asset carries more risk
# than the market; values below 1.0 indicate relative calm.

# %% [markdown]
# ## Multi-Asset Regime Aggregation (Illustrative)
#
# Individual regime classifications (from `11_hmm_regimes`) gain
# power when aggregated across the universe. We approximate per-asset
# regime probabilities from realized volatility percentiles to
# demonstrate the aggregation mechanics. This is exploratory — the
# full-history ranking introduces look-ahead. In production, use
# HMM-filtered regime probabilities estimated within walk-forward folds.

# %%
# Illustrative regime probabilities from rolling volatility (proxy for HMM crisis probs).
# WARNING: This ranks each asset's volatility over its FULL history — a look-ahead
# shortcut for demonstrating the aggregation mechanics. In production, use per-asset
# HMM filtered probabilities estimated inside walk-forward folds (see 11_hmm_regimes).
vol_panel = panel.select(["timestamp", "symbol", "vol_60d"]).with_columns(
    crisis_prob=(
        pl.col("vol_60d").rank().over(["symbol"]) / pl.col("vol_60d").count().over(["symbol"])
    ),
)

# Universe-level aggregations at each timestamp
regime_agg = (
    vol_panel.group_by("timestamp")
    .agg(
        universe_crisis_prob=pl.col("crisis_prob").mean(),
        crisis_breadth=(pl.col("crisis_prob") > 0.7).mean(),  # fraction above 70th pct
        regime_dispersion=pl.col("crisis_prob").std(),
    )
    .sort("timestamp")
)

# Plot against SPY returns
spy_ret = (
    etf_data.filter(pl.col("symbol") == "SPY")
    .select(["timestamp", "close"])
    .sort("timestamp")
    .with_columns(spy_cumret=(pl.col("close") / pl.col("close").first() - 1) * 100)
)

agg_pd = regime_agg.join(
    spy_ret.select(["timestamp", "spy_cumret"]), on="timestamp", how="inner"
).to_pandas()

# %%
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

axes[0].plot(agg_pd["timestamp"], agg_pd["spy_cumret"], linewidth=0.8, color="gray")
axes[0].set_ylabel("SPY Cumulative Return (%)")
axes[0].set_title("Market Context")

axes[1].plot(agg_pd["timestamp"], agg_pd["crisis_breadth"], linewidth=0.8, color="#d62728")
axes[1].fill_between(agg_pd["timestamp"], 0, agg_pd["crisis_breadth"], alpha=0.2, color="#d62728")
axes[1].set_ylabel("Crisis Breadth")
axes[1].set_title("Fraction of Assets in High-Volatility Regime")

axes[2].plot(agg_pd["timestamp"], agg_pd["regime_dispersion"], linewidth=0.8, color="#1f77b4")
axes[2].set_ylabel("Regime Dispersion")
axes[2].set_title("Cross-Sectional Dispersion of Regime Probabilities")

plt.tight_layout()
plt.show()

# %% [markdown]
# Crisis breadth spikes during market selloffs (visible as SPY drawdowns),
# confirming that aggregating individual regime indicators produces
# a useful system-wide stress signal. Regime dispersion is low during
# broad crises (all assets move together) and high during sector
# rotations (differentiated behavior) — making it a complement to
# crisis breadth rather than a substitute.

# %% [markdown]
# ## Key Takeaways
#
# 1. **Cross-sectional ranking normalizes temporal features** — ranks are
#    scale-invariant, outlier-robust, and stationary, making heterogeneous
#    assets comparable for panel models
# 2. **Relative features decompose alpha from beta** — asset momentum minus
#    sector momentum isolates idiosyncratic signal from systematic exposure
# 3. **Two cointegration tests for robustness** — Engle-Granger and Johansen
#    may disagree; treat pairs with caution when they do
# 4. **Kalman filter adapts** — dynamic hedge ratios track structural shifts
#    that a static OLS estimate misses
# 5. **Half-life guides lookback** — use 2x half-life for rolling z-score
#    windows; short half-lives are more practical for trading
# 6. **Screen multiple pairs** — economic relatedness does not guarantee
#    cointegration; the `ml4t-engineer` library scales this analysis
#    with Polars-native expressions
#
# **Previous**: `13_regime_as_feature` for regime features in ML pipelines.
# **Book**: Chapter 9, Section 9.6 discusses cross-sectional ranking,
# relative features, and panel aggregation in depth.