多重検定を考慮したCME先物の単変量特徴量評価

# Feature Evaluation - CME Futures


# Feature Evaluation - CME Futures

One screening pass over every candidate feature the two upstream notebooks built,
against the forward return the case study is trying to predict. For each feature it
asks three questions in order. Can the value be trusted at the moment a decision is
made? Does it sort next week's returns better than chance would? And is that sorting
the same in every year of the walk-forward test, or the product of one good period?

The answers become one record per feature. This notebook chooses nothing: it does not
fit a model, it does not pick the feature set the training notebooks use, and it says
nothing about whether a strategy built on these features would make money. Chapter 7
is explicit that screening features one at a time is necessary and not sufficient.

**What it reads**: `features/financial.parquet` (built by `03_financial_features`),
`features/model_based.parquet` (built by `04_model_based_features`), the forward-return
labels under `labels/` (built by `02_labels`), and `config/setup.yaml` for the label
names, the walk-forward geometry and the date the holdout begins.

**What it writes**: `evaluation/triage_ledger.parquet`, one row per feature carrying
its decision and the evidence behind it, and `evaluation/ic_timeseries.parquet`, the
per-date correlation series the figures here are drawn from.

**Learning objectives**:

- Measure how well a feature sorts next week's returns, by correlating the rank of the
  feature with the rank of the return across the products quoted on each decision date,
  then averaging that correlation over dates.
- Widen the error bars on that average to allow for consecutive dates being computed
  from price windows that overlap, so the daily measurements are not independent of
  each other and there are fewer of them than there appear to be.
- Correct for having tested every candidate feature at once rather than one, using
  Benjamini-Hochberg control of the share of announced findings that are false.
- Split a feature's results by walk-forward validation year, to tell one that works
  throughout the sample from one carried by a single period.
- Record a keep, revisit or drop decision per feature - the book calls these PROCEED,
  REVISE and STOP - together with the evidence that produced it.

**Book reference**: Chapter 7, Section 7.3 (Univariate feature-label evaluation) and
Section 7.4 (Search accounting and multiple testing). Section 8.6 covers controlling
the search once features are combined rather than screened one at a time.

**Prerequisites**: run `03_financial_features` and `04_model_based_features` first.

```python
"""Feature Evaluation - CME Futures.

Screens the financial and model-based features against the primary forward-return
label and writes one triage record per feature.
"""

import re
from datetime import date

import numpy as np
import pandas as pd
import plotly.graph_objects as go
import polars as pl
from IPython.display import display
from ml4t.diagnostic.evaluation.stats import benjamini_hochberg_fdr
from ml4t.diagnostic.metrics import (
    compute_ic_hac_stats,
    compute_ic_uncertainty,
    cross_sectional_ic_series,
)
from plotly.subplots import make_subplots

from case_studies.utils.feature_engineering import quantile_profile
from utils.artifact_specs import load_setup_config, resolve_label_buffer
from utils.cv_splits import generate_cv_splits
from utils.data_quality import top_entities, validate_modeling_inputs
from utils.paths import display_path, get_case_study_dir
from utils.style import COLORS, GRAY_FILLS, show_plotly_with_alt
```

```python
MAX_SYMBOLS = 0
```

## Configuration

The label, the horizon it looks forward over and the date the holdout begins are read
from `config/setup.yaml`, which is also what the training notebooks read, so the two
cannot drift apart on which return is being predicted. The thresholds below are
choices this notebook makes. Each is bound once, and the value is not repeated in
prose or in a figure title, because a threshold written twice is a decision with two
sources of truth. What each one decides:

- **Coverage** and **staleness** are the two correctness gates. A feature must have a
  value on most of the rows it is screened on, and must not simply repeat the previous
  session's number on most of them. A feature that is mostly absent cannot be ranked
  on, and one that rarely moves cannot separate this week's products from last week's.
- **The false-discovery level** is the share of the features this notebook calls
  significant that we are willing to have be noise.
- **The exploration bar** is a pair of conditions. A feature that is not significant on
  its own is still carried forward if it points the same way in most of the validation
  years and its average correlation clears a floor. A rank correlation is unitless and
  does not convert into a return, so the floor is a judgement about how weak an
  association can be and still be worth carrying into a strategy that pays commission,
  bid-ask spread and roll slippage. It is a stated choice rather than an estimate, and
  Section 6 says what follows from that.
- **The redundancy cut** is the correlation above which two features are reported as
  one piece of evidence rather than two.
- **The minimum cross-section** is how many products must carry both a feature value
  and a label on a date before that date's correlation is used at all. A rank
  correlation over a handful of products is mostly noise. It is capped by the universe
  actually loaded, so a run over a reduced universe narrows the gate with it rather
  than silently discarding every date.
- **The extreme-value limits** are what the quality gate treats as impossible rather
  than merely large: a weekly futures return above the return limit is a failure of
  the back-adjustment that stitches consecutive contracts into one price series, not a
  market move.

```python
CASE_STUDY_ID = "cme_futures"
CASE_DIR = get_case_study_dir(CASE_STUDY_ID)
EVAL_DIR = CASE_DIR / "evaluation"
EVAL_DIR.mkdir(exist_ok=True)

SETUP = load_setup_config(CASE_STUDY_ID)
eval_config = SETUP["evaluation"]
SECTORS = {
    product: sector
    for sector, products in SETUP["universe"]["product_groups"].items()
    for product in products
}

PRIMARY_LABEL = SETUP["labels"]["primary"]
LABEL_BUFFER = resolve_label_buffer(CASE_STUDY_ID, PRIMARY_LABEL, SETUP)
assert LABEL_BUFFER, f"No label buffer configured for {PRIMARY_LABEL}"


def horizon_sessions(buffer_spec: str) -> int:
    """Sessions a label looks forward over, from its `setup.yaml` buffer."""
    return int(re.match(r"^(\d+)", buffer_spec).group(1))


PRIMARY_HORIZON = horizon_sessions(LABEL_BUFFER)
# The bandwidth of the serial-correlation correction is the label horizon, because that
# is the overlap a daily-sampled forward return induces. Deriving it here rather than
# typing it keeps the correction and the label from drifting apart.
HAC_MAXLAGS = PRIMARY_HORIZON
# Every forward return the case study declares, primary first, for the horizon profile.
LABEL_HORIZONS = {PRIMARY_LABEL: PRIMARY_HORIZON} | {
    name: horizon_sessions(spec) for name, spec in SETUP["labels"]["variant_buffers"].items()
}

FDR_ALPHA = 0.05  # Benjamini-Hochberg level
NAIVE_T = 1.96  # two-sided normal critical value, both significance tiers are read against
REDUNDANCY_CUT = 0.7  # |rho| above which two features are reported as one piece of evidence
MIN_SIGN_CONSISTENCY = 0.60  # fold-sign agreement the exploration arm requires
IC_THRESHOLD = 0.008  # |IC| the exploration arm requires, at a weekly horizon
N_QUANTILES = 5
MIN_COVERAGE = 0.70  # non-null fraction the correctness gate requires
MAX_STALENESS = 0.50  # unchanged-from-prior-date fraction the correctness gate allows
MAX_ABS_RETURN = 1.0  # a weekly futures return past this is a back-adjustment failure
MAX_ABS_FEATURE = 1e6  # feature magnitude past this is a construction failure, not a value
CROSS_SECTION_TARGET = 10  # products a date needs before its rank correlation is used
MIN_IC_DATES = 20  # dates a feature needs before its IC series is summarized at all
MIN_FOLD_DATES = 5  # dates a fold must contribute before its mean IC is read
LEADING_FOR_SERIES = 3  # features the IC-through-time figure draws, by |IC|
ROLLING_DAYS = 126  # sessions in the rolling mean drawn over the IC series, about half a year
RANKED_SHOWN = 25  # features the ranking figure draws
FOLDS_SHOWN = 12  # features the fold-stability figure draws
HORIZON_SHOWN = 8  # features the horizon profile draws
SHAPE_SHOWN = 6  # features the quantile-profile figure draws
TOP_PAIRS = 15  # correlated pairs the redundancy figure ranks
```

## The data this screen runs on

Three upstream artifacts are combined into one table with a row per product per
session: the price-derived features from `features/financial.parquet`, the
model-derived features from `features/model_based.parquet`, and the forward return
from `labels/`. Every candidate feature is a column, and the label is the column each
one is measured against.

**One contract per product.** A futures product does not trade as a single
instrument. On any day the exchange lists several delivery months of the same
underlying, and the pipeline tracks three of them: the *front month*, the contract
nearest to delivery and the most heavily traded, and the two behind it, the first and
second *deferred*. This screen keeps the front month alone. The three contracts on one
product are three views of the same market, so keeping all of them would present the
rank correlation with roughly ninety entities per date when there are only thirty
independent ones, and a correlation across near-duplicate rows looks more reliable
than it is.

**The holdout stays unread.** Everything computed below is an input to which features
the later notebooks work with, which is a form of selection. `config/setup.yaml` names
the date the holdout starts; no row on or after it is loaded into any statistic here,
so the holdout is still untouched when it is used, once, to confirm or disconfirm the
final result.

**Which rows the model-derived features are read on.** `model_based.parquet` carries one
row per `(timestamp, product, position)` and no `fold` column. Both fitted families in it
are re-estimated on the schedule `setup.yaml` declares, so a value on a date was produced
by an estimate made from sessions strictly earlier than it -
`04_model_based_features` section A sets out why a walk-forward period does not do that
job on its own - and the same value is out of sample wherever it is read. There is no
vintage to choose between and nothing to select by fold id.

The screen still runs on the union of the five validation windows rather than on the
whole span before the holdout, and the price-derived features are measured on the same
rows. A correlation measured over the whole span is not comparable with one measured over
the validation windows alone, so a ranking that mixed the two would rank the window as
much as the feature - and the training sessions are the ones the estimates behind a
column read, which is the half a screen must not report on.

```python
features = pl.read_parquet(CASE_DIR / "features" / "financial.parquet").filter(
    pl.col("position") == 0
)
temporal = pl.read_parquet(CASE_DIR / "features" / "model_based.parquet").filter(
    pl.col("position") == 0
)
label_df = pl.read_parquet(CASE_DIR / "labels" / f"{PRIMARY_LABEL}.parquet").filter(
    pl.col("position") == 0
)
label_col = [c for c in label_df.columns if c not in ("timestamp", "product", "position")][0]

HOLDOUT_START = date(*map(int, eval_config["holdout_start"].split("-")))

# position is fixed at 0 throughout, and kept in the key so every join below is on the
# same three columns the artifacts are written with.
JOIN_COLS = ["timestamp", "product", "position"]
DATE_COL = "timestamp"
```

### The universe

The thirty products the screen ranks against each other, grouped the way the exchange
and the trading desk group them. This is what "the cross-section" means in every
statistic below: on a given session, each product contributes one feature value and
one forward return, and the correlation is taken across the rows of this table.

Two things in it matter for what follows. The sectors differ in how many products they
contain, so a sector with more products moves the cross-sectional ranking more than
one with fewer, whatever the economics. And the first session differs across products,
so the early years of the panel rank fewer products than the late ones.

```python
universe = (
    label_df.with_columns(pl.col("product").replace_strict(SECTORS).alias("sector"))
    .group_by("sector", "product")
    .agg(pl.col(DATE_COL).min().alias("start"))
    .group_by("sector")
    .agg(
        pl.len().alias("n_products"),
        pl.col("product").sort().str.join(" ").alias("products"),
        pl.col("start").min().alias("earliest_start"),
        pl.col("start").max().alias("latest_start"),
    )
    .sort("sector")
)
with pl.Config(tbl_rows=universe.height, tbl_width_chars=170, fmt_str_lengths=40):
    display(universe)
```

### The walk-forward folds

*Walk-forward* means the sample is cut into consecutive blocks, each with an earlier
training window and a later validation window that the training window never reaches
into, so a feature is always judged on dates after the ones any estimator saw. The
boundaries come from `generate_cv_splits`, which reads them out of `config/setup.yaml`
and is the same call every other notebook in this pipeline makes. They are derived
here rather than copied, so the periods printed below are the ones
`04_model_based_features` screens over too.

The gap between a training window's end and its validation window's start is the label
horizon. Without it, the last training dates would carry a return that resolves after
validation begins, and the estimator would have seen part of what it is about to be
tested on. The same reasoning ends each validation window short of the holdout: a
decision made on the last validation date has to have its return observable before the
holdout period starts, so the last few sessions of the final fold are dropped.

```python
splits = generate_cv_splits(
    label_df.select(DATE_COL).unique().sort(DATE_COL),
    case_study_id=CASE_STUDY_ID,
    label_buffer=LABEL_BUFFER,
)


def _as_date(value) -> date:
    return pd.Timestamp(value).date()


display(
    pl.DataFrame(
        [
            {
                "fold": split["fold"],
                "train_start": _as_date(split["train_start"]),
                "train_end": _as_date(split["train_end"]),
                "val_start": _as_date(split["val_start"]),
                "val_end": _as_date(split["val_end"]),
            }
            for split in splits
        ]
    ).sort("fold")
)
```

### Checking the artifact is keyed the way this notebook believes it is

Under the previous artifact this section had real work to do: the hidden Markov model was
fitted once per fold, so its value for a date depended on which fold was asking, and each
fitted column had to be taken from its own fold's validation window while the unfitted
ones were asserted identical across folds and read once.

None of that survives the conversion. Every fitted column is now produced by a refit
schedule, so it is out of sample wherever it is read, and the artifact holds one row per
key. What is left is the check that this is actually true of the file in front of us: a
`fold` column here, or a duplicated key, means the writer changed and this cell did not.

```python
assert "fold" not in temporal.columns, (
    "model_based.parquet carries a fold column; this notebook reads the fold-free artifact "
    "04_model_based_features writes and would count a session once per fold"
)
assert temporal.select(JOIN_COLS).is_duplicated().sum() == 0, (
    "model_based.parquet is not one row per (timestamp, product, position)"
)
# The artifact as written. The quality gate below reads this rather than the screened
# frame, because the training notebooks read every row of it and the screened frame keeps
# only the validation windows.
temporal_artifact = temporal
temporal_feature_cols = [c for c in temporal.columns if c not in JOIN_COLS]
val_windows = {
    int(sp["fold"]): (_as_date(sp["val_start"]), _as_date(sp["val_end"])) for sp in splits
}
IN_VALIDATION = pl.any_horizontal(
    [(pl.col(DATE_COL) >= start) & (pl.col(DATE_COL) <= end) for start, end in val_windows.values()]
)
print(
    f"Model-based artifact: {len(temporal):,} rows on {JOIN_COLS}, "
    f"{temporal[DATE_COL].min()} to {temporal[DATE_COL].max()}, no fold column."
)
```

### Assembling the table

`ls_signal` and `risk_adj_score` are dropped: they are the case study's own composite
trading signals, built downstream from several of the columns beside them, so scoring
them here would measure the combination rather than a candidate input to it.

```python
COMPOSITE_COLS = ["ls_signal", "risk_adj_score"]
METADATA_COLS = {*JOIN_COLS, *COMPOSITE_COLS}
financial_cols = [c for c in features.columns if c not in METADATA_COLS]
temporal_cols = [c for c in temporal.columns if c not in JOIN_COLS]
all_feature_cols = financial_cols + temporal_cols

eval_panel = (
    features.drop(COMPOSITE_COLS, strict=False)
    .join(temporal, on=JOIN_COLS, how="left")
    .join(label_df, on=JOIN_COLS, how="inner")
    # The holdout first, then the validation windows, which is what every screen reads.
    .filter(pl.col(DATE_COL) < HOLDOUT_START)
    .filter(IN_VALIDATION)
)
assert eval_panel[DATE_COL].max() < HOLDOUT_START

# Reduce the universe for a fast development run. Left at zero, every product is kept.
# `top_entities` breaks a tie on the product code, which a local sort did not: on a padded
# grid every product quoting the whole window carries the same row count, and the winner
# then came from frame order. Two callers reducing the same panel to the same size have to
# choose the same universe, or a product one of them kept joins to null features in the
# other and the run answers cleanly and wrongly.
if MAX_SYMBOLS > 0:
    kept = top_entities(eval_panel, MAX_SYMBOLS, entity_col="product")
    eval_panel = eval_panel.filter(pl.col("product").is_in(kept))

n_rows = len(eval_panel)
n_symbols = eval_panel["product"].n_unique()
n_dates = eval_panel[DATE_COL].n_unique()
# A rank correlation over a handful of products is mostly noise, so a date has to carry
# a minimum cross-section to count. Capping it by the universe actually loaded means a
# reduced run measures something rather than discarding every date.
MIN_CROSS_SECTION = min(CROSS_SECTION_TARGET, n_symbols)

print(f"Screening {len(all_feature_cols)} features against {label_col}")
print(
    f"  {n_rows:,} product-sessions, {n_symbols} products, {n_dates} sessions, "
    f"{eval_panel[DATE_COL].min()} to {eval_panel[DATE_COL].max()}"
)
print(f"  {len(financial_cols)} price-derived, {len(temporal_cols)} model-derived")
print(f"  a session counts once at least {MIN_CROSS_SECTION} products carry both columns")
```

### What each fold's training window can actually supply

A feature can be present in the declared set, present in every model's recorded
specification, and still be almost entirely absent from the rows a given fold trains on.
Anything with a long warm-up does this by construction: a rolling statistic over a year
needs a year before it produces anything, and a feature derived from a model fitted on
history cannot exist before that model has enough history to fit.

Nothing downstream will complain. The imputer fills a missing value with the training
median, the scaler standardises that median to zero, and the fit proceeds with a column
that is present, inert, and indistinguishable in the registry from a column carrying real
variation. The only case that raises is a feature missing from *every* training row, which
is a narrow escape rather than a safety net: it fires at zero coverage and says nothing at
one per cent.

So the rate is reported here, before any model is fitted, for both halves of the declared
set. Both are read from the artifacts as written rather than from the screened frame
above, which keeps only the validation windows - and the training rows are exactly what
this question is about. For the model-derived columns the rate now varies across periods
only because the periods have different training windows, not because a different fit
produced each one.

**This table does not repair anything, and it is not meant to.** Whether a fold should be
fitted at all when one of its declared features barely exists there is a question about
what this case study is teaching, and it is answered by a reader looking at the numbers,
not by a notebook quietly imputing its way past them.

```python
train_windows = {
    int(split["fold"]): (_as_date(split["train_start"]), _as_date(split["train_end"]))
    for split in splits
}


def _coverage_rows(
    frame: pl.DataFrame, columns: list[str], source: str, fold_id: int
) -> list[dict]:
    height = frame.height
    return [
        {
            "fold": fold_id,
            "source": source,
            "feature": column,
            "train_rows": height,
            "non_null_pct": (frame.get_column(column).drop_nulls().len() / height * 100)
            if height
            else 0.0,
        }
        for column in columns
    ]


fold_coverage_rows: list[dict] = []
for fold_id, (start, end) in sorted(train_windows.items()):
    in_window = pl.col(DATE_COL).is_between(start, end)
    fold_coverage_rows += _coverage_rows(
        features.filter(in_window), financial_cols, "price-derived", fold_id
    )
    fold_coverage_rows += _coverage_rows(
        temporal_artifact.filter(in_window),
        temporal_feature_cols,
        "model-derived",
        fold_id,
    )
fold_coverage = pl.DataFrame(fold_coverage_rows)

SPARSE_PCT = 5.0
sparse = fold_coverage.filter(pl.col("non_null_pct") < SPARSE_PCT).sort("fold", "feature")
print(f"{fold_coverage.height} (fold, feature) pairs across {len(train_windows)} folds")
print(f"{sparse.height} are below {SPARSE_PCT:.0f}% non-null in their fold's training window")
```

```python
# Every feature whose fold sees almost none of it, and the per-fold rate for the columns
# that vary most across folds. A feature absent here is imputed, not refused.
fold_coverage.filter(
    pl.col("feature").is_in(sparse.get_column("feature").unique().to_list())
).pivot(on="fold", index=["source", "feature"], values="non_null_pct").sort("source", "feature")
```

## 0. Data quality gate

Before asking whether a feature predicts anything, check that the numbers are numbers:
no infinities, no column that is entirely absent, no return so large it can only be a
construction failure. The failures this catches are silent ones. A division by zero in
a derived feature or a break in the back-adjustment that stitches consecutive delivery
months into one continuous price series does not raise anything here; it surfaces
several notebooks later as an unexplained loss or a model that refuses to fit.

**This gate reads more rows than the screens below do, deliberately.** The screens are
a selection decision and run on the validation windows, where every candidate exists.
This gate asks whether the artifact is sound at all, and the rows the training
notebooks read include every fold's training window back to the start of the panel, so
a broken value there reaches a model whether or not this notebook screened it. It
therefore runs on the whole span before the holdout. It stops there for two reasons:
its counts are printed, so reading holdout rows would put a description of the holdout
into this notebook's output, and it is allowed to halt the run, which would make
whether this notebook executes depend on data it is not permitted to see.

```python
sealed_features = features.filter(pl.col(DATE_COL) < HOLDOUT_START)
# The fold-bearing artifact rather than the resolved frame, for the reason in the
# markdown above: every fold's value is a value a training notebook can read.
sealed_temporal = temporal_artifact.filter(pl.col(DATE_COL) < HOLDOUT_START)
sealed_labels = label_df.filter(pl.col(DATE_COL) < HOLDOUT_START)

print("Price-derived features and the label:")
quality_result = validate_modeling_inputs(
    features_df=sealed_features,
    label_df=sealed_labels,
    feature_cols=[c for c in features.columns if c not in set(JOIN_COLS)],
    label_col=label_col,
    join_cols=JOIN_COLS,
    asset_col="product",
    max_abs_return=MAX_ABS_RETURN,
    max_abs_feature=MAX_ABS_FEATURE,
    fail_on_critical=True,
)

print("Model-derived features, every session before the holdout, and the label:")
temporal_quality = validate_modeling_inputs(
    features_df=sealed_temporal,
    label_df=sealed_labels,
    feature_cols=temporal_cols,
    label_col=label_col,
    join_cols=JOIN_COLS,
    asset_col="product",
    max_abs_return=MAX_ABS_RETURN,
    max_abs_feature=MAX_ABS_FEATURE,
    fail_on_critical=True,
)
```

## 1. Can the value be trusted at decision time?

Two properties are checked per feature, and both are about whether the column can
carry a ranking at all rather than about whether it predicts anything.

**Coverage** is the share of rows where the feature has a value. A product with no
value on a date cannot be placed in that date's ranking, so a sparse feature is
ranking a shifting subset of the universe and its correlation is not comparable with
one measured on the full cross-section.

**Staleness** is the share of rows where the value repeats the same product's previous
session. A feature that rarely changes ranks this week's products almost exactly as it
ranked last week's, so whatever it appears to predict it predicts by persistence.
Slow-moving is not the same as broken - a quarterly quantity is legitimately flat
between releases - but at a weekly decision cadence a column that almost never moves
cannot be the thing that distinguishes one week from the next.

A feature failing either is recorded as STOP and takes no further part.

The book names four correctness checks at this point in the pipeline. These are two of
them. The other two - whether a feature's timestamp is the moment its information was
available, and whether the feature and the label are restricted to the same eligible
rows - are settled where the columns are built, in `03_financial_features` and
`04_model_based_features`, and cannot be re-derived from the artifacts they write.

Both are measured on the validation windows, so a feature fitted per fold is judged
over the rows where it is out of sample rather than over a span where it does not
exist.

```python
coverage = {}
staleness = {}
for feat in all_feature_cols:
    coverage[feat] = eval_panel[feat].drop_nulls().len() / n_rows
    ordered = eval_panel.select(["product", DATE_COL, feat]).sort(["product", DATE_COL])
    unchanged = ordered.with_columns(
        (pl.col(feat) == pl.col(feat).shift(1).over("product")).alias("unchanged")
    )["unchanged"].sum()
    # One row per product has no predecessor, so those rows are out of the denominator.
    staleness[feat] = float(unchanged) / max(n_rows - n_symbols, 1)

correctness = {
    feat: coverage[feat] >= MIN_COVERAGE and staleness[feat] <= MAX_STALENESS
    for feat in all_feature_cols
}

n_pass = sum(correctness.values())
n_fail = len(correctness) - n_pass
print(f"Correctness gate: {n_pass} pass, {n_fail} fail")

if n_fail > 0:
    fail_df = (
        pl.DataFrame(
            {
                "feature": [f for f, ok in correctness.items() if not ok],
                "coverage": [coverage[f] for f, ok in correctness.items() if not ok],
                "staleness": [staleness[f] for f, ok in correctness.items() if not ok],
            }
        )
        .with_columns(
            pl.when(pl.col("coverage") < MIN_COVERAGE)
            .then(pl.lit("coverage"))
            .otherwise(pl.lit("staleness"))
            .alias("failed_on")
        )
        .sort(["failed_on", "feature"])
    )
    with pl.Config(tbl_rows=fail_df.height, fmt_str_lengths=40):
        display(fail_df)
```

## 2. Does the feature sort next week's returns?

The measure is the **information coefficient**, or IC: on each session, rank the
products by the feature, rank the same products by the return they went on to earn,
and take the correlation between the two rankings. That gives one number per session,
and the series of those numbers is the object this section summarizes. Ranks rather
than raw values, because what a long-short book acts on is the order of the products,
not the distance between them, and because a single outsized futures return would
otherwise dominate the statistic.

The series is computed by `cross_sectional_ic_series`, which pairs the feature and the
label on both the session and the product, so a session's correlation is taken across
that session's cross-section and nothing else. A session where too few products carry
both columns is dropped, as is one where the feature is the same for every product,
because a correlation needs the feature to vary across the things being ranked.

### Why the average of that series needs a wider error bar than it looks like it needs

The label looks forward over several sessions, and it is measured on every session, so
two consecutive observations are computed from price windows that overlap almost
entirely. They are close to the same measurement made twice. A standard error that
treats them as independent therefore counts far more information than the sample
holds. Newey-West standard errors correct for exactly this, by estimating how much the
series correlates with itself at each lag up to a chosen bandwidth and inflating the
error bar accordingly:

$$N_{\text{eff}} \approx \frac{N}{1 + 2\sum_{k=1}^{q} w_k \hat{\rho}_k}$$

where $w_k$ are Bartlett kernel weights and $\hat{\rho}_k$ are the estimated
autocorrelations. With a label spanning $h$ overlapping sessions the effective sample
is roughly $N / h$, so an uncorrected standard error understates the true one by about
$\sqrt{h}$. The bandwidth is the label horizon, taken from the same configuration
entry that decides the label, so the two cannot drift apart.

The three counts printed below are a different quantity from the error bars: they say
how many features clear a fixed t-statistic under each correction. They fall as the
correction gets stricter, but not in proportion to it.

```python
evaluable_features = [f for f in all_feature_cols if correctness[f]]

ic_results = {}
ic_timeseries = {}
for feat in evaluable_features:
    valid = eval_panel.select([DATE_COL, "product", feat, label_col]).drop_nulls()
    ic_df = (
        cross_sectional_ic_series(
            valid.select([DATE_COL, "product", feat]),
            valid.select([DATE_COL, "product", label_col]),
            pred_col=feat,
            ret_col=label_col,
            date_col=DATE_COL,
            entity_col="product",
            method="spearman",
            min_obs=MIN_CROSS_SECTION,
        )
        # A thin session comes back with a null IC; a session where the feature is the
        # same for every product comes back NaN, which drop_nulls does not remove.
        .drop_nulls("ic")
        .filter(pl.col("ic").is_not_nan())
        .sort(DATE_COL)
    )
    if len(ic_df) >= MIN_IC_DATES:
        # The estimator reads the series in stored order, so the sort above is the
        # caller's responsibility and is kept visible at the call site.
        ic_results[feat] = compute_ic_hac_stats(ic_df, ic_col="ic", maxlags=HAC_MAXLAGS)
        ic_timeseries[feat] = ic_df

print(f"Measured an IC series for {len(ic_results)} of {len(evaluable_features)} features")
```

### Is it the same in every validation year?

A feature whose IC is positive in one fold and negative in the next has an average
that describes no period in the sample. The measure of that is **sign consistency**:
the share of validation windows whose mean IC points the same way as the feature's own
overall mean.

Agreement with its *own* direction is the point. A feature that is negative in every
window is exactly as stable as one that is positive in every window, and just as
usable - a long-short book takes the short side of it. The strongest associations in
this panel are negative, so a rule that scored agreement with a positive sign would
make them unpromotable at any threshold.

```python
def fold_mean_ics(ts: pl.DataFrame) -> list[float]:
    """Mean of an IC series inside each fold's validation window."""
    out = []
    for split in splits:
        window = ts.filter(
            (pl.col(DATE_COL) >= _as_date(split["val_start"]))
            & (pl.col(DATE_COL) <= _as_date(split["val_end"]))
        )
        if len(window) >= MIN_FOLD_DATES:
            out.append(window["ic"].mean())
    return out


fold_stats = {}
for feat in ic_results:
    fold_ics = fold_mean_ics(ic_timeseries[feat])
    if fold_ics:
        pooled_sign = np.sign(ic_results[feat]["mean_ic"])
        agreeing = sum(1 for ic in fold_ics if np.sign(ic) == pooled_sign and pooled_sign != 0)
        fold_stats[feat] = {
            "n_folds": len(fold_ics),
            "sign_consistency": agreeing / len(fold_ics),
            "worst_fold_ic": min(fold_ics),
            "best_fold_ic": max(fold_ics),
            "median_fold_ic": float(np.median(fold_ics)),
        }

n_consistent = sum(1 for s in fold_stats.values() if s["sign_consistency"] >= MIN_SIGN_CONSISTENCY)
print(
    f"{n_consistent} features keep their own sign in at least "
    f"{MIN_SIGN_CONSISTENCY:.0%} of the validation windows"
)
```

## 3. Correcting for having tested every feature at once

A significance level of five percent means one test in twenty clears it on noise. Run
dozens of tests and a handful of apparent discoveries is the expected outcome even if
no feature carries anything. The Benjamini-Hochberg procedure controls the **false
discovery rate**: of the features this notebook announces as significant, the expected
share that are noise is held to the configured level. It does that by ranking the
p-values and comparing each against a threshold that loosens with rank, which is less
severe than demanding that no false positive occur at all and is the right trade when
the output is a shortlist to investigate rather than a single claim.

**What was searched has to be stated, or the corrected p-value means nothing.** The
search here is every column the two upstream notebooks wrote, screened once, against
one label. The price-derived columns are the term-structure, carry, momentum,
volatility, technical and calendar families built in `03_financial_features`; the
model-derived ones are the three families in `04_model_based_features`. Neither
notebook generated its columns by scanning parameter grids for association with this
label, so the candidate set is not itself the product of a search - which is what makes
the count below the right denominator. The search that this correction cannot see is
the one across horizons, and Section 3.4 is where that becomes visible.

Three tiers are printed. **Uncorrected** ignores both problems. **Newey-West** allows
for the overlapping label but not for the number of tests. **False-discovery** allows
for both.

```python
feature_names = list(ic_results.keys())
p_values = [ic_results[f]["p_value"] for f in feature_names]

fdr_result = benjamini_hochberg_fdr(p_values, alpha=FDR_ALPHA, return_details=True)

# Build evaluation summary
eval_summary = pl.DataFrame(
    {
        "feature": feature_names,
        "source": ["temporal" if f in temporal_cols else "financial" for f in feature_names],
        "ic_mean": [ic_results[f]["mean_ic"] for f in feature_names],
        "hac_se": [ic_results[f]["hac_se"] for f in feature_names],
        "hac_t": [ic_results[f]["t_stat"] for f in feature_names],
        "hac_p": p_values,
        "fdr_p": list(fdr_result["adjusted_p_values"]),
        "fdr_sig": list(fdr_result["rejected"]),
        "naive_t": [ic_results[f]["naive_t_stat"] for f in feature_names],
    },
    schema_overrides={
        "ic_mean": pl.Float64,
        "hac_se": pl.Float64,
        "hac_t": pl.Float64,
        "hac_p": pl.Float64,
        "fdr_p": pl.Float64,
        "fdr_sig": pl.Boolean,
        "naive_t": pl.Float64,
    },
).sort(pl.col("ic_mean").cast(pl.Float64, strict=False).abs(), descending=True)

# The uncorrected tier reads the uncorrected t-statistic; `p_values` above holds the
# Newey-West one, and reading it here would make the two tiers the same test.
n_significant_naive = sum(1 for f in feature_names if abs(ic_results[f]["naive_t_stat"]) > NAIVE_T)
n_significant_hac = sum(1 for f in feature_names if abs(ic_results[f]["t_stat"]) > NAIVE_T)
n_significant_fdr = int(fdr_result["n_rejected"])


def shrinkage(numerator: int, denominator: int) -> str:
    """How many times more features the uncorrected test announces, or why there is no ratio."""
    if denominator == 0:
        return "undefined - the corrected test announces nothing"
    return f"{numerator / denominator:.2f}x"


print(f"Searched set: {len(feature_names)} features, one label")
print(f"  uncorrected, |t| > {NAIVE_T}:      {n_significant_naive}")
print(f"  Newey-West, |t| > {NAIVE_T}:       {n_significant_hac}")
print(f"  false-discovery, q < {FDR_ALPHA}:   {n_significant_fdr}")
print(f"  uncorrected / Newey-West:      {shrinkage(n_significant_naive, n_significant_hac)}")
print(f"  uncorrected / false-discovery: {shrinkage(n_significant_naive, n_significant_fdr)}")
```

### The series behind the average

Every number above is a one-line summary of the same object, the per-session IC
series, and two ways of arriving at a respectable average are invisible in the
summary. One is an association confined to a single episode, with the series flat
either side of it. The other is a series that changes sign between validation windows
and averages to something no period actually earned. So the series is drawn before it
is reduced, for the features with the largest absolute average.

The thin line is the session-by-session IC and it is mostly noise at this width; the
heavy line is its rolling mean over about half a year, which is what makes an episode
visible. The horizontal lines are the full-sample average and the Newey-West interval
around it - the interval covers the uncertainty in the average, not the spread of the
series, and the gap between the two is the point of the figure.

```python
series_features = eval_summary.head(LEADING_FOR_SERIES)["feature"].to_list()

fig = make_subplots(
    rows=len(series_features), cols=1, shared_xaxes=True, subplot_titles=series_features
)
for row, feat in enumerate(series_features, start=1):
    series = ic_timeseries[feat].sort(DATE_COL)
    bands = compute_ic_uncertainty(series, horizon=HAC_MAXLAGS, ic_col="ic")
    dates = series[DATE_COL].to_list()
    fig.add_trace(
        go.Scatter(
            x=dates,
            y=series["ic"].to_list(),
            mode="lines",
            line=dict(color=GRAY_FILLS["muted"], width=0.6),
            name="One session",
            legendgroup="session",
            showlegend=row == 1,
        ),
        row=row,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=dates,
            y=series["ic"].rolling_mean(ROLLING_DAYS, min_samples=ROLLING_DAYS).to_list(),
            mode="lines",
            line=dict(color=COLORS["blue"], width=1.4),
            name="Rolling mean",
            legendgroup="rolling",
            showlegend=row == 1,
        ),
        row=row,
        col=1,
    )
    for value, dash, label, group in (
        (bands["mean_ic"], "solid", "Full-sample average", "mean"),
        (bands["ci_hac_lower"], "dot", "Its Newey-West interval", "band"),
        (bands["ci_hac_upper"], "dot", None, None),
    ):
        fig.add_hline(y=value, line=dict(color=COLORS["amber"], width=1, dash=dash), row=row, col=1)
        if label and row == 1:
            fig.add_trace(
                go.Scatter(
                    x=[dates[0]],
                    y=[value],
                    mode="lines",
                    line=dict(color=COLORS["amber"], width=1, dash=dash),
                    name=label,
                    legendgroup=group,
                ),
                row=row,
                col=1,
            )
    fig.add_hline(y=0, line=dict(color=GRAY_FILLS["border"], width=0.8), row=row, col=1)
fig.update_yaxes(title_text="Rank correlation")
fig.update_layout(
    template="ml4t",
    height=200 * len(series_features) + 140,
    width=900,
    title_text="The daily IC swings far wider than the mean it averages to",
    legend=dict(orientation="h", y=-0.12),
)
show_plotly_with_alt(
    fig,
    "One stacked time-series panel per feature. In each, the per-session rank correlation is a "
    "pale line filling most of the vertical range, swinging between roughly plus and minus 0.5 "
    "from one session to the next. Over it sit a dark rolling mean that stays within a narrow "
    "band near zero, a flat full-sample average line, and a dashed Newey-West interval around "
    "that average. The session-by-session scatter is an order of magnitude wider than the mean "
    "it averages to, which is the comparison the figure is making.",
)
```

### Every feature ranked, with the correction attached

The bars are the average IC, ordered by size, and the number on each is its
Newey-West t-statistic. Colour is one thing throughout this section and the next: the
feature cleared false-discovery control, or it did not. Reading the two together is
the point - a bar can be among the longest and still be uncoloured, which says the
association is large relative to the others and small relative to its own uncertainty.

```python
top = eval_summary.head(min(RANKED_SHOWN, len(eval_summary))).sort("ic_mean")

SURVIVES, DOES_NOT = COLORS["blue"], GRAY_FILLS["muted"]

fig = go.Figure(
    go.Bar(
        x=top["ic_mean"].to_list(),
        y=top["feature"].to_list(),
        orientation="h",
        marker_color=[SURVIVES if sig else DOES_NOT for sig in top["fdr_sig"].to_list()],
        text=[f"t={value:.1f}" for value in top["hac_t"].to_list()],
        textposition="outside",
        showlegend=False,
    )
)
for colour, label in ((SURVIVES, "Clears false-discovery control"), (DOES_NOT, "Does not")):
    fig.add_trace(go.Bar(x=[None], y=[None], orientation="h", marker_color=colour, name=label))
fig.add_vline(x=0, line=dict(color=GRAY_FILLS["border"], width=1))
# The t-statistic is drawn past the end of its bar, so the axis needs room beyond the
# longest one or the leftmost label runs into the feature names.
ic_span = max(abs(value) for value in top["ic_mean"].to_list()) * 1.35
fig.update_layout(
    template="ml4t",
    height=620,
    width=900,
    title_text="Almost nothing here survives false-discovery control",
    xaxis_title="Mean cross-sectional rank correlation with the forward return",
    xaxis_range=[-ic_span, ic_span],
    yaxis_title="Feature",
    margin=dict(l=170),
    legend=dict(orientation="h", y=-0.14),
)
show_plotly_with_alt(
    fig,
    "Horizontal bar chart of the mean cross-sectional rank correlation with the forward return, "
    "one bar per feature, sorted from the most positive at the top through zero to the most "
    "negative at the bottom, with each bar annotated by its t-statistic. The volatility family "
    "occupies the positive end and the carry-momentum family the negative end, but every bar is "
    "small, within a few hundredths of zero. Bars are shaded by whether they clear "
    "false-discovery control, and almost the whole chart is in the colour that does not.",
)
```

### The same features, one window at a time

The average above is drawn over the whole span. This is the same quantity computed
inside each validation window separately: one dot per window, with the middle one
marked. A feature whose dots sit on one side of zero is one the average describes; a
feature whose dots straddle zero has an average that is a compromise between periods
that disagreed. The dots and the marker are produced by the same function the
sign-consistency screen reads, so the marked value is always the middle of the windows
drawn beside it.

```python
fold_features = [f for f in eval_summary["feature"].to_list() if f in fold_stats][:FOLDS_SHOWN]
fig = go.Figure()
for feat in fold_features:
    per_fold = fold_mean_ics(ic_timeseries[feat])
    fig.add_trace(
        go.Scatter(
            x=per_fold,
            y=[feat] * len(per_fold),
            mode="markers",
            marker=dict(color=GRAY_FILLS["muted"], size=7),
            showlegend=False,
        )
    )
    fig.add_trace(
        go.Scatter(
            x=[fold_stats[feat]["median_fold_ic"]],
            y=[feat],
            mode="markers",
            marker=dict(color=COLORS["amber"], size=11, symbol="diamond"),
            showlegend=False,
        )
    )
fig.add_vline(x=0, line=dict(color=GRAY_FILLS["border"], width=1))
fig.update_layout(
    template="ml4t",
    height=520,
    width=900,
    title_text="Most leading features change sign in one window, and a couple never do",
    xaxis_title="Mean rank correlation within one validation window (diamond marks the middle one)",
    margin=dict(l=170),
)
show_plotly_with_alt(
    fig,
    "Dot plot with one row per leading feature. Each row shows the feature's mean rank "
    "correlation in each validation window as a small grey dot, with an amber diamond marking "
    "the middle window. A vertical rule marks zero. Most rows have dots on both sides of that "
    "rule, so the feature's direction reverses from one window to the next; only a couple of "
    "rows keep every dot on one side. Within a row the dots are spread across much of the "
    "chart's width, so the disagreement between windows is the dominant feature of the "
    "picture rather than the position of any one of them.",
)
```

### What the overlap correction costs

One point per feature: its t-statistic before the correction against its t-statistic
after. The dashed line is where the two agree. A point below it in the upper half, or
above it in the lower half, is a feature whose evidence was overstated by treating
overlapping observations as independent.

```python
fig = go.Figure(
    go.Scatter(
        x=eval_summary["naive_t"].to_list(),
        y=eval_summary["hac_t"].to_list(),
        mode="markers",
        marker=dict(
            color=[SURVIVES if s else DOES_NOT for s in eval_summary["fdr_sig"].to_list()],
            size=7,
        ),
        text=eval_summary["feature"].to_list(),
        showlegend=False,
    )
)
max_t = (
    max(
        eval_summary["naive_t"].cast(pl.Float64, strict=False).abs().max() or 1.0,
        eval_summary["hac_t"].cast(pl.Float64, strict=False).abs().max() or 1.0,
    )
    * 1.1
)
fig.add_trace(
    go.Scatter(
        x=[-max_t, max_t],
        y=[-max_t, max_t],
        mode="lines",
        line=dict(dash="dash", color=GRAY_FILLS["border"]),
        showlegend=False,
    )
)
fig.update_layout(
    template="ml4t",
    height=480,
    width=760,
    title_text="Overlapping labels pull every t-statistic toward zero",
    xaxis_title="t-statistic treating each session as independent",
    yaxis_title="t-statistic after the Newey-West correction",
)
show_plotly_with_alt(
    fig,
    "Scatter plot of each feature's t-statistic after the Newey-West correction against the "
    "same t-statistic computed as if sessions were independent, with a dashed 45-degree "
    "reference line. Every point sits between that line and the horizontal axis rather than on "
    "it, so each corrected value is smaller in magnitude than the uncorrected one, and the gap "
    "widens with distance from the origin. Points are shaded by whether the feature clears "
    "false-discovery control, and almost all of them are in the colour that does not.",
)
```

The correction moves points toward the horizontal axis, which is what a positive
autocorrelation in the IC series does: the same evidence, counted once instead of
several times. False-discovery control then applies a second penalty for the number of
features tested, and the features drawn in the survivor colour are the ones that clear
both.

### The same features against the other forward horizon

Everything above reads one label. The case study also builds a longer forward return,
and the difference between the two matters for what a strategy could do with a
feature: an association that exists at the weekly horizon and has gone by the monthly
one cannot be held for a month, whatever its average says. Repeating the measurement
at each declared horizon is also what makes the size of the search visible - screening
two horizons and reporting the better one is a search the false-discovery correction
above cannot see, because it was given one set of p-values.

Each horizon is held out of the holdout on its own terms. A longer forward return
settles later, so a session safely inside the development period for the weekly label
can have its monthly label realized inside the holdout. The sessions within one
horizon's reach of the holdout are therefore dropped for that horizon separately.

The count runs over the sessions each *product* trades, not over the sessions the market
trades. The products here do not quote the same calendar: the equity-index contract
quotes 2,170 sessions in this window and the grains over 3,750, so twenty-one of the
market's sessions can be a good deal fewer than twenty-one of a given product's, and a
row sealed on the market's calendar can still have its outcome fall inside the holdout.
Counting within the product is what makes the seal mean the same thing for every one.

Two quantities per feature per horizon. The **average IC** is the level. The **IC
information ratio** is that average divided by how much it varies across the
validation windows, which is a per-unit-of-instability reading of the same thing. The
ratio is computed across windows rather than across sessions on purpose: at the longer
horizon consecutive sessions overlap far more, so a session-level ratio would flatter
the longer horizon for a reason that has nothing to do with the feature.

```python
leaders = eval_summary.head(HORIZON_SHOWN)["feature"].to_list()
horizon_rows = []
for variant_label, variant_horizon in LABEL_HORIZONS.items():
    variant = (
        pl.read_parquet(CASE_DIR / "labels" / f"{variant_label}.parquet")
        .filter(pl.col("position") == 0)
        .select([*JOIN_COLS, variant_label])
    )
    # A label advances over the sessions its own product trades, and the products do not
    # trade the same sessions: RTY quotes 2,170 of the 3,853 in this window while the
    # grains quote over 3,750. Counting the horizon on the market's calendar therefore
    # lands short for a product that skips sessions, and the row it leaves in has its
    # outcome realized inside the holdout. So each row's settling session is found within
    # its own product, and a row is kept only when that session falls before the holdout.
    sealed = (
        variant.sort(["product", DATE_COL])
        .with_columns(pl.col(DATE_COL).shift(-variant_horizon).over("product").alias("_settles"))
        .filter(pl.col("_settles").is_not_null() & (pl.col("_settles") < HOLDOUT_START))
        .drop("_settles")
    )
    joined = eval_panel.select([*JOIN_COLS, *leaders]).join(sealed, on=JOIN_COLS, how="inner")
    assert joined.is_empty() or joined[DATE_COL].max() < HOLDOUT_START
    for feat in leaders:
        valid = joined.select([DATE_COL, "product", feat, variant_label]).drop_nulls()
        series = (
            cross_sectional_ic_series(
                valid.select([DATE_COL, "product", feat]),
                valid.select([DATE_COL, "product", variant_label]),
                pred_col=feat,
                ret_col=variant_label,
                date_col=DATE_COL,
                entity_col="product",
                method="spearman",
                min_obs=MIN_CROSS_SECTION,
            )
            .drop_nulls("ic")
            .filter(pl.col("ic").is_not_nan())
            .sort(DATE_COL)
        )
        if len(series) < MIN_IC_DATES:
            continue
        window_means = fold_mean_ics(series)
        dispersion = float(np.std(window_means, ddof=1)) if len(window_means) > 1 else float("nan")
        horizon_rows.append(
            {
                "feature": feat,
                "horizon": variant_horizon,
                "ic_mean": float(series["ic"].mean()),
                "icir": float(np.mean(window_means)) / dispersion if dispersion else float("nan"),
            }
        )

horizon_ic = pl.DataFrame(horizon_rows)
print(
    f"Horizon profile: {len(leaders)} features across "
    f"{horizon_ic['horizon'].n_unique()} declared forward horizons"
)
```

```python
fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("Average rank correlation", "Average divided by its spread across windows"),
    horizontal_spacing=0.12,
)
shown_direction = set()
for feat in leaders:
    profile = horizon_ic.filter(pl.col("feature") == feat).sort("horizon")
    if not len(profile):
        continue
    at_primary = profile.filter(pl.col("horizon") == PRIMARY_HORIZON)
    if not len(at_primary):
        continue
    positive = float(at_primary["ic_mean"][0]) > 0
    group = "Positive at the primary horizon" if positive else "Negative at the primary horizon"
    line = dict(color=COLORS["blue"] if positive else COLORS["copper"], width=1.6)
    for column, column_index in (("ic_mean", 1), ("icir", 2)):
        fig.add_trace(
            go.Scatter(
                x=profile["horizon"].to_list(),
                y=profile[column].to_list(),
                mode="lines+markers",
                line=line,
                opacity=0.85,
                name=group,
                legendgroup=group,
                showlegend=column_index == 1 and group not in shown_direction,
                hovertext=feat,
            ),
            row=1,
            col=column_index,
        )
    shown_direction.add(group)
fig.add_hline(y=0, line=dict(color=GRAY_FILLS["border"], width=0.8))
fig.update_layout(
    template="ml4t",
    height=460,
    width=980,
    title_text="Every leading feature keeps its direction at the longer horizon",
    legend=dict(orientation="h", y=-0.22),
)
for column_index in (1, 2):
    fig.update_xaxes(
        title_text="Forward horizon (sessions)",
        tickmode="array",
        tickvals=sorted(LABEL_HORIZONS.values()),
        row=1,
        col=column_index,
    )
fig.update_yaxes(title_text="Mean rank correlation", row=1, col=1)
fig.update_yaxes(title_text="Mean divided by its spread across windows", row=1, col=2)
show_plotly_with_alt(
    fig,
    "Two side-by-side slope panels, each with one line per leading feature joining its value at "
    "the 5-session horizon to its value at the 21-session horizon. Lines are coloured by their "
    "sign at the primary horizon. In the left panel, measuring the mean rank correlation, the "
    "positive and negative groups stay cleanly on their own sides of zero at both horizons and "
    "no line crosses over. The right panel divides each mean by its spread across windows; the "
    "same separation holds there, and the negative group sits further from zero on that scale "
    "than the positive group does.",
)
```

## 4. What shape is the relationship?

A correlation says a feature and the return move together; it does not say the
relationship is one a model can use in a straight line. This section splits the
products into five buckets by the feature and reads off the average return in each,
from the lowest bucket to the highest. A profile that rises or falls all the way
across is one a single coefficient can carry. One that turns in the middle says the
information is there but a linear model cannot reach it, and a tree-based model in a
later notebook could.

**The buckets are formed inside each session and every session counts once.** The
correlation beside them ranks products against each other on a given day, so the
profile has to be built the same way. Cutting the buckets once over every row of the
panel would sort an observation from one year against one from another and let the
feature's movement over time stand in for its spread across products, which is a
different question with a different answer - and one whose profile could contradict
the correlation it sits next to, with nothing in the notebook to reconcile them.

Two summaries are drawn per bucket. The **mean** is what a book holding every product
in the bucket would earn. The **median** is the typical product. Where they disagree,
the difference is a few large returns rather than the shape of the relationship, and a
rank correlation follows the median.

The panels show the features that cleared false-discovery control first, then the next
largest by absolute correlation, so the figure stays informative when few features
clear it.

```python
fdr_shape = eval_summary.filter(pl.col("fdr_sig").fill_null(False))["feature"].to_list()
top_features_for_shape = (
    fdr_shape + [f for f in eval_summary["feature"].to_list() if f not in fdr_shape]
)[:SHAPE_SHOWN]

QUANTILE_LABELS = [f"Q{i + 1}" for i in range(N_QUANTILES)]

monotonicity_scores = {}
quantile_spreads = {}
for feat in top_features_for_shape:
    profile = quantile_profile(
        eval_panel,
        feat,
        label_col,
        date_col=DATE_COL,
        n_quantiles=N_QUANTILES,
        min_cross_section=MIN_CROSS_SECTION,
    )
    if profile is None or profile.periods_used < MIN_IC_DATES:
        continue
    quantile_spreads[feat] = {
        "q_means": profile.means,
        "q_medians": profile.medians,
        "spread": profile.spread,
    }
    # The recorded shape score is the rank correlation between the bucket's position
    # and its mean return: one at a profile that rises all the way across, minus one at
    # one that falls all the way, near zero at one that turns in the middle.
    monotonicity_scores[feat] = profile.monotonicity
```

```python
if quantile_spreads:
    n_show = min(SHAPE_SHOWN, len(quantile_spreads))
    feats_to_show = list(quantile_spreads.keys())[:n_show]
    n_rows_fig = (n_show + 2) // 3
    fig = make_subplots(rows=n_rows_fig, cols=3, subplot_titles=feats_to_show, shared_yaxes=True)
    for idx, feat in enumerate(feats_to_show):
        r, c = divmod(idx, 3)
        fig.add_trace(
            go.Bar(
                x=QUANTILE_LABELS,
                y=quantile_spreads[feat]["q_means"],
                marker_color=COLORS["blue"],
                name="mean",
                legendgroup="mean",
                showlegend=idx == 0,
            ),
            row=r + 1,
            col=c + 1,
        )
        fig.add_trace(
            go.Scatter(
                x=QUANTILE_LABELS,
                y=quantile_spreads[feat]["q_medians"],
                mode="markers",
                marker=dict(color=COLORS["amber"], size=9, symbol="diamond"),
                name="median",
                legendgroup="median",
                showlegend=idx == 0,
            ),
            row=r + 1,
            col=c + 1,
        )
    # One y range across the panels: the point of the panels side by side is that their
    # heights are comparable, which independent auto-scaling would destroy.
    span = max(
        abs(value)
        for feat in feats_to_show
        for key in ("q_means", "q_medians")
        for value in quantile_spreads[feat][key]
    )
    fig.update_yaxes(range=[-1.15 * span, 1.15 * span])
    fig.update_yaxes(title_text="Mean forward return", col=1)
    fig.update_layout(
        template="ml4t",
        height=260 * n_rows_fig + 60,
        width=900,
        title_t