scripts/backtest_range_upgrades.py

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