#!/opt/anaconda3/bin/python
"""pn-001 runner — see 2026-08-22-pn001-orb-vxn.md (frozen design)."""
import math, glob
from collections import defaultdict
import pandas as pd

GLBX = "<data-root>/fp-data/glbx"
FRICTION = {"NQ": 0.5, "ES": 0.25}  # per side

# VXN prior-day close, tercile cuts frozen on the 2010-2015 subsample
vx = pd.read_csv("<data-root>/fp-data/fred/VXNCLS.csv")
vx.columns = [c.strip().lower() for c in vx.columns]
dcol = "date" if "date" in vx.columns else vx.columns[0]
vcol = [c for c in vx.columns if c != dcol][0]
vx[vcol] = pd.to_numeric(vx[vcol], errors="coerce")
vx = vx.dropna()
vx[dcol] = pd.to_datetime(vx[dcol]).dt.strftime("%Y-%m-%d")
vx = vx.sort_values(dcol).reset_index(drop=True)
vx["prior"] = vx[vcol].shift(1)
prior_map = dict(zip(vx[dcol], vx["prior"]))
calib = vx[(vx[dcol] >= "2010-01-01") & (vx[dcol] <= "2015-12-31")][vcol].dropna()
t1_cut, t2_cut = calib.quantile(1/3), calib.quantile(2/3)
print(f"VXN tercile cuts (2010-2015 calibration): {t1_cut:.1f} / {t2_cut:.1f}")

def quad_friday(day):
    d = pd.Timestamp(day)
    if d.month not in (3, 6, 9, 12): return False
    fridays = [x for x in pd.date_range(d.replace(day=1), periods=31) if x.month == d.month and x.weekday() == 4]
    return len(fridays) >= 3 and d.normalize() == fridays[2].normalize()

def grid5(hm): return int(hm[3:]) % 5 == 4

res = defaultdict(list)  # (sym, orw, tercile, exit) -> list of R
pts_res = defaultdict(list)
for sym in ("NQ", "ES"):
    for path in sorted(glob.glob(f"{GLBX}/{sym}/1m/20*.csv")):
        df = pd.read_csv(path, usecols=["ts_event", "open", "high", "low", "close"])
        ts = pd.to_datetime(df["ts_event"], utc=True).dt.tz_convert("America/New_York")
        df["day"] = ts.dt.strftime("%Y-%m-%d"); df["hm"] = ts.dt.strftime("%H:%M")
        df = df[(df["hm"] >= "09:30") & (df["hm"] <= "15:59")]
        for day, sd in df.groupby("day"):
            if quad_friday(day): continue
            v = prior_map.get(day)
            if v is None or not (v == v): continue
            ter = "T1-low" if v <= t1_cut else "T2-mid" if v <= t2_cut else "T3-high"
            bars = sd.sort_values("hm").reset_index(drop=True)
            if len(bars) < 200: continue
            hm = bars["hm"].tolist(); close = bars["close"].tolist()
            hi = bars["high"].tolist(); lo = bars["low"].tolist()
            last_i = max(i for i, h in enumerate(hm) if h <= "15:55")
            for orw, or_end in (("15m", "09:44"), ("30m", "09:59")):
                ors = [i for i, h in enumerate(hm) if h <= or_end]
                if not ors: continue
                or_hi = max(hi[i] for i in ors); or_lo = min(lo[i] for i in ors)
                width = or_hi - or_lo
                if width <= 0: continue
                ent = None
                for i, h in enumerate(hm):
                    if h <= or_end or h > "11:59" or not grid5(h): continue
                    if close[i] > or_hi: ent = (i, 1); break
                    if close[i] < or_lo: ent = (i, -1); break
                if not ent: continue
                i0, sgn = ent
                j60 = min(i0 + 60, last_i)
                for ex, j in (("60m", j60), ("1555", last_i)):
                    pts = sgn * (close[j] - close[i0])
                    res[(sym, orw, ter, ex)].append(pts / width)
                    pts_res[(sym, orw, ter, ex)].append(pts)

def tstat(xs):
    n = len(xs)
    if n < 3: return float("nan")
    mu = sum(xs)/n; sd = math.sqrt(sum((x-mu)**2 for x in xs)/(n-1))
    return mu / (sd / math.sqrt(n)) if sd > 0 else float("nan")

print("\nsym or  tercile  exit :   n    meanR   meanPts  net/trade($)  win%   t")
mult = {"NQ": 20, "ES": 50}
for (sym, orw, ter, ex), rs in sorted(res.items()):
    ps = pts_res[(sym, orw, ter, ex)]
    mu_r = sum(rs)/len(rs); mu_p = sum(ps)/len(ps)
    net = (mu_p - 2*FRICTION[sym]) * mult[sym]
    wr = 100*sum(1 for x in ps if x > 0)/len(ps)
    print(f"{sym} {orw} {ter} {ex:>4}: {len(rs):4d}  {mu_r:+.3f}R  {mu_p:+7.2f}  {net:+9.0f}  {wr:4.0f}%  t={tstat(ps):+.2f}")

import itertools
print("\nH1 (T3-high vs T1-low, Welch t on points):")
for sym, orw, ex in itertools.product(("NQ", "ES"), ("15m", "30m"), ("60m", "1555")):
    a = pts_res.get((sym, orw, "T3-high", ex), []); b = pts_res.get((sym, orw, "T1-low", ex), [])
    if len(a) > 2 and len(b) > 2:
        ma, mb = sum(a)/len(a), sum(b)/len(b)
        va = sum((x-ma)**2 for x in a)/(len(a)-1); vb = sum((x-mb)**2 for x in b)/(len(b)-1)
        se = math.sqrt(va/len(a)+vb/len(b))
        print(f"  {sym} {orw} {ex:>4}: high {ma:+6.2f} vs low {mb:+6.2f} pts · t={((ma-mb)/se if se>0 else float('nan')):+.2f}")