spb/anomaly-atlas Public License

1# =============================================================================2#  Project   : anomaly-atlas3#  File      : experiments/micro/expF_multiple_testing/benchmark.py4#  Purpose   : Survival battery: naive -> FDR -> RC/SPA -> DSR over C/D/E rules5#  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"""Experiment F — the survival curve (protocol pre-specified in hypothesis.md).1415Rules are built mechanically from EVERYTHING the C/D/E scans searched (both16signs), evaluated on the same TRAIN data (in-sample by design — OOS is expH),17and pushed through naive-t -> BH-FDR -> White RC / Hansen SPA -> DSR.18"""1920from __future__ import annotations2122import json23import sys24from collections import defaultdict25from datetime import UTC, datetime26from pathlib import Path2728import numpy as np29from scipy.stats import norm3031REPO_ROOT = Path(__file__).resolve().parents[3]32sys.path.insert(0, str(REPO_ROOT / "benchmarks"))33sys.path.insert(0, str(REPO_ROOT / "src"))3435from hardware_manifest import collect_manifest  # noqa: E4023637from anomaly_atlas.data.cleaning import RTH_SLOTS, rth_day_grids  # noqa: E40238from anomaly_atlas.data.hf_client import HFMarketDataClient  # noqa: E40239from anomaly_atlas.data.universe import TRAIN, TRAIN_SUBPERIODS, core_universe  # noqa: E40240from anomaly_atlas.stats.multiple_testing import benjamini_hochberg  # noqa: E40241from anomaly_atlas.stats.spa import deflated_sharpe, reality_check, spa_test  # noqa: E4024243ADJ = "adj_split"44N_BOOT, MEAN_BLOCK, SEED = 500, 5.0, 4245LIQUID = ["AAPL", "MSFT", "NVDA", "AMZN", "GOOGL", "META", "TSLA", "JPM", "XOM", "UNH", "QQQ"]46SECTORS = ["XLF", "XLE", "XLK", "XLV", "XLI", "XLY", "XLP", "XLU", "XLB"]47D_WINDOWS = {"2006-2007": ("2006-01-01", "2008-01-01"),48             "2014-2015": ("2014-01-01", "2016-01-01")}49ONE_MIN_WINDOW = ("2014-01-01", "2016-01-01")505152def contrarian_daily(bars: list[dict], timeframe: str) -> dict[str, float]:53    """day -> contrarian rule return at the cell's timeframe (trade-time RTH)."""54    by_day: dict[str, list[float]] = defaultdict(list)55    for b in bars:56        dt = b["datetime"]57        if timeframe == "1day":58            by_day[dt[:10]].append(np.log(b["close"]))59        elif "09:30" <= dt[11:16] < "16:00":60            by_day[dt[:10]].append(np.log(b["close"]))61    days = sorted(by_day)62    out: dict[str, float] = {}63    if timeframe == "1day":64        closes = np.array([by_day[d][0] for d in days])65        r = np.diff(closes)66        for i in range(1, len(r)):67            out[days[i + 1]] = float(-np.sign(r[i - 1]) * r[i])68        return out69    for d in days:70        p = np.array(by_day[d])71        if len(p) < 3:72            continue73        r = np.diff(p)74        out[d] = float(np.sum(-np.sign(r[:-1]) * r[1:]))75    return out767778def leadlag_daily(gx: dict, gy: dict, day_list: list[str]) -> dict[str, float]:79    """day -> sign(leader_{t-1}) * follower_t summed over both-fresh minutes."""80    out: dict[str, float] = {}81    for day in day_list:82        px, py = gx.get(day), gy.get(day)83        if px is None or py is None:84            continue85        ox, oy = np.isfinite(px), np.isfinite(py)86        fx = np.zeros(RTH_SLOTS - 1, dtype=bool)87        fy = np.zeros(RTH_SLOTS - 1, dtype=bool)88        rx = np.zeros(RTH_SLOTS - 1)89        ry = np.zeros(RTH_SLOTS - 1)90        fpx, fpy = px.copy(), py.copy()91        for t in range(1, RTH_SLOTS):92            if not ox[t]:93                fpx[t] = fpx[t - 1]94            if not oy[t]:95                fpy[t] = fpy[t - 1]96        rx[:] = np.diff(fpx)97        ry[:] = np.diff(fpy)98        fx[:] = ox[1:] & ox[:-1]99        fy[:] = oy[1:] & oy[:-1]100        keep = fx[:-1] & fy[1:] & np.isfinite(rx[:-1]) & np.isfinite(ry[1:])101        if keep.sum() < 30:102            continue103        out[day] = float(np.sum(np.sign(rx[:-1][keep]) * ry[1:][keep]))104    return out105106107def build_matrix(rules: dict[str, dict[str, float]]) -> tuple[np.ndarray, list[str], list[str]]:108    """(days x 2N signed rules) matrix; missing day = 0 (idle), as declared."""109    day_set = sorted({d for r in rules.values() for d in r})110    names, cols = [], []111    for name, series in rules.items():112        v = np.array([series.get(d, 0.0) for d in day_set])113        names += [f"+{name}", f"-{name}"]114        cols += [v, -v]115    return np.column_stack(cols), names, day_set116117118def battery(x: np.ndarray, names: list[str]) -> dict:119    t_len = x.shape[0]120    mu = x.mean(axis=0)121    sd = x.std(axis=0, ddof=1)122    sd = np.maximum(sd, 1e-12)123    t = np.sqrt(t_len) * mu / sd124    p_two = 2 * (1 - norm.cdf(np.abs(t)))125    rc = reality_check(x, n_boot=N_BOOT, mean_block=MEAN_BLOCK, seed=SEED)126    sp = spa_test(x, n_boot=N_BOOT, mean_block=MEAN_BLOCK, seed=SEED)127    srs = mu / sd128    best = int(np.argmax(srs))129    r = x[:, best]130    dsr = deflated_sharpe(131        sr=float(srs[best]), t_len=t_len,132        skew=float(((r - r.mean()) ** 3).mean() / r.std() ** 3),133        kurt=float(((r - r.mean()) ** 4).mean() / r.std() ** 4),134        n_trials=x.shape[1], sr_variance=float(srs.var(ddof=1)),135    )136    return {137        "days": int(t_len), "rules": len(names),138        "naive_t196": int((np.abs(t) > 1.96).sum()),139        "p_two_sided": p_two, "t": t,140        "rc_p": rc["p"], "spa_p": sp["p"],141        "spa_step1_survivors": [names[i] for i, tv in enumerate(sp["rule_t"])142                                if tv >= sp["t95"]],143        "best_rule": names[best], "best_daily_sharpe": round(float(srs[best]), 4),144        "best_annualized_sharpe": round(float(srs[best] * np.sqrt(252)), 3),145        "dsr_sr0": round(dsr["sr0"], 4), "dsr": round(dsr["dsr"], 4),146    }147148149def main() -> None:150    run_utc = datetime.now(UTC)151    client = HFMarketDataClient()152    universe = core_universe(client.tickers("stock", timeframe="1min", adjustment=ADJ))153154    blocks: dict[str, dict[str, dict[str, float]]] = defaultdict(dict)155156    # ---- R-family (expC universe)157    for asset, ticker, _ in universe:158        for sub, (s, e) in TRAIN_SUBPERIODS.items():159            day_bars = client.get_bars(asset, ticker, "1day", ADJ, s, e)160            if len(day_bars) < 200:161                continue162            for tf in ("1day", "30min", "5min"):163                bars = day_bars if tf == "1day" else client.get_bars(asset, ticker, tf, ADJ, s, e)164                series = contrarian_daily(bars, tf)165                if len(series) >= 150:166                    blocks[f"expC {sub}"][f"R:{ticker}:{tf}"] = series167        print(f"R {ticker} done")168    for asset, ticker in [("stock", t) for t in LIQUID if t != "QQQ"] + \169                         [("etf", t) for t in ("SPY", "QQQ")]:170        bars = client.get_bars(asset, ticker, "1min", ADJ, *ONE_MIN_WINDOW)171        series = contrarian_daily(bars, "1min")172        if len(series) >= 150:173            blocks["expC 1min 2014-2015"][f"R:{ticker}:1min"] = series174175    # ---- L-family (expD universe)176    random10 = sorted(np.random.default_rng(42).choice(177        sorted(client.tickers("stock", timeframe="1min", adjustment=ADJ)), 30,178        replace=False))[:10]179    for window, (s, e) in D_WINDOWS.items():180        spec = ({"SPY": ("etf", "SPY", ADJ)}181                | {t: ("stock", t, ADJ) for t in LIQUID if t != "QQQ"}182                | {"QQQ": ("etf", "QQQ", ADJ)}183                | {t: ("etf", t, ADJ) for t in SECTORS}184                | {t: ("stock", t, ADJ) for t in random10})185        if window == "2014-2015":186            spec |= {"SPX": ("index", "SPX", None),187                     "ES": ("futures", "ES", "contin_adj_ratio")}188        grids = {}189        for name, (asset, ticker, adj) in spec.items():190            g = rth_day_grids(client.get_bars(asset, ticker, "1min", adj, s, e))191            if len(g) >= 200:192                grids[name] = g193        day_list = sorted(grids["SPY"].keys())194        for name, g in grids.items():195            if name == "SPY":196                continue197            if name in ("SPX", "ES"):198                series = leadlag_daily(g, grids["SPY"], day_list)  # x leads SPY199                key = f"L:{name}->SPY"200            else:201                series = leadlag_daily(grids["SPY"], g, day_list)202                key = f"L:SPY->{name}"203            if len(series) >= 150:204                blocks[f"expD {window}"][key] = series205        print(f"L {window} done ({len(blocks[f'expD {window}'])} pairs)")206207    # ---- C-family (expE classes, drift-adjusted)208    bars = client.get_bars("etf", "SPY", "1day", "adj_splitdiv", TRAIN[0], TRAIN[1])209    dates = [b["datetime"][:10] for b in bars][1:]210    rets = np.diff(np.log([b["close"] for b in bars]))211    mu = rets.mean()212    sys.path.insert(0, str(REPO_ROOT / "experiments" / "micro" / "expE_calendar_scan"))213    from benchmark import class_masks  # noqa: E402214215    masks = class_masks(dates)216    for cname, m in masks.items():217        blocks["expE train"][f"C:{cname}"] = {218            d: float(rets[i] - mu) for i, d in enumerate(dates) if m[i]219        }220221    # ---- battery per block + global funnel222    per_block: dict[str, dict] = {}223    all_p, all_index = [], []224    for bname, rules in blocks.items():225        x, names, _ = build_matrix(rules)226        res = battery(x, names)227        all_p.extend(res.pop("p_two_sided").tolist())228        res.pop("t")229        all_index.extend([(bname, n) for n in names])230        per_block[bname] = res231        print(f"{bname}: rules={res['rules']} naive={res['naive_t196']} "232              f"rc_p={res['rc_p']} spa_p={res['spa_p']} dsr={res['dsr']}")233234    fdr_mask = benjamini_hochberg(np.array(all_p), alpha=0.05)235    total_rules = len(all_p)236    funnel = {237        "universe_rules": total_rules,238        "naive_t196": int(sum(per_block[b]["naive_t196"] for b in per_block)),239        "fdr_survivors": int(fdr_mask.sum()),240        "spa_step1_survivors": sorted({n for b in per_block241                                       for n in per_block[b]["spa_step1_survivors"]}),242        "blocks_spa_significant": [b for b in per_block if per_block[b]["spa_p"] < 0.05],243        "dsr_by_block": {b: per_block[b]["dsr"] for b in per_block},244    }245    funnel["survival_rate"] = {246        "naive": round(funnel["naive_t196"] / total_rules, 4),247        "fdr": round(funnel["fdr_survivors"] / total_rules, 4),248        "spa_step1": round(len(funnel["spa_step1_survivors"]) / total_rules, 4),249    }250251    results = {252        "experiment": "expF_multiple_testing",253        "run_utc": run_utc.isoformat(),254        "author": "Simon-Pierre Boucher",255        "contact": "contact@spboucher.ai",256        "data_source": "hfmarketdata.io",257        "confidence_level": 0,258        "protocol": {"n_boot": N_BOOT, "mean_block": MEAN_BLOCK, "seed": SEED,259                     "note": "in-sample search survival on TRAIN; OOS = expH"},260        "per_block": per_block,261        "funnel": funnel,262        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},263        "manifest": collect_manifest(),264    }265    out_dir = REPO_ROOT / "results" / "expF_multiple_testing" / run_utc.strftime("%Y%m%dT%H%M%SZ")266    out_dir.mkdir(parents=True)267    (out_dir / "results.json").write_text(json.dumps(results, indent=2) + "\n")268    print(f"\nwrote {out_dir.relative_to(REPO_ROOT)}/results.json")269    print(json.dumps(funnel, indent=1))270271272if __name__ == "__main__":273    main()274