| """Four candidate upgrades to the live price-range forecaster (model B). |
| |
| Same walk-forward as backtest_improvements.py (refit every quarter on the two |
| years ending 72h before it, test the quarter, Jan 2023 - Sep 2026) and the same |
| scoring, so every row is comparable with the live model: |
| |
| B (live) past-week volatility x hour-of-week x Fed/CPI calendar, |
| empirical |z| quantiles |
| HAR blend the per-hour variance forecast is a fitted mix of the last |
| day, week and month (Corsi 2009), in logs, per horizon |
| Range estimator the past week's Parkinson high/low variance replaces the |
| close-to-close one: the same hours, less noise |
| HAR + range both |
| Asymmetric B's scale with separate lower and upper quantiles |
| Adaptive (ACI) B, with each level's width nudged after every forecast |
| whose outcome is already known (Gibbs & Candes 2021) |
| |
| Seasonality is divided out before any variance is averaged, and multiplied |
| back over the forecast window, exactly as in B. Diebold-Mariano tests compare |
| each row's 80% interval score with B's (Newey-West, h-1 lags). |
| |
| Usage: |
| python scripts/backtest_range_upgrades.py |
| """ |
| from __future__ import annotations |
|
|
| import json |
| import math |
| import sys |
| from pathlib import Path |
|
|
| import numpy as np |
| import pandas as pd |
|
|
| sys.path.insert(0, str(Path(__file__).resolve().parent)) |
| from backtest_improvements import GAP, STEP, TEST_START, WINDOW, fit_factors, prepare |
| from backtest_price_range import HORIZONS, LEVELS |
|
|
| ROOT = Path(__file__).resolve().parents[1] |
| RULES = ["B (live)", "HAR blend", "Range estimator", "HAR + range", "Asymmetric", "Adaptive (ACI)"] |
| ACI_GAMMA = 0.005 |
|
|
|
|
| def high_low(index: pd.DatetimeIndex) -> pd.Series: |
| b = pd.read_csv(ROOT / "outputs" / "volatility" / "btc_1h.csv", usecols=["timestamp", "high", "low"]) |
| b["timestamp"] = pd.to_datetime(b["timestamp"], utc=True) |
| b = b.drop_duplicates("timestamp").set_index("timestamp").reindex(index) |
| return (np.log(b["high"] / b["low"]) ** 2) / (4 * math.log(2)) |
|
|
|
|
| def fwd_sum(x: pd.Series, h: int) -> np.ndarray: |
| return x.rolling(h).sum().shift(-h).to_numpy() |
|
|
|
|
| def log_ols(y, X, train): |
| """Fit log(y) on log(X) over the training rows; return a predictor.""" |
| ok = train & np.isfinite(y) & (y > 0) & np.isfinite(X).all(1) & (X > 0).all(1) |
| A = np.column_stack([np.ones(ok.sum()), np.log(X[ok])]) |
| beta, *_ = np.linalg.lstsq(A, np.log(y[ok]), rcond=None) |
| smear = float(np.mean(np.exp(np.log(y[ok]) - A @ beta))) |
|
|
| def predict(Xn): |
| with np.errstate(divide="ignore", invalid="ignore"): |
| return np.exp(np.column_stack([np.ones(len(Xn)), np.log(Xn)]) @ beta) * smear |
| return predict |
|
|
|
|
| def scores(z_lo, z_hi, scale, move, lo_q, hi_q): |
| """Coverage flags per level and the 80% interval score, as % of price.""" |
| out = {} |
| for l in LEVELS: |
| out[f"in{int(l * 100)}"] = (move >= lo_q[l] * scale) & (move <= hi_q[l] * scale) |
| lo, hi = np.exp(lo_q[0.8] * scale) - 1, np.exp(hi_q[0.8] * scale) - 1 |
| y = np.exp(move) - 1 |
| out["width"] = hi - lo |
| out["score"] = (hi - lo) + 10 * (np.clip(lo - y, 0, None) + np.clip(y - hi, 0, None)) |
| return out |
|
|
|
|
| def dm(a: np.ndarray, b: np.ndarray, h: int) -> dict: |
| d = a - b |
| u = d - d.mean() |
| var = u @ u / len(d) |
| for lag in range(1, h): |
| var += 2 * (1 - lag / h) * (u[lag:] @ u[:-lag]) / len(d) |
| t = float(d.mean() / math.sqrt(var / len(d))) |
| return {"t": t, "p": float(math.erfc(abs(t) / math.sqrt(2)))} |
|
|
|
|
| def main() -> None: |
| base = prepare() |
| idx = base["ret"].index |
| ret2 = base["ret"] ** 2 |
| park = high_low(idx) |
| logp = base["logp"] |
| moves = {h: (logp.shift(-h) - logp).to_numpy() for h in HORIZONS} |
| rows = {h: {r: [] for r in RULES} for h in HORIZONS} |
| meta = {h: [] for h in HORIZONS} |
|
|
| start = pd.Timestamp(TEST_START, tz="UTC") |
| while start < idx.max(): |
| end = start + STEP |
| train = ((idx >= start - WINDOW) & (idx <= start - GAP)) |
| test = (idx >= start) & (idx < end) |
| f = fit_factors(base, train) |
| season = f["var"]["+ Fed/CPI calendar"] |
| e2 = ret2 / season |
| p2 = park / season |
| X_har = np.column_stack([e2.rolling(w, min_periods=int(w * .9)).mean().to_numpy() for w in (24, 168, 720)]) |
| X_rng = np.column_stack([p2.rolling(w, min_periods=int(w * .9)).mean().to_numpy() for w in (24, 168, 720)]) |
| for h in HORIZONS: |
| fs = fwd_sum(season, h) |
| target = fwd_sum(e2, h) / h |
| scale_b = base["sigma"].to_numpy() * np.sqrt(fs) |
| preds = { |
| "HAR blend": log_ols(target, X_har, train)(X_har), |
| "Range estimator": X_rng[:, 1], |
| "HAR + range": log_ols(target, np.hstack([X_har, X_rng]), train)(np.hstack([X_har, X_rng])), |
| } |
| scales = {"B (live)": scale_b, **{k: np.sqrt(v * fs) for k, v in preds.items()}} |
| mv = moves[h] |
| for rule in RULES: |
| sc = scales["B (live)" if rule in ("Asymmetric", "Adaptive (ACI)") else rule] |
| z = mv / sc |
| ok_tr = train & np.isfinite(z) |
| ok_te = test & np.isfinite(z) & np.isfinite(scale_b) |
| if rule == "Asymmetric": |
| lo_q = {l: float(np.quantile(z[ok_tr], (1 - l) / 2)) for l in LEVELS} |
| hi_q = {l: float(np.quantile(z[ok_tr], 1 - (1 - l) / 2)) for l in LEVELS} |
| else: |
| q = {l: float(np.quantile(np.abs(z[ok_tr]), l)) for l in LEVELS} |
| lo_q, hi_q = {l: -v for l, v in q.items()}, q |
| rows[h][rule].append(pd.DataFrame({ |
| "z": z[ok_te], "scale": sc[ok_te], "move": mv[ok_te], |
| **{f"lo{l}": lo_q[l] for l in LEVELS}, **{f"hi{l}": hi_q[l] for l in LEVELS}, |
| "release": base["covered"][h].to_numpy()[ok_te]}, index=idx[ok_te])) |
| start = end |
|
|
| out, table = {}, {} |
| for h in HORIZONS: |
| frames = {r: pd.concat(rows[h][r]) for r in RULES} |
| common = frames["B (live)"].index |
| for r in RULES: |
| common = common.intersection(frames[r].index) |
| for r in RULES: |
| fr = frames[r].loc[common] |
| lo_q = {l: fr[f"lo{l}"].to_numpy() for l in LEVELS} |
| hi_q = {l: fr[f"hi{l}"].to_numpy() for l in LEVELS} |
| if r == "Adaptive (ACI)": |
| |
| |
| |
| zz = fr["z"].to_numpy() |
| for l in LEVELS: |
| theta, th = 0.0, np.zeros(len(zz)) |
| for i in range(len(zz)): |
| th[i] = theta |
| j = i - h |
| if j >= 0: |
| miss = not (lo_q[l][j] * math.exp(th[j]) <= zz[j] <= hi_q[l][j] * math.exp(th[j])) |
| theta += ACI_GAMMA * ((1.0 if miss else 0.0) - (1 - l)) |
| lo_q[l] = lo_q[l] * np.exp(th) |
| hi_q[l] = hi_q[l] * np.exp(th) |
| s = scores(None, None, fr["scale"].to_numpy(), fr["move"].to_numpy(), |
| {l: lo_q[l] for l in LEVELS}, {l: hi_q[l] for l in LEVELS}) |
| rel = fr["release"].to_numpy() |
| table[(h, r)] = {"score": s["score"], |
| "coverage": {f"{int(l * 100)}%": float(s[f"in{int(l * 100)}"].mean()) for l in LEVELS}, |
| "coverage_80_release": float(s["in80"][rel].mean()), |
| "mean_width_80_pct": float(s["width"].mean() * 100), |
| "interval_score_80_pct": float(s["score"].mean() * 100), "n": int(len(fr))} |
| b = table[(h, "B (live)")]["score"] |
| for r in RULES: |
| t = table[(h, r)] |
| if r != "B (live)": |
| t["dm_vs_B"] = dm(t["score"], b, h) |
| out[f"{h}h | {r}"] = {k: v for k, v in t.items() if k != "score"} |
|
|
| (ROOT / "outputs" / "price_range_upgrades.json").write_text(json.dumps({ |
| "protocol": "walk-forward as price_backtest_fixes.json; identical test hours for every rule", |
| "results": out}, indent=2)) |
| print(f"{'':4}{'rule':17s} {'50%':>6} {'80%':>6} {'95%':>6} {'release':>8} {'width':>7} " |
| f"{'score':>7} {'vs B':>7} {'DM t':>6} {'p':>6}") |
| for h in HORIZONS: |
| bs = out[f"{h}h | B (live)"]["interval_score_80_pct"] |
| for r in RULES: |
| t = out[f"{h}h | {r}"] |
| c = t["coverage"] |
| d = t.get("dm_vs_B") |
| print(f"{h:>2}h {r:17s} {c['50%']:>6.1%} {c['80%']:>6.1%} {c['95%']:>6.1%} " |
| f"{t['coverage_80_release']:>8.1%} {t['mean_width_80_pct']:>6.2f}% " |
| f"{t['interval_score_80_pct']:>7.3f} {t['interval_score_80_pct'] / bs - 1:>+7.1%}" |
| + (f" {d['t']:>+6.2f} {d['p']:>6.3f}" if d else "")) |
|
|
|
|
| if __name__ == "__main__": |
| main() |
|
|
