spb/anomaly-atlas Public License

1# =============================================================================2#  Project   : anomaly-atlas3#  File      : experiments/micro/expD_leadlag_scan/benchmark.py4#  Purpose   : Lead-lag scan — raw LOCF vs both-fresh, artifact share, FDR5#  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 D — lead-lag scan (protocol pre-specified in hypothesis.md).1415The design measures T3 directly: every pair reports its raw-LOCF lead, its16both-fresh lead, and the difference (artifact share). Level 0 throughout.17"""1819from __future__ import annotations2021import json22import sys23from datetime import UTC, datetime24from pathlib import Path2526import numpy as np2728REPO_ROOT = Path(__file__).resolve().parents[3]29sys.path.insert(0, str(REPO_ROOT / "benchmarks"))30sys.path.insert(0, str(REPO_ROOT / "src"))3132from hardware_manifest import collect_manifest  # noqa: E4023334from anomaly_atlas.data.cleaning import (  # noqa: E40235    RTH_SLOTS,36    both_fresh,37    nan_xcorr,38    rth_day_grids,39)40from anomaly_atlas.data.hf_client import HFMarketDataClient  # noqa: E40241from anomaly_atlas.data.universe import random_stock_universe  # noqa: E40242from anomaly_atlas.stats.multiple_testing import (  # noqa: E40243    benjamini_hochberg,44    bootstrap_pvalue,45)4647ADJ = "adj_split"48LIQUID = ["AAPL", "MSFT", "NVDA", "AMZN", "GOOGL", "META", "TSLA", "JPM", "XOM", "UNH", "QQQ"]49SECTORS = ["XLF", "XLE", "XLK", "XLV", "XLI", "XLY", "XLP", "XLU", "XLB"]50WINDOWS = {"2006-2007": ("2006-01-01", "2008-01-01"),51           "2014-2015": ("2014-01-01", "2016-01-01")}52MAX_LAG = 353BOOT_N, BOOT_SEED = 200, 4254MIN_DAYS, MIN_FRESH = 200, 10_000555657def day_matrix(days: dict[str, np.ndarray], day_list: list[str]) -> tuple[np.ndarray, np.ndarray]:58    """(n_days, 389) LOCF returns + interval-fresh masks, day-aligned."""59    rets = np.full((len(day_list), RTH_SLOTS - 1), np.nan)60    fresh = np.zeros((len(day_list), RTH_SLOTS - 1), dtype=bool)61    for i, day in enumerate(day_list):62        p = days.get(day)63        if p is None:64            continue65        observed = np.isfinite(p)66        filled = p.copy()67        for t in range(1, RTH_SLOTS):68            if not observed[t]:69                filled[t] = filled[t - 1]70        rets[i] = np.diff(filled)71        fresh[i] = observed[1:] & observed[:-1]72    return rets, fresh737475def day_block_boot_xcorr(76    rx: np.ndarray, fx: np.ndarray, ry: np.ndarray, fy: np.ndarray, lag: int77) -> np.ndarray:78    """Bootstrap distribution of the both-fresh xcorr at `lag`, resampling days."""79    rng = np.random.default_rng(BOOT_SEED)80    n_days = rx.shape[0]81    out = np.empty(BOOT_N)82    for b in range(BOOT_N):83        idx = rng.integers(0, n_days, n_days)84        bx, by = both_fresh(rx[idx].ravel(), fx[idx].ravel(), ry[idx].ravel(), fy[idx].ravel())85        out[b] = nan_xcorr(bx, by, abs(lag)).get(lag, np.nan)86    return out878889def analyze_pair(rx, fx, ry, fy) -> dict | None:90    bx, by = both_fresh(rx.ravel(), fx.ravel(), ry.ravel(), fy.ravel())91    n_fresh = int(np.isfinite(bx).sum())92    if n_fresh < MIN_FRESH:93        return None94    raw = nan_xcorr(rx.ravel(), ry.ravel(), MAX_LAG)95    fresh_xc = nan_xcorr(bx, by, MAX_LAG)96    boot_p1 = day_block_boot_xcorr(rx, fx, ry, fy, 1)97    boot_m1 = day_block_boot_xcorr(rx, fx, ry, fy, -1)98    return {99        "n_fresh_pairs": n_fresh,100        "fresh_fraction": round(n_fresh / rx.size, 4),101        "raw_xcorr": {str(k): round(v, 5) if np.isfinite(v) else None for k, v in raw.items()},102        "fresh_xcorr": {str(k): round(v, 5) if np.isfinite(v) else None103                        for k, v in fresh_xc.items()},104        "artifact_share_+1": round(raw[1] - fresh_xc[1], 5)105        if np.isfinite(raw[1]) and np.isfinite(fresh_xc[1]) else None,106        "p_fresh_+1": bootstrap_pvalue(boot_p1, 0.0),107        "p_fresh_-1": bootstrap_pvalue(boot_m1, 0.0),108    }109110111def main() -> None:112    run_utc = datetime.now(UTC)113    client = HFMarketDataClient()114    random10 = random_stock_universe(115        client.tickers("stock", timeframe="1min", adjustment=ADJ))[:10]116117    series_spec: dict[str, tuple[str, str, str | None]] = (118        {"SPY": ("etf", "SPY", ADJ), "SPX": ("index", "SPX", None)}119        | {t: ("stock", t, ADJ) for t in LIQUID if t != "QQQ"}120        | {"QQQ": ("etf", "QQQ", ADJ)}121        | {t: ("etf", t, ADJ) for t in SECTORS}122        | {t: ("stock", t, ADJ) for t in random10}123        | {f"ES[{a}]": ("futures", "ES", a)124           for a in ("contin_adj_ratio", "contin_adj_absolute", "contin_UNadj")}125    )126127    pairs: list[dict] = []128    for window, (s, e) in WINDOWS.items():129        grids: dict[str, dict] = {}130        for name, (asset, ticker, adj) in series_spec.items():131            if window == "2006-2007" and (name == "SPX" or name.startswith("ES[")):132                continue133            bars = client.get_bars(asset, ticker, "1min", adj, s, e)134            g = rth_day_grids(bars)135            if len(g) >= MIN_DAYS:136                grids[name] = g137            print(f"{window} {name}: {len(g)} days")138        if "SPY" not in grids:139            continue140        day_list = sorted(grids["SPY"].keys())141        mats = {n: day_matrix(g, day_list) for n, g in grids.items()}142        rx, fx = mats["SPY"]143        for name in mats:144            if name == "SPY":145                continue146            ry, fy = mats[name]147            # convention: x = the hypothesized leader.148            # ES/SPX pairs: x = ES or SPX, y = SPY. Others: x = SPY, y = ticker.149            if name.startswith("ES[") or name == "SPX":150                m = analyze_pair(ry, fy, rx, fx)151            else:152                m = analyze_pair(rx, fx, ry, fy)153            if m is None:154                continue155            bucket = ("es" if name.startswith("ES[") else "index" if name == "SPX"156                      else "sector" if name in SECTORS157                      else "liquid" if name in LIQUID else "random")158            m |= {"pair": f"{name}->SPY" if bucket in ("es", "index") else f"SPY->{name}",159                  "bucket": bucket, "window": window}160            pairs.append(m)161        print(f"{window}: {len(pairs)} pair-cells so far")162163    # FDR over all fresh ±1 tests jointly164    tests = [(i, "p_fresh_+1") for i in range(len(pairs))] + \165            [(i, "p_fresh_-1") for i in range(len(pairs))]166    mask = benjamini_hochberg(np.array([pairs[i][k] for i, k in tests]), alpha=0.05)167    for (i, k), rej in zip(tests, mask, strict=True):168        pairs[i][f"fdr_{k[2:]}"] = bool(rej)169170    fdr_survivors = [171        {"pair": p["pair"], "window": p["window"], "bucket": p["bucket"],172         "fresh_+1": p["fresh_xcorr"]["1"], "fresh_-1": p["fresh_xcorr"]["-1"],173         "which": [k for k in ("fresh_+1", "fresh_-1") if p[f"fdr_{k}"]]}174        for p in pairs if p.get("fdr_fresh_+1") or p.get("fdr_fresh_-1")175    ]176    summary = {177        "pair_cells": len(pairs),178        "tests": len(tests),179        "fdr_alpha": 0.05,180        "fdr_survivors": fdr_survivors,181        "median_artifact_share_by_bucket": {182            b: round(float(np.median(183                [p["artifact_share_+1"] for p in pairs184                 if p["bucket"] == b and p["artifact_share_+1"] is not None])), 5)185            for b in ("liquid", "sector", "random")186        },187    }188189    results = {190        "experiment": "expD_leadlag_scan",191        "run_utc": run_utc.isoformat(),192        "author": "Simon-Pierre Boucher",193        "contact": "contact@spboucher.ai",194        "data_source": "hfmarketdata.io",195        "confidence_level": 0,196        "protocol": {"windows": WINDOWS, "adjustment": ADJ, "max_lag": MAX_LAG,197                     "boot": [BOOT_N, BOOT_SEED], "random10": random10,198                     "min_days": MIN_DAYS, "min_fresh": MIN_FRESH},199        "summary": summary,200        "pairs": pairs,201        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},202        "manifest": collect_manifest(),203    }204    out_dir = REPO_ROOT / "results" / "expD_leadlag_scan" / run_utc.strftime("%Y%m%dT%H%M%SZ")205    out_dir.mkdir(parents=True)206    (out_dir / "results.json").write_text(json.dumps(results, indent=2) + "\n")207    print(f"\nwrote {out_dir.relative_to(REPO_ROOT)}/results.json")208    print(json.dumps(summary, indent=1))209210211if __name__ == "__main__":212    main()213