"""Three fixes to the price-range forecast, learned on 2019-2022, tested on 2023-2026. backtest_price_range.py found the plain forecast honest at 80% overall but wrong in three systematic ways. Each fix targets one, and is added in turn so its contribution can be read separately: 1. Fat tails. The range assumed a normal distribution, so 50% ranges held ~59% and 95% ranges ~92%. The fix replaces normal quantiles with the empirical quantiles of standardised outcomes in the training years. 2. Hour of the week. A week of volatility mixes busy weekday afternoons into quiet weekend nights. Each hour of the week gets a variance factor, and an h-hour range sums the factors of the hours it actually covers. 3. The release calendar. Fed decision and CPI hours, and the two after each, get a further variance factor. Their dates are published weeks ahead, so using them is not look-ahead. Everything is fitted on 2019-2022 and frozen. The 2023-2026 numbers are therefore what each rule would have delivered live. Scoring: coverage at 50/80/95%, split by weekday, weekend and release hours, and the 80% interval score (Gneiting & Raftery 2007): range width plus a penalty of 2/alpha times any miss. Lower is better; it rewards a range that is narrow where it can be and wide where it must be, which coverage alone does not. Usage: python scripts/backtest_improvements.py """ from __future__ import annotations import json import math import sys from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd sys.path.insert(0, str(Path(__file__).resolve().parent)) from backtest_price_range import (GRID, HORIZONS, INK, INK_2, INK_3, LEVELS, # noqa: E402 LOOKBACK, SURFACE, load, z_for) from macro_event_study import EXTRA_RELEASES, RELEASES, release_times # noqa: E402 ROOT = Path(__file__).resolve().parents[1] TRAIN_END, TEST_START = "2022-12-31 23:00", "2023-01-01" CONFIGS = ["plain", "+ fat tails", "+ hour of week", "+ Fed/CPI calendar"] # Scored and reported, not drawn: does the jobs report earn a place in the calendar? CANDIDATE = "+ jobs report" COLORS = ["#8b8a85", "#eb6834", "#1baf7a", "#2a78d6"] # grey baseline, then slots 2, 3, 1 GAP = pd.Timedelta(hours=max(HORIZONS)) # a training outcome must end before testing starts WINDOW = pd.DateOffset(years=2) # walk-forward: refit on the previous two years STEP = pd.DateOffset(months=3) # ... every quarter def prepare() -> dict: close = load() logp = np.log(close) ret = logp.diff() sigma = ret.rolling(LOOKBACK, min_periods=int(LOOKBACK * 0.9)).std() rel_hours = set() for dates, clock in RELEASES.values(): rel_hours |= {t.floor("h") for t in release_times(dates, clock)} release = ret.index.isin(list(rel_hours)) after = ret.index.isin([t + pd.Timedelta(hours=k) for t in rel_hours for k in (1, 2)]) dates, clock = EXTRA_RELEASES["Jobs report (NFP)"] nfp_hours = {t.floor("h") for t in release_times(dates, clock)} - rel_hours nfp = ret.index.isin(list(nfp_hours)) nfp_after = ret.index.isin([t + pd.Timedelta(hours=k) for t in nfp_hours for k in (1, 2)]) & ~release return {"logp": logp, "ret": ret, "sigma": sigma, "release": release, "after": after, "nfp": nfp, "nfp_after": nfp_after, "how": ret.index.dayofweek * 24 + ret.index.hour, "u2": (ret / sigma) ** 2, "weekend": ret.index.dayofweek >= 5, "covered": {h: pd.Series(release, index=ret.index).astype(float) .rolling(h).max().shift(-h).fillna(0).astype(bool) for h in HORIZONS}} def fit_factors(base: dict, train: np.ndarray) -> dict: """Hour-of-week and release variance factors from the training hours only.""" u2, how = base["u2"], base["how"] normal = train & ~base["release"] & ~base["after"] & u2.notna().to_numpy() f2 = u2[normal].groupby(how[normal]).mean() f2 = f2 / f2.mean() # a week averages to 1 week = pd.Series(how, index=u2.index).map(f2).astype(float).fillna(1.0) expected = u2 / week k_rel = float(expected[train & base["release"]].mean()) k_after = float(expected[train & base["after"]].mean()) cal = np.where(base["release"], k_rel, np.where(base["after"], k_after, 1.0)) k_nfp = float(expected[train & base["nfp"]].mean()) k_nfp_after = float(expected[train & base["nfp_after"] & ~base["after"]].mean()) nfp_cal = np.where(base["nfp"], k_nfp, np.where(base["nfp_after"] & ~base["after"], k_nfp_after, 1.0)) ones = pd.Series(1.0, index=u2.index) # Variance per future hour under each rule. Known at forecast time: the # week pattern comes from the past and the calendar is published ahead. return {"var": {"plain": ones, "+ fat tails": ones, "+ hour of week": week, "+ Fed/CPI calendar": week * cal, CANDIDATE: week * cal * nfp_cal}, "k_nfp": k_nfp, "k_release": k_rel, "k_after": k_after, "weekend_factor": float(f2[f2.index >= 120].mean())} def score_rows(base: dict, factors: dict, train: np.ndarray, test: np.ndarray) -> pd.DataFrame: """Per-forecast outcomes on the test hours, every horizon and rule.""" out = [] index = base["ret"].index for h in HORIZONS: move = base["logp"].shift(-h) - base["logp"] for name in CONFIGS + [CANDIDATE]: # sigma_h^2 = sigma^2 * the sum of the next h hours' variance factors. v = factors["var"][name].rolling(h).sum().shift(-h) scale = base["sigma"] * np.sqrt(v) z = move / scale ok = z.notna().to_numpy() if name == "plain": q = {l: z_for(l) for l in LEVELS} else: # Fat tails: the training hours' own quantiles of |z|. q = {l: float(z[train & ok].abs().quantile(l)) for l in LEVELS} sel = test & ok zt, st, mt = z[sel], scale[sel], move[sel] half = q[0.8] * st lo_r, hi_r = np.exp(-half) - 1, np.exp(half) - 1 # range as % of price y = np.exp(mt) - 1 alpha = 0.2 frame = pd.DataFrame({ "h": h, "rule": name, **{f"in{int(l * 100)}": (zt.abs() <= q[l]).to_numpy() for l in LEVELS}, "weekend": base["weekend"][sel], "release": base["covered"][h].to_numpy()[sel], "width": (hi_r - lo_r).to_numpy(), "score": ((hi_r - lo_r) + (2 / alpha) * ((lo_r - y).clip(lower=0) + (y - hi_r).clip(lower=0))).to_numpy(), }, index=index[sel]) out.append(frame) return pd.concat(out) def summarise(rows: pd.DataFrame) -> dict: res = {} for (h, name), g in rows.groupby(["h", "rule"]): res[(h, name)] = { "coverage": {f"{int(l * 100)}%": float(g[f"in{int(l * 100)}"].mean()) for l in LEVELS}, "coverage_80_weekday": float(g.loc[~g["weekend"], "in80"].mean()), "coverage_80_weekend": float(g.loc[g["weekend"], "in80"].mean()), "coverage_80_release": float(g.loc[g["release"], "in80"].mean()), "n_release": int(g["release"].sum()), "interval_score_80_pct": float(g["score"].mean() * 100), "mean_width_80_pct": float(g["width"].mean() * 100), "n": int(len(g)), } return res def single_split(base: dict) -> tuple[dict, dict]: index = base["ret"].index train = (index <= pd.Timestamp(TRAIN_END, tz="UTC") - GAP) test = (index >= pd.Timestamp(TEST_START, tz="UTC")) factors = fit_factors(base, train) return summarise(score_rows(base, factors, train, test)), factors def walk_forward(base: dict) -> tuple[dict, list]: """Refit every quarter on the two years ending GAP before it, test the quarter.""" index = base["ret"].index start = pd.Timestamp(TEST_START, tz="UTC") rows, fits = [], [] while start < index.max(): end = start + STEP train = (index >= start - WINDOW) & (index <= start - GAP) test = (index >= start) & (index < end) factors = fit_factors(base, train) rows.append(score_rows(base, factors, train, test)) fits.append({"quarter": f"{start:%Y-%m}", "k_release": factors["k_release"], "k_nfp": factors["k_nfp"], "weekend_factor": factors["weekend_factor"]}) start = end return summarise(pd.concat(rows)), fits def run() -> dict: base = prepare() frozen, factors = single_split(base) walked, fits = walk_forward(base) return {"results": walked, "frozen": frozen, "fits": fits, "k_release": float(np.mean([f["k_release"] for f in fits])), "k_after": factors["k_after"], "weekend_factor": float(np.mean([f["weekend_factor"] for f in fits])), "weekend_factor_range": (float(min(f["weekend_factor"] for f in fits)), float(max(f["weekend_factor"] for f in fits)))} def draw(out: dict) -> None: res = out["results"] plt.rcParams.update({"font.family": "DejaVu Sans", "text.color": INK, "axes.labelcolor": INK_2}) fig = plt.figure(figsize=(13.5, 9.5), facecolor=SURFACE) grid = fig.add_gridspec(2, 2, hspace=0.5, wspace=0.28, left=0.07, right=0.98, top=0.81, bottom=0.08) def style(ax): ax.set_facecolor(SURFACE) for s in ("top", "right"): ax.spines[s].set_visible(False) for s in ("left", "bottom"): ax.spines[s].set_color(GRID) ax.tick_params(colors=INK_2, labelsize=8.5, length=0) ax.grid(color=GRID, linewidth=0.8, axis="y") ax.set_axisbelow(True) # A: 80% coverage where the plain range failed --------------------------------- for col, h in enumerate((1, 24)): ax = fig.add_subplot(grid[0, col]); style(ax) groups = [("all hours", lambda r: r["coverage"]["80%"]), ("weekdays", lambda r: r["coverage_80_weekday"]), ("weekends", lambda r: r["coverage_80_weekend"]), ("Fed/CPI hours" if h == 1 else "windows with Fed/CPI", lambda r: r["coverage_80_release"])] width = 0.2 for i, name in enumerate(CONFIGS): vals = [fn(res[(h, name)]) * 100 for _, fn in groups] ax.bar(np.arange(len(groups)) + (i - 1.5) * width, vals, width * 0.92, color=COLORS[i], label=name, zorder=2) ax.axhline(80, color=INK, linewidth=1.1, linestyle=(0, (4, 2)), zorder=3) ax.set_xticks(range(len(groups))) ax.set_xticklabels([g for g, _ in groups], fontsize=9) ax.set_ylim(40, 100) ax.set_ylabel("real price inside the 80% range, %", fontsize=8.5) ax.set_title(f"{h}h ahead: does the 80% range hold everywhere?", loc="left", fontsize=11, fontweight="bold") if col == 0: # One legend for the whole figure, above the panels, where it covers # no bars; all four bar panels share the same colours. handles, labels = ax.get_legend_handles_labels() fig.legend(handles, labels, loc="upper left", bbox_to_anchor=(0.07, 0.875), ncol=4, frameon=False, fontsize=9) # B: calibration at 24h ----------------------------------------------------------- ax = fig.add_subplot(grid[1, 0]); style(ax) ax.plot([45, 100], [45, 100], color=INK_3, linewidth=1, linestyle=(0, (4, 2))) for i, name in enumerate((CONFIGS[0], CONFIGS[-1])): cov = res[(24, name)]["coverage"] ax.plot([int(k[:-1]) for k in cov], [v * 100 for v in cov.values()], color=COLORS[0 if i == 0 else -1], linewidth=2.2, marker="o", markersize=6, markeredgecolor=SURFACE, label=f"{name}" if i == 0 else "all three fixes") ax.set_xlim(45, 100); ax.set_ylim(45, 100) ax.set_xticks([50, 80, 95]); ax.set_yticks([50, 60, 70, 80, 90, 95, 100]) ax.grid(axis="x", color=GRID, linewidth=0.8) ax.set_xlabel("stated confidence of the range, %", fontsize=8.5) ax.set_ylabel("how often the real price landed inside, %", fontsize=8.5) ax.text(97, 48, "on the dashed line = honest", ha="right", fontsize=8, color=INK_2) ax.legend(loc="upper left", frameon=False, fontsize=8.5) ax.set_title("24h ahead: is each confidence level honest?", loc="left", fontsize=11, fontweight="bold") # C: interval score ------------------------------------------------------------------- ax = fig.add_subplot(grid[1, 1]); style(ax) width = 0.2 for i, name in enumerate(CONFIGS): vals = [(res[(h, name)]["interval_score_80_pct"] / res[(h, "plain")]["interval_score_80_pct"] - 1) * 100 for h in HORIZONS] bars = ax.bar(np.arange(len(HORIZONS)) + (i - 1.5) * width, vals, width * 0.92, color=COLORS[i], zorder=2) if name == CONFIGS[-1]: for b, v in zip(bars, vals): ax.text(b.get_x() + b.get_width() / 2, v - 0.25, f"{v:+.1f}%", ha="center", va="top", fontsize=8, color=INK_2) ax.axhline(0, color=INK_2, linewidth=0.9) ax.set_xticks(range(len(HORIZONS))) ax.set_xticklabels([f"{h}h ahead" for h in HORIZONS]) ax.set_ylabel("80% interval score vs plain, % (lower = better)", fontsize=8.5) ax.set_title("Overall forecast quality", loc="left", fontsize=11, fontweight="bold") n = res[(24, "plain")]["n"] fig.suptitle("Three fixes to the price-range forecast, tested on data they never saw", x=0.07, ha="left", fontsize=15, fontweight="bold", y=0.975) lo_w, hi_w = out["weekend_factor_range"] fig.text(0.07, 0.9, f"Walk-forward: every quarter each fix is refitted on the previous two years " f"(ending 72h before the quarter) and scored on the quarter, 2023 – Sep 2026 " f"({n:,} forecasts per horizon).\nFixes added in turn, left to right. Weekend " f"hours carried {lo_w:.2f}–{hi_w:.2f}× average variance across refits; a Fed/CPI " f"release hour {out['k_release']:.1f}× its usual variance on average.", fontsize=9.3, color=INK_2, ha="left") fig.savefig(ROOT / "outputs" / "figures" / "fig12_forecast_fixes.png", dpi=150, facecolor=SURFACE) plt.close(fig) print("written: outputs/figures/fig12_forecast_fixes.png") def main() -> None: out = run() draw(out) res = out["results"] (ROOT / "outputs" / "price_backtest_fixes.json").write_text(json.dumps({ "protocol": ("walk-forward: refit quarterly on the previous 2 years ending 72h before " "each quarter; 'frozen' = single fit on 2019-2022 for comparison"), "test": f"{TEST_START} to Sep 2026", "refits": out["fits"], "walk_forward": {f"{h}h | {name}": v for (h, name), v in res.items()}, "frozen_split": {f"{h}h | {name}": v for (h, name), v in out["frozen"].items()}}, indent=2)) for label, table in (("WALK-FORWARD", res), ("FROZEN SPLIT (fit 2019-2022)", out["frozen"])): print(f"\n=== {label} ===") for h in HORIZONS: print(f"{h}h {'rule':20s} {'50%':>6} {'80%':>6} {'95%':>6} {'wkday':>6} " f"{'wkend':>6} {'release':>8} {'width':>7} {'score':>7}") for name in CONFIGS + [CANDIDATE]: r = table[(h, name)] c = r["coverage"] print(f" {name:20s} {c['50%']:>6.1%} {c['80%']:>6.1%} {c['95%']:>6.1%} " f"{r['coverage_80_weekday']:>6.1%} {r['coverage_80_weekend']:>6.1%} " f"{r['coverage_80_release']:>8.1%} {r['mean_width_80_pct']:>6.2f}% " f"{r['interval_score_80_pct']:>6.3f}") print("\nweekend factor per refit:", " ".join(f"{f['quarter']}:{f['weekend_factor']:.2f}" for f in out["fits"])) if __name__ == "__main__": main()