Honesty doctrine. Every candidate anomaly is an artifact until proven otherwise; in-sample results are never findings; past statistical regularity does not imply future returns. This is research on statistical properties of market data — not investment advice, not a trading system.

Code / experiments/micro/expB_artifact_baselines/benchmark.py

experiments/micro/expB_artifact_baselines/benchmark.py 271 lines
# =============================================================================
#  Project   : anomaly-atlas
#  File      : experiments/micro/expB_artifact_baselines/benchmark.py
#  Purpose   : Measure the artifact nulls: bounce, staleness, LOCF lead-lag
#  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 B — artifact baselines on real data (pre-specified protocol in
hypothesis.md; detectors gated on synthetic ground truth first, §8.1).

Every number produced here is a NULL LEVEL (Level 0 by construction): the
fake-signal magnitude that later experiments must exceed before claiming
anything. Universe, window, seeds are pre-specified; all data flows through
the cached hf_client.
"""

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.hf_client import HFMarketDataClient  # noqa: E402
from anomaly_atlas.stats.bootstrap import moving_block_bootstrap, percentile_ci  # noqa: E402
from anomaly_atlas.stats.reversion import ac1, variance_ratio  # noqa: E402
from anomaly_atlas.validation.artifacts import roll_spread  # noqa: E402

LIQUID_STOCK = ["AAPL", "MSFT", "NVDA", "AMZN", "GOOGL", "META", "TSLA", "JPM", "XOM", "UNH"]
LIQUID_ETF = ["SPY", "QQQ"]
N_RANDOM, RANDOM_SEED = 30, 42
START, END = "2024-01-02", "2024-04-01"
ADJ = "adj_split"
BOOT_N, BOOT_SEED = 300, 42
MAX_LAG = 3

RTH_MINUTES = [f"{h:02d}:{m:02d}" for h in range(9, 16) for m in range(60)]
RTH_MINUTES = [t for t in RTH_MINUTES if "09:30" <= t < "16:00"]  # 390 slots
SLOT = {t: i for i, t in enumerate(RTH_MINUTES)}


def rth_day_grids(bars: list[dict]) -> dict[str, np.ndarray]:
    """day -> 390-slot array of log close prices (NaN where no print)."""
    days: dict[str, np.ndarray] = {}
    for b in bars:
        dt = b["datetime"]
        t = dt[11:16]
        if not ("09:30" <= t < "16:00"):
            continue
        grid = days.setdefault(dt[:10], np.full(390, np.nan))
        grid[SLOT[t]] = np.log(b["close"])
    return days


def trade_time_returns(days: dict[str, np.ndarray]) -> np.ndarray:
    """Within-day log returns between consecutive PRINTS (no grid, no LOCF)."""
    out = []
    for day in sorted(days):
        p = days[day]
        obs = p[np.isfinite(p)]
        if len(obs) >= 2:
            out.append(np.diff(obs))
    return np.concatenate(out) if out else np.array([])


def locf_grid_returns(days: dict[str, np.ndarray], day_list: list[str]) -> np.ndarray:
    """Concatenated per-day LOCF grid returns, NaN before first print and at
    day boundaries — the join that MANUFACTURES the stale-price artifact."""
    out = []
    for day in day_list:
        p = days.get(day)
        if p is None:
            out.append(np.full(389, np.nan))
            continue
        filled = p.copy()
        for i in range(1, 390):
            if not np.isfinite(filled[i]):
                filled[i] = filled[i - 1]
        out.append(np.diff(filled))  # NaN propagates before first print
    return np.concatenate(out)


def nan_xcorr(x: np.ndarray, y: np.ndarray, max_lag: int) -> dict[int, float]:
    """corr(x_{t-k}, y_t) over finite pairs only; k>0 = x leads y."""
    n = min(len(x), len(y))
    x, y = x[:n], y[:n]
    out: dict[int, float] = {}
    for k in range(-max_lag, max_lag + 1):
        a = x[: n - k] if k >= 0 else x[-k:]
        b = y[k:] if k >= 0 else y[: n + k]
        m = np.isfinite(a) & np.isfinite(b)
        if m.sum() < 100 or a[m].std() == 0 or b[m].std() == 0:
            out[k] = float("nan")
            continue
        out[k] = float(np.corrcoef(a[m], b[m])[0, 1])
    return out


def analyze_ticker(days: dict[str, np.ndarray], day_list: list[str]) -> dict | None:
    present = (
        np.concatenate([np.isfinite(days[d]) for d in day_list if d in days])
        if any(d in days for d in day_list)
        else np.array([])
    )
    n_days_covered = sum(d in days for d in day_list)
    if n_days_covered < 30:
        return None
    r = trade_time_returns(days)
    if len(r) < 2_000:
        return None
    block = max(50, len(r) // max(n_days_covered, 1))
    boot = moving_block_bootstrap(r, ac1, block=block, n_boot=BOOT_N, seed=BOOT_SEED)
    lo, hi = percentile_ci(boot)
    spread = roll_spread(r)
    return {
        "days_covered": n_days_covered,
        "staleness": round(1.0 - present.mean() * len(present) / (390 * n_days_covered), 4)
        if n_days_covered
        else None,
        "rth_fill_ratio": round(present.sum() / (390 * n_days_covered), 4),
        "n_trade_returns": int(len(r)),
        "ac1": round(ac1(r), 5),
        "ac1_ci95": [round(lo, 5), round(hi, 5)],
        "roll_rel_spread": round(spread, 6) if np.isfinite(spread) else None,
        "vr5": round(variance_ratio(r, 5), 4),
        "vr30": round(variance_ratio(r, 30), 4),
    }


def main() -> None:
    run_utc = datetime.now(UTC)
    client = HFMarketDataClient()

    # deterministic random universe (seed pre-specified)
    all_stock = client.tickers("stock", timeframe="1min", adjustment=ADJ)
    rng = np.random.default_rng(RANDOM_SEED)
    random_universe = sorted(rng.choice(sorted(all_stock), N_RANDOM, replace=False))
    universe = (
        [("stock", t, "liquid") for t in LIQUID_STOCK]
        + [("etf", t, "liquid") for t in LIQUID_ETF]
        + [("stock", t, "random") for t in random_universe]
    )

    # fetch + grid everything
    grids: dict[str, dict[str, np.ndarray]] = {}
    for asset, ticker, _ in universe:
        bars = client.get_bars(asset, ticker, "1min", ADJ, START, END)
        grids[ticker] = rth_day_grids(bars)
        print(f"{ticker}: {sum(len(v[np.isfinite(v)]) for v in grids[ticker].values())} RTH bars")
    day_list = sorted(grids["SPY"].keys())  # trading calendar := SPY days

    # B1-B3: per-ticker artifact levels
    per_ticker: dict[str, dict] = {}
    for asset, ticker, bucket in universe:
        m = analyze_ticker(grids[ticker], day_list)
        if m is not None:
            m["bucket"] = bucket
            m["asset"] = asset
            per_ticker[ticker] = m

    # B4: LOCF lead-lag vs SPY
    spy_r = locf_grid_returns(grids["SPY"], day_list)
    for ticker, m in per_ticker.items():
        if ticker == "SPY":
            continue
        r = locf_grid_returns(grids[ticker], day_list)
        xc = nan_xcorr(spy_r, r, MAX_LAG)
        m["xcorr_vs_spy"] = {str(k): round(v, 5) if np.isfinite(v) else None for k, v in xc.items()}
        m["spy_leads_+1"] = round(xc[1], 5) if np.isfinite(xc[1]) else None

    # SPX (index) vs SPY — the non-synchronous-session case
    spx_bars = client.get_bars("index", "SPX", "1min", None, START, END)
    spx_grid = rth_day_grids(spx_bars)
    spx_r = locf_grid_returns(spx_grid, day_list)
    spx_xc = nan_xcorr(spx_r, spy_r, MAX_LAG)

    # staleness -> artifact monotonicity (Spearman)
    pairs = [
        (m["staleness"], m["spy_leads_+1"])
        for m in per_ticker.values()
        if m.get("spy_leads_+1") is not None and m["staleness"] is not None
    ]
    xs = np.array([p[0] for p in pairs])
    ys = np.array([p[1] for p in pairs])
    rx = np.argsort(np.argsort(xs)).astype(float)
    ry = np.argsort(np.argsort(ys)).astype(float)
    spearman = float(np.corrcoef(rx, ry)[0, 1]) if len(pairs) > 5 else float("nan")

    # aggregates by staleness tercile
    stale_vals = sorted(m["staleness"] for m in per_ticker.values())
    t1, t2 = np.percentile(stale_vals, [33.3, 66.7])

    def tercile(s: float) -> str:
        return "fresh" if s <= t1 else "mid" if s <= t2 else "stale"

    agg: dict[str, dict] = {}
    for name in ("fresh", "mid", "stale"):
        rows = [m for m in per_ticker.values() if tercile(m["staleness"]) == name]
        if not rows:
            continue
        agg[name] = {
            "n": len(rows),
            "median_staleness": round(float(np.median([m["staleness"] for m in rows])), 4),
            "median_ac1": round(float(np.median([m["ac1"] for m in rows])), 5),
            "median_roll_spread": round(
                float(np.median([m["roll_rel_spread"] for m in rows if m["roll_rel_spread"]])), 6
            ),
            "median_vr5": round(float(np.median([m["vr5"] for m in rows])), 4),
            "median_vr30": round(float(np.median([m["vr30"] for m in rows])), 4),
            "median_spy_leads_+1": round(
                float(
                    np.median(
                        [m["spy_leads_+1"] for m in rows if m.get("spy_leads_+1") is not None]
                    )
                ),
                5,
            ),
        }

    results = {
        "experiment": "expB_artifact_baselines",
        "run_utc": run_utc.isoformat(),
        "author": "Simon-Pierre Boucher",
        "contact": "contact@spboucher.ai",
        "data_source": "hfmarketdata.io",
        "protocol": {
            "window": [START, END],
            "adjustment": ADJ,
            "rth": "09:30-16:00",
            "liquid": LIQUID_STOCK + LIQUID_ETF,
            "random_universe": list(random_universe),
            "random_seed": RANDOM_SEED,
            "boot": [BOOT_N, BOOT_SEED],
            "confidence_level": 0,
            "note": "artifact NULL levels — descriptive, in-sample by design",
        },
        "per_ticker": per_ticker,
        "terciles": {"cuts": [round(float(t1), 4), round(float(t2), 4)], "agg": agg},
        "staleness_vs_spy_lead_spearman": round(spearman, 4),
        "spx_vs_spy_xcorr": {
            str(k): round(v, 5) if np.isfinite(v) else None for k, v in spx_xc.items()
        },
        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},
        "manifest": collect_manifest(),
    }

    out_dir = REPO_ROOT / "results" / "expB_artifact_baselines" / 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("terciles:", json.dumps(agg, indent=1))
    print("spearman(staleness, SPY leads +1):", round(spearman, 4))
    print("SPX vs SPY xcorr:", results["spx_vs_spy_xcorr"])


if __name__ == "__main__":
    main()