"""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 # noqa: E402 from backtest_price_range import HORIZONS, LEVELS # noqa: E402 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)) # Parkinson per-hour variance 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))) # Duan's retransformation 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"] # per-hour factor e2 = ret2 / season # deseasonalised 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) # sum of next h factors target = fwd_sum(e2, h) / h # realised deseasonalised var 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)": # Multiply each level's width by exp(theta); after every hour, # update theta with the outcome of the forecast made h hours # earlier, the newest one whose result is already known. 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()