Blame for scripts/run_volatility_selective_forecasting.py

"""Selective forecasting for the volatility-regime model: accuracy vs coverage.

The model answers every hour, including the hours where its score sits on the
decision boundary and it plainly does not know. Abstaining on those and
reporting accuracy only on the hours it does qualify is the one honest route to
a higher headline figure -- and it is the method specified in
notebooks/BTC_1Y_Regime_MetaLabel_Selective_Forecasting.ipynb, which has never
been executed.

It also fills a real hole in the application. ForecastContextAgent resolves
`validated_reliability`, `meta_reliability` and `reliability`, and nothing in the
codebase ever writes them, so every forecast falls through to UNKNOWN and signal
coverage is zero. This produces the thresholds that gate would need.

Two rules keep the result honest.

First, the comparison must be on identical hours. Abstaining on hard cases makes
*any* rule look better on what remains, persistence included, so measuring a
reduced-coverage model against the full-sample 65.61% baseline would be
meaningless. Persistence is therefore rescored on exactly the hours the model
qualifies.

Second, the operating point is chosen on a validation slice carved out of each
training fold, never on the test fold. The full curve is reported as well, so
nothing depends on one chosen point.

The deployed artifact cannot be unpickled on this machine -- its scikit-learn
1.6.1 pickle needs a module 1.9.1 renamed -- so models are refit per fold from
the same configuration. That configuration reproduces the artifact closely:
AUC 0.7875 here against 0.7895 recorded in its metadata.

Usage:
    python scripts/run_volatility_selective_forecasting.py
    python scripts/run_volatility_selective_forecasting.py --folds 5 --target-coverage 0.6
"""
from __future__ import annotations

import argparse
import json
import sys
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import matthews_corrcoef, roc_auc_score

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))
warnings.filterwarnings("ignore")

from src.volatility import add_volatility_target, build_volatility_features, purge_split  # noqa: E402
from scripts.run_volatility_experiment import fetch_klines                                 # noqa: E402

# identical to the deployed configuration
MODEL = dict(max_iter=400, learning_rate=0.06, max_depth=6,
             l2_regularization=1.0, random_state=42)
COVERAGES = [1.00, 0.90, 0.80, 0.70, 0.60, 0.50, 0.40, 0.30, 0.20]


def fit_fold(train: pd.DataFrame, features: list[str], horizon: int):
    """Fit, and choose the decision threshold on a tail of the training fold."""
    split = int(len(train) * 0.85)
    fit_part, validation = train.iloc[: split - horizon], train.iloc[split:]
    probe = HistGradientBoostingClassifier(**MODEL).fit(
        fit_part[features], fit_part["target_volatility_regime"].astype(int))
    val_probs = probe.predict_proba(validation[features])[:, 1]
    val_truth = validation["target_volatility_regime"].astype(int).to_numpy()
    threshold = max(np.linspace(0.30, 0.70, 41),
                    key=lambda t: matthews_corrcoef(val_truth, (val_probs >= t).astype(int)))
    model = HistGradientBoostingClassifier(**MODEL).fit(
        train[features], train["target_volatility_regime"].astype(int))
    return model, float(threshold), val_probs, val_truth, validation


def selective_curve(probs, truth, persistence, threshold, coverages) -> list[dict]:
    """Accuracy at each coverage, with persistence rescored on the same hours.

    Confidence is distance from the decision threshold. Keeping the most
    confident share means abstaining on the boundary, which is where the model's
    errors concentrate.
    """
    confidence = np.abs(probs - threshold)
    order = np.argsort(-confidence)             # most confident first
    predictions = (probs >= threshold).astype(int)
    rows = []
    for coverage in coverages:
        keep = max(1, int(round(len(order) * coverage)))
        idx = order[:keep]
        model_hits = int((predictions[idx] == truth[idx]).sum())
        persistence_hits = int((persistence[idx] == truth[idx]).sum())
        rows.append(dict(
            coverage=coverage, n=keep,
            accuracy=model_hits / keep,
            persistence=persistence_hits / keep,
            margin=(model_hits - persistence_hits) / keep,
            min_confidence=float(confidence[idx].min()),
            base_rate=float(truth[idx].mean())))
    return rows


def validation_operating_point(val_probs, val_truth, threshold,
                              target_coverage: float) -> float:
    """The confidence cut-off that yields `target_coverage` on validation only.

    Chosen before the test fold is touched, so the figure it produces there is
    genuinely out of sample.
    """
    confidence = np.abs(val_probs - threshold)
    return float(np.quantile(confidence, 1.0 - target_coverage))


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--folds", type=int, default=5)
    ap.add_argument("--horizon", type=int, default=24)
    ap.add_argument("--target-coverage", type=float, default=0.60)
    ap.add_argument("--cache", default="outputs/volatility/btc_1h.csv")
    ap.add_argument("--start", default="2019-01-01")
    args = ap.parse_args()

    raw = fetch_klines(args.start, ROOT / args.cache)
    featured, features = build_volatility_features(raw)
    featured = add_volatility_target(featured, args.horizon)
    data = featured.dropna(
        subset=features + ["forward_vol", "threshold_causal", "target_volatility_regime"]
    ).reset_index(drop=True)
    print(f"hourly rows {len(data):,}   {data.timestamp.min().date()} -> {data.timestamp.max().date()}")
    print(f"operating point chosen on validation for {args.target_coverage:.0%} coverage\n")

    total = len(data)
    first = int(total * 0.40)
    step = (total - first) // args.folds
    per_fold, chosen_rows = [], []

    for fold in range(args.folds):
        train_end = first + fold * step
        train, test = purge_split(data, train_end, min(train_end + step, total), args.horizon)
        if len(test) < 500:
            continue
        model, threshold, val_probs, val_truth, _ = fit_fold(train, features, args.horizon)
        probs = model.predict_proba(test[features])[:, 1]
        truth = test["target_volatility_regime"].astype(int).to_numpy()
        persistence = (test["vol_24"] > test["threshold_causal"]).astype(int).to_numpy()

        curve = selective_curve(probs, truth, persistence, threshold, COVERAGES)
        for row in curve:
            row.update(fold=fold + 1, start=str(test.timestamp.min().date()),
                       end=str(test.timestamp.max().date()))
        per_fold += curve

        # the pre-committed operating point
        cut = validation_operating_point(val_probs, val_truth, threshold, args.target_coverage)
        keep = np.abs(probs - threshold) >= cut
        if keep.sum() >= 50:
            predictions = (probs >= threshold).astype(int)
            chosen_rows.append(dict(
                fold=fold + 1, cut=cut, n=int(keep.sum()),
                realised_coverage=float(keep.mean()),
                accuracy=float((predictions[keep] == truth[keep]).mean()),
                persistence=float((persistence[keep] == truth[keep]).mean()),
                full_accuracy=float((predictions == truth).mean()),
                full_persistence=float((persistence == truth).mean()),
                auc=float(roc_auc_score(truth, probs)) if len(set(truth)) > 1 else None))
        print(f"  fold {fold + 1}  {test.timestamp.min().date()} -> {test.timestamp.max().date()}"
              f"   n {len(test):,}   threshold {threshold:.2f}   AUC "
              f"{roc_auc_score(truth, probs):.4f}")

    df = pd.DataFrame(per_fold)
    agg = (df.groupby("coverage")
             .agg(folds=("fold", "count"), n=("n", "sum"), accuracy=("accuracy", "mean"),
                  persistence=("persistence", "mean"), margin=("margin", "mean"),
                  beat=("margin", lambda s: int((s > 0).sum())))
             .reset_index().sort_values("coverage", ascending=False))

    print("\n" + "=" * 96)
    print("ACCURACY vs COVERAGE   (persistence rescored on exactly the hours the model keeps)")
    print("=" * 96)
    print(f"{'coverage':>9}{'hours kept':>12}{'accuracy':>11}{'persistence':>13}"
          f"{'margin':>10}{'beats':>8}")
    for _, r in agg.iterrows():
        print(f"{r['coverage'] * 100:>8.0f}%{int(r['n']):>12,}{r['accuracy'] * 100:>10.2f}%"
              f"{r['persistence'] * 100:>12.2f}%{r['margin'] * 100:>+9.2f}pp"
              f"{int(r['beat'])}/{int(r['folds']):>3}")

    print("\n  Full coverage reproduces the deployed figure. Reading down the table shows")
    print("  what abstention buys -- and the margin column shows whether it buys anything")
    print("  the baseline does not also get for free.")

    if chosen_rows:
        chosen = pd.DataFrame(chosen_rows)
        print("\n" + "=" * 96)
        print(f"PRE-COMMITTED OPERATING POINT  (cut-off set on validation for "
              f"{args.target_coverage:.0%} coverage, applied to the test fold)")
        print("=" * 96)
        print(f"{'fold':>5}{'coverage':>10}{'hours':>8}{'accuracy':>11}{'persistence':>13}"
              f"{'margin':>10}{'all-hours acc':>15}")
        for _, r in chosen.iterrows():
            print(f"{int(r['fold']):>5}{r['realised_coverage'] * 100:>9.1f}%{int(r['n']):>8,}"
                  f"{r['accuracy'] * 100:>10.2f}%{r['persistence'] * 100:>12.2f}%"
                  f"{(r['accuracy'] - r['persistence']) * 100:>+9.2f}pp"
                  f"{r['full_accuracy'] * 100:>14.2f}%")
        gain = (chosen["accuracy"] - chosen["full_accuracy"]).to_numpy()
        edge = (chosen["accuracy"] - chosen["persistence"]).to_numpy()
        t_gain, p_gain = stats.ttest_1samp(gain, 0) if len(gain) > 2 else (np.nan, np.nan)
        t_edge, p_edge = stats.ttest_1samp(edge, 0) if len(edge) > 2 else (np.nan, np.nan)
        print(f"\n  mean coverage          {chosen['realised_coverage'].mean() * 100:.1f}%")
        print(f"  mean accuracy          {chosen['accuracy'].mean() * 100:.2f}%  "
              f"(all hours {chosen['full_accuracy'].mean() * 100:.2f}%)")
        print(f"  accuracy gained        {gain.mean() * 100:+.2f}pp   p = {p_gain:.3f}")
        print(f"  edge over persistence  {edge.mean() * 100:+.2f}pp   p = {p_edge:.3f}")
        print(f"  beat persistence in    {int((edge > 0).sum())}/{len(edge)} folds")

    out = ROOT / "outputs" / "volatility_selective_forecasting.json"
    out.write_text(json.dumps({
        "target_coverage": args.target_coverage,
        "curve": per_fold,
        "curve_mean": agg.to_dict("records"),
        "operating_point": chosen_rows,
    }, indent=2, default=str), encoding="utf-8")
    print(f"\nwritten: {out.relative_to(ROOT)}")


if __name__ == "__main__":
    main()