bitcoin-model/btcmodel/models/drift.py

"""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]