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/expC_reversion_scan/benchmark.py

experiments/micro/expC_reversion_scan/benchmark.py 220 lines
# =============================================================================
#  Project   : anomaly-atlas
#  File      : experiments/micro/expC_reversion_scan/benchmark.py
#  Purpose   : Mean-reversion scan net of the EDGE bounce null (TRAIN only)
#  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 C — reversion scan (protocol pre-specified in hypothesis.md).

Every output cell is Level 0. The scan's job is triage: which cells show
excess reversion beyond the bounce null after FDR — those go to the expF
correction battery, nothing goes to the atlas from here.
"""

from __future__ import annotations

import json
import sys
from collections import defaultdict
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.data.universe import (  # noqa: E402
    LIQUID_ETF,
    LIQUID_STOCK,
    TRAIN_SUBPERIODS,
    core_universe,
)
from anomaly_atlas.stats.bootstrap import moving_block_bootstrap, percentile_ci  # noqa: E402
from anomaly_atlas.stats.multiple_testing import (  # noqa: E402
    benjamini_hochberg,
    bootstrap_pvalue,
)
from anomaly_atlas.stats.reversion import ac1, half_life, variance_ratio  # noqa: E402
from anomaly_atlas.validation.artifacts import edge_spread  # noqa: E402

ADJ = "adj_split"
BOOT_N, BOOT_SEED = 200, 42
TIMEFRAMES = ["1day", "30min", "5min"]
BARS_PER_DAY = {"5min": 78, "30min": 13, "1day": 1}
ONE_MIN_WINDOW = ("2014-01-01", "2016-01-01")  # liquid 12 only (declared)
MIN_RETURNS = {"1day": 350, "30min": 2_000, "5min": 5_000, "1min": 5_000}


def rth_returns_by_day(bars: list[dict], timeframe: str) -> np.ndarray:
    """Within-day log returns on RTH bars only (no overnight, no LOCF)."""
    days: dict[str, list[float]] = defaultdict(list)
    for b in bars:
        dt = b["datetime"]
        if timeframe == "1day":
            days[dt[:10]].append(np.log(b["close"]))
            continue
        if "09:30" <= dt[11:16] < "16:00":
            days[dt[:10]].append(np.log(b["close"]))
    if timeframe == "1day":
        allp = [v[0] for _, v in sorted(days.items())]
        return np.diff(allp) if len(allp) > 2 else np.array([])
    out = [np.diff(v) for _, v in sorted(days.items()) if len(v) >= 2]
    return np.concatenate(out) if out else np.array([])


def daily_ohlc(bars: list[dict]) -> tuple[np.ndarray, ...]:
    o = np.array([b["open"] for b in bars])
    h = np.array([b["high"] for b in bars])
    lo = np.array([b["low"] for b in bars])
    c = np.array([b["close"] for b in bars])
    return o, h, lo, c


def analyze_cell(r: np.ndarray, timeframe: str, spread: float) -> dict | None:
    if len(r) < MIN_RETURNS[timeframe]:
        return None
    var = r.var()
    if var == 0:
        return None
    bounce_ac1 = -(spread**2) / 4.0 / var if np.isfinite(spread) else 0.0
    block = max(20, BARS_PER_DAY.get(timeframe, 390) * 5)
    boot = moving_block_bootstrap(r, ac1, block=block, n_boot=BOOT_N, seed=BOOT_SEED)
    excess_boot = boot - bounce_ac1
    a = ac1(r)
    lo, hi = percentile_ci(excess_boot)
    boot_vr5 = moving_block_bootstrap(
        r, lambda x: variance_ratio(x, 5), block=block, n_boot=BOOT_N, seed=BOOT_SEED
    )
    boot_vr30 = moving_block_bootstrap(
        r, lambda x: variance_ratio(x, 30), block=block, n_boot=BOOT_N, seed=BOOT_SEED
    )
    return {
        "n": int(len(r)),
        "ac1": round(a, 5),
        "edge_spread": round(spread, 6) if np.isfinite(spread) else None,
        "bounce_ac1": round(bounce_ac1, 5),
        "excess_ac1": round(a - bounce_ac1, 5),
        "excess_ci95": [round(lo, 5), round(hi, 5)],
        "p_excess": bootstrap_pvalue(excess_boot, 0.0),
        "vr5": round(variance_ratio(r, 5), 4),
        "p_vr5": bootstrap_pvalue(boot_vr5, 1.0),
        "vr30": round(variance_ratio(r, 30), 4),
        "p_vr30": bootstrap_pvalue(boot_vr30, 1.0),
    }


def main() -> None:
    run_utc = datetime.now(UTC)
    client = HFMarketDataClient()
    universe = core_universe(client.tickers("stock", timeframe="1min", adjustment=ADJ))

    cells: list[dict] = []
    for asset, ticker, bucket in universe:
        # daily bars per sub-period: returns for 1day cells + EDGE spread input
        for sub, (s, e) in TRAIN_SUBPERIODS.items():
            day_bars = client.get_bars(asset, ticker, "1day", ADJ, s, e)
            if len(day_bars) < 200:
                continue
            o, h, lo, c = daily_ohlc(day_bars)
            spread = edge_spread(o, h, lo, c)
            hl = half_life(np.log(c))
            for tf in TIMEFRAMES:
                bars = day_bars if tf == "1day" else client.get_bars(asset, ticker, tf, ADJ, s, e)
                r = rth_returns_by_day(bars, tf)
                m = analyze_cell(r, tf, spread)
                if m is None:
                    continue
                m |= {"ticker": ticker, "bucket": bucket, "timeframe": tf, "period": sub}
                if tf == "1day":
                    m["half_life_days"] = round(hl, 1) if np.isfinite(hl) else None
                cells.append(m)
            print(f"{ticker} {sub}: done ({len(cells)} cells)")

    # 1min cells: liquid 12, declared window
    for asset, ticker in [("stock", t) for t in LIQUID_STOCK] + [("etf", t) for t in LIQUID_ETF]:
        s, e = ONE_MIN_WINDOW
        day_bars = client.get_bars(asset, ticker, "1day", ADJ, s, e)
        if len(day_bars) < 200:
            continue
        spread = edge_spread(*daily_ohlc(day_bars))
        bars = client.get_bars(asset, ticker, "1min", ADJ, s, e)
        r = rth_returns_by_day(bars, "1min")
        m = analyze_cell(r, "1min", spread)
        if m is not None:
            m |= {"ticker": ticker, "bucket": "liquid", "timeframe": "1min",
                  "period": "2014-2015"}
            cells.append(m)
        print(f"{ticker} 1min: done")

    # FDR within each statistic family, all cells jointly
    for key, pkey in [("excess_ac1", "p_excess"), ("vr5", "p_vr5"), ("vr30", "p_vr30")]:
        mask = benjamini_hochberg(np.array([c[pkey] for c in cells]), alpha=0.05)
        for c, rej in zip(cells, mask, strict=True):
            c[f"fdr_{key}"] = bool(rej)

    def survivors(key: str, sign_key: str, negative: bool) -> list[dict]:
        out = []
        for c in cells:
            if not c[f"fdr_{key}"]:
                continue
            v = c[sign_key] - (1.0 if sign_key.startswith("vr") else 0.0)
            if (v < 0) == negative:
                out.append({k: c[k] for k in ("ticker", "bucket", "timeframe", "period",
                                              sign_key, "edge_spread")})
        return out

    summary = {
        "cells_total": len(cells),
        "families_tested": 3,
        "fdr_alpha": 0.05,
        "excess_ac1_negative_survivors": survivors("excess_ac1", "excess_ac1", True),
        "excess_ac1_positive_survivors": survivors("excess_ac1", "excess_ac1", False),
        "vr30_below_1_survivors": len(survivors("vr30", "vr30", True)),
        "vr30_above_1_survivors": len(survivors("vr30", "vr30", False)),
        "liquid_2008_2015_intraday_negative": [
            c["ticker"] for c in cells
            if c["bucket"] == "liquid" and c["period"] == "2008-2015"
            and c["timeframe"] in ("5min", "30min")
            and c["fdr_excess_ac1"] and c["excess_ac1"] < 0
        ],
    }

    results = {
        "experiment": "expC_reversion_scan",
        "run_utc": run_utc.isoformat(),
        "author": "Simon-Pierre Boucher",
        "contact": "contact@spboucher.ai",
        "data_source": "hfmarketdata.io",
        "confidence_level": 0,
        "protocol": {
            "train_subperiods": TRAIN_SUBPERIODS, "adjustment": ADJ,
            "one_min_window": ONE_MIN_WINDOW, "boot": [BOOT_N, BOOT_SEED],
            "min_returns": MIN_RETURNS,
        },
        "summary": summary,
        "cells": cells,
        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},
        "manifest": collect_manifest(),
    }
    out_dir = REPO_ROOT / "results" / "expC_reversion_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({k: (v if not isinstance(v, list) else len(v))
                      for k, v in summary.items()}, indent=1))


if __name__ == "__main__":
    main()