| """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, |
| LOOKBACK, SURFACE, load, z_for) |
| from macro_event_study import EXTRA_RELEASES, RELEASES, release_times |
|
|
| 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"] |
| |
| CANDIDATE = "+ jobs report" |
| COLORS = ["#8b8a85", "#eb6834", "#1baf7a", "#2a78d6"] |
|
|
|
|
| GAP = pd.Timedelta(hours=max(HORIZONS)) |
| WINDOW = pd.DateOffset(years=2) |
| STEP = pd.DateOffset(months=3) |
|
|
|
|
| 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() |
| 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) |
| |
| |
| 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]: |
| |
| 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: |
| |
| 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 |
| 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) |
|
|
| |
| 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: |
| |
| |
| 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) |
|
|
| |
| 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") |
|
|
| |
| 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() |
|
|
