diff --git a/README.md b/README.md index f123895..6ae433c 100644 --- a/README.md +++ b/README.md @@ -59,20 +59,45 @@ the same way: - `cycle`: the 2024 model's one real idea. Expected return depends on the day of the halving cycle, estimated from past cycles, with recent cycles weighted more. +- `powerlaw`: log price grows linearly in log time since genesis, so growth + keeps slowing. Fitted walk-forward; the exponent has stayed between 5.4 and + 6.0 in every fit since 2014. ### Findings so far (development data) -- Nothing beats the random walk with any confidence at any horizon. At 2-4 - years there are only 3-5 independent outcomes in the whole history. -- `drift_rw` leads at 2-4 years (+13-18% skill, but the intervals span zero). -- `cycle` loses to both at every horizon, and so does every setting tried - (recency half-life 0.25-2 cycles, smoothing bandwidth 15-60 days). The cycle - *shape* costs accuracy. The level is the problem: each cycle has grown less - than the last (log return 4.0, 2.6, 2.0, i.e. roughly ×55, ×13, ×7), so any - average of past cycles overshoots. +- `cycle` loses to the random walk at every horizon, and so does every setting + tried (recency half-life 0.25-2 cycles, smoothing bandwidth 15-60 days). + The level is the problem: each cycle has grown less than the last (log + return 4.0, 2.6, 2.0, i.e. roughly ×55, ×13, ×7), so any average of past + cycles overshoots. +- `powerlaw` models exactly that, and it is the first model to beat the random + walk with some confidence: +53% and +63% skill at 3 and 4 years, with + unbiased outcomes (mean PIT 0.51). Only 3-4 independent windows back that + up, the functional form is famous *because* it fits Bitcoin's history, and + the holdout hasn't been run yet. +- Its intervals are too wide at long horizons (the 80% interval held every + 3-year outcome), because it treats deviations from the trend as permanent. -That points at diminishing returns as the structure worth modelling, e.g. a -power-law trend, which is next. +### A/B tests of the 2024 ideas + +`just ab` runs the experiments in `btcmodel/experiments.py`: ideas salvaged +from the old branches (catalogued in [docs/2024-ideas.md](docs/2024-ideas.md)), +each a control plus variants that change one component. Hypotheses and the +verdict rule were written down before anything ran. With ~100 comparisons, +expect a few flukes either way. + +| Experiment | Idea | Verdict | +|---|---|---| +| shrink-cycle | scale the cycle drift by 0.25/0.5/0.75 | better, all three: it fixes the level crudely | +| diminishing-returns | power-law trend, optionally reverting to it, or the cycle shape rescaled to it | better (both power-law variants); cycle shape on the power law inconclusive | +| vol-window | EWMA blends, shorter or longer windows | worse: the plain 365-day window wins | +| vol-reversion | volatility reverting to a long-run level or falling trend | worse | +| cycle-vol | volatility by cycle position | inconclusive (no effect) | +| tails | Student-t, or empirical horizon-level shape | Student-t worse; empirical +5% at 1 month only | +| cycle-phase | align cycles by fraction elapsed, not days | inconclusive | + +Next: a power law whose deviations revert with bounded variance, and a +registered test of whether the cycle's timing adds anything on top of it. ## Usage @@ -83,6 +108,7 @@ nix develop just update # fetch new daily prices from Coinbase just backtest # score models on development data -> output/backtest/ just forecast # forecast from the latest price -> output/forecast/ +just ab # run A/B experiments -> output/ab/ just test just holdout # score on held-out outcomes; sparingly ``` @@ -98,7 +124,8 @@ btcmodel/ halving.py halving calendar, position in cycle forecast.py Forecast (quantiles of log price), CRPS, PIT evaluate.py walk-forward backtest and summary - models/ one file per model family; register new ones in __init__.py + experiments.py A/B tests: hypothesis, control, variants, verdict rule + models/ drift, volatility and shape components; register models in __init__.py plots.py fan chart, skill and calibration charts ``` diff --git a/btcmodel/__main__.py b/btcmodel/__main__.py index 54a4f67..069f3ae 100644 --- a/btcmodel/__main__.py +++ b/btcmodel/__main__.py @@ -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) diff --git a/btcmodel/evaluate.py b/btcmodel/evaluate.py index c2efc61..8d02a0b 100644 --- a/btcmodel/evaluate.py +++ b/btcmodel/evaluate.py @@ -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) ], diff --git a/btcmodel/experiments.py b/btcmodel/experiments.py new file mode 100644 index 0000000..71c8db9 --- /dev/null +++ b/btcmodel/experiments.py @@ -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()),), + ), + ) +} diff --git a/btcmodel/halving.py b/btcmodel/halving.py index f0f58b9..118c530 100644 --- a/btcmodel/halving.py +++ b/btcmodel/halving.py @@ -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) diff --git a/btcmodel/models/__init__.py b/btcmodel/models/__init__.py index 54657c2..3c9d6b8 100644 --- a/btcmodel/models/__init__.py +++ b/btcmodel/models/__init__.py @@ -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()), + ) +} diff --git a/btcmodel/models/base.py b/btcmodel/models/base.py new file mode 100644 index 0000000..12224bf --- /dev/null +++ b/btcmodel/models/base.py @@ -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) diff --git a/btcmodel/models/baselines.py b/btcmodel/models/baselines.py deleted file mode 100644 index 3ec3db1..0000000 --- a/btcmodel/models/baselines.py +++ /dev/null @@ -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)) diff --git a/btcmodel/models/cycle.py b/btcmodel/models/cycle.py deleted file mode 100644 index e225dc2..0000000 --- a/btcmodel/models/cycle.py +++ /dev/null @@ -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)) diff --git a/btcmodel/models/drift.py b/btcmodel/models/drift.py new file mode 100644 index 0000000..ce1c00d --- /dev/null +++ b/btcmodel/models/drift.py @@ -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] diff --git a/btcmodel/models/shape.py b/btcmodel/models/shape.py new file mode 100644 index 0000000..b0fe317 --- /dev/null +++ b/btcmodel/models/shape.py @@ -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 diff --git a/btcmodel/models/volatility.py b/btcmodel/models/volatility.py new file mode 100644 index 0000000..5d4f148 --- /dev/null +++ b/btcmodel/models/volatility.py @@ -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]) diff --git a/btcmodel/plots.py b/btcmodel/plots.py index 95c49b9..23437f6 100644 --- a/btcmodel/plots.py +++ b/btcmodel/plots.py @@ -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, diff --git a/docs/2024-ideas.md b/docs/2024-ideas.md new file mode 100644 index 0000000..dcdb4c3 --- /dev/null +++ b/docs/2024-ideas.md @@ -0,0 +1,65 @@ +# Ideas in the 2024 model's history + +A catalogue of the modeling ideas in the old `model.py` and its branches +(September 2026). Labels are referenced from `btcmodel/experiments.py`. +"Fudge" means the constant or rule existed to hit backtest coverage/MAPE +targets rather than to express a hypothesis. + +## Drift + +- **D1. Mean return by halving-cycle position.** The core idea throughout + (initial commit, `Switch to simpler log-based projection`, `Implement trend + smoothing`, `tuning-b`). Smoothing went Savitzky-Golay → 60-day centred mean + → Gaussian kernel with recency weights. The old code stretched every cycle + to 1460 days (phase = fraction of the halving interval). Bugs: simple + returns compounded as log returns (+~25%/yr bias) in the initial commit; + genesis date off by a year; `tuning-3` regressed log price on row index + across cycles. The "blend ratio of 0 is best?" commit was blending the + simple- and log-return versions of the same quantity. +- **D2. Drift damping.** Constants everywhere (×0.6–0.9, asymmetric, era- and + cycle-position-keyed, a 3%/day cap, the fundamentals ×0.65–0.75). Fudge as + implemented; the hypothesis underneath (the cycle drift is overfit) is real. +- **D3. Diminishing returns.** Initial commit decayed drift by 0.9^(t/365); + `old/backtests-trend-enhancement-1` regressed per-cycle returns on cycle + number, but added that on top of the cycle drift (double counting) and + clipped paths. +- **D4. "Skew".** `loc += sign(μ)·0.087σ`, tuned so 68% coverage hit 68.1%. + A fudge; its honest cousin is momentum. +- **D5. Stock-to-flow.** A 30% blend of the daily % change in S2F into the + drift. Broken: mismatched units, and S2F "halvings" on the wrong dates + created a large fake post-halving drift. + +## Volatility + +- **V1. EWMA blend.** 30/90/180-day spans, weights .5/.3/.2 then .2/.5/.3, + times a 1.2 fudge. +- **V2. Adaptive windows / regime weights.** Short/long vol ratio rescales + window lengths and blend weights; over-parameterised, and after the + fundamentals rewrite the adaptive windows were computed but unused. +- **V3. Cycle-position volatility.** `old/backtest-vol-cycle`. +- **V4. Market maturity.** Volume growth, inverse vol, autocorrelation and a + post-futures dummy, min-max normalised in-sample. Removed as "complexity + without clear benefit". +- **V5. "Fundamentals".** Supply growth, volume/supply and "depth" combined + and clipped to [0.65, 0.75]: effectively a constant, grid-searched to match + the era constants it replaced. +- **V6. Market conditions.** Vol ratio, MA50−MA200 trend strength (≈30× too + large from a units bug) and drawdown widening uncertainty. + +## Calibration (all fudges) + +- **C1. Era scaling keyed to the backtest's training start date**, with stacked + factors (net ≈0.45× vol) and hand-picked era boundaries. +- **C2. Horizon uncertainty multipliers**, caps and floors on daily vol. +- **C3. Changing the quantile levels**: a "95%" band read off the ~81% + quantiles past a year. + +## Shape + +- **S1. Skewed innovations**, with skew set per era (`tuning-1`). + +## Also + +Plumbing (backtest framework, multiprocessing, output handling, cycle and +CDPR plots). The old NOTES.md described "machine learning" and "macro +indicator" experiments; no such code ever existed. diff --git a/justfile b/justfile index a1b7a81..9910ef3 100644 --- a/justfile +++ b/justfile @@ -32,3 +32,7 @@ lint: clean: rm -rf output + +# Run A/B experiments on development data (default: all) +ab *names: + python -m btcmodel ab {{ names }} diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_components.py b/tests/test_components.py new file mode 100644 index 0000000..af439db --- /dev/null +++ b/tests/test_components.py @@ -0,0 +1,66 @@ +import numpy as np +import pandas as pd +import pytest + +from btcmodel.experiments import EXPERIMENTS +from btcmodel.halving import GENESIS +from btcmodel.models.drift import CycleDrift, PowerLawDrift, ShrunkDrift +from btcmodel.models.shape import StudentT +from btcmodel.models.volatility import ReversionVol, TrailingVol + +from .test_models import synthetic_prices + +ALL_VARIANTS = {m.name: m for e in EXPERIMENTS.values() for m in e.models} + + +@pytest.mark.parametrize("name", list(ALL_VARIANTS)) +def test_experiment_models_produce_valid_forecasts(name): + prices = synthetic_prices(lambda day: 0.001 + 0 * day, noise=0.03) + horizons = np.array([1, 30, 365, 1460]) + f = ALL_VARIANTS[name].forecast(prices, horizons) + assert np.isfinite(f.log_quantiles).all() + assert np.all(np.diff(f.log_quantiles, axis=1) >= 0) + lo, hi = f.interval(0.8) + assert np.all(np.diff(hi - lo) > 0) + + +def test_shrunk_drift_scales_linearly(): + prices = synthetic_prices(lambda day: np.where(day < 700, 0.002, -0.001)) + horizons = np.array([100, 1000]) + full = CycleDrift().expected_log_return(prices, horizons) + np.testing.assert_allclose( + ShrunkDrift(CycleDrift(), 0.5).expected_log_return(prices, horizons), full / 2 + ) + np.testing.assert_allclose( + ShrunkDrift(CycleDrift(), 0.0).expected_log_return(prices, horizons), 0 + ) + + +def test_power_law_recovers_its_exponent(): + dates = pd.date_range("2011-01-01", "2024-11-26", freq="D", name="date") + t = (dates - GENESIS).days.to_numpy() + prices = pd.DataFrame({"close": 1e-17 * t**5.8}, index=dates) + _, slope, _ = PowerLawDrift().fit(prices) + assert slope == pytest.approx(5.8, rel=1e-6) + expected = 5.8 * np.log((t[-1] + 365) / t[-1]) + assert PowerLawDrift().expected_log_return(prices, np.array([365]))[0] == pytest.approx( + expected + ) + + +def test_reversion_vol_matches_trailing_when_already_at_long_run(): + rng = np.random.default_rng(0) + dates = pd.date_range("2011-01-01", "2024-11-26", freq="D", name="date") + prices = pd.DataFrame( + {"close": 100 * np.exp(np.cumsum(0.03 * rng.standard_normal(len(dates))))}, index=dates + ) + horizons = np.array([30, 365, 1460]) + reverting = ReversionVol(now_span=1460, long_window=len(dates)).sd(prices, horizons) + flat = TrailingVol(window=len(dates)).sd(prices, horizons) + np.testing.assert_allclose(reverting, flat, rtol=0.05) + + +def test_student_t_shape_has_unit_variance_and_fatter_tails(): + q = StudentT(4).standard_quantiles(None, np.array([1]))[0] + assert q[-1] > 2.576 # beyond the normal 99.5% quantile + assert np.interp(0.8413, np.linspace(0.005, 0.995, 100), q) < 1.0 # thinner shoulders diff --git a/tests/test_models.py b/tests/test_models.py index 5f8b378..c049f4a 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -5,7 +5,7 @@ import pytest from btcmodel import data from btcmodel.evaluate import backtest from btcmodel.halving import HALVINGS, cycle_position -from btcmodel.models import MODELS, CycleModel +from btcmodel.models import MODELS, CycleDrift def synthetic_prices(daily_return, start="2011-01-01", end="2024-11-26", noise=0.0, seed=0): @@ -31,7 +31,7 @@ def test_cycle_model_recovers_a_cycle_shaped_drift(): return np.where(day < 700, 0.002, -0.001) prices = synthetic_prices(shape) - drift = CycleModel(prior_days=0).drift_by_cycle_day(prices) + drift = CycleDrift(prior_days=0).by_cycle_day(prices) assert drift[300] == pytest.approx(0.002, abs=2e-4) assert drift[1100] == pytest.approx(-0.001, abs=2e-4) @@ -42,8 +42,8 @@ def test_cycle_model_weights_recent_cycles_more(): cycle, _ = cycle_position(prices.index) r = 0.004 / 2.0**cycle prices["close"] = 100 * np.exp(np.cumsum(r)) - drift = CycleModel(recency_half_life=0.25).drift_by_cycle_day(prices) - equal = CycleModel(recency_half_life=1e9).drift_by_cycle_day(prices) + drift = CycleDrift(recency_half_life=0.25).by_cycle_day(prices) + equal = CycleDrift(recency_half_life=1e9).by_cycle_day(prices) latest_complete = 0.004 / 2.0**3 assert abs(drift[900] - latest_complete) < abs(equal[900] - latest_complete)