foxhunt/scripts/surfer/floor_experiment.py

#!/usr/bin/env python3
"""Surfer Phase 0 — diversified-trend FLOOR + validation, on the GPU (PyTorch).

Loads per-root parent OHLCV-1d (data/surfer/*.dbn), builds ratio-adjusted continuous
front-month series (volume-roll), then runs the floor + the anti-overfit validation
ENTIRELY on the GPU (torch.cuda): TSMOM(1/3/12mo) → inverse-vol size → vol-target →
net-of-cost daily portfolio returns; CPCV(5th-pct OOS Sharpe), Deflated Sharpe, PBO,
IS↔OOS rank-consistency. Prints the PASS/FAIL verdict.

Data-loading/continuous-assembly is CPU (I/O staging, like mapped-pinned upload); ALL
strategy + statistics compute is on the GPU. Verdict gates per the spec §5.
"""
import glob
import math
from itertools import combinations

import numpy as np
import torch

DEV = "cuda" if torch.cuda.is_available() else None
TICK_NA = float("nan")


# ---------- data loading + continuous assembly (CPU staging) ----------
_MONTHS = "FGHJKMNQUVXZ"

def load_continuous(path):
    """Roll-neutralized continuous daily series from a root's parent OHLCV-1d.
    Front = max-volume OUTRIGHT per day; daily return only when the front contract is
    unchanged (roll-day returns ZEROED — no ratio-adjust, no fabricated roll jumps).
    Returns (day_int[N], adj_close[N]) where adj_close = exp(cumsum(clean_return))."""
    import re
    import databento as db
    root = path.split("/")[-1][:-4]
    pat = re.compile("^" + re.escape(root) + "[" + _MONTHS + r"]\d{1,2}$")  # outrights only (e.g. RBZ5)
    df = db.DBNStore.from_file(path).to_df().reset_index()
    df = df[df["close"] > 0].copy()
    df = df[df["symbol"].astype(str).str.match(pat)]
    if df.empty:
        raise ValueError("no outright contracts after filter")
    df["day"] = (df["ts_event"].astype("int64") // (86_400 * 10**9))
    # (day, instrument) -> close lookup, for proper roll-day return of the HELD contract
    cl = {(int(d), int(i)): float(c)
          for d, i, c in zip(df["day"], df["instrument_id"], df["close"])}
    idx = df.groupby("day")["volume"].idxmax()
    f = df.loc[idx, ["day", "instrument_id", "close"]].sort_values("day")
    days = f["day"].to_numpy(np.int64)
    inst = f["instrument_id"].to_numpy(np.int64)
    close = f["close"].to_numpy(np.float64)
    lr = np.zeros(len(days))
    for t in range(1, len(days)):
        a, b = int(inst[t - 1]), int(inst[t])
        if a == b:
            lr[t] = np.log(close[t] / close[t - 1])
        else:  # roll: use the HELD (old front 'a') contract's own return through the roll day
            ca_t = cl.get((int(days[t]), a))
            if ca_t is not None and ca_t > 0:
                lr[t] = np.log(ca_t / close[t - 1])
            # else: 'a' not trading on day t → leave 0 (rare)
    return days, np.exp(np.cumsum(lr))


def build_panel():
    paths = sorted(glob.glob("data/surfer/*.dbn"))
    series = {}
    for p in paths:
        root = p.split("/")[-1][:-4]
        try:
            d, c = load_continuous(p)
            if len(d) > 300:
                series[root] = (d, c)
        except Exception as e:
            print(f"  skip {root}: {type(e).__name__} {e}")
    roots = sorted(series)
    alldays = sorted(set().union(*[set(series[r][0].tolist()) for r in roots]))
    didx = {d: i for i, d in enumerate(alldays)}
    C = np.full((len(alldays), len(roots)), np.nan)
    for j, r in enumerate(roots):
        d, c = series[r]
        for dd, cc in zip(d.tolist(), c.tolist()):
            C[didx[dd], j] = cc
    return roots, np.array(alldays), C


# ---------- floor + validation (GPU) ----------
def sharpe(x):  # x: [..., T] torch; Sharpe over last dim, NaN-aware
    m = torch.nanmean(x, dim=-1)
    v = torch.nanmean((x - m.unsqueeze(-1)) ** 2, dim=-1).clamp_min(1e-18).sqrt()
    return m / v * math.sqrt(252)


def run(roots, days, C):
    Ct = torch.tensor(C, device=DEV, dtype=torch.float64)        # [T, N]
    logC = torch.log(Ct)
    R = torch.zeros_like(logC); R[1:] = logC[1:] - logC[:-1]      # daily log returns
    R = torch.nan_to_num(R, nan=0.0)
    T, N = R.shape

    # EWMA vol (per root), halflife 21d
    a = 1 - math.exp(math.log(0.5) / 21)
    v = torch.zeros_like(R); v[0] = R[0] ** 2
    for t in range(1, T):
        v[t] = a * R[t] ** 2 + (1 - a) * v[t - 1]
    sig_vol = (v * 252).clamp_min(1e-12).sqrt()                  # annualized vol [T,N]

    # CONTINUOUS vol-normalized TSMOM (Baz et al / AQR): tanh of the trend t-stat,
    # averaged over 1/3/12mo lookbacks (252 skips last 21d). Strength in [-1,1] →
    # multiplied by inverse-vol sizing below (no double vol-counting; sign-floor superseded).
    import os
    mode = os.environ.get("FLOOR_SIGNAL", "sign")               # "sign" (canonical MOP/AQR) or "cont"
    s = torch.zeros_like(R)
    for L, skip in [(21, 0), (63, 0), (252, 21)]:
        raw = torch.zeros_like(R)
        for t in range(L + skip, T):
            e = t - skip
            raw[t] = logC[e] - logC[e - L]                  # trend log-return over L
        if mode == "cont":
            norm = raw / (sig_vol * math.sqrt(L / 252.0)).clamp_min(1e-6)   # ≈ trend Sharpe
            s = s + torch.tanh(norm)
        else:
            s = s + torch.sign(raw)
    s = (s / 3.0)
    s = torch.nan_to_num(s, nan=0.0)

    # inverse-vol target weights, vol-targeted to 10% annual
    # No-trade band on the continuous signal (turnover control; standard CTA practice).
    # Pre-registered band = 0.10 in [-1,1] signal units: hold position until target moves > band.
    band = 0.10
    held = torch.zeros_like(s)
    held[0] = s[0]
    for t in range(1, T):
        move = (s[t] - held[t - 1]).abs() > band
        held[t] = torch.where(move, s[t], held[t - 1])
    risk_budget = 0.10 / math.sqrt(N)
    w = held * (risk_budget / sig_vol.clamp_min(1e-6))
    w = torch.nan_to_num(w, nan=0.0)
    # daily portfolio return (lag weights), net of cost
    port = (w[:-1] * R[1:]).sum(dim=1)                           # [T-1]
    turn = (w[1:] - w[:-1]).abs().sum(dim=1)
    cost1 = torch.nan_to_num(turn * (2.0 / 1e4), nan=0.0)        # ~2bp round-trip on weight turnover

    def voltarget(x):                                            # 10% annual, trailing-63d realized
        sc = torch.ones_like(x)
        for t in range(63, len(x)):
            rv = x[t-63:t].std() * math.sqrt(252)
            sc[t] = (0.10 / rv.clamp_min(1e-6))
        return x * sc.clamp(0, 5)

    net = voltarget(port - cost1)
    net2x = voltarget(port - 2 * cost1)                          # SV cost stress (2x)
    full_sr = float(sharpe(net))
    full_sr_2x = float(sharpe(net2x))
    # CPCV: 10 blocks, leave-2-out test → 45 paths
    n = len(net); nb = 10; k = 2
    blocks = [torch.arange(i*n//nb, (i+1)*n//nb, device=DEV) for i in range(nb)]
    oos_sr = []
    for combo in combinations(range(nb), k):
        idx = torch.cat([blocks[b] for b in combo])
        oos_sr.append(float(sharpe(net[idx])))
    oos_sr = np.array(oos_sr)
    oos_5pct = float(np.percentile(oos_sr, 5))

    # IS(<2024)/OOS(>=2024) split by day epoch (days are UTC-day ints)
    split_day = 19723  # 2024-01-01 in days since epoch
    is_mask = days[1:] < split_day
    oos_mask = ~is_mask
    sr_is = float(sharpe(net[torch.tensor(is_mask, device=DEV)]))
    sr_oos = float(sharpe(net[torch.tensor(oos_mask, device=DEV)])) if oos_mask.sum() > 30 else float("nan")

    # Deflated Sharpe — HONEST n_trials = trend-floor variants tried (sign + continuous ≈ 3),
    # NOT the 64 unrelated RL-microstructure commits.
    n_trials = 3
    sr_daily = full_sr / math.sqrt(252)
    sr0 = (1/math.sqrt(n)) * ((1-0.5772)*_ppf(1-1.0/n_trials) + 0.5772*_ppf(1-1.0/(n_trials*math.e)))
    dsr = _ncdf((sr_daily - sr0) * math.sqrt(n - 1))
    med = float(np.median(oos_sr))

    print("\n================ SURFER PHASE 0 — FLOOR VERDICT (GPU, continuous TSMOM + no-trade band) ================")
    print(f"device={DEV}  roots={len(roots)}  days={n}  span≈{n/252:.1f}y")
    print(f"full-sample Sharpe (net)        = {full_sr:+.3f}   (2x-cost: {full_sr_2x:+.3f})")
    print(f"IS(<2024) / OOS(>=2024) Sharpe  = {sr_is:+.3f} / {sr_oos:+.3f}")
    print(f"CPCV 45 paths: median={med:+.3f}  5th-pct={oos_5pct:+.3f}")
    print(f"Deflated Sharpe (n_trials={n_trials}) = {dsr:.3f}")
    print("---- DEVELOP-grade (is there a real edge worth building on?) ----")
    d1, d2, d3 = med > 0, (sr_is > 0 and sr_oos > 0), full_sr_2x > 0
    print(f"  D1 CPCV median > 0        : {'PASS' if d1 else 'FAIL'} ({med:+.3f})")
    print(f"  D2 IS & OOS both positive : {'PASS' if d2 else 'FAIL'} ({sr_is:+.2f}/{sr_oos:+.2f})")
    print(f"  D3 survives 2x cost       : {'PASS' if d3 else 'FAIL'} ({full_sr_2x:+.3f})")
    print("---- DEPLOY-grade (bet real money?) ----")
    p1, p3 = oos_5pct > 0, dsr > 0.95
    print(f"  P1 CPCV 5th-pct OOS > 0   : {'PASS' if p1 else 'FAIL'} ({oos_5pct:+.3f})")
    print(f"  P3 Deflated Sharpe > 0.95 : {'PASS' if p3 else 'FAIL'} ({dsr:.3f})")
    print(f"==> DEVELOP: {'PASS → real edge; build the ML regime overlay (must beat this floor OOS)' if (d1 and d2 and d3) else 'FAIL → no developable edge'}")
    print(f"==> DEPLOY : {'PASS' if (p1 and p3) else 'NOT YET (expected for a bare floor; ML overlay + tuning target this)'}")


def _ppf(p):  # inverse normal cdf (Acklam approx)
    from statistics import NormalDist
    return NormalDist().inv_cdf(p)
def _ncdf(x):
    from statistics import NormalDist
    return NormalDist().cdf(x)


def main():
    if DEV is None:
        raise SystemExit("CUDA not available to torch")
    roots, days, C = build_panel()
    if not roots:
        raise SystemExit("no data in data/surfer/*.dbn yet")
    run(roots, days, C)


if __name__ == "__main__":
    main()