Warrior_EA/research/test_wyckoff.py

"""Does volume-at-price structure do anything measurable? (Villahermosa, "Wyckoff 2.0")

This is a different family from everything tested so far in this project. Every previous
test asked a DIRECTIONAL question - can something predict the sign of the next move - and
the answer was consistently "a little, and less than the spread". The volume-profile claims
are not directional. They are claims about the PATH:

    HVN / VPOC   agreement -> a magnet. Price is drawn to it and lingers.  ("good targets")
    LVN          rejection -> price refuses to trade there and turns.      ("good stops")
    Value area   the 80% rule: re-entry and acceptance implies traversal.

That distinction is why this is worth testing after six negative families. A path claim can
be true while every directional claim is false, and if HVN/LVN really are non-uniform then
STOP AND TARGET PLACEMENT carries edge with no forecast at all - which is exactly what was
asked for: where to enter, where to exit, where to put the stop.

THE NULL - and why the obvious one is degenerate
------------------------------------------------
The tempting test is "how often does price reach the naked VPOC". That number is
meaningless on its own, and worse, it cannot be fixed by comparing against a control level
at the SAME distance, because at the same distance the control IS the same price. The
probability of touching a level is a function of its distance and nothing else.

So the null has to destroy the level's IDENTITY while preserving its GEOMETRY:

    null offset  =  ATR_i  x  (d / ATR)_permuted

Distances in volatility units are shuffled across events. The marginal distribution of
distance is preserved exactly, the local volatility scaling is preserved exactly (a naked
VPOC in a quiet market stays near, in a wild one stays far), the market conditions are the
real ones - and the only thing removed is that the level sits where the market previously
agreed on value. Anything the real levels do beyond the permuted ones is attributable to
volume structure.

The 80% rule gets a better null still: an EXACT martingale benchmark. Entering a band at
price p between a near edge and a far edge, a driftless process reaches the far edge first
with probability (p - near)/(far - near). No simulation needed, and it correctly punishes
the fact that acceptance usually happens close to the edge you came in through.
"""
import numpy as np, sys, os, datetime as dt
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
from vplevels import (load, day_index, session_levels, composite_nodes, naked_vpocs,
                      broker_day)

BARS = 'c:/Users/admin/Documents/Workspaces/Market Data/bars/'
SYMS = ('EURUSD', 'USDJPY', 'XAUUSD', 'SP500')


def load_bars(sym, tf='M5'):
    z = np.load(f"{BARS}{sym}_{tf}_ticks.npz", allow_pickle=True)
    a = z['bars']
    return a, {str(c): k for k, c in enumerate(z['columns'])}


def atr_of(h, l, c, n=14):
    pc = np.roll(c, 1); pc[0] = c[0]
    tr = np.maximum(h - l, np.maximum(np.abs(h - pc), np.abs(l - pc)))
    out = np.convolve(tr, np.ones(n) / n, mode='full')[:len(tr)]
    out[:n] = tr[:n].mean()
    return out


def forward_cummax(h, l, idx, H):
    """Running max-high and min-low over bars idx+1 .. idx+k, for every k <= H.

    This is the whole permutation test in one array. 'Does price reach a level at distance
    d' is exactly 'd <= the forward excursion', so once the excursion is known every
    permuted level is a single comparison instead of another 288-pass scan - the difference
    between 5 hours and a few seconds. Monotone in k, so the first-touch BAR is a
    searchsorted rather than a scan too.
    """
    W = np.lib.stride_tricks.sliding_window_view
    keep = idx + H < len(h)
    idx = idx[keep]
    cu = np.maximum.accumulate(W(h, H)[idx + 1], axis=1).astype(np.float32)
    cl = np.minimum.accumulate(W(l, H)[idx + 1], axis=1).astype(np.float32)
    return idx, keep, cu, cl


def touched(cu, cl, up, price):
    """Was `price` reached within the window? Pure comparison against the final excursion."""
    return np.where(up, price <= cu[:, -1], price >= cl[:, -1])


def touch_bar(cu, cl, up, price, idx):
    """Bar index of first touch, or -1.

    argmax over the boolean 'reached by bar k' matrix, not a per-row searchsorted: rows
    differ in their target so searchsorted cannot be batched, and a Python loop over
    100k rows inside a permutation loop is the whole runtime.
    """
    reach = np.where(up[:, None], cu >= price[:, None], cl <= price[:, None])
    k = reach.argmax(axis=1)
    return np.where(reach[np.arange(len(price)), k], idx + 1 + k, -1)


def first_touch(h, l, idx, price, up, H):
    """First bar within (idx, idx+H] whose range covers `price`. -1 if never.

    Vectorised over events by stepping the OFFSET, not the event: every event shares the
    same k = 1..H, so this is H vector passes rather than one Python loop per event. On
    100k events x H=288 that is ~30M element ops, about a second.
    """
    n = len(h)
    out = np.full(len(idx), -1, np.int64)
    live = np.ones(len(idx), bool)
    for k in range(1, H + 1):
        j = idx + k
        ok = live & (j < n)
        if not ok.any():
            break
        jj = j[ok]
        hit = np.where(up[ok], h[jj] >= price[ok], l[jj] <= price[ok])
        w = np.nonzero(ok)[0][hit]
        out[w] = j[w]
        live[w] = False
    return out


def race(h, l, idx, lo_px, hi_px, H):
    """Which of two barriers is touched first: +1 upper, -1 lower, 0 neither in H bars.
    A bar spanning both counts as the LOWER first - the pessimistic convention used
    everywhere else in this project, so results stay comparable."""
    n = len(h)
    out = np.zeros(len(idx), np.int8)
    live = np.ones(len(idx), bool)
    for k in range(1, H + 1):
        j = idx + k
        ok = live & (j < n)
        if not ok.any():
            break
        jj = j[ok]
        dn = l[jj] <= lo_px[ok]
        up = h[jj] >= hi_px[ok]
        res = np.where(dn, -1, np.where(up, 1, 0)).astype(np.int8)
        w = np.nonzero(ok)[0]
        got = res != 0
        out[w[got]] = res[got]
        live[w[got]] = False
    return out


def build_levels(sym, lookback=20):
    day, b, tk, dw, binsize = load(sym)
    days, cells = day_index(day, b, tk)
    sess = session_levels(days, cells)
    nodes = composite_nodes(days, cells, lookback=lookback)
    naked = naked_vpocs(sess)
    return sess, nodes, naked, binsize


def nearest(levels, price_bin, above):
    """Nearest level strictly above / below a price bin, or None."""
    if levels is None or len(levels) == 0:
        return None
    d = levels - price_bin
    m = d > 0 if above else d < 0
    if not m.any():
        return None
    return int(levels[m][np.argmin(np.abs(d[m]))])


#----------------------------------------------------------------------------------------
def collect(sym, kind, H=288, every=12, lookback=20, seed=7):
    """Reference bars -> (offset in ATR, real level price, market state) for one level type.

    `every` subsamples the M5 grid (12 = hourly) so that consecutive reference bars are not
    near-duplicates of each other; overlapping events inflate significance without adding
    information.
    """
    a, I = load_bars(sym)
    g = lambda c: a[:, I[c]]
    t, o, h, l, c = g('time'), g('open'), g('high'), g('low'), g('close')
    atr = atr_of(h, l, c, 14)
    atr = np.concatenate([[atr[0]], atr[:-1]])            # causal
    bd = broker_day(t.astype(np.int64))
    sess, nodes, naked, binsize = build_levels(sym, lookback)
    sday = {int(d): i for i, d in enumerate(sess[:, 0])}

    idx, offs, ups = [], [], []
    n = len(c)
    step = np.arange(14 + 1, n - H - 1, every)
    for i in step:
        d = int(bd[i])
        si = sday.get(d)
        if si is None or si < 1:
            continue
        pb = c[i] / binsize
        if kind == 'naked':
            lv = naked.get(d)
        elif kind in ('hvn', 'lvn'):
            nd = nodes.get(d)
            if nd is None:
                continue
            lv = nd[0] if kind == 'hvn' else nd[1]
        elif kind == 'vpoc':
            lv = np.array([sess[si - 1, 1]])
        elif kind == 'vaedge':
            lv = np.array([sess[si - 1, 2], sess[si - 1, 3]])
        else:
            raise ValueError(kind)
        if lv is None or len(lv) == 0:
            continue
        for up in (True, False):
            v = nearest(lv, pb, up)
            if v is None:
                continue
            off = v * binsize - c[i]
            if atr[i] <= 0:
                continue
            r = abs(off) / atr[i]
            if not (0.2 <= r <= 8.0):        # a level 30 ATR away is not an operational level
                continue
            idx.append(i); offs.append(off); ups.append(up)
    return (np.array(idx, np.int64), np.array(offs), np.array(ups, bool),
            atr, h, l, c, len(c))


def permute_within(r, up, rng):
    """Shuffle distances WITHIN each direction.

    Pooling the two would let an 'up' event inherit a 'down' event's distance. Up and down
    distances have different distributions (drift, and asymmetric ranges), and up and down
    touch rates differ, so a pooled shuffle mixes a direction effect into what is supposed
    to be a pure level-identity test.
    """
    out = np.empty_like(r)
    for m in (up, ~up):
        if m.any():
            out[m] = rng.permutation(r[m])
    return out


def magnet_test(sym, kind, H=288, every=12, nperm=2000, seed=7, quiet=False, pre=None):
    """Does price reach a real level more often than a distance-matched placebo?"""
    idx, off, up, atr, h, l, c, n = pre if pre else collect(sym, kind, H, every)
    if len(idx) < 500:
        print(f"  {sym} {kind}: only {len(idx)} events, skipping")
        return None
    idx2, keep, cu, cl = forward_cummax(h, l, idx, H)
    off, up = off[keep], up[keep]
    a = atr[idx2]; base = c[idx2]
    rate = touched(cu, cl, up, base + off).mean()

    rng = np.random.default_rng(seed)
    r = np.abs(off) / a
    sgn = np.sign(off)
    pr = np.empty(nperm)
    for k in range(nperm):
        pr[k] = touched(cu, cl, up, base + sgn * permute_within(r, up, rng) * a).mean()
    mu, sd = pr.mean(), max(pr.std(ddof=1), 1e-12)
    z = (rate - mu) / sd
    if not quiet:
        print(f"  {sym:>7} {kind:<7} n={len(idx2):>7,}  touch {100*rate:6.2f}%  "
              f"null {100*mu:6.2f}%  diff {100*(rate-mu):+6.2f}pp  z {z:+7.2f}")
    return rate, mu, sd, z, len(idx2), (pr - mu) / sd


def reaction_test(sym, kind, H=288, K=48, w=0.30, every=12, nperm=200, seed=11,
                  quiet=False, pre=None):
    """CONDITIONAL ON REACHING the level, does anything happen there?

    At the touch bar, put symmetric barriers w*ATR either side of the level and ask which
    is hit first. Rejection means price turns back the way it came. Under no effect this is
    ~0.5 by symmetry, so the test needs no model of drift - and the same permuted-offset
    placebo controls for the fact that arriving anywhere after a directional run is not a
    neutral state.
    """
    idx, off, up, atr, h, l, c, n = pre if pre else collect(sym, kind, H, every)
    if len(idx) < 500:
        return None
    idx2, keep, cu, cl = forward_cummax(h, l, idx, H)
    off, up = off[keep], up[keep]
    a = atr[idx2]; base = c[idx2]

    def run(offsets):
        px = base + offsets
        j = touch_bar(cu, cl, up, px, idx2)
        m = (j >= 0) & (j + K < n)
        if m.sum() < 100:
            return np.nan, 0
        res = race(h, l, j[m], px[m] - w * a[m], px[m] + w * a[m], K)
        dec = res != 0
        if dec.sum() < 100:
            return np.nan, 0
        #--- "reject" = the barrier AWAY from the direction of approach is hit first
        away = np.where(up[m], -1, 1)[dec]
        return float((res[dec] == away).mean()), int(dec.sum())

    rate, nn = run(off)
    if not np.isfinite(rate):
        return None
    rng = np.random.default_rng(seed)
    r = np.abs(off) / a; sgn = np.sign(off)
    pr = []
    for k in range(nperm):
        v, _ = run(sgn * permute_within(r, up, rng) * a)
        if np.isfinite(v):
            pr.append(v)
    pr = np.array(pr)
    mu, sd = pr.mean(), max(pr.std(ddof=1), 1e-12)
    z = (rate - mu) / sd
    if not quiet:
        print(f"  {sym:>7} {kind:<7} touches={nn:>7,}  reject {100*rate:6.2f}%  "
              f"null {100*mu:6.2f}%  diff {100*(rate-mu):+6.2f}pp  z {z:+7.2f}")
    return rate, mu, sd, z, nn, (pr - mu) / sd


def eighty_rule(sym, accept_bars=12, H=288):
    """Market Profile's 80% rule, with an exact martingale benchmark.

    Setup: the session OPENS outside the previous session's value area, later trades back
    inside, and is ACCEPTED (`accept_bars` consecutive M5 closes inside - 12 = two 30-minute
    periods, which is the classic formulation). Claim: ~80% chance of traversing the whole
    value area to the far edge.

    Benchmark: from the acceptance price p, a driftless process reaches the far edge before
    the near one with probability (distance to near edge)/(width). Quoted per event and
    averaged, so the comparison is against what geometry alone already delivers.
    """
    a, I = load_bars(sym)
    g = lambda c: a[:, I[c]]
    t, o, h, l, c = g('time'), g('open'), g('high'), g('low'), g('close')
    bd = broker_day(t.astype(np.int64))
    sess, nodes, naked, binsize = build_levels(sym)
    sday = {int(d): i for i, d in enumerate(sess[:, 0])}
    starts = np.concatenate(([0], np.flatnonzero(np.diff(bd)) + 1))
    ends = np.concatenate((starts[1:], [len(bd)]))

    #--- three readings of the same rule, because they answer different questions:
    #---   acc  = tradeable, after acceptance   (far edge BEFORE the near edge)
    #---   re   = tradeable, on bare re-entry   (does 'acceptance' add anything?)
    #---   lit  = literal, after acceptance     (far edge touched AT ALL in the session)
    res = {k: ([], []) for k in ('acc', 're')}
    lit = []
    for s, e in zip(starts, ends):
        d = int(bd[s]); si = sday.get(d)
        if si is None or si < 1 or e - s < 60:
            continue
        val = sess[si - 1, 2] * binsize
        vah = sess[si - 1, 3] * binsize
        if vah - val <= 0:
            continue
        op = o[s]
        if val <= op <= vah:
            continue                      # opened inside value - not this setup
        from_above = op > vah
        inside = (c[s:e] >= val) & (c[s:e] <= vah)
        near, far = (vah, val) if from_above else (val, vah)

        run = 0; acc = -1; re = -1
        for k in range(len(inside)):
            if inside[k] and re < 0:
                re = s + k
            run = run + 1 if inside[k] else 0
            if run >= accept_bars:
                acc = s + k
                break
        for tag, at in (('acc', acc), ('re', re)):
            if at < 0 or at + 1 >= e:
                continue
            r = race(h, l, np.array([at]), np.array([min(near, far)]),
                     np.array([max(near, far)]), H)[0]
            if r == 0:
                continue
            res[tag][0].append((r == -1) if from_above else (r == 1))
            res[tag][1].append(abs(c[at] - near) / (vah - val))
        if acc >= 0 and acc + 1 < e:
            seg = slice(acc, e)
            lit.append(bool((l[seg] <= far).any() if from_above else (h[seg] >= far).any()))

    for tag, label in (('acc', 'after acceptance'), ('re', 'on bare re-entry')):
        hits = np.array(res[tag][0], bool); bench = np.array(res[tag][1])
        if len(hits) < 30:
            print(f"  {sym}: only {len(hits)} qualifying sessions ({label})")
            continue
        var = (bench * (1 - bench)).sum()
        z = (hits.sum() - bench.sum()) / max(np.sqrt(var), 1e-9)
        print(f"  {sym:>7} {label:<18} n={len(hits):>5}  traverse {100*hits.mean():6.2f}%  "
              f"martingale {100*bench.mean():6.2f}%  "
              f"diff {100*(hits.mean()-bench.mean()):+6.2f}pp  z {z:+6.2f}")
    if lit:
        print(f"  {sym:>7} {'literal (no stop)':<18} n={len(lit):>5}  far edge touched at "
              f"some point in the session: {100*np.mean(lit):5.2f}%   "
              f"-- no stop, so no benchmark and no money in it")
    return res


if __name__ == '__main__':
    which = sys.argv[1] if len(sys.argv) > 1 else 'all'
    syms = [s for s in sys.argv[2:] if s in SYMS] or list(SYMS)
    KINDS = ('naked', 'vpoc', 'vaedge', 'hvn', 'lvn')

    #--- collect() re-decodes the bars and rebuilds every level, ~45 s a time. Do it once
    #--- per (symbol, kind) and hand the same event set to both tests.
    PRE = {}
    for sym in syms:
        for kind in KINDS:
            try:
                PRE[(sym, kind)] = collect(sym, kind)
            except FileNotFoundError:
                pass

    def family(title, fn):
        print(title)
        acc = []
        for sym in syms:
            for kind in KINDS:
                p = PRE.get((sym, kind))
                if p is None:
                    continue
                r = fn(sym, kind, pre=p)
                if r:
                    acc.append(np.abs(r[5]))
        if acc:
            m = min(len(x) for x in acc)
            crit = float(np.quantile(np.maximum.reduce([x[:m] for x in acc]), 0.95))
            print(f"  family-wise |z| bar over {len(acc)} tests: {crit:.2f}")

    if which in ('all', 'magnet'):
        family("\n=== 1. MAGNET: is a real level reached more often than a "
               "distance-matched placebo? (M5, H=288 bars = 1 session) ===", magnet_test)

    if which in ('all', 'reaction'):
        family("\n=== 2. REACTION: conditional on touching it, does price reject the "
               "level? (barriers +/-0.30 ATR, K=48 bars) ===", reaction_test)

    if which in ('all', 'eighty'):
        print("\n=== 3. THE 80% RULE: open outside value, accept back inside, "
              "traverse to the far edge? ===")
        for sym in syms:
            eighty_rule(sym)