ForexBot2/backtest/research/burst_precursors.py

"""Are there signs BEFORE a violent bar?  Analysis tooling only — never imported by the live code.

For every bar t, features are computed from bars <= t-1 only. Targets:
  next_bar  : bar t itself is a burst (range >= K x prior Wilder ATR14)
  next_hour : any burst in bars t .. t+H-1
Reports per-quintile burst rates (lift vs base rate) in both halves, the same
within/outside the scheduled-data hours, and an out-of-sample logistic model
(fit on half 1, scored on half 2): AUC + lift in the top 5% of scores.

    python backtest/research/burst_precursors.py PARQUET BAR_MINUTES [--k 3] [--start ... --end ...]
"""
import argparse, sys
import numpy as np, pandas as pd


def wilder_atr(df, period=14):
    pc = df["Close"].shift(1)
    tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
    return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()


def build(df, bar_minutes, k):
    d = df.copy()
    rng = d["High"] - d["Low"]
    atr_prev = wilder_atr(d).shift(1)                       # ATR known at t-1
    d["burst"] = (rng >= k * atr_prev).astype(int)
    H = max(1, 60 // bar_minutes)
    # burst within next H bars (t..t+H-1)
    d["next_hour"] = d["burst"][::-1].rolling(H, min_periods=1).max()[::-1].astype(int)
    d["next_bar"] = d["burst"]
    close_t = d.index + pd.Timedelta(minutes=bar_minutes)
    d["hour"] = close_t.hour
    r1 = rng.shift(1)
    d["squeeze6"] = r1.rolling(6).mean() / atr_prev
    d["squeeze12"] = r1.rolling(12).mean() / atr_prev
    d["last_bar_range"] = r1 / atr_prev
    sma = d["Close"].shift(1).rolling(20).mean(); sd = d["Close"].shift(1).rolling(20).std(ddof=0)
    bbw = (4 * sd / sma)
    d["bbw_pct500"] = bbw.rolling(500).rank(pct=True)
    v1 = d["Volume"].shift(1)
    d["vol_ratio"] = v1.rolling(3).mean() / v1.rolling(50).mean()
    d["drift12"] = (d["Close"].shift(1) - d["Close"].shift(13)).abs() / atr_prev
    d["atr_pct500"] = atr_prev.rolling(500).rank(pct=True)
    b = d["burst"].shift(1).fillna(0).astype(int)
    since = np.zeros(len(d)); cnt = 200
    bv = b.values
    for i in range(len(d)):
        cnt = 0 if bv[i] == 1 else min(200, cnt + 1)
        since[i] = cnt
    d["bars_since_burst"] = since
    return d.dropna(subset=["squeeze6", "squeeze12", "bbw_pct500", "vol_ratio", "drift12", "atr_pct500"])


FEATS = ["squeeze6", "squeeze12", "last_bar_range", "bbw_pct500", "vol_ratio", "drift12", "atr_pct500", "bars_since_burst"]


def quintile_table(d, feat, target, half):
    q = pd.qcut(d[feat], 5, labels=False, duplicates="drop")
    base1 = d.loc[half == 1, target].mean(); base2 = d.loc[half == 2, target].mean()
    rows = []
    for b in sorted(q.dropna().unique()):
        m = q == b
        lo, hi = d.loc[m, feat].min(), d.loc[m, feat].max()
        r1 = d.loc[m & (half == 1), target].mean(); r2 = d.loc[m & (half == 2), target].mean()
        rows.append((int(b), lo, hi, int(m.sum()), 100 * r1, r1 / base1, 100 * r2, r2 / base2))
    t = pd.DataFrame(rows, columns=["bin", "lo", "hi", "n", "h1_rate%", "h1_lift", "h2_rate%", "h2_lift"])
    return t.round(3)


def main():
    ap = argparse.ArgumentParser(); ap.add_argument("parquet"); ap.add_argument("bar_minutes", type=int)
    ap.add_argument("--k", type=float, default=3.0); ap.add_argument("--start"); ap.add_argument("--end")
    a = ap.parse_args()
    df = pd.read_parquet(a.parquet)
    if a.start: df = df[df.index >= pd.Timestamp(a.start, tz="UTC")]
    if a.end: df = df[df.index < pd.Timestamp(a.end, tz="UTC")]
    d = build(df, a.bar_minutes, a.k)
    half = pd.Series(np.where(np.arange(len(d)) < len(d) // 2, 1, 2), index=d.index)
    print(f"bars {len(d)}  K={a.k}  base rate next_bar {100*d.next_bar.mean():.3f}%  next_hour {100*d.next_hour.mean():.2f}%")
    print(f"half 1 ends {d.index[len(d)//2 - 1]}")
    # 1. the clock
    hr = d.groupby("hour")[["next_bar", "next_hour"]].mean() * 100
    hr["lift_next_bar"] = hr["next_bar"] / (100 * d.next_bar.mean())
    print("\n== THE CLOCK (bar-close hour UTC): burst rate % and lift ==")
    print(hr.round(2).T.to_string())
    sched = d["hour"].isin([1, 12, 13])
    print(f"\nscheduled hours (01,12,13) hold {100*d.loc[sched,'next_bar'].sum()/d.next_bar.sum():.1f}% of bursts in {100*sched.mean():.1f}% of bars")
    # 2. per-feature quintiles, both halves
    for target in ("next_bar", "next_hour"):
        print(f"\n== TARGET {target}: quintile burst rate (lift vs base) — ALL hours ==")
        for f in FEATS:
            t = quintile_table(d, f, target, half)
            print(f"-- {f}: lifts h1 {list(t.h1_lift)}  h2 {list(t.h2_lift)}  (bins lo..hi {t.lo.iloc[0]:.2f}..{t.hi.iloc[-1]:.2f})")
    # 3. same, OUTSIDE the scheduled hours (does anything predict the unscheduled ones?)
    u = d[~sched]; hu = half[~sched]
    print(f"\n== next_bar lifts OUTSIDE scheduled hours (base {100*u.next_bar.mean():.3f}%) ==")
    for f in FEATS:
        t = quintile_table(u, f, "next_bar", hu)
        print(f"-- {f}: h1 {list(t.h1_lift)}  h2 {list(t.h2_lift)}")
    # 4. OOS logistic model: fit half 1, score half 2
    try:
        from sklearn.linear_model import LogisticRegression
        from sklearn.metrics import roc_auc_score
        X = d[FEATS].copy()
        for h in range(24): X[f"h{h}"] = (d["hour"] == h).astype(float)
        Xc = d[FEATS].copy()   # clock-free variant
        for name, XX in (("features+clock", X), ("features only", Xc), ("clock only", X[[f"h{h}" for h in range(24)]])):
            mu, sd = XX[half == 1].mean(), XX[half == 1].std().replace(0, 1)
            Z = (XX - mu) / sd
            for target in ("next_bar", "next_hour"):
                m = LogisticRegression(C=0.1, max_iter=2000, class_weight="balanced").fit(Z[half == 1], d.loc[half == 1, target])
                p = m.predict_proba(Z[half == 2])[:, 1]; y = d.loc[half == 2, target]
                auc = roc_auc_score(y, p)
                top = p >= np.quantile(p, 0.95)
                print(f"OOS {name:15s} {target:9s}: AUC {auc:.3f} | top-5% of scores: burst rate {100*y[top].mean():.2f}% vs base {100*y.mean():.2f}% (lift {y[top].mean()/y.mean():.2f}x), captures {100*y[top].sum()/y.sum():.1f}% of bursts")
    except ImportError:
        print("sklearn missing — OOS model skipped")
    # 5. direction: does the burst continue the prior 12-bar drift?
    bt = d[d.burst == 1]
    prior = np.sign(bt["Close"].shift(1) - bt["Close"].shift(13)) if False else np.sign(d["Close"].shift(1) - d["Close"].shift(13))[d.burst == 1]
    bdir = np.sign(bt["Close"] - bt["Open"])
    agree = (prior == bdir)[(prior != 0) & (bdir != 0)]
    print(f"\nburst direction == prior-hour drift direction: {100*agree.mean():.1f}% of {len(agree)} bursts (50% = no information)")


if __name__ == "__main__":
    main()