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

experiments/micro/expG_cost_frontier/benchmark.py 285 lines
# =============================================================================
#  Project   : anomaly-atlas
#  File      : experiments/micro/expG_cost_frontier/benchmark.py
#  Purpose   : Cost frontier: net-of-spread sweep over the double-filtered pool
#  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 G — the cost frontier (protocol pre-specified in hypothesis.md).

Pool = mechanical intersection of committed expC/expD/expF outputs. For each
rule: gross daily stream + daily turnover, then net = gross − κ·(EDGE/2)·
turnover across the declared κ sweep, with the analytic breakeven κ*.
"""

from __future__ import annotations

import glob
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.cleaning import RTH_SLOTS, rth_day_grids  # noqa: E402
from anomaly_atlas.data.hf_client import HFMarketDataClient  # noqa: E402
from anomaly_atlas.data.universe import TRAIN_SUBPERIODS  # noqa: E402
from anomaly_atlas.stats.bootstrap import moving_block_bootstrap, percentile_ci  # noqa: E402
from anomaly_atlas.validation.artifacts import edge_spread  # noqa: E402

ADJ = "adj_split"
KAPPAS = [0.0, 0.1, 0.25, 0.5, 1.0, 2.0]
ONE_MIN_WINDOW = ("2014-01-01", "2016-01-01")
D_WINDOW = ("2014-01-01", "2016-01-01")


def latest(pattern: str) -> dict:
    return json.loads(Path(sorted(glob.glob(pattern))[-1]).read_text())


def pool_from_committed_results() -> tuple[list[dict], list[dict]]:
    """Rebuild the double-filtered pool mechanically (no hand-picking)."""
    rc = latest(str(REPO_ROOT / "results/expC_reversion_scan/*/results.json"))
    rf = latest(str(REPO_ROOT / "results/expF_multiple_testing/*/results.json"))
    rd = latest(str(REPO_ROOT / "results/expD_leadlag_scan/*/results.json"))

    cells = rc["cells"]
    for c in cells:
        c["vr30_excess"] = c["vr30"] - (1 + 2 * c["ac1"] * (1 - 1 / 30))
    triage = {(c["ticker"], c["timeframe"], c["period"]) for c in cells
              if c.get("fdr_vr30") and c["vr30_excess"] < -0.05}
    spa_r = {tuple(s[3:].split(":")) for s in rf["funnel"]["spa_step1_survivors"]
             if s[1] == "R"}  # (ticker, timeframe)
    rev_pool = [{"ticker": t, "timeframe": tf, "period": per,
                 "asset": "etf" if t in ("SPY", "QQQ") else "stock"}
                for (t, tf, per) in sorted(triage) if (t, tf) in spa_r]

    ll_pool = []
    for p in rd["pairs"]:
        if p["window"] != "2014-2015" or p["bucket"] == "index":
            continue
        if not (p.get("fdr_fresh_+1") or p.get("fdr_fresh_-1")):
            continue
        follower = p["pair"].split("->")[1] if p["pair"].startswith("SPY->") else "SPY"
        leader = "SPY" if follower != "SPY" else p["pair"].split("->")[0].split("[")[0]
        sign = np.sign(p["fresh_xcorr"]["1"]) if p.get("fdr_fresh_+1") else np.sign(
            p["fresh_xcorr"]["-1"])
        ll_pool.append({"pair": p["pair"], "leader": leader, "follower": follower,
                        "sign": int(sign), "bucket": p["bucket"]})
    return rev_pool, ll_pool


def contrarian_stream(bars: list[dict], timeframe: str) -> dict[str, tuple[float, float]]:
    """day -> (gross rule return, turnover) for the +contrarian rule."""
    by_day: dict[str, list[float]] = defaultdict(list)
    for b in bars:
        dt = b["datetime"]
        if timeframe == "1day" or "09:30" <= dt[11:16] < "16:00":
            by_day[dt[:10]].append(np.log(b["close"]))
    days = sorted(by_day)
    out: dict[str, tuple[float, float]] = {}
    if timeframe == "1day":
        closes = np.array([by_day[d][0] for d in days])
        r = np.diff(closes)
        pos_prev = 0.0
        for i in range(1, len(r)):
            pos = -np.sign(r[i - 1])
            out[days[i + 1]] = (float(pos * r[i]), float(abs(pos - pos_prev)))
            pos_prev = pos
        return out
    for d in days:
        p = np.array(by_day[d])
        if len(p) < 3:
            continue
        r = np.diff(p)
        pos = -np.sign(r[:-1])
        gross = float(np.sum(pos * r[1:]))
        turnover = float(abs(pos[0]) + np.sum(np.abs(np.diff(pos))) + abs(pos[-1]))
        out[d] = (gross, turnover)
    return out


def leadlag_stream(gx: dict, gy: dict, day_list: list[str], sign: int
                   ) -> dict[str, tuple[float, float]]:
    """day -> (gross, turnover) for pos_t = sign * sign(x_{t-1}) on fresh minutes."""
    out: dict[str, tuple[float, float]] = {}
    for day in day_list:
        px, py = gx.get(day), gy.get(day)
        if px is None or py is None:
            continue
        ox, oy = np.isfinite(px), np.isfinite(py)
        fpx, fpy = px.copy(), py.copy()
        for t in range(1, RTH_SLOTS):
            if not ox[t]:
                fpx[t] = fpx[t - 1]
            if not oy[t]:
                fpy[t] = fpy[t - 1]
        rx, ry = np.diff(fpx), np.diff(fpy)
        fx = ox[1:] & ox[:-1]
        fy = oy[1:] & oy[:-1]
        keep = fx[:-1] & fy[1:] & np.isfinite(rx[:-1]) & np.isfinite(ry[1:])
        if keep.sum() < 30:
            continue
        pos = np.zeros(len(ry))
        pos[1:][keep] = sign * np.sign(rx[:-1][keep])
        gross = float(np.sum(pos[1:] * ry[1:]))
        turnover = float(np.sum(np.abs(np.diff(np.concatenate([[0.0], pos, [0.0]])))))
        out[day] = (gross, turnover)
    return out


def half_spread(client: HFMarketDataClient, asset: str, ticker: str,
                adj: str | None, start: str, end: str) -> float:
    bars = client.get_bars(asset, ticker, "1day", adj, start, end)
    if len(bars) < 100:
        return float("nan")
    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])
    s = edge_spread(o, h, lo, c)
    return s / 2.0 if np.isfinite(s) else float("nan")


def sanitize(obj):
    """Replace NaN/inf with None recursively — Python json emits bare NaN,
    which is invalid JSON for every other consumer (incl. the web platform)."""
    if isinstance(obj, dict):
        return {k: sanitize(v) for k, v in obj.items()}
    if isinstance(obj, list):
        return [sanitize(v) for v in obj]
    if isinstance(obj, float) and not np.isfinite(obj):
        return None
    return obj


def sweep(stream: dict[str, tuple[float, float]], hs: float) -> dict:
    days = sorted(d for d in stream if np.isfinite(stream[d][0]) and np.isfinite(stream[d][1]))
    gross = np.array([stream[d][0] for d in days])
    turn = np.array([stream[d][1] for d in days])
    res = {
        "n_days": len(days),
        "gross_mean_daily_bp": round(float(gross.mean()) * 1e4, 3),
        "turnover_per_day": round(float(turn.mean()), 2),
        "half_spread_bp": round(hs * 1e4, 3) if np.isfinite(hs) else None,
        "net": {},
    }
    if not np.isfinite(hs) or hs <= 0:
        res["kappa_star"] = None
        return res
    for k in KAPPAS:
        net = gross - k * hs * turn
        boot = moving_block_bootstrap(net, lambda x: float(np.mean(x)),
                                      block=21, n_boot=300, seed=42)
        lo_ci, hi_ci = percentile_ci(boot)
        res["net"][str(k)] = {
            "mean_daily_bp": round(float(net.mean()) * 1e4, 3),
            "ci95_bp": [round(lo_ci * 1e4, 3), round(hi_ci * 1e4, 3)],
            "positive": bool(net.mean() > 0),
        }
    denom = hs * turn.mean()
    res["kappa_star"] = round(float(gross.mean() / denom), 4) if denom > 0 else None
    return res


def main() -> None:
    run_utc = datetime.now(UTC)
    client = HFMarketDataClient()
    rev_pool, ll_pool = pool_from_committed_results()
    print("reversion pool:", [(r["ticker"], r["timeframe"], r["period"]) for r in rev_pool])
    print("leadlag pool:", [(p["pair"], p["sign"]) for p in ll_pool])

    items = []
    for r in rev_pool:
        if r["period"] in TRAIN_SUBPERIODS:
            s, e = TRAIN_SUBPERIODS[r["period"]]
        else:
            s, e = ONE_MIN_WINDOW
        bars = client.get_bars(r["asset"], r["ticker"], r["timeframe"], ADJ, s, e)
        stream = contrarian_stream(bars, r["timeframe"])
        hs = half_spread(client, r["asset"], r["ticker"], ADJ, s, e)
        item = {"rule": f"R:{r['ticker']}:{r['timeframe']}:{r['period']}",
                "family": "reversion"} | sweep(stream, hs)
        items.append(item)
        print(item["rule"], "kappa* =", item["kappa_star"])

    # lead-lag pool (2014-2015 window)
    s, e = D_WINDOW
    needed = {"SPY"} | {p["leader"] for p in ll_pool} | {p["follower"] for p in ll_pool}
    grids = {}
    for name in sorted(needed):
        if name == "ES":
            g = rth_day_grids(client.get_bars("futures", "ES", "1min",
                                              "contin_adj_ratio", s, e))
        else:
            asset = "etf" if name in ("SPY", "QQQ", "XLF", "XLE", "XLK", "XLV", "XLI",
                                      "XLY", "XLP", "XLU", "XLB") else "stock"
            g = rth_day_grids(client.get_bars(asset, name, "1min", ADJ, s, e))
        grids[name] = g
    day_list = sorted(grids["SPY"].keys())
    for p in ll_pool:
        stream = leadlag_stream(grids[p["leader"]], grids[p["follower"]],
                                day_list, p["sign"])
        traded = p["follower"]
        asset = ("futures" if traded == "ES" else
                 "etf" if traded in ("SPY", "QQQ") or traded.startswith("XL") else "stock")
        adj = "contin_adj_ratio" if traded == "ES" else ADJ
        hs = half_spread(client, asset, traded, adj, s, e)
        item = {"rule": f"L:{p['pair']}:{'+' if p['sign'] > 0 else '-'}",
                "family": "leadlag", "traded": traded} | sweep(stream, hs)
        items.append(item)
        print(item["rule"], "kappa* =", item["kappa_star"])

    ks = [i["kappa_star"] for i in items
          if i["kappa_star"] is not None and np.isfinite(i["kappa_star"])]
    summary = {
        "pool_size": len(items),
        "n_kappa_star_defined": len(ks),
        "kappa_star_median": round(float(np.median(ks)), 4) if ks else None,
        "kappa_star_max": round(float(max(ks)), 4) if ks else None,
        "survivors_at": {str(k): [i["rule"] for i in items
                                  if i["net"].get(str(k), {}).get("positive")]
                         for k in (0.1, 0.25, 0.5, 1.0)},
        "intraday_survivor_at_1.0": [
            i["rule"] for i in items
            if ":1day:" not in i["rule"] and i["net"].get("1.0", {}).get("positive")],
    }

    results = {
        "experiment": "expG_cost_frontier",
        "run_utc": run_utc.isoformat(),
        "author": "Simon-Pierre Boucher",
        "contact": "contact@spboucher.ai",
        "data_source": "hfmarketdata.io",
        "confidence_level": 0,
        "protocol": {"kappas": KAPPAS, "cost_model": "net = gross - k*(EDGE/2)*turnover",
                     "pool": "mechanical intersection of committed expC/expD/expF results"},
        "items": items,
        "summary": summary,
        "client_stats": vars(client.stats) | {"refreshes": list(client.stats.refreshes)},
        "manifest": collect_manifest(),
    }
    out_dir = REPO_ROOT / "results" / "expG_cost_frontier" / run_utc.strftime("%Y%m%dT%H%M%SZ")
    out_dir.mkdir(parents=True)
    (out_dir / "results.json").write_text(
        json.dumps(sanitize(results), indent=2, allow_nan=False) + "\n")
    print(f"\nwrote {out_dir.relative_to(REPO_ROOT)}/results.json")
    print(json.dumps(summary, indent=1))


if __name__ == "__main__":
    main()