בניית מאפייני קארי, בין־נכסיים ומאופציות

# 03_structural_cross_instrument_features.py


```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]
# # Structural and Cross-Instrument Features
#
# **Chapter 8: Feature Engineering**
# **Section Reference**: 8.3 - Structural and Cross-Instrument Features
# **Docker image**: `ml4t`
#
# ## Purpose
#
# This notebook demonstrates features that require **data beyond a single asset's
# price series**: term structures, cross-instrument relationships, and
# derivatives-implied quantities. These encode information invisible in any
# individual price history.
#
# ## Learning Objectives
#
# 1. Compute carry (roll yield) from futures term structure data
# 2. Build cross-asset beta and lead-lag correlation features
# 3. Extract options-implied features: ATM IV, risk reversal, IV term structure, VRP
# 4. Understand surface stability requirements for options features
#
# ## Feature Families Covered
#
# | Family | Representative Features | Data Source |
# |--------|------------------------|-------------|
# | **Carry & Term Structure** | Roll yield, curve slope/curvature | CME futures |
# | **Cross-Asset Structure** | Rolling beta, lead-lag, relative value | ETFs |
# | **Options-Implied** | ATM IV, risk reversal, VRP, term slope | S&P 500 options |
#
# ## Data Policy
#
# All examples use **real data**: CME futures, ETFs, and S&P 500 options.

# %%
"""Structural and Cross-Instrument Features: carry, cross-asset, and options-implied families."""

from __future__ import annotations

from datetime import datetime
from typing import cast

import numpy as np
import plotly.graph_objects as go
import polars as pl
from plotly.subplots import make_subplots

from utils.paths import get_chapter_dir

# Importing utils.style registers and activates the ML4T Plotly template
# (house palette, fonts, gridlines) as the repo-wide default, so every Plotly
# figure below inherits the book style; matplotlib is styled via matplotlibrc.
from utils.style import COLORS, show_plotly_with_alt, show_with_alt

# %% tags=["parameters"]
SEED = 42
FUTURES_START_DATE = ""
ETF_START_DATE = ""
OPT_SYMBOLS = ["AAPL", "MSFT", "GOOGL", "AMZN", "JPM", "BA", "XOM", "KO"]
OPTIONS_START_DATE = "2019-01-01"


# %%
def ensure_df(df: pl.DataFrame | pl.LazyFrame) -> pl.DataFrame:
    """Normalize loader outputs to eager DataFrames for notebook use."""
    return cast(pl.DataFrame, df.collect() if isinstance(df, pl.LazyFrame) else df)


# %% [markdown]
# ---
#
# # Part 1: Carry and Term Structure
#
# Carry features encode **risk compensation**: the return from holding a position
# absent price change. Roll yield from futures is the canonical example.
#
# Key formula:
#
# $$\text{carry}_{t,c} = \frac{F_{t,c}^{\text{near}} - F_{t,c}^{\text{far}}}{F_{t,c}^{\text{near}}} \times \frac{365}{\Delta T}$$
#
# - Positive carry (backwardation): near > far, so a holder earns by rolling
# - Negative carry (contango): near < far, so a holder pays to roll

# %% [markdown]
# ## Load Futures Data

# %%
from data import load_cme_futures

# Focus on 4 liquid products spanning asset classes
PRODUCTS = ["ES", "CL", "GC", "ZN"]

futures = ensure_df(load_cme_futures(products=PRODUCTS, tenors=[0, 1, 2]))

if FUTURES_START_DATE:
    futures = futures.filter(
        pl.col("session_date") >= datetime.strptime(FUTURES_START_DATE, "%Y-%m-%d")
    )

print(f"Futures: {len(futures):,} rows")
print(f"Products: {futures['product'].unique().to_list()}")
print(f"Tenors: {futures['tenor'].unique().sort().to_list()}")
print(f"Date range: {futures['session_date'].min()} to {futures['session_date'].max()}")

# %% [markdown]
# ## Roll Yield (Manual Implementation)
#
# Roll yield is the annualized price difference between near and far contracts.
# We pivot tenors to compute the spread.


# %%
def compute_roll_yield(df: pl.DataFrame) -> pl.DataFrame:
    """
    Compute annualized roll yield from tenor 0 (front) and tenor 1 (deferred).

    Roll yield = (F_near - F_far) / F_near * (365 / DeltaT)

    We approximate DeltaT ≈ 30 days (typical monthly roll).

    Reads ``raw_close`` rather than the adjusted series. Term structure is a spread
    between contemporaneous tenor *levels*, and the ratio-adjusted series bakes each
    tenor's accumulated roll history into its level, which is not the curve.
    """
    front = df.filter(pl.col("tenor") == 0).select(
        ["session_date", "product", pl.col("raw_close").alias("close_front")]
    )
    deferred = df.filter(pl.col("tenor") == 1).select(
        ["session_date", "product", pl.col("raw_close").alias("close_deferred")]
    )

    # Join and compute carry
    carry = front.join(deferred, on=["session_date", "product"], how="inner")
    carry = carry.with_columns(
        [
            # Raw roll yield
            ((pl.col("close_front") - pl.col("close_deferred")) / pl.col("close_front")).alias(
                "roll_yield_raw"
            ),
            # Annualized (assuming ~30 day roll period)
            (
                (pl.col("close_front") - pl.col("close_deferred"))
                / pl.col("close_front")
                * (365 / 30)
            ).alias("roll_yield_ann"),
        ]
    )

    return carry.sort(["product", "session_date"])


carry_df = compute_roll_yield(futures)

print("Roll yield (last rows per product):")
for product in PRODUCTS:
    row = carry_df.filter(pl.col("product") == product).tail(1)
    if len(row) > 0:
        ry = row["roll_yield_ann"][0]
        print(f"  {product}: {ry:+.4f} ({'backwardation' if ry > 0 else 'contango'})")

# %% [markdown]
# ## Term Structure Slope and Curvature
#
# With three tenors we can also extract slope and curvature:
#
# - **Slope**: $(F_0 - F_2) / F_0$, the overall term structure direction
# - **Curvature**: $F_0 - 2 \cdot F_1 + F_2$, the butterfly shape


# %%
def compute_term_structure_features(df: pl.DataFrame) -> pl.DataFrame:
    """Compute slope and curvature from 3-tenor term structure."""
    # Contemporaneous tenor levels → raw_close (see compute_roll_yield note).
    t0 = df.filter(pl.col("tenor") == 0).select(
        ["session_date", "product", pl.col("raw_close").alias("c0")]
    )
    t1 = df.filter(pl.col("tenor") == 1).select(
        ["session_date", "product", pl.col("raw_close").alias("c1")]
    )
    t2 = df.filter(pl.col("tenor") == 2).select(
        ["session_date", "product", pl.col("raw_close").alias("c2")]
    )

    ts = t0.join(t1, on=["session_date", "product"], how="inner").join(
        t2, on=["session_date", "product"], how="inner"
    )

    return ts.with_columns(
        [
            # Slope: normalized front-back spread
            ((pl.col("c0") - pl.col("c2")) / pl.col("c0")).alias("ts_slope"),
            # Curvature: butterfly (positive = humped)
            (pl.col("c0") - 2 * pl.col("c1") + pl.col("c2")).alias("ts_curvature"),
        ]
    )


ts_df = compute_term_structure_features(futures)

print("\nTerm structure features:")
ts_df.select(["session_date", "product", "ts_slope", "ts_curvature"]).tail(8)

# %% [markdown]
# ## Visualize Carry Across Products

# %%
# Smooth carry with 21-day EMA for visualization
carry_smooth = carry_df.sort(["product", "session_date"]).with_columns(
    pl.col("roll_yield_ann")
    .ewm_mean(span=21, ignore_nulls=True)
    .over("product")
    .alias("carry_ema21")
)

fig = make_subplots(
    rows=2,
    cols=2,
    subplot_titles=[f"{p}: annualized roll yield" for p in PRODUCTS],
    vertical_spacing=0.12,
    horizontal_spacing=0.08,
)

for idx, product in enumerate(PRODUCTS):
    row, col = idx // 2 + 1, idx % 2 + 1
    subset = carry_smooth.filter(pl.col("product") == product)
    fig.add_trace(
        go.Scatter(
            x=subset["session_date"].to_list(),
            y=subset["carry_ema21"].to_list(),
            name=product,
            fill="tozeroy",
        ),
        row=row,
        col=col,
    )
    fig.add_hline(y=0, line_dash="dash", line_color=COLORS["neutral"], row=row, col=col)

fig.update_layout(height=500, title="Carry (Roll Yield) Across Products", showlegend=False)
fig.update_yaxes(title_text="Roll Yield (ann.)", row=1, col=1)
fig.update_yaxes(title_text="Roll Yield (ann.)", row=2, col=1)
show_plotly_with_alt(
    fig,
    (
        "A two-by-two grid of annualized roll yield, one panel per product, over roughly "
        "fifteen years to the mid-2020s, each drawn against a dashed zero line. Each panel "
        "has its own vertical scale and they differ greatly, so only the position relative "
        "to zero is comparable between panels. ES oscillates tightly above zero for the "
        "first stretch, sits near zero in the late 2010s, spikes upward around 2020 and "
        "then settles into a sustained band below zero for the final years. CL is close to "
        "zero for most of the span apart from dips below it in the mid-2010s and one "
        "isolated plunge in 2020 that reaches far beyond anything else in the grid and "
        "sets that panel's scale on its own. GC hugs zero early, then swings increasingly "
        "far below it from the early 2020s with tall recoveries above zero in between. ZN "
        "crosses zero rapidly and symmetrically across the whole span, with its widest "
        "swings in the mid-2010s and a visibly calmer stretch after 2020."
    ),
)

# %% [markdown]
# The four panels carry very different magnitudes, so each has its own vertical scale and
# only the sign is comparable across them. Count the signs rather than reading them off
# the panel heights.

# %%
backwardation = (
    carry_df.drop_nulls("roll_yield_ann")
    .group_by("product")
    .agg(
        [
            (pl.col("roll_yield_ann") > 0).mean().alias("share_backwardated"),
            pl.col("roll_yield_ann").median().alias("median_roll_yield"),
            pl.len().alias("n_sessions"),
        ]
    )
    .sort("share_backwardated", descending=True)
)
print("share of sessions in backwardation (positive roll yield), by product:")
print(backwardation)

# %% [markdown]
# **Interpretation**: read the shares printed above rather than the panel heights, which
# are on four different scales. Backwardation means the front contract trades above the
# deferred one, so a holder who rolls earns the difference; contango is the reverse and
# the holder pays. The table says how often each product sat on each side of that line
# over this sample.
#
# The received explanations are worth holding next to the numbers rather than in place of
# them. Backwardation in an energy contract like CL is usually read as physical tightness,
# a market paying up for a barrel today. Gold is the standard contango example, because
# storage and financing dominate and there is no convenience yield to offset them. A
# note like ZN carries the coupon against a financing rate, so its sign follows that
# spread. ES sits on the cost-of-carry relation between dividends and financing. Whether
# this sample agrees with each of those is what the table is for, and a sample this short
# settles none of them.
#
# **Crypto funding rates** operate on a distinct clock (8-hour settlements) with
# much higher volatility. See the `crypto_perps_funding` case study for the
# full treatment.
#
# **Full pipeline**: The `cme_futures` case study applies these features to all
# 30 CME products with proper roll handling and cross-sectional ranks.

# %% [markdown]
# ### Book Figure: Cross-Sectional Roll Yield Snapshot
#
# Load all CME products and plot a single-date snapshot of annualized roll yield
# sorted by magnitude, which is the canonical carry signal of this section.

# %%
import matplotlib.pyplot as plt

# Load all products for the cross-section
all_futures = ensure_df(load_cme_futures(tenors=[0, 1]))

# Compute roll yield for all products
all_carry = compute_roll_yield(all_futures)

# Pick a representative date with good coverage (latest date with most products)
date_counts = all_carry.group_by("session_date").len().sort("len", descending=True)
len_max = date_counts["len"].max() if len(date_counts) > 0 else 0
max_products = len_max if isinstance(len_max, int) else 0
min_products = min(15, max_products) if max_products > 0 else 1
snapshot_date = date_counts.filter(pl.col("len") >= min_products)["session_date"].max()

snapshot = (
    all_carry.filter(pl.col("session_date") == snapshot_date)
    .select(["product", "roll_yield_ann"])
    .drop_nulls()
    .sort("roll_yield_ann")
)

print(f"Snapshot date: {snapshot_date}")
print(f"Products: {len(snapshot)}")

# %%
# Sector mapping for fill patterns
SECTORS = {
    "CL": "Energy",
    "HO": "Energy",
    "NG": "Energy",
    "RB": "Energy",
    "GC": "Gold/Silver",
    "SI": "Gold/Silver",
    "HG": "Metals",
    "PA": "Metals",
    "PL": "Metals",
    "ZC": "Agriculture",
    "ZS": "Agriculture",
    "ZW": "Agriculture",
    "ZM": "Agriculture",
    "ZL": "Agriculture",
    "KC": "Agriculture",
    "CT": "Agriculture",
    "SB": "Agriculture",
    "CC": "Agriculture",
    "LC": "Agriculture",
    "LH": "Agriculture",
    "FC": "Agriculture",
    "ES": "Financials",
    "NQ": "Financials",
    "YM": "Financials",
    "RTY": "Financials",
    "ZN": "Rates",
    "ZB": "Rates",
    "ZF": "Rates",
    "ZT": "Rates",
    "GE": "Rates",
    "6E": "FX",
    "6J": "FX",
    "6B": "FX",
    "6A": "FX",
    "6C": "FX",
    "6S": "FX",
}
SECTOR_GRAYS = {
    "Energy": "0.15",
    "Gold/Silver": "0.30",
    "Metals": "0.40",
    "Agriculture": "0.55",
    "Financials": "0.70",
    "Rates": "0.80",
    "FX": "0.90",
}

# %% [markdown]
# The bars, the regime annotations and the sector legend are built in one cell because the
# inline backend flushes a figure at the end of the cell that created it; splitting the
# legend into a later cell would display the bars before it is attached.

# %%
products = snapshot["product"].to_list()
yields = snapshot["roll_yield_ann"].to_list()
colors = [SECTOR_GRAYS.get(SECTORS.get(p, "Other"), "0.50") for p in products]

fig_mpl, ax = plt.subplots(figsize=(12, 5))
ax.bar(range(len(products)), yields, color=colors, edgecolor="black", linewidth=0.5)
ax.set_xticks(range(len(products)))
ax.set_xticklabels(products, fontsize=7, rotation=45, ha="right")
ax.axhline(y=0, color="black", linewidth=0.8)
ax.set_ylabel("Annualized Roll Yield")
ax.set_title(f"Annualized Roll Yield Across CME Products ({snapshot_date})")

# Regime annotations
y_max, y_min = (max(yields), min(yields)) if yields else (0.1, -0.1)
text_kw = dict(ha="center", va="center", fontsize=8, style="italic", color="0.3")
ax.text(len(products) * 0.75, y_max * 0.7, "Backwardation\n(holder earns)", **text_kw)
ax.text(len(products) * 0.25, y_min * 0.7, "Contango\n(holder pays)", **text_kw)

# Compact sector legend for the grayscale fills
unique_sectors = sorted(set(SECTORS.get(p, "Other") for p in products))
legend_handles = [
    plt.Rectangle((0, 0), 1, 1, fc=SECTOR_GRAYS.get(s, "0.50"), ec="black", lw=0.5)
    for s in unique_sectors
]
ax.legend(legend_handles, unique_sectors, loc="upper left", fontsize=7, frameon=False, ncol=4)

show_with_alt(
    fig_mpl,
    (
        "A bar chart of annualized roll yield across the CME product set on a single "
        "snapshot date, sorted from most negative on the left to most positive on the "
        "right, with a solid line at zero. Bars are filled in shades of grey keyed to a "
        "sector legend across the top. Most products cluster in a shallow band just "
        "either side of zero. The left tail is led by a single agriculture contract "
        "reaching well below the rest, followed by two more agriculture bars. The right "
        "tail is dominated by one energy contract whose bar towers over everything else "
        "in the chart, several times the height of the next largest. Italic annotations "
        "label the region above zero as backwardation, where the holder earns, and the "
        "region below as contango, where the holder pays."
    ),
)

# %%
# Persisted so the book figure script can re-render at print resolution.
_FIG_8_4_ARTIFACT = (
    get_chapter_dir(8) / "output" / "book_figure_artifacts" / "figure_8_4_carry_roll_yield.parquet"
)
_FIG_8_4_ARTIFACT.parent.mkdir(parents=True, exist_ok=True)
snapshot.with_columns(pl.lit(snapshot_date).alias("snapshot_date")).write_parquet(_FIG_8_4_ARTIFACT)

# %% [markdown]
# ---
#
# # Part 2: Cross-Asset Structure and Relative Value
#
# Cross-asset features model relationships *between* instruments:
# - **Rolling beta**: Sensitivity to a common factor (market, sector)
# - **Lead-lag**: Temporal structure across markets
# - **Relative value**: Deviation from peer mean (z-scored)

# %% [markdown]
# ## Load ETF Data

# %%
from data import load_etfs

etfs = load_etfs().sort(["symbol", "timestamp"])

if ETF_START_DATE:
    etfs = etfs.filter(pl.col("timestamp") >= datetime.strptime(ETF_START_DATE, "%Y-%m-%d"))

# Use a subset for cross-asset analysis: market + sectors + bonds + gold
CA_SYMBOLS = ["SPY", "QQQ", "XLF", "XLE", "XLK", "XLV", "TLT", "GLD"]
ca_etfs = etfs.filter(pl.col("symbol").is_in(CA_SYMBOLS)).sort(["symbol", "timestamp"])

print(f"Cross-asset universe: {ca_etfs['symbol'].n_unique()} assets, {len(ca_etfs):,} rows")

# %% [markdown]
# ## Rolling Beta to Market
#
# Rolling beta measures time-varying sensitivity to the market factor:
#
# $$\beta_{t,a} = \frac{\text{Cov}(r_a, r_m)_t}{\text{Var}(r_m)_t}$$

# %%
from ml4t.engineer.features.cross_asset import beta_to_market

# Compute returns
ca_returns = ca_etfs.with_columns(pl.col("close").pct_change().over("symbol").alias("ret"))

# Get SPY returns as market factor
spy_ret = ca_returns.filter(pl.col("symbol") == "SPY").select(
    pl.col("timestamp"), pl.col("ret").alias("market_ret")
)

# Join market returns to all ETFs
ca_with_market = ca_returns.join(spy_ret, on="timestamp", how="inner")

# Rolling 63-day beta
beta_df = ca_with_market.with_columns(
    beta_to_market("ret", "market_ret", window=63).alias("beta_63d")
)

# Show latest cross-section
latest = beta_df["timestamp"].max()
print(f"Rolling 63-day Beta to SPY ({latest}):")
(
    beta_df.filter(pl.col("timestamp") == latest)
    .select(["symbol", "beta_63d"])
    .sort("beta_63d", descending=True)
)

# %% [markdown]
# ## Beta-Adjusted Residual Momentum
#
# Isolate stock-specific momentum by subtracting market contribution:
#
# $$r^{\text{resid}}_{t,a} = r_{t,a} - \beta_{t,a} \cdot r_{t,m}$$

# %%
# Compute residual returns and residual momentum
resid_mom_df = (
    beta_df.with_columns(
        (pl.col("ret") - pl.col("beta_63d") * pl.col("market_ret")).alias("resid_ret")
    )
    .sort(["symbol", "timestamp"])
    .with_columns(
        [
            pl.col("ret").rolling_sum(21).over("symbol").alias("raw_mom_21d"),
            pl.col("resid_ret").rolling_sum(21).over("symbol").alias("resid_mom_21d"),
        ]
    )
)

# Compare raw vs residual momentum for the latest date
print(f"\nRaw vs Residual 21d Momentum ({latest}):")
(
    resid_mom_df.filter(pl.col("timestamp") == latest)
    .select(["symbol", "raw_mom_21d", "resid_mom_21d", "beta_63d"])
    .sort("resid_mom_21d", descending=True)
)

# %% [markdown]
# **Interpretation**: High-beta ETFs (XLK, QQQ) often show strong raw momentum
# that is largely market-driven. Residual momentum isolates the ETF-specific
# component, which is what a rotation strategy wants.

# %% [markdown]
# ## Lead-Lag Correlations

# %% [markdown]
# Does SPY lead the sector ETFs? The loop below correlates the SPY return at t with each
# sector's return at t plus a lag, so a non-zero correlation at a positive lag is SPY
# moving first.

# %%
sectors = ["XLF", "XLE", "XLK", "XLV"]

for sector in sectors:
    sector_ret = ca_returns.filter(pl.col("symbol") == sector).select(
        pl.col("timestamp"), pl.col("ret").alias("sector_ret")
    )
    paired = spy_ret.join(sector_ret, on="timestamp", how="inner").drop_nulls()

    if len(paired) > 252:
        # Lagged correlation: SPY_t vs Sector_{t+1}
        corr_0 = paired.select(pl.corr("market_ret", "sector_ret")).item()
        # Lag 1: shift sector returns
        paired_lag = paired.with_columns(pl.col("sector_ret").shift(-1).alias("sector_ret_lag1"))
        corr_1 = paired_lag.drop_nulls().select(pl.corr("market_ret", "sector_ret_lag1")).item()
        print(f"  SPY → {sector}: contemporaneous={corr_0:.3f}, lag-1={corr_1:.3f}")

# %% [markdown]
# **Interpretation**: If lag-1 correlation (SPY_t vs Sector_{t+1}) is meaningfully
# positive, SPY may lead the sector -- a tradeable signal if the lead is genuine
# and not a staleness artifact (e.g., less liquid sector ETFs closing before SPY).
# Always verify that the leading instrument actually traded at the "lead" timestamp.

# %% [markdown]
# ## Deviation from Peer Mean (Relative Value Z-Score)
#
# For each ETF, compute how far its return deviates from the cross-sectional mean:
#
# $$z_{t,a} = \frac{r_{t,a} - \bar{r}_t}{\sigma_t}$$

# %%
# Cross-sectional z-score of 21-day momentum
cs_mom = ca_returns.sort(["symbol", "timestamp"]).with_columns(
    pl.col("ret").rolling_sum(21).over("symbol").alias("mom_21d")
)

cs_zscore = cs_mom.with_columns(
    (
        (pl.col("mom_21d") - pl.col("mom_21d").mean().over("timestamp"))
        / pl.col("mom_21d").std().over("timestamp")
    ).alias("mom_21d_zscore")
)

print(f"\nRelative value z-scores ({latest}):")
(
    cs_zscore.filter(pl.col("timestamp") == latest)
    .select(["symbol", "mom_21d", "mom_21d_zscore"])
    .sort("mom_21d_zscore", descending=True)
)

# %% [markdown]
# **Interpretation**: Relative-value z-scores identify assets that have diverged
# from their peer group. A z-score of +2 means the asset's 21-day momentum is
# two standard deviations above the cross-sectional mean -- a candidate for
# mean-reversion or confirmation of breakout depending on the strategy.

# %% [markdown]
# **Full pipelines**: The `etfs` case study (`case_studies/etfs/03_financial_features`)
# applies cross-asset features to the full 100-ETF universe. The `fx_pairs` case
# study uses analogous relative value features for currency pairs.

# %% [markdown]
# ---
#
# # Part 3: Options-Implied Features
#
# Options encode the market's probabilistic assessment of future outcomes.
# Four key features capture distinct information:
#
# | Feature | What It Measures | Role |
# |---------|-----------------|------|
# | ATM IV | Expected future volatility | State |
# | Risk reversal (25$\delta$ skew) | Directional fear/greed | State |
# | IV term structure slope | Near-term vs long-term stress | State |
# | VRP (IV minus RV) | Variance risk premium | Signal/state |

# %% [markdown]
# ## Load Options and Equity Data

# %%
from data import load_sp500_daily_bars, load_sp500_options_eda

options: pl.DataFrame = ensure_df(
    load_sp500_options_eda(
        symbols=OPT_SYMBOLS,
        start_date=OPTIONS_START_DATE,
        include_greeks=True,
    )
)
# Quality filter: converged IV only
options = options.filter(pl.col("iv_convergence") == "Converged")

equities: pl.DataFrame = ensure_df(
    load_sp500_daily_bars(
        symbols=OPT_SYMBOLS,
        start_date=OPTIONS_START_DATE,
    )
)

print(f"Options: {len(options):,} rows, {options['symbol'].n_unique()} assets")
print(f"Date range: {options['timestamp'].min()} to {options['timestamp'].max()}")
print(f"Equities: {len(equities):,} rows")

# %% [markdown]
# ## ATM Implied Volatility
#
# ATM IV is the market's expectation of future volatility over the option's life. The
# selection below is on **moneyness**, not on delta: `compute_atm_iv` keeps calls whose
# strike sits inside `moneyness_range` of the underlying and whose maturity sits inside
# `dte_range`, then takes the strike closest to the money on each day. Near the money and
# near a one-month maturity the closest-to-the-money call has a delta in the region of a
# half, which is where "the fifty-delta" as a name for this contract comes from, but the
# delta is never read and a contract is never selected by it. The defaults are in the
# function signature below.


# %%
def compute_atm_iv(
    df: pl.DataFrame,
    dte_range: tuple[int, int] = (25, 35),
    moneyness_range: tuple[float, float] = (0.98, 1.02),
) -> pl.DataFrame:
    """
    Compute ATM implied volatility per symbol per day.

    Selects calls near ATM with ~30 days to expiry.
    """
    return (
        df.with_columns((pl.col("strike") / pl.col("underlying_price")).alias("moneyness"))
        .filter(pl.col("moneyness").is_between(*moneyness_range))
        .filter(pl.col("days_to_maturity").is_between(*dte_range))
        .filter(pl.col("call_put") == "C")
        # Pick strike closest to ATM per day/asset
        .with_columns((pl.col("moneyness") - 1.0).abs().alias("atm_dist"))
        .sort(["timestamp", "symbol", "atm_dist"])
        .group_by(["timestamp", "symbol"])
        .first()
        .select(["timestamp", "symbol", "underlying_price", "implied_vol"])
        .rename({"implied_vol": "iv_atm_30"})
        .sort(["symbol", "timestamp"])
    )


atm_iv = compute_atm_iv(options)
print(f"ATM IV: {len(atm_iv):,} rows")
atm_iv.head(5)

# %% [markdown]
# ## Risk Reversal (25-Delta Skew)
#
# The risk reversal measures directional skew: the price of downside protection
# relative to upside:
#
# $$\mathrm{RR}_{25\delta} = IV_{25\delta,\,\mathrm{put}} - IV_{25\delta,\,\mathrm{call}}$$
#
# Positive values indicate puts are more expensive (crash fear).


# %%
def compute_risk_reversal(
    df: pl.DataFrame,
    dte_range: tuple[int, int] = (25, 35),
    otm_put_range: tuple[float, float] = (0.93, 0.97),
    otm_call_range: tuple[float, float] = (1.03, 1.07),
) -> pl.DataFrame:
    """
    Compute 25-delta risk reversal (put IV - call IV).

    Positive = puts more expensive (bearish skew / crash fear).
    """
    df = df.with_columns((pl.col("strike") / pl.col("underlying_price")).alias("moneyness")).filter(
        pl.col("days_to_maturity").is_between(*dte_range)
    )

    # OTM puts
    otm_puts = (
        df.filter((pl.col("call_put") == "P") & pl.col("moneyness").is_between(*otm_put_range))
        .group_by(["timestamp", "symbol"])
        .agg(pl.col("implied_vol").mean().alias("iv_put_25d"))
    )

    # OTM calls
    otm_calls = (
        df.filter((pl.col("call_put") == "C") & pl.col("moneyness").is_between(*otm_call_range))
        .group_by(["timestamp", "symbol"])
        .agg(pl.col("implied_vol").mean().alias("iv_call_25d"))
    )

    return (
        otm_puts.join(otm_calls, on=["timestamp", "symbol"], how="inner")
        .with_columns((pl.col("iv_put_25d") - pl.col("iv_call_25d")).alias("risk_reversal_25d"))
        .sort(["symbol", "timestamp"])
    )


rr_df = compute_risk_reversal(options)
print(f"Risk reversal: {len(rr_df):,} rows")
rr_df.head(5)

# %% [markdown]
# ## IV Term Structure Slope
#
# The ratio of short-dated to long-dated IV captures near-term stress:
#
# $$\text{Term Slope} = \frac{IV_{\text{short}}}{IV_{\text{long}}}$$
#
# - $> 1$: Inverted (near-term stress, fear of imminent event)
# - $< 1$: Normal contango (calm)


# %%
def compute_iv_term_slope(
    df: pl.DataFrame,
    short_dte: tuple[int, int] = (8, 30),
    long_dte: tuple[int, int] = (60, 180),
    moneyness_range: tuple[float, float] = (0.98, 1.02),
) -> pl.DataFrame:
    """
    Compute IV term structure slope = short-dated ATM IV / long-dated ATM IV.
    """
    atm = (
        df.with_columns((pl.col("strike") / pl.col("underlying_price")).alias("moneyness"))
        .filter(pl.col("moneyness").is_between(*moneyness_range))
        .filter(pl.col("call_put") == "C")
    )

    short_iv = (
        atm.filter(pl.col("days_to_maturity").is_between(*short_dte))
        .group_by(["timestamp", "symbol"])
        .agg(pl.col("implied_vol").mean().alias("iv_short"))
    )

    long_iv = (
        atm.filter(pl.col("days_to_maturity").is_between(*long_dte))
        .group_by(["timestamp", "symbol"])
        .agg(pl.col("implied_vol").mean().alias("iv_long"))
    )

    EPSILON = 1e-10
    return (
        short_iv.join(long_iv, on=["timestamp", "symbol"], how="inner")
        .with_columns(
            (pl.col("iv_short") / pl.col("iv_long").clip(EPSILON, None)).alias("iv_term_slope")
        )
        .sort(["symbol", "timestamp"])
    )


term_df = compute_iv_term_slope(options)
print(f"IV term slope: {len(term_df):,} rows")
term_df.head(5)

# %% [markdown]
# ## Variance Risk Premium (IV - RV)
#
# The VRP as computed here is the difference between a forward-looking implied volatility
# and a trailing realized volatility:
#
# $$\text{VRP}_t = IV_{30,\text{atm},t} - RV_{20,t}$$
#
# where $RV_{20,t}$ looks back over the twenty sessions ending at $t$. The realized leg is
# built from split-adjusted closes; the note in `compute_vrp` says why that matters.
#
# Be precise about what this quantity is before reading anything off it. It subtracts a
# **trailing** twenty-day realized volatility from a **forward-looking** thirty-day
# implied volatility, so the two terms describe different and barely overlapping windows.
# It is a spread between two volatility measurements taken on the same day, and it is not
# a return: nothing here buys or sells an option, and no position is held to expiry. The
# quantity a volatility seller actually earns compares the implied volatility quoted at t
# with the realized volatility over the *following* thirty days, which this feature does
# not compute.
#
# With that established, the cell below describes the spread's distribution. Read the
# share of days on which it is positive together with the median and the mean. A spread
# that is positive on most days while its mean sits below its median is left-skewed: small
# positive values most of the time, and occasional negative values large enough to move
# the average on their own. The most positive and most negative values printed say how
# lopsided, and the bottom panel of the figure further down shows the same thing for one
# name: a series near zero whose vertical extent is set by a single deep excursion.
#
# The skew is a property of the spread worth carrying forward, because it is what makes a
# summary statistic misleading here: a mean and a median disagree about the sign, and
# quoting either alone describes a different series. It also suggests treating a feature
# built from this spread as a state variable rather than reading it as a signal at face
# value. Whether selling volatility is profitable is a separate question that needs
# the horizon-matched calculation above, not this table. Note too that the share of
# positive days is not uniform across names: check whether every symbol is above a half
# before treating "the spread is positive" as something to lean on.


# %%
def compute_vrp(
    iv_df: pl.DataFrame,
    equity_df: pl.DataFrame,
    rv_window: int = 20,
) -> pl.DataFrame:
    """
    Compute variance risk premium = ATM IV - realized vol.

    iv_df must have columns: timestamp, symbol, iv_atm_30
    equity_df must have columns: timestamp, symbol, close, adj_factor

    Returns are taken from ``close * adj_factor``, not from ``close``. The bundled
    daily bars carry the as-traded price and the cumulative split factor in separate
    columns, so a raw close steps down by the split ratio on the ex-date: AAPL goes
    from 499.23 to 129.04 on 2020-08-31, which enters a return series as a 74% fall
    and then inflates every rolling window that contains it. Multiplying by
    ``adj_factor`` first restores the continuous series, and that day becomes the
    3% rise it was.
    """
    if "adj_factor" not in equity_df.columns:
        raise KeyError(
            "equity_df needs adj_factor to build a split-continuous return series; "
            f"got {sorted(equity_df.columns)}"
        )

    rv = (
        equity_df.sort(["symbol", "timestamp"])
        .with_columns((pl.col("close") * pl.col("adj_factor")).alias("adj_close"))
        .with_columns(
            (
                pl.col("adj_close").pct_change().over("symbol").rolling_std(rv_window)
                * np.sqrt(252)
            ).alias("rv_20")
        )
        .select(["timestamp", "symbol", "rv_20"])
    )

    return (
        iv_df.join(rv, on=["timestamp", "symbol"], how="inner")
        .with_columns((pl.col("iv_atm_30") - pl.col("rv_20")).alias("vrp"))
        .sort(["symbol", "timestamp"])
    )


vrp_df = compute_vrp(atm_iv, equities)
print(f"VRP: {len(vrp_df):,} rows")
print(vrp_df.select(["timestamp", "symbol", "iv_atm_30", "rv_20", "vrp"]).tail(5))

# The received claim is that this premium is usually positive. Check it on this sample
# rather than repeating it.
print()
print("VRP sign, by symbol:")
print(
    vrp_df.drop_nulls("vrp")
    .group_by("symbol")
    .agg(
        [
            (pl.col("vrp") > 0).mean().alias("share_positive"),
            pl.col("vrp").median().alias("median_vrp"),
            pl.col("iv_atm_30").median().alias("median_iv"),
            pl.col("rv_20").median().alias("median_rv"),
            pl.len().alias("n"),
        ]
    )
    .sort("share_positive", descending=True)
)
_v = vrp_df.drop_nulls("vrp")["vrp"]
print(f"pooled share positive: {(_v > 0).mean():.3f}")
print(f"pooled median {_v.median():+.4f} against a mean of {_v.mean():+.4f}")
print(f"most positive {_v.max():+.3f}, most negative {_v.min():+.3f}")

# %% [markdown]
# ## Visualize Options Features

# %%
# Pick one symbol for visualization
viz_symbol = "AAPL" if "AAPL" in OPT_SYMBOLS else OPT_SYMBOLS[0]

# Combine features for the selected symbol
viz_iv = atm_iv.filter(pl.col("symbol") == viz_symbol)
viz_rr = rr_df.filter(pl.col("symbol") == viz_symbol).select(["timestamp", "risk_reversal_25d"])
viz_term = term_df.filter(pl.col("symbol") == viz_symbol).select(["timestamp", "iv_term_slope"])
viz_vrp = vrp_df.filter(pl.col("symbol") == viz_symbol).select(["timestamp", "vrp"])

viz = (
    viz_iv.join(viz_rr, on="timestamp", how="left")
    .join(viz_term, on="timestamp", how="left")
    .join(viz_vrp, on="timestamp", how="left")
)

# %%
fig = make_subplots(
    rows=4,
    cols=1,
    shared_xaxes=True,
    subplot_titles=[
        f"{viz_symbol}: ATM IV (30d)",
        "Risk Reversal (25δ)",
        "IV Term Slope",
        "Variance Risk Premium (IV − RV)",
    ],
    vertical_spacing=0.06,
)

fig.add_trace(
    go.Scatter(x=viz["timestamp"].to_list(), y=viz["iv_atm_30"].to_list(), name="ATM IV"),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=viz["timestamp"].to_list(),
        y=viz["risk_reversal_25d"].to_list(),
        name="RR 25δ",
        fill="tozeroy",
    ),
    row=2,
    col=1,
)
fig.add_hline(y=0, line_dash="dash", line_color=COLORS["neutral"], row=2, col=1)

fig.add_trace(
    go.Scatter(
        x=viz["timestamp"].to_list(),
        y=viz["iv_term_slope"].to_list(),
        name="Term Slope",
    ),
    row=3,
    col=1,
)
fig.add_hline(y=1, line_dash="dash", line_color=COLORS["neutral"], row=3, col=1)
fig.add_trace(
    go.Scatter(
        x=viz["timestamp"].to_list(),
        y=viz["vrp"].to_list(),
        name="VRP",
        fill="tozeroy",
    ),
    row=4,
    col=1,
)
fig.add_hline(y=0, line_dash="dash", line_color=COLORS["neutral"], row=4, col=1)

fig.update_layout(height=700, title=f"Options-implied features for {viz_symbol}", showlegend=False)
fig.update_yaxes(title_text="IV", row=1, col=1)
fig.update_yaxes(title_text="RR 25d", row=2, col=1)
fig.update_yaxes(title_text="Short/Long", row=3, col=1)
fig.update_yaxes(title_text="VRP", row=4, col=1)
show_plotly_with_alt(
    fig,
    (
        "Four stacked panels of options-implied features for one symbol, sharing a date "
        "axis across 2019 and 2020. The top panel plots at-the-money implied volatility, "
        "flat at a low level through 2019, spiking to its highest point of the window in "
        "the spring of 2020, falling back, then rising to a second lower peak in the "
        "autumn. The second fills the twenty-five delta risk reversal around a dashed "
        "zero line: positive for most of 2019 and early 2020, peaking as implied "
        "volatility does, then turning persistently negative from late summer 2020 with "
        "its deepest trough near the end. The third plots the implied volatility term "
        "slope, the ratio of short-dated to long-dated volatility, crossing a dashed "
        "reference line at one repeatedly across the window. The bottom panel fills the "
        "variance risk premium, implied volatility minus trailing realized volatility, "
        "around a dashed zero line: it sits above the line through much of 2019 and "
        "again from the autumn of 2020, dips below it in several short stretches, and "
        "carries one pronounced downward excursion in the spring of 2020 that reaches "
        "several times deeper than any other move in the panel."
    ),
)

# %% [markdown]
# **Interpretation**:
# - ATM IV spikes during market stress (the COVID crash, for one), which makes it a state variable
# - Positive risk reversal indicates elevated put demand (crash fear)
# - Term slope > 1 = inverted term structure (near-term event risk)
# - VRP > 0 is normal; extreme VRP signals attractive vol-selling opportunities
#
# **OI-based features** (dealer gamma, positioning) are conceptually valuable
# but require granular open interest data not available in this dataset. The
# aggregate Greeks approach in `22_sp500_equities_option_analytics_features`
# provides a workable proxy.
#
# **Surface stability**: If the quote selection policy (delta convention,
# interpolation method, maturity mapping) changes mid-backtest, all
# options-derived features become contaminated. Version your construction
# choices alongside the features.

# %% [markdown]
# ---
#
# # Summary
#
# ## Feature Families
#
# | Family | Key Features | Data Required | Primary Role |
# |--------|-------------|---------------|-------------|
# | **Carry** | Roll yield, funding rate | Multi-tenor futures, perps | Signal |
# | **Term Structure** | Slope, curvature | 3+ tenors | Signal/state |
# | **Cross-Asset** | Beta, lead-lag, relative value | Multi-asset prices | Signal |
# | **Options-Implied** | ATM IV, skew, VRP, term slope | Options chain | State |
#
# ## Implementation Rules
#
# 1. **Carry requires paired tenors**: Missing deferred contracts produce NaN carry
# 2. **Beta window matters**: 21d for fast regime, 63d for stable estimate
# 3. **Lead-lag can be spurious**: Validate that leads aren't staleness artifacts
# 4. **Options features need surface policy**: Version construction choices
# 5. **VRP is usually positive**: Negative VRP is extreme and informative
#
# ## Practical Takeaways
#
# 1. **Carry varies by asset class**: in the snapshot above, CL trades closer to
#    backwardation than GC does, which the printed shares either bear out on this sample
#    or do not; read them rather than this sentence, and read them as one sample's
#    behaviour rather than as the asset-class rule the received explanations describe
# 2. **Residual momentum isolates idiosyncratic return**: subtracting
#    $\beta_{t,a} r_{t,m}$ removes the market component from each ETF's return.
#    Whether residual momentum has higher IC than raw momentum is evaluated in
#    the `etfs` case study; this notebook demonstrates only the construction
# 3. **Options-implied features are state, not signal**: use them to condition
#    faster signals
# 4. **Surface stability is non-negotiable**: changing quote conventions
#    mid-backtest invalidates all options features
#
# ## Next Notebooks
#
# - `04_fundamentals_macro_calendar`: fundamentals, macro, calendar encodings
# - `case_studies/cme_futures`: full 30-product futures pipeline
# - `case_studies/sp500_equity_option_analytics`: full equity and options pipeline
# - `case_studies/sp500_options`: options straddle strategy pipeline

```