Chapter 51: Time Series Cross-Validation

import numpy as np
from sklearn.linear_model import Ridge

from kaggle_companion.activities._common import clean

SEED = 51
DRAWS = 160                      # independent series; every estimate is a mean over draws (40 was too few to rank the recipes)
DAYS, ROWS, FEATURES = 420, 10, 5
BLOCK = 30                       # each assessment block covers the 30 days after its forecast origin
ORIGINS = [100, 130, 160, 190, 220]
DRIFT = 0.04                     # daily pull of the coefficients toward a fresh random draw (about 25 days of memory)
SHORT_WINDOW = 30                # the "recent history only" recipe trains on the last 30 usable days


def generate(rng, drift):
    """Ten rows a day. The target is linear in five features whose coefficients drift slowly (a stationary process), plus noise of sd 1."""
    beta = np.zeros((DAYS, FEATURES))
    beta[0] = rng.normal(size=FEATURES)
    for t in range(1, DAYS):
        beta[t] = (1 - drift) * beta[t - 1] + np.sqrt(1 - (1 - drift) ** 2) * rng.normal(size=FEATURES)
    X = rng.normal(size=(DAYS, ROWS, FEATURES))
    return X, (beta[:, None, :] * X).sum(2) + rng.normal(size=(DAYS, ROWS))


def backtest_mae(X, y, delay, window, respect_availability):
    """Mean absolute error over the assessment blocks. A row observed on day d has its label available on day d + delay.
    Honest rule: train on rows with d + delay < origin. Naive rule: train on every row observed before the origin,
    including rows whose labels would still be on the way. `window` keeps only the most recent usable days (None keeps all)."""
    errors = []
    for origin in ORIGINS:
        last = origin - delay if respect_availability else origin
        first = 0 if window is None else max(0, last - window)
        model = Ridge(alpha=1.0).fit(X[first:last].reshape(-1, FEATURES), y[first:last].reshape(-1))
        block = slice(origin, origin + BLOCK)
        errors.append(np.mean(np.abs(model.predict(X[block].reshape(-1, FEATURES)) - y[block].reshape(-1))))
    return float(np.mean(errors))


def run(delay):
    rng = np.random.default_rng(SEED)
    recipes = {"short": SHORT_WINDOW, "expanding": None}
    rows = {f"{rule}_{name}": [] for rule in ("naive", "honest") for name in recipes}
    stable = []
    for _ in range(DRAWS):
        X, y = generate(rng, DRIFT)
        for name, window in recipes.items():
            for rule in ("naive", "honest"):
                rows[f"{rule}_{name}"].append(backtest_mae(X, y, delay, window, rule == "honest"))
        X, y = generate(rng, 0.0)                              # control world: coefficients never change
        stable.append(backtest_mae(X, y, delay, None, True) - backtest_mae(X, y, delay, None, False))
    mean = {k: float(np.mean(v)) for k, v in rows.items()}
    # Optimism of the naive backtest: how much lower its error is than the honest one, per recipe (paired by series).
    optimism = {name: np.array(rows[f"honest_{name}"]) - np.array(rows[f"naive_{name}"]) for name in recipes}
    # Difference between recipes under each rule; negative means the short-window recipe has the lower error.
    naive_diff = np.array(rows["naive_short"]) - np.array(rows["naive_expanding"])
    honest_diff = np.array(rows["honest_short"]) - np.array(rows["honest_expanding"])
    se = lambda a: np.std(a, ddof=1) / np.sqrt(DRAWS)
    return clean({
        "delay": delay, **mean,
        **{f"optimism_{name}": optimism[name].mean() for name in recipes},
        **{f"optimism_{name}_se": se(optimism[name]) for name in recipes},
        "naive_diff": naive_diff.mean(), "naive_diff_se": se(naive_diff),
        "honest_diff": honest_diff.mean(), "honest_diff_se": se(honest_diff),
        "stable_gap": np.mean(stable), "stable_gap_se": se(stable),
        "draws": DRAWS, "origins": len(ORIGINS), "noise_floor": 0.7979,
    })

import pandas as pd
from kaggle_companion.activities.ch51 import explain

VALUES = [0, 10, 30, 60]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch51 import draw, NCOLS, HEIGHT

DEFAULT = 30
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()