Code / experiments/micro/expD_leadlag_scan/benchmark.py
experiments/micro/expD_leadlag_scan/benchmark.py
213 lines
# =============================================================================
# Project : anomaly-atlas
# File : experiments/micro/expD_leadlag_scan/benchmark.py
# Purpose : Lead-lag scan — raw LOCF vs both-fresh, artifact share, FDR
# Author : Simon-Pierre Boucher
# Contact : contact@spboucher.ai
# Data src : hfmarketdata.io (sole data source)
# Created : 2026-08-12
# Modified : 2026-08-12
# Platform : macOS / Apple Silicon (arm64)
# License : All rights reserved (research code)
# =============================================================================
"""Experiment D — lead-lag scan (protocol pre-specified in hypothesis.md).
The design measures T3 directly: every pair reports its raw-LOCF lead, its
both-fresh lead, and the difference (artifact share). Level 0 throughout.
"""
from __future__ import annotations
import json
import sys
from datetime import UTC, datetime
from pathlib import Path
import numpy as np
REPO_ROOT = Path(__file__).resolve().parents[3]
sys.path.insert(0, str(REPO_ROOT / "benchmarks"))
sys.path.insert(0, str(REPO_ROOT / "src"))
from hardware_manifest import collect_manifest # noqa: E402
from anomaly_atlas.data.cleaning import ( # noqa: E402
RTH_SLOTS,
both_fresh,
nan_xcorr,
rth_day_grids,
)
from anomaly_atlas.data.hf_client import HFMarketDataClient # noqa: E402
from anomaly_atlas.data.universe import random_stock_universe # noqa: E402
from anomaly_atlas.stats.multiple_testing import ( # noqa: E402
benjamini_hochberg,
bootstrap_pvalue,
)
ADJ = "adj_split"
LIQUID = ["AAPL", "MSFT", "NVDA", "AMZN", "GOOGL", "META", "TSLA", "JPM", "XOM", "UNH", "QQQ"]
SECTORS = ["XLF", "XLE", "XLK", "XLV", "XLI", "XLY", "XLP", "XLU", "XLB"]
WINDOWS = {"2006-2007": ("2006-01-01", "2008-01-01"),
"2014-2015": ("2014-01-01", "2016-01-01")}
MAX_LAG = 3
BOOT_N, BOOT_SEED = 200, 42
MIN_DAYS, MIN_FRESH = 200, 10_000
def day_matrix(days: dict[str, np.ndarray], day_list: list[str]) -> tuple[np.ndarray, np.ndarray]:
"""(n_days, 389) LOCF returns + interval-fresh masks, day-aligned."""
rets = np.full((len(day_list), RTH_SLOTS - 1), np.nan)
fresh = np.zeros((len(day_list), RTH_SLOTS - 1), dtype=bool)
for i, day in enumerate(day_list):
p = days.get(day)
if p is None:
continue
observed = np.isfinite(p)
filled = p.copy()
for t in range(1, RTH_SLOTS):
if not observed[t]:
filled[t] = filled[t - 1]
rets[i] = np.diff(filled)
fresh[i] = observed[1:] & observed[:-1]
return rets, fresh
def day_block_boot_xcorr(
rx: np.ndarray, fx: np.ndarray, ry: np.ndarray, fy: np.ndarray, lag: int
) -> np.ndarray:
"""Bootstrap distribution of the both-fresh xcorr at `lag`, resampling days."""
rng = np.random.default_rng(BOOT_SEED)
n_days = rx.shape[0]
out = np.empty(BOOT_N)
for b in range(BOOT_N):
idx = rng.integers(0, n_days, n_days)
bx, by = both_fresh(rx[idx].ravel(), fx[idx].ravel(), ry[idx].ravel(), fy[idx].ravel())
out[b] = nan_xcorr(bx, by, abs(lag)).get(lag, np.nan)
return out
def analyze_pair(rx, fx, ry, fy) -> dict | None:
bx, by = both_fresh(rx.ravel(), fx.ravel(), ry.ravel(), fy.ravel())
n_fresh = int(np.isfinite(bx).sum())
if n_fresh < MIN_FRESH:
return None
raw = nan_xcorr(rx.ravel(), ry.ravel(), MAX_LAG)
fresh_xc = nan_xcorr(bx, by, MAX_LAG)
boot_p1 = day_block_boot_xcorr(rx, fx, ry, fy, 1)
boot_m1 = day_block_boot_xcorr(rx, fx, ry, fy, -1)
return {
"n_fresh_pairs": n_fresh,
"fresh_fraction": round(n_fresh / rx.size, 4),
"raw_xcorr": {str(k): round(v, 5) if np.isfinite(v) else None for k, v in raw.items()},
"fresh_xcorr": {str(k): round(v, 5) if np.isfinite(v) else None
for k, v in fresh_xc.items()},
"artifact_share_+1": round(raw[1] - fresh_xc[1], 5)
if np.isfinite(raw[1]) and np.isfinite(fresh_xc[1]) else None,
"p_fresh_+1": bootstrap_pvalue(boot_p1, 0.0),
"p_fresh_-1": bootstrap_pvalue(boot_m1, 0.0),
}
def main() -> None:
run_utc = datetime.now(UTC)
client = HFMarketDataClient()
random10 = random_stock_universe(
client.tickers("stock", timeframe="1min", adjustment=ADJ))[:10]
series_spec: dict[str, tuple[str, str, str | None]] = (
{"SPY": ("etf", "SPY", ADJ), "SPX": ("index", "SPX", None)}
| {t: ("stock", t, ADJ) for t in LIQUID if t != "QQQ"}
| {"QQQ": ("etf", "QQQ", ADJ)}
| {t: ("etf", t, ADJ) for t in SECTORS}
| {t: ("stock", t, ADJ) for t in random10}
| {f"ES[{a}]": ("futures", "ES", a)
for a in ("contin_adj_ratio", "contin_adj_absolute", "contin_UNadj")}
)
pairs: list[dict] = []
for window, (s, e) in WINDOWS.items():
grids: dict[str, dict] = {}
for name, (asset, ticker, adj) in series_spec.items():
if window == "2006-2007" and (name == "SPX" or name.startswith("ES[")):
continue
bars = client.get_bars(asset, ticker, "1min", adj, s, e)
g = rth_day_grids(bars)
if len(g) >= MIN_DAYS:
grids[name] = g
print(f"{window} {name}: {len(g)} days")
if "SPY" not in grids:
continue
day_list = sorted(grids["SPY"].keys())
mats = {n: day_matrix(g, day_list) for n, g in grids.items()}
rx, fx = mats["SPY"]
for name in mats:
if name == "SPY":
continue
ry, fy = mats[name]
# convention: x = the hypothesized leader.
# ES/SPX pairs: x = ES or SPX, y = SPY. Others: x = SPY, y = ticker.
if name.startswith("ES[") or name == "SPX":
m = analyze_pair(ry, fy, rx, fx)
else:
m = analyze_pair(rx, fx, ry, fy)
if m is None:
continue
bucket = ("es" if name.startswith("ES[") else "index" if name == "SPX"
else "sector" if name in SECTORS
else "liquid" if name in LIQUID else "random")
m |= {"pair": f"{name}->SPY" if bucket in ("es", "index") else f"SPY->{name}",
"bucket": bucket, "window": window}
pairs.append(m)
print(f"{window}: {len(pairs)} pair-cells so far")
# FDR over all fresh ±1 tests jointly
tests = [(i, "p_fresh_+1") for i in range(len(pairs))] + \
[(i, "p_fresh_-1") for i in range(len(pairs))]
mask = benjamini_hochberg(np.array([pairs[i][k] for i, k in tests]), alpha=0.05)
for (i, k), rej in zip(tests, mask, strict=True):
pairs[i][f"fdr_{k[2:]}"] = bool(rej)
fdr_survivors = [
{"pair": p["pair"], "window": p["window"], "bucket": p["bucket"],
"fresh_+1": p["fresh_xcorr"]["1"], "fresh_-1": p["fresh_xcorr"]["-1"],
"which": [k for k in ("fresh_+1", "fresh_-1") if p[f"fdr_{k}"]]}
for p in pairs if p.get("fdr_fresh_+1") or p.get("fdr_fresh_-1")
]
summary = {
"pair_cells": len(pairs),
"tests": len(tests),
"fdr_alpha": 0.05,
"fdr_survivors": fdr_survivors,
"median_artifact_share_by_bucket": {
b: round(float(np.median(
[p["artifact_share_+1"] for p in pairs
if p["bucket"] == b and p["artifact_share_+1"] is not None])), 5)
for b in ("liquid", "sector", "random")
},
}
results = {
"experiment": "expD_leadlag_scan",
"run_utc": run_utc.isoformat(),
"author": "Simon-Pierre Boucher",
"contact": "contact@spboucher.ai",
"data_source": "hfmarketdata.io",
"confidence_level": 0,
"protocol": {"windows": WINDOWS, "adjustment": ADJ, "max_lag": MAX_LAG,
"boot": [BOOT_N, BOOT_SEED], "random10": random10,
"min_days": MIN_DAYS, "min_fresh": MIN_FRESH},
"summary": summary,
"pairs": pairs,
"client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},
"manifest": collect_manifest(),
}
out_dir = REPO_ROOT / "results" / "expD_leadlag_scan" / run_utc.strftime("%Y%m%dT%H%M%SZ")
out_dir.mkdir(parents=True)
(out_dir / "results.json").write_text(json.dumps(results, indent=2) + "\n")
print(f"\nwrote {out_dir.relative_to(REPO_ROOT)}/results.json")
print(json.dumps(summary, indent=1))
if __name__ == "__main__":
main()