ٹریڈنگ کیس اسٹڈیز میں لکیری ماڈلز کا موازنہ

# 07_case_study_insights.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] tags=[]
# # Case Study Insights: Linear Models
#
# **Docker image**: `ml4t`
#
# **Purpose**: assemble the cross-case-study view of the linear family
# (OLS, Ridge, Lasso, ElasticNet). Per-case-study deep dives live in
# `case_studies/{cs}/13_model_analysis.py`; this notebook is the
# comparative view across all nine case studies.
#
# **Learning objectives**
#
# - For each case study, read the highest-IC complete linear configuration's
#   mean daily cross-sectional Spearman IC with its HAC interval on the primary label
# - Compare regularization families (OLS, Ridge, Lasso, ElasticNet) at
#   each case study's primary label
# - Inspect per-fold IC distributions and the ICIR diagnostic for
#   stability
# - Extend the comparison across labels (horizon view) and across the
#   classification ↔ regression metric symmetry
# - Inspect coefficient sign consistency, Lasso sparsity, and the overlap
#   of selected features across case studies
#
# **Book reference**: Section 11.6 - Case study insights.
#
# **Prerequisites**: each case study's `06_linear.py` pipeline has populated
# `run_log/registry.db` for the linear family. Teaching notebooks NB01-NB06
# cover the underlying techniques.

# %% tags=[]
"""Case Study Insights: Linear - cross-case-study aggregation from the registry."""

import sqlite3
import warnings

import joblib
import matplotlib.pyplot as plt
import numpy as np
import polars as pl

# ml4t.diagnostic dlopens cudart; load torch first so its bundled CUDA
# runtime wins. Same precedence pattern as case_studies/utils/model_analysis.py.
import torch  # noqa: F401
from IPython.display import Markdown, display
from matplotlib.colors import LinearSegmentedColormap
from scipy import stats
from sklearn.exceptions import InconsistentVersionWarning
from sklearn.metrics import roc_auc_score

from case_studies.research.population import retired_prediction_hashes
from case_studies.utils.analytics import (
    CASE_STUDY_IDS,
    PRIMARY_LABELS,
    SHORT_NAMES,
)
from case_studies.utils.insight_chapter import (
    SYMMETRY_TABLE_SCHEMA,
    discover_symmetry_pairs,
    plot_cross_cs_forest,
    plot_rolling_daily_ic,
)
from case_studies.utils.model_analysis import (
    load_predictions,
)
from utils.paths import get_case_study_dir, registry_readonly_uri
from utils.style import COLORS, show_with_alt

# %% tags=["parameters"]
FAMILY = "linear"

# %% [markdown] tags=[]
# ## Complete-result guard
#
# A model result is eligible for comparison only when no fold has an undefined
# IC and its validation-day coverage equals the maximum for that label. The
# local loader excludes predictions with any null fold IC and rejects shorter,
# partially populated results before ranking.


# %% tags=[]
METRICS_QUERY = """
    SELECT t.training_hash, t.family, t.config_name, t.label,
           p.prediction_hash, p.checkpoint_value, p.checkpoint_kind,
           pm.ic_mean, pm.ic_std, pm.ic_mean_daily, pm.ic_se_hac,
           pm.ic_ci_lo, pm.ic_ci_hi, pm.ic_t_hac, pm.ic_p_hac,
           pm.ic_n_days, pm.ic_hac_lag,
           pm.auc_mean_daily, pm.auc_se_hac, pm.auc_ci_lo, pm.auc_ci_hi,
           pm.auc_n_days,
           (SELECT COUNT(*) FROM fold_metrics fm
            WHERE fm.prediction_hash = p.prediction_hash AND fm.ic IS NULL) AS n_null_folds
    FROM training_runs t
    JOIN prediction_sets p ON p.training_hash = t.training_hash
    JOIN prediction_metrics pm ON pm.prediction_hash = p.prediction_hash
    WHERE t.family = ? AND p.split = 'validation'
"""
FINITE_METRIC_FIELDS = [
    "ic_mean_daily",
    "ic_se_hac",
    "ic_ci_lo",
    "ic_ci_hi",
    "ic_t_hac",
    "ic_p_hac",
    "ic_n_days",
    "ic_hac_lag",
]


# %% [markdown] tags=[]
# The loader preserves exact training and prediction identities so every later
# statistic and coefficient view can be traced to one registry-selected artifact.
#
# It also drops the generations a refit has retired. A case study that refits under
# a corrected input artifact publishes new prediction identities and a snapshot that
# supersedes the previous one; the retired rows stay in the registry, because the
# record of what was superseded is evidence rather than litter. That record lives in
# `official_populations`, one layer above the tables joined below, so a query that
# stops at `prediction_metrics` sees two rows where the study has one. It cannot
# separate them on the numbers either: a refit that changes only a declared input
# leaves the computation, and therefore every metric, bit-identical, so a `max` or an
# `ORDER BY created_at` would pick one arbitrarily and look like it had worked.


# %% tags=[]
def load_complete_metrics(
    case_study: str,
    *,
    label: str | None = None,
) -> pl.DataFrame:
    """Load nondegenerate linear metrics with maximum coverage per label."""
    db_path = get_case_study_dir(case_study) / "run_log" / "registry.db"
    query = METRICS_QUERY
    params: list[str] = [FAMILY]
    if label is not None:
        query += " AND t.label = ?"
        params.append(label)
    with sqlite3.connect(registry_readonly_uri(db_path), uri=True) as connection:
        connection.row_factory = sqlite3.Row
        rows = [dict(row) for row in connection.execute(query, params).fetchall()]
        retired = retired_prediction_hashes(connection)
    if not rows:
        return pl.DataFrame()
    df = pl.DataFrame(rows, infer_schema_length=None)
    return (
        df.filter(
            ~pl.col("prediction_hash").is_in(list(retired))
            & (pl.col("n_null_folds") == 0)
            & pl.all_horizontal(
                pl.col(column).is_not_null() & pl.col(column).is_finite()
                for column in FINITE_METRIC_FIELDS
            )
        )
        .with_columns(_max_n_days=pl.col("ic_n_days").max().over("label"))
        .filter(pl.col("ic_n_days") == pl.col("_max_n_days"))
        .drop("_max_n_days", "n_null_folds")
    )


# %% [markdown] tags=[]
# ### Resolving a tie for rank one
#
# Two rows can tie on mean daily IC for two different reasons, and only one of them is
# an ambiguity. `config_name` alone does not identify a configuration: a name can be
# reused across generations while the parameters behind it move, so collapsing two such
# rows would pick one by hash order. A refit that only re-declared an input, on the other
# hand, leaves the row bit-identical apart from its two hashes. Requiring agreement on
# every column except those two hashes separates the duplicate from the genuine
# ambiguity without reading the training spec here.
#
# There is a third case, and it is the one this chapter is about. Ridge at a shrinkage
# small enough to be numerically inert reproduces OLS, and the registry then holds two
# or three rows whose every recorded statistic agrees to the last digit and whose only
# difference is the name of the configuration. `crypto_perps_funding` at `fwd_ret_24h`
# carries exactly that: `ols`, `ridge_a0.001` and `ridge_a0.01` at one IC, one standard
# error and one day count. Nothing is ambiguous about the measurement, so the rows
# collapse to the first name alphabetically and the others are carried in
# `tied_configs` and printed, rather than the notebook refusing to draw the section. A
# tie that differs on any measured column is still a real ambiguity and still raises.
#
# The excluded set is written as "everything except the two identities a refit moves"
# rather than as a list of fields to compare, because a list is a second copy of
# `METRICS_QUERY` to keep in step: the first draft of one omitted `ic_mean` and
# `ic_std` while claiming to cover every recorded statistic.

# %% tags=[]
TIE_IDENTITY_EXCLUDED = frozenset({"training_hash", "prediction_hash"})
TIE_NAME_ONLY = frozenset({"config_name"})


def resolve_generation_tie(tied: pl.DataFrame, what: str) -> pl.DataFrame:
    """Collapse rows recording one measurement; raise when the tie is a real ambiguity.

    Adds `tied_configs`, which names the tied configurations when several of them were
    collapsed and is null otherwise.
    """
    no_tie = pl.lit(None, dtype=pl.Utf8).alias("tied_configs")
    if tied.height == 1:
        return tied.with_columns(no_tie)
    fields = [column for column in tied.columns if column not in TIE_IDENTITY_EXCLUDED]
    differing = [column for column in fields if tied[column].n_unique() > 1]
    names = sorted(tied["config_name"].unique().to_list()) if "config_name" in tied.columns else []
    if not differing:
        return tied.sort("prediction_hash").head(1).with_columns(no_tie)
    if set(differing) <= TIE_NAME_ONLY:
        return (
            tied.sort("config_name")
            .head(1)
            .with_columns(pl.lit(", ".join(names)).alias("tied_configs"))
        )
    raise RuntimeError(
        f"{what}: daily-IC rank one is ambiguous"
        + (f" between {', '.join(names)}" if names else "")
        + f"; the tied rows differ on {', '.join(differing)}"
    )


def report_name_only_ties(frame: pl.DataFrame, what: str) -> None:
    """Print the selections that were decided by name order because nothing else differed."""
    if "tied_configs" not in frame.columns:
        return
    ties = frame.filter(pl.col("tied_configs").is_not_null())
    if ties.is_empty():
        return
    print(f"{what}: tied on every recorded statistic, carried under the first name:")
    for row in ties.iter_rows(named=True):
        who = row.get("short_name") or row.get("case_study") or ""
        print(f"  {who} {row.get('label', '')}: {row['tied_configs']}")


# %% [markdown] tags=[]
# The primary-label collector applies the same completeness rule independently
# to every case study before selecting the highest mean daily IC.


# %% tags=[]
def collect_complete_rank1(case_studies: list[str]) -> pl.DataFrame:
    """Return the highest-IC complete linear result per primary label."""
    frames = []
    for cs in case_studies:
        label = PRIMARY_LABELS[cs]
        df = load_complete_metrics(cs, label=label)
        if df.is_empty():
            raise RuntimeError(f"{cs}: no complete linear result for {label}")
        best_ic = df["ic_mean_daily"].max()
        winner = df.filter(pl.col("ic_mean_daily") == best_ic)
        # A tie between two generations of one configuration is a duplicate, not an ambiguity:
        # a refit writes a second prediction set under a new training hash with bit-identical
        # scores, and every fx_pairs linear configuration carries two such generations.
        winner = resolve_generation_tie(winner, f"{cs}: primary linear")
        frames.append(
            winner.with_columns(
                case_study=pl.lit(cs),
                short_name=pl.lit(SHORT_NAMES[cs]),
            )
        )
    if len(frames) != len(case_studies):
        raise RuntimeError("primary linear coverage is incomplete")
    return pl.concat(frames, how="diagonal_relaxed")


# %% [markdown] tags=[]
# HAC inference is rebuilt from the exact selected prediction artifact. Daily
# correlations are sorted before the Bartlett Newey-West calculation, avoiding
# the fold-order dependence in several stored registry intervals.


# %% tags=[]
def _daily_ic_values(row: dict, target_column: str | None = None) -> np.ndarray:
    """Load chronological daily IC values against an explicit artifact target."""
    parquet = (
        get_case_study_dir(row["case_study"])
        / "run_log"
        / "predictions"
        / row["prediction_hash"]
        / "predictions.parquet"
    )
    schema = pl.scan_parquet(parquet).collect_schema()
    score = "y_score" if "y_score" in schema else "prediction"
    actual = target_column or ("y_true" if "y_true" in schema else "actual")
    if actual not in schema:
        raise RuntimeError(
            f"{row['case_study']}: prediction artifact lacks target column {actual!r}"
        )
    daily = (
        pl.scan_parquet(parquet)
        .select("timestamp", score, actual)
        .group_by("timestamp")
        .agg(
            pl.corr(score, actual, method="spearman").alias("ic"),
            (pl.col(score).is_finite() & pl.col(actual).is_finite()).sum().alias("n_obs"),
        )
        .filter((pl.col("n_obs") >= 5) & pl.col("ic").is_finite())
        .sort("timestamp")
        .collect()
    )
    return daily["ic"].to_numpy()


# %% tags=[]
def _bartlett_hac(values: np.ndarray, lag: int, case_study: str) -> dict[str, float]:
    """Return mean and Bartlett Newey-West inference for ordered observations."""
    n_days = len(values)
    if n_days < 3:
        raise RuntimeError(f"{case_study}: fewer than three daily IC observations")
    mean_ic = float(values.mean())
    residuals = values - mean_ic
    lag = min(lag, n_days - 1)
    long_run_sum = float(residuals @ residuals)
    for offset in range(1, lag + 1):
        weight = 1.0 - offset / (lag + 1.0)
        long_run_sum += 2.0 * weight * float(residuals[offset:] @ residuals[:-offset])
    variance = long_run_sum / n_days**2 * n_days / (n_days - 1)
    if variance < -1e-15:
        raise RuntimeError(f"{case_study}: negative HAC variance {variance}")
    variance = max(variance, 0.0)
    standard_error = float(np.sqrt(variance))
    if not np.isfinite(standard_error) or standard_error == 0.0:
        raise RuntimeError(f"{case_study}: degenerate HAC standard error")
    t_hac = mean_ic / standard_error
    critical = float(stats.t.ppf(0.975, n_days - 1))
    return {
        "ic_mean_daily": mean_ic,
        "ic_n_days": float(n_days),
        "ic_se_hac": standard_error,
        "ic_ci_lo": mean_ic - critical * standard_error,
        "ic_ci_hi": mean_ic + critical * standard_error,
        "ic_t_hac": t_hac,
        "ic_p_hac": float(2.0 * stats.t.sf(abs(t_hac), n_days - 1)),
        "ic_hac_lag": float(lag),
    }


# %% tags=[]
def chronological_hac(row: dict, target_column: str | None = None) -> dict[str, float]:
    """Recompute daily IC and Bartlett Newey-West inference in time order."""
    values = _daily_ic_values(row, target_column=target_column)
    return _bartlett_hac(values, int(row["ic_hac_lag"]), row["case_study"])


# %% [markdown] tags=[]
# Recomputed inference replaces only the selected rows, preserving their exact
# training and prediction identities for later coefficient and fold diagnostics.


# %% tags=[]
def replace_with_chronological_hac(
    selected: pl.DataFrame, target_column: str | None = None
) -> pl.DataFrame:
    """Bind selected rows to inference recomputed against an explicit target."""
    rows = []
    for row in selected.iter_rows(named=True):
        row.update(chronological_hac(row, target_column=target_column))
        rows.append(row)
    return pl.DataFrame(rows, infer_schema_length=None)


# %% [markdown] tags=[]
# Every grouped comparison fails closed when two artifacts tie for the highest daily
# IC, so row order never decides which configuration a section reports.


# %% tags=[]
def select_unique_best(frame: pl.DataFrame, groups: list[str]) -> pl.DataFrame:
    """Select one highest daily-IC row in every requested group.

    Raises on a tie whose rows differ on anything measured, because row order must never
    decide which configuration a section reports. Two ties resolve instead, and neither is
    a judgement call: two generations of one configuration, where a refit wrote a second
    prediction set under a new training hash with bit-identical scores, and two or three
    configurations whose every recorded statistic agrees and which differ only in name.
    `resolve_generation_tie` names the second kind in `tied_configs`.
    """
    winners = []
    for group in frame.partition_by(groups, maintain_order=True):
        best_ic = group["ic_mean_daily"].max()
        best = group.filter(pl.col("ic_mean_daily") == best_ic)
        identity = {column: group[0, column] for column in groups}
        best = resolve_generation_tie(best, str(identity))
        winners.append(best)
    return pl.concat(winners, how="diagonal_relaxed") if winners else pl.DataFrame()


# %% [markdown] tags=[]
# Fold diagnostics use the exact complete configurations selected above, rather
# than independently ranking a second time.


# %% tags=[]
def collect_selected_fold_ic(selected: pl.DataFrame) -> pl.DataFrame:
    """Load fold metrics for the exact primary-label selections."""
    frames = []
    for row in selected.iter_rows(named=True):
        db_path = get_case_study_dir(row["case_study"]) / "run_log" / "registry.db"
        with sqlite3.connect(registry_readonly_uri(db_path), uri=True) as connection:
            records = connection.execute(
                """
                SELECT fold_id, ic, ic_std, n_entities, rmse, mae
                FROM fold_metrics WHERE prediction_hash = ? ORDER BY fold_id
                """,
                (row["prediction_hash"],),
            ).fetchall()
        folds = pl.DataFrame(
            records,
            schema=["fold_id", "ic", "ic_std", "n_entities", "rmse", "mae"],
            orient="row",
        )
        if folds.is_empty() or folds["ic"].null_count() > 0:
            raise RuntimeError(f"{row['case_study']}: selected result has missing fold IC")
        frames.append(
            folds.with_columns(
                case_study=pl.lit(row["case_study"]),
                short_name=pl.lit(row["short_name"]),
            )
        )
    return pl.concat(frames, how="diagonal_relaxed") if frames else pl.DataFrame()


# %% [markdown] tags=[]
# Multi-label comparisons repeat the complete-coverage selection within every
# case-study and label pair.


# %% tags=[]
def collect_complete_multi_label(
    case_studies: list[str],
    label_resolver,
) -> pl.DataFrame:
    """Return the highest-IC complete result for each requested label."""
    frames = []
    for cs in case_studies:
        for label in label_resolver(cs):
            df = load_complete_metrics(cs, label=label)
            if df.is_empty():
                continue
            winner = select_unique_best(df, ["label"])
            frames.append(
                winner.with_columns(
                    case_study=pl.lit(cs),
                    short_name=pl.lit(SHORT_NAMES[cs]),
                )
            )
    return pl.concat(frames, how="diagonal_relaxed") if frames else pl.DataFrame()


# %% [markdown] tags=[]
# ## 1. Scope and Coverage
#
# The linear family is the baseline against which Ch12 (gradient boosting),
# Ch13 (deep learning for time series), and Ch14 (latent factors) are
# compared. Each case study's primary label is fixed in
# `PRIMARY_LABELS`; the primary metric is mean daily cross-sectional Spearman IC.
# Selection uses registry coverage and IC fields, while HAC inference is recomputed
# from each selected prediction artifact after sorting daily IC chronologically.

# %% tags=[]
coverage_rows = []
for cs in CASE_STUDY_IDS:
    primary = PRIMARY_LABELS[cs]
    df = load_complete_metrics(cs)
    if df.is_empty():
        coverage_rows.append(
            {
                "case_study": SHORT_NAMES[cs],
                "primary_label": primary,
                "n_labels": 0,
                "n_complete_configs_primary": 0,
            }
        )
        continue
    n_labels = df["label"].n_unique()
    n_configs_primary = df.filter(pl.col("label") == primary)["config_name"].n_unique()
    coverage_rows.append(
        {
            "case_study": SHORT_NAMES[cs],
            "primary_label": primary,
            "n_labels": n_labels,
            "n_complete_configs_primary": n_configs_primary,
        }
    )

coverage_df = pl.DataFrame(coverage_rows)
print("Complete linear coverage per case study (primary label and labels trained):")
coverage_df

# %% [markdown] tags=[]
# ## 2. Cross-CS Forest of Highest-IC Linear Configurations
#
# For each case study, the linear configuration with the highest
# mean daily IC on the primary label is plotted with its HAC interval.
# Filled markers indicate $|t_{HAC}| > 2$ (CI excludes zero); open
# markers indicate the CI overlaps zero.

# %% tags=[]
rank1 = replace_with_chronological_hac(collect_complete_rank1(CASE_STUDY_IDS))
report_name_only_ties(rank1, "Primary label")
print(
    "Highest-IC complete linear configuration per case study "
    "(primary label, mean daily IC ± HAC 95 % CI):"
)
rank1.select(
    "short_name",
    "label",
    "config_name",
    pl.col("ic_mean_daily").round(4).alias("ic"),
    pl.col("ic_ci_lo").round(4).alias("ci_lo"),
    pl.col("ic_ci_hi").round(4).alias("ci_hi"),
    pl.col("ic_t_hac").round(2).alias("t_hac"),
    pl.col("ic_n_days").cast(pl.Int64).alias("n_days"),
)

# %% tags=[]
fig, ax = plot_cross_cs_forest(
    rank1,
    family=FAMILY,
    title="Mean daily IC by case study, with HAC intervals",
)
ax.set_xlabel("Mean daily IC (HAC 95 % CI)")
show_with_alt(
    fig,
    "Forest plot of mean daily information coefficient with HAC intervals, one row "
    "per case study on its primary label, against a vertical line at zero.",
)

# %% tags=[]
significant = rank1.filter((pl.col("ic_ci_lo") > 0) | (pl.col("ic_ci_hi") < 0))
overlapping = rank1.filter(~((pl.col("ic_ci_lo") > 0) | (pl.col("ic_ci_hi") < 0)))
sig_names = ", ".join(significant.sort("short_name")["short_name"].to_list())
overlap_names = ", ".join(overlapping.sort("short_name")["short_name"].to_list())
display(
    Markdown(
        f"The HAC interval excludes zero for **{significant.height} of {rank1.height}** primary "
        f"labels ({sig_names}). It overlaps zero for {overlap_names}. The ETF result now selects "
        "the full-coverage Ridge fit; shorter or degenerate L1 fits are not eligible for this "
        "comparison. These are validation estimates read from the curated registries, not a new "
        "cross-case-study model selection exercise."
    )
)

# %% [markdown] tags=[]
# ## 3. Regularization Comparison
#
# How much does the choice of regularization family change the IC at
# the primary label? For each case study we take the highest-IC OLS,
# Ridge, Lasso, and ElasticNet configurations and compare their daily-
# mean daily IC with HAC intervals, then trace the Ridge regularization path
# with its HAC band. The axes carry the confidence level.


# %% tags=[]
def family_from_config(s: pl.Expr) -> pl.Expr:
    return (
        pl.when(s.str.starts_with("ols"))
        .then(pl.lit("ols"))
        .when(s.str.starts_with("ridge"))
        .then(pl.lit("ridge"))
        .when(s.str.starts_with("lasso"))
        .then(pl.lit("lasso"))
        .when(s.str.starts_with("enet"))
        .then(pl.lit("elastic_net"))
        .otherwise(pl.lit("other"))
    )


# %% tags=[]
reg_rows = []
for cs in CASE_STUDY_IDS:
    primary = PRIMARY_LABELS[cs]
    df = load_complete_metrics(cs, label=primary)
    if df.is_empty():
        continue
    df = df.filter(pl.col("ic_mean_daily").is_not_null())
    df = df.with_columns(
        reg_family=family_from_config(pl.col("config_name")),
        short_name=pl.lit(SHORT_NAMES[cs]),
        case_study=pl.lit(cs),
    )
    reg_rows.append(df)

reg_df = pl.concat(reg_rows, how="diagonal_relaxed") if reg_rows else pl.DataFrame()

# Rank-1 per (CS, regularization family)
reg_best = replace_with_chronological_hac(
    select_unique_best(reg_df, ["case_study", "reg_family"])
).select(
    "case_study",
    "short_name",
    "reg_family",
    "config_name",
    "prediction_hash",
    "ic_mean_daily",
    "ic_ci_lo",
    "ic_ci_hi",
    "ic_t_hac",
)

# Wide pivot for table
ic_pivot = reg_best.pivot(on="reg_family", index="short_name", values="ic_mean_daily").sort(
    "short_name"
)
print("Mean daily IC by regularization family (complete primary-label results):")
ic_pivot

# %% [markdown] tags=[]
# ### Figure: Regularization Family Comparison
#
# Missing families remain gaps rather than zeros in the grouped chart.


# %% tags=[]
def regularization_metric(
    frame: pl.DataFrame,
    case_studies: list[str],
    column: str,
) -> np.ndarray:
    """Align one regularization metric to the chart's case-study order."""
    lookup = dict(zip(frame["short_name"].to_list(), frame[column].to_list(), strict=True))
    return np.array([lookup.get(case_study, np.nan) for case_study in case_studies], dtype=float)


# %% tags=[]
families_present = ["ols", "ridge", "lasso", "elastic_net"]
plot_data = reg_best.filter(pl.col("reg_family").is_in(families_present)).sort(
    ["short_name", "reg_family"]
)

fig, ax = plt.subplots(figsize=(11, 5))
cs_order = sorted(plot_data["short_name"].unique().to_list())
x = np.arange(len(cs_order))
width = 0.2
fam_colors = {
    "ols": COLORS["neutral"],
    "ridge": COLORS["blue"],
    "lasso": COLORS["amber"],
    "elastic_net": COLORS["copper"],
}

for i, fam in enumerate(families_present):
    sub = plot_data.filter(pl.col("reg_family") == fam)
    ic = regularization_metric(sub, cs_order, "ic_mean_daily")
    lo = regularization_metric(sub, cs_order, "ic_ci_lo")
    hi = regularization_metric(sub, cs_order, "ic_ci_hi")
    err = np.vstack([ic - lo, hi - ic])
    ax.bar(
        x + (i - 1.5) * width,
        ic,
        width=width,
        yerr=err,
        capsize=2,
        color=fam_colors[fam],
        alpha=0.9,
        label=fam.replace("_", " ").title(),
    )

ax.set_xticks(x)
ax.set_xticklabels(cs_order, rotation=35, ha="right")
ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
ax.set_ylabel("Mean daily IC (HAC 95 % CI)")
ax.set_title("Mean daily IC by regularization family and case study")
ax.legend(frameon=False, fontsize=8, loc="best")
show_with_alt(
    fig,
    "Grouped bars of mean daily information coefficient by regularizer, one group per "
    "case study, with HAC error bars and a horizontal line at zero.",
)

# %% tags=[]
ridge_ols = reg_best.filter(pl.col("reg_family").is_in(["ridge", "ols"]))
ridge_ols_wide = ridge_ols.pivot(
    on="reg_family", index="short_name", values="ic_mean_daily"
).drop_nulls(["ridge", "ols"])
ridge_wins = ridge_ols_wide.filter(pl.col("ridge") >= pl.col("ols")).height
l1_cases = reg_best.filter(pl.col("reg_family").is_in(["lasso", "elastic_net"]))[
    "short_name"
].n_unique()
display(
    Markdown(
        f"Ridge matches or exceeds OLS in **{ridge_wins} of {ridge_ols_wide.height}** comparable "
        f"case studies. Complete Lasso or ElasticNet results are available for {l1_cases} case "
        "studies. Most family-to-family point-estimate shifts remain smaller than the displayed "
        "uncertainty intervals, so regularization choice rarely changes the primary conclusion."
    )
)

# %% [markdown] tags=[]
# ### Ridge regularization path
#
# How does mean daily IC change as Ridge $\alpha$ increases? A flat
# path means the conditioning of the feature matrix is already benign;
# a pronounced peak indicates an optimal shrinkage strength.

# %% tags=[]
ridge_path = replace_with_chronological_hac(
    reg_df.filter(pl.col("reg_family") == "ridge")
    .with_columns(
        alpha=pl.col("config_name").str.extract(r"ridge_a([\d.e+]+)", 1).cast(pl.Float64),
    )
    .filter(pl.col("alpha").is_not_null())
    .sort(["case_study", "alpha"])
)

palette = [
    COLORS["blue"],
    COLORS["copper"],
    COLORS["amber"],
    COLORS["positive"],
    COLORS["negative"],
    COLORS["slate"],
    COLORS["amber_light"],
    COLORS["neutral"],
    COLORS["blue_light"],
]
ridge_markers = ["o", "s", "D", "^", "v", "P", "X", "*", "h"]
ridge_styles = ["-", "--", "-.", ":", "-", "--", "-.", ":", "-"]

# %% tags=[]
fig, ax = plt.subplots(figsize=(10, 5))
cs_list_sorted = sorted(ridge_path["case_study"].unique().to_list())
for i, cs in enumerate(cs_list_sorted):
    sub = ridge_path.filter(pl.col("case_study") == cs).sort("alpha")
    if sub.height < 2:
        continue
    a = sub["alpha"].to_numpy()
    ic = sub["ic_mean_daily"].to_numpy()
    lo = sub["ic_ci_lo"].to_numpy()
    hi = sub["ic_ci_hi"].to_numpy()
    color = palette[i % len(palette)]
    name = SHORT_NAMES.get(cs, cs)
    ax.fill_between(a, lo, hi, color=color, alpha=0.10)
    ax.plot(
        a,
        ic,
        marker=ridge_markers[i],
        linestyle=ridge_styles[i],
        color=color,
        label=name,
        linewidth=1.5,
        markersize=4,
        alpha=0.9,
    )

ax.set_xscale("log")
ax.set_xlabel(r"Ridge $\alpha$ (log scale)")
ax.set_ylabel("Mean daily IC (HAC 95 % CI band)")
ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
ax.set_title("Mean daily IC against Ridge alpha, by case study")
ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), frameon=False, fontsize=8)
show_with_alt(
    fig,
    "Mean daily information coefficient against Ridge alpha on a log axis, one line per "
    "case study with a shaded HAC band, against a horizontal line at zero.",
)

# %% [markdown] tags=[]
# The confidence bands put the path curvature in context. An apparent optimum
# matters only when its change is large relative to uncertainty; across these
# panels, most neighboring Ridge settings remain statistically difficult to
# distinguish.

# %% [markdown] tags=[]
# ## 4. Persistence and Stability
#
# A positive mean daily IC averaged over the whole validation period can
# still come from a few good stretches. The primary persistence view is a
# three-month rolling mean of the daily IC, which traces how ranking quality
# evolves over time. Because case studies cover different validation windows,
# the comparison is clipped to the common validation period for ETFs and FX.
# (Figure 11.8 in the chapter.) Per-fold IC and ICIR below are secondary stability
# diagnostics.

# %% tags=[]
fig, _ = plot_rolling_daily_ic(
    ["etfs", "fx_pairs"],
    window=63,
    common_window=True,
    selected_prediction_hashes={
        row["case_study"]: row["prediction_hash"]
        for row in rank1.filter(pl.col("case_study").is_in(["etfs", "fx_pairs"])).iter_rows(
            named=True
        )
    },
    title="63-day rolling daily IC for ETFs and FX",
)
show_with_alt(
    fig,
    "Rolling 63-day information coefficient for the ETF and FX selected fits over their "
    "shared window, against a horizontal line at zero.",
)

# %% [markdown] tags=[]
# The rolling view tests whether a full-period average is broadly persistent or
# concentrated in a few stretches. Both paths cross zero here, so neither
# full-period average describes every validation regime. Both are validation
# estimates: the holdout period is untouched by this chapter.

# %% tags=[]
fold_df = collect_selected_fold_ic(rank1)
fold_summary = (
    fold_df.group_by(["case_study", "short_name"])
    .agg(
        n_folds=pl.col("ic").count(),
        mean=pl.col("ic").mean(),
        median=pl.col("ic").median(),
        std=pl.col("ic").std(),
        pct_positive=(pl.col("ic") > 0).mean(),
    )
    .sort("median", descending=True)
)
print("Per-fold IC summary for the highest-IC linear config (primary label):")
fold_summary

# %% [markdown] tags=[]
# ### Figure: Per-Fold IC Distribution

# %% tags=[]
order = rank1.sort("ic_mean_daily", descending=True)["short_name"].to_list()
present = [c for c in order if c in fold_df["short_name"].unique().to_list()]

fig, ax = plt.subplots(figsize=(11, 4.5))
data = [fold_df.filter(pl.col("short_name") == cs)["ic"].to_numpy() for cs in present]
positions = np.arange(len(present))
ax.boxplot(data, positions=positions, widths=0.55, showfliers=True)
for i, arr in enumerate(data):
    if len(arr):
        ax.scatter(np.full(len(arr), i), arr, alpha=0.5, s=14, color=COLORS["blue"])
ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
ax.set_xticks(positions)
ax.set_xticklabels(present, rotation=30, ha="right")
ax.set_ylabel("Per-fold Spearman IC")
ax.set_title("Per-fold Spearman IC of the selected linear fits")
show_with_alt(
    fig,
    "Box plots of per-fold Spearman information coefficient, one box per case study with "
    "the individual folds overlaid as points, against a horizontal line at zero.",
)

# %% tags=[]
most_positive = fold_summary.sort("pct_positive", descending=True).row(0, named=True)
least_variable = fold_summary.sort("std").row(0, named=True)
display(
    Markdown(
        f"**{most_positive['short_name']}** has the largest positive-fold share "
        f"({most_positive['pct_positive']:.0%}), while **{least_variable['short_name']}** has the "
        f"smallest fold-level standard deviation ({least_variable['std']:.4f}). Fold counts range "
        f"from {fold_summary['n_folds'].min()} to {fold_summary['n_folds'].max()}, so the boxplots "
        "are stability diagnostics rather than equal-precision estimates."
    )
)

# %% [markdown] tags=[]
# ### ICIR cross-CS bar
#
# The information ratio $\text{ICIR} = |\overline{\text{IC}}| / \sigma(\text{IC})$
# is a stability-adjusted signal-strength summary. With only a handful of folds per
# case study it is sample-starved and should be read alongside the
# mean daily IC and per-fold distribution, not as a standalone score.

# %% tags=[]
icir = fold_summary.with_columns(icir=(pl.col("mean").abs() / pl.col("std")).round(3)).sort(
    "icir", descending=True, nulls_last=True
)

fig, ax = plt.subplots(figsize=(9, 3.6))
y = np.arange(icir.height)
ax.barh(y, icir["icir"].to_numpy(), color=COLORS["blue"], alpha=0.9, height=0.55)
ax.set_yticks(y)
ax.set_yticklabels(icir["short_name"].to_list())
ax.invert_yaxis()
ax.set_xlabel("ICIR (|mean fold IC| / fold standard deviation)")
ax.set_title("ICIR of the selected linear fit, by case study")
show_with_alt(
    fig,
    "Horizontal bars of ICIR, the absolute mean fold information coefficient divided by "
    "its fold standard deviation, one bar per case study.",
)

# %% tags=[]
icir_leader = icir.row(0, named=True)
display(
    Markdown(
        f"**{icir_leader['short_name']}** has the largest fold ICIR "
        f"({icir_leader['icir']:.3f}). With only "
        f"{icir_leader['n_folds']} folds for that estimate, ICIR remains a descriptive "
        "stability diagnostic and does not replace the daily HAC interval."
    )
)

# %% [markdown] tags=[]
# ## 5. Multi-Label Horizon and Metric Symmetry

# %% [markdown] tags=[]
# ### 5a. Linear IC across regression labels per case study
#
# Several case studies trained the linear pipeline on multiple regression
# horizons. For each (case study, label), we take the highest complete
# mean daily IC across linear configurations and compare across horizons -
# this diagnoses whether the cross-sectional signal strengthens or
# weakens with the prediction window.


# %% tags=[]
def regression_labels(cs: str) -> list[str]:
    df = load_complete_metrics(cs)
    if df.is_empty():
        return []
    return [
        lbl
        for lbl in df["label"].unique().to_list()
        if lbl is not None
        and lbl.startswith("fwd_ret_")
        and "spot" not in lbl
        and "_win" not in lbl
        and "risk_adj" not in lbl
    ]


horizon_df = replace_with_chronological_hac(
    collect_complete_multi_label(CASE_STUDY_IDS, regression_labels)
)
report_name_only_ties(horizon_df, "Regression labels")
print(
    f"Highest-IC complete linear configuration per (case study, regression label): "
    f"{horizon_df.height} rows"
)
horizon_df.select(
    "short_name",
    "label",
    pl.col("ic_mean_daily").round(4).alias("ic"),
    pl.col("ic_ci_lo").round(4).alias("ci_lo"),
    pl.col("ic_ci_hi").round(4).alias("ci_hi"),
    pl.col("ic_t_hac").round(2).alias("t_hac"),
)

# %% tags=[]
HORIZON_DAYS = {
    "fwd_ret_5m": 5 / (6.5 * 60),
    "fwd_ret_15m": 15 / (6.5 * 60),
    "fwd_ret_60m": 60 / (6.5 * 60),
    "fwd_ret_8h": 1.0 / 3,
    "fwd_ret_24h": 1.0,
    "fwd_ret_1d": 1.0,
    "fwd_ret_5d": 5.0,
    "fwd_ret_10d": 10.0,
    "fwd_ret_21d": 21.0,
    "fwd_ret_1m": 21.0,
    "fwd_ret_3m": 63.0,
    "fwd_ret_1m_win": 21.0,
    "fwd_ret_risk_adj_5d": 5.0,
}

plot_horizon = horizon_df.with_columns(
    horizon_days=pl.col("label").replace_strict(HORIZON_DAYS, default=None).cast(pl.Float64),
).filter(pl.col("horizon_days").is_not_null())

multi_cs = (
    plot_horizon.group_by("short_name").len().filter(pl.col("len") >= 2)["short_name"].to_list()
)
plot_horizon = plot_horizon.filter(pl.col("short_name").is_in(multi_cs))

# %% tags=[]
if plot_horizon.height > 0:
    palette = [
        COLORS["blue"],
        COLORS["copper"],
        COLORS["amber"],
        COLORS["positive"],
        COLORS["negative"],
        COLORS["slate"],
        COLORS["amber_light"],
        COLORS["neutral"],
    ]

# %% tags=[]
if plot_horizon.height > 0:
    fig, ax = plt.subplots(figsize=(10, 5))
    cs_sorted = sorted(plot_horizon["short_name"].unique().to_list())
    markers = ["o", "s", "D", "^", "v", "P", "X", "*"]
    linestyles = ["-", "--", "-.", ":", "-", "--", "-.", ":"]
    for idx, cs in enumerate(cs_sorted):
        sub = plot_horizon.filter(pl.col("short_name") == cs).sort("horizon_days")
        if sub.height < 2:
            continue
        x = sub["horizon_days"].to_numpy()
        ic = sub["ic_mean_daily"].to_numpy()
        lo = sub["ic_ci_lo"].to_numpy()
        hi = sub["ic_ci_hi"].to_numpy()
        color = palette[idx % len(palette)]
        ax.fill_between(x, lo, hi, color=color, alpha=0.12)
        ax.plot(
            x,
            ic,
            marker=markers[idx % len(markers)],
            linestyle=linestyles[idx % len(linestyles)],
            color=color,
            label=cs,
            linewidth=1.6,
            markersize=6,
            alpha=0.9,
        )
    ax.set_xscale("log")
    ax.set_xlabel("Horizon (trading days, log scale)")
    ax.set_ylabel("Mean daily IC (HAC 95 % CI band)")
    ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_title("Mean daily IC against label horizon, by case study")
    ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), frameon=False, fontsize=8)
    show_with_alt(
        fig,
        "Mean daily information coefficient against label horizon in trading days on a log "
        "axis, one line per case study with a shaded HAC band, against a horizontal line at "
        "zero.",
    )

# %% tags=[]
display(
    Markdown(
        f"The horizon view covers **{plot_horizon['short_name'].n_unique()}** case studies with "
        "at least two comparable regression labels. Direction and magnitude vary by panel, and "
        "the HAC bands show that many within-panel orderings remain uncertain. Each horizon is a "
        "separate estimand; the chart does not treat a longer horizon as more observations."
    )
)

# %% [markdown] tags=[]
# ### 5b. Classification ↔ regression metric symmetry
#
# A linear model trained on a continuous return label has a regression
# score; the same modelling apparatus can be trained on the binary
# direction label of the same horizon. Across the case studies that
# carry binary direction labels, two questions are symmetric:
#
# - **Direction A** - the classification model's score, evaluated as
#   IC against the *continuous* return, asks whether the directional
#   classifier is also a useful *cross-sectional ranker*. This number
#   is on the registry as `prediction_metrics.ic_mean_daily` for
#   `task_type='classification'` rows.
# - **Direction B** - the regression model's score, evaluated as mean daily AUC
#   against the *binary* direction, asks whether the continuous
#   regression score is also a useful *cross-sectional binary classifier*. This number
#   is computed on the fly here from raw OOF predictions, since the
#   registry stores AUC only for classification training runs.
#
# We restrict to case studies with a binary direction label paired to a
# regression label for the same horizon; ternary direction labels (e.g.
# NASDAQ-100 `fwd_dir_15m`, Crypto `fwd_dir_8h_3c`) need a multi-class
# AUC and are out of scope for this brief subsection.

# %% tags=[]
# Discovered from the registry, not listed. A hand-written list cannot report what has
# gone missing from it: Chapter 12 carried the same literal, lost its
# `us_firm_characteristics` entry in a 2026-07-31 chapter-tree restore, and went on
# reporting a full count of itself. A ternary direction label is excluded by measuring
# its domain and naming it, which is also how the skip becomes visible here.
SYMMETRY_PAIRS, SYMMETRY_SKIPS = discover_symmetry_pairs(CASE_STUDY_IDS, FAMILY)
for line in SYMMETRY_SKIPS:
    print(f"  skipped {line}")
print(
    f"{sum(len(pairs) for pairs in SYMMETRY_PAIRS.values())} matched label pairs across "
    f"{len(SYMMETRY_PAIRS)} case studies"
)


# %% tags=[]
def _load_binary_label(cs: str, dir_label: str) -> pl.DataFrame:
    """Load one canonical binary direction label surface."""
    p = get_case_study_dir(cs) / "labels" / f"{dir_label}.parquet"
    if not p.exists():
        return pl.DataFrame()
    df = pl.read_parquet(p)
    return df.rename({dir_label: "y_dir"}).select("timestamp", "symbol", "y_dir")


# %% tags=[]
def _aligned_direction_scores(cs: str, reg_label: str, dir_label: str, best: dict) -> pl.DataFrame:
    """Join one exact regression prediction artifact to its binary labels."""
    preds = load_predictions(
        cs,
        family=FAMILY,
        label=reg_label,
        config_name=best["config_name"],
        checkpoint_value=best["checkpoint_value"],
        split="validation",
    )
    if preds.height == 0:
        return pl.DataFrame()
    preds = preds.filter(pl.col("prediction_hash") == best["prediction_hash"])
    if preds.is_empty():
        raise RuntimeError(f"{cs}: selected regression prediction artifact is missing")
    dir_df = _load_binary_label(cs, dir_label)
    if dir_df.is_empty():
        return pl.DataFrame()

    # Align timestamp + symbol dtypes for the join.
    if preds["timestamp"].dtype != dir_df["timestamp"].dtype:
        preds = preds.with_columns(pl.col("timestamp").cast(pl.Datetime("ms")))
        dir_df = dir_df.with_columns(pl.col("timestamp").cast(pl.Datetime("ms")))
    if preds["symbol"].dtype != dir_df["symbol"].dtype:
        preds = preds.with_columns(pl.col("symbol").cast(pl.Utf8))
        dir_df = dir_df.with_columns(pl.col("symbol").cast(pl.Utf8))

    return preds.join(dir_df, on=["timestamp", "symbol"], how="inner").filter(
        pl.col("y_dir").is_in([0, 1]) & pl.col("y_score").is_not_null()
    )


# %% tags=[]
def _mean_daily_auc(merged: pl.DataFrame) -> tuple[float, int] | None:
    """Compute unweighted mean daily AUC on dates containing both classes."""
    if merged.height == 0 or merged["y_dir"].n_unique() < 2:
        return None
    daily_auc = []
    for group in merged.sort("timestamp").partition_by("timestamp", maintain_order=True):
        if group["y_dir"].n_unique() < 2:
            continue
        daily_auc.append(roc_auc_score(group["y_dir"].to_numpy(), group["y_score"].to_numpy()))
    if not daily_auc:
        return None
    return float(np.mean(daily_auc)), len(daily_auc)


# %% tags=[]
def _direction_b_auc(cs: str, reg_label: str, dir_label: str) -> dict | None:
    """Highest-IC complete regression config to mean daily binary AUC."""
    reg_metrics = load_complete_metrics(cs, label=reg_label)
    if reg_metrics.is_empty():
        return None
    best = select_unique_best(reg_metrics, ["label"]).row(0, named=True)
    merged = _aligned_direction_scores(cs, reg_label, dir_label, best)
    auc_result = _mean_daily_auc(merged)
    if auc_result is None:
        return None
    auc, n_days = auc_result
    return {
        "case_study": cs,
        "short_name": SHORT_NAMES[cs],
        "reg_label": reg_label,
        "dir_label": dir_label,
        "reg_config": best["config_name"],
        "reg_ic_mean_daily": best["ic_mean_daily"],
        "reg_score_auc": auc,
        "n": merged.height,
        "n_days": n_days,
    }


# %% tags=[]
sym_rows = []
for cs, pairs in SYMMETRY_PAIRS.items():
    for reg_lbl, dir_lbl in pairs:
        # Direction A: classification model score to IC vs continuous return
        cls_metrics = load_complete_metrics(cs, label=dir_lbl)
        cls_ic = cls_ic_lo = cls_ic_hi = cls_t = None
        cls_config = None
        cls_auc = cls_auc_lo = cls_auc_hi = None
        if not cls_metrics.is_empty():
            cls_metrics = cls_metrics.filter(pl.col("ic_mean_daily").is_not_null())
            if not cls_metrics.is_empty():
                selected_cls = select_unique_best(cls_metrics, ["label"]).with_columns(
                    case_study=pl.lit(cs)
                )
                top = replace_with_chronological_hac(selected_cls, target_column="eval_actual").row(
                    0, named=True
                )
                cls_ic = top["ic_mean_daily"]
                cls_ic_lo = top.get("ic_ci_lo")
                cls_ic_hi = top.get("ic_ci_hi")
                cls_t = top.get("ic_t_hac")
                cls_config = top["config_name"]
                # The classifier's AUC against its own direction label - the book's
                # "Native AUC" column. It is in the same prediction_metrics row as the
                # IC above; until 2026-09-18 METRICS_QUERY did not select it, so this
                # notebook published a table it could not reproduce from its own output.
                cls_auc = top.get("auc_mean_daily")
                cls_auc_lo = top.get("auc_ci_lo")
                cls_auc_hi = top.get("auc_ci_hi")
        # Direction B: regression model score to AUC vs binary direction
        b = _direction_b_auc(cs, reg_lbl, dir_lbl)
        sym_rows.append(
            {
                "short_name": SHORT_NAMES[cs],
                "reg_label": reg_lbl,
                "dir_label": dir_lbl,
                "cls_config": cls_config,
                "cls_score_ic": cls_ic,
                "cls_score_ic_lo": cls_ic_lo,
                "cls_score_ic_hi": cls_ic_hi,
                "cls_score_ic_t": cls_t,
                "cls_score_auc": cls_auc,
                "cls_score_auc_lo": cls_auc_lo,
                "cls_score_auc_hi": cls_auc_hi,
                "reg_config": (b or {}).get("reg_config"),
                "reg_score_auc": (b or {}).get("reg_score_auc"),
                "n_b": (b or {}).get("n"),
                "n_b_days": (b or {}).get("n_days"),
            }
        )


# %% tags=[]

# A full schema, not schema_overrides. Discovery can legitimately return nothing - every
# declared pair skipped, or a registry with no classification runs for this family - and a
# frame built from an empty list with partial overrides has no string columns at all, so
# the selection below raises ColumnNotFoundError and the skip reasons this section exists
# to print never reach the reader. An empty frame with the right columns renders as an
# empty table, which is the correct answer.
sym_df = pl.DataFrame(sym_rows, schema=SYMMETRY_TABLE_SCHEMA)
print(
    "Native AUC (classification score on its own label), Direction A (classification "
    "score to IC) and Direction B (regression score to AUC):"
)
sym_df.select(
    "short_name",
    "reg_label",
    "dir_label",
    pl.col("cls_score_auc").round(4).alias("native_auc"),
    pl.col("cls_score_auc_lo").round(4).alias("native_lo"),
    pl.col("cls_score_auc_hi").round(4).alias("native_hi"),
    pl.col("cls_score_ic").round(4).alias("A_ic"),
    pl.col("cls_score_ic_lo").round(4).alias("A_lo"),
    pl.col("cls_score_ic_hi").round(4).alias("A_hi"),
    pl.col("cls_score_ic_t").round(2).alias("A_t"),
    pl.col("reg_score_auc").round(4).alias("B_auc"),
    "n_b_days",
)

# %% [markdown] tags=[]
# ### Figure: Direction A and Direction B side by side

# %% tags=[]
if sym_df.is_empty():
    raise RuntimeError(
        "no declared regression/direction pair qualified, so there is nothing to compare. "
        f"Skipped: {'; '.join(SYMMETRY_SKIPS) or 'nothing'}"
    )

fig, axes = plt.subplots(1, 2, figsize=(13, 4.0))
labels_y = [f"{r['short_name']} · {r['dir_label']}" for r in sym_df.iter_rows(named=True)]
y = np.arange(sym_df.height)

# Panel (a): Direction A, classification score IC vs continuous return
ax = axes[0]
ic = sym_df["cls_score_ic"].to_numpy()
lo = sym_df["cls_score_ic_lo"].to_numpy()
hi = sym_df["cls_score_ic_hi"].to_numpy()
ax.errorbar(
    ic,
    y,
    xerr=[ic - lo, hi - ic],
    fmt="o",
    color=COLORS["blue"],
    capsize=3,
    lw=1,
)
ax.axvline(0, color=COLORS["neutral"], lw=0.7, linestyle="--")
ax.set_yticks(y)
ax.set_yticklabels(labels_y)
ax.invert_yaxis()
ax.set_xlabel("Mean daily IC (HAC 95 % CI)")
ax.set_title("Classification score against the continuous return")

# Panel (b): Direction B, regression score AUC vs binary direction
ax = axes[1]
auc = sym_df["reg_score_auc"].to_numpy()
ax.scatter(auc, y, color=COLORS["copper"], s=60, zorder=3)
ax.axvline(0.5, color=COLORS["neutral"], lw=0.7, linestyle="--")
ax.set_yticks(y)
ax.set_yticklabels([])
ax.invert_yaxis()
ax.set_xlabel("Mean daily AUC (regression score)")
ax.set_title("Regression score against the binary direction")
ax.set_xlim(0.47, 0.53)

show_with_alt(
    fig,
    "Two panels sharing one row of case-study labels. Left: the classification score's "
    "mean daily information coefficient against the continuous return, with HAC error "
    "bars and a line at zero. Right: the regression score's mean daily AUC against the "
    "binary direction, with a line at one half.",
)

# %% tags=[]
direction_a_clear = sym_df.filter((pl.col("cls_score_ic_lo") > 0) | (pl.col("cls_score_ic_hi") < 0))
max_auc_gap = float((sym_df["reg_score_auc"] - 0.5).abs().max())
display(
    Markdown(
        f"Direction A's HAC interval excludes zero in **{direction_a_clear.height} of "
        f"{sym_df.height}** comparisons. Direction B's mean daily AUC remains within "
        f"**{max_auc_gap:.3f}** of chance across the same pairs. Cross-sectional rank "
        "correlation "
        "therefore does not collapse to binary-direction discrimination. The Ch16 equal-weight "
        "baseline uses ranked predictions, while direction-aware classifiers can still inform "
        "long-short construction."
    )
)

# %% [markdown] tags=[]
# ## 6. Coefficient Analysis
#
# The linear training pipeline stores a fitted pipeline per fold under the
# training run's `models/` directory, so the coefficients below are read out of
# those fold models. They are coefficients of the standardized features: the
# saved preprocessor imputes and scales and leaves the feature order alone, so
# they are comparable across features within a fit.
#
# A stored model was pickled by the scikit-learn that trained it, which is older
# than the one running this notebook, so unpickling one raises
# `InconsistentVersionWarning` once per fold. The loader silences that warning and
# nothing else, because it reads `coef_` and `feature_names` and never calls
# `predict`: an array attribute and a tuple of strings survive the version gap that
# the warning is about. A fit whose coefficient count stops matching its feature
# names raises instead of being read.


# %% tags=[]
def load_coefficients(cs: str, training_hash: str, config_name: str) -> pl.DataFrame | None:
    """Read the per-fold coefficients of one registry-tracked training run.

    One row per (fold, feature). Returns None when the run kept no fold models, which is
    what a run log holding only prediction sets looks like.
    """
    models_dir = get_case_study_dir(cs) / "run_log" / "training" / training_hash / "models"
    fold_paths = sorted(models_dir.glob("fold_*.joblib")) if models_dir.exists() else []
    rows = []
    for fold_path in fold_paths:
        with warnings.catch_warnings():
            warnings.filterwarnings(
                "ignore", category=InconsistentVersionWarning, module="sklearn.base"
            )
            payload = joblib.load(fold_path)
        names = list(payload["feature_names"])
        coef = np.asarray(payload["model"].coef_, dtype=float).ravel()
        if coef.size != len(names):
            raise RuntimeError(
                f"{cs} {training_hash}: {fold_path.name} has {coef.size} coefficients against "
                f"{len(names)} feature names. The stored preprocessor changed the column set, "
                "so the two cannot be paired."
            )
        fold = int(fold_path.stem.removeprefix("fold_"))
        rows.extend(
            {"fold": fold, "feature": name, "coefficient": float(value), "config_name": config_name}
            for name, value in zip(names, coef, strict=True)
        )
    if not rows:
        return None
    return pl.DataFrame(
        rows,
        schema={
            "fold": pl.Int64,
            "feature": pl.Utf8,
            "coefficient": pl.Float64,
            "config_name": pl.Utf8,
        },
    )


# %% [markdown] tags=[]
# ### 6a. Sign consistency across folds
#
# For each case study's highest-IC linear configuration on the primary label,
# sign consistency is the share of a feature's non-zero folds that agree on its
# sign. A regularized fit sets many coefficients to exactly zero, and a zero has
# no sign, so the zero folds are excluded from the fraction and a feature that is
# zero in every fold is excluded from the statistic altogether. Counting those
# zeros as agreement would score a fit that selected almost nothing as the most
# stable one in the table; how much a fit zeroes out is section 6b's subject.
# `n_active` is how many features the mean is taken over, and
# `SIGN_STABLE_MIN` is the consistency a feature has to reach to count as
# sign-stable.

# %% tags=[]
ZERO_TOL = 1e-10
SIGN_STABLE_MIN = 0.8

sign_rows = []
for row in rank1.iter_rows(named=True):
    cs = row["case_study"]
    cfg = row["config_name"]
    coef_df = load_coefficients(cs, row["training_hash"], cfg)
    if coef_df is None:
        continue
    coefs = coef_df.filter(pl.col("config_name") == cfg)
    if coefs.is_empty():
        continue
    stats = (
        coefs.group_by("feature")
        .agg(
            n_pos=(pl.col("coefficient") > ZERO_TOL).sum(),
            n_neg=(pl.col("coefficient") < -ZERO_TOL).sum(),
            n_folds=pl.col("coefficient").count(),
        )
        .with_columns(n_signed=pl.col("n_pos") + pl.col("n_neg"))
    )
    active = stats.filter(pl.col("n_signed") > 0).with_columns(
        consistency=pl.max_horizontal("n_pos", "n_neg") / pl.col("n_signed"),
    )
    if active.is_empty():
        continue
    sign_rows.append(
        {
            "short_name": row["short_name"],
            "config_name": cfg,
            "n_features": stats.height,
            "n_active": active.height,
            "mean_consistency": float(active["consistency"].mean()),
            "pct_active_sign_stable": float(
                (active["consistency"] >= SIGN_STABLE_MIN).mean() * 100
            ),
        }
    )

sign_df = pl.DataFrame(
    sign_rows,
    schema={
        "short_name": pl.Utf8,
        "config_name": pl.Utf8,
        "n_features": pl.Int64,
        "n_active": pl.Int64,
        "mean_consistency": pl.Float64,
        "pct_active_sign_stable": pl.Float64,
    },
).sort("mean_consistency", descending=True)
print("Coefficient sign consistency for the highest-IC linear configuration (primary label):")
sign_df

# %% tags=[]
if sign_df.is_empty():
    display(
        Markdown(
            "No chart: no selected linear fit in this run log has stored fold models, so "
            "there are no coefficients to count signs over."
        )
    )
else:
    fig, ax = plt.subplots(figsize=(10, 3.5))
    plot_sign = sign_df.sort("mean_consistency", descending=True)
    y = np.arange(plot_sign.height)
    ax.barh(
        y,
        plot_sign["mean_consistency"].to_numpy(),
        color=COLORS["blue"],
        alpha=0.9,
        height=0.5,
    )
    ax.set_yticks(y)
    ax.set_yticklabels(plot_sign["short_name"].to_list())
    ax.invert_yaxis()
    ax.set_xlim(0.5, 1.0)
    ax.axvline(SIGN_STABLE_MIN, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_xlabel("Mean sign consistency over a fit's active features")
    ax.set_title("Sign consistency of the selected linear fits")
    show_with_alt(
        fig,
        "Horizontal bars of mean fold-level coefficient sign consistency, one per case "
        "study, against a dashed line at the sign-stable threshold.",
    )

# %% tags=[]
# Counted on both paths, so the closing takeaways quote a zero rather than an absent name.
stable_sign_count = sign_df.filter(pl.col("mean_consistency") >= SIGN_STABLE_MIN).height
if sign_df.is_empty():
    display(
        Markdown(
            "No selected linear fit has stored fold models in this run log, so sign "
            "consistency cannot be measured here. The diagnostic reads the coefficients out "
            "of the per-fold models the training pipeline saves, which a run log holding "
            "only prediction sets does not carry."
        )
    )
else:
    sign_leader = sign_df.row(0, named=True)
    display(
        Markdown(
            f"Mean sign consistency reaches {SIGN_STABLE_MIN:.1f} in **{stable_sign_count} of "
            f"{sign_df.height}** selected fits, taken over the features each fit leaves "
            f"non-zero. **{sign_leader['short_name']}** is highest at "
            f"{sign_leader['mean_consistency']:.2f}. Sign stability and predictive ranking are "
            "different diagnostics: a stable coefficient direction does not by itself establish "
            "a nonzero out-of-sample IC."
        )
    )

# %% [markdown] tags=[]
# ### 6b. Lasso sparsity
#
# Lasso drives coefficients exactly to zero. The fraction of zero
# coefficients per case study summarizes how aggressively Lasso prunes
# features at the highest-IC $\alpha$ for the Lasso family.

# %% tags=[]
sparsity_rows = []
for cs in CASE_STUDY_IDS:
    label = PRIMARY_LABELS.get(cs)
    if label is None:
        continue
    lasso_metrics = load_complete_metrics(cs, label=label).filter(
        pl.col("config_name").str.starts_with("lasso")
    )
    if lasso_metrics.is_empty():
        continue
    lasso_selected = select_unique_best(lasso_metrics, ["label"]).row(0, named=True)
    lasso_config = lasso_selected["config_name"]
    lasso = load_coefficients(cs, lasso_selected["training_hash"], lasso_config)
    if lasso is None:
        continue
    n_total = lasso.height
    n_zero = lasso.filter(pl.col("coefficient").abs() <= ZERO_TOL).height
    always_zero = (
        lasso.group_by("feature")
        .agg(all_zero=(pl.col("coefficient").abs() <= ZERO_TOL).all())
        .filter(pl.col("all_zero"))
        .height
    )
    n_features = lasso["feature"].n_unique()
    sparsity_rows.append(
        {
            "short_name": SHORT_NAMES[cs],
            "config_name": lasso_config,
            "n_features": n_features,
            "zero_fraction": round(n_zero / n_total, 3),
            "features_always_zero": always_zero,
            "pct_always_zero": round(always_zero / n_features * 100, 1),
        }
    )


# %% tags=[]
sparsity_df = pl.DataFrame(
    sparsity_rows,
    schema={
        "short_name": pl.Utf8,
        "config_name": pl.Utf8,
        "n_features": pl.Int64,
        "zero_fraction": pl.Float64,
        "features_always_zero": pl.Int64,
        "pct_always_zero": pl.Float64,
    },
)
sparsity_df = sparsity_df.sort("zero_fraction", descending=True)
print("Lasso sparsity per case study (at the highest-IC Lasso $\\alpha$):")
sparsity_df

# %% tags=[]
if sparsity_df.is_empty():
    display(
        Markdown(
            "No case study in this run log pairs a complete Lasso result on its primary label "
            "with a stored coefficient artifact, so there is nothing to compare here. Sparsity "
            "is read off the coefficients themselves, not off any metric t