Caractéristiques de volatilité GJR-GARCH à un instant donné pour l’analyse des actions et des options

# 04_model_based_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]
# # S&P 500 Equity + Option Analytics: Model-Based Features
#
# Every feature in [`03_financial_features`](03_financial_features.ipynb) is a rule written in
# advance: a window, a difference, a rank. The feature built here is not. It is the output of a
# model whose parameters are **estimated from the data**, and that changes what has to be checked.
# A rule cannot leak by being fitted, because there is nothing to fit. A fitted parameter can, and
# the only thing that decides whether it did is which rows the estimation window contained.
#
# The model is a **GJR-GARCH(1,1)**, an equation for how a share's variance evolves from one
# session to the next. What it gives back is a **conditional volatility**: how much the share is
# expected to move over a session, worked out from what was known at the close of the session
# before it. Following the usual convention, the value is stamped on the session it describes, so
# `garch_cond_vol` at date $t$ is the square root of the variance the model gives session $t$ from
# the return and variance at $t-1$, annualized, and the recursion producing it reads nothing dated
# $t$ or later.
#
# **The recursion is only half of what decides whether a value is a forecast, and section C.1 is
# exact about the other half.** The parameters it runs on were estimated over a window, and on the
# dates inside that window a value is a retrospective transform rather than something that could
# have been computed on the day. On the dates after it the parameters predate the value entirely,
# which is why section F scores those dates alone and F2 draws them alone.
#
# That is the difference this notebook is built on. Stage 03's variance risk premium subtracts a
# *realized* volatility - an average of the last twenty sessions - from the option market's
# implied volatility, so it compares a forecast against a memory. Replacing the second half with
# a forecast makes both halves forecasts.
#
# ## Learning objectives
#
# - Re-estimate a volatility model on a fixed cadence over a share's own history, so that every
#   emitted value carries parameters fitted strictly before the session it speaks for
# - Declare that schedule where the feature is defined rather than inside the notebook, and assert
#   that an emitted value does not move when later observations are deleted rather than saying so
# - Pin every quantity the recursion starts from - the scaling, the starting variance and the
#   bounds it is clipped to - to the estimation window, because a library that recomputes them
#   over whatever series it is handed will read later data through them
# - Rebuild the option market's variance premium against a forecast rather than against a memory,
#   and test whether the change moves what the premium says about future returns
# - Test a difference between two dependent measurements directly, instead of inferring it from
#   two tests that were each about something else
#
# ## Book reference, prerequisites and artifacts
#
# Chapter 9, Section 9.3. Reads the daily share bars through `load_sp500_daily_bars()`, the
# implied volatility level and the variance risk premium from
# [`03_financial_features`](03_financial_features.ipynb)'s `features/financial.parquet`, and the
# primary label from [`02_labels`](02_labels.ipynb). Writes
# `features/model_based.parquet` with a `.digest.json` sidecar, read by every downstream model
# notebook through `load_modeling_dataset`, and by [`05_evaluation`](05_evaluation.ipynb).
# [`08_garch_volatility`](../../09_model_based_features/08_garch_volatility.ipynb) introduces the
# model itself.

# %%
"""S&P 500 Equity + Option Analytics: Model-Based Features."""

import numpy as np
import pandas as pd
import plotly.graph_objects as go
import polars as pl
from arch import arch_model
from ml4t.diagnostic.evaluation.stats import benjamini_hochberg_fdr
from ml4t.diagnostic.metrics import compute_ic_hac_stats, cross_sectional_ic_series

from case_studies.utils.artifact_digest import read_digest, value_digest
from case_studies.utils.artifact_quality import (
    label_universe,
    quality_report,
    render_quality_report,
)
from case_studies.utils.cv_window import fold_boundary_date, modeling_fold_boundaries
from case_studies.utils.temporal import (
    garch11_conditional_volatility,
    refit_boundaries,
    walk_forward_feature,
    write_model_based,
)
from case_studies.utils.warning_policy import apply_notebook_warning_policy
from data import load_sp500_daily_bars, load_sp500_options_surface
from utils.artifact_specs import load_setup_config, resolve_label_horizon
from utils.cv_splits import load_evaluation_config
from utils.data_quality import top_entities
from utils.paths import display_path, get_case_study_dir
from utils.style import ml4t_palette, show_plotly_with_alt

apply_notebook_warning_policy()

# %% [markdown]
# `MIN_OBS` is the shortest estimation window a fit is attempted on. A GJR-GARCH has four
# parameters and its likelihood is flat in the asymmetry term unless the window contains enough
# down sessions to identify it, so a year of sessions is the floor below which the notebook
# declines to fit rather than returning an unstable estimate. `MAX_SYMBOLS` caps the universe;
# `None` runs all of it. Both are read below and neither is re-assigned outside this cell, so the
# value a test injects is the value the notebook runs on.

# %% tags=["parameters"]
CASE_STUDY_ID = "sp500_equity_option_analytics"
MIN_OBS = 252
MAX_SYMBOLS = None
START_DATE = "2017-01-01"
END_DATE = "2021-12-31"

# %%
CASE_DIR = get_case_study_dir(CASE_STUDY_ID)
FEATURES_DIR = CASE_DIR / "features"
PANEL_KEY = ["timestamp", "symbol"]
FINANCIAL_PATH = FEATURES_DIR / "financial.parquet"
LABELS_DIR = CASE_DIR / "labels"

# %% [markdown]
# ## Configuration
#
# Three things are bound here and each decides something below: the label the section F diagnostic
# scores against and how far ahead it resolves, the refit schedule the estimation windows come
# from, and the date the holdout begins.
#
# **The schedule is declared in `config/setup.yaml`, not here.** An estimation window is part of a
# fitted feature's definition, so it lives where the definition lives. `model_based.gjr_garch`
# gives the burn-in and the refit cadence, and the comments there say what each decides and what
# it costs.
#
# **The folds are read for one purpose and it is presentational.** `modeling_fold_boundaries`
# still supplies the validation spans that section E's chart rules mark and that section F scores
# over, because a reader comparing this feature against a downstream model wants the same dates.
# No parameter is estimated per fold and no emitted value depends on one - that is the whole
# subject of section B.

# %%
SETUP = load_setup_config(CASE_STUDY_ID)
PRIMARY_LABEL = SETUP["labels"]["primary"]
# The outcome horizon, not the fold buffer: `labels.buffer` is 10D here because the purge has to
# clear the longest variant, while this label resolves in 5 sessions. The horizon is what sets
# section F's HAC bandwidth, because it is what makes consecutive daily IC values overlap.
_horizon = resolve_label_horizon(CASE_STUDY_ID, PRIMARY_LABEL, SETUP)
if not _horizon or not str(_horizon).endswith("D"):
    raise ValueError(f"Expected a daily label horizon for {PRIMARY_LABEL}, found {_horizon!r}")
LABEL_HORIZON = int(str(_horizon)[:-1])
LABEL_VARIANTS = [PRIMARY_LABEL, *SETUP["labels"].get("variants", [])]

canonical_splits = modeling_fold_boundaries(CASE_STUDY_ID, PRIMARY_LABEL)
if not canonical_splits:
    raise RuntimeError(f"Canonical {PRIMARY_LABEL} modeling folds are unavailable")
evaluation_config = load_evaluation_config(CASE_STUDY_ID)
HOLDOUT_START = fold_boundary_date(evaluation_config["holdout_start"])
HOLDOUT_END = fold_boundary_date(evaluation_config["holdout_end"])

_schedule = SETUP["model_based"]["gjr_garch"]
GARCH_BURNIN = int(_schedule["burnin"])
GARCH_REFIT_EVERY = int(_schedule["refit_every"])

print(
    f"{len(canonical_splits)} walk-forward folds come from the {PRIMARY_LABEL} label file. "
    "They mark the spans sections E and F report over; no parameter is fitted per fold."
)
print(
    f"GJR-GARCH schedule: {GARCH_BURNIN}-session burn-in, then re-estimated every "
    f"{GARCH_REFIT_EVERY} sessions on that security's whole history to date."
)
print(
    f"That label resolves {LABEL_HORIZON} sessions after the decision, so consecutive daily "
    f"scores in section F share {LABEL_HORIZON - 1} sessions of outcome and the standard errors "
    "there are corrected for it."
)
print(
    f"A fit is attempted only where the estimation window holds at least {MIN_OBS} sessions of "
    "returns."
)
print(f"Everything from {HOLDOUT_START} to {HOLDOUT_END} is the holdout.")

# %% [markdown]
# ## A. Why a fitted feature is different
#
# A feature from stage 03 is a function of past prices. Ask for the 20-session realized volatility
# at some date and the answer depends on twenty numbers and nothing else; run the calculation on a
# longer history and the value at that date does not move.
#
# A model-based feature is a function of *parameters estimated from* prices, and the parameters
# depend on every row in the estimation window. So the estimation window is part of the feature's
# information set, in a way that no window in stage 03 was. If the window used to estimate the
# parameters reaches past the date the feature is stamped with, then the feature at that date was
# computed with knowledge of what happened afterwards - and it will look far better than it can
# be, because the model was told the answer.
#
# The discipline that removes it has three parts, and section B, section C and section E carry one
# each: declare the refit schedule before any model runs, estimate on the sessions before each
# block and never on the block itself, and emit one value per date and security so there is
# nothing left for a consumer to join wrongly.
#
# There is one more subtlety, and it is the one a careless implementation leaves in place. Fixing
# the *parameters* is not enough. A recursion also needs a value to start from and a range to stay
# inside, and a library asked to run a recursion will derive both from whatever series it is
# handed. Section C measures what that costs here.

# %% [markdown]
# ## B. The schedule contract
#
# The parameters are re-estimated on a fixed cadence over each security's own history, and every
# emitted value comes from the last estimate made strictly before it. There is no estimation
# window per fold and no fold column on the artifact: a session gets one conditional volatility,
# and whichever fold selects that session later reads the same number.
#
# **The burn-in is paid out of the run-up, not out of training data.** seoa's panel opens
# 2017-01-03 and the oldest fold's training window opens 2018-01-04, so about 250 sessions
# precede the first fold. The declared burn-in of 252 is very nearly that, which is why it is
# 252 and not the 504 the longer-history case studies get - `config/setup.yaml::model_based`
# carries the measurement behind that choice.
#
# **The holdout is covered without being estimated on.** `freeze_after` stops re-estimation at
# the last session before `HOLDOUT_START`, so holdout sessions receive values from parameters
# fitted wholly before the seal. That is sound for the reason this stage rests on throughout:
# the model reads prices, not labels.
#
# **One schedule serves all five labels.** The old arrangement needed an argument that the
# primary label's folds were usable for the other four, because the estimation windows were the
# folds. A cadence over a security's own history has no fold in it, so that argument is no
# longer needed and the assertion that carried it is gone with the thing it protected.

# %%
# The block structure, on a security holding the full panel. `refit_boundaries` returns
# (fit_end_exclusive, emit_end_exclusive) pairs: the first block spends the burn-in, and each
# subsequent block is emitted from parameters fitted on everything before it.
_full_history = 1259
_blocks = refit_boundaries(_full_history, GARCH_BURNIN, GARCH_REFIT_EVERY)
print(f"A security with the full {_full_history}-session panel is fitted {len(_blocks)} times.")
print(f"  Burn-in, no value emitted: sessions 0 to {GARCH_BURNIN - 1}")
print(f"  First emitted session: {GARCH_BURNIN}")
print(f"  Sessions per block after the burn-in: {GARCH_REFIT_EVERY}")
pl.DataFrame(
    {
        "block": list(range(len(_blocks))),
        "fitted on sessions before": [b[0] for b in _blocks],
        "emits through session": [b[1] for b in _blocks],
    }
).head(8)

# %% [markdown]
# The spans below are the primary label's validation windows. Nothing is estimated per fold any
# more, so they decide no parameter; sections E and F report over them so that a reader can line
# this feature up against the downstream models that are scored on the same dates.

# %%
REPORTING_SPANS = [
    {
        "fold": s["fold"],
        "infer_start": fold_boundary_date(s["val_start"]),
        "infer_end": fold_boundary_date(s["val_end"]),
    }
    for s in canonical_splits
]
CHRONOLOGICAL = sorted(REPORTING_SPANS, key=lambda row: row["infer_start"])
print(
    f"{len(REPORTING_SPANS)} validation spans, {CHRONOLOGICAL[0]['infer_start']} to "
    f"{CHRONOLOGICAL[-1]['infer_end']}, and the holdout runs {HOLDOUT_START} to {HOLDOUT_END}."
)
pl.DataFrame(CHRONOLOGICAL)

# %% [markdown]
# ### F1. The schedule contract
#
# One row per refit block, for the first several blocks of a security holding the full panel. The
# dark bar is what the parameters were estimated on and the gold bar is the sessions those
# parameters emit. Every dark bar ends exactly where its own gold bar begins and never reaches
# past it, which is the property the whole revision is for. The dark bars grow because the walk
# expands: each fit sees everything before it, not a fixed window.

# %%
fit_colour, emit_colour = ml4t_palette(2, categorical=True)
_shown = _blocks[: min(8, len(_blocks))]
fig = go.Figure()
for row_index, (fit_end, emit_end) in enumerate(_shown):
    label = f"block {row_index}"
    emit_start = GARCH_BURNIN if row_index == 0 else _shown[row_index - 1][1]
    for span, colour, name in (
        ((0, fit_end), fit_colour, "estimated on"),
        ((emit_start, emit_end), emit_colour, "emits"),
    ):
        fig.add_trace(
            go.Scatter(
                x=list(span),
                y=[label, label],
                mode="lines",
                line=dict(color=colour, width=16),
                name=name,
                showlegend=row_index == 0,
            )
        )
fig.add_vline(x=GARCH_BURNIN, line_dash="dash", line_color=ml4t_palette(3)[2])
fig.update_layout(
    title="No parameter comes from the right of the sessions it speaks for",
    xaxis_title=f"Session index within the security; the dashed rule marks the "
    f"{GARCH_BURNIN}-session burn-in",
    height=340,
)
fig.update_yaxes(autorange="reversed")
show_plotly_with_alt(
    fig,
    "Timeline with one row per refit block, earliest at the top, drawn against session index "
    "within a single security rather than calendar dates. Each row carries a dark estimation bar "
    "starting at session zero and a gold emission bar immediately to its right. The dark bars "
    "lengthen row by row because the walk expands, while the gold bars stay the same short "
    "length, one refit interval each. A dashed vertical rule marks the end of the burn-in, where "
    "the first gold bar begins. No dark bar extends past the left edge of its own gold bar.",
)

# %% [markdown]
# ## C. The conditional volatility model
#
# GJR-GARCH(1,1) says next session's variance is a weighted sum of three things: a constant, how
# large the last session's move was, and what the variance already was.
#
# $$\sigma^2_t = \omega + (\alpha + \gamma \mathbb{1}_{r_{t-1}<0}) r^2_{t-1} + \beta \sigma^2_{t-1}$$
#
# Read the subscripts: $\sigma^2_t$ is built from $r_{t-1}$ and $\sigma^2_{t-1}$ and from nothing
# dated $t$ or later, which is what makes it a forecast of session $t$ rather than a description
# of it. Every column below inherits that. `garch_ivrv_spread` at $t$ compares the option market's
# implied volatility with a forecast of $t$ made before $t$ opened, and `garch_vol_surprise` at
# $t$ divides what the session actually did by what had been predicted for it - which is only a
# surprise because the denominator was fixed in advance.
#
# The $\gamma$ term is what the GJR variant adds to a plain GARCH: it lets a *negative* return
# raise the variance by more than a positive return of the same size. That asymmetry is well
# documented for equities and it is the reason the variant is worth the extra parameter here
# (Glosten, Jagannathan and Runkle, 1993). $\alpha + \gamma/2 + \beta$ is the model's
# **persistence** - how much of today's variance is still present tomorrow - and section D shows
# what it comes out at.
#
# Returns are multiplied by a hundred before fitting. The optimizer works on the variance of the
# series it is handed, and a daily equity return is a fraction of a percent whose square is four
# orders of magnitude smaller again; on that scale $\omega$ is small enough that the likelihood is
# flat in it to the optimizer's tolerance. Working in percent puts every parameter within a couple
# of orders of magnitude of one.
#
# ### C.1 Fit on the estimation window, then run forward without re-estimating
#
# The returns are taken on the **adjusted** close, `close` multiplied by `adj_factor`, for the
# reason `03_financial_features` gives at length: a four-for-one split reads as a three-quarter
# loss on the printed close, and a single such session is a squared return four orders of
# magnitude above anything else in the window. It would move $\omega$, it would be scored as an
# enormous surprise by the third feature below, and neither would be about the share.
#
# **Fixing the parameters is not enough on its own, and this is where a convenience method costs
# a value.** `arch`'s `fix` runs the recursion with parameters you supply, but it derives two
# other quantities from the series it is handed: the value the recursion starts from, and a pair
# of bounds every variance is clipped to. Handed the estimation window and the inference span
# together, it derives both from all of it - so a variance at a date inside the estimation window
# depends on returns that came after it, through the clipping range.
#
# So the recursion is run by `garch11_conditional_volatility`, the shared helper every case
# study that fits a volatility model here calls, and the three things it needs beyond the
# parameters - the scaling, the starting variance and the bounds - are all derived from the
# estimation window alone, using `arch`'s own functions for each. The recursion itself reads
# only $r_{t-1}$ and $\sigma^2_{t-1}$, so once it has started it never reads forward.
#
# **Every emitted value is now the same kind of quantity, and that is the change.** Under the
# refit schedule the parameters, the starting variance and the bounds behind a value all come
# from sessions strictly before the block that value sits in - on a training row exactly as much
# as on a validation row. There is no longer a retrospective half: the old design emitted values
# across the estimation window itself, and those carried parameters fitted from up to two years
# of their own future, which is precisely what made a training row and a validation row
# incomparable.
#
# The burn-in is where that guarantee is paid for. A security's first `burnin` sessions carry no
# value at all, because there is nothing to fit on yet, and section E reports how many securities
# that silences.
#
# C.2 asserts the property rather than describing it, by deleting later observations and checking
# that no retained value moves.


# %%
def training_filter_state(result, train_returns: pd.Series) -> tuple[float, float, tuple]:
    """Scaling, starting variance and clipping bounds, from the estimation window alone."""
    scale = float(result.model.scale)
    mu = float(result.params.get("mu", result.params.get("Const", 0.0)))
    residuals = train_returns.to_numpy(dtype=float) * scale - mu
    backcast = float(result.model.volatility.backcast(residuals))
    bounds = result.model.volatility.variance_bounds(residuals)
    return scale, backcast, (float(bounds[:, 0].min()), float(bounds[:, 1].max()))


# %% [markdown]
# One walk per security over that security's whole return history. `walk_forward_feature` spends
# the burn-in, fits on everything before the block it is about to emit, hands those parameters to
# `garch11_conditional_volatility`, and keeps only the block's own rows. `freeze_after` is the count of
# sessions strictly before the holdout, so the last estimate made is fitted on every pre-holdout
# session and on no holdout session.
#
# `fit` returns the four things the recursion needs and every one of them comes from the
# estimation prefix alone: the parameters, `arch`'s scaling, the backcast that seeds the first
# variance, and the clipping range. Deriving the last two from the emitted block instead is
# exactly the leak section C measures, so they are taken here and passed through.
#
# **A degenerate fit is recorded, not silently dropped.** A GJR-GARCH estimated on a short window
# can come back with `alpha + gamma < 0`, which says a larger down-shock *lowers* next session's
# variance. That is not a volatility model, and on this universe it is not rare: 19.0% of
# securities on their first block against 1.8% of fits across the whole walk, which section C.3
# reports. The recursion still runs - the clipping range keeps it on the positive reals - and the
# rate is published rather than hidden, because the burn-in that produces it is a declared
# schedule choice a reader is entitled to weigh.


# %%
ANNUALIZE = SETUP["evaluation"]["periods_per_year"] ** 0.5
# The specification is read, not written here, for the reason the schedule is: `o=1` is this case
# study's declared deviation from the shared GARCH(1,1)-Normal default, and a deviation stated in
# a notebook is one a reader has to find by reading the notebook. `config/setup.yaml::model_based`
# carries it and says what property of the data justifies it.
GARCH_KW = {k: _schedule[k] for k in ("mean", "vol", "p", "o", "q", "dist")}


def garch_walk(returns: pd.Series, freeze_after: int | None) -> tuple[pd.Series, list[dict]]:
    """Filter one security on the refit schedule; return the path and one record per fit."""
    diagnostics: list[dict] = []
    observations = returns.to_numpy(dtype=float).reshape(-1, 1)

    def fit(train: np.ndarray) -> dict:
        # The walk expands, so the prefix always starts at this security's first return and its
        # length is the index one past the last session the fit is allowed to see.
        train_returns = returns.iloc[: len(train)]
        result = arch_model(train_returns, **GARCH_KW).fit(disp="off", show_warning=False)
        scale, backcast, bounds = training_filter_state(result, train_returns)
        alpha = float(result.params.get("alpha[1]", 0.0))
        gamma = float(result.params.get("gamma[1]", 0.0))
        beta = float(result.params.get("beta[1]", 0.0))
        diagnostics.append(
            {
                # The first session these parameters speak for. The walk schedules on the
                # security's own observations, so a security that lists late refits on different
                # dates than its neighbours and this is the only honest key for section D.
                "emit_start": returns.index[len(train)],
                "n_train": len(train),
                "alpha": alpha,
                "gamma": gamma,
                "beta": beta,
                "persistence": alpha + gamma / 2 + beta,
                "degenerate": alpha + gamma < 0,
            }
        )
        return {
            "mu": float(result.params.get("mu", result.params.get("Const", 0.0))),
            "omega": float(result.params["omega"]),
            "alpha": alpha,
            "gamma": gamma,
            "beta": beta,
            "scale": scale,
            "backcast": backcast,
            "bounds": bounds,
        }

    def apply(model: dict, prefix: np.ndarray) -> np.ndarray:
        # `arch` estimates on returns it has scaled, so mu, omega and the bounds are all in the
        # scaled units and the recursion has to run there too. The volatility comes back divided
        # by the same factor, which returns it to percent, and the annualization follows.
        scale = model["scale"]
        sigma = (
            garch11_conditional_volatility(
                prefix[:, 0] * scale,
                mu=model["mu"],
                omega=model["omega"],
                alpha=model["alpha"],
                gamma=model["gamma"],
                beta=model["beta"],
                backcast=model["backcast"],
                bounds=model["bounds"],
            )
            / scale
        )
        return (sigma * ANNUALIZE / 100.0).reshape(-1, 1)

    values = walk_forward_feature(
        observations,
        timestamps=returns.index.to_numpy(),
        burnin=GARCH_BURNIN,
        refit_every=GARCH_REFIT_EVERY,
        fit=fit,
        apply=apply,
        n_features=1,
        window=None,
        freeze_after=freeze_after,
        on_fit_error="skip",
    )
    return pd.Series(values[:, 0], index=returns.index), diagnostics


def count_before(index: pd.Index, boundary: pd.Timestamp) -> int:
    """How many entries fall strictly before *boundary*.

    This is what ``walk_forward_feature`` wants for ``freeze_after``: the largest exclusive fit
    end it may use, so the last estimate before the seal sees every session before it and none
    after.
    """
    return int(index.searchsorted(boundary, side="left"))


# %% [markdown]
# The securities and their returns.
#
# **The series a model is estimated on is one security's, not one ticker's.** A ticker is
# reassigned after a merger or a spin-off and `adj_factor` restarts with the new security, so a
# return taken across the changeover divides one company's price by another's - and a single such
# session is a return large enough to move every parameter estimated from a window containing it.
# `03_financial_features` takes every one of its windows inside `sec_id` for the same reason, and
# this notebook does the same. The output is written on the ticker, because that is the key the
# downstream join uses, so the security is what the model is estimated inside and the ticker is
# what the result is stamped with.
#
# A row is emitted for a security in a given fold when the estimation window holds enough sessions
# and the optimizer returns; the count of those is reported below, because a security missing from
# one fold's rows is a security a downstream model has no conditional volatility for on those
# dates, and that is worth knowing before it turns up as a gap.
#
# The swap from security to ticker happens on an inner join against the bar panel, so every
# emitted value is named by the ticker its security actually traded under on that session. The
# row-count assertion beside it is what would catch a filtered value with no bar to name it,
# which would be a value nothing downstream could reach.
#
# **The universe is bounded before any fit, and the bound is the one `setup.yaml::universe`
# declares.** `eligibility_rule` is `sp500_with_options`, so a name qualifies by carrying an
# option surface. The share-bar extract is wider than that and the loader takes a `symbols=`
# argument nothing was passing, so this artifact carried a fitted conditional volatility for
# names the strategy can never hold - they have no implied volatility to rank on, so
# `garch_ivrv_spread` was null for every one of their rows while the other two columns were not.
# The roster is derived from the surface rather than typed out, and every roster name is checked
# against the share bars. How many names that comes to is reported rather than asserted:
# `n_assets` describes the production extract, and the reduced one CI runs on carries a handful by
# design. `tests/test_eoa_universe_roster.py` holds the declaration to the production data.

# %%
ENTITY = "sec_id"
# The surface is loaded for its roster alone; no column of it is read here. Both extracts are read
# whole rather than over the requested window, because which names have listed options is a
# property of the dataset. The dates bound the panel below.
_surface = load_sp500_options_surface()
full_bars = load_sp500_daily_bars()
ROSTER = sorted(_surface["symbol"].unique().to_list())
assert ROSTER, "the option-surface extract carries no names to rank"
priced = set(full_bars["symbol"].unique().to_list())
assert not set(ROSTER) - priced, f"no share bars for {sorted(set(ROSTER) - priced)}"
outside = sorted(priced - set(ROSTER))
DECLARED = SETUP["universe"]["n_assets"]
print(
    f"Universe {len(ROSTER)} names with an option surface"
    + ("" if len(ROSTER) == DECLARED else f" (a reduced extract; n_assets declares {DECLARED})")
    + f"; {len(outside)} priced names carry none and are excluded ({', '.join(outside)})"
)

window = pl.col("timestamp").is_between(
    pl.lit(START_DATE).str.to_date(), pl.lit(END_DATE).str.to_date()
)
bars = full_bars.filter(window & pl.col("symbol").is_in(ROSTER)).sort([ENTITY, "timestamp"])
bars = bars.with_columns((pl.col("close") * pl.col("adj_factor")).alias("adj_close"))
# One security trades under one ticker on one session. Without that, `sec_id` does not identify a
# price series and every window below would be taken across companies - which is the failure this
# whole section is arranged to avoid, so it is checked rather than assumed.
assert bars.select([ENTITY, "timestamp"]).is_duplicated().sum() == 0, (
    "a security carries more than one ticker on the same session, so sec_id does not identify a "
    "price series in this data"
)
securities = bars[ENTITY].unique().sort().to_list()
if MAX_SYMBOLS is not None:
    # Reduce by the rule every other stage reduces by, on the column they reduce on. This used to
    # take a prefix of the sorted sec_id list, which is the lexically first securities and has
    # nothing to do with which ones carry data: `05_evaluation` reduces to the most-observed
    # symbols, so a reduced 04 fitted volatility models for one universe while a reduced 05 scored
    # another, and the symbols only one side chose carried null features that ran clean.
    #
    # The selection is made on `symbol` and then mapped to `sec_id`, rather than counting sec_ids
    # directly, because `symbol` is what the artifact is keyed on and what every consumer reduces
    # by. The assertion above - one security per ticker per session - is what makes the mapping
    # one-to-one, so the two counts agree.
    kept_symbols = top_entities(bars, MAX_SYMBOLS, entity_col="symbol")
    securities = bars.filter(pl.col("symbol").is_in(kept_symbols))[ENTITY].unique().sort().to_list()
    print(f"Reduced run: {len(securities)} of {bars[ENTITY].n_unique()} securities")

log_returns: dict[int, pd.Series] = {}
for key, frame in (
    bars.select([ENTITY, "timestamp", "adj_close"])
    .filter(pl.col(ENTITY).is_in(securities))
    .partition_by(ENTITY, as_dict=True, maintain_order=True)
    .items()
):
    series = (
        frame.sort("timestamp")
        .with_columns((pl.col("adj_close").log().diff() * 100).alias("ret_pct"))
        .drop_nulls("ret_pct")
        .select(["timestamp", "ret_pct"])
        .to_pandas()
        .set_index("timestamp")["ret_pct"]
    )
    if not series.empty:
        log_returns[key[0]] = series

TICKER_OF = bars.select([ENTITY, "timestamp", "symbol"])
print(
    f"{len(log_returns)} securities under {bars['symbol'].n_unique()} tickers, "
    f"{bars['timestamp'].min()} to {bars['timestamp'].max()}"
)

# %% [markdown]
# ### C.2 What the convenience method would have cost
#
# The claim above is that handing `fix` the estimation window and the inference span together lets
# later returns reach an earlier date through the clipping range. The range is wide - the bounds
# sit six orders of magnitude either side of a local variance estimate - so it is entirely
# possible that it never binds and the whole concern is theoretical. That is a measurement, not an
# argument, so it is measured.
#
# The test isolates the clipping range from everything else. For a sample of securities on the
# first fold, `fix` is called twice with the same estimated parameters: once on the estimation
# window alone, once on the estimation window followed by the inference span. Both series start at
# the same session, so the starting variance - which `arch` reads from the first sessions only -
# is identical between them. The clipping range is then the one thing that differs, and any
# difference on the sessions the two calls share is what it did.

# %%
TRUNCATION_SAMPLE = 40
truncation_effect = []
for security in securities[:TRUNCATION_SAMPLE]:
    returns = log_returns.get(security)
    if returns is None or len(returns) <= GARCH_BURNIN + GARCH_REFIT_EVERY:
        continue
    blocks = refit_boundaries(len(returns), GARCH_BURNIN, GARCH_REFIT_EVERY)
    if len(blocks) < 2:
        continue
    # Cut INSIDE a block, never at a boundary. Truncating at a boundary leaves every retained
    # block with exactly the prefix it already had, so the two walks agree whatever the recursion
    # does inside a block and the check passes without testing anything. The cut below lands
    # mid-block by construction, and the assertion after this loop is what holds it there.
    boundaries = {end for _, end in blocks}
    # The last block can be a single observation, when the series ends one session past a refit.
    # It then has no session strictly inside it, and a cut clamped into it lands back on the
    # boundary that opens it - which is the case this assertion exists to refuse. Every earlier
    # block is a full `GARCH_REFIT_EVERY` wide by construction, so step back one.
    block_start, block_end = blocks[-1]
    if block_end - block_start < 2:
        block_start, block_end = blocks[-2]
    cut = min(block_start + GARCH_REFIT_EVERY // 2, block_end - 1)
    assert cut not in boundaries, "the truncation point fell on a refit boundary and tests nothing"
    full, _ = garch_walk(returns, None)
    short, _ = garch_walk(returns.iloc[:cut], None)
    shared = full.iloc[: len(short)].to_numpy()
    both = np.isfinite(shared) & np.isfinite(short.to_numpy())
    if not both.any():
        continue
    denom = np.maximum(np.abs(shared[both]), 1e-12)
    truncation_effect.append(float(np.max(np.abs(shared[both] - short.to_numpy()[both]) / denom)))

if truncation_effect:
    print(f"{len(truncation_effect)} securities truncated mid-block and re-walked")
    print(
        "largest relative move in a retained value when every later observation is deleted: "
        f"{max(truncation_effect):.3e}"
    )
    assert max(truncation_effect) < 1e-9, (
        "deleting later observations moved an earlier emitted value, so something in the walk "
        "reads past the row it is emitting"
    )
else:
    print("No security in the sample was long enough to truncate mid-block.")

# %% [markdown]
# The property is asserted rather than described: an emitted value does not move when
# observations after it are deleted. The cut is taken inside a block rather than at a refit
# boundary, because a boundary cut leaves every retained block with exactly the prefix it had and
# would agree no matter what the recursion did.
#
# What the convenience method would cost is worth keeping in view beside it: handing `arch`'s
# `fix` the estimation window and the inference span together moves the volatility it reports on
# a *shared* session by up to **0.9%** of that security's own largest value, measured across 38
# securities here. That is what `garch11_conditional_volatility` exists to avoid, and it is why
# the backcast and the clipping range are derived from the estimation prefix and passed in rather
# than left for the library to infer from whatever series it is handed.

# %%
filtered_parts: list[pl.DataFrame] = []
parameter_rows: list[dict] = []
too_short = 0

for security in securities:
    returns = log_returns.get(security)
    if returns is None or len(returns) <= GARCH_BURNIN:
        # A security listed too late to clear its own burn-in emits nothing. That is a statement
        # about the security rather than a failure, and it is counted so the coverage report
        # below can say how many rather than leaving a reader to infer it from a gap.
        too_short += 1
        continue
    path, records = garch_walk(returns, count_before(returns.index, pd.Timestamp(HOLDOUT_START)))
    emitted = path.dropna()
    if emitted.empty:
        too_short += 1
        continue
    filtered_parts.append(
        pl.DataFrame({"timestamp": emitted.index.values, "garch_cond_vol": emitted.to_numpy()})
        .with_columns(pl.col("timestamp").cast(pl.Date))
        .with_columns(pl.lit(security).alias(ENTITY))
    )
    parameter_rows.extend({ENTITY: security, **record} for record in records)

_filtered = pl.concat(filtered_parts)
garch = _filtered.join(TICKER_OF, on=[ENTITY, "timestamp"], how="inner").select(
    ["timestamp", "symbol", ENTITY, "garch_cond_vol"]
)
assert garch.height == _filtered.height, "a filtered value had no bar to name it"
parameters = pl.DataFrame(parameter_rows)
print(f"{len(filtered_parts):,} of {len(securities):,} securities emitted a value")
print(f"  never cleared the {GARCH_BURNIN}-session burn-in: {too_short}")
print(f"  fits across every security and refit: {parameters.height:,}")
print(
    f"{garch.height:,} filtered values over {garch[ENTITY].n_unique()} securities, "
    f"written under {garch['symbol'].n_unique()} tickers"
)

# %% [markdown]
# ## D. Fit stability across the schedule
#
# The method rests on re-estimating every `refit_every` sessions, and that is only worth its cost
# if the estimates actually move. If the same parameters come back block after block, the cadence
# bought nothing and one estimation would have done. If they swing, the feature carries whatever
# the expanding window happened to contain as much as it carries anything about the share, and a
# reader should treat it accordingly.
#
# What a consumer no longer has to do is join on a fold. Each session's value comes from the one
# estimate that preceded it, so the date and the security are the whole key.
#
# Persistence is the summary to read, and it is the sum of the other three: $\alpha + \gamma/2 +
# \beta$. A value near one means a shock to volatility decays slowly, so today's variance says
# something about tomorrow's - which is the property that makes a conditional volatility worth
# forecasting at all. A value near zero would mean the model has nothing to carry forward and the
# forecast is close to a constant.
#
# ### F3. Fit stability across the schedule
#
# One line per estimated parameter against the date the estimate first speaks for, drawn at the
# median across securities with a band spanning the middle half. The x-axis is calendar time
# rather than a fold index because the blocks are calendar time: a security refits every
# `refit_every` of its own sessions, so securities that list at different dates refit on different
# dates and the median is taken over whoever is estimated at each. All four parameters are
# dimensionless and bounded the same way, so they share an axis.

# %%
parameter_colours = ml4t_palette(4, categorical=True)
styles = [
    ("alpha", parameter_colours[0], "solid"),
    ("gamma", parameter_colours[1], "solid"),
    ("beta", parameter_colours[2], "dash"),
    ("persistence", parameter_colours[3], "solid"),
]
by_date = (
    parameters.with_columns(pl.col("emit_start").cast(pl.Date))
    .group_by("emit_start")
    .agg(
        pl.len().alias("securities"),
        *[
            expr
            for name, _, _ in styles
            for expr in (
                pl.col(name).median().alias(f"{name}_med"),
                pl.col(name).quantile(0.25).alias(f"{name}_lo"),
                pl.col(name).quantile(0.75).alias(f"{name}_hi"),
            )
        ],
    )
    .sort("emit_start")
)
x = by_date["emit_start"].to_list()
fig = go.Figure()
for name, colour, dash in styles:
    fig.add_scatter(
        x=x + x[::-1],
        y=by_date[f"{name}_hi"].to_list() + by_date[f"{name}_lo"].to_list()[::-1],
        fill="toself",
        fillcolor=colour,
        opacity=0.15,
        line=dict(width=0),
        hoverinfo="skip",
        showlegend=False,
    )
    fig.add_scatter(
        x=x,
        y=by_date[f"{name}_med"].to_list(),
        mode="lines",
        name=name,
        line=dict(color=colour, width=2, dash=dash),
    )
fig.update_layout(
    title="Every parameter moves with the window, so every refit moves the feature",
    xaxis_title="Date the estimate first speaks for; line is the median across securities, "
    "band the middle half",
    yaxis_title="Estimated parameter value",
    height=430,
)
show_plotly_with_alt(
    fig,
    "Four parameter series against calendar time, one per estimated GJR-GARCH coefficient plus "
    "persistence, each drawn as a median line across securities with a shaded band spanning the "
    "middle half of them. Persistence and beta run highest and track each other; gamma and alpha "
    "run lower. The bands are widest early, when few securities have cleared the burn-in and the "
    "estimation windows are shortest, and narrow as the expanding windows lengthen. No series is "
    "flat, which is the point: the refit cadence changes the estimates it produces.",
)

# %%
stability = (
    parameters.with_columns(pl.col("emit_start").cast(pl.Date).dt.year().alias("year"))
    .group_by("year")
    .agg(
        pl.len().alias("fits"),
        pl.col(ENTITY).n_unique().alias("securities"),
        pl.col("alpha").median().round(4).alias("median alpha"),
        pl.col("gamma").median().round(4).alias("median gamma"),
        pl.col("beta").median().round(4).alias("median beta"),
        pl.col("persistence").median().round(4).alias("median persistence"),
        (pl.col("degenerate").mean() * 100).round(1).alias("alpha+gamma < 0 %"),
    )
    .sort("year")
)
stability

# %% [markdown] tags=["results"]
# Read the table by year and the figure by shape. What the schedule buys is visible in whether
# median persistence moves between refits: if it were flat, one estimation would have done and
# the cadence would be waste. The parts are worth reading separately from the sum, because
# $\alpha$ and $\gamma$ can trade against each other while persistence barely moves - a share of
# the response to a shock shifting between the symmetric and asymmetric terms is a change in what
# the model says about *down* sessions specifically, and persistence alone will not show it.
#
# The last column is the one to weigh against the schedule. It is the share of fits in each year
# that came back with $\alpha + \gamma < 0$, which says a larger down-shock lowers next session's
# variance. That is not a volatility model, and it is a burn-in artifact: measured across this
# universe it runs about 19% on a security's first block against about 1.8% of fits over the
# whole walk, falling as the expanding window lengthens. `config/setup.yaml::model_based` records
# why the burn-in is 252 anyway.
#
# A downstream model reading this feature is reading a quantity whose generating parameters were
# re-estimated on a schedule, not a fixed transform - but every session carries exactly one such
# value, so the date and the symbol are the whole key.

# %% [markdown]
# ### F2. What the model inferred, over the sessions it was run over
#
# The feature itself, drawn against the quantity it is meant to improve on. The line is the median
# across securities of the conditional volatility the model reports, and beside it the median of
# stage 03's 20-session realized volatility on the same dates. Two restrictions apply and both
# matter. Only the **inference** spans are drawn, because the estimation windows are the sessions
# the parameters came from and the reader is owed a picture of the feature where it is out of
# sample. And only the **development** folds are drawn: the holdout fold's values are in the
# artifact because a later stage needs them, and a figure read during development does not get to
# look at them.
#
# What to look for is not that the two agree - they are different quantities, one an estimate of
# the next session and the other an average over the last twenty - but *when* they part company.
# A conditional volatility responds to a shock the session after it, while a twenty-session
# average takes a month to absorb it and another month to forget it.

# %%
financial = pl.read_parquet(FINANCIAL_PATH).select(
    [*PANEL_KEY, "iv_30_atm", "rv_20", "ivrv_spread"]
)
# The validation spans of the development folds. Nothing about the feature depends on them now;
# they are the dates the downstream models are scored on, so reporting over them is what lets a
# reader line this chart up against those. The holdout is excluded, as it is everywhere in
# development.
inference_only = (
    pl.concat(
        [
            garch.filter(
                (pl.col("timestamp") >= pl.lit(row["infer_start"], dtype=pl.Date))
                & (pl.col("timestamp") <= pl.lit(row["infer_end"], dtype=pl.Date))
            )
            for row in REPORTING_SPANS
        ]
    )
    .unique(subset=[*PANEL_KEY])
    .sort(PANEL_KEY)
)

daily_medians = (
    inference_only.join(financial, on=PANEL_KEY, how="left")
    .group_by("timestamp")
    .agg(
        pl.col("garch_cond_vol").median().alias("conditional"),
        pl.col("rv_20").median().alias("realized"),
    )
    .sort("timestamp")
)

series_colours = ml4t_palette(2, categorical=True)
fig = go.Figure()
for column, name, colour in (
    ("conditional", "GJR-GARCH conditional volatility", series_colours[0]),
    ("realized", "20-session realized volatility", series_colours[1]),
):
    fig.add_scatter(
        x=daily_medians["timestamp"].to_list(),
        y=daily_medians[column].to_list(),
        mode="lines",
        name=name,
        line=dict(color=colour, width=1.6),
    )
for row in REPORTING_SPANS:
    fig.add_vline(x=row["infer_start"].isoformat(), line_dash="dot", line_color=ml4t_palette(3)[2])
fig.update_layout(
    title="The forecast turns within days of a shock where the average takes weeks",
    xaxis_title="Median across securities, validation spans only; dotted rules mark span starts",
    yaxis_title="Annualized volatility",
    height=420,
)
show_plotly_with_alt(
    fig,
    "Two lines of median annualized volatility across securities over the two inference spans, "
    "2019 and 2020, with dotted rules at each validation span start. Through 2019 both run "
    "together between "
    "roughly 0.18 and 0.30. In late February 2020 both climb steeply; the conditional volatility "
    "peaks first, near 1.0, and falls back below 0.4 within about six weeks, while the "
    "twenty-session realized volatility peaks higher, near 1.2, holds that level for several "
    "weeks and then declines in steps through the summer. From July 2020 the realized line stays "
    "above the conditional one for most of the rest of the year, including a separate rise to "
    "about 0.42 in the autumn that the conditional line barely registers.",
)

# %% [markdown]
# ## E. Combine and emit
#
# Three columns go out, and the second is the one this notebook was written for.
#
# - `garch_cond_vol` - the conditional volatility itself, annualized.
# - `garch_ivrv_spread` - the implied volatility of the coming month less that forecast. Stage 03's
#   `ivrv_spread` subtracts a *realized* volatility from the same implied level, so the two differ
#   in the denominator alone and in nothing else.
# - `garch_vol_surprise` - the size of the session's actual move divided by the volatility that
#   had been forecast for it. The denominator is a standard deviation, and most sessions move less
#   than one, so a well-calibrated forecast puts the typical value below one rather than at it;
#   what the column is for is the tail, where a value of three or four says the session moved far
#   more than the model had any reason to expect.
#
# The return in `garch_vol_surprise` is the same adjusted log return the model was estimated on,
# taken inside the security for the reason C.1 gives and annualized to the denominator's units so
# the ratio compares like with like.

# %%
returns_panel = (
    bars.select([ENTITY, "timestamp", "adj_close"])
    .sort([ENTITY, "timestamp"])
    .with_columns((pl.col("adj_close").log().diff().over(ENTITY) * ANNUALIZE).alias("_ann_ret"))
    .select([ENTITY, "timestamp", "_ann_ret"])
)

model_based = (
    garch.join(returns_panel, on=[ENTITY, "timestamp"], how="left")
    .join(financial.select([*PANEL_KEY, "iv_30_atm"]), on=PANEL_KEY, how="left")
    .with_columns(
        (pl.col("_ann_ret").abs() / pl.col("garch_cond_vol")).alias("garch_vol_surprise"),
        (pl.col("iv_30_atm") - pl.col("garch_cond_vol")).alias("garch_ivrv_spread"),
    )
    .drop(["_ann_ret", "iv_30_atm", ENTITY])
    .select([*PANEL_KEY, "garch_cond_vol", "garch_ivrv_spread", "garch_vol_surprise"])
    .filter(pl.col("timestamp") <= pl.lit(HOLDOUT_END, dtype=pl.Date))
    .sort(PANEL_KEY)
)
FEATURE_COLS = ["garch_cond_vol", "garch_ivrv_spread", "garch_vol_surprise"]
assert model_based.select(PANEL_KEY).is_duplicated().sum() == 0, "duplicate panel key"
assert model_based["timestamp"].max() <= HOLDOUT_END, "a value was emitted past the holdout"

# %% [markdown]
# ### E.1 What each column covers
#
# The panel key is `timestamp` + `symbol`, and that is the whole key. A date and security carry
# one conditional volatility, from the last parameters estimated before that date, so there is
# nothing for a consumer to disambiguate and no fold column to join on. That is the change this
# revision makes and the assertion above is what holds it.
#
# The walk covers each security's whole history, which runs past the holdout's end for nothing's
# benefit, so the artifact is bounded at `holdout_end` above. Coverage is reported over the
# development spans and the holdout separately, which keeps every distribution printed during
# development to sessions development is allowed to see.

# %%
coverage = (
    model_based.with_columns(
        pl.when(pl.col("timestamp") >= pl.lit(HOLDOUT_START, dtype=pl.Date))
        .then(pl.lit("holdout"))
        .otherwise(pl.lit("development"))
        .alias("span")
    )
    .group_by("span")
    .agg(
        pl.len().alias("rows"),
        pl.col("symbol").n_unique().alias("symbols"),
        *[(pl.col(c).is_not_null().mean() * 100).round(1).alias(f"{c} %") for c in FEATURE_COLS],
    )
    .sort("span")
)
coverage

# %%
development = model_based.filter(pl.col("timestamp") < pl.lit(HOLDOUT_START, dtype=pl.Date))
pl.DataFrame(
    {
        "feature": FEATURE_COLS,
        "median": [development[c].median() for c in FEATURE_COLS],
        "std": [development[c].std() for c in FEATURE_COLS],
        "1st pct": [development[c].quantile(0.01) for c in FEATURE_COLS],
        "99th pct": [development[c].quantile(0.99) for c in FEATURE_COLS],
    }
)

# %% [markdown]
# ### E.2 Write the artifact
#
# Beside the parquet goes a small JSON file with the same name and a `.digest.json` suffix, the
# same way stage 03 writes its matrix. What it buys is provenance a file name cannot give: a model
# trained on these values records the hash, and a later run can tell whether the feature it is
# reading is the feature it trained on.
#
# It holds a hash computed over the values in the file; the number of rows; the columns that
# identify a row; and a hash of each input the values were built from. The hash is over the content
# rather than the file's bytes, so a re-write in a different row order leaves it alone and any
# changed value moves it.
#
# The two inputs recorded here are the share bars the model was estimated on and stage 03's matrix,
# which supplied the implied volatility the second feature subtracts from.

# %%
record = write_model_based(
    model_based,
    FEATURES_DIR / "model_based.parquet",
    keys=PANEL_KEY,
    feature_columns=FEATURE_COLS,
    time_column="timestamp",
    # No fold column. The parameters come off a refit schedule over the whole panel, so a row is
    # identified by its keys alone and there is no fold for a value to belong to. Passing None
    # rather than leaving the default is what makes `expected_folds` refused instead of ignored.
    fold_column=None,
    written_by=f"case_studies/{CASE_STUDY_ID}/04_model_based_features.py",
    inputs={
        # `sec_id` belongs in this digest: it decides which rows form one series, so a changed
        # security mapping moves every fitted value while the prices themselves stand still.
        "load_sp500_daily_bars": value_digest(
            bars.select([*PANEL_KEY, ENTITY, "close", "adj_factor"])
        ),
        "features/financial.parquet": read_digest(FINANCIAL_PATH)["digest"],
    },
)
print(f"Wrote {display_path(FEATURES_DIR / 'model_based.parquet')}, digest {record['digest']}")
print(f"{record['n_rows']:,} rows on {record['keys']}, read through load_modeling_dataset")

# %% [markdown]
# ## F. Incremental evaluation
#
# One question, asked on validation rows only: does the option market's variance premium say
# something different about future returns when its denominator is a forecast instead of a memory?
#
# The measure is the **information coefficient**: on each decision date, the rank correlation
# across securities between the feature and the return that follows. That gives one number per
# date; the average over dates is the feature's IC, and its standard error has to account for
# consecutive dates sharing outcome sessions, which is what the label horizon bound above is for.
#
# **This section selects nothing.** All three features were written to the artifact in E.2 before
# it ran, and the holdout is excluded from every number below. `05_evaluation` is where features
# are screened, and it runs that screen over the whole matrix.

# %%
label = pl.read_parquet(LABELS_DIR / f"{PRIMARY_LABEL}.parquet").rename(
    {PRIMARY_LABEL: "forward_return"}
)
scored = (
    model_based.join(inference_only.select(PANEL_KEY), on=PANEL_KEY, how="semi")
    .join(label.select([*PANEL_KEY, "forward_return"]), on=PANEL_KEY, how="inner")
    .sort(PANEL_KEY)
)
assert scored.filter(pl.col("timestamp") >= pl.lit(HOLDOUT_START, dtype=pl.Date)).is_empty(), (
    "the validation scoring frame reaches into the holdout"
)
print(f"{scored.height:,} scored rows over {scored['timestamp'].n_unique():,} decision dates")


# %%
# Two floors, and they count different things. `MIN_CROSS_SECTION` is a number of securities on
# one date, below which a rank correlation is noise rather than a measurement. `MIN_IC_DATES` is
# a number of dates, below which the HAC standard error over those correlations means nothing.
# One constant used to serve both, which read as deliberate and was not: the second comparison
# was measuring a series length against a universe size.
#
# The cross-section floor is clamped to the universe this run actually loaded. At full width
# `MAX_SYMBOLS` is None and the clamp does nothing, so production is unchanged; under a reduced
# run it is what keeps this section measuring something. A fixed floor of ten against a five-name
# reduction does not shrink the section, it empties it - every date falls short, every feature
# reports no IC, and the run stays green over a measurement it never made.
MIN_CROSS_SECTION = min(10, scored["symbol"].n_unique())
MIN_IC_DATES = 10


def ic_series(frame: pl.DataFrame, column: str) -> pl.DataFrame:
    """Per-date cross-sectional rank correlation of *column* against the forward return.

    A date with fewer than ``MIN_CROSS_SECTION`` securities comes back as a null rather
    than as a correlation over three names, and is dropped here rather than averaged in.
    """
    scored_frame = frame.select([*PANEL_KEY, column, "forward_return"]).drop_nulls()
    scored_frame = scored_frame.rename({column: "prediction"})
    series = cross_sectional_ic_series(
        scored_frame,
        scored_frame,
        date_col="timestamp",
        entity_col="symbol",
        method="spearman",
        min_obs=MIN_CROSS_SECTION,
    )
    return series.filter(pl.col("ic").is_not_null() & pl.col("ic").is_not_nan()).sort("timestamp")


ic_rows = []
for column in FEATURE_COLS:
    series = ic_series(scored, column)
    if len(series) < MIN_IC_DATES:
        print(
            f"{column}: only {len(series)} dates carry a cross-section of at least "
            f"{MIN_CROSS_SECTION} securities, and {MIN_IC_DATES} are needed for a standard "
            "error, so no IC is reported for it"
        )
        continue
    stats = compute_ic_hac_stats(series, label_horizon=LABEL_HORIZON)
    ic_rows.append(
        {
            "feature": column,
            "ic": stats["mean_ic"],
            "hac_t": stats["t_stat"],
            "p": stats["p_value"],
            "dates": len(series),
        }
    )

if ic_rows:
    fdr = benjamini_hochberg_fdr([r["p"] for r in ic_rows], alpha=0.05, return_details=True)
    for row, survives in zip(ic_rows, fdr["rejected"], strict=True):
        row["survives FDR"] = bool(survives)
    ic_table = pl.DataFrame(ic_rows).with_columns(
        pl.col("ic").round(4), pl.col("hac_t").round(2), pl.col("p").round(4)
    )
else:
    ic_table = pl.DataFrame(schema={"feature": pl.Utf8, "ic": pl.Float64})
ic_table

# %% [markdown]
# ### F4. Incremental IC by feature
#
# Each bar is one feature's mean IC with its HAC t-statistic printed past the bar end. A
# t-statistic of about two is the conventional threshold for calling a mean distinguishable from
# zero at all; the Benjamini-Hochberg correction applied above lowers the bar each individual
# feature has to clear, because three tests give three chances to find something by accident. A
# bar is drawn in the darker colour where its feature clears that corrected threshold, so a chart
# in which every bar is the lighter colour is one in which none of them did.

# %%
if ic_table.is_empty():
    print("No IC series survived the cross-section floor, so there is nothing to draw here.")
else:
    bar_colour, faded = ml4t_palette(2, categorical=True)
    ordered = ic_table.sort(pl.col("ic").abs(), descending=True)
    fig = go.Figure()
    fig.add_bar(
        x=ordered["ic"].to_list(),
        y=ordered["feature"].to_list(),
        orientation="h",
        marker_color=[bar_colour if s else faded for s in ordered["survives FDR"].to_list()],
        text=[f"HAC t {t:+.2f}" for t in ordered["hac_t"].to_list()],
        textposition="outside",
        cliponaxis=False,
    )
    fig.add_vline(x=0.0, line_width=1, line_color=ml4t_palette(3)[2])
    _span = max(abs(v) for v in ordered["ic"].to_list()) * 1.7
    fig.update_layout(
        title="All three sit near zero and none is distinguishable from it",
        xaxis_title=(f"Mean cross-sectional IC against the {LABEL_HORIZON}-session forward return"),
        xaxis=dict(range=[-_span, _span]),
        yaxis=dict(autorange="reversed"),
        margin=dict(l=170, r=70),
        height=340,
    )
    show_plotly_with_alt(