Evaluación de señales transversales con IC, cuantiles y rotación

# Signal Evaluation: IC, Quantiles, and Spreads


# Signal Evaluation: IC, Quantiles, and Spreads

**Chapter 7: Defining the Learning Task**
**Section Reference**: 7.3 - Feature and Label Evaluation as Triage

**Docker image**: `ml4t`

## Purpose

This notebook demonstrates **single-factor evaluation** using Information Coefficient
(IC) analysis, quantile returns, and spread metrics. We answer: "Is this factor
predictive in the cross-section, and what horizon does it live on?"

## Learning Objectives

1. Compute cross-sectional IC and understand its time series properties
2. Interpret IC, ICIR, and HAC-adjusted significance
3. Analyze quantile returns, spread, and the monotonicity of the quantile ladder
4. Understand horizon comparison with proper overlap warnings
5. Measure turnover and signal half-life

## Data Policy

All examples use **real ETF data** from the case study store.
NO synthetic data is used in this notebook.

## Prerequisites

- `02_preprocessing_pipeline` - for split-aware preprocessing concepts that
  underpin fold-aware IC evaluation under **Fold-Aware Evaluation**.
- `03_label_methods` - supplies the forward-return labels used as `y_true`.
- Familiarity with rank correlations (Spearman) and walk-forward CV.

```python
"""Signal Evaluation - IC analysis, quintile spreads, and classification diagnostics for alpha signals."""

from __future__ import annotations

import json

import numpy as np
import plotly.graph_objects as go
import polars as pl
from IPython.display import display
from ml4t.diagnostic.evaluation.binary_metrics import (
    binary_classification_report,
    wilson_score_interval,
)
from ml4t.diagnostic.metrics import cross_sectional_ic, pooled_ic
from ml4t.diagnostic.signal import analyze_signal
from plotly.subplots import make_subplots
from scipy.stats import rankdata, spearmanr
from sklearn.metrics import (
    auc,
    confusion_matrix,
    precision_recall_curve,
    roc_auc_score,
    roc_curve,
)

from data import load_etfs
from utils.paths import get_chapter_dir
from utils.reproducibility import set_global_seeds
from utils.style import (  # importing utils.style activates the ml4t Plotly template
    COLORS,
    show_plotly_with_alt,
)
```

```python
SEED = 42
# Resolved from the chapter, not the working directory: the runner sets cwd to the chapter
# dir, so a repo-relative literal writes the publication artifact one level too deep and
# the book figure pipeline keeps reading an older copy at the intended path.
OUTPUT_DIR = get_chapter_dir(7) / "output"
START_DATE = "2006-01-01"
MAX_SYMBOLS = 0
N_PERMUTATIONS = 1000
DECAY_HORIZONS = (1, 2, 3, 5, 7, 10, 15, 21, 42)
N_SPLITS = 8
```

```python
set_global_seeds(SEED)
```

## Data Contract

Signal analysis requires two DataFrames with specific schemas:

**Factor Panel**:
- `timestamp`: Decision date
- `symbol`: Asset identifier
- `factor` or signal column(s): Factor value(s)

**Prices Panel**:
- `timestamp`: Same as factor panel
- `symbol`: Same as factor panel
- `close` or `price`: Closing price

The `ml4t-diagnostic` library aligns these using ASOF joins to compute
forward returns at each decision point.

```python
# Load real ETF data
etfs = load_etfs()
print(f"ETF universe: {etfs['symbol'].n_unique()} symbols, {len(etfs):,} rows")
print(f"Date range: {etfs['timestamp'].min()} to {etfs['timestamp'].max()}")
```

## Preparing Factor and Price Panels

We compute a simple momentum factor (21-day return) and prepare the data
in the format required by `analyze_signal()`.

The factor below is a 21-day return, used here as a teaching example. Production factors
come from the Chapter 8 feature pipelines.

```python
if START_DATE != "2006-01-01":
    etfs = etfs.filter(pl.col("timestamp") >= pl.lit(START_DATE).str.to_date())
if MAX_SYMBOLS > 0:
    keep = sorted(etfs["symbol"].unique().to_list())[:MAX_SYMBOLS]
    etfs = etfs.filter(pl.col("symbol").is_in(keep))

# Compute 21-day momentum per asset
factor_df = (
    etfs.sort(["symbol", "timestamp"])
    .with_columns(
        [(pl.col("close") / pl.col("close").shift(21).over("symbol") - 1).alias("factor")]
    )
    .filter(pl.col("factor").is_not_null())
    .select(["timestamp", "symbol", "factor"])
)

# Price panel
prices_df = etfs.select(["timestamp", "symbol", "close"]).rename({"close": "price"})

print(f"Factor panel: {factor_df.shape}")
print(f"Price panel: {prices_df.shape}")
print("Factor summary:")
display(factor_df.select("factor").describe())
```

```python
# Forward returns, reused by the fold-aware and binary sections below.
eval_df = (
    factor_df.join(prices_df, on=["timestamp", "symbol"], how="inner")
    .sort(["symbol", "timestamp"])
    .with_columns(
        fwd_21d=(pl.col("price").shift(-21).over("symbol") / pl.col("price") - 1),
    )
    .filter(pl.col("fwd_21d").is_not_null())
)
print(f"\nEvaluation panel: {eval_df.shape} (factor + 21D forward returns)")
```

### Correctness Screens

Before evaluating predictive power, verify that the factor is usable under the
stated protocol. The chapter's *Correctness screens* section prescribes four checks;
we demonstrate coverage
and staleness here. Timing/lag consistency and mask alignment become critical
with fundamental or third-party data (Chapters 8-10) but are trivially satisfied
for a price-derived momentum signal.

**Coverage** is the fraction of (date, asset) pairs carrying a non-null factor value,
reported overall and per date.

```python
all_pairs = prices_df.select("timestamp", "symbol").unique()
factor_pairs = factor_df.select("timestamp", "symbol").unique()
coverage = len(factor_pairs) / len(all_pairs)

# Per-date coverage (assets with factor / total assets)
daily_coverage = (
    all_pairs.join(factor_df, on=["timestamp", "symbol"], how="left")
    .group_by("timestamp")
    .agg(
        total=pl.len(),
        has_factor=pl.col("factor").is_not_null().sum(),
    )
    .with_columns(coverage=(pl.col("has_factor") / pl.col("total")))
    .sort("timestamp")
)

print("=== Correctness Screen: Coverage ===\n")
print(f"Overall coverage: {coverage:.1%}")
print(
    f"Per-date coverage - min: {daily_coverage['coverage'].min():.1%}, "
    f"median: {daily_coverage['coverage'].median():.1%}, "
    f"max: {daily_coverage['coverage'].max():.1%}"
)
print("\nNote: Coverage < 100% is expected - momentum requires 21 days of history,")
print("so new listings lack factor values during their first 21 trading days.")
```

**Staleness** asks whether the factor updates as often as its definition implies. A
price-derived momentum signal should take a new value every session, so a change rate
materially below one points at a data gap rather than at the signal.

```python
staleness = (
    factor_df.sort(["symbol", "timestamp"])
    .with_columns(
        factor_change=(pl.col("factor") != pl.col("factor").shift(1).over("symbol")).cast(pl.Int32)
    )
    .group_by("symbol")
    .agg(
        n_obs=pl.len(),
        n_changes=pl.col("factor_change").sum(),
    )
    .with_columns(change_rate=(pl.col("n_changes") / pl.col("n_obs")))
)

median_change_rate = staleness["change_rate"].median()
min_change_rate = staleness["change_rate"].min()

print("\n=== Correctness Screen: Staleness ===\n")
print(f"Median daily change rate: {median_change_rate:.1%}")
print(f"Min change rate (worst asset): {min_change_rate:.1%}")
if median_change_rate > 0.9:
    print("[PASS] Factor updates daily as expected for a price-derived signal.")
else:
    print("[WARNING] Some assets show stale factor values - investigate data gaps.")
```

## Information Coefficient (IC) Analysis

IC measures the **cross-sectional** rank correlation between signals and forward returns:

$$IC_t = \text{Spearman}(\text{signal}_{t}, \text{return}_{t \to t+h})$$

Where the correlation is computed across assets at each time $t$.

```python
# Run signal analysis
PERIODS = (1, 5, 21)  # Forward return horizons (days)
QUANTILES = 5  # Quintile analysis

result = analyze_signal(
    factor_df,
    prices_df,
    periods=PERIODS,
    quantiles=QUANTILES,
    ic_method="spearman",  # Rank correlation (robust to outliers)
    date_col="timestamp",
    asset_col="symbol",
)

print("Signal analysis complete")
print(f"Assets: {result.n_assets}, Dates: {result.n_dates}")
```

```python
# Information Coefficient by Horizon
ic_rows = []
for period in PERIODS:
    period_key = f"{period}D"
    ic_mean = result.ic.get(period_key, float("nan"))
    icir = result.ic_ir.get(period_key, float("nan"))
    t_stat = result.ic_t_stat.get(period_key, float("nan"))
    p_value = result.ic_p_value.get(period_key, float("nan"))
    sig = "***" if p_value < 0.01 else "**" if p_value < 0.05 else "*" if p_value < 0.10 else ""
    ic_rows.append(
        {
            "horizon": period_key,
            "mean_ic": round(ic_mean, 4),
            "icir": round(icir, 3),
            "t_stat": round(t_stat, 2),
            "p_value": round(p_value, 4),
            "sig": sig,
        }
    )

ic_summary = pl.DataFrame(ic_rows)
ic_summary
```

### Two IC Conventions: Pooled vs Cross-Sectional

The library exposes both, and the distinction matters for ranking strategies:

- **`pooled_ic`** - one global Spearman correlation across **all (date, asset)**
  observations.  Conflates *which days were good* with *which assets ranked
  correctly within a day*; sensitive to time-series mean shifts in returns.
- **`cross_sectional_ic`** - Spearman per date, then mean across dates.  Measures
  only the daily ranking skill that a long-short strategy can monetise, and exposes
  IC IR / t-stat / p-value on the per-date series.

Chapter 14 standardises on the cross-sectional convention. The two can disagree
materially on the same data - pooled may inflate or deflate magnitude depending
on the regime structure of returns.

```python
# Build a (date, symbol, y_pred, y_true) frame from eval_df at horizon 21D
ic_frame = eval_df.select(
    [
        pl.col("timestamp").alias("date"),
        pl.col("symbol"),
        pl.col("factor").alias("y_pred"),
        pl.col("fwd_21d").alias("y_true"),
    ]
).drop_nulls()

ic_pooled = pooled_ic(ic_frame["y_pred"], ic_frame["y_true"], method="spearman")
ic_xs = cross_sectional_ic(
    ic_frame,
    ic_frame,
    pred_col="y_pred",
    ret_col="y_true",
    date_col="date",
    entity_col="symbol",
    method="spearman",
    min_obs=5,
)

print(f"pooled_ic           : {ic_pooled:.4f}  (one global Spearman)")
print(
    f"cross_sectional_ic  : {ic_xs['ic_mean']:.4f}  "
    f"(mean of {ic_xs['n_periods']} daily Spearmans; "
    f"t={ic_xs['ic_t']:.2f}, p={ic_xs['p_value']:.4f})"
)
```

### Interpreting an IC magnitude

The right anchor for interpreting a mean IC is not the headline value but the standard
error of that mean, which is set by the number of periods $T$ in the daily-IC series and
by the dispersion of that series:

$$\text{SE}(\bar{\text{IC}}) \approx \frac{\sigma_{\text{IC}}}{\sqrt{T}}$$

A single point estimate therefore carries very different evidence depending on
$\sigma_{\text{IC}}$ and $T$. The cell below works the confidence interval for one
fixed $\bar{\text{IC}}$ under three sample-size and dispersion combinations, so the
arithmetic can be checked and re-run rather than read.

```python
IC_SCENARIO_MEAN = 0.02  # one headline IC, held fixed across the three scenarios
IC_SCENARIOS = (
    ("ten years, tight daily IC", 2_500, 0.05),
    ("ten years, wide daily IC", 2_500, 0.30),
    ("one year, wide daily IC", 250, 0.30),
)
Z_95 = 1.96

print(f"Mean IC held at {IC_SCENARIO_MEAN} in every row; only T and sigma move.\n")
print(f"{'scenario':<28}{'T':>7}{'sigma':>8}{'SE':>9}{'95% CI':>20}")
print("-" * 72)
for label, T, sigma in IC_SCENARIOS:
    se = sigma / np.sqrt(T)
    lo, hi = IC_SCENARIO_MEAN - Z_95 * se, IC_SCENARIO_MEAN + Z_95 * se
    print(f"{label:<28}{T:>7,}{sigma:>8.2f}{se:>9.4f}{f'[{lo:+.3f}, {hi:+.3f}]':>20}")
```

The same headline number is comfortably above zero in the first row, above zero but
uninformative about the signal's tail behaviour in the second, and indistinguishable
from zero in the third. Reporting a daily-mean IC therefore requires the CI (or the
$t$-statistic) alongside, and ideally the dispersion of the daily series as well.

The **ICIR** $= \bar{\text{IC}} / \sigma_{\text{IC}}$ is the signal-level analog of an
information ratio: $t \approx \text{ICIR} \times \sqrt{T}$ for serially uncorrelated
daily IC, and a HAC-adjusted $t$ for the realistic correlated case.

The bands printed next are typical magnitudes from the equity-factor literature on
multi-year daily-rebalanced cross-sectional studies. They are starting points for the SE
calculation above, not standalone readings, and they live in a declared table so the
code below can score against the same numbers the prose refers to.

```python
IC_BANDS = (
    (
        0.02,
        "at or below the daily-IC noise floor on multi-year samples; the SE alone often spans it",
    ),
    (0.04, "detectable on 5-to-10-year samples at the dispersions seen in published studies"),
    (
        0.06,
        "comparable to documented equity-factor effects such as momentum or short-term reversal",
    ),
    (
        0.10,
        "above most factor-zoo benchmarks; the question is whether it survives expanding-window evaluation",
    ),
    (
        float("inf"),
        "outside the published academic range; the prior is leakage or label corruption until ruled out",
    ),
)


def ic_band(ic: float) -> str:
    """Return the literature band an absolute mean IC falls into."""
    for upper, description in IC_BANDS:
        if abs(ic) < upper:
            return description
    return IC_BANDS[-1][1]


print("Reference bands for |mean IC| (equity-factor literature):")
lower = 0.0
for upper, description in IC_BANDS:
    edge = "and above" if upper == float("inf") else f"to {upper:.2f}"
    print(f"  {lower:.2f} {edge:<12} {description}")
    lower = upper

print(f"\nThis factor's 21-day mean IC: {result.ic.get('21D', float('nan')):.4f}")
print(f"  falls in: {ic_band(result.ic.get('21D', float('nan')))}")
```

The daily IC series is plotted at low opacity behind two rolling means. At this sample
size the raw series is a solid band and carries no readable structure on its own; what a
reader can act on is whether the smoothed level drifts away from zero, and the two
windows show whether an apparent drift is still there under a longer average.

```python
fig = make_subplots(
    rows=1, cols=2, subplot_titles=["Daily IC, 21-day horizon", "Distribution of daily IC"]
)

# Get 21D IC series
ic_21d = result.ic_series.get("21D", [])
if ic_21d:
    # Raw series as context only: 5,000 overlapping daily values plot as a solid band.
    fig.add_trace(
        go.Scatter(
            y=ic_21d,
            mode="lines",
            name="Daily IC",
            line=dict(color=COLORS["neutral"], width=0.4),
            opacity=0.25,
        ),
        row=1,
        col=1,
    )

    ic_series = pl.Series(ic_21d)
    for window, color, width in ((63, COLORS["amber"], 1.5), (252, COLORS["blue"], 2.0)):
        fig.add_trace(
            go.Scatter(
                y=ic_series.rolling_mean(window_size=window).to_list(),
                mode="lines",
                name=f"{window}-day mean",
                line=dict(color=color, width=width),
            ),
            row=1,
            col=1,
        )

    fig.add_hline(y=0, line_dash="dash", line_color=COLORS["neutral"], row=1, col=1)

    # Distribution
    fig.add_trace(
        go.Histogram(x=ic_21d, nbinsx=30, name="IC Distribution", marker_color=COLORS["blue"]),
        row=1,
        col=2,
    )

    mean_ic = np.mean(ic_21d)
    fig.add_vline(x=mean_ic, line_dash="dash", line_color=COLORS["negative"], row=1, col=2)

fig.update_xaxes(title_text="Trading day (index, chronological)", row=1, col=1)
fig.update_yaxes(title_text="IC (Spearman)", row=1, col=1)
fig.update_xaxes(title_text="IC", row=1, col=2)
fig.update_yaxes(title_text="Count", row=1, col=2)
fig.update_layout(
    height=350, showlegend=True, title_text="Daily cross-sectional IC of the momentum factor"
)
show_plotly_with_alt(
    fig,
    alt=(
        "Two panels. The left panel plots the daily cross-sectional IC against a "
        "chronological index of roughly five thousand trading days. The raw series is a "
        "faint grey band filling the range from about minus 0.9 to plus 0.9 with no "
        "visible trend, and two rolling means are drawn over it. The shorter amber mean "
        "swings within roughly a fifth of a correlation unit either side of zero; the "
        "longer navy mean is flatter still and stays inside about half that. Neither "
        "settles at a level away from zero anywhere in the sample. The right panel is a "
        "histogram of the same "
        "daily values: a broad, roughly symmetric bell centred on zero and running out to "
        "about plus and minus 0.8, with a dashed red line marking the mean sitting "
        "essentially on the zero gridline."
    ),
)
```

### Publication Figure Artifact

The book IC time-series figure reads a compact NumPy artifact so formatting
changes do not reload the ETF panel or recompute daily cross-sectional ICs.

The pairs below are sorted on the native timestamp rather than on its string form.
Lexicographic sorting of stringified dates is correct only for zero-padded ISO output
and would silently reorder the series under any other rendering.

```python
ic_pairs: list[tuple[object, float]] = []
for date_df in eval_df.partition_by("timestamp"):
    if len(date_df) < 20:
        continue
    corr, _ = spearmanr(date_df["factor"].to_numpy(), date_df["fwd_21d"].to_numpy())
    if not np.isnan(corr):
        ic_pairs.append((date_df["timestamp"][0], float(corr)))

ic_pairs.sort(key=lambda t: t[0])
ic_dates_for_figure = [str(ts) for ts, _ in ic_pairs]
ic_values_arr = np.array([v for _, v in ic_pairs])

rolling_window = 63
rolling_ic_for_figure = np.full_like(ic_values_arr, np.nan)
for i in range(rolling_window, len(ic_values_arr)):
    rolling_ic_for_figure[i] = np.mean(ic_values_arr[i - rolling_window : i])

min_train = int(len(ic_dates_for_figure) * 0.3)
test_size = 252
fold_boundaries_for_figure = ic_dates_for_figure[min_train::test_size]

OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
figure_7_3_artifact = OUTPUT_DIR / "figure_7_3_ic_time_series_with_folds.npz"
np.savez(
    figure_7_3_artifact,
    ic_dates=np.array(ic_dates_for_figure),
    ic_values=ic_values_arr,
    rolling_ic=rolling_ic_for_figure,
    fold_boundaries=np.array(fold_boundaries_for_figure),
)
print(f"Wrote publication figure artifact: {figure_7_3_artifact}")
```

## Quantile Analysis

Examine returns by signal quantile to assess **monotonicity** (do higher signal
values lead to higher returns?) and **spread** (what's the return difference
between top and bottom quantiles?).

```python
# Quantile returns
print("\n=== Mean Returns by Quantile ===\n")

for period in PERIODS:
    period_key = f"{period}D"
    quantile_rets = result.quantile_returns.get(period_key, {})

    if quantile_rets:
        print(f"{period_key} Forward Returns:")
        for quantile in sorted(quantile_rets.keys()):
            ret = quantile_rets[quantile]
            bar = "█" * int(abs(ret) * 500)  # Simple text bar
            sign = "+" if ret >= 0 else ""
            print(f"  Q{quantile}: {sign}{ret:.4%} {bar}")
        print()
```

```python
# Visualize quantile returns
fig = make_subplots(
    rows=1,
    cols=len(PERIODS),
    subplot_titles=[f"{p}D Forward Returns" for p in PERIODS],
    horizontal_spacing=0.12,
)

for i, period in enumerate(PERIODS, 1):
    period_key = f"{period}D"
    period_returns = result.quantile_returns.get(period_key, {})

    if period_returns:
        quantiles = [f"Q{q}" for q in sorted(period_returns.keys())]
        returns = [period_returns[q] for q in sorted(period_returns.keys())]

        fig.add_trace(
            go.Bar(
                x=quantiles,
                y=returns,
                name=period_key,
                showlegend=False,
                marker_color=[COLORS["amber"] if v < 0 else COLORS["blue"] for v in returns],
            ),
            row=1,
            col=i,
        )

for i in range(1, len(PERIODS) + 1):
    fig.update_xaxes(title_text="Quantile", row=1, col=i)
    fig.update_yaxes(tickformat=".2%", row=1, col=i)
# Label the shared y-axis only on the first panel to avoid the title overlapping
# the neighbouring panel's bars; each panel keeps its own (per-horizon) scale.
fig.update_yaxes(title_text="Mean forward return", row=1, col=1)
fig.update_layout(title="Mean forward return by signal quantile, three horizons", height=350)
show_plotly_with_alt(
    fig,
    alt=(
        "Three bar panels, one per forward-return horizon, each showing the mean forward "
        "return of the five signal quantiles from Q1 to Q5. Every bar in every panel is "
        "positive, so all five quantiles earned money over the sample. The ordering does "
        "not follow the signal: in the one-day panel Q1 is the tallest bar and Q4 the "
        "shortest; in the five-day panel the bars fall from Q1 to Q4 and tick up at Q5; in "
        "the twenty-one-day panel the five bars are nearly level with Q5 the lowest. "
        "Each panel carries its own vertical scale, so the heights are comparable within a "
        "panel and not across panels."
    ),
)
```

The bars do not rise with the quantile in any of the three panels, and in the shortest
horizon they fall - the bottom quantile of the signal earns the most. A long-short book
built the obvious way round would be short the better-performing leg. The spread and
monotonicity numbers below put figures on that, and the sign is the part to read first.

```python
# Spread and monotonicity analysis
spread_rows = []
for period in PERIODS:
    period_key = f"{period}D"
    spread = result.spread.get(period_key, float("nan"))
    t_stat = result.spread_t_stat.get(period_key, float("nan"))
    mono = result.monotonicity.get(period_key, float("nan"))
    spread_rows.append(
        {
            "horizon": period_key,
            "spread_pct": round(spread * 100, 4),
            "t_stat": round(t_stat, 2),
            "monotonicity_rho": round(mono, 3),
        }
    )

spread_summary = pl.DataFrame(spread_rows)
spread_summary
```

### Monotonicity Interpretation

`analyze_signal` reports monotonicity as the **Spearman rank correlation between the
quantile ranks and their mean returns**, so it runs from $-1$ to $+1$: $+1$ is a
perfectly increasing ladder of quantile returns, $-1$ a perfectly decreasing one, and
$0$ no rank association at all. It is a correlation, not a percentage, and is reported
as such below.

That is not the same statistic as the one the fold section computes further down, which
is the *fraction of adjacent quantile steps that move upward* and runs from $0$ to $1$.
The two answer the same question differently and do not share a scale; where both appear
they are named apart.

The bands below score the **strength** of the ordering, so they read the magnitude. The
sign is a separate fact and the more important one here: a strongly ordered signal
pointing the wrong way is not a weak signal, it is an inverted one, and the two call for
opposite actions. The bands are kept in a declared table so the reading the prose gives
is the one the code applies.

```python
MONOTONICITY_BANDS = (
    (0.60, "weak ordering; consider a non-linear model"),
    (0.80, "moderate ordering; may work as a long-short book"),
    (1.01, "strong, consistent ordering"),
)


def monotonicity_band(value: float) -> str:
    """Return the strength convention the magnitude of a monotonicity rho falls into."""
    for upper, description in MONOTONICITY_BANDS:
        if abs(value) < upper:
            return description
    return MONOTONICITY_BANDS[-1][1]


lower = 0.0
print("Reference bands for |rho| (ordering strength; the sign is read separately):")
for upper, description in MONOTONICITY_BANDS:
    print(f"  {lower:.2f} to {min(upper, 1.0):.2f}  {description}")
    lower = upper

print()
for period in PERIODS:
    mono = result.monotonicity.get(f"{period}D", float("nan"))
    direction = "as signalled" if mono > 0 else "inverted" if mono < 0 else "flat"
    print(f"  {period}D: {mono:+.2f}  {direction}, {monotonicity_band(mono)}")
```

## Horizon Comparison

Compare IC across different forward return horizons to identify the optimal
holding period for the signal.

```python
fig = make_subplots(rows=1, cols=2, subplot_titles=["Mean IC by horizon", "ICIR by horizon"])

horizons = list(PERIODS)
ics = [result.ic.get(f"{h}D", float("nan")) for h in horizons]
icirs = [result.ic_ir.get(f"{h}D", float("nan")) for h in horizons]

# IC plot (colorblind-safe blue/amber)
fig.add_trace(
    go.Bar(
        x=[f"{h}D" for h in horizons],
        y=ics,
        name="Mean IC",
        marker_color=[COLORS["blue"] if v > 0 else COLORS["amber"] for v in ics],
    ),
    row=1,
    col=1,
)

# ICIR plot
icir_colors = [
    COLORS["blue"] if v > 0.5 else COLORS["slate"] if v > 0 else COLORS["amber"] for v in icirs
]
fig.add_trace(
    go.Bar(x=[f"{h}D" for h in horizons], y=icirs, name="ICIR", marker_color=icir_colors),
    row=1,
    col=2,
)

fig.add_hline(y=0.5, line_dash="dash", line_color=COLORS["neutral"], row=1, col=2)

fig.update_xaxes(title_text="Horizon", row=1, col=1)
fig.update_yaxes(title_text="Mean IC", row=1, col=1)
fig.update_xaxes(title_text="Horizon", row=1, col=2)
fig.update_yaxes(title_text="ICIR", row=1, col=2)
fig.update_layout(
    height=350, showlegend=False, title_text="Mean IC and ICIR across the three evaluated horizons"
)
show_plotly_with_alt(
    fig,
    alt=(
        "Two bar panels over the same three horizons, one day, five days and twenty-one "
        "days. The left panel plots mean IC: the one-day and five-day bars hang below "
        "zero in amber, the five-day the longer of the two, and the twenty-one-day bar "
        "rises just above zero in navy. The whole vertical range spans less than one "
        "hundredth of a correlation unit. The right panel plots ICIR on an axis whose top "
        "is a dashed reference line; all three bars are so close to zero that they read as "
        "a flat line along the baseline, nowhere near that reference."
    ),
)
```

Read the two panels together. The left panel's bars change sign across horizons, which
invites a story about reversal giving way to momentum; the right panel says the whole
left panel is inside the noise, because an ICIR of essentially zero means the daily IC
series has a standard deviation orders of magnitude larger than its mean. The sign
pattern is not a finding at these magnitudes.

### Overlapping Returns Warning

When comparing IC across horizons, be aware that **longer horizons have overlapping
forward returns**, which introduces autocorrelation in the IC series.

For example, 21-day forward returns on consecutive days share 20 days of overlap.
This means:

1. **IC series are autocorrelated** at longer horizons
2. **Standard errors are understated** without HAC adjustment
3. **Comparing IC across horizons requires caution** - higher IC at longer
   horizons may reflect overlap, not better predictability

**Best practice**: Use HAC-adjusted t-statistics (as provided by `analyze_signal`)
and be skeptical of IC that increases monotonically with horizon.

### IC Decay Analysis

IC decay determines the optimal rebalancing frequency: a signal whose IC halves within a
week should not be held for a month. We compute IC on a finer horizon grid to estimate
the signal's useful life.

The half-peak marker below is drawn only for a crossing that comes **after** the peak.
Searching from the shortest horizon instead puts the marker on the first grid point
whenever the profile rises, which is the opposite of decay and reads to a reader as
decay. Where IC never falls back within the grid - the peak being the longest horizon on
it - no marker is drawn at all.

```python
# Compute IC across a finer horizon grid (single call - batches all horizons)
decay_horizons = DECAY_HORIZONS

decay_result = analyze_signal(
    factor_df,
    prices_df,
    periods=decay_horizons,
    quantiles=3,
    ic_method="spearman",
    date_col="timestamp",
    asset_col="symbol",
)
decay_ics = [decay_result.ic.get(f"{h}D", float("nan")) for h in decay_horizons]
```

```python
# Plot IC decay curve
fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x=decay_horizons,
        y=decay_ics,
        mode="lines+markers",
        name="Mean IC",
        line=dict(width=2),
    )
)
fig.add_hline(y=0, line_dash="dash", line_color=COLORS["neutral"])

# The crossing search starts after the peak; see the markdown above this figure.
peak_idx = int(np.argmax(decay_ics))
peak_ic = decay_ics[peak_idx]
half_ic = peak_ic / 2
crossing = None
if peak_ic > 0:
    for i in range(peak_idx + 1, len(decay_ics)):
        if decay_ics[i] < half_ic:
            crossing = decay_horizons[i]
            break
if crossing is not None:
    fig.add_vline(
        x=crossing,
        line_dash="dot",
        line_color=COLORS["negative"],
        annotation_text=f"IC below half its peak by day {crossing}",
    )

fig.update_layout(
    title="Mean IC against forward-return horizon",
    xaxis_title="Forward return horizon (days)",
    yaxis_title="Mean IC (Spearman)",
    height=350,
)
show_plotly_with_alt(
    fig,
    alt=(
        "A line with markers showing mean IC against forward-return horizon over a grid "
        "running from one to forty-two days. The whole vertical range spans less than "
        "one hundredth of a correlation unit. The line starts negative at the shortest "
        "horizons, reaches its lowest point around five days, climbs through the dashed "
        "zero line between ten and fifteen days, dips slightly at twenty-one and ends at "
        "its highest point at forty-two days. There is no peak followed by decay: the "
        "profile is still rising where the grid stops."
    ),
)
```

**Interpretation**: The horizon-IC curve is the tool for choosing a rebalancing
frequency - for a cleanly decaying signal you rebalance near where IC peaks and stop
before it fades. This factor does *not* show that decay. Every horizon IC printed above
sits within a few thousandths of zero, none is statistically distinct from zero, and the
profile is still rising at the longest horizon on the grid rather than falling away from
a peak. The half-peak marker is drawn only when a crossing exists *after* the peak, so
on this profile there is nothing to mark; a marker placed at the first horizon would
have described decay running backwards. The honest read is a factor with no exploitable
cross-sectional horizon structure on this universe.

## Turnover Analysis

High turnover erodes returns through transaction costs. A signal with high IC
but excessive turnover may not be profitable after costs.

```python
# Turnover metrics
print("\n=== Turnover Analysis ===\n")

if result.turnover:
    print("Mean Turnover by Period:")
    for period_key, turnover in result.turnover.items():
        print(f"  {period_key}: {turnover:.1%}")
else:
    print("Turnover: Not computed")

if result.autocorrelation:
    print(
        f"\nSignal Autocorrelation (lag 1-5): {[f'{ac:.3f}' for ac in result.autocorrelation[:5]]}"
    )
else:
    print("\nAutocorrelation: Not computed")

if result.half_life:
    print(f"\nSignal Half-Life: {result.half_life:.1f} periods")

    # Interpretation
    if result.half_life < 5:
        print("  Interpretation: Fast decay - requires frequent rebalancing")
    elif result.half_life < 20:
        print("  Interpretation: Moderate decay - weekly rebalancing appropriate")
    else:
        print("  Interpretation: Slow decay - monthly rebalancing sufficient")
```

### Turnover and Costs

A simple cost-adjusted IC (Grinold approximation):

$$IC_{net} \approx IC - \frac{c \times \text{turnover}}{E[r]}$$

Where $c$ is round-trip transaction cost and $E[r]$ is expected return. For most equity
strategies, turning the whole book over once a month or more erodes alpha materially.

### Break-Even Cost Analysis

A feasibility check asks: **could this signal survive transaction costs?**

We compare the expected spread (top-bottom quantile return) to the cost of
achieving that spread. If round-trip costs exceed the expected edge, the
signal is not tradeable at the given rebalancing frequency.

```python
# Break-even cost analysis
print("\n=== Break-Even Cost Analysis ===\n")

# Get the 21-day spread and turnover
spread_21d = result.spread.get("21D", float("nan"))
turnover_21d = result.turnover.get("21D", 0.5) if result.turnover else 0.5  # default 50%

# Define cost assumptions (conservative for US equities)
# These should match the trading setup's cost model
COST_ASSUMPTIONS = {
    "spread_bps": 10,  # Half-spread in basis points
    "commission_bps": 5,  # Commission per side
    "market_impact_bps": 10,  # Expected market impact
}

round_trip_cost = (
    2 * COST_ASSUMPTIONS["spread_bps"]
    + 2 * COST_ASSUMPTIONS["commission_bps"]
    + 2 * COST_ASSUMPTIONS["market_impact_bps"]
) / 10000  # Convert to decimal

print("Cost Assumptions (per leg):")
print(f"  Spread:        {COST_ASSUMPTIONS['spread_bps']} bps")
print(f"  Commission:    {COST_ASSUMPTIONS['commission_bps']} bps")
print(f"  Market impact: {COST_ASSUMPTIONS['market_impact_bps']} bps")
print(f"  Round-trip:    {round_trip_cost * 10000:.0f} bps ({round_trip_cost:.2%})")
```

```python
# Compute cost-adjusted spread
expected_turnover_per_period = turnover_21d * 2  # Both legs
cost_drag = round_trip_cost * expected_turnover_per_period

cost_adjusted_spread = spread_21d - cost_drag

print("\n21D Signal Analysis:")
print(f"  Raw spread:            {spread_21d:.2%}")
print(f"  Expected turnover:     {expected_turnover_per_period:.0%}")
print(f"  Cost drag per period:  {cost_drag:.2%}")
print(f"  Cost-adjusted spread:  {cost_adjusted_spread:.2%}")

# Break-even calculation
if spread_21d > 0:
    break_even_cost = spread_21d / expected_turnover_per_period
    print(f"  Break-even cost:       {break_even_cost * 10000:.0f} bps")

    if cost_adjusted_spread > 0:
        print("\n[PASS] Signal survives cost assumptions")
    else:
        print("\n[WARNING] Signal does not survive cost assumptions at this turnover")
        print("   Consider: longer horizon, lower turnover, or reduced position sizing")
else:
    print("\n[WARNING] Negative spread - signal direction may be inverted")
```

### Feasibility Guidelines

Three checks for signal feasibility, from the chapter's *Preliminary feasibility
checks* section:

1. **Turnover proxies**: Measure entry/exit rates in the top-k set (see above)
2. **Break-even cost checks**: Compare spread to conservative cost estimates
3. **Capacity warnings**: Recompute IC by liquidity bucket (deferred to Ch8)

**Note**: Liquidity-bucket analysis requires actual market microstructure data
(average volume, bid-ask spreads) which is not available for this placeholder
feature. Chapter 8 demonstrates this check with real case study data.

## Fold-Aware Evaluation

**Critical**: The IC computed above pools all dates into a single statistic. However,
real trading strategies are evaluated on **out-of-sample** data using walk-forward
validation. This section demonstrates fold-aware IC computation.

### Why Folds Matter

Global IC can be misleading because:
1. **Regime dependence**: IC may be high in some periods and zero in others
2. **Lookahead contamination**: Parameters tuned on full data leak future information
3. **Overfitting detection**: Consistent IC across folds suggests robust signal

The text emphasizes computing IC **per fold** and reporting the distribution of
fold-level statistics, not just their pooled mean.

The splits below use an expanding window: train on everything up to the split point,
test on the period that follows.

```python
def create_walk_forward_splits(
    dates: list, n_splits: int = 5, min_train_pct: float = 0.2, test_periods: int = 63
) -> list[tuple[list, list]]:
    """Create walk-forward cross-validation splits.

    Args:
        dates: Sorted unique dates
        n_splits: Number of test folds
        min_train_pct: Minimum training data as fraction of total
        test_periods: Number of periods per test fold

    Returns:
        List of (train_dates, test_dates) tuples
    """
    n_dates = len(dates)
    min_train = int(n_dates * min_train_pct)

    splits = []
    for i in range(n_splits):
        # Test window
        test_start = min_train + i * test_periods
        test_end = min(test_start + test_periods, n_dates)

        if test_start >= n_dates:
            break

        train_dates = dates[:test_start]
        test_dates = dates[test_start:test_end]

        if len(test_dates) > 0:
            splits.append((train_dates, test_dates))

    return splits
```

```python
# Create splits for our data
unique_dates = sorted(factor_df["timestamp"].unique().to_list())
n_splits = N_SPLITS
test_periods = 63  # ~3 months per fold

splits = create_walk_forward_splits(
    unique_dates, n_splits=n_splits, min_train_pct=0.3, test_periods=test_periods
)

print(f"Created {len(splits)} walk-forward splits")
print("\nSplit structure:")
for i, (train_dates, test_dates) in enumerate(splits):
    print(f"  Fold {i + 1}: Train {len(train_dates)} days, Test {len(test_dates)} days")
    print(f"           Train: {train_dates[0]} to {train_dates[-1]}")
    print(f"           Test:  {test_dates[0]} to {test_dates[-1]}")
```

### Compute Per-Fold IC

For each fold, we compute IC on the **test period only**. This mimics how the
signal would perform in live trading, where we only see future returns after
making predictions.

```python
# Slice eval_df per fold; the forward returns were computed once above.
fold_results = []
evaluated_dates: list = []  # every date inside a test window - needed for a like-for-like IC

for fold_idx, (train_dates, test_dates) in enumerate(splits):
    test_data = eval_df.filter(pl.col("timestamp").is_in(test_dates))

    if len(test_data) < 100:
        continue

    # Cross-sectional IC per date, then average
    ic_per_date = []
    for date_df in test_data.partition_by("timestamp"):
        if len(date_df) < 10:  # Need enough assets for meaningful correlation
            continue
        corr, _ = spearmanr(date_df["factor"].to_numpy(), date_df["fwd_21d"].to_numpy())
        if not np.isnan(corr):
            ic_per_date.append(corr)

    if not ic_per_date:
        continue

    fold_ic = np.mean(ic_per_date)

    # Quantile spread, and the adjacent-step ordering fraction (NOT the Spearman rho
    # analyze_signal reports; see "Monotonicity Interpretation" above).
    test_with_q = test_data.with_columns(
        quantile=pl.col("factor")
        .rank()
        .over("timestamp")
        .qcut(5, labels=[str(i) for i in range(1, 6)])
        .over("timestamp")
    )
    q_rets = test_with_q.group_by("quantile").agg(pl.col("fwd_21d").mean()).sort("quantile")
    q_vals = q_rets["fwd_21d"].to_list()
    fold_spread = q_vals[-1] - q_vals[0] if len(q_vals) >= 2 else float("nan")
    # Monotonicity: fraction of consecutive quantiles in correct order
    if len(q_vals) >= 2:
        diffs = [q_vals[i + 1] - q_vals[i] for i in range(len(q_vals) - 1)]
        fold_step_frac = sum(1 for d in diffs if d > 0) / len(diffs)
    else:
        fold_step_frac = float("nan")

    evaluated_dates.extend(test_dates)
    fold_results.append(
        {
            "fold": fold_idx + 1,
            "test_start": str(test_dates[0]),
            "test_end": str(test_dates[-1]),
            "n_obs": len(test_data),
            "ic": fold_ic,
            "spread": fold_spread,
            "up_step_fraction": fold_step_frac,
        }
    )
```

```python
fold_df = pl.DataFrame(fold_results)
print("\n=== Per-Fold IC Results (21D horizon) ===\n")
print(fold_df)
```

```python
# Summarize fold-level statistics
fold_ic_mean = fold_df["ic"].mean()
fold_ic_std = fold_df["ic"].std()
fold_ic_min = fold_df["ic"].min()
fold_ic_max = fold_df["ic"].max()
pct_positive = (fold_df["ic"] > 0).mean() * 100

print("\n=== Fold-Level Summary ===")
print(f"Mean IC:       {fold_ic_mean:.4f}")
print(f"Std IC:        {fold_ic_std:.4f}")
print(f"IC Range:      [{fold_ic_min:.4f}, {fold_ic_max:.4f}]")
print(f"% Folds > 0:   {pct_positive:.0f}%")
print(f"Fold ICIR:     {fold_ic_mean / fold_ic_std:.3f}" if fold_ic_std > 0 else "N/A")
```

```python
fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=["IC by test fold", "Distribution of fold IC"],
    horizontal_spacing=0.15,
)

# Fold test-start dates as x-axis labels, for temporal context
fold_labels = [r["test_start"][:7] for r in fold_results]  # YYYY-MM format

# Colorblind-safe: blue for positive, amber for negative
fig.add_trace(
    go.Bar(
        x=fold_labels,
        y=[r["ic"] for r in fold_results],
        marker_color=[COLORS["blue"] if r["ic"] > 0 else COLORS["amber"] for r in fold_results],
        name="Fold IC",
    ),
    row=1,
    col=1,
)

# Add global mean line
fig.add_hline(
    y=fold_ic_mean,
    line_dash="dash",
    line_color=COLORS["blue"],
    line_width=1.5,
    annotation_text=f"Mean={fold_ic_mean:.3f}",
    row=1,
    col=1,
)
fig.add_hline(y=0, line_dash="dot", line_color=COLORS["neutral"], row=1, col=1)

# Histogram of IC values
fig.add_trace(
    go.Histogram(x=[r["ic"] for r in fold_results], nbinsx=10, marker_color=COLORS["blue"]),
    row=1,
    col=2,
)
fig.add_vline(x=0, line_dash="dot", line_color=COLORS["neutral"], row=1, col=2)
fig.add_vline(x=fold_ic_mean, line_dash="dash", line_color=COLORS["blue"], row=1, col=2)

fig.update_xaxes(title_text="Test Fold Start", row=1, col=1)
fig.update_yaxes(title_text="IC (Spearman)", row=1, col=1)
fig.update_xaxes(title_text="IC", row=1, col=2)
fig.update_yaxes(title_text="Count", row=1, col=2)

fig.update_layout(
    height=400,
    showlegend=False,
    font=dict(size=12),
    title_text="Out-of-sample IC across the walk-forward test folds",
)
show_plotly_with_alt(
    fig,
    alt=(
        "Two panels covering the eight walk-forward test folds, whose start dates run "
        "from early 2012 to late 2013. The left panel is a bar per fold: seven bars stand "
        "above zero in navy and one amber bar hangs below it, the negative fold being much "
        "the shortest of the eight. A dashed line marks the fold mean, well above zero, "
        "and a dotted line marks zero itself. The right panel is a histogram of those "
        "eight values, sparse by construction, with most of the mass just above zero and a "
        "single isolated bar far to the right."
    ),
)
```

### Interpretation: Full-Sample vs Fold-Level IC

The obvious move is to read the full-sample IC against the fold-level mean and
call the difference an aggregation effect - the folds are out of sample, so a
gap looks like the honest shrinkage of an optimistic full-sample number.

That reading is only available if both statistics cover the **same dates**.
Ours do not. The full-sample IC averages every date in the panel. The fold-level
mean averages only dates inside a test window, and with the `min_train_pct` and fold
count set above those windows are a few hundred consecutive dates near the front of the
sample - the printout below gives the span and the share. The two numbers describe
different periods, so their
difference cannot be attributed to aggregation - or to leakage, which the folds
structurally cannot produce, since each fold's IC only ever touches its own test
dates.

The way out is a third number: the same cross-sectional IC, restricted to
exactly the dates the folds evaluated. It splits the gap into two parts we can
name separately:

| Comparison | Isolates |
|------------|----------|
| Fold-date IC vs full-sample IC | **Period** - is this window unrepresentative? |
| Fold-level mean vs fold-date IC | **Aggregation** - does averaging folds change anything? |

```python
# The like-for-like statistic: same per-date Spearman, restricted to fold dates
fold_window = eval_df.filter(pl.col("timestamp").is_in(evaluated_dates))
fold_window_ics = []
for date_df in fold_window.partition_by("timestamp"):
    if len(date_df) < 10:
        continue
    corr, _ = spearmanr(date_df["factor"].to_numpy(), date_df["fwd_21d"].to_numpy())
    if not np.isnan(corr):
        fold_window_ics.append(corr)
fold_window_ic = float(np.mean(fold_window_ics))

full_sample_ic = result.ic.get("21D", float("nan"))
n_all_dates = eval_df["timestamp"].n_unique()
n_fold_dates = len(fold_window_ics)

print("\n=== Like-for-Like IC Comparison (21D) ===")
print(f"Full-sample cross-sectional IC ({n_all_dates:,} dates): {full_sample_ic:>8.4f}")
print(f"Same statistic, fold dates only ({n_fold_dates:,} dates): {fold_window_ic:>8.4f}")
print(f"Fold-level mean IC              ({n_fold_dates:,} dates): {fold_ic_mean:>8.4f}")
print("\n--- Decomposition of the full-sample vs fold-level gap ---")
print(
    f"Period effect      (fold-date IC - full-sample IC): {fold_window_ic - full_sample_ic:>8.4f}"
)
print(f"Aggregation effect (fold mean    - fold-date IC):   {fold_ic_mean - fold_window_ic:>8.4f}")

fold_span = f"{min(evaluated_dates)} to {max(evaluated_dates)}"
coverage_pct = 100 * n_fold_dates / n_all_dates
print(f"\nFolds evaluate {fold_span} - {coverage_pct:.0f}% of the panel's dates.")
if abs(fold_window_ic - full_sample_ic) > abs(fold_ic_mean - fold_window_ic):
    print("[READ] The gap is a period effect: this window is not the average window.")
    print("       It says nothing about overfitting, and nothing about leakage.")
else:
    print("[READ] The gap survives on identical dates, so it is an aggregation effect.")
```

### Within-Time Permutation Test

The chapter recommends a **within-time permutation test** as a
null-distribution benchmark: break the feature-label pairing while preserving the
structure of the data, then ask how often chance alone reproduces the observed IC.
Two details decide whether the answer means anything.

**The null must cover the same dates as the observed statistic.** A null built
from every date in the panel is a distribution of a mean over ~5,000 dates; our
observed statistic is a mean over ~500. The mean of a larger sample is mechanically
less dispersed, so such a null is too narrow by roughly $\sqrt{5000/500} \approx 3$
and would reject almost anything. We therefore permute only the fold dates.

**The null must preserve temporal dependence.** Shuffling assets independently on
each date implies the per-date ICs are independent, so the null mean's standard
error shrinks like $\sigma/\sqrt{n_{\text{dates}}}$. Our labels are 21-day forward
returns sampled daily: consecutive dates share 20 of 21 days of return, so both the
returns and the per-date ICs are strongly autocorrelated, and the true standard
error is much larger. Independent within-date shuffling would understate it by
roughly an order of magnitude.

The repair is a **block permutation**. We draw one asset relabeling per block of 21
sessions - the label horizon - and hold it fixed across the block. Within a block
each asset's return path stays intact, so the autocorrelation the observed IC
inherits also lives in the null; across blocks the relabelings are independent.
The null then answers the right question: *given how persistent this data is, how
often does a random assignment rank this well over this window?*

```python
# Block permutation: one asset relabeling per BLOCK_SESSIONS, restricted to fold dates
rng = np.random.default_rng(SEED)
n_permutations = N_PERMUTATIONS
BLOCK_SESSIONS = 21  # label horizon - blocks must be at least as long as the overlap

perm_df = fold_window.sort(["timestamp", "symbol"])
dates_arr = perm_df["timestamp"].to_numpy()
factors_arr = perm_df["factor"].to_numpy()
returns_arr = perm_df["fwd_21d"].to_numpy()
symbol_codes = perm_df["symbol"].cast(pl.Categorical).to_physical().to_numpy()
n_symbols = int(symbol_codes.max()) + 1
unique_dates_perm = np.unique(dates_arr)

# Group rows by date; within a date, rows are already in symbol order
date_groups = [np.where(dates_arr == d)[0] for d in unique_dates_perm]
date_groups = [idx for idx in date_groups if len(idx) >= 10]

# Ranks are invariant to relabeling, so pre-compute both sides once
factor_ranks_by_group = [rankdata(factors_arr[idx]) for idx in date_groups]
return_ranks_by_group = [rankdata(returns_arr[idx]) for idx in date_groups]
symbols_by_group = [symbol_codes[idx] for idx in date_groups]
```

```python
permuted_ics = []
for _ in range(n_permutations):
    ic_per_date = []
    block_keys = None
    for i, f_ranks in enumerate(factor_ranks_by_group):
        # New relabeling only when a block boundary is crossed
        if i % BLOCK_SESSIONS == 0:
            block_keys = rng.permutation(n_symbols)
        # Same key vector across the block => same asset->asset mapping across the block
        order = np.argsort(block_keys[symbols_by_group[i]], kind="stable")
        r_ranks = return_ranks_by_group[i][order]
        corr = np.corrcoef(f_ranks, r_ranks)[0, 1]
        if not np.isnan(corr):
            ic_per_date.append(corr)
    if ic_per_date:
        permuted_ics.append(np.mean(ic_per_date))

permuted_ics = np.array(permuted_ics)
# Observed statistic and null are now the same estimator on the same dates
observed_ic_mean = fold_window_ic

# (r + 1) / (B + 1): a permutation p-value can never be exactly zero - the observed
# assignment is itself one of the arrangements under the null
n_at_least = int(np.sum(permuted_ics >= observed_ic_mean))
p_value_perm = (n_at_least + 1) / (len(permuted_ics) + 1)

print("=== Within-Time Block Permutation Test ===\n")
print(
    f"Dates: {len(factor_ranks_by_group):,} (fold windows only, block = {BLOCK_SESSIONS} sessions)"
)
print(f"Observed mean IC: {observed_ic_mean:.4f}")
print(f"Permutation null - mean: {permuted_ics.mean():.4f}, std: {permuted_ics.std():.4f}")
print(f"Permutation p-value (one-sided): {p_value_perm:.4f}")
print(f"Permutations: {n_permutations} (finest resolvable p = {1 / (n_permutations + 1):.4f})")
```

```python
# Visualize permutation distribution vs observed IC
fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=permuted_ics,
        nbinsx=30,
        name="Permuted IC",
        marker_color=COLORS["slate"],
        opacity=0.8,
    )
)
fig.add_vline(
    x=observed_ic_mean,
    line_dash="solid",
    line_color=COLORS["negative"],
    line_width=2,
    annotation_text=f"Observed IC={observed_ic_mean:.4f}",
)
fig.update_layout(
    title="Block-permutation null for mean IC, against the observed value",
    xaxis_title="Mean IC (block-permuted, fold dates only)",
    yaxis_title="Count",
    height=350,
    showlegend=False,
)
show_plotly_with_alt(
    fig,
    alt=(
        "A histogram of the block-permutation null for mean IC over the fold dates, drawn "
        "in slate. It is a symmetric bell centred on zero, running out to roughly plus and "
        "minus 0.045. A solid red vertical line marks the observed mean IC, standing well "
        "clear of the right-hand edge of the null with no permuted value anywhere near it."
    ),
)
```

The claim that the block null is the honest one is checkable: rerun the same permutation
with a block of a single session, which is exactly the independent within-date shuffle,
and compare the two null widths on identical dates.

```python
naive_ics = []
for _ in range(n_permutations):
    ic_per_date = []
    for i, f_ranks in enumerate(factor_ranks_by_group):
        keys = rng.permutation(n_symbols)  # a fresh relabeling every date
        order = np.argsort(keys[symbols_by_group[i]], kind="stable")
        corr = np.corrcoef(f_ranks, return_ranks_by_group[i][order])[0, 1]
        if not np.isnan(corr):
            ic_per_date.append(corr)
    if ic_per_date:
        naive_ics.append(np.mean(ic_per_date))
naive_ics = np.array(naive_ics)

naive_p = (int(np.sum(naive_ics >= observed_ic_mean)) + 1) / (len(naive_ics) + 1)
print("=== Null width: independent within-date shuffle vs block permutation ===\n")
print(f"  independent shuffle  std {naive_ics.std():.4f}   one-sided p {naive_p:.4f}")
print(
    f"  block of {BLOCK_SESSIONS} sessions  std {permuted_ics.std():.4f}   one-sided p {p_value_perm:.4f}"
)
print(f"\n  ratio of null widths: {permuted_ics.std() / naive_ics.std():.1f}x, on identical dates")
```

**Interpretation**: over the fold window, the observed IC sits outside every block
permutation, so the factor ranked ETFs better than chance *in the years those folds
cover*. That is the only claim the test supports, and it is worth being precise about
why.

Read it against the full-sample IC printed under **Information Coefficient (IC)
Analysis**, which is indistinguishable from zero with a HAC $t$ well under one. These
are not opposing readings of one quantity - they are two different windows, as the
decomposition above showed: the entire full-sample-versus-fold gap was a period effect,
and the aggregation term rounded to nothing. The factor worked in the folds' two years
and does nothing over the full panel. The folds say something about those years, not
about the signal.

Two lessons generalize past this factor:

- **A window that flatters the signal is the default outcome, not a surprise.** The
  `min_train_pct` and fold count set above evaluate only the front of the panel and
  never score its last decade - the coverage share is printed with the decomposition. A
  fold scheme that leaves most of the sample unevaluated cannot support a claim about
  the sample.
- **A null must be as dependent as the data.** The comparison printed just above puts a
  number on it: shuffling assets independently each date gives a null several times
  narrower than blocking at the label horizon, on identical dates. Both nulls reject
  here, and at this permutation count both p-values sit on the floor, so the p-value is
  not what separates them - the width is. A null understated by that factor is the
  difference between rejecting and not rejecting for any factor whose IC lands near the
  boundary rather than far outside it, which is where most candidate factors land.

The rest of the evidence points the other way, and the quantile analysis already showed
it: the full-sample quintile spread is **negative** at every horizon, and monotonicity is
negative at every horizon too - strongly so at the short ones, where the ordering is
nearly perfect with its sign reversed. The bottom quintile out-earns the top.
Cross-sectional 21-day ETF momentum is a **reversal** signal over this panel. A single
favorable window does not overturn that; it illustrates how easily a fold scheme can
hide it.

## Factor Scorecard Output

Export a structured summary for downstream use, including both global and
fold-level statistics.

```python
# Build factor scorecard with fold-level stats
scorecard = {
    "factor_name": "momentum_21d",
    "n_assets": result.n_assets,
    "n_dates": result.n_dates,
    "horizons": {},
}

for period in PERIODS:
    period_key = f"{period}D"
    scorecard["horizons"][period_key] = {
        "ic_mean": round(result.ic.get(period_key, float("nan")), 4),
        "ic_std": round(np.std(result.ic_series.get(period_key, [])), 4)
        if result.ic_series.get(period_key)
        else None,
        "icir": round(result.ic_ir.get(period_key, float("nan")), 3),
        "t_stat": round(result.ic_t_stat.get(period_key, float("nan")), 2),
        "p_value": round(result.ic_p_value.get(period_key, float("nan")), 4),
        "spread": round(result.spread.get(period_key, float("nan")), 4),
        "monotonicity": round(result.monotonicity.get(period_key, float("nan")), 3),
    }

# Add fold-level statistics (21D only)
if fold_results:
    scorecard["fold_evaluation"] = {
        "n_folds": len(fold_results),
        "ic_mean": round(fold_ic_mean, 4),
        "ic_std": round(fold_ic_std, 4),
        "ic_min": round(fold_ic_min, 4),
        "ic_max": round(fold_ic_max, 4),
        "pct_positive_folds": round(pct_positive, 1),
        "fold_icir": round(fold_ic_mean / fold_ic_std, 3) if fold_ic_std > 0 else None,
    }

if result.turnover:
    scorecard["turnover"] = {k: round(v, 3) for k, v in result.turnover.items()}
if result.half_life:
    scorecard["half_life_periods"] = round(result.half_life, 1)

# Display scorecard
print("\n=== Factor Scorecard ===\n")
print(json.dumps(scorecard, indent=2))
```

## Binary Label Evaluation

When labels are binary (e.g., "positive return" vs "negative return"), we evaluate
using classification metrics rather than IC. The feature acts as a **score** that
separates positives from negatives.

**Key metrics:**
- **ROC AUC**: Area under ROC curve (threshold-free ranking metric)
- **PR AUC**: Area under Precision-Recall curve (better for imbalanced data)
- **Confusion matrix**: TP, FP, TN, FN at a chosen threshold

```python
# Positive = return > 0, negative = return <= 0.
binary_df = eval_df.with_columns(
    pl.when(pl.col("fwd_21d") > 0).then(1).otherwise(0).alias("binary_label")
)

# Use factor as score (higher = predict positive)
y_true = binary_df["binary_label"].to_numpy()
y_score = binary_df["factor"].to_numpy()

# Handle NaN values; sklearn accepts the native int32/float64 dtypes.
mask = ~(np.isnan(y_true) | np.isnan(y_score))
y_true = y_true[mask]
y_score = y_score[mask]

print(f"Binary evaluation: {len(y_true):,} samples")
print(f"Class balance: {y_true.mean():.1%} positive")
```

The ROC and precision-recall sweeps below are built by hand rather than taken from
sklearn's curve helpers, for an environment reason recorded in the code comment: the
library path fails under the pin this project carries. The manual sweep is deterministic
and agrees with the library to well inside plotting resolution.

```python
roc_auc = roc_auc_score(y_true, y_score)

# sklearn's _binary_clf_curve raises IndexError on state-dependent runs of this array
# under the sklearn pin econml imposes; the failure mode is uncharacterized.
_order = np.argsort(-y_score, kind="mergesort")
# int64 promotion guards np.cumsum from int32 overflow at 466k samples.
_yt_sorted = y_true[_order].astype(np.int64)
_score_sorted = y_score[_order]
_n_pos = int(_yt_sorted.sum())
_n_neg = len(_yt_sorted) - _n_pos
_tps_cum = np.cumsum(_yt_sorted)
_fps_cum = np.cumsum(1 - _yt_sorted)
fpr = np.concatenate([[0.0], _fps_cum / _n_neg])
tpr = np.concatenate([[0.0], _tps_cum / _n_pos])
# Thresholds aligned with fpr/tpr via the same _order index; no re-sort.
# Ties are not collapsed (unlike sklearn) - at 466k unique-or-near-unique
# float scores the AUC difference is sub-1e-6.
roc_thresholds = np.concatenate([[np.inf], _score_sorted])

# Precision-Recall (manual sweep omits sklearn's trailing (precision=1, recall=0)
# sentinel; the trapezoidal-AUC difference at 466k points is sub-1e-6)
_pred_pos = np.arange(1, len(_yt_sorted) + 1)
precision = _tps_cum / _pred_pos
recall = _tps_cum / _n_pos
pr_auc = auc(recall, precision)

print("\nThreshold-Free Metrics:")
print(f"  ROC AUC: {roc_auc:.3f}")
print(f"  PR AUC:  {pr_auc:.3f}")

# Interpretation
if roc_auc > 0.55:
    print("  Interpretation: AUC above 0.55 - score ranks positives above negatives")
elif roc_auc > 0.52:
    print(
        "  Interpretation: AUC in (0.52, 0.55] - small ranking signal; tradeability not evaluated here"
    )
else:
    print("  Interpretation: AUC at or near 0.5 - score does not separate the two classes")
```

```python
fig = make_subplots(rows=1, cols=2, subplot_titles=["ROC curve", "Precision-recall curve"])

# ROC curve
fig.add_trace(
    go.Scatter(
        x=fpr,
        y=tpr,
        mode="lines",
        name=f"ROC (AUC={roc_auc:.3f})",
        line=dict(color=COLORS["blue"], width=2),
    ),