"""Open-factor sweep — pre-registered 2026-09-08 (studies/2026-09-08-open-factors-prereg.md).
Univariate quintile tests + walk-forward logistic vs the reconstructed 3-factor baseline."""
import os, sys, json, math, warnings
import numpy as np, pandas as pd
from scipy.stats import norm
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, log_loss
warnings.filterwarnings("ignore")
FP="<data-root>/fp-data"; CACHE="data"; os.makedirs(CACHE, exist_ok=True)
YEARS=range(2011, 2027)

def load_1m(sym):
    p=f"{CACHE}/{sym}_1m_et.parquet"
    if os.path.exists(p): return pd.read_parquet(p)
    fr=[]
    for y in YEARS:
        f=f"{FP}/glbx/{sym}/1m/{y}.csv"
        if not os.path.exists(f): continue
        d=pd.read_csv(f, usecols=["ts_event","open","high","low","close","volume"])
        fr.append(d)
    d=pd.concat(fr, ignore_index=True)
    d["ts"]=pd.to_datetime(d["ts_event"], utc=True).dt.tz_convert("America/New_York")
    d=d.drop(columns=["ts_event"]).sort_values("ts").reset_index(drop=True)
    # session date: bars from 18:00 belong to the NEXT calendar day's session
    t=d["ts"]; d["sday"]=(t + pd.Timedelta(hours=6)).dt.date  # 18:00 → next day
    d["min"]=t.dt.hour*60+t.dt.minute
    d.to_parquet(p); return d

def quad_dates(years):
    out=set()
    for y in years:
        for m in (3,6,9,12):
            d=pd.Timestamp(year=y, month=m, day=1); fr=[x for x in pd.date_range(d, d+pd.Timedelta(days=27)) if x.weekday()==4]
            out.add(fr[2].date())
    return out
QUAD=quad_dates(range(2010,2028))
def opex_week(day):
    d=pd.Timestamp(day); fr=[x for x in pd.date_range(d.replace(day=1), d.replace(day=1)+pd.Timedelta(days=27)) if x.weekday()==4][2]
    return abs((d - fr).days) <= 4 and d.weekday() <= fr.weekday() and (fr - d).days >= 0

def value_area(rth):
    if rth.empty: return (np.nan,)*3
    px=((rth["open"]+rth["close"])/2/5).round()*5; vol=rth["volume"].values
    prof=pd.Series(vol).groupby(px.values).sum().sort_index()
    if prof.empty: return (np.nan,)*3
    poc=prof.idxmax(); total=prof.sum(); acc=prof[poc]; lo=hi=list(prof.index).index(poc); idx=list(prof.index)
    while acc < 0.7*total and (lo>0 or hi<len(idx)-1):
        up = prof.iloc[hi+1] if hi<len(idx)-1 else -1; dn = prof.iloc[lo-1] if lo>0 else -1
        if up>=dn: hi+=1; acc+=up
        else: lo-=1; acc+=dn
    return poc, idx[hi], idx[lo]

def sessions(sym, need_vol=False):
    d=load_1m(sym); rows=[]
    for sday, g in d.groupby("sday"):
        rth=g[(g["min"]>=570)&(g["min"]<960)]; on=g[(g["min"]>=1080)|(g["min"]<570)]
        if len(rth)<300: continue
        r={"sday":sday, "o":rth.iloc[0]["open"], "c":rth.iloc[-1]["close"], "h":rth["high"].max(), "l":rth["low"].min(),
           "on_first": on.iloc[0]["open"] if len(on) else np.nan, "on_last": on.iloc[-1]["close"] if len(on) else np.nan,
           "on_h": on["high"].max() if len(on) else np.nan, "on_l": on["low"].min() if len(on) else np.nan}
        if sym=="NQ":
            def at(m):
                b=rth[rth["min"]==m]; return b.iloc[0]["close"] if len(b) else np.nan
            r.update({"c5":at(574),"c15":at(584),"c30":at(599)})
            last30=rth[rth["min"]>=930]; r["last30"]= (last30.iloc[-1]["close"]/last30.iloc[0]["open"]-1) if len(last30)>5 else np.nan
            r["poc"],r["vah"],r["val"]=value_area(rth)
        rows.append(r)
    return pd.DataFrame(rows).set_index("sday").sort_index()

print("loading…", flush=True)
NQ=sessions("NQ"); ES=sessions("ES"); ZN=sessions("ZN"); E6=sessions("6E")
vxn=pd.read_csv(f"{FP}/fred/VXNCLS.csv"); vxn.columns=["d","v"]; vxn["d"]=pd.to_datetime(vxn["d"]).dt.date; vxn["v"]=pd.to_numeric(vxn["v"],errors="coerce"); vxn=vxn.set_index("d")["v"].dropna()
vix=pd.read_csv(f"{FP}/fred/VIXCLS.csv"); vix.columns=["d","v"]; vix["d"]=pd.to_datetime(vix["d"]).dt.date; vix["v"]=pd.to_numeric(vix["v"],errors="coerce"); vix=vix.set_index("d")["v"].dropna()

df=NQ.copy(); df["prev_c"]=df["c"].shift(1); df["prev_h"]=df["h"].shift(1); df["prev_l"]=df["l"].shift(1); df["prev_poc"]=df["poc"].shift(1); df["prev_vah"]=df["vah"].shift(1); df["prev_val"]=df["val"].shift(1); df["prev_last30"]=df["last30"].shift(1)
df["gapdays"]=pd.Series(df.index).diff().dt.days.values
df["rng"]=df["h"]-df["l"]; df["atr20"]=df["rng"].shift(1).rolling(20).mean()
# targets
df["up5"]=(df["c5"]>df["o"]).astype(float).where(df["c5"].notna()); df["up15"]=(df["c15"]>df["o"]).astype(float).where(df["c15"].notna()); df["up30"]=(df["c30"]>df["o"]).astype(float).where(df["c30"].notna())
df["expand"]=(df["rng"]/df["atr20"]>=1.0).astype(float); df["trend"]=((df["c"]-df["o"]).abs()/df["rng"]>=0.6).astype(float)
# baseline
vx=vxn.reindex(df.index).ffill(); df["B1_vxn_shift"]=(vx.shift(1)/vx.shift(2)-1).values
df["B2_on_rates"]=(ZN["on_last"]/ZN["on_first"]-1).reindex(df.index).values
df["B3_on_dollar"]=-(E6["on_last"]/E6["on_first"]-1).reindex(df.index).values
# candidates
df["F1_gap"]=(df["o"]-df["prev_c"])/df["atr20"]
df["F2_on_accept"]=np.where(df["on_last"]>df["prev_vah"],1,np.where(df["on_last"]<df["prev_val"],-1,0)).astype(float)
df["F3_nq_vs_es"]=(df["on_last"]/df["on_first"]-1)-(ES["on_last"]/ES["on_first"]-1).reindex(df.index)
df["F4_vol_ratio"]=(vx.shift(1)/vix.reindex(df.index).ffill().shift(1)).values
df["F5_close_loc"]=(df["prev_c"]-df["prev_l"])/(df["prev_h"]-df["prev_l"])
df["F6_last30"]=df["prev_last30"]
df["F7_on_compress"]=(df["on_h"]-df["on_l"])/df["atr20"]
dow=pd.Series(pd.to_datetime(pd.Series(df.index)).dt.weekday.values, index=df.index); df["F8_dow"]=dow.astype(float)
df["F8_opex"]=pd.Series([1.0 if opex_week(d) else 0.0 for d in df.index], index=df.index)
nqret=df["c"]/df["prev_c"]-1; znret=(ZN["c"]/ZN["c"].shift(1)-1).reindex(df.index); df["F9_rates_corr"]=nqret.rolling(20).corr(znret).shift(1)
df["F10_on_poc_dist"]=(df["on_last"]-df["prev_poc"])/df["atr20"]
# exclusions
df=df[(df.index>=pd.Timestamp("2012-01-03").date())]; df=df[~pd.Series(df.index).isin(QUAD).values]; df=df[df["gapdays"]<=4]
df=df.dropna(subset=["atr20","up5","up15","up30"])
print(f"sessions after exclusions: {len(df)}  {df.index[0]} → {df.index[-1]}", flush=True)
TARGETS=["up5","up15","up30","expand","trend"]; BASE=["B1_vxn_shift","B2_on_rates","B3_on_dollar"]
CANDS=["F1_gap","F2_on_accept","F3_nq_vs_es","F4_vol_ratio","F5_close_loc","F6_last30","F7_on_compress","F8_dow","F8_opex","F9_rates_corr","F10_on_poc_dist"]
ALPHA=0.05/50; out=[]
def ztest(p1,n1,p2,n2):
    p=(p1*n1+p2*n2)/(n1+n2); se=math.sqrt(p*(1-p)*(1/n1+1/n2)) if p*(1-p)>0 else np.nan
    return 2*(1-norm.cdf(abs(p1-p2)/se)) if se and se>0 else np.nan
half=pd.Timestamp("2019-01-01").date()
def buckets(s, f):
    if f in ("F2_on_accept","F8_dow","F8_opex"): return s.astype(int).astype(str)
    return pd.qcut(s.rank(method="first"), 5, labels=["Q1","Q2","Q3","Q4","Q5"]).astype(str)
def wf_auc(feats, target):
    d=df.dropna(subset=feats+[target]); aucs=[]; lls=[]; n=0
    for Y in range(2018, 2027):
        tr=d[pd.Series(d.index).apply(lambda x:x.year).values<Y]; te=d[pd.Series(d.index).apply(lambda x:x.year).values==Y]
        if len(te)<30 or tr[target].nunique()<2: continue
        m=LogisticRegression(C=1.0, max_iter=500); mu=tr[feats].mean(); sd=tr[feats].std().replace(0,1)
        m.fit(((tr[feats]-mu)/sd).values, tr[target].values); p=m.predict_proba(((te[feats]-mu)/sd).values)[:,1]
        if te[target].nunique()<2: continue
        aucs.append(roc_auc_score(te[target], p)*len(te)); lls.append(log_loss(te[target], p)*len(te)); n+=len(te)
    return (sum(aucs)/n if n else np.nan, sum(lls)/n if n else np.nan, n)
lines=[]; lines.append(f"# Open-factor sweep — results · {pd.Timestamp.now():%Y-%m-%d %H:%M}\n\nSessions: {len(df)} ({df.index[0]} → {df.index[-1]}), quad Fridays excluded, α = {ALPHA:.4f}.\n")
lines.append("## Base rates\n"+" · ".join(f"{t}: {df[t].mean()*100:.1f}%" for t in TARGETS)+"\n")
base_auc={t: wf_auc(BASE,t) for t in TARGETS}
lines.append("## Baseline walk-forward AUC (B1–B3, test years 2018–2026)\n"+" · ".join(f"{t}: {base_auc[t][0]:.3f} (n={base_auc[t][2]})" for t in TARGETS)+"\n")
lines.append("## Univariate + walk-forward, per factor × target\n\n| factor | target | Q-low P | Q-high P | Δpp | p | H1 sign | H2 sign | WF AUC base→+f | Δ | verdict |\n|---|---|---|---|---|---|---|---|---|---|---|")
verdicts={}
for f in CANDS:
    for t in TARGETS:
        d=df.dropna(subset=[f,t]); b=buckets(d[f],f); g=d.groupby(b)[t].agg(["mean","count"])
        keys=list(g.index); lo,hi=keys[0],keys[-1]
        if f=="F8_dow": lo,hi=g["mean"].idxmin(),g["mean"].idxmax()
        p=ztest(g.loc[hi,"mean"],g.loc[hi,"count"],g.loc[lo,"mean"],g.loc[lo,"count"]); dpp=(g.loc[hi,"mean"]-g.loc[lo,"mean"])*100
        def sgn(sub):
            gg=sub.groupby(buckets(sub[f],f))[t].mean(); 
            try: return np.sign(gg.loc[hi]-gg.loc[lo])
            except Exception: return np.nan
        s1=sgn(d[pd.Series(d.index).values<half]); s2=sgn(d[pd.Series(d.index).values>=half])
        a=wf_auc(BASE+[f],t); lift=a[0]-base_auc[t][0] if not np.isnan(a[0]) else np.nan
        ok = abs(dpp)>=6 and p<ALPHA and s1==s2 and not np.isnan(s1) and lift>=0.02
        v="VALIDATED" if ok else ("suggestive" if (abs(dpp)>=4 and p<0.01 and s1==s2) else "refuted")
        verdicts[(f,t)]=v
        lines.append(f"| {f} | {t} | {g.loc[lo,'mean']*100:.1f}% ({lo}) | {g.loc[hi,'mean']*100:.1f}% ({hi}) | {dpp:+.1f} | {p:.4f} | {s1:+.0f} | {s2:+.0f} | {base_auc[t][0]:.3f}→{a[0]:.3f} | {lift:+.3f} | {v} |")
surv=[f for f in CANDS if any(verdicts[(f,t)]=="VALIDATED" for t in TARGETS)]
lines.append("\n## Survivors (VALIDATED on ≥1 target)\n"+(", ".join(surv) if surv else "none")+"\n")
if surv:
    lines.append("## Baseline + all survivors, walk-forward\n"+" · ".join(f"{t}: {base_auc[t][0]:.3f}→{wf_auc(BASE+surv,t)[0]:.3f}" for t in TARGETS)+"\n")
lines.append("## Suggestive (fails the frozen rule, worth a forward shadow)\n"+"\n".join(f"- {f} × {t}" for (f,t),v in verdicts.items() if v=="suggestive")+"\n")
res="\n".join(lines); open(f"{FP}/studies/2026-09-08-open-factors-results.md","w").write(res); print(res)