spb/anomaly-atlas Public License

1# =============================================================================2#  Project   : anomaly-atlas3#  File      : src/anomaly_atlas/stats/spa.py4#  Purpose   : White Reality Check & Hansen SPA over a rule-return matrix5#  Author    : Simon-Pierre Boucher6#  Contact   : contact@spboucher.ai7#  Data src  : hfmarketdata.io (sole data source)8#  Created   : 2026-08-129#  Modified  : 2026-08-1210#  Platform  : macOS / Apple Silicon (arm64)11#  License   : All rights reserved (research code)12# =============================================================================13"""Data-snooping corrections over a searched universe of rules.1415Inputs are a (T days × N rules) matrix of rule returns. H0: no rule has16positive expected return — max_k E[f_k] <= 0.1718* White (2000) Reality Check: max-statistic over the centered stationary19  bootstrap (Politis-Romano 1994).20* Hansen (2005) SPA: studentized statistic with the recentering threshold,21  less sensitive to poor/irrelevant rules in the universe.2223Validated on synthetic ground truth (§8.1 gate: pure noise must not24survive; a planted profitable rule must).25"""2627from __future__ import annotations2829import numpy as np303132def stationary_bootstrap_indices(33    n: int, mean_block: float, n_boot: int, seed: int = 4234) -> np.ndarray:35    """(n_boot, n) index matrix from the Politis-Romano stationary bootstrap.3637    Geometric block lengths with mean `mean_block`, circular wrapping —38    resamples preserve short-range dependence in expectation.39    """40    rng = np.random.default_rng(seed)41    p = 1.0 / mean_block42    idx = np.empty((n_boot, n), dtype=np.int64)43    for b in range(n_boot):44        t = 045        while t < n:46            start = rng.integers(0, n)47            length = min(int(rng.geometric(p)), n - t)48            idx[b, t : t + length] = (start + np.arange(length)) % n49            t += length50    return idx515253def reality_check(54    x: np.ndarray, n_boot: int = 500, mean_block: float = 5.0, seed: int = 4255) -> dict:56    """White's Reality Check p-value for max_k mean(x_k) > 0.5758    x: (T, N) rule-return matrix (NaN rows dropped listwise).59    """60    x = np.asarray(x, dtype=float)61    x = x[np.isfinite(x).all(axis=1)]62    t_len, n_rules = x.shape63    if t_len < 30 or n_rules == 0:64        return {"p": float("nan"), "best_rule": None, "v_stat": float("nan")}65    means = x.mean(axis=0)66    v = np.sqrt(t_len) * means.max()67    idx = stationary_bootstrap_indices(t_len, mean_block, n_boot, seed)68    centered = x - means  # White: bootstrap distribution of centered means69    v_boot = np.empty(n_boot)70    for b in range(n_boot):71        v_boot[b] = np.sqrt(t_len) * centered[idx[b]].mean(axis=0).max()72    p = float((np.sum(v_boot >= v) + 1) / (n_boot + 1))73    return {"p": p, "best_rule": int(means.argmax()), "v_stat": float(v),74            "best_mean_daily": float(means.max())}757677def spa_test(78    x: np.ndarray, n_boot: int = 500, mean_block: float = 5.0, seed: int = 4279) -> dict:80    """Hansen's SPA p-value (consistent variant) for max_k mean(x_k) > 0."""81    x = np.asarray(x, dtype=float)82    x = x[np.isfinite(x).all(axis=1)]83    t_len, n_rules = x.shape84    if t_len < 30 or n_rules == 0:85        return {"p": float("nan"), "best_rule": None}86    means = x.mean(axis=0)87    idx = stationary_bootstrap_indices(t_len, mean_block, n_boot, seed)88    boot_means = np.empty((n_boot, n_rules))89    for b in range(n_boot):90        boot_means[b] = x[idx[b]].mean(axis=0)91    omega = np.sqrt(t_len) * boot_means.std(axis=0, ddof=1)92    omega = np.maximum(omega, 1e-12)93    t_stat = float((np.sqrt(t_len) * means / omega).max())94    # Hansen recentering: rules with sufficiently negative means contribute 095    thresh = -omega / np.sqrt(t_len) * np.sqrt(2.0 * np.log(np.log(max(t_len, 3))))96    center = np.where(means >= thresh, means, 0.0)97    t_boot = np.empty(n_boot)98    for b in range(n_boot):99        z = np.sqrt(t_len) * (boot_means[b] - center) / omega100        t_boot[b] = max(z.max(), 0.0)101    p = float((np.sum(t_boot >= max(t_stat, 0.0)) + 1) / (n_boot + 1))102    rule_t = np.sqrt(t_len) * means / omega103    t95 = float(np.percentile(t_boot, 95))104    return {"p": p, "best_rule": int(rule_t.argmax()), "t_stat": t_stat,105            "rule_t": rule_t.tolist(), "t95": t95,106            "n_step1_survivors": int((rule_t >= t95).sum()) if np.isfinite(t95) else 0}107108109def deflated_sharpe(110    sr: float, t_len: int, skew: float, kurt: float,111    n_trials: int, sr_variance: float,112) -> dict:113    """Bailey & López de Prado (2014) Deflated Sharpe Ratio.114115    `sr` is the per-period (e.g. daily) Sharpe of the BEST rule; `sr_variance`116    the variance of Sharpe estimates across the searched universe; `kurt` is117    Pearson kurtosis (normal = 3). Returns the expected max Sharpe under118    pure selection (`sr0`) and DSR = P[true SR > 0 | selection].119    """120    from math import sqrt121122    from scipy.stats import norm123124    if n_trials < 2 or sr_variance <= 0 or t_len < 10:125        return {"sr0": float("nan"), "dsr": float("nan")}126    gamma = 0.5772156649015329127    z1 = norm.ppf(1.0 - 1.0 / n_trials)128    z2 = norm.ppf(1.0 - 1.0 / (n_trials * np.e))129    sr0 = sqrt(sr_variance) * ((1.0 - gamma) * z1 + gamma * z2)130    denom = sqrt(max(1.0 - skew * sr + (kurt - 1.0) / 4.0 * sr**2, 1e-12))131    dsr = float(norm.cdf((sr - sr0) * sqrt(t_len - 1.0) / denom))132    return {"sr0": float(sr0), "dsr": dsr}133