spb/anomaly-atlas Public License

1# =============================================================================2#  Project   : anomaly-atlas3#  File      : experiments/micro/expC_reversion_scan/benchmark.py4#  Purpose   : Mean-reversion scan net of the EDGE bounce null (TRAIN only)5#  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 C — reversion scan (protocol pre-specified in hypothesis.md).1415Every output cell is Level 0. The scan's job is triage: which cells show16excess reversion beyond the bounce null after FDR — those go to the expF17correction battery, nothing goes to the atlas from here.18"""1920from __future__ import annotations2122import json23import sys24from collections import defaultdict25from datetime import UTC, datetime26from pathlib import Path2728import numpy as np2930REPO_ROOT = Path(__file__).resolve().parents[3]31sys.path.insert(0, str(REPO_ROOT / "benchmarks"))32sys.path.insert(0, str(REPO_ROOT / "src"))3334from hardware_manifest import collect_manifest  # noqa: E4023536from anomaly_atlas.data.hf_client import HFMarketDataClient  # noqa: E40237from anomaly_atlas.data.universe import (  # noqa: E40238    LIQUID_ETF,39    LIQUID_STOCK,40    TRAIN_SUBPERIODS,41    core_universe,42)43from anomaly_atlas.stats.bootstrap import moving_block_bootstrap, percentile_ci  # noqa: E40244from anomaly_atlas.stats.multiple_testing import (  # noqa: E40245    benjamini_hochberg,46    bootstrap_pvalue,47)48from anomaly_atlas.stats.reversion import ac1, half_life, variance_ratio  # noqa: E40249from anomaly_atlas.validation.artifacts import edge_spread  # noqa: E4025051ADJ = "adj_split"52BOOT_N, BOOT_SEED = 200, 4253TIMEFRAMES = ["1day", "30min", "5min"]54BARS_PER_DAY = {"5min": 78, "30min": 13, "1day": 1}55ONE_MIN_WINDOW = ("2014-01-01", "2016-01-01")  # liquid 12 only (declared)56MIN_RETURNS = {"1day": 350, "30min": 2_000, "5min": 5_000, "1min": 5_000}575859def rth_returns_by_day(bars: list[dict], timeframe: str) -> np.ndarray:60    """Within-day log returns on RTH bars only (no overnight, no LOCF)."""61    days: dict[str, list[float]] = defaultdict(list)62    for b in bars:63        dt = b["datetime"]64        if timeframe == "1day":65            days[dt[:10]].append(np.log(b["close"]))66            continue67        if "09:30" <= dt[11:16] < "16:00":68            days[dt[:10]].append(np.log(b["close"]))69    if timeframe == "1day":70        allp = [v[0] for _, v in sorted(days.items())]71        return np.diff(allp) if len(allp) > 2 else np.array([])72    out = [np.diff(v) for _, v in sorted(days.items()) if len(v) >= 2]73    return np.concatenate(out) if out else np.array([])747576def daily_ohlc(bars: list[dict]) -> tuple[np.ndarray, ...]:77    o = np.array([b["open"] for b in bars])78    h = np.array([b["high"] for b in bars])79    lo = np.array([b["low"] for b in bars])80    c = np.array([b["close"] for b in bars])81    return o, h, lo, c828384def analyze_cell(r: np.ndarray, timeframe: str, spread: float) -> dict | None:85    if len(r) < MIN_RETURNS[timeframe]:86        return None87    var = r.var()88    if var == 0:89        return None90    bounce_ac1 = -(spread**2) / 4.0 / var if np.isfinite(spread) else 0.091    block = max(20, BARS_PER_DAY.get(timeframe, 390) * 5)92    boot = moving_block_bootstrap(r, ac1, block=block, n_boot=BOOT_N, seed=BOOT_SEED)93    excess_boot = boot - bounce_ac194    a = ac1(r)95    lo, hi = percentile_ci(excess_boot)96    boot_vr5 = moving_block_bootstrap(97        r, lambda x: variance_ratio(x, 5), block=block, n_boot=BOOT_N, seed=BOOT_SEED98    )99    boot_vr30 = moving_block_bootstrap(100        r, lambda x: variance_ratio(x, 30), block=block, n_boot=BOOT_N, seed=BOOT_SEED101    )102    return {103        "n": int(len(r)),104        "ac1": round(a, 5),105        "edge_spread": round(spread, 6) if np.isfinite(spread) else None,106        "bounce_ac1": round(bounce_ac1, 5),107        "excess_ac1": round(a - bounce_ac1, 5),108        "excess_ci95": [round(lo, 5), round(hi, 5)],109        "p_excess": bootstrap_pvalue(excess_boot, 0.0),110        "vr5": round(variance_ratio(r, 5), 4),111        "p_vr5": bootstrap_pvalue(boot_vr5, 1.0),112        "vr30": round(variance_ratio(r, 30), 4),113        "p_vr30": bootstrap_pvalue(boot_vr30, 1.0),114    }115116117def main() -> None:118    run_utc = datetime.now(UTC)119    client = HFMarketDataClient()120    universe = core_universe(client.tickers("stock", timeframe="1min", adjustment=ADJ))121122    cells: list[dict] = []123    for asset, ticker, bucket in universe:124        # daily bars per sub-period: returns for 1day cells + EDGE spread input125        for sub, (s, e) in TRAIN_SUBPERIODS.items():126            day_bars = client.get_bars(asset, ticker, "1day", ADJ, s, e)127            if len(day_bars) < 200:128                continue129            o, h, lo, c = daily_ohlc(day_bars)130            spread = edge_spread(o, h, lo, c)131            hl = half_life(np.log(c))132            for tf in TIMEFRAMES:133                bars = day_bars if tf == "1day" else client.get_bars(asset, ticker, tf, ADJ, s, e)134                r = rth_returns_by_day(bars, tf)135                m = analyze_cell(r, tf, spread)136                if m is None:137                    continue138                m |= {"ticker": ticker, "bucket": bucket, "timeframe": tf, "period": sub}139                if tf == "1day":140                    m["half_life_days"] = round(hl, 1) if np.isfinite(hl) else None141                cells.append(m)142            print(f"{ticker} {sub}: done ({len(cells)} cells)")143144    # 1min cells: liquid 12, declared window145    for asset, ticker in [("stock", t) for t in LIQUID_STOCK] + [("etf", t) for t in LIQUID_ETF]:146        s, e = ONE_MIN_WINDOW147        day_bars = client.get_bars(asset, ticker, "1day", ADJ, s, e)148        if len(day_bars) < 200:149            continue150        spread = edge_spread(*daily_ohlc(day_bars))151        bars = client.get_bars(asset, ticker, "1min", ADJ, s, e)152        r = rth_returns_by_day(bars, "1min")153        m = analyze_cell(r, "1min", spread)154        if m is not None:155            m |= {"ticker": ticker, "bucket": "liquid", "timeframe": "1min",156                  "period": "2014-2015"}157            cells.append(m)158        print(f"{ticker} 1min: done")159160    # FDR within each statistic family, all cells jointly161    for key, pkey in [("excess_ac1", "p_excess"), ("vr5", "p_vr5"), ("vr30", "p_vr30")]:162        mask = benjamini_hochberg(np.array([c[pkey] for c in cells]), alpha=0.05)163        for c, rej in zip(cells, mask, strict=True):164            c[f"fdr_{key}"] = bool(rej)165166    def survivors(key: str, sign_key: str, negative: bool) -> list[dict]:167        out = []168        for c in cells:169            if not c[f"fdr_{key}"]:170                continue171            v = c[sign_key] - (1.0 if sign_key.startswith("vr") else 0.0)172            if (v < 0) == negative:173                out.append({k: c[k] for k in ("ticker", "bucket", "timeframe", "period",174                                              sign_key, "edge_spread")})175        return out176177    summary = {178        "cells_total": len(cells),179        "families_tested": 3,180        "fdr_alpha": 0.05,181        "excess_ac1_negative_survivors": survivors("excess_ac1", "excess_ac1", True),182        "excess_ac1_positive_survivors": survivors("excess_ac1", "excess_ac1", False),183        "vr30_below_1_survivors": len(survivors("vr30", "vr30", True)),184        "vr30_above_1_survivors": len(survivors("vr30", "vr30", False)),185        "liquid_2008_2015_intraday_negative": [186            c["ticker"] for c in cells187            if c["bucket"] == "liquid" and c["period"] == "2008-2015"188            and c["timeframe"] in ("5min", "30min")189            and c["fdr_excess_ac1"] and c["excess_ac1"] < 0190        ],191    }192193    results = {194        "experiment": "expC_reversion_scan",195        "run_utc": run_utc.isoformat(),196        "author": "Simon-Pierre Boucher",197        "contact": "contact@spboucher.ai",198        "data_source": "hfmarketdata.io",199        "confidence_level": 0,200        "protocol": {201            "train_subperiods": TRAIN_SUBPERIODS, "adjustment": ADJ,202            "one_min_window": ONE_MIN_WINDOW, "boot": [BOOT_N, BOOT_SEED],203            "min_returns": MIN_RETURNS,204        },205        "summary": summary,206        "cells": cells,207        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},208        "manifest": collect_manifest(),209    }210    out_dir = REPO_ROOT / "results" / "expC_reversion_scan" / run_utc.strftime("%Y%m%dT%H%M%SZ")211    out_dir.mkdir(parents=True)212    (out_dir / "results.json").write_text(json.dumps(results, indent=2) + "\n")213    print(f"\nwrote {out_dir.relative_to(REPO_ROOT)}/results.json")214    print(json.dumps({k: (v if not isinstance(v, list) else len(v))215                      for k, v in summary.items()}, indent=1))216217218if __name__ == "__main__":219    main()220