"""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()