Composable models and A/B tests of the 2024 ideas; add powerlaw.
Models are now a Composite of drift, volatility and (optional) shape components, so an experiment can swap one part against a fixed control. btcmodel/experiments.py holds seven experiments built from the ideas in the old branches (catalogued in docs/2024-ideas.md), each with its hypothesis and source, and a verdict rule fixed before anything ran. `just ab` runs them on development data. Results: - Shrinking the cycle drift, and a power-law trend (plain or reverting), beat their controls. The power law beats the random walk by 53-63% at 3-4 years with unbiased outcomes, so it is promoted to MODELS. - Every alternative volatility estimate (EWMA blends, other windows, reversion to a level or trend) is worse than the trailing 365-day window. Cycle-dependent volatility, heavy tails and stretched cycle phase show no reliable effect.
This commit is contained in:
+43
-1
@@ -5,6 +5,7 @@ Command line entry point.
|
||||
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
|
||||
python -m btcmodel ab [NAME ...] run A/B experiments on development data
|
||||
"""
|
||||
|
||||
import argparse
|
||||
@@ -13,7 +14,7 @@ from pathlib import Path
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from . import data, evaluate, plots
|
||||
from . import data, evaluate, experiments, plots
|
||||
from .models import MODELS
|
||||
|
||||
FORECAST_REPORT_HORIZONS = (182, 365, 730, 1095, 1460)
|
||||
@@ -33,7 +34,11 @@ def main() -> None:
|
||||
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")
|
||||
ab = commands.add_parser("ab", help="run A/B experiments on development data")
|
||||
ab.add_argument("names", nargs="*", help="experiments to run (default: all)")
|
||||
args = parser.parse_args()
|
||||
if args.command == "ab" and (unknown := set(args.names) - set(experiments.EXPERIMENTS)):
|
||||
parser.error(f"unknown experiments: {', '.join(sorted(unknown))}")
|
||||
models = [MODELS[name] for name in args.models]
|
||||
|
||||
if args.command == "update":
|
||||
@@ -43,6 +48,8 @@ def main() -> None:
|
||||
run_backtest(models, args.output, args.holdout)
|
||||
elif args.command == "forecast":
|
||||
run_forecast(models, args.output)
|
||||
elif args.command == "ab":
|
||||
run_ab(args.names or list(experiments.EXPERIMENTS), args.output)
|
||||
|
||||
|
||||
def run_backtest(models, output: Path, holdout: bool) -> None:
|
||||
@@ -70,6 +77,41 @@ def run_backtest(models, output: Path, holdout: bool) -> None:
|
||||
print(f"\nwrote {out}/")
|
||||
|
||||
|
||||
def run_ab(names: list[str], output: Path) -> None:
|
||||
prices = data.load_prices(until=data.DEV_CUTOFF)
|
||||
verdicts = []
|
||||
for name in names:
|
||||
experiment = experiments.EXPERIMENTS[name]
|
||||
scores, summary = experiments.run(experiment, prices)
|
||||
out = output / "ab" / name
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
control = experiment.control.name
|
||||
report = (
|
||||
f"{name}: {experiment.hypothesis}\n(from {experiment.source})\n\n"
|
||||
+ evaluate.format_summary(summary, baseline=control)
|
||||
)
|
||||
print(report + "\n")
|
||||
(out / "report.txt").write_text(report + "\n")
|
||||
summary.to_csv(out / "summary.csv", index=False)
|
||||
plots.skill_chart(summary, out / "skill.png", f"{name}: skill vs {control}", control)
|
||||
plots.calibration_chart(summary, out / "calibration.png", f"{name}: interval coverage")
|
||||
for variant in experiment.variants:
|
||||
v = summary[summary.model == variant.name].set_index("horizon")
|
||||
verdicts.append(
|
||||
{
|
||||
"experiment": name,
|
||||
"variant": variant.name,
|
||||
"control": control,
|
||||
**{evaluate.horizon_label(h): f"{s:+.0%}" for h, s in v["skill"].items()},
|
||||
"verdict": experiments.verdict(v),
|
||||
}
|
||||
)
|
||||
table = pd.DataFrame(verdicts).to_string(index=False)
|
||||
print(table)
|
||||
(output / "ab").mkdir(parents=True, exist_ok=True)
|
||||
(output / "ab" / "verdicts.txt").write_text(table + "\n")
|
||||
|
||||
|
||||
def run_forecast(models, output: Path) -> None:
|
||||
prices = data.load_prices()
|
||||
horizons = np.arange(1, max(FORECAST_REPORT_HORIZONS) + 1)
|
||||
|
||||
@@ -55,9 +55,11 @@ def backtest(
|
||||
return pd.concat(rows, ignore_index=True)
|
||||
|
||||
|
||||
def summarize(scores: pd.DataFrame, n_boot: int = 2000, seed: int = 0) -> pd.DataFrame:
|
||||
def summarize(
|
||||
scores: pd.DataFrame, baseline: str = BASELINE, n_boot: int = 2000, seed: int = 0
|
||||
) -> pd.DataFrame:
|
||||
"""
|
||||
Per model and horizon: mean CRPS, skill relative to the baseline, and coverage.
|
||||
Per model and horizon: mean CRPS, skill relative to `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
|
||||
@@ -68,7 +70,7 @@ def summarize(scores: pd.DataFrame, n_boot: int = 2000, seed: int = 0) -> pd.Dat
|
||||
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()
|
||||
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)
|
||||
@@ -99,14 +101,14 @@ def _block_bootstrap_indices(n, block, n_boot, rng) -> np.ndarray:
|
||||
return ((starts[:, :, None] + np.arange(block)) % n).reshape(n_boot, -1)[:, :n]
|
||||
|
||||
|
||||
def format_summary(summary: pd.DataFrame) -> str:
|
||||
def format_summary(summary: pd.DataFrame, baseline: str = BASELINE) -> 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"skill vs {baseline} [90%]": [
|
||||
f"{s:+.0%} [{lo:+.0%}, {hi:+.0%}]"
|
||||
for s, lo, hi in zip(summary.skill, summary.skill_lo, summary.skill_hi, strict=True)
|
||||
],
|
||||
|
||||
@@ -0,0 +1,142 @@
|
||||
"""
|
||||
A/B tests of individual ideas, mostly salvaged from the 2024 model's branches.
|
||||
|
||||
Each experiment states its hypothesis and pairs a control with variants that
|
||||
change one component. Experiments are written down before they are run, and
|
||||
every result is reported, including the failures; with a few dozen
|
||||
comparisons over a handful of independent windows, some "wins" will be luck.
|
||||
Verdicts use one rule, fixed in advance (see `verdict`), and anything that
|
||||
passes still has to hold up on the holdout.
|
||||
"""
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import pandas as pd
|
||||
|
||||
from .evaluate import backtest, summarize
|
||||
from .models import MODELS
|
||||
from .models.base import Composite
|
||||
from .models.drift import (
|
||||
CycleDrift,
|
||||
PowerLawDrift,
|
||||
PowerLawScaledCycleDrift,
|
||||
ShrunkDrift,
|
||||
TrailingMeanDrift,
|
||||
)
|
||||
from .models.shape import Empirical, StudentT
|
||||
from .models.volatility import CycleVol, EwmaVol, ReversionVol, TrailingVol
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Experiment:
|
||||
name: str
|
||||
hypothesis: str
|
||||
source: str # where in the old history the idea came from
|
||||
control: Composite
|
||||
variants: tuple[Composite, ...]
|
||||
|
||||
@property
|
||||
def models(self) -> list[Composite]:
|
||||
return [self.control, *self.variants]
|
||||
|
||||
|
||||
def run(experiment: Experiment, prices: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]:
|
||||
scores = backtest(experiment.models, prices)
|
||||
return scores, summarize(scores, baseline=experiment.control.name)
|
||||
|
||||
|
||||
def verdict(variant_summary: pd.DataFrame) -> str:
|
||||
"""
|
||||
Decide from skill vs the control, per horizon, with 90% intervals:
|
||||
|
||||
- "worse" if the interval is entirely below zero at any horizon;
|
||||
- "better" if the interval is entirely above zero at two or more horizons
|
||||
and the point estimate is non-negative at every horizon;
|
||||
- "inconclusive" otherwise.
|
||||
"""
|
||||
if (variant_summary["skill_hi"] < 0).any():
|
||||
return "worse"
|
||||
if (variant_summary["skill_lo"] > 0).sum() >= 2 and (variant_summary["skill"] >= 0).all():
|
||||
return "better"
|
||||
return "inconclusive"
|
||||
|
||||
|
||||
# Idea labels (D1, V2, ...) refer to docs/2024-ideas.md.
|
||||
CYCLE = MODELS["cycle"]
|
||||
DRIFT_RW = MODELS["drift_rw"]
|
||||
# Volatility and shape experiments use drift_rw as the control: the zero-drift
|
||||
# random walk is biased low at long horizons, so anything that merely widened
|
||||
# its intervals would look like an improvement.
|
||||
_VOL = dict(drift=TrailingMeanDrift())
|
||||
|
||||
EXPERIMENTS: dict[str, Experiment] = {
|
||||
e.name: e
|
||||
for e in (
|
||||
Experiment(
|
||||
"shrink-cycle",
|
||||
"the cycle drift is overfit; pulling it toward zero improves it",
|
||||
"D2: damping constants throughout the 2024 model",
|
||||
CYCLE,
|
||||
tuple(
|
||||
Composite(f"cycle_x{f}", ShrunkDrift(CycleDrift(), f), TrailingVol())
|
||||
for f in (0.25, 0.5, 0.75)
|
||||
),
|
||||
),
|
||||
Experiment(
|
||||
"diminishing-returns",
|
||||
"each cycle grows less than the last; a power-law trend captures that",
|
||||
"D3: initial commit, old/backtests-trend-enhancement-1",
|
||||
DRIFT_RW,
|
||||
(
|
||||
Composite("powerlaw", PowerLawDrift(), TrailingVol()),
|
||||
Composite("powerlaw_revert", PowerLawDrift(revert=True), TrailingVol()),
|
||||
Composite("cycle_on_powerlaw", PowerLawScaledCycleDrift(), TrailingVol()),
|
||||
),
|
||||
),
|
||||
Experiment(
|
||||
"vol-window",
|
||||
"recent volatility predicts future volatility better than a flat 365-day window",
|
||||
"V1: 'Add improved vol calculation'",
|
||||
DRIFT_RW,
|
||||
(
|
||||
Composite("ewma_blend", volatility=EwmaVol(), **_VOL),
|
||||
Composite("ewma_90", volatility=EwmaVol(spans=(90,), weights=(1.0,)), **_VOL),
|
||||
Composite("trailing_730", volatility=TrailingVol(730), **_VOL),
|
||||
),
|
||||
),
|
||||
Experiment(
|
||||
"vol-reversion",
|
||||
"volatility shocks fade toward a long-run level, which itself is falling",
|
||||
"V2, V4, C2: adaptive windows, market maturity, horizon multipliers",
|
||||
DRIFT_RW,
|
||||
(
|
||||
Composite("revert_level", volatility=ReversionVol(), **_VOL),
|
||||
Composite("revert_trend", volatility=ReversionVol(long_run="trend"), **_VOL),
|
||||
),
|
||||
),
|
||||
Experiment(
|
||||
"cycle-vol",
|
||||
"volatility depends on position in the halving cycle",
|
||||
"V3: old/backtest-vol-cycle",
|
||||
DRIFT_RW,
|
||||
(Composite("cycle_vol", volatility=CycleVol(), **_VOL),),
|
||||
),
|
||||
Experiment(
|
||||
"tails",
|
||||
"log returns are heavier-tailed than normal, even at long horizons",
|
||||
"S1: skewed innovations in tuning-1",
|
||||
DRIFT_RW,
|
||||
(
|
||||
Composite("student_t4", volatility=TrailingVol(), shape=StudentT(4), **_VOL),
|
||||
Composite("empirical", volatility=TrailingVol(), shape=Empirical(), **_VOL),
|
||||
),
|
||||
),
|
||||
Experiment(
|
||||
"cycle-phase",
|
||||
"cycles align better by fraction elapsed than by days since the halving",
|
||||
"D1: the 2024 model stretched every cycle to 1460 days",
|
||||
CYCLE,
|
||||
(Composite("cycle_fraction", CycleDrift(phase="fraction"), TrailingVol()),),
|
||||
),
|
||||
)
|
||||
}
|
||||
@@ -32,3 +32,10 @@ def cycle_position(dates) -> tuple[np.ndarray, np.ndarray]:
|
||||
index = CYCLE_STARTS.searchsorted(dates, side="right") - 1
|
||||
days = (dates - CYCLE_STARTS[index]).days
|
||||
return np.asarray(index), np.asarray(days)
|
||||
|
||||
|
||||
def cycle_fraction(dates) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""For each date, return (cycle index, fraction of that cycle elapsed, in [0, 1))."""
|
||||
index, days = cycle_position(dates)
|
||||
lengths = (CYCLE_STARTS[index + 1] - CYCLE_STARTS[index]).days
|
||||
return index, days / np.asarray(lengths)
|
||||
|
||||
@@ -5,11 +5,26 @@ 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.
|
||||
Most models are a Composite of a drift and a volatility component.
|
||||
"""
|
||||
|
||||
from .baselines import DriftRandomWalk, RandomWalk
|
||||
from .cycle import CycleModel
|
||||
from .base import Composite
|
||||
from .drift import CycleDrift, PowerLawDrift, TrailingMeanDrift, ZeroDrift
|
||||
from .volatility import TrailingVol
|
||||
|
||||
BASELINE = "random_walk"
|
||||
|
||||
# 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
|
||||
MODELS = {
|
||||
m.name: m
|
||||
for m in (
|
||||
# "It stays about here, give or take."
|
||||
Composite(BASELINE, ZeroDrift(), TrailingVol()),
|
||||
# "It keeps doing what it did last cycle."
|
||||
Composite("drift_rw", TrailingMeanDrift(), TrailingVol()),
|
||||
# The 2024 model, distilled.
|
||||
Composite("cycle", CycleDrift(), TrailingVol()),
|
||||
# "Growth keeps slowing, like it always has." Passed the diminishing-returns A/B.
|
||||
Composite("powerlaw", PowerLawDrift(), TrailingVol()),
|
||||
)
|
||||
}
|
||||
|
||||
@@ -0,0 +1,50 @@
|
||||
"""
|
||||
Models built from parts.
|
||||
|
||||
Most ideas about Bitcoin prices say something about either the expected return
|
||||
(drift) or the size of the uncertainty (volatility). A Composite pairs one of
|
||||
each, so an A/B test can swap exactly one part and hold the other fixed.
|
||||
"""
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import Protocol
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from ..forecast import Forecast
|
||||
from .shape import Normal
|
||||
|
||||
|
||||
class Drift(Protocol):
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
"""Expected cumulative log return from the origin to each horizon."""
|
||||
...
|
||||
|
||||
|
||||
class Volatility(Protocol):
|
||||
def sd(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
"""Standard deviation of the cumulative log return at each horizon."""
|
||||
...
|
||||
|
||||
|
||||
class Shape(Protocol):
|
||||
def standard_quantiles(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
"""Quantiles at LEVELS of a mean-0, sd-1 distribution, shape (H, N_LEVELS)."""
|
||||
...
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Composite:
|
||||
"""Log price: mean from `drift`, spread from `volatility`, normal unless `shape` says."""
|
||||
|
||||
name: str
|
||||
drift: Drift
|
||||
volatility: Volatility
|
||||
shape: Shape = Normal()
|
||||
|
||||
def forecast(self, history: pd.DataFrame, horizons: np.ndarray) -> Forecast:
|
||||
mean = np.log(history["close"].iloc[-1]) + self.drift.expected_log_return(history, horizons)
|
||||
sd = self.volatility.sd(history, horizons)
|
||||
z = self.shape.standard_quantiles(history, horizons)
|
||||
return Forecast(history.index[-1], horizons, mean[:, None] + sd[:, None] * z)
|
||||
@@ -1,47 +0,0 @@
|
||||
"""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))
|
||||
@@ -1,73 +0,0 @@
|
||||
"""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,175 @@
|
||||
"""Expected-return components."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from ..data import log_returns
|
||||
from ..halving import GENESIS, cycle_fraction, cycle_position
|
||||
|
||||
# Longer than any cycle so far (the longest, cycle 0, is 1425 days).
|
||||
MAX_CYCLE_DAYS = 1500
|
||||
# Length every cycle is stretched to when phase="fraction".
|
||||
NOMINAL_CYCLE_DAYS = 1440
|
||||
|
||||
|
||||
def mean_by_cycle_day(values, day, weight, bandwidth_days, prior_days) -> np.ndarray:
|
||||
"""
|
||||
Weighted mean of `values` around each day of the cycle, shape (MAX_CYCLE_DAYS,).
|
||||
|
||||
Neighbouring days are pooled with a Gaussian kernel. Where the data is thin,
|
||||
the mean shrinks toward the overall weighted mean, as if `prior_days` extra
|
||||
observations sat at that value.
|
||||
"""
|
||||
sum_wv = np.bincount(day, weights=weight * values, minlength=MAX_CYCLE_DAYS)
|
||||
sum_w = np.bincount(day, weights=weight, minlength=MAX_CYCLE_DAYS)
|
||||
|
||||
# Peak-1 kernel, so smoothed weights count (weighted) days of data.
|
||||
half_width = int(np.ceil(4 * bandwidth_days))
|
||||
offsets = np.arange(-half_width, half_width + 1)
|
||||
kernel = np.exp(-0.5 * (offsets / bandwidth_days) ** 2)
|
||||
smooth_wv = np.convolve(sum_wv, kernel, mode="same")
|
||||
smooth_w = np.convolve(sum_w, kernel, mode="same")
|
||||
|
||||
overall = sum_wv.sum() / sum_w.sum()
|
||||
return (smooth_wv + prior_days * overall) / (smooth_w + prior_days)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class ZeroDrift:
|
||||
"""No expected change in log price."""
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
return np.zeros(len(horizons))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TrailingMeanDrift:
|
||||
"""Mean daily log return over the trailing `window` days, extrapolated."""
|
||||
|
||||
window: int = 1460
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
return log_returns(history).iloc[-self.window :].mean() * horizons
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class CycleDrift:
|
||||
"""
|
||||
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 mean (see
|
||||
`mean_by_cycle_day`).
|
||||
|
||||
`phase="days"` aligns cycles by days since the halving; `phase="fraction"`
|
||||
stretches every cycle to NOMINAL_CYCLE_DAYS, as the 2024 model did.
|
||||
"""
|
||||
|
||||
bandwidth_days: float = 30.0
|
||||
recency_half_life: float = 1.0
|
||||
prior_days: float = 10.0
|
||||
phase: str = "days"
|
||||
|
||||
def position(self, dates) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""(cycle index, cycle day) under this model's phase convention."""
|
||||
if self.phase == "days":
|
||||
return cycle_position(dates)
|
||||
if self.phase == "fraction":
|
||||
index, fraction = cycle_fraction(dates)
|
||||
return index, np.round(fraction * NOMINAL_CYCLE_DAYS).astype(int)
|
||||
raise ValueError(f"unknown phase {self.phase!r}")
|
||||
|
||||
def 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 = self.position(returns.index)
|
||||
current_cycle = self.position(history.index[-1:])[0][0]
|
||||
return mean_by_cycle_day(
|
||||
returns.to_numpy(),
|
||||
day,
|
||||
weight=0.5 ** ((current_cycle - cycle) / self.recency_half_life),
|
||||
bandwidth_days=self.bandwidth_days,
|
||||
prior_days=self.prior_days,
|
||||
)
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
drift = self.by_cycle_day(history)
|
||||
future = history.index[-1] + pd.to_timedelta(np.arange(1, horizons.max() + 1), unit="D")
|
||||
_, future_day = self.position(future)
|
||||
return np.cumsum(drift[future_day])[horizons - 1]
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class ShrunkDrift:
|
||||
"""Another drift scaled by `factor`: 0 is no drift, 1 is the original."""
|
||||
|
||||
inner: CycleDrift | TrailingMeanDrift
|
||||
factor: float
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
return self.factor * self.inner.expected_log_return(history, horizons)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PowerLawDrift:
|
||||
"""
|
||||
Log price grows linearly in log time since genesis: ln P = a + b ln t.
|
||||
|
||||
Growth therefore slows like b/t, which is the diminishing returns the
|
||||
cycle-by-cycle numbers show. The line is fitted by least squares to the
|
||||
history. With `revert=True` the gap between price and line also closes,
|
||||
at the rate of an AR(1) fitted to the daily residuals.
|
||||
"""
|
||||
|
||||
revert: bool = False
|
||||
|
||||
def fit(self, history: pd.DataFrame) -> tuple[float, float, float]:
|
||||
"""Return (intercept, slope, daily AR(1) coefficient of the residuals)."""
|
||||
log_t = np.log((history.index - GENESIS).days.to_numpy())
|
||||
log_p = np.log(history["close"].to_numpy())
|
||||
slope, intercept = np.polyfit(log_t, log_p, 1)
|
||||
resid = log_p - (intercept + slope * log_t)
|
||||
phi = resid[1:] @ resid[:-1] / (resid[:-1] @ resid[:-1])
|
||||
return intercept, slope, phi
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
intercept, slope, phi = self.fit(history)
|
||||
t0 = (history.index[-1] - GENESIS).days
|
||||
trend = slope * np.log((t0 + horizons) / t0)
|
||||
if not self.revert:
|
||||
return trend
|
||||
gap = np.log(history["close"].iloc[-1]) - (intercept + slope * np.log(t0))
|
||||
return trend + (phi**horizons - 1) * gap
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PowerLawScaledCycleDrift:
|
||||
"""
|
||||
The cycle shape, rescaled so its average over a cycle equals the power-law
|
||||
trend's growth rate over the next cycle: keep the timing, fix the level.
|
||||
"""
|
||||
|
||||
cycle: CycleDrift = CycleDrift()
|
||||
|
||||
def expected_log_return(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
shape = self.cycle.by_cycle_day(history)
|
||||
cycle_mean = shape[:NOMINAL_CYCLE_DAYS].mean()
|
||||
level = PowerLawDrift().expected_log_return(history, np.array([NOMINAL_CYCLE_DAYS]))[0]
|
||||
level /= NOMINAL_CYCLE_DAYS
|
||||
if cycle_mean <= 1e-6:
|
||||
return level * horizons
|
||||
future = history.index[-1] + pd.to_timedelta(np.arange(1, horizons.max() + 1), unit="D")
|
||||
_, future_day = self.cycle.position(future)
|
||||
return np.cumsum(shape[future_day] * level / cycle_mean)[horizons - 1]
|
||||
@@ -0,0 +1,52 @@
|
||||
"""Distribution shapes: standardised quantiles (mean 0, sd 1) at each horizon."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from scipy.stats import norm, t
|
||||
|
||||
from ..forecast import LEVELS
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Normal:
|
||||
def standard_quantiles(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
return np.broadcast_to(norm.ppf(LEVELS), (len(horizons), len(LEVELS)))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class StudentT:
|
||||
"""Student's t with `df` degrees of freedom, rescaled to unit variance."""
|
||||
|
||||
df: float = 4.0
|
||||
|
||||
def standard_quantiles(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
q = t.ppf(LEVELS, self.df) / np.sqrt(self.df / (self.df - 2))
|
||||
return np.broadcast_to(q, (len(horizons), len(LEVELS)))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Empirical:
|
||||
"""
|
||||
Filtered historical simulation at the horizon level: the shape of past
|
||||
h-day log returns, each divided by the trailing volatility at its start.
|
||||
|
||||
Overlapping h-day returns are far from independent, so a horizon falls back
|
||||
to normal unless the history spans at least `min_windows` of them.
|
||||
"""
|
||||
|
||||
vol_window: int = 365
|
||||
min_windows: int = 3
|
||||
|
||||
def standard_quantiles(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
log_price = np.log(history["close"])
|
||||
sigma = log_price.diff().rolling(self.vol_window).std()
|
||||
out = np.empty((len(horizons), len(LEVELS)))
|
||||
for i, h in enumerate(horizons):
|
||||
z = ((log_price.shift(-h) - log_price) / (sigma * np.sqrt(h))).dropna()
|
||||
if len(z) < self.min_windows * h:
|
||||
out[i] = norm.ppf(LEVELS)
|
||||
else:
|
||||
out[i] = np.quantile((z - z.mean()) / z.std(), LEVELS)
|
||||
return out
|
||||
@@ -0,0 +1,106 @@
|
||||
"""Uncertainty components."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from ..data import log_returns
|
||||
from ..halving import cycle_position
|
||||
from .drift import mean_by_cycle_day
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TrailingVol:
|
||||
"""Standard deviation of daily log returns over the trailing `window` days, scaled by √h."""
|
||||
|
||||
window: int = 365
|
||||
|
||||
def sd(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
return log_returns(history).iloc[-self.window :].std() * np.sqrt(horizons)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class EwmaVol:
|
||||
"""Weighted blend of exponentially weighted standard deviations, scaled by √h."""
|
||||
|
||||
spans: tuple[int, ...] = (30, 90, 180)
|
||||
weights: tuple[float, ...] = (0.2, 0.5, 0.3)
|
||||
|
||||
def sd(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
returns = log_returns(history)
|
||||
sigma = sum(
|
||||
w * returns.ewm(span=s).std().iloc[-1]
|
||||
for s, w in zip(self.spans, self.weights, strict=True)
|
||||
)
|
||||
return sigma * np.sqrt(horizons)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class ReversionVol:
|
||||
"""
|
||||
Volatility starts at its current (short EWMA) level and decays toward a
|
||||
long-run level with a `half_life_days` half-life.
|
||||
|
||||
`long_run="level"` uses the trailing `long_window` standard deviation;
|
||||
`long_run="trend"` extrapolates a log-linear trend in 90-day realised
|
||||
volatility over the same window, since volatility has fallen for years.
|
||||
"""
|
||||
|
||||
now_span: int = 30
|
||||
half_life_days: float = 90.0
|
||||
long_run: str = "level"
|
||||
long_window: int = 1460
|
||||
|
||||
def sd(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
returns = log_returns(history)
|
||||
now = returns.ewm(span=self.now_span).std().iloc[-1]
|
||||
days = np.arange(1, horizons.max() + 1)
|
||||
trailing = returns.iloc[-self.long_window :]
|
||||
if self.long_run == "level":
|
||||
long_run = np.full(len(days), trailing.std())
|
||||
elif self.long_run == "trend":
|
||||
# Non-overlapping 90-day windows, anchored at the origin.
|
||||
realised = trailing.rolling(90).std().iloc[::-90].dropna()
|
||||
age = (realised.index - history.index[-1]).days.to_numpy()
|
||||
slope, intercept = np.polyfit(age, np.log(realised.to_numpy()), 1)
|
||||
long_run = np.exp(intercept + slope * days)
|
||||
else:
|
||||
raise ValueError(f"unknown long_run {self.long_run!r}")
|
||||
phi = 0.5 ** (1 / self.half_life_days)
|
||||
daily_var = long_run**2 + (now**2 - long_run**2) * phi**days
|
||||
return np.sqrt(np.cumsum(daily_var)[horizons - 1])
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class CycleVol:
|
||||
"""
|
||||
Trailing volatility, modulated by how volatile each day of the halving
|
||||
cycle has been relative to the rest of the cycle.
|
||||
|
||||
The per-day ratio is estimated like CycleDrift's drift (kernel-smoothed,
|
||||
recency-weighted), then pulled `shrink` of the way back toward 1. The
|
||||
trailing estimate is first divided by the ratio it was measured under, so
|
||||
the cycle effect isn't counted twice.
|
||||
"""
|
||||
|
||||
window: int = 365
|
||||
bandwidth_days: float = 30.0
|
||||
recency_half_life: float = 1.0
|
||||
shrink: float = 0.5
|
||||
|
||||
def sd(self, history: pd.DataFrame, horizons: np.ndarray) -> np.ndarray:
|
||||
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)
|
||||
squared = returns.to_numpy() ** 2
|
||||
variance = mean_by_cycle_day(squared, day, weight, self.bandwidth_days, prior_days=10.0)
|
||||
overall = (weight * squared).sum() / weight.sum()
|
||||
ratio = 1 + (1 - self.shrink) * (np.sqrt(variance / overall) - 1)
|
||||
|
||||
past_ratio = ratio[day[-self.window :]]
|
||||
base = returns.iloc[-self.window :].std() / np.sqrt(np.mean(past_ratio**2))
|
||||
future = history.index[-1] + pd.to_timedelta(np.arange(1, horizons.max() + 1), unit="D")
|
||||
_, future_day = cycle_position(future)
|
||||
return base * np.sqrt(np.cumsum(ratio[future_day] ** 2)[horizons - 1])
|
||||
+17
-10
@@ -56,8 +56,12 @@ plt.rcParams.update(
|
||||
)
|
||||
|
||||
|
||||
def model_color(name: str) -> str:
|
||||
return SERIES[list(MODELS).index(name) % len(SERIES)]
|
||||
def model_colors(names) -> dict[str, str]:
|
||||
"""Registered models keep their MODELS slot; others follow in order of appearance."""
|
||||
names = list(dict.fromkeys(names))
|
||||
if all(n in MODELS for n in names):
|
||||
return {n: SERIES[list(MODELS).index(n)] for n in names}
|
||||
return {n: SERIES[i % len(SERIES)] for i, n in enumerate(names)}
|
||||
|
||||
|
||||
def price_formatter(x, _=None) -> str:
|
||||
@@ -74,8 +78,9 @@ def fan_chart(history: pd.DataFrame, forecasts: dict[str, Forecast], path: Path)
|
||||
)
|
||||
axes = np.atleast_1d(axes)
|
||||
shown = history[history.index >= history.index[-1] - pd.Timedelta(days=6 * 365)]
|
||||
colors = model_colors(forecasts)
|
||||
for ax, (name, f) in zip(axes, forecasts.items(), strict=True):
|
||||
color = model_color(name)
|
||||
color = colors[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)
|
||||
@@ -108,15 +113,16 @@ def fan_chart(history: pd.DataFrame, forecasts: dict[str, Forecast], path: Path)
|
||||
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."""
|
||||
def skill_chart(summary: pd.DataFrame, path: Path, title: str, baseline: str = BASELINE) -> None:
|
||||
"""CRPS skill vs `baseline`, by horizon, with bootstrap intervals."""
|
||||
horizons = sorted(summary["horizon"].unique())
|
||||
x = np.arange(len(horizons))
|
||||
colors = model_colors(summary["model"])
|
||||
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):
|
||||
ax.axhline(0, color=colors[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)
|
||||
color = colors[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(
|
||||
@@ -129,7 +135,7 @@ def skill_chart(summary: pd.DataFrame, path: Path, title: str) -> None:
|
||||
)
|
||||
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.set_ylabel(f"CRPS skill vs {baseline} (higher is better)")
|
||||
ax.yaxis.set_major_formatter(PercentFormatter(1.0, decimals=0))
|
||||
ax.set_title(title)
|
||||
ax.legend(loc="lower left")
|
||||
@@ -142,6 +148,7 @@ 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))
|
||||
colors = model_colors(summary["model"])
|
||||
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)
|
||||
@@ -153,7 +160,7 @@ def calibration_chart(summary: pd.DataFrame, path: Path, title: str) -> None:
|
||||
ax.plot(
|
||||
x,
|
||||
g[f"cov{c:.0%}"],
|
||||
color=model_color(name),
|
||||
color=colors[name],
|
||||
marker="o",
|
||||
ms=6,
|
||||
mec=SURFACE,
|
||||
|
||||
Reference in New Issue
Block a user