Spaces:
Sleeping
Sleeping
| """Backtest the live price forecast against what the price actually did. | |
| The forecast being tested is the one used live: | |
| centre = the current price (a random walk) | |
| range = current price x exp(+/- z * sigma * sqrt(h)), where sigma is the | |
| standard deviation of hourly log returns over the previous 168 | |
| hours | |
| It is issued at every hour of the test period and scored against the real | |
| price h hours later. Nothing about the test period is used to set the range: | |
| sigma at hour t comes only from hours before t. | |
| A probabilistic forecast is right when its stated confidence matches the | |
| outcome: an 80% range should contain the real price 80% of the time, a 50% | |
| range half the time. Too often means the range is wider than it needs to be; | |
| too rarely means it understates the risk. | |
| Usage: | |
| python scripts/backtest_price_range.py | |
| """ | |
| from __future__ import annotations | |
| import json | |
| import math | |
| from pathlib import Path | |
| import matplotlib | |
| matplotlib.use("Agg") | |
| import matplotlib.dates as mdates | |
| import matplotlib.pyplot as plt | |
| import numpy as np | |
| import pandas as pd | |
| ROOT = Path(__file__).resolve().parents[1] | |
| HORIZONS = (1, 6, 24, 72) | |
| LEVELS = (0.5, 0.8, 0.95) | |
| LOOKBACK = 168 | |
| START = "2023-01-01" | |
| SURFACE, INK, INK_2, INK_3, GRID = "#fcfcfb", "#0b0b0b", "#52514e", "#8b8a85", "#e6e5e1" | |
| SERIES = ["#2a78d6", "#eb6834", "#1baf7a", "#4a3aa7"] # validated order, slots 1,2,3,7 | |
| def z_for(level: float) -> float: | |
| """Two-sided normal quantile without scipy (blocked on some machines).""" | |
| p = 0.5 + level / 2 | |
| lo, hi = 0.0, 10.0 | |
| for _ in range(80): | |
| mid = (lo + hi) / 2 | |
| if 0.5 * math.erfc(-mid / math.sqrt(2)) < p: | |
| lo = mid | |
| else: | |
| hi = mid | |
| return (lo + hi) / 2 | |
| def load() -> pd.Series: | |
| btc = pd.read_csv(ROOT / "outputs" / "volatility" / "btc_1h.csv", | |
| usecols=["timestamp", "close"]) | |
| btc["timestamp"] = pd.to_datetime(btc["timestamp"], utc=True) | |
| close = btc.set_index("timestamp")["close"].sort_index() | |
| return close[~close.index.duplicated()].asfreq("1h") | |
| def backtest(close: pd.Series) -> pd.DataFrame: | |
| logp = np.log(close) | |
| ret = logp.diff() | |
| sigma = ret.rolling(LOOKBACK, min_periods=int(LOOKBACK * 0.9)).std() # known at t | |
| rows = [] | |
| for h in HORIZONS: | |
| future = logp.shift(-h) | |
| frame = pd.DataFrame({"logp": logp, "future": future, "sigma": sigma}).dropna() | |
| frame = frame[frame.index >= START] | |
| z_score = (frame["future"] - frame["logp"]) / (frame["sigma"] * math.sqrt(h)) | |
| price, actual = np.exp(frame["logp"]), np.exp(frame["future"]) | |
| for t, zs, p, a, s in zip(frame.index, z_score, price, actual, frame["sigma"]): | |
| rows.append((h, t, zs, p, a, s)) | |
| return pd.DataFrame(rows, columns=["h", "t", "z", "price", "actual", "sigma"]) | |
| def summarise(bt: pd.DataFrame) -> dict: | |
| out = {} | |
| for h, g in bt.groupby("h"): | |
| absz = g["z"].abs() | |
| err = (g["actual"] - g["price"]).abs() | |
| z80 = z_for(0.8) | |
| width = (np.exp(z80 * g["sigma"] * math.sqrt(h)) - np.exp(-z80 * g["sigma"] * math.sqrt(h))) * g["price"] | |
| by_year = {int(y): float((gy["z"].abs() <= z80).mean()) | |
| for y, gy in g.groupby(g["t"].dt.year)} | |
| weekend = g["t"].dt.dayofweek >= 5 | |
| out[f"{h}h"] = { | |
| "forecasts": int(len(g)), | |
| "coverage": {f"{int(l * 100)}%": float((absz <= z_for(l)).mean()) for l in LEVELS}, | |
| "coverage_80_by_year": by_year, | |
| "coverage_80_weekday": float((absz[~weekend] <= z80).mean()), | |
| "coverage_80_weekend": float((absz[weekend] <= z80).mean()), | |
| "mean_abs_error_usd": float(err.mean()), | |
| "median_abs_error_pct": float((err / g["price"]).median() * 100), | |
| "mean_range_width_usd": float(width.mean()), | |
| "share_up": float((g["actual"] > g["price"]).mean()), | |
| } | |
| return out | |
| def draw(close: pd.Series, bt: pd.DataFrame, summary: dict) -> None: | |
| plt.rcParams.update({"font.family": "DejaVu Sans", "text.color": INK, | |
| "axes.labelcolor": INK_2}) | |
| fig = plt.figure(figsize=(13, 9.5), facecolor=SURFACE) | |
| grid = fig.add_gridspec(2, 2, height_ratios=[1.2, 1], hspace=0.42, wspace=0.25, | |
| left=0.07, right=0.98, top=0.86, bottom=0.08) | |
| def style(ax): | |
| ax.set_facecolor(SURFACE) | |
| for side in ("top", "right"): | |
| ax.spines[side].set_visible(False) | |
| for side in ("left", "bottom"): | |
| ax.spines[side].set_color(GRID) | |
| ax.tick_params(colors=INK_2, labelsize=8.5, length=0) | |
| ax.grid(color=GRID, linewidth=0.8) | |
| ax.set_axisbelow(True) | |
| # --- the last 60 days: forecast range issued 24h earlier vs the real price ---- | |
| ax = fig.add_subplot(grid[0, :]); style(ax) | |
| g = bt[bt["h"] == 24].set_index("t") | |
| end = g.index.max() | |
| g = g[g.index >= end - pd.Timedelta(days=60)] | |
| z80 = z_for(0.8) | |
| shift = pd.Timedelta(hours=24) | |
| when = g.index + shift # the hour each forecast is about | |
| lo = g["price"] * np.exp(-z80 * g["sigma"] * math.sqrt(24)) | |
| hi = g["price"] * np.exp(z80 * g["sigma"] * math.sqrt(24)) | |
| ax.fill_between(when, lo, hi, color=SERIES[0], alpha=0.16, linewidth=0, | |
| label="80% range, forecast 24h earlier") | |
| ax.plot(when, g["actual"], color=INK, linewidth=1.3, label="actual price") | |
| miss = (g["actual"] < lo) | (g["actual"] > hi) | |
| ax.scatter(when[miss], g["actual"][miss], color=SERIES[1], s=10, zorder=3, | |
| label=f"outside the range ({miss.mean():.0%} of hours)") | |
| ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda v, _: f"${v / 1000:,.0f}k")) | |
| ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) | |
| ax.legend(loc="upper left", frameon=False, fontsize=8.5) | |
| ax.set_title("The last 60 days of the test: 24-hour forecast range vs the real price", | |
| loc="left", fontsize=11, fontweight="bold") | |
| # --- calibration: stated confidence vs how often it held ------------------------ | |
| ax = fig.add_subplot(grid[1, 0]); style(ax) | |
| ax.plot([40, 100], [40, 100], color=INK_3, linewidth=1, linestyle=(0, (4, 2))) | |
| for i, h in enumerate(HORIZONS): | |
| cov = summary[f"{h}h"]["coverage"] | |
| xs = [int(k[:-1]) for k in cov] | |
| ys = [v * 100 for v in cov.values()] | |
| ax.plot(xs, ys, color=SERIES[i], linewidth=2, marker="o", markersize=5, | |
| markeredgecolor=SURFACE, label=f"{h}h ahead") | |
| 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.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, 49, "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("Are the ranges honest?", loc="left", fontsize=11, fontweight="bold") | |
| # --- 80% coverage by year ----------------------------------------------------------- | |
| ax = fig.add_subplot(grid[1, 1]); style(ax) | |
| years = sorted(summary["24h"]["coverage_80_by_year"]) | |
| width = 0.2 | |
| for i, h in enumerate(HORIZONS): | |
| vals = [summary[f"{h}h"]["coverage_80_by_year"].get(y, np.nan) * 100 for y in years] | |
| ax.bar(np.arange(len(years)) + (i - 1.5) * width, vals, width * 0.9, color=SERIES[i], | |
| label=f"{h}h") | |
| ax.axhline(80, color=INK_2, linewidth=1, linestyle=(0, (4, 2))) | |
| ax.set_xticks(range(len(years))) | |
| ax.set_xticklabels([str(y) + (" (to Sep)" if y == years[-1] else "") for y in years]) | |
| ax.set_ylim(60, 100) | |
| ax.set_ylabel("80% range: share of real prices inside", fontsize=8.5) | |
| ax.legend(loc="upper center", frameon=False, fontsize=8, ncol=4) | |
| ax.set_title("Does it hold every year?", loc="left", fontsize=11, fontweight="bold") | |
| n = summary["24h"]["forecasts"] | |
| fig.suptitle("Backtest: the live price forecast against the real price", x=0.07, ha="left", | |
| fontsize=15, fontweight="bold", y=0.975) | |
| fig.text(0.07, 0.905, | |
| f"Issued every hour, {START[:4]} – Sep 2026 ({n:,} forecasts per horizon). Centre = " | |
| "current price; range from the previous 7 days' hourly volatility only, so no " | |
| "forecast sees its own future.", | |
| fontsize=9.3, color=INK_2, ha="left") | |
| fig.savefig(ROOT / "outputs" / "figures" / "fig11_price_backtest.png", dpi=150, | |
| facecolor=SURFACE) | |
| plt.close(fig) | |
| print("written: outputs/figures/fig11_price_backtest.png") | |
| def main() -> None: | |
| close = load() | |
| bt = backtest(close) | |
| summary = summarise(bt) | |
| draw(close, bt, summary) | |
| (ROOT / "outputs" / "price_backtest.json").write_text(json.dumps({ | |
| "method": f"centre = current price; range = price*exp(+/-z*sigma*sqrt(h)), sigma = std " | |
| f"of the previous {LOOKBACK} hourly log returns; scored every hour from {START}", | |
| "horizons": summary}, indent=2)) | |
| for h, s in summary.items(): | |
| cov = " ".join(f"{k}->{v:.1%}" for k, v in s["coverage"].items()) | |
| print(f"{h:>4}: {s['forecasts']:,} forecasts | {cov} | weekday {s['coverage_80_weekday']:.1%} " | |
| f"weekend {s['coverage_80_weekend']:.1%} | median miss {s['median_abs_error_pct']:.2f}% " | |
| f"| up {s['share_up']:.1%}") | |
| print(f" 80% by year: " + " ".join(f"{y}:{v:.0%}" for y, v in s["coverage_80_by_year"].items())) | |
| if __name__ == "__main__": | |
| main() | |