Warrior_EA/research/test_cause_effect.py

"""Wyckoff's law of cause and effect: does a bigger cause really produce a bigger effect?

This is the load-bearing claim of the whole method and the one that decides where a TAKE
PROFIT goes. Villahermosa states it plainly (book 1, ch. 8): "the effect produced by the
cause will always be directly proportional to that cause... the longer the market spends in
a range building a campaign, the further the subsequent trend move will travel." He rejects
point-and-figure counting as too subjective and recommends the objective version instead -
project the range's height 1:1 from the breakout.

Two testable halves, neither tested anywhere else in this project:

  1. Is the projected target reachable at a profit?  Enter on the breakout, stop at the far
     side of the range, target k x risk. A driftless market gives P(win) = 1/(1+k) EXACTLY,
     so expected R is exactly zero for every k and no permutation null is needed - the
     benchmark is analytic. Any k with positive expectancy after spread is a real finding.
  2. Is the effect proportional to the cause?  Does a LONGER or TALLER range produce more
     R-multiples of travel? This is the actual Wyckoff claim, and it is the one that would
     tell the EA how far to aim.

WHY THE UNITS MATTER
--------------------
Measuring the move in ATR is nearly useless here: a fixed horizon in a volatile market
produces a big number whatever the setup was. Everything below is in RANGE HEIGHTS (R),
i.e. multiples of the risk actually taken, which is the only unit in which "where should
the take profit go" has an answer.

Two biases are deliberately left in place, both against the hypothesis, so a positive result
cannot be an artifact of generous accounting:
  - a bar that spans stop and target books the LOSS,
  - the spread is charged on the entry and on both barriers.

RANGE DETECTION - objective, causal, volatility-normalised
----------------------------------------------------------
For a random walk the span of L bars grows like ATR*sqrt(L), so

    compression = span(L) / (ATR * sqrt(L))

is ~1 for ordinary wandering and well below 1 when price is genuinely contained. Taking the
LONGEST L whose compression stays under a threshold finds the range with no threshold in
price units - it means the same thing on gold and on EURUSD - and never looks at a bar
after the one being evaluated.
"""
import numpy as np, sys, os
sys.stdout.reconfigure(encoding='utf-8', errors='replace')

BARS = 'c:/Users/admin/Documents/Workspaces/Market Data/bars/'
SYMS = ('EURUSD', 'USDJPY', 'XAUUSD', 'SP500')
KS = (0.5, 1.0, 1.5, 2.0, 3.0, 5.0)


def load_bars(sym, tf):
    z = np.load(f"{BARS}{sym}_{tf}_ticks.npz", allow_pickle=True)
    return z['bars'], {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 rolling_span(h, l, L):
    """(max high, min low) over the L bars ENDING AT i-1. Never includes bar i."""
    from numpy.lib.stride_tricks import sliding_window_view
    n = len(h)
    hi = np.full(n, np.nan); lo = np.full(n, np.nan)
    if n > L:
        hi[L:] = sliding_window_view(h, L).max(axis=1)[:-1]
        lo[L:] = sliding_window_view(l, L).min(axis=1)[:-1]
    return hi, lo


def find_ranges(h, l, c, atr, Ls=(12, 18, 24, 36, 48, 72, 96, 144), theta=0.60):
    """For every bar, the LONGEST candidate window that qualifies as a range.

    Longest-first so a 24-bar consolidation is not recorded as a 12-bar one; the duration
    test needs the real duration and truncating it would flatten the very effect under test.
    """
    n = len(c)
    L_out = np.zeros(n, np.int32)
    hi_out = np.full(n, np.nan); lo_out = np.full(n, np.nan)
    for L in sorted(Ls, reverse=True):
        hi, lo = rolling_span(h, l, L)
        comp = (hi - lo) / np.maximum(atr * np.sqrt(L), 1e-12)
        ok = (L_out == 0) & np.isfinite(comp) & (comp <= theta)
        L_out[ok] = L; hi_out[ok] = hi[ok]; lo_out[ok] = lo[ok]
    return L_out, hi_out, lo_out


def race_first(h, l, e, d, stop_px, targ_px, H):
    """+1 target first, -1 stop first, 0 unresolved. Stop wins a bar that spans both.

    H is PER EVENT, so a 12-bar range is not silently given the 144-bar range's horizon -
    that alone would fabricate a duration effect out of nothing.
    """
    n = len(h)
    H = np.asarray(H)
    out = np.zeros(len(e), np.int8)
    live = np.ones(len(e), bool)
    for k in range(0, int(H.max()) + 1):
        j = e + k
        ok = live & (j < n) & (k <= H)
        if not ok.any():
            continue
        jj = j[ok]
        lose = np.where(d[ok] > 0, l[jj] <= stop_px[ok], h[jj] >= stop_px[ok])
        winb = np.where(d[ok] > 0, h[jj] >= targ_px[ok], l[jj] <= targ_px[ok])
        res = np.where(lose, -1, np.where(winb, 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 events(sym, tf, theta=0.60, hk=6, hcap=600):
    """Range breakouts, with a horizon PROPORTIONAL to the cause.

    A 144-bar consolidation given the same 240 bars to work with as a 12-bar one would be
    judged mostly on the horizon rather than on the setup, and that alone would manufacture
    a negative duration effect. Horizon = hk x range length, capped.
    """
    a, I = load_bars(sym, tf)
    g = lambda k: a[:, I[k]]
    o, h, l, c = g('open'), g('high'), g('low'), g('close')
    spb = g('spread_mean'); spb = np.concatenate([[spb[0]], spb[:-1]])
    atr = atr_of(h, l, c, 14); atr = np.concatenate([[atr[0]], atr[:-1]])
    n = len(c)
    L, hi, lo = find_ranges(h, l, c, atr, theta=theta)
    up = (L > 0) & (c > hi); dn = (L > 0) & (c < lo)
    fire = np.nonzero(up | dn)[0]
    H = np.minimum(L[fire] * hk, hcap)
    fire_ok = (fire > 200) & (fire + H < n - 2)
    fire, H = fire[fire_ok], H[fire_ok]
    if not len(fire):
        return None
    d = np.where(up[fire], 1, -1)
    #--- non-overlapping: the next event may not start before this one's horizon expires
    keep, busy = [], -1
    for k, i in enumerate(fire):
        if i <= busy:
            continue
        keep.append(k); busy = i + H[k]
    keep = np.array(keep, int)
    f, d, H = fire[keep], d[keep], H[keep]
    e = f + 1                                   # enter at the NEXT open, never the signal close
    return dict(sym=sym, tf=tf, f=f, e=e, d=d, H=H, hi=hi[f], lo=lo[f], L=L[f].astype(float),
                atr=atr[f], ent=o[e], sp=spb[e], h=h, l=l, c=c, n=n)


def trade(E, k):
    """The breakout trade at a target of k x risk. Returns (R per trade, resolved, ok).

    An unresolved trade is CLOSED AT THE HORIZON, not discarded. Discarding it is the
    subtle bias that made the first version of this test look worse than reality at large
    k: stops resolve fast and stay in the sample, while distant targets resolve slowly and
    a truncated horizon quietly deletes exactly the winners. Marking to market at the last
    bar keeps every trade in the accounting where it belongs.
    """
    d, e, ent, sp = E['d'], E['e'], E['ent'], E['sp']
    far = np.where(d > 0, E['lo'], E['hi'])       # stop at the OTHER side of the range
    risk = np.abs(ent - far) + sp
    ok = risk > 0
    stop = ent - d * risk
    targ = ent + d * k * risk
    res = race_first(E['h'], E['l'], e, d, stop, targ, E['H'])
    R = np.where(res > 0, k, np.where(res < 0, -1.0, 0.0))
    un = res == 0
    if un.any():
        j = np.minimum(e[un] + E['H'][un], len(E['c']) - 1)
        R[un] = (E['c'][j] - ent[un]) * d[un] / np.maximum(risk[un], 1e-12) \
                - sp[un] / np.maximum(risk[un], 1e-12)
    return R, res, ok


def run(sym, tf, theta=0.60, verbose=True):
    E = events(sym, tf, theta)
    if E is None or len(E['e']) < 80:
        print(f"  {sym:>7} {tf}: too few range breakouts")
        return None
    rows = []
    for k in KS:
        R, res, ok = trade(E, k)
        if ok.sum() < 50:
            continue
        R = R[ok]
        m = R.mean()
        #--- t on the realised R series: the analytic benchmark is exactly 0 expected R
        t = m / max(R.std(ddof=1) / np.sqrt(len(R)), 1e-12)
        dec = res[ok] != 0
        w = (res[ok] > 0).sum() / max(dec.sum(), 1)
        rows.append((k, int(ok.sum()), 100 * w, 100 / (1 + k), m, t,
                     100 * (1 - dec.mean())))
    if verbose and rows:
        print(f"  {sym:>7} {tf}  {len(E['e'])} breakouts")
        print(f"          {'k':>5}{'n':>7}{'win%':>8}{'martingale%':>13}{'expR':>9}"
              f"{'t':>8}{'unres%':>8}")
        for r in rows:
            print(f"          {r[0]:>5.1f}{r[1]:>7}{r[2]:>8.2f}{r[3]:>13.2f}"
                  f"{r[4]:>+9.3f}{r[5]:>+8.2f}{r[6]:>8.1f}")
    return E, rows


def proportionality(sym, tf, theta=0.60, _cache={}):
    """The Wyckoff claim itself: is travel (in range-heights) larger after a bigger cause?

    Reported as median MFE / R by duration and by height. If the law holds these columns
    rise. If they are flat, the projection is a constant rule and 'proportional to the
    cause' adds nothing to it.
    """
    E = events(sym, tf, theta)
    if E is None or len(E['e']) < 80:
        return None
    d, e, ent = E['d'], E['e'], E['ent']
    far = np.where(d > 0, E['lo'], E['hi'])
    risk = np.abs(ent - far) + E['sp']
    Hm = int(E['H'].max())
    from numpy.lib.stride_tricks import sliding_window_view
    keep = e + Hm < E['n']
    e = e[keep]; d = d[keep]; ent = ent[keep]; risk = risk[keep]
    E = dict(E, e=e, d=d, ent=ent, H=E['H'][keep], L=E['L'][keep],
             lo=E['lo'][keep], hi=E['hi'][keep], atr=E['atr'][keep])
    hh = sliding_window_view(E['h'], Hm)[e]
    ll = sliding_window_view(E['l'], Hm)[e]
    #--- mask each event to its own horizon
    cols = np.arange(Hm)[None, :]
    mask = cols <= E['H'][:, None]
    hh = np.where(mask, hh, -np.inf)
    ll = np.where(mask, ll, np.inf)
    mfe = np.where(d > 0, hh.max(axis=1) - ent, ent - ll.min(axis=1)) / np.maximum(risk, 1e-12)
    return E, mfe


def exponent(sym, tf, theta=0.60):
    """How does the effect actually scale with the cause?

    Wyckoff's law says effect ~ cause, i.e. an exponent of 1 in

        log(MFE / ATR)  =  a  +  b * log(height / ATR)

    b = 1 means the 1:1 projection is the right rule at every size. b = 0 means the size of
    the range says nothing about how far price then travels. Anything in between says the
    projection over-reaches for big ranges and under-reaches for small ones - which is a
    more useful answer than a yes/no, because b is directly the exponent a target rule
    should use.
    """
    pr = proportionality(sym, tf, theta)
    if pr is None:
        return None
    E, mfe_r = pr
    d = E['d']; ent = E['ent']
    far = np.where(d > 0, E['lo'], E['hi'])
    risk = np.abs(ent - far)
    #--- ATR must be a FREE regressor, not a shared denominator. Dividing both sides by the
    #--- same noisy ATR correlates the errors and biases b TOWARDS 1 - i.e. towards the
    #--- hypothesis - so the constrained form cannot be used to argue against it.
    x = np.log(np.maximum(risk, 1e-12))
    z = np.log(np.maximum(E['atr'], 1e-12))
    y = np.log(np.maximum(mfe_r * risk, 1e-12))
    A = np.column_stack([np.ones(len(x)), x, z])
    beta, *_ = np.linalg.lstsq(A, y, rcond=None)
    resid = y - A @ beta
    s2 = resid @ resid / max(len(x) - 3, 1)
    se = np.sqrt(s2 * np.linalg.inv(A.T @ A)[1, 1])
    return float(beta[1]), float(se), len(x)


def by_cause_size(sym, tf, theta=0.60):
    """Sublinearity would be monetisable only if the RIGHT k differs by cause size and the
    expectancy at that k is positive somewhere. Split the breakouts by range height and
    price the whole k grid inside each half."""
    E = events(sym, tf, theta)
    if E is None or len(E['e']) < 200:
        return
    far = np.where(E['d'] > 0, E['lo'], E['hi'])
    hgt = np.abs(E['ent'] - far) / E['atr']
    med = np.median(hgt)
    for lbl, m in (('small cause', hgt <= med), ('big cause', hgt > med)):
        best, cells = None, []
        for k in KS:
            R, res, ok = trade(E, k)
            sel = ok & m
            if sel.sum() < 60:
                cells.append('    n/a'); continue
            r = R[sel]
            t = r.mean() / max(r.std(ddof=1) / np.sqrt(len(r)), 1e-12)
            cells.append(f"{r.mean():+7.3f}")
            if best is None or r.mean() > best[1]:
                best = (k, r.mean(), t, int(sel.sum()))
        print(f"  {sym:>7} {tf} {lbl:<12} expR by k: " + " ".join(cells) +
              (f"   best k={best[0]} expR {best[1]:+.3f} t {best[2]:+.2f} n={best[3]}"
               if best else ""))


if __name__ == '__main__':
    syms = [s for s in sys.argv[1:] if s in SYMS] or list(SYMS)
    TFS = ('M15', 'H1')
    print("=== 1. THE PROJECTED TARGET AS A TRADE ===")
    print("    enter on the range breakout, stop at the far side of the range, target k x risk.")
    print("    A driftless market gives win% = 100/(1+k) and expR = 0 at EVERY k.\n")
    pool = {}
    for tf in TFS:
        print(f"--- {tf} ---")
        for sym in syms:
            r = run(sym, tf)
            if r:
                pool[(sym, tf)] = r

    print("\n=== 2. IS THE EFFECT PROPORTIONAL TO THE CAUSE? ===")
    print("    travel before the stop, in units of the risk taken (MFE / R).")
    print("    Wyckoff says these columns should RISE with the size of the cause.\n")
    for tf in TFS:
        for sym in syms:
            pr = proportionality(sym, tf)
            if pr is None:
                continue
            E, mfe = pr
            L = E['L']; height = np.abs(E['ent'] - np.where(E['d'] > 0, E['lo'], E['hi'])) / E['atr']
            print(f"  {sym:>7} {tf}  n={len(mfe)}")
            for lbl, key in (('duration', L), ('height/ATR', height)):
                qs = np.quantile(key, [0, .25, .5, .75, 1.0])
                cells = []
                for a, b in zip(qs[:-1], qs[1:]):
                    m = (key >= a) & (key <= b)
                    if m.sum() < 20:
                        cells.append('   n/a'); continue
                    cells.append(f"{np.median(mfe[m]):6.2f}")
                print(f"            {lbl:<11} Q1..Q4 median MFE/R: " + "  ".join(cells))

    print("\n=== 3. THE SCALING EXPONENT ===")
    print("    log(MFE/ATR) = a + b.log(height/ATR).  Wyckoff's law is b = 1.")
    print("    b = 0 would mean the size of the cause says nothing at all.\n")
    print(f"  {'symbol':>8} {'tf':>4}{'n':>7}{'b':>9}{'se':>7}   {'b=1?':>12}   {'b=0?':>10}")
    for tf in TFS:
        for sym in syms:
            r = exponent(sym, tf)
            if r is None:
                continue
            b, se, n = r
            print(f"  {sym:>8} {tf:>4}{n:>7}{b:>9.3f}{se:>7.3f}   "
                  f"{(b-1)/se:>+9.2f} sd   {b/se:>+7.2f} sd")

    print("\n=== 4. CAN THE SUBLINEARITY BE TRADED? expR by cause size and target ===")
    print(f"    k = {', '.join(str(k) for k in KS)}\n")
    for tf in TFS:
        for sym in syms:
            by_cause_size(sym, tf)