Rewrite as a probabilistic model with walk-forward evaluation.
Replace the 2024 model (model.py, ~2000 lines) with the btcmodel package, the baseline for future work: - Forecasts are quantiles of log price at each horizon, scored with CRPS in a walk-forward backtest (origins every 30 days from 2014, horizons 1 month to 4 years). Skill is relative to a zero-drift random walk, with circular block-bootstrap intervals and a count of independent windows. - Development data stops at 2024-11-26, the last day the 2024 model saw. Later outcomes are a holdout, scored only by `backtest --holdout`. - Models: random_walk, drift_rw, and cycle (the 2024 model's cycle-position drift, now kernel-smoothed and recency-weighted). On development data nothing beats the random walk with confidence; cycle loses at every horizon. - Prices: the Investing.com archive moves to data/ (cut at 2024-11-26; its last row was intraday) and is extended with Coinbase daily closes by `update`. Also: Nix flake dev shell (Python 3.13, pandas 3), ruff in place of black, pytest suite, and a rewritten README. NOTES.md is removed as inaccurate, and poetry is dropped.
This commit is contained in:
@@ -0,0 +1,100 @@
|
||||
"""
|
||||
Command line entry point.
|
||||
|
||||
python -m btcmodel update fetch new daily prices from Coinbase
|
||||
python -m btcmodel backtest score models on development data
|
||||
python -m btcmodel backtest --holdout score models on outcomes after DEV_CUTOFF
|
||||
python -m btcmodel forecast forecast from the latest price
|
||||
"""
|
||||
|
||||
import argparse
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from . import data, evaluate, plots
|
||||
from .models import MODELS
|
||||
|
||||
FORECAST_REPORT_HORIZONS = (182, 365, 730, 1095, 1460)
|
||||
FORECAST_REPORT_LEVELS = (0.05, 0.25, 0.5, 0.75, 0.95)
|
||||
|
||||
|
||||
def main() -> None:
|
||||
parser = argparse.ArgumentParser(prog="btcmodel", description="Bitcoin price model")
|
||||
parser.add_argument("-o", "--output", type=Path, default=Path("output"))
|
||||
parser.add_argument("-m", "--models", nargs="+", choices=list(MODELS), default=list(MODELS))
|
||||
commands = parser.add_subparsers(dest="command", required=True)
|
||||
commands.add_parser("update", help="fetch new daily prices")
|
||||
backtest = commands.add_parser("backtest", help="walk-forward evaluation")
|
||||
backtest.add_argument(
|
||||
"--holdout",
|
||||
action="store_true",
|
||||
help=f"score outcomes after {data.DEV_CUTOFF:%Y-%m-%d} (don't use while developing)",
|
||||
)
|
||||
commands.add_parser("forecast", help="forecast from the latest price")
|
||||
args = parser.parse_args()
|
||||
models = [MODELS[name] for name in args.models]
|
||||
|
||||
if args.command == "update":
|
||||
added = data.update_coinbase()
|
||||
print(f"added {added} days; latest {data.load_prices().index[-1]:%Y-%m-%d}")
|
||||
elif args.command == "backtest":
|
||||
run_backtest(models, args.output, args.holdout)
|
||||
elif args.command == "forecast":
|
||||
run_forecast(models, args.output)
|
||||
|
||||
|
||||
def run_backtest(models, output: Path, holdout: bool) -> None:
|
||||
if holdout:
|
||||
name, prices = "holdout", data.load_prices()
|
||||
scores = evaluate.backtest(models, prices, score_after=data.DEV_CUTOFF)
|
||||
else:
|
||||
name, prices = "backtest", data.load_prices(until=data.DEV_CUTOFF)
|
||||
scores = evaluate.backtest(models, prices)
|
||||
out = output / name
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
summary = evaluate.summarize(scores)
|
||||
report = (
|
||||
f"{name}: origins every {evaluate.ORIGIN_STEP_DAYS} days from "
|
||||
f"{evaluate.FIRST_ORIGIN:%Y-%m-%d}, outcomes through {prices.index[-1]:%Y-%m-%d}\n\n"
|
||||
+ evaluate.format_summary(summary)
|
||||
)
|
||||
print(report)
|
||||
(out / "report.txt").write_text(report + "\n")
|
||||
scores.to_csv(out / "scores.csv", index=False)
|
||||
summary.to_csv(out / "summary.csv", index=False)
|
||||
plots.skill_chart(summary, out / "skill.png", f"{name}: skill by horizon")
|
||||
plots.calibration_chart(summary, out / "calibration.png", f"{name}: interval coverage")
|
||||
print(f"\nwrote {out}/")
|
||||
|
||||
|
||||
def run_forecast(models, output: Path) -> None:
|
||||
prices = data.load_prices()
|
||||
horizons = np.arange(1, max(FORECAST_REPORT_HORIZONS) + 1)
|
||||
forecasts = {m.name: m.forecast(prices, horizons) for m in models}
|
||||
out = output / "forecast"
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
rows = []
|
||||
for name, f in forecasts.items():
|
||||
for h in FORECAST_REPORT_HORIZONS:
|
||||
row = {"model": name, "date": f.dates[h - 1].date(), "horizon": h}
|
||||
for level in FORECAST_REPORT_LEVELS:
|
||||
row[f"p{level * 100:02.0f}"] = np.exp(f.quantile(level)[h - 1])
|
||||
rows.append(row)
|
||||
table = pd.DataFrame(rows)
|
||||
table.to_csv(out / "forecast.csv", index=False)
|
||||
|
||||
shown = table.copy()
|
||||
for column in shown.columns[3:]:
|
||||
shown[column] = shown[column].map(plots.price_formatter)
|
||||
print(f"from {prices.index[-1]:%Y-%m-%d} at {plots.price_formatter(prices.close.iloc[-1])}\n")
|
||||
print(shown.to_string(index=False))
|
||||
plots.fan_chart(prices, forecasts, out / "fan.png")
|
||||
print(f"\nwrote {out}/")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,117 @@
|
||||
"""Daily BTC-USD closing prices.
|
||||
|
||||
Two sources, stitched at ARCHIVE_END:
|
||||
|
||||
- data/investing.csv: the original Investing.com download (2010-07-18 onward).
|
||||
Its last row (2024-11-27) was an intraday snapshot, so it is cut a day early.
|
||||
- data/coinbase.csv: Coinbase Exchange daily candles (UTC days), appended by
|
||||
`python -m btcmodel update`.
|
||||
|
||||
Everything up to ARCHIVE_END is development data. Everything after it is the
|
||||
holdout: outcomes nobody had seen while the 2024 model was being built, and which
|
||||
model development here must not look at (see evaluate.py).
|
||||
"""
|
||||
|
||||
import datetime as dt
|
||||
import json
|
||||
import urllib.request
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
DATA_DIR = Path(__file__).resolve().parent.parent / "data"
|
||||
ARCHIVE_CSV = DATA_DIR / "investing.csv"
|
||||
COINBASE_CSV = DATA_DIR / "coinbase.csv"
|
||||
|
||||
ARCHIVE_END = pd.Timestamp("2024-11-26")
|
||||
DEV_CUTOFF = ARCHIVE_END
|
||||
|
||||
# 2010 has four distinct prices and no change on 89% of days; it is noise.
|
||||
DATA_START = pd.Timestamp("2011-01-01")
|
||||
|
||||
COINBASE_URL = "https://api.exchange.coinbase.com/products/BTC-USD/candles"
|
||||
COINBASE_MAX_CANDLES = 300
|
||||
|
||||
|
||||
def load_prices(
|
||||
until: pd.Timestamp | str | None = None, start: pd.Timestamp | str = DATA_START
|
||||
) -> pd.DataFrame:
|
||||
"""
|
||||
Load daily closes as a frame indexed by date, with a single `close` column.
|
||||
|
||||
Models receive a prefix of this frame, so additional data sources can be
|
||||
joined in as extra columns later without changing the model interface.
|
||||
"""
|
||||
archive = _read_archive()
|
||||
parts = [archive[archive.index <= ARCHIVE_END]]
|
||||
if COINBASE_CSV.exists():
|
||||
coinbase = pd.read_csv(COINBASE_CSV, index_col="date", parse_dates=["date"])
|
||||
parts.append(coinbase.loc[coinbase.index > ARCHIVE_END, ["close"]])
|
||||
df = pd.concat(parts).sort_index()
|
||||
|
||||
df = df[df.index >= pd.Timestamp(start)]
|
||||
if until is not None:
|
||||
df = df[df.index <= pd.Timestamp(until)]
|
||||
|
||||
expected = pd.date_range(df.index[0], df.index[-1], freq="D")
|
||||
missing = expected.difference(df.index)
|
||||
if len(missing):
|
||||
raise ValueError(f"{len(missing)} missing days, first {missing[0].date()}")
|
||||
if not (df["close"] > 0).all():
|
||||
raise ValueError("non-positive closing price")
|
||||
return df
|
||||
|
||||
|
||||
def log_returns(df: pd.DataFrame) -> pd.Series:
|
||||
"""Daily log returns of the close, without the leading NaN."""
|
||||
return np.log(df["close"]).diff().iloc[1:]
|
||||
|
||||
|
||||
def _read_archive() -> pd.DataFrame:
|
||||
raw = pd.read_csv(ARCHIVE_CSV, encoding="utf-8-sig", thousands=",")
|
||||
dates = pd.to_datetime(raw["Date"], format="%m/%d/%Y")
|
||||
return pd.DataFrame({"close": raw["Price"].astype(float).values}, index=dates.rename("date"))
|
||||
|
||||
|
||||
def update_coinbase() -> int:
|
||||
"""Append completed daily candles since the last stored day. Returns rows added."""
|
||||
if COINBASE_CSV.exists():
|
||||
existing = pd.read_csv(COINBASE_CSV, index_col="date", parse_dates=["date"])
|
||||
first = existing.index.max() + pd.Timedelta(days=1)
|
||||
else:
|
||||
existing = None
|
||||
first = ARCHIVE_END + pd.Timedelta(days=1)
|
||||
# Today's candle is still forming; only take finished UTC days.
|
||||
last = pd.Timestamp(dt.datetime.now(dt.UTC).date()) - pd.Timedelta(days=1)
|
||||
|
||||
rows = []
|
||||
chunk_start = first
|
||||
while chunk_start <= last:
|
||||
chunk_end = min(chunk_start + pd.Timedelta(days=COINBASE_MAX_CANDLES - 1), last)
|
||||
rows.extend(_fetch_candles(chunk_start, chunk_end))
|
||||
chunk_start = chunk_end + pd.Timedelta(days=1)
|
||||
if not rows:
|
||||
return 0
|
||||
|
||||
new = pd.DataFrame(rows, columns=["date", "close"]).set_index("date").sort_index()
|
||||
new = new[(new.index >= first) & (new.index <= last)]
|
||||
combined = new if existing is None else pd.concat([existing, new])
|
||||
combined = combined[~combined.index.duplicated(keep="last")].sort_index()
|
||||
combined.to_csv(COINBASE_CSV, date_format="%Y-%m-%d")
|
||||
return len(new)
|
||||
|
||||
|
||||
def _fetch_candles(start: pd.Timestamp, end: pd.Timestamp) -> list[tuple[pd.Timestamp, float]]:
|
||||
url = (
|
||||
f"{COINBASE_URL}?granularity=86400"
|
||||
f"&start={start:%Y-%m-%d}T00:00:00Z&end={end:%Y-%m-%d}T00:00:00Z"
|
||||
)
|
||||
# Coinbase rejects requests without a User-Agent.
|
||||
request = urllib.request.Request(url, headers={"User-Agent": "btcmodel"})
|
||||
with urllib.request.urlopen(request, timeout=30) as response:
|
||||
candles = json.load(response)
|
||||
# Each candle is [time, low, high, open, close, volume].
|
||||
return [
|
||||
(pd.Timestamp(dt.datetime.fromtimestamp(c[0], dt.UTC).date()), float(c[4])) for c in candles
|
||||
]
|
||||
@@ -0,0 +1,123 @@
|
||||
"""
|
||||
Walk-forward evaluation.
|
||||
|
||||
From each origin (every ORIGIN_STEP_DAYS from FIRST_ORIGIN), each model sees the
|
||||
data up to that day only and forecasts every horizon. Each forecast whose target
|
||||
date has been observed is scored against what happened.
|
||||
|
||||
Development runs load data only up to DEV_CUTOFF, so outcomes after it cannot
|
||||
influence model design. The holdout run scores only targets after DEV_CUTOFF.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from .forecast import crps, pit
|
||||
from .models import BASELINE
|
||||
|
||||
HORIZONS = np.array([30, 91, 182, 365, 730, 1095, 1460])
|
||||
FIRST_ORIGIN = pd.Timestamp("2014-01-01")
|
||||
ORIGIN_STEP_DAYS = 30
|
||||
COVERAGES = (0.5, 0.8, 0.95)
|
||||
|
||||
|
||||
def backtest(
|
||||
models, data: pd.DataFrame, horizons=HORIZONS, score_after: pd.Timestamp | None = None
|
||||
) -> pd.DataFrame:
|
||||
"""One row per (model, origin, horizon) with an observed outcome."""
|
||||
log_close = np.log(data["close"])
|
||||
last = data.index[-1]
|
||||
rows = []
|
||||
for origin in pd.date_range(FIRST_ORIGIN, last, freq=f"{ORIGIN_STEP_DAYS}D"):
|
||||
targets = origin + pd.to_timedelta(horizons, unit="D")
|
||||
scored = targets <= last
|
||||
if score_after is not None:
|
||||
scored &= targets > score_after
|
||||
if not scored.any():
|
||||
continue
|
||||
history = data.loc[:origin]
|
||||
outcome = log_close.loc[targets[scored]].to_numpy()
|
||||
for model in models:
|
||||
f = model.forecast(history, horizons[scored])
|
||||
row = {
|
||||
"model": model.name,
|
||||
"origin": origin,
|
||||
"horizon": f.horizons,
|
||||
"outcome": outcome,
|
||||
"median": f.quantile(0.5),
|
||||
"crps": crps(f.log_quantiles, outcome),
|
||||
"pit": pit(f.log_quantiles, outcome),
|
||||
}
|
||||
for c in COVERAGES:
|
||||
lo, hi = f.interval(c)
|
||||
row[f"in{c:.0%}"] = (lo <= outcome) & (outcome <= hi)
|
||||
rows.append(pd.DataFrame(row))
|
||||
return pd.concat(rows, ignore_index=True)
|
||||
|
||||
|
||||
def summarize(scores: pd.DataFrame, n_boot: int = 2000, seed: int = 0) -> pd.DataFrame:
|
||||
"""
|
||||
Per model and horizon: mean CRPS, skill relative to the baseline, and coverage.
|
||||
|
||||
Skill is 1 - CRPS / baseline CRPS (positive = better than the baseline),
|
||||
with a 90% moving-block bootstrap interval over origins. Forecasts from
|
||||
nearby origins overlap heavily, so `windows` (the span covered divided by
|
||||
the horizon) is the honest count of independent outcomes. Treat intervals
|
||||
with fewer than ~5 windows as optimistic.
|
||||
"""
|
||||
rng = np.random.default_rng(seed)
|
||||
rows = []
|
||||
for horizon, at_h in scores.groupby("horizon"):
|
||||
base = at_h[at_h.model == BASELINE].set_index("origin")["crps"].sort_index()
|
||||
span = (base.index[-1] - base.index[0]).days + horizon
|
||||
block = max(1, min(int(np.ceil(horizon / ORIGIN_STEP_DAYS)), len(base) // 2))
|
||||
boot_index = _block_bootstrap_indices(len(base), block, n_boot, rng)
|
||||
for model, g in at_h.groupby("model", sort=False):
|
||||
m = g.set_index("origin")["crps"].reindex(base.index).to_numpy()
|
||||
boot = 1 - m[boot_index].mean(axis=1) / base.to_numpy()[boot_index].mean(axis=1)
|
||||
rows.append(
|
||||
{
|
||||
"model": model,
|
||||
"horizon": horizon,
|
||||
"forecasts": len(g),
|
||||
"windows": span / horizon,
|
||||
"crps": g["crps"].mean(),
|
||||
"skill": 1 - m.mean() / base.mean(),
|
||||
"skill_lo": np.quantile(boot, 0.05),
|
||||
"skill_hi": np.quantile(boot, 0.95),
|
||||
**{f"cov{c:.0%}": g[f"in{c:.0%}"].mean() for c in COVERAGES},
|
||||
"mean_pit": g["pit"].mean(),
|
||||
}
|
||||
)
|
||||
return pd.DataFrame(rows)
|
||||
|
||||
|
||||
def _block_bootstrap_indices(n, block, n_boot, rng) -> np.ndarray:
|
||||
"""Circular block bootstrap, so the first and last origins aren't under-sampled."""
|
||||
n_blocks = int(np.ceil(n / block))
|
||||
starts = rng.integers(0, n, size=(n_boot, n_blocks))
|
||||
return ((starts[:, :, None] + np.arange(block)) % n).reshape(n_boot, -1)[:, :n]
|
||||
|
||||
|
||||
def format_summary(summary: pd.DataFrame) -> str:
|
||||
table = pd.DataFrame(
|
||||
{
|
||||
"model": summary["model"],
|
||||
"horizon": summary["horizon"].map(horizon_label),
|
||||
"windows": summary["windows"].map("{:.1f}".format),
|
||||
"crps": summary["crps"].map("{:.3f}".format),
|
||||
"skill vs rw [90%]": [
|
||||
f"{s:+.0%} [{lo:+.0%}, {hi:+.0%}]"
|
||||
for s, lo, hi in zip(summary.skill, summary.skill_lo, summary.skill_hi, strict=True)
|
||||
],
|
||||
**{f"in {c:.0%}": summary[f"cov{c:.0%}"].map("{:.0%}".format) for c in COVERAGES},
|
||||
"mean pit": summary["mean_pit"].map("{:.2f}".format),
|
||||
}
|
||||
)
|
||||
return table.to_string(index=False)
|
||||
|
||||
|
||||
def horizon_label(days: int) -> str:
|
||||
if days < 365:
|
||||
return f"{round(days / 30.4)}mo"
|
||||
return f"{round(days / 365)}y"
|
||||
@@ -0,0 +1,76 @@
|
||||
"""
|
||||
The one output format every model produces, and how it is scored.
|
||||
|
||||
A forecast is a set of quantiles of the natural-log price at each horizon.
|
||||
Quantiles work for any model (closed-form, simulated, bootstrapped, mixtures)
|
||||
and make scoring simple. Working in log price makes errors relative: a CRPS of
|
||||
0.1 is roughly "typically 10% off", in 2013 or in 2026.
|
||||
"""
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from scipy.stats import norm
|
||||
|
||||
N_LEVELS = 100
|
||||
# Midpoints of 100 equal-probability bins: 0.005, 0.015, ..., 0.995.
|
||||
LEVELS = (np.arange(N_LEVELS) + 0.5) / N_LEVELS
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Forecast:
|
||||
origin: pd.Timestamp
|
||||
horizons: np.ndarray # days after origin, shape (H,)
|
||||
log_quantiles: np.ndarray # log price at LEVELS, shape (H, N_LEVELS), rows nondecreasing
|
||||
|
||||
@classmethod
|
||||
def normal(cls, origin, horizons, mean, sd) -> "Forecast":
|
||||
"""Normal distribution in log price (i.e. lognormal price) at each horizon."""
|
||||
mean = np.broadcast_to(np.asarray(mean, dtype=float), np.shape(horizons))
|
||||
sd = np.broadcast_to(np.asarray(sd, dtype=float), np.shape(horizons))
|
||||
q = mean[:, None] + sd[:, None] * norm.ppf(LEVELS)[None, :]
|
||||
return cls(pd.Timestamp(origin), np.asarray(horizons), q)
|
||||
|
||||
@classmethod
|
||||
def from_samples(cls, origin, horizons, samples) -> "Forecast":
|
||||
"""Empirical quantiles of simulated log prices, shape (n_samples, H)."""
|
||||
q = np.quantile(np.asarray(samples), LEVELS, axis=0).T
|
||||
return cls(pd.Timestamp(origin), np.asarray(horizons), q)
|
||||
|
||||
@property
|
||||
def dates(self) -> pd.DatetimeIndex:
|
||||
return self.origin + pd.to_timedelta(self.horizons, unit="D")
|
||||
|
||||
def quantile(self, level: float) -> np.ndarray:
|
||||
"""Log price at an arbitrary level, interpolated between grid levels."""
|
||||
return np.array([np.interp(level, LEVELS, row) for row in self.log_quantiles])
|
||||
|
||||
def interval(self, coverage: float) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Central interval in log price holding `coverage` probability."""
|
||||
tail = (1 - coverage) / 2
|
||||
return self.quantile(tail), self.quantile(1 - tail)
|
||||
|
||||
|
||||
def crps(log_quantiles: np.ndarray, outcome: np.ndarray) -> np.ndarray:
|
||||
"""
|
||||
Continuous ranked probability score, from quantiles, in log-price units.
|
||||
|
||||
CRPS is twice the pinball loss integrated over all quantile levels; the
|
||||
quantile grid gives the integral directly. It rewards sharpness and
|
||||
calibration together and has no free parameters to game. Lower is better.
|
||||
"""
|
||||
outcome = np.asarray(outcome, dtype=float)
|
||||
u = outcome[..., None] - log_quantiles
|
||||
pinball = u * (LEVELS - (u < 0))
|
||||
return 2 * pinball.mean(axis=-1)
|
||||
|
||||
|
||||
def pit(log_quantiles: np.ndarray, outcome: np.ndarray) -> np.ndarray:
|
||||
"""Probability integral transform: forecast CDF evaluated at the outcome."""
|
||||
return np.array(
|
||||
[
|
||||
np.interp(y, q, LEVELS, left=0.0, right=1.0)
|
||||
for q, y in zip(np.atleast_2d(log_quantiles), np.atleast_1d(outcome), strict=True)
|
||||
]
|
||||
)
|
||||
@@ -0,0 +1,34 @@
|
||||
"""Halving calendar and position within the halving cycle."""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
GENESIS = pd.Timestamp("2009-01-03")
|
||||
|
||||
# Block heights 210k, 420k, 630k, 840k (UTC dates).
|
||||
HALVINGS = pd.DatetimeIndex(["2012-11-28", "2016-07-09", "2020-05-11", "2024-04-20"])
|
||||
|
||||
# Later halvings are projected at the length of the last cycle. Block times drift
|
||||
# by weeks per cycle, which is noise at the resolution this is used.
|
||||
_LAST_CYCLE = HALVINGS[-1] - HALVINGS[-2]
|
||||
_PROJECTED = pd.DatetimeIndex([HALVINGS[-1] + k * _LAST_CYCLE for k in range(1, 6)])
|
||||
|
||||
# Genesis starts cycle 0.
|
||||
CYCLE_STARTS = pd.DatetimeIndex([GENESIS]).append(HALVINGS).append(_PROJECTED)
|
||||
|
||||
|
||||
def cycle_position(dates) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""
|
||||
For each date, return (cycle index, days since that cycle began).
|
||||
|
||||
Cycle 0 runs from genesis to the first halving; a halving day is day 0 of
|
||||
the cycle it starts.
|
||||
"""
|
||||
dates = pd.DatetimeIndex(dates)
|
||||
if (dates < GENESIS).any():
|
||||
raise ValueError("date before genesis")
|
||||
if (dates >= CYCLE_STARTS[-1]).any():
|
||||
raise ValueError("date beyond projected halvings")
|
||||
index = CYCLE_STARTS.searchsorted(dates, side="right") - 1
|
||||
days = (dates - CYCLE_STARTS[index]).days
|
||||
return np.asarray(index), np.asarray(days)
|
||||
@@ -0,0 +1,15 @@
|
||||
"""
|
||||
Candidate models.
|
||||
|
||||
A model is any object with a `name` and a
|
||||
`forecast(history: pd.DataFrame, horizons: np.ndarray) -> Forecast` method.
|
||||
`history` holds every row up to and including the forecast origin and nothing
|
||||
after it; the harness guarantees that, so models can use all of it freely.
|
||||
"""
|
||||
|
||||
from .baselines import DriftRandomWalk, RandomWalk
|
||||
from .cycle import CycleModel
|
||||
|
||||
# Order is fixed: it sets each model's colour in every chart.
|
||||
MODELS = {m.name: m for m in (RandomWalk(), DriftRandomWalk(), CycleModel())}
|
||||
BASELINE = RandomWalk.name
|
||||
@@ -0,0 +1,47 @@
|
||||
"""Reference forecasts every other model has to beat."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import ClassVar
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from ..data import log_returns
|
||||
from ..forecast import Forecast
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class RandomWalk:
|
||||
"""
|
||||
Zero-drift random walk in log price: "it stays about here, give or take".
|
||||
|
||||
Volatility is the trailing standard deviation of daily log returns.
|
||||
"""
|
||||
|
||||
name: ClassVar[str] = "random_walk"
|
||||
vol_window: int = 365
|
||||
|
||||
def forecast(self, history: pd.DataFrame, horizons: np.ndarray) -> Forecast:
|
||||
sigma = log_returns(history).iloc[-self.vol_window :].std()
|
||||
mean = np.log(history["close"].iloc[-1])
|
||||
return Forecast.normal(history.index[-1], horizons, mean, sigma * np.sqrt(horizons))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class DriftRandomWalk:
|
||||
"""
|
||||
Random walk whose drift is the mean daily log return over the trailing
|
||||
`drift_window` days (one halving cycle by default): "it keeps doing what it
|
||||
did last cycle".
|
||||
"""
|
||||
|
||||
name: ClassVar[str] = "drift_rw"
|
||||
drift_window: int = 1460
|
||||
vol_window: int = 365
|
||||
|
||||
def forecast(self, history: pd.DataFrame, horizons: np.ndarray) -> Forecast:
|
||||
returns = log_returns(history)
|
||||
mu = returns.iloc[-self.drift_window :].mean()
|
||||
sigma = returns.iloc[-self.vol_window :].std()
|
||||
mean = np.log(history["close"].iloc[-1]) + mu * horizons
|
||||
return Forecast.normal(history.index[-1], horizons, mean, sigma * np.sqrt(horizons))
|
||||
@@ -0,0 +1,73 @@
|
||||
"""The 2024 model, distilled."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import ClassVar
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from ..data import log_returns
|
||||
from ..forecast import Forecast
|
||||
from ..halving import cycle_position
|
||||
|
||||
# Longer than any cycle so far (the longest, cycle 0, is 1425 days).
|
||||
MAX_CYCLE_DAYS = 1500
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class CycleModel:
|
||||
"""
|
||||
Expected return depends on how many days it has been since the last halving.
|
||||
|
||||
The drift for day d of the cycle is a weighted mean of the daily log returns
|
||||
observed around day d of every past cycle. Of the 2024 model's ~2000 lines,
|
||||
this idea did all the work. It differs from that model in two ways:
|
||||
|
||||
- Neighbouring cycle days are pooled with a Gaussian kernel. The 2024 model
|
||||
averaged each day separately and then took a rolling mean.
|
||||
- Past cycles are down-weighted, halving each `recency_half_life` cycles, so
|
||||
the 10-100x cycles of 2011-2017 don't set the level. The 2024 model
|
||||
averaged all cycles equally, then scaled by ~0.7; it overshot the
|
||||
2025 peak by ~60%. (This is the idea on the old `tuning-b` branch.)
|
||||
|
||||
Where the data is thin, the drift shrinks toward the overall weighted mean,
|
||||
as if `prior_days` extra observations sat at that value. Noise is a
|
||||
constant-volatility random walk, so the distribution is closed-form.
|
||||
"""
|
||||
|
||||
name: ClassVar[str] = "cycle"
|
||||
bandwidth_days: float = 30.0
|
||||
recency_half_life: float = 1.0
|
||||
prior_days: float = 10.0
|
||||
vol_window: int = 365
|
||||
|
||||
def drift_by_cycle_day(self, history: pd.DataFrame) -> np.ndarray:
|
||||
"""Expected daily log return for each day of the cycle, shape (MAX_CYCLE_DAYS,)."""
|
||||
returns = log_returns(history)
|
||||
cycle, day = cycle_position(returns.index)
|
||||
current_cycle = cycle_position(history.index[-1:])[0][0]
|
||||
weight = 0.5 ** ((current_cycle - cycle) / self.recency_half_life)
|
||||
|
||||
sum_wr = np.bincount(day, weights=weight * returns.values, minlength=MAX_CYCLE_DAYS)
|
||||
sum_w = np.bincount(day, weights=weight, minlength=MAX_CYCLE_DAYS)
|
||||
|
||||
# Peak-1 kernel, so smoothed weights count (recency-weighted) days of data.
|
||||
half_width = int(np.ceil(4 * self.bandwidth_days))
|
||||
offsets = np.arange(-half_width, half_width + 1)
|
||||
kernel = np.exp(-0.5 * (offsets / self.bandwidth_days) ** 2)
|
||||
smooth_wr = np.convolve(sum_wr, kernel, mode="same")
|
||||
smooth_w = np.convolve(sum_w, kernel, mode="same")
|
||||
|
||||
overall = sum_wr.sum() / sum_w.sum()
|
||||
return (smooth_wr + self.prior_days * overall) / (smooth_w + self.prior_days)
|
||||
|
||||
def forecast(self, history: pd.DataFrame, horizons: np.ndarray) -> Forecast:
|
||||
drift = self.drift_by_cycle_day(history)
|
||||
origin = history.index[-1]
|
||||
future = origin + pd.to_timedelta(np.arange(1, horizons.max() + 1), unit="D")
|
||||
_, future_day = cycle_position(future)
|
||||
cumulative = np.cumsum(drift[future_day])
|
||||
|
||||
sigma = log_returns(history).iloc[-self.vol_window :].std()
|
||||
mean = np.log(history["close"].iloc[-1]) + cumulative[horizons - 1]
|
||||
return Forecast.normal(origin, horizons, mean, sigma * np.sqrt(horizons))
|
||||
@@ -0,0 +1,172 @@
|
||||
"""Charts. One colour per model, fixed by its position in MODELS."""
|
||||
|
||||
from pathlib import Path
|
||||
|
||||
import matplotlib
|
||||
|
||||
matplotlib.use("Agg")
|
||||
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
import numpy as np # noqa: E402
|
||||
import pandas as pd # noqa: E402
|
||||
from matplotlib.ticker import FuncFormatter, LogLocator, PercentFormatter # noqa: E402
|
||||
|
||||
from .evaluate import COVERAGES, horizon_label # noqa: E402
|
||||
from .forecast import Forecast # noqa: E402
|
||||
from .halving import CYCLE_STARTS # noqa: E402
|
||||
from .models import BASELINE, MODELS # noqa: E402
|
||||
|
||||
SURFACE = "#fcfcfb"
|
||||
INK = "#0b0b0b"
|
||||
INK_SECONDARY = "#52514e"
|
||||
MUTED = "#898781"
|
||||
GRID = "#e1e0d9"
|
||||
AXIS = "#c3c2b7"
|
||||
SERIES = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#4a3aa7", "#e34948"]
|
||||
|
||||
plt.rcParams.update(
|
||||
{
|
||||
"figure.facecolor": SURFACE,
|
||||
"axes.facecolor": SURFACE,
|
||||
"savefig.facecolor": SURFACE,
|
||||
"axes.edgecolor": AXIS,
|
||||
"axes.linewidth": 0.8,
|
||||
"axes.labelcolor": INK_SECONDARY,
|
||||
"axes.titlecolor": INK,
|
||||
"axes.titlesize": 11,
|
||||
"axes.titleweight": "semibold",
|
||||
"axes.titlelocation": "left",
|
||||
"axes.spines.top": False,
|
||||
"axes.spines.right": False,
|
||||
"axes.grid": True,
|
||||
"grid.color": GRID,
|
||||
"grid.linewidth": 0.8,
|
||||
"xtick.color": MUTED,
|
||||
"ytick.color": MUTED,
|
||||
"xtick.labelcolor": INK_SECONDARY,
|
||||
"ytick.labelcolor": INK_SECONDARY,
|
||||
"legend.frameon": False,
|
||||
"legend.labelcolor": INK_SECONDARY,
|
||||
"lines.linewidth": 2,
|
||||
"lines.solid_capstyle": "round",
|
||||
"lines.solid_joinstyle": "round",
|
||||
"font.family": "sans-serif",
|
||||
"font.size": 9,
|
||||
}
|
||||
)
|
||||
|
||||
|
||||
def model_color(name: str) -> str:
|
||||
return SERIES[list(MODELS).index(name) % len(SERIES)]
|
||||
|
||||
|
||||
def price_formatter(x, _=None) -> str:
|
||||
for scale, suffix in ((1e9, "B"), (1e6, "M"), (1e3, "K")):
|
||||
if x >= scale:
|
||||
return f"${x / scale:.3g}{suffix}"
|
||||
return f"${x:.3g}"
|
||||
|
||||
|
||||
def fan_chart(history: pd.DataFrame, forecasts: dict[str, Forecast], path: Path) -> None:
|
||||
"""History plus 50/80/95% intervals, one panel per model on shared axes."""
|
||||
fig, axes = plt.subplots(
|
||||
len(forecasts), 1, figsize=(10, 3.2 * len(forecasts)), sharex=True, sharey=True
|
||||
)
|
||||
axes = np.atleast_1d(axes)
|
||||
shown = history[history.index >= history.index[-1] - pd.Timedelta(days=6 * 365)]
|
||||
for ax, (name, f) in zip(axes, forecasts.items(), strict=True):
|
||||
color = model_color(name)
|
||||
for c in sorted(COVERAGES, reverse=True):
|
||||
lo, hi = f.interval(c)
|
||||
ax.fill_between(f.dates, np.exp(lo), np.exp(hi), color=color, alpha=0.1, lw=0)
|
||||
median = np.exp(f.quantile(0.5))
|
||||
ax.plot(shown.index, shown["close"], color=INK, lw=1.2)
|
||||
ax.plot(f.dates, median, color=color)
|
||||
ax.annotate(
|
||||
f"median {price_formatter(median[-1])}",
|
||||
(f.dates[-1], median[-1]),
|
||||
xytext=(6, 0),
|
||||
textcoords="offset points",
|
||||
va="center",
|
||||
color=INK_SECONDARY,
|
||||
)
|
||||
for start in CYCLE_STARTS:
|
||||
if shown.index[0] <= start <= f.dates[-1]:
|
||||
ax.axvline(start, color=AXIS, lw=0.8, zorder=0)
|
||||
ax.set_yscale("log")
|
||||
ax.yaxis.set_major_locator(LogLocator(base=10, subs=(1, 2, 5)))
|
||||
ax.yaxis.set_major_formatter(FuncFormatter(price_formatter))
|
||||
ax.yaxis.set_minor_formatter(FuncFormatter(lambda *_: ""))
|
||||
ax.set_title(f"{name}: forecast from {f.origin:%Y-%m-%d}")
|
||||
axes[0].fill_between([], [], color=MUTED, alpha=0.3, label="50% interval")
|
||||
axes[0].fill_between([], [], color=MUTED, alpha=0.2, label="80% interval")
|
||||
axes[0].fill_between([], [], color=MUTED, alpha=0.1, label="95% interval")
|
||||
axes[0].legend(loc="upper left")
|
||||
axes[-1].set_xlabel("vertical lines: halvings (future ones projected)", color=MUTED)
|
||||
fig.tight_layout()
|
||||
fig.savefig(path, dpi=150)
|
||||
plt.close(fig)
|
||||
|
||||
|
||||
def skill_chart(summary: pd.DataFrame, path: Path, title: str) -> None:
|
||||
"""CRPS skill vs the random walk, by horizon, with bootstrap intervals."""
|
||||
horizons = sorted(summary["horizon"].unique())
|
||||
x = np.arange(len(horizons))
|
||||
fig, ax = plt.subplots(figsize=(8, 4.5))
|
||||
ax.axhline(0, color=model_color(BASELINE), lw=2, label=BASELINE)
|
||||
for name, g in summary[summary.model != BASELINE].groupby("model", sort=False):
|
||||
g = g.set_index("horizon").reindex(horizons)
|
||||
color = model_color(name)
|
||||
ax.fill_between(x, g.skill_lo, g.skill_hi, color=color, alpha=0.1, lw=0)
|
||||
ax.plot(x, g.skill, color=color, marker="o", ms=6, mec=SURFACE, mew=2, label=name)
|
||||
ax.annotate(
|
||||
name,
|
||||
(x[-1], g.skill.iloc[-1]),
|
||||
xytext=(8, 0),
|
||||
textcoords="offset points",
|
||||
va="center",
|
||||
color=INK_SECONDARY,
|
||||
)
|
||||
ax.set_xticks(x, [horizon_label(h) for h in horizons])
|
||||
ax.set_xlabel("forecast horizon")
|
||||
ax.set_ylabel("CRPS skill vs random walk (higher is better)")
|
||||
ax.yaxis.set_major_formatter(PercentFormatter(1.0, decimals=0))
|
||||
ax.set_title(title)
|
||||
ax.legend(loc="lower left")
|
||||
fig.tight_layout()
|
||||
fig.savefig(path, dpi=150)
|
||||
plt.close(fig)
|
||||
|
||||
|
||||
def calibration_chart(summary: pd.DataFrame, path: Path, title: str) -> None:
|
||||
"""How often each nominal interval contained the outcome, by horizon."""
|
||||
horizons = sorted(summary["horizon"].unique())
|
||||
x = np.arange(len(horizons))
|
||||
fig, axes = plt.subplots(1, len(COVERAGES), figsize=(12, 4), sharey=True)
|
||||
for ax, c in zip(axes, COVERAGES, strict=True):
|
||||
ax.axhline(c, color=INK_SECONDARY, lw=1)
|
||||
ax.annotate(
|
||||
"target", (x[-1], c), xytext=(0, 4), textcoords="offset points", ha="right", color=MUTED
|
||||
)
|
||||
for name, g in summary.groupby("model", sort=False):
|
||||
g = g.set_index("horizon").reindex(horizons)
|
||||
ax.plot(
|
||||
x,
|
||||
g[f"cov{c:.0%}"],
|
||||
color=model_color(name),
|
||||
marker="o",
|
||||
ms=6,
|
||||
mec=SURFACE,
|
||||
mew=2,
|
||||
label=name,
|
||||
)
|
||||
ax.set_xticks(x, [horizon_label(h) for h in horizons])
|
||||
ax.set_title(f"{c:.0%} interval")
|
||||
ax.set_ylim(0, 1.02)
|
||||
ax.yaxis.set_major_formatter(PercentFormatter(1.0, decimals=0))
|
||||
axes[0].set_ylabel("share of outcomes inside")
|
||||
axes[-1].legend(loc="lower left")
|
||||
fig.suptitle(title, x=0.01, ha="left", color=INK, fontsize=11, fontweight="semibold")
|
||||
fig.tight_layout()
|
||||
fig.savefig(path, dpi=150)
|
||||
plt.close(fig)
|
||||
Reference in New Issue
Block a user