用于超参数选择偏差评估的嵌套滚动验证

# Hyperparameter Selection and Validation Bias


# Hyperparameter Selection and Validation Bias

**Docker image**: `ml4t`

**Purpose**: develop a systematic approach to hyperparameter selection for
regularized models. The notebook answers two questions:

1. How does regularization affect Ridge performance on the ETF panel? (alpha
   grid analysis)
2. How much does HPO selection bias inflate single-loop CV estimates relative
   to nested CV? (Cawley and Talbot 2010)

**Learning objectives**

- Map performance stability across a wide regularization range
- Distinguish a robust alpha plateau from a noisy landscape
- Implement nested walk-forward CV with proper inner/outer purging and embargo
- Quantify the inflation from optimizing hyperparameters and reporting on the
  same data
- Apply the 3σ overfitting diagnostic to a hyperparameter search

**Book reference**: Section 11.2 - Regularized Regression (nested CV +
selection bias).

**Prerequisites**

- Ch7 21-day forward return labels at `case_studies/etfs/labels/fwd_ret_21d.parquet`
- Ch8 ETF features at `case_studies/etfs/features/financial.parquet`
- `02_regularization_paths` (coefficient paths context for Ridge)

**Downstream**: Ch12 extends the same nested-CV machinery to gradient boosting's
higher-dimensional hyperparameter space.

## 1. Setup and Imports

```python
"""Hyperparameter Selection and Validation Bias - nested CV for unbiased performance evaluation."""

import joblib
import matplotlib.pyplot as plt
import numpy as np
import optuna
import pandas as pd
import polars as pl
from matplotlib.colors import LinearSegmentedColormap
from ml4t.diagnostic.splitters import WalkForwardCV
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler

from utils.modeling import cross_sectional_ic_mean
from utils.paths import display_path, get_case_study_dir, get_chapter_dir
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, ml4t_diverging, show_with_alt

optuna.logging.set_verbosity(optuna.logging.WARNING)
```

```python
SEED = 42
RETRAIN = False
MAX_SYMBOLS = 0
N_SPLITS = 5
N_TRIALS = 20
CV_MAX_SPLITS = 0
CV_MAX_TRIALS = 0
ALPHA_GRID_POINTS = 23
TEST_SIZE = 200
ARTIFACT_TAG = ""
```

```python
# Configuration
RANDOM_SEED = SEED
set_global_seeds(SEED)
N_SPLITS_EFFECTIVE = N_SPLITS if CV_MAX_SPLITS <= 0 else min(N_SPLITS, CV_MAX_SPLITS)
N_TRIALS_EFFECTIVE = N_TRIALS if CV_MAX_TRIALS <= 0 else min(N_TRIALS, CV_MAX_TRIALS)
```

## 2. Load ETF Features

We use the pre-computed ETF features from Chapter 8 and labels from Chapter 7.

```python
CASE_DIR = get_case_study_dir("etfs")
FEATURES_PATH = CASE_DIR / "features" / "financial.parquet"
LABELS_PATH = CASE_DIR / "labels" / "fwd_ret_21d.parquet"

assert FEATURES_PATH.exists(), (
    f"Features not found: {FEATURES_PATH}\nRun Ch8 feature engineering first."
)
assert LABELS_PATH.exists(), f"Labels not found: {LABELS_PATH}\nRun Ch7 label engineering first."

features_df = pl.read_parquet(FEATURES_PATH).with_columns(pl.col("timestamp").cast(pl.Date))
labels_df = pl.read_parquet(LABELS_PATH).with_columns(pl.col("timestamp").cast(pl.Date))

# Join features and labels
TARGET_COL = "fwd_ret_21d"
ASSET_COL = "symbol"
dataset = features_df.join(labels_df, on=["timestamp", ASSET_COL], how="inner").drop_nulls(
    subset=[TARGET_COL]
)
dataset = dataset.sort(["timestamp", ASSET_COL])

if MAX_SYMBOLS > 0:
    assets = sorted(dataset[ASSET_COL].unique().to_list())[:MAX_SYMBOLS]
    dataset = dataset.filter(pl.col(ASSET_COL).is_in(assets))

print(f"Dataset: {len(dataset):,} rows")
print(f"Assets: {dataset[ASSET_COL].n_unique()}")
print(f"Date range: {dataset['timestamp'].min()} to {dataset['timestamp'].max()}")
```

## 3. Prepare Features and Target

A feature that could not be computed arrives as a null, and every row carrying
one is dropped. Zero is not a neutral filler for any of these columns: it is a
flat 252-day return, a zero volatility, and a mid-range oscillator reading, so
filling with it would put invented observations into the training set. The rows
this removes are the early history of a handful of symbols, where a long-window
feature has not had its warm-up period yet.

```python
meta_cols = {"timestamp", ASSET_COL}
label_cols = {c for c in dataset.columns if c.startswith("fwd_ret")}
feature_cols = [c for c in dataset.columns if c not in meta_cols and c not in label_cols]

df = dataset.select(feature_cols + [TARGET_COL, "timestamp", ASSET_COL]).with_columns(
    [
        pl.when(pl.col(c).is_nan() | pl.col(c).is_infinite())
        .then(None)
        .otherwise(pl.col(c))
        .alias(c)
        for c in feature_cols
    ]
)
rows_before = df.height
df = df.drop_nulls(subset=feature_cols)
print(f"Dropped {rows_before - df.height:,} rows with an incomplete feature vector")

X = df.select(feature_cols).to_numpy()
y = df[TARGET_COL].to_numpy()

# Create DataFrame with datetime index for WalkForwardCV
dates = df["timestamp"].to_list()
dates_arr = df["timestamp"].to_numpy()
symbols_arr = df[ASSET_COL].to_numpy()
df_for_cv = pd.DataFrame(X, columns=feature_cols)
df_for_cv.index = pd.DatetimeIndex(dates).tz_localize("UTC")


print(f"Final dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Target: {TARGET_COL}")
```

---
# Part A: Alpha Grid Analysis

Before optimizing hyperparameters, we should understand the **landscape** of
performance across regularization levels. This reveals:

- Is there a stable plateau where performance is robust?
- Or is performance noisy/sensitive to alpha choice?
- Should we use a single alpha or an ensemble?

This analysis informs our HPO strategy before we even run Optuna.

## 4. Define Alpha Grid and Evaluation Function

```python
ALPHAS = np.logspace(-2, 9, ALPHA_GRID_POINTS)
LABEL_HORIZON = 21  # sessions the label looks forward, and the purge each split needs
```

### Evaluate Alpha Grid

Walk-forward CV across the full alpha range, returning per-fold IC for each alpha.

```python
def evaluate_alpha_grid(
    X: np.ndarray,
    y: np.ndarray,
    df_cv: pd.DataFrame,
    alphas: np.ndarray,
    n_splits: int,
    test_size: int,
    label_horizon: int,
) -> pd.DataFrame:
    """
    Evaluate Ridge regression across a grid of alphas for each CV fold.

    Returns DataFrame with columns: fold, alpha, ic
    """
    results = []

    cv = WalkForwardCV(
        n_splits=n_splits,
        test_size=test_size,
        label_horizon=label_horizon,
        embargo_size=10,
        expanding=True,
    )

    for fold_idx, (train_idx, test_idx) in enumerate(cv.split(df_cv)):
        if len(train_idx) < 10:
            continue

        X_train_raw, X_test_raw = X[train_idx], X[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]

        # Scale inside fold
        scaler = StandardScaler()
        X_train = scaler.fit_transform(X_train_raw)
        X_test = scaler.transform(X_test_raw)

        dates_test = dates_arr[test_idx]
        symbols_test = symbols_arr[test_idx]

        # Evaluate each alpha
        for alpha in alphas:
            model = Ridge(alpha=alpha, random_state=RANDOM_SEED)
            model.fit(X_train, y_train)
            pred = model.predict(X_test)
            ic = cross_sectional_ic_mean(y_test, pred, dates_test, symbols_test)
            results.append({"fold": fold_idx + 1, "alpha": alpha, "ic": ic})

    return pd.DataFrame(results)
```

## 5. Run Alpha Grid Analysis

```python
MODELS_DIR = get_chapter_dir(11) / "models" / "04_nested_cv_hpo"
is_reduced_run = CV_MAX_SPLITS > 0 or CV_MAX_TRIALS > 0 or ALPHA_GRID_POINTS < 23 or MAX_SYMBOLS > 0
if is_reduced_run and not ARTIFACT_TAG:
    ARTIFACT_TAG = "_fast"

GRID_PATH = MODELS_DIR / f"grid_results{ARTIFACT_TAG}.joblib"
CV_PATH = MODELS_DIR / f"cv_results{ARTIFACT_TAG}.joblib"

NEED_GRID = RETRAIN or not GRID_PATH.exists()

if not NEED_GRID:
    grid_results = joblib.load(GRID_PATH)
    cached_alphas = np.sort(grid_results["alpha"].unique())
    if (
        cached_alphas.min() > ALPHAS.min()
        or cached_alphas.max() < ALPHAS.max()
        or len(cached_alphas) < len(ALPHAS)
    ):
        print("Cached grid uses a narrower alpha range; recomputing with expanded sweep.")
        NEED_GRID = True

if NEED_GRID:
    print("Evaluating alpha grid across CV folds...")
    grid_results = evaluate_alpha_grid(
        X, y, df_for_cv, ALPHAS, N_SPLITS_EFFECTIVE, TEST_SIZE, LABEL_HORIZON
    )
    MODELS_DIR.mkdir(parents=True, exist_ok=True)
    joblib.dump(grid_results, GRID_PATH)

# Pivot to get fold × alpha matrix
ic_matrix = grid_results.pivot(index="fold", columns="alpha", values="ic")
print(f"\nGrid shape: {ic_matrix.shape[0]} folds × {ic_matrix.shape[1]} alphas")
```

## 6. Visualize Alpha Landscape

The colour scale diverges around zero because zero is the meaningful point for
an IC: it separates a ranking that helped from one that hurt. The scale is made
symmetric so that equal distances either side of zero read as equal.

```python
ml4t_div = LinearSegmentedColormap.from_list("ml4t_div", ml4t_diverging())
fig, ax = plt.subplots(figsize=(12, 4))
ic_limit = float(np.nanmax(np.abs(ic_matrix.values)))
im = ax.imshow(
    ic_matrix.values,
    aspect="auto",
    cmap=ml4t_div,
    vmin=-ic_limit,
    vmax=ic_limit,
)
# Compact scientific labels (every other tick) keep the alpha grid legible.
tick_idx = list(range(0, len(ic_matrix.columns), 2))
ax.set_xticks(tick_idx)
ax.set_xticklabels(
    [f"{ic_matrix.columns[i]:.0e}" for i in tick_idx], rotation=45, ha="right", fontsize=8
)
ax.set_yticks(range(len(ic_matrix.index)))
ax.set_yticklabels([f"Fold {i}" for i in ic_matrix.index])
ax.set_xlabel(r"Alpha (log spacing, $10^{-2}$ to $10^{9}$)")
ax.set_title("IC by alpha and fold")
fig.colorbar(im, ax=ax, label="IC", shrink=0.8)
show_with_alt(
    fig,
    "Heatmap of information coefficient with one row per fold and one column per "
    "alpha, on a scale diverging around zero.",
)
```

```python
alpha_stats = (
    grid_results.groupby("alpha")
    .agg(
        mean_ic=("ic", "mean"),
        std_ic=("ic", "std"),
        min_ic=("ic", "min"),
        max_ic=("ic", "max"),
    )
    .reset_index()
)

# Plot mean IC with error bands
fig, ax = plt.subplots(figsize=(10, 5))
ax.fill_between(
    alpha_stats["alpha"],
    alpha_stats["mean_ic"] - alpha_stats["std_ic"],
    alpha_stats["mean_ic"] + alpha_stats["std_ic"],
    alpha=0.2,
    color=COLORS["amber"],
)
ax.semilogx(
    alpha_stats["alpha"],
    alpha_stats["mean_ic"],
    "o-",
    color=COLORS["amber"],
    lw=2,
    markersize=6,
)
ax.set_xlabel("Alpha (log scale)")
ax.set_ylabel("Information Coefficient")
ax.set_title("Mean IC against alpha, with the spread across folds")
show_with_alt(
    fig,
    "Mean information coefficient against alpha on a log axis, with a shaded band "
    "one standard deviation wide either side.",
)
```

## 7. Stability Analysis

```python
alpha_stats["cv"] = alpha_stats["std_ic"] / alpha_stats["mean_ic"].abs()

# Find the most stable region
best_alpha_idx = alpha_stats["mean_ic"].idxmax()
best_alpha = alpha_stats.loc[best_alpha_idx, "alpha"]
best_ic = alpha_stats.loc[best_alpha_idx, "mean_ic"]
best_std = alpha_stats.loc[best_alpha_idx, "std_ic"]

# Find stable plateau (alphas within 90% of best performance)
threshold = best_ic * 0.9 if best_ic > 0 else best_ic * 1.1
stable_alphas = alpha_stats[alpha_stats["mean_ic"] >= threshold]["alpha"]

# Assess stability
overall_cv = alpha_stats["cv"].median()
if overall_cv < 0.5:
    stability_label = "stable"
elif overall_cv < 1.0:
    stability_label = "moderate"
else:
    stability_label = "unstable"
```

```python
stability_df = pl.DataFrame(
    {
        "Metric": [
            "Best alpha",
            "Mean IC",
            "Std IC",
            "Plateau size",
            "Plateau range",
            "Median CV",
            "Stability",
        ],
        "Value": [
            f"{best_alpha:.3f}",
            f"{best_ic:.4f}",
            f"{best_std:.4f}",
            f"{len(stable_alphas)} alphas",
            f"[{stable_alphas.min():.3f}, {stable_alphas.max():.3f}]",
            f"{overall_cv:.2f}",
            stability_label,
        ],
    }
)
stability_df
```

**How to read the table.** Three of its rows describe the landscape and one
describes how much you should trust it.

The plateau size and range say whether the alpha that scores highest is
meaningfully better than its neighbours or merely first past the post. A
plateau spanning orders of magnitude means the choice within it hardly
matters, and picking the exact maximum is picking noise.

The median coefficient of variation is the ratio of fold-to-fold dispersion to
the size of the mean, taken across alphas. Above one, the IC at a given alpha
varies between folds by more than the mean IC itself is worth, and no ranking
of alphas computed from these folds is reliable.

That ratio is what makes Part B necessary. Where fold-level noise dwarfs the
mean, a search that takes the highest-scoring alpha per fold is taking whichever alpha
suited that fold's noise, and reporting its score as performance measures the
search rather than the model.

---
# Part B: Nested vs Single-Loop Cross-Validation

Now that we understand the alpha landscape, we address the **measurement problem**:
when we select hyperparameters and evaluate on the same data, we get optimistic
estimates due to selection bias.

**Key insight**: Single-loop CV (select best alpha, report that performance)
overestimates because we're measuring "how well did we fit the test set during HPO"
rather than "how well will we generalize."

A fold can end with nothing to select from. Every trial in it is pruned when no
validation date in that fold carries a defined IC - too few entities priced on
any one date, for instance - and Optuna then has no completed trial to report a
best parameter from. That is a fold the protocol could not evaluate, not a fold
that scored badly, so both runners record it as missing rather than dropping it:
the two protocols are compared fold by fold, and deleting a fold from one of
them would misalign the pair.

```python
def _selected_alpha(study) -> float:
    """The alpha of the best completed trial, or NaN when the fold produced none."""
    completed = study.get_trials(deepcopy=False, states=(optuna.trial.TrialState.COMPLETE,))
    return float(study.best_params["alpha"]) if completed else float("nan")
```

## 8. Single-Loop CV (BIASED)

This approach performs HPO and evaluation on the same data splits.
The selected hyperparameters are optimized for the test set, leading to
overly optimistic results.

```python
def _make_single_loop_objective(
    X_train: np.ndarray,
    y_train: np.ndarray,
    X_test: np.ndarray,
    y_test: np.ndarray,
    dates_test: np.ndarray,
    symbols_test: np.ndarray,
):
    """Create objective function that optimizes on test set (biased!)."""

    def objective(trial):
        alpha = trial.suggest_float("alpha", 0.01, 1e9, log=True)
        model = Ridge(alpha=alpha, random_state=RANDOM_SEED)
        model.fit(X_train, y_train)
        pred = model.predict(X_test)
        ic = cross_sectional_ic_mean(y_test, pred, dates_test, symbols_test)
        if np.isnan(ic):
            raise optuna.TrialPruned("no validation date carried a defined IC")
        return ic

    return objective
```

### Single-Loop CV Runner

Runs HPO and evaluation on the same data splits - the biased baseline.

```python
def single_loop_cv(
    X: np.ndarray,
    y: np.ndarray,
    df_cv: pd.DataFrame,
    n_splits: int = 5,
    n_trials: int = 15,
    label_horizon: int = 21,
    test_size: int = 200,
) -> dict:
    """
    Single-loop CV: HPO and evaluation on same splits (BIASED).

    This is the WRONG way - hyperparameters are selected based on
    the same test set used for final evaluation.
    """
    results = {"ic": [], "best_alpha": []}

    cv = WalkForwardCV(
        n_splits=n_splits,
        test_size=test_size,
        label_horizon=label_horizon,
        embargo_size=10,
        expanding=True,
    )

    for fold_idx, (train_idx, test_idx) in enumerate(cv.split(df_cv)):
        if len(train_idx) < 10:
            continue

        X_train_raw, X_test_raw = X[train_idx], X[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        dates_test = dates_arr[test_idx]
        symbols_test = symbols_arr[test_idx]

        # Scale inside fold
        scaler = StandardScaler()
        X_train = scaler.fit_transform(X_train_raw)
        X_test = scaler.transform(X_test_raw)

        # HPO on test set (this is the bias!)
        objective = _make_single_loop_objective(
            X_train, y_train, X_test, y_test, dates_test, symbols_test
        )
        study = optuna.create_study(
            direction="maximize",
            sampler=optuna.samplers.TPESampler(seed=RANDOM_SEED),
        )
        study.optimize(objective, n_trials=n_trials, show_progress_bar=False)

        # Evaluate with best params on SAME test set
        best_alpha = _selected_alpha(study)
        if np.isnan(best_alpha):
            results["ic"].append(float("nan"))
            results["best_alpha"].append(best_alpha)
            continue
        model = Ridge(alpha=best_alpha, random_state=RANDOM_SEED)
        model.fit(X_train, y_train)
        pred = model.predict(X_test)
        ic = cross_sectional_ic_mean(y_test, pred, dates_test, symbols_test)

        results["ic"].append(ic)
        results["best_alpha"].append(best_alpha)

    return results
```

## 9. Nested CV (UNBIASED)

This approach uses an **inner loop** for HPO and **outer loop** for evaluation.
The test set is never seen during hyperparameter selection.

**Both loops use WalkForwardCV** to maintain proper temporal structure
with purging and embargo.

```python
def _make_nested_objective(
    X_outer_train: np.ndarray,
    y_outer_train: np.ndarray,
    dates_outer_train: np.ndarray,
    symbols_outer_train: np.ndarray,
    df_inner: pd.DataFrame,
    n_inner: int,
    label_horizon: int,
):
    """
    Create objective function for nested CV inner loop.

    Uses WalkForwardCV for inner loop to maintain proper temporal
    structure with purging and embargo.
    """
    # WalkForwardCV sizes its windows in sessions; this is a panel with many rows
    # per session, so the inner test size is counted in sessions too.
    inner_sessions = int(pd.DatetimeIndex(df_inner.index).nunique())
    inner_cv = WalkForwardCV(
        n_splits=n_inner,
        test_size=max(label_horizon, inner_sessions // (n_inner + 2)),
        label_horizon=label_horizon,
        embargo_size=label_horizon,  # Inner embargo must also respect label horizon
        expanding=True,
    )

    def objective(trial):
        alpha = trial.suggest_float("alpha", 0.01, 1e9, log=True)
        inner_scores = []

        for inner_train_idx, inner_val_idx in inner_cv.split(df_inner):
            if len(inner_train_idx) < 10:
                continue

            X_inner_train = X_outer_train[inner_train_idx]
            X_inner_val = X_outer_train[inner_val_idx]
            y_inner_train = y_outer_train[inner_train_idx]
            y_inner_val = y_outer_train[inner_val_idx]
            dates_inner_val = dates_outer_train[inner_val_idx]
            symbols_inner_val = symbols_outer_train[inner_val_idx]

            # Scale inside inner fold
            inner_scaler = StandardScaler()
            X_inner_train_scaled = inner_scaler.fit_transform(X_inner_train)
            X_inner_val_scaled = inner_scaler.transform(X_inner_val)

            model = Ridge(alpha=alpha, random_state=RANDOM_SEED)
            model.fit(X_inner_train_scaled, y_inner_train)
            pred = model.predict(X_inner_val_scaled)
            ic = cross_sectional_ic_mean(y_inner_val, pred, dates_inner_val, symbols_inner_val)
            if not np.isnan(ic):
                inner_scores.append(ic)

        if not inner_scores:
            raise optuna.TrialPruned("no inner fold carried a defined IC")
        return float(np.mean(inner_scores))

    return objective
```

### Nested CV Runner

Uses an inner WalkForwardCV loop for HPO, keeping the outer test set
completely hidden from hyperparameter selection.

```python
def nested_cv(
    X: np.ndarray,
    y: np.ndarray,
    df_cv: pd.DataFrame,
    n_outer: int = 5,
    n_inner: int = 3,
    n_trials: int = 15,
    label_horizon: int = 21,
    test_size: int = 200,
) -> dict:
    """
    Nested CV: HPO isolated from evaluation (UNBIASED).

    Outer loop: WalkForwardCV for temporal validation
    Inner loop: WalkForwardCV for HPO (within training data only)

    This is the CORRECT approach - hyperparameters are selected using only
    training data, and final evaluation is on truly held-out test data.
    """
    results = {"ic": [], "best_alpha": []}

    outer_cv = WalkForwardCV(
        n_splits=n_outer,
        test_size=test_size,
        label_horizon=label_horizon,
        embargo_size=10,
        expanding=True,
    )

    for fold_idx, (outer_train_idx, outer_test_idx) in enumerate(outer_cv.split(df_cv)):
        if len(outer_train_idx) < 10:
            continue

        X_outer_train_raw = X[outer_train_idx]
        X_outer_test_raw = X[outer_test_idx]
        y_outer_train = y[outer_train_idx]
        y_outer_test = y[outer_test_idx]
        dates_outer_train = dates_arr[outer_train_idx]
        symbols_outer_train = symbols_arr[outer_train_idx]
        dates_outer_test = dates_arr[outer_test_idx]
        symbols_outer_test = symbols_arr[outer_test_idx]

        # Create inner DataFrame for WalkForwardCV
        df_inner = df_cv.iloc[outer_train_idx].copy()

        # Inner loop: HPO on outer training data ONLY
        objective = _make_nested_objective(
            X_outer_train_raw,
            y_outer_train,
            dates_outer_train,
            symbols_outer_train,
            df_inner,
            n_inner,
            label_horizon,
        )

        study = optuna.create_study(
            direction="maximize",
            sampler=optuna.samplers.TPESampler(seed=RANDOM_SEED),
        )
        study.optimize(objective, n_trials=n_trials, show_progress_bar=False)

        # Outer evaluation with best params (UNBIASED)
        best_alpha = _selected_alpha(study)
        if np.isnan(best_alpha):
            results["ic"].append(float("nan"))
            results["best_alpha"].append(best_alpha)
            continue

        # Scale outer fold for final evaluation
        outer_scaler = StandardScaler()
        X_outer_train = outer_scaler.fit_transform(X_outer_train_raw)
        X_outer_test = outer_scaler.transform(X_outer_test_raw)

        model = Ridge(alpha=best_alpha, random_state=RANDOM_SEED)
        model.fit(X_outer_train, y_outer_train)
        pred = model.predict(X_outer_test)
        ic = cross_sectional_ic_mean(y_outer_test, pred, dates_outer_test, symbols_outer_test)

        results["ic"].append(ic)
        results["best_alpha"].append(best_alpha)

    return results
```

## 10. Run Comparison

```python
N_INNER = 3

NEED_CV = RETRAIN or not CV_PATH.exists()

if NEED_CV:
    print("=" * 60)
    print("Running Single-Loop CV (BIASED)...")
    print("=" * 60)
    single_results = single_loop_cv(
        X,
        y,
        df_for_cv,
        n_splits=N_SPLITS_EFFECTIVE,
        n_trials=N_TRIALS_EFFECTIVE,
        label_horizon=LABEL_HORIZON,
        test_size=TEST_SIZE,
    )

    print("\n" + "=" * 60)
    print("Running Nested CV (UNBIASED)...")
    print("=" * 60)
    nested_results = nested_cv(
        X,
        y,
        df_for_cv,
        n_outer=N_SPLITS_EFFECTIVE,
        n_inner=N_INNER,
        n_trials=N_TRIALS_EFFECTIVE,
        label_horizon=LABEL_HORIZON,
        test_size=TEST_SIZE,
    )

    MODELS_DIR.mkdir(parents=True, exist_ok=True)
    joblib.dump({"single_results": single_results, "nested_results": nested_results}, CV_PATH)
    print(f"Saved results to {display_path(CV_PATH)}")
```

### Load Cached Results

If CV results are already cached, load them and verify completeness.

```python
if not NEED_CV:
    _cached = joblib.load(CV_PATH)
    single_results = _cached["single_results"]
    nested_results = _cached["nested_results"]
    del _cached
    if (
        len(single_results["ic"]) < N_SPLITS_EFFECTIVE
        or len(nested_results["ic"]) < N_SPLITS_EFFECTIVE
    ):
        print("Cached CV results appear incomplete; recomputing.")
        single_results = single_loop_cv(
            X,
            y,
            df_for_cv,
            n_splits=N_SPLITS_EFFECTIVE,
            n_trials=N_TRIALS_EFFECTIVE,
            label_horizon=LABEL_HORIZON,
            test_size=TEST_SIZE,
        )
        nested_results = nested_cv(
            X,
            y,
            df_for_cv,
            n_outer=N_SPLITS_EFFECTIVE,
            n_inner=N_INNER,
            n_trials=N_TRIALS_EFFECTIVE,
            label_horizon=LABEL_HORIZON,
            test_size=TEST_SIZE,
        )
        MODELS_DIR.mkdir(parents=True, exist_ok=True)
        joblib.dump({"single_results": single_results, "nested_results": nested_results}, CV_PATH)
```

## 11. Results Analysis

```python
def _alpha_range(results: dict) -> str:
    """The span of the alphas actually selected, or a note when none were."""
    chosen = np.asarray(results["best_alpha"], dtype=float)
    if not np.isfinite(chosen).any():
        return "no fold selected one"
    return f"[{np.nanmin(chosen):.2f}, {np.nanmax(chosen):.2f}]"


single_mean_ic = np.nanmean(single_results["ic"])
single_std_ic = np.nanstd(single_results["ic"])
nested_mean_ic = np.nanmean(nested_results["ic"])
nested_std_ic = np.nanstd(nested_results["ic"])

# Calculate inflation
if nested_mean_ic != 0:
    inflation_pct = (single_mean_ic - nested_mean_ic) / abs(nested_mean_ic) * 100
else:
    inflation_pct = 0.0

comparison_df = pl.DataFrame(
    {
        "Method": ["Single-Loop (Biased)", "Nested (Unbiased)", "Inflation"],
        "Mean IC": [f"{single_mean_ic:.4f}", f"{nested_mean_ic:.4f}", f"{inflation_pct:+.1f}%"],
        "Std IC": [f"{single_std_ic:.4f}", f"{nested_std_ic:.4f}", ""],
        "Alpha Range": [_alpha_range(single_results), _alpha_range(nested_results), ""],
    }
)
comparison_df
```

**How to read the table.** The two protocols ran on identical outer splits and
differ in one respect: where the alpha came from. Single-loop CV chose it by
maximizing IC on the very rows it then reports, so its number answers "how
well can this model be made to fit this test fold". Nested CV chose it inside
the training window and never showed the outer test fold to the search, so its
number answers "how well does the procedure generalize". The gap between them
is the selection bias, and it is a property of the search, not of Ridge.

The alpha column is the more direct evidence. Selections that jump across
orders of magnitude from fold to fold mean the search is tracking each fold's
noise rather than a stable property of the data. Selections that cluster mean
the inner-loop average is finding something the folds agree on. Compare the
spread of the two columns before reading either mean IC.

### Validation Overfitting Diagnostic (3σ Heuristic)

The text proposes a concrete diagnostic: if the top configuration's IC exceeds the
median by more than 3× the inter-configuration standard deviation, the search is
likely overfitting to validation noise. We apply this to the alpha grid from Part A,
where each alpha's mean IC across folds serves as a "trial" value.

```python
trial_ics = alpha_stats["mean_ic"].values
best_trial_ic = trial_ics.max()
median_trial_ic = np.median(trial_ics)
std_trial_ic = np.std(trial_ics)

gap_sigma = (best_trial_ic - median_trial_ic) / std_trial_ic if std_trial_ic > 0 else 0.0
overfitting_flag = "WARNING: likely overfitting" if gap_sigma > 3.0 else "OK"

diagnostic_df = pl.DataFrame(
    {
        "Metric": ["Best IC", "Median IC", "Std IC", "Gap (σ)", "Overfit flag"],
        "Value": [
            f"{best_trial_ic:.4f}",
            f"{median_trial_ic:.4f}",
            f"{std_trial_ic:.4f}",
            f"{gap_sigma:.1f}σ",
            overfitting_flag,
        ],
    }
)
diagnostic_df
```

This diagnostic and the stability classification above measure dispersion in
two different directions, and they can disagree without either being wrong.
Stability looks *within* an alpha, across folds: it asks whether the same
setting scores consistently on different periods. The three-sigma gap looks
*across* alphas, at the fold-averaged scores: it asks whether the top setting
stands out from the field.

A landscape can be smooth across alphas and noisy across folds at the same
time, and that combination is the one to watch for. It means the alpha-by-alpha
curve looks orderly enough to invite a confident choice, while the evidence
behind each point on it is thin.

## 12. Visualization

```python
folds = list(range(1, len(single_results["ic"]) + 1))
x = np.arange(len(folds))
bar_w = 0.35

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# (a) IC comparison
axes[0].bar(
    x - bar_w / 2, single_results["ic"], bar_w, label="Single-Loop (Biased)", color=COLORS["blue"]
)
axes[0].bar(
    x + bar_w / 2, nested_results["ic"], bar_w, label="Nested (Unbiased)", color=COLORS["amber"]
)
axes[0].set_xticks(x)
axes[0].set_xticklabels([f"Fold {f}" for f in folds])
axes[0].set_ylabel("Information Coefficient")
axes[0].set_title("IC by fold, single-loop against nested")
axes[0].legend(frameon=False)

# (b) Alpha selection
axes[1].semilogy(
    folds, single_results["best_alpha"], "o-", color=COLORS["blue"], label="Single-Loop"
)
axes[1].semilogy(folds, nested_results["best_alpha"], "s-", color=COLORS["amber"], label="Nested")
axes[1].set_xlabel("Fold")
axes[1].set_ylabel("Best Alpha (log)")
axes[1].set_title("Alpha chosen in each fold, by procedure")
axes[1].legend(frameon=False)

fig.suptitle("Single-loop and nested selection, fold by fold", fontsize=13)
show_with_alt(
    fig,
    "Two panels: paired bars of information coefficient per fold for the single-loop "
    "and nested protocols, and the alpha each protocol selected per fold on a log axis.",
)
```

The left panel is the cost of selecting and scoring on the same data. Where the two bars
differ, the single-loop score is the more favourable of the pair - less negative where the
fold is negative, higher where it is positive - and the nested score is what the procedure
is worth when the alpha it picked was chosen without seeing the fold it is scored on. Where
the bars are equal, the two procedures happened to select the same alpha.

The gap has one cause and it is not in the right panel: the single-loop procedure picks its
alpha using the same fold it is then scored on, so its score is partly a report on how well
the choice fitted that fold. Nested selection picks on an inner split and scores on data
the choice never saw, which is what makes its number an estimate of the procedure rather
than of one fold.

The right panel is worth reading alongside that and does not establish it. It shows the
single-loop choice swinging across orders of magnitude while the nested choice moves
smoothly, which is a fact about these five folds. Choosing a different alpha for different
training data is what an adaptive procedure does and is not itself a defect - nested CV
evaluates that procedure honestly whether its choice is stable or not, and a procedure that
picked the same alpha every time would be no less exposed to selection bias if it picked it
on the evaluation fold.

## 13. Key Takeaways

1. **Report the number the protocol earns.** A hyperparameter chosen on a set
   of rows makes any score computed on those same rows a statement about the
   search. Nested CV separates the two by keeping the outer test fold out of
   the selection entirely, and the difference between the two mean ICs above
   is what the separation costs on paper and buys in honesty.

2. **The spread of the selected values is the tell.** Where a search chases
   noise, its choice moves from fold to fold with nothing in the data driving
   it. Plotting the selected alpha per fold, as panel (b) does, diagnoses that
   faster than comparing scores.

3. **Two dispersion diagnostics, in different directions.** Fold-to-fold
   variation at a fixed alpha and alpha-to-alpha variation of the fold average
   answer different questions, and a landscape that is smooth in one can be
   dominated by noise in the other.

4. **A wide grid before an expensive search is cheap insurance.** Sweeping
   alpha across many orders of magnitude costs one pass and tells you whether
   the landscape has a plateau worth searching in at all.

5. **Both loops need the purge.** The inner split has the same overlapping
   label as the outer one, so it needs the same gap; an inner loop that
   silently produces no folds leaves the search with nothing to optimize and
   hands back whatever its sampler tried first.

**Reference**: Cawley and Talbot (2010) document the same bias mechanism
across multiple datasets and model classes.

```python
print("Hyperparameter selection analysis complete")
print(
    f"  Alpha grid: {len(ALPHAS)} levels evaluated, mean IC range "
    f"{alpha_stats['mean_ic'].min():.4f} to {alpha_stats['mean_ic'].max():.4f}"
)
print(f"  Single-loop mean IC: {single_mean_ic:+.4f}  |  Nested mean IC: {nested_mean_ic:+.4f}")
print(f"  HPO inflation:       {inflation_pct:+.1f}%")
print(f"  Stability label:     {stability_label}   |  Cross-alpha gap: {gap_sigma:.1f}σ")
```
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)
![notebook output](figures/p1_3.png)