fithia2/apps/intraday_bt/scripts/diag_orb_rsi14.py

"""
V36 Diagnostic: RSI-14 Pre-Breakout Signal

Hypothesis: RSI-14 in the days before an ORB gap-up predicts follow-through quality.
Two competing sub-hypotheses:
  (A) Momentum: high RSI (>60) = confirmed uptrend, clean breakout. Positive Pearson.
  (B) Mean-reversion: high RSI = overbought, gap sells off. Negative Pearson.
  (C) Sweet spot: moderate RSI (40-60) = breakout from accumulation, not extended.

Features:
  rsi_14     : raw RSI-14 [0, 100] on last prev_bar
  rsi_momentum: (RSI - 50) / 50, centered at midline [-1, +1]

Source: V24 400d run JSON, daily bars from parquet cache.
Uses Cutler's RSI (simple average, same formula as libs/features/market_features.py).
"""
from __future__ import annotations

import concurrent.futures
import datetime as dt
import json
import os
import sys
from pathlib import Path

import pyarrow.parquet as pq

sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../../..")))

from zoneinfo import ZoneInfo

_ET = ZoneInfo("America/New_York")
_MKT_OPEN = dt.time(9, 30)
_MKT_CLOSE = dt.time(16, 0)

INTRADAY_CACHE_DIR = "data/cache/intraday"
V24_400D_RUN = "runs/intraday_orb/intraday_20260422_011012_06f59ede.json"


def _parse_ts(ts_raw: object) -> dt.datetime:
    s = str(ts_raw)
    if s.endswith("Z"):
        s = s[:-1] + "+00:00"
    return dt.datetime.fromisoformat(s).astimezone(_ET)


def _build_daily_bar_from_parquet(path: Path, date: str) -> dict | None:
    try:
        table = pq.read_table(str(path))
        rows = table.to_pydict()
    except Exception:
        return None
    opens, highs, lows, closes, vols = [], [], [], [], []
    for i, ts_raw in enumerate(rows.get("timestamp", [])):
        try:
            ts = _parse_ts(ts_raw)
        except Exception:
            continue
        if _MKT_OPEN <= ts.time() < _MKT_CLOSE:
            opens.append(float(rows["open"][i] or 0))
            highs.append(float(rows["high"][i] or 0))
            lows.append(float(rows["low"][i] or 0))
            closes.append(float(rows["close"][i] or 0))
            vols.append(float(rows["volume"][i] or 0))
    if not opens or closes[-1] <= 0:
        return None
    return {
        "date": date, "open": opens[0], "high": max(highs),
        "low": min(lows), "close": closes[-1], "volume": sum(vols),
    }


def load_ticker_daily_bars(ticker: str, needed_dates: list[str]) -> list[dict]:
    root = Path(INTRADAY_CACHE_DIR) / ticker
    if not root.is_dir():
        return []
    bars = []
    for date in needed_dates:
        p = root / f"{date}.parquet"
        if not p.exists():
            continue
        bar = _build_daily_bar_from_parquet(p, date)
        if bar:
            bars.append(bar)
    return sorted(bars, key=lambda b: b["date"])


def compute_rsi_14(bars: list[dict], period: int = 14) -> float | None:
    """Cutler's RSI from daily close prices. bars sorted oldest→newest."""
    if len(bars) < period + 2:
        return None
    tail = bars[-(period + 1):]
    gains, losses = [], []
    for i in range(period):
        change = tail[i + 1]["close"] - tail[i]["close"]
        if change >= 0:
            gains.append(change)
            losses.append(0.0)
        else:
            gains.append(0.0)
            losses.append(abs(change))
    avg_gain = sum(gains) / period
    avg_loss = sum(losses) / period
    if avg_loss == 0:
        return 100.0
    return 100.0 - (100.0 / (1.0 + avg_gain / avg_loss))


def pearson(xs: list[float], ys: list[float]) -> float | None:
    n = len(xs)
    if n < 2:
        return None
    xm, ym = sum(xs) / n, sum(ys) / n
    num = sum((xs[i] - xm) * (ys[i] - ym) for i in range(n))
    dx = sum((x - xm) ** 2 for x in xs) ** 0.5
    dy = sum((y - ym) ** 2 for y in ys) ** 0.5
    if dx <= 0 or dy <= 0:
        return None
    return num / (dx * dy)


def tercile_stats(vals: list[float], rs: list[float]) -> dict:
    if len(vals) < 9:
        return {}
    pairs = sorted(zip(vals, rs), key=lambda p: p[0])
    n = len(pairs)
    t = n // 3
    def stats(sub):
        ys = [p[1] for p in sub]
        return {"n": len(ys), "wr": sum(1 for y in ys if y > 0) / len(ys), "avg_r": sum(ys) / len(ys)}
    return {"low": stats(pairs[:t]), "mid": stats(pairs[t:2*t]), "high": stats(pairs[2*t:])}


def report_feature(label: str, vals: list[float], rs: list[float], obv20: list[float] | None = None) -> None:
    n = len(vals)
    p = pearson(vals, rs)
    ts = tercile_stats(vals, rs)
    rho_obv20 = pearson(vals, obv20) if obv20 else None
    if not ts or p is None:
        print(f"  {label}: insufficient data")
        return
    low, mid, high = ts["low"], ts["mid"], ts["high"]
    avg_r_gap = abs(high["avg_r"] - low["avg_r"])
    wr_gap = abs(high["wr"] - low["wr"])
    g1 = abs(p) >= 0.07 and n >= 120
    g2 = avg_r_gap >= 0.30
    g3 = wr_gap >= 0.05
    g5a = rho_obv20 is None or abs(rho_obv20) < 0.70
    print(f"\n  [{label}]  n={n}  Pearson={p:+.3f}")
    print(f"    G1: {'PASS' if g1 else 'FAIL'} (|{abs(p):.3f}| {'≥' if abs(p)>=0.07 else '<'} 0.07, n={n} {'≥' if n>=120 else '<'} 120)")
    print(f"    G2: {'PASS' if g2 else 'FAIL'} (avg_R gap = {avg_r_gap:.3f}R  [threshold 0.30R])")
    print(f"    G3: {'PASS' if g3 else 'FAIL'} (WR gap = {wr_gap*100:.1f}pp  [threshold 5pp])")
    if rho_obv20 is not None:
        print(f"    G5a: {'PASS' if g5a else 'FAIL'} (|ρ(feature, obv_slope_20)| = {abs(rho_obv20):.3f}  [threshold 0.70])")
    print(f"    Bottom tercile: WR {low['wr']*100:.1f}%  avg_R {low['avg_r']:+.3f}  n={low['n']}")
    print(f"    Middle tercile: WR {mid['wr']*100:.1f}%  avg_R {mid['avg_r']:+.3f}  n={mid['n']}")
    print(f"    Top tercile:    WR {high['wr']*100:.1f}%  avg_R {high['avg_r']:+.3f}  n={high['n']}")
    all_pass = g1 and g2 and g3 and g5a
    print(f"    → {'ALL GATES PASS ✓' if all_pass else 'FAIL'}")


def main() -> None:
    print("=== V36 RSI-14 Pre-Breakout Diagnostic ===\n")

    with open(V24_400D_RUN) as f:
        run_data = json.load(f)
    trades = run_data.get("trades", [])
    m = run_data.get("metrics", {})
    print(f"Loaded V24 400d run: {m.get('total_trades')} trades, "
          f"{m.get('start_date')} → {m.get('end_date')}")
    print(f"Return: {m.get('total_return_pct', 0)*100:.2f}%  "
          f"DD: {m.get('max_drawdown_pct', 0)*100:.2f}%  "
          f"Sharpe: {m.get('sharpe_ratio', 0):.3f}\n")

    trade_records = [
        {"ticker": t["ticker"], "date": t["date"][:10],
         "r_multiple": float(t["r_multiple_at_exit"])}
        for t in trades
        if t.get("r_multiple_at_exit") is not None
    ]
    print(f"Trades with r_multiple: {len(trade_records)}")

    # Gather all unique tickers and dates needed (extended window for RSI lookback)
    tickers_needed = sorted(set(r["ticker"] for r in trade_records))
    dates_by_ticker: dict[str, set[str]] = {t: set() for t in tickers_needed}
    for rec in trade_records:
        dates_by_ticker[rec["ticker"]].add(rec["date"])

    print(f"Loading daily bars for {len(tickers_needed)} tickers...")

    # Determine needed calendar range per ticker (need 30+ prior trading days)
    import datetime as dt
    min_date = min(r["date"] for r in trade_records)
    max_date = max(r["date"] for r in trade_records)
    # Build full calendar range with buffer
    start_cal = (dt.date.fromisoformat(min_date) - dt.timedelta(days=90)).isoformat()
    all_dates = []
    d = dt.date.fromisoformat(start_cal)
    end_d = dt.date.fromisoformat(max_date)
    while d <= end_d:
        all_dates.append(d.isoformat())
        d += dt.timedelta(days=1)

    def _load(ticker: str) -> tuple[str, list[dict]]:
        return ticker, load_ticker_daily_bars(ticker, all_dates)

    ticker_bars: dict[str, list[dict]] = {}
    with concurrent.futures.ThreadPoolExecutor(max_workers=8) as ex:
        for ticker, bars in ex.map(_load, tickers_needed):
            if bars:
                ticker_bars[ticker] = bars

    print(f"Loaded bars for {len(ticker_bars)} / {len(tickers_needed)} tickers\n")

    # Compute features per trade
    rsi_vals, rsi_mom_vals, obv20_vals, r_mults = [], [], [], []
    missing = 0

    for rec in trade_records:
        ticker = rec["ticker"]
        date = rec["date"]
        r = rec["r_multiple"]

        bars = ticker_bars.get(ticker, [])
        prev_bars = [b for b in bars if b["date"] < date]

        rsi = compute_rsi_14(prev_bars, period=14) if len(prev_bars) >= 16 else None
        if rsi is None:
            missing += 1
            continue

        # OBV slope 20d for redundancy check
        from libs.intraday.features import compute_obv_slope_approx
        obv20 = compute_obv_slope_approx(prev_bars, lookback=20) if len(prev_bars) >= 22 else None

        rsi_vals.append(rsi)
        rsi_mom_vals.append((rsi - 50.0) / 50.0)
        obv20_vals.append(obv20 if obv20 is not None else 0.0)
        r_mults.append(r)

    n_valid = len(r_mults)
    print(f"Valid trades (RSI computable): {n_valid} / {len(trade_records)}  (missing: {missing})")

    print("\n" + "=" * 60)
    print("GATE RESULTS (G1: |P|≥0.07 & n≥120; G2: avg_R≥0.30R; G3: WR≥5pp; G5a: ρ<0.70 vs obv_slope_20)")

    report_feature("rsi_14", rsi_vals, r_mults, obv20_vals)
    report_feature("rsi_momentum (RSI-50)/50", rsi_mom_vals, r_mults, obv20_vals)

    # Distribution stats
    mean_rsi = sum(rsi_vals) / len(rsi_vals) if rsi_vals else 0
    print(f"\n  RSI distribution: mean={mean_rsi:.1f}  min={min(rsi_vals):.1f}  max={max(rsi_vals):.1f}")

    print("\n=== Summary ===")
    p = pearson(rsi_vals, r_mults)
    if p is not None:
        direction = "Momentum (high RSI = better)" if p > 0 else "Mean-reversion (low RSI = better)"
        print(f"Signal direction: {direction}  (Pearson={p:+.3f})")
    print("Context: V24's only passing axis: OBV-slope G2=0.394R. Target: G2 ≥ 0.30R.")


if __name__ == "__main__":
    main()