fithia2/apps/intraday_bt/scripts/diag_orb_bb_pctb.py

"""
V37 Diagnostic: Bollinger Band %B and BB Width pre-breakout signal

Hypothesis A (compression): stocks with narrowing Bollinger Bands (low BB width) before a
gap-up may be in a coiling phase — explosive breakout, better follow-through.

Hypothesis B (position): stocks near the upper band (%B > 0.8) have confirmed momentum
and continue higher after the gap.

Features:
  bb_pct_b  : (close - lower_band) / (upper_band - lower_band), last prev_bar [0-1+]
  bb_width  : (upper_band - lower_band) / close — normalized band width (compression)
  bb_width_pct: percentile rank of bb_width vs own 60d history (0=tightest compression)

Source: V24 400d run JSON + daily parquet cache.
"""
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(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))
            closes.append(float(rows["close"][i] or 0))
    if not closes or closes[-1] <= 0:
        return None
    return {"date": date, "open": opens[0] if opens else 0, "close": closes[-1]}


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


def compute_bb(bars: list[dict], window: int = 20, num_std: float = 2.0) -> tuple[float, float] | None:
    """Returns (pct_b, width) from last `window` bars. bars sorted oldest→newest."""
    if len(bars) < window:
        return None
    tail = bars[-window:]
    closes = [b["close"] for b in tail]
    mean = sum(closes) / window
    variance = sum((c - mean) ** 2 for c in closes) / window
    std = variance ** 0.5
    if std == 0:
        return None
    upper = mean + num_std * std
    lower = mean - num_std * std
    last_close = closes[-1]
    band_width = upper - lower
    if band_width <= 0:
        return None
    pct_b = (last_close - lower) / band_width
    width = band_width / last_close  # normalized
    return pct_b, width


def compute_bb_width_percentile(bars: list[dict], window: int = 20, history: int = 60) -> float | None:
    """Percentile rank of current BB width vs own last `history` days."""
    if len(bars) < window + history:
        return None
    widths = []
    for i in range(history):
        end_idx = len(bars) - history + i + 1
        tail = bars[max(0, end_idx - window):end_idx]
        if len(tail) < window:
            continue
        closes = [b["close"] for b in tail]
        mean = sum(closes) / len(closes)
        std = (sum((c - mean) ** 2 for c in closes) / len(closes)) ** 0.5
        if std == 0:
            continue
        widths.append((mean + 2 * std - (mean - 2 * std)) / closes[-1])
    if len(widths) < 20:
        return None
    current = widths[-1]
    rank = sum(1 for w in widths if w <= current) / len(widths)
    return rank


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 n={n}")
        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})")
    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']}")
    print(f"    → {'ALL GATES PASS ✓' if (g1 and g2 and g3 and g5a) else 'FAIL'}")


def main() -> None:
    print("=== V37 BB %B and Width 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}%  DD: {m.get('max_drawdown_pct', 0)*100:.2f}%  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)}")

    tickers_needed = sorted(set(r["ticker"] for r in trade_records))
    min_date = min(r["date"] for r in trade_records)
    max_date = max(r["date"] for r in trade_records)
    start_cal = (dt.date.fromisoformat(min_date) - dt.timedelta(days=150)).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)

    print(f"Loading bars for {len(tickers_needed)} tickers ({start_cal} → {max_date})...")
    ticker_bars: dict[str, list[dict]] = {}
    with concurrent.futures.ThreadPoolExecutor(max_workers=8) as ex:
        def _load(ticker: str) -> tuple[str, list[dict]]:
            return ticker, load_ticker_bars(ticker, all_dates)
        for ticker, bars in ex.map(_load, tickers_needed):
            if bars:
                ticker_bars[ticker] = bars
    print(f"Loaded: {len(ticker_bars)} / {len(tickers_needed)} tickers\n")

    pctb_vals, width_vals, width_pct_vals, obv20_vals, r_mults = [], [], [], [], []
    missing_bb, missing_wp = 0, 0

    from libs.intraday.features import compute_obv_slope_approx
    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]

        bb = compute_bb(prev_bars, window=20) if len(prev_bars) >= 20 else None
        if bb is None:
            missing_bb += 1
            continue

        pct_b, width = bb
        wp = compute_bb_width_percentile(prev_bars, window=20, history=60)
        if wp is None:
            missing_wp += 1

        obv20 = compute_obv_slope_approx(prev_bars, lookback=20) if len(prev_bars) >= 22 else None

        pctb_vals.append(pct_b)
        width_vals.append(width)
        if wp is not None:
            width_pct_vals.append(wp)
        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 (BB computable): {n_valid} / {len(trade_records)}")
    print(f"Missing BB: {missing_bb}  Missing width_pct: {missing_wp}  Valid width_pct: {len(width_pct_vals)}")

    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)")

    obv20_aligned_pctb = obv20_vals[:len(pctb_vals)]
    report_feature("bb_pct_b", pctb_vals, r_mults, obv20_aligned_pctb)
    report_feature("bb_width (compression)", width_vals, r_mults, obv20_aligned_pctb)
    if len(width_pct_vals) >= 120:
        report_feature("bb_width_percentile_60d", width_pct_vals, r_mults[:len(width_pct_vals)], obv20_aligned_pctb[:len(width_pct_vals)])
    else:
        print(f"\n  bb_width_percentile_60d: insufficient n={len(width_pct_vals)} (need 120)")

    p_pctb_width = pearson(pctb_vals, width_vals)
    print(f"\n  Inter-feature: ρ(bb_pct_b, bb_width) = {p_pctb_width:.3f}")
    mean_pctb = sum(pctb_vals)/len(pctb_vals) if pctb_vals else 0
    print(f"  %B distribution: mean={mean_pctb:.2f}  (0=lower, 0.5=mid, 1=upper, >1=above band)")

    print("\n=== Summary ===")
    print(f"G2 target: ≥ 0.30R. OBV-slope = 0.394R (only axis that cleared).")


if __name__ == "__main__":
    main()