#!/usr/bin/env python3
"""
vzd_reproduce.py -- reproduce the Volatility-Normalized Zone Displacement (VZD) result.

VZD = |price - AOI| / ATR  : how far price has displaced from the zone it is
reacting to, in units of current volatility.

This script takes a list of trades, each with its VZD at signal time and its
realized R-multiple, and shows that expectancy rises monotonically with VZD --
and that the relationship survives a strict out-of-sample (chronological) split.

NO third-party libraries required. Just Python 3.

USAGE
-----
    python vzd_reproduce.py                 # uses the bundled vzd_sample.csv
    python vzd_reproduce.py mytrades.csv    # uses your own data

INPUT CSV FORMAT
----------------
A header row, then one row per FILLED trade, with (at minimum) columns:
    vzd          float   the trade's VZD at signal time
    realized_r   float   the trade's realized R-multiple (profit in units of risk)
Optional:
    seq          int     a monotonic sequence/time key; if present, rows are
                         sorted by it so the out-of-sample split is chronological.

If you have raw trades instead, compute VZD yourself first:
    vzd = abs(price_at_signal - aoi_level) / atr
"""
import csv
import sys


def load(path):
    rows = []
    with open(path, newline="") as f:
        for r in csv.DictReader(f):
            try:
                v = float(r["vzd"]); rr = float(r["realized_r"])
            except (KeyError, ValueError):
                continue
            seq = float(r["seq"]) if r.get("seq") not in (None, "") else len(rows)
            rows.append((seq, v, rr))
    rows.sort(key=lambda x: x[0])           # chronological for the OOS split
    return [(v, rr) for _, v, rr in rows]


def quantiles(xs, ps):
    s = sorted(xs)
    out = []
    for p in ps:
        i = p * (len(s) - 1)
        lo = int(i)
        out.append(s[lo] if lo + 1 >= len(s) else s[lo] + (i - lo) * (s[lo + 1] - s[lo]))
    return out


def mean(xs):
    return sum(xs) / len(xs) if xs else 0.0


def ranks(xs):
    # average-rank for ties
    order = sorted(range(len(xs)), key=lambda i: xs[i])
    rk = [0.0] * len(xs)
    i = 0
    while i < len(xs):
        j = i
        while j + 1 < len(xs) and xs[order[j + 1]] == xs[order[i]]:
            j += 1
        avg = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            rk[order[k]] = avg
        i = j + 1
    return rk


def spearman(a, b):
    ra, rb = ranks(a), ranks(b)
    ma, mb = mean(ra), mean(rb)
    num = sum((x - ma) * (y - mb) for x, y in zip(ra, rb))
    da = sum((x - ma) ** 2 for x in ra) ** 0.5
    db = sum((y - mb) ** 2 for y in rb) ** 0.5
    return num / (da * db) if da and db else 0.0


def profit_factor(rs):
    gp = sum(r for r in rs if r > 0)
    gl = -sum(r for r in rs if r <= 0)
    return gp / gl if gl else float("inf")


def main():
    path = sys.argv[1] if len(sys.argv) > 1 else "vzd_sample.csv"
    data = load(path)
    if not data:
        print(f"No usable rows in {path} (need 'vzd' and 'realized_r' columns)."); return
    V = [v for v, _ in data]; R = [r for _, r in data]
    print(f"\nLoaded {len(data):,} filled trades from {path}")
    print(f"Base expectancy: {mean(R):+.4f} R/trade   profit factor: {profit_factor(R):.2f}\n")

    print("=== Expectancy by VZD quartile (full sample) ===")
    qs = quantiles(V, [0, .25, .5, .75, 1.0])
    for i in range(4):
        lo, hi = qs[i], qs[i + 1]
        bucket = [r for v, r in data if (v >= lo and (v <= hi if i == 3 else v < hi))]
        bar = "#" * max(0, int(mean(bucket) * 100))
        print(f"  Q{i+1}  VZD [{lo:5.2f}, {hi:5.2f}]   n={len(bucket):5d}   "
              f"avgR={mean(bucket):+.3f}  {bar}")

    cut = int(len(data) * 0.70)
    tr, te = data[:cut], data[cut:]
    print("\n=== Out-of-sample check (train = first 70%, test = last 30%) ===")
    print(f"  train Spearman(VZD, R) = {spearman([v for v,_ in tr], [r for _,r in tr]):+.4f}")
    print(f"  TEST  Spearman(VZD, R) = {spearman([v for v,_ in te], [r for _,r in te]):+.4f}")
    te_sorted = sorted(te, key=lambda x: x[0])
    d = max(1, len(te_sorted) // 10)
    d1 = mean([r for _, r in te_sorted[:d]]); d10 = mean([r for _, r in te_sorted[-d:]])
    print(f"  TEST  top-decile minus bottom-decile R = {d10 - d1:+.4f}R "
          f"(low-VZD decile {d1:+.3f} -> high-VZD decile {d10:+.3f})")

    print("\n=== Gate it: keep trades with VZD >= threshold (full sample) ===")
    for thr in (0.0, 0.20, 0.30, 0.45, 0.60):
        sub = [r for v, r in data if v >= thr]
        if sub:
            print(f"  VZD >= {thr:.2f}   n={len(sub):5d}   avgR={mean(sub):+.3f}   "
                  f"PF={profit_factor(sub):.2f}")
    print()


if __name__ == "__main__":
    main()