| """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 |
| from scripts.run_volatility_experiment import fetch_klines |
|
|
| |
| 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) |
| 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 |
|
|
| |
| 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() |
|
|
