# ar2: arrival-pattern features for the arb/retail call (Taylor's open problem 3)

Credit: Taylor (gtmcknight), github.com/gtmcknight/amm-sim @ 4812a6382311624e1506aca131ca35c378f930d9, research/solved. All code and tables are his; only the classifier hook below is new. Replica only (no engine run).

## Verdict: NOT worth porting
Report-set gain is **+0.011, se 0.009** (threshold was > +0.3 with se < 0.1). The three arrival-pattern features carry almost no information beyond the filter's own likelihood ratio. The true flag is worth **+3.72** on this plain-table arm, so the gap is real, but this idea closes about 0.3% of it.

## Gate P: PASS
Unmodified `python3 -I gates.py` printed gaussian+plain 443.455598 / 571.463409 / 525.168028 for seeds 10000..10002, equal to the README block (also the other columns matched). BASE arm mean over 10000..12047 = 532.386 (README: 532.39).

## What the idea computes
In `RealPolicy.on_trade`, on the first own trade of a step, real.py gets pi = L_arb/(L_arb+L_ret). TEST replaces it with pi' = sigmoid(c0 + cl*lgt + cg1*[dt==1] + cg2*[dt==2] + cp*prev_pi + cd*ln(pside)), where lgt = ln(L_arb/L_ret) (clipped to +-30), dt = step gap since the last own swap, prev_pi = pi' the previous step's first trade got (the inferred previous type), pside = mass of the pre-trade belief N(m,v) that lies beyond the quote band on the side this trade pushed (the filter's predicted direction of mispricing; floor 1e-4). pi' then feeds real.py's unchanged hard/soft belief update and tables (soft=1, pol_mid.json, plain tables).
Selected setting (choose seeds): c0=0, cl=1, cg1=cg2=cd=0, **cp=0.2** (i.e. only the previous-type feature). The other settings are listed below; all were no better.

## Choose set, seeds 60000..61023 (TEST minus BASE, edge)
Logistic fit of all features on 284k sampled first-trade records (true flag as label) improved log-likelihood by only 556 (-106468 -> -105912) and classification error 16.74% -> 16.71%.

| setting | mean_diff | se |
|---|---|---|
| all features, fitted (c0 -0.6662, cl 0.8821, cg1 0.4921, cg2 0.1307, cp 0.1178, cd -0.1867) | -0.4175 | 0.0189 |
| gap only (cg1 0.4, cg2 0.1) | -0.0957 | 0.0182 |
| direction only (cd -0.12) | -0.1340 | 0.0173 |
| bias +0.3 | -0.0958 | 0.0181 |
| bias -0.3 | -0.1397 | 0.0168 |
| prev only cp 0.5 | -0.0333 | 0.0165 |
| **prev only cp 0.2 (chosen)** | **+0.0063** | **0.0119** |

Chosen setting on choose seeds: n 1024, mean_diff +0.0063, se 0.0119, wins 516, worst5 (60725,-1.182) (60211,-1.171) (60225,-1.160) (60086,-1.065) (60092,-1.055).

## Report set, seeds 10000..12047 (scored once, same seeds both arms)
- TEST minus BASE: n 2048, mean_diff **+0.0110**, se **0.0094**, wins 1048, worst5 (11701,-3.676) (11510,-1.882) (11021,-1.602) (10359,-1.493) (10949,-1.323).
- Means: BASE 532.386, TEST 532.397.
- Diagnostic, true flag (`oraclearb=1`, same plain tables, soft=1) minus BASE: n 2048, mean_diff **+3.7226**, se 0.0237, wins 2044. (README quotes +6.98 total vs HA with qm tables and true flag, and 2.7 as the worth; this plain arm's gap is larger.) Oracle mean 536.108.

## Confusion matrix, report seeds (first own trade of each step; 14,218,237 calls)
| | called arb | called retail |
|---|---|---|
| BASE true arb | 6,315,679 | 1,287,954 |
| BASE true retail | 1,071,919 | 5,542,685 |
| TEST true arb | 6,399,972 | 1,203,661 |
| TEST true retail | 1,159,050 | 5,455,554 |

Error rate BASE 16.60%, TEST 16.62%. TEST trades arb misses for retail false alarms about one for one, so the call is no better; the mispricing direction and the gap are already inside the filter's likelihood.

## Time
BASE about 0.12 CPU-s per seed (0.06 s/seed wall on 2 workers, 2048 seeds in about 2 min); TEST about the same (0.07-0.11 s wall; the collecting run logs records).

## Complete changed file (real2.py = real.py plus the ar2 block and per-record logging)
sha256 of the file as run: `cc348f59b4e1b486b71dfe122f19924150399b17adb6612ffc130752c2e85261`. The `collect`/`rec` lines are logging only. Run via an `opt` dict with ar=1,soft=1,c0=0,cl=1,cg1=0,cg2=0,cp=0.2,cd=0 (RealPolicy(T, opt)); it needs the repo's replica2.py and evalpol.py beside it.

```python
"""Solved policy driven ONLY by what afterSwap sees (side, amounts, timestamp, reserves).

Filter: Gaussian belief over ln p (mean, variance).
  - new step: variance grows by sigma^2 a step; steps with no trade of ours truncate the belief to our no-arb interval
  - first trade of a step: arb or retail by likelihood ratio. Arb => p is exact (post-trade spot through the fee in force).
    Retail => belief truncated to the no-arb interval.
Rival: absolute spot dead-reckoned (router identity on own retail trades, band clip against the price estimate).
Policy: table lookup at the estimates.  Oracle arguments are used for diagnostics only.

usage: python3 -I real.py <pol.json> [variant ...]      variants: ce | qm:<qmdp.json> ...
"""
import sys, os, math, json, time
from multiprocessing import Pool
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
import replica2 as R
from evalpol import DPPolicy, load, stats, CALIB

G2 = 1.0 - R.NF
SIG2 = 0.000945 ** 2
RHO = 0.26                       # P(a no-arb step's first own trade is retail on a given side)
LD_M, LD_S = math.log(11.12e-4), 1.176   # own-fill spot move of a retail first trade: lognormal fit
SQ2PI = 2.5066282746310002
H0 = 0.55
LD_FULL = math.log(20.0) - 0.72 + math.log(2.0 / 10000.0)   # ln median of a full order's displacement


def phi(x): return math.exp(-0.5 * x * x) / SQ2PI
def Phi(x): return 0.5 * (1.0 + math.erf(x / 1.4142135623730951))


def trunc(m, v, lo, hi):
    """moments of N(m, v) restricted to (lo, hi); returns (mass, mean, var)"""
    if v <= 0.0:
        return (1.0 if lo <= m <= hi else 0.0), m, 0.0
    sd = math.sqrt(v); a = (lo - m) / sd; b = (hi - m) / sd
    Z = Phi(b) - Phi(a)
    if Z < 1e-12:
        # belief is outside the interval: snap to the nearer edge
        return Z, (lo if m < lo else hi), v * 0.25
    pa, pb = phi(a), phi(b)
    r = (pa - pb) / Z
    mean = m + sd * r
    var = v * (1.0 + (a * pa - b * pb) / Z - r * r)
    if var < 0.0: var = 0.0
    return Z, mean, var


class RealPolicy(DPPolicy):
    wants_j = True
    def __init__(self, T, opt=None):
        DPPolicy.__init__(self, T)
        self.o = opt or {}
        self.Q = load(self.o['qm']) if 'qm' in self.o else None
        self.synced = False
        self.m = math.log(100.0); self.v = 0.0
        self.tlast = 0; self.cnt = 0
        self.s2h = 100.0
        self.fb = self.fa = 0.003
        self.sp = 100.0               # own spot after our last callback
        # diagnostics
        self.prev_pi = 0.0; self.rec = []
        self.nmis = 0; self.ncls = 0; self.se = 0.0; self.nse = 0; self.se2 = 0.0
    def init(self, *a):
        fb, fa = self._fees(0.0, 0.0, 0)
        self.fb, self.fa = R.clampf(fb), R.clampf(fa)
        return fb, fa
    def lookup(self, w, w2, phs, sd):
        Q = self.Q
        if Q is None:
            return self._fees(w, w2, phs)
        sd = sd * self.o.get('sdscale', 1.0)
        S = Q['S']; k = 0
        while k < len(S) - 2 and sd >= S[k + 1]: k += 1
        tt = (sd - S[k]) / (S[k + 1] - S[k])
        if tt < 0.0: tt = 0.0
        if tt > 1.0: tt = 1.0
        tag = '0c' if phs == 0 else (str(phs) + ('s' if self.synced else 'c'))
        A, B = Q['A' + tag], Q['B' + tag]
        pw = w / self.h + self.W0
        if pw < 0: pw = 0.0
        if pw > self.NW - 1: pw = float(self.NW - 1)
        l2 = w2 / self.h
        if phs == 0:
            if l2 < -self.FB: l2 = -self.FB
            if l2 > self.FB: l2 = self.FB
            p2 = l2 + self.FB; n2 = 2 * self.FB + 1
        else:
            p2 = self._pos2(l2); n2 = self.N2
        i = int(pw); j = int(p2)
        if i > self.NW - 2: i = self.NW - 2
        if j > n2 - 2: j = n2 - 2
        f = pw - i; g = p2 - j
        def bl(Tb):
            return (Tb[i][j] * (1 - f) + Tb[i + 1][j] * f) * (1 - g) + (Tb[i][j + 1] * (1 - f) + Tb[i + 1][j + 1] * f) * g
        a = bl(A[k]) * (1 - tt) + bl(A[k + 1]) * tt
        b = bl(B[k]) * (1 - tt) + bl(B[k + 1]) * tt
        fa = 1.0 - math.exp(-(a - w) * 1e-4) if a > w else 0.0
        fb = 1.0 - math.exp(-(b + w) * 1e-4) if b > -w else 0.0
        return fb, fa
    def on_trade(self, is_buy, ax, ay, t, x, y, arb, p, s2, j):
        s = y / x
        dt = t - self.tlast
        is_arb = False
        if dt >= 1:
            self.cnt = 0
            lo = math.log(self.sp * (1.0 - self.fb)); hi = math.log(self.sp / (1.0 - self.fa))   # no-arb interval of the quotes in force
            m, v = self.m, self.v
            for _ in range(min(dt - 1, 12)):          # fully silent steps
                v += SIG2
                _, m, v = trunc(m, v, lo, hi)
            v += SIG2                                  # this step's move
            # arb hypothesis: p exactly at the post-trade band edge
            lp_arb = math.log(s * (1.0 - self.fb)) if is_buy else math.log(s / (1.0 - self.fa))
            sd = math.sqrt(v)
            o = self.o
            l_arb = phi((lp_arb - m) / sd) / sd
            eps = o.get('tail', 0.0)
            if eps > 0.0:
                l_arb = (1.0 - eps) * l_arb + eps * phi((lp_arb - m) / 40e-4) / 40e-4
            Z, mt, vt = trunc(m, v, lo, hi)
            d = abs(math.log(s / self.sp))
            if o.get('fill', 0):
                # retail fill density from the order-size law and how our quote compared with the rival's (estimates)
                pe0 = math.exp(mt)
                r2 = self.s2h
                if r2 < pe0 * G2: r2 = pe0 * G2
                elif r2 > pe0 / G2: r2 = pe0 / G2
                if is_buy:   # trader sold: compare bids (log)
                    delta = math.log(r2 * G2) - math.log(self.sp * (1.0 - self.fb))      # >0: rival bid better
                else:
                    delta = math.log(self.sp / (1.0 - self.fa)) - math.log(r2 / G2)      # >0: rival ask better
                if d > -delta:
                    D = 2.0 * d + delta; jac = 2.0
                else:
                    D = d; jac = 1.0
                D = max(D, 1e-9)
                l_ret = Z * 0.5 * H0 * jac * phi((math.log(D) - LD_FULL) / 1.2) / (D * 1.2)
            else:
                l_ret = Z * RHO * phi((math.log(max(d, 1e-9)) - LD_M) / LD_S) / (max(d, 1e-9) * LD_S)
            pi = l_arb / (l_arb + l_ret) if (l_arb + l_ret) > 0 else 0.0
            # --- ar2: arrival-pattern features
            if is_buy: pside = Phi((lo - m) / sd)
            else: pside = 1.0 - Phi((hi - m) / sd)
            lgt = math.log(max(l_arb, 1e-300)) - math.log(max(l_ret, 1e-300))
            lgt = max(-30.0, min(30.0, lgt))
            if o.get('collect', 0):
                self.rec.append((1 if arb else 0, lgt, dt, self.prev_pi, pside, 1 if is_buy else 0))
            if o.get('ar', 0):
                ps = min(max(pside, 1e-4), 1.0)
                x = (o['c0'] + o['cl'] * lgt + o['cg1'] * (1.0 if dt == 1 else 0.0) + o['cg2'] * (1.0 if dt == 2 else 0.0)
                     + o['cp'] * self.prev_pi + o['cd'] * math.log(ps))
                x = max(-40.0, min(40.0, x))
                pi = 1.0 / (1.0 + math.exp(-x))
            self.prev_pi = pi
            if o.get('oraclearb', 0): pi = 1.0 if arb else 0.0
            is_arb = pi > 0.5
            self.ncls += 1
            if is_arb != arb: self.nmis += 1
            if o.get('soft', 0) and pi < o.get('hard', 0.995) and pi > 1.0 - o.get('hard', 0.995):
                m2 = pi * lp_arb + (1.0 - pi) * mt
                v2 = pi * (lp_arb - m2) ** 2 + (1.0 - pi) * (vt + (mt - m2) ** 2)
                self.m, self.v = m2, v2
            elif is_arb:
                self.m, self.v = lp_arb, 0.0
            else:
                self.m, self.v = mt, vt
            # rival was arbed at the start of this step: band clip against the estimate
            pe = math.exp(self.m)
            if self.s2h < pe * G2: self.s2h = pe * G2
            elif self.s2h > pe / G2: self.s2h = pe / G2
            self.synced = False
        pe = math.exp(self.m)
        if not is_arb:
            if is_buy:
                cand = s * (1.0 - self.fb) / G2
                if cand < self.s2h: self.s2h = cand; self.synced = True
            else:
                cand = s * G2 / (1.0 - self.fa)
                if cand > self.s2h: self.s2h = cand; self.synced = True
            self.cnt += 1
            phs = 1 if self.cnt <= 1 else 2
        else:
            phs = 0
        err = (self.m - math.log(p)) * 1e4
        self.se += err * err; self.nse += 1
        if not is_arb: self.se2 += err * err
        self.tlast = t; self.sp = s
        w = (math.log(s) - self.m) * 1e4; w2 = (math.log(self.s2h) - self.m) * 1e4
        fb, fa = self.lookup(w, w2, phs, math.sqrt(self.v) * 1e4)
        self.fb, self.fa = R.clampf(fb), R.clampf(fa)
        return fb, fa


def ev(a):
    path, opt, seeds = a
    T = load(path); out = []
    for s in seeds:
        pol = RealPolicy(T, opt)
        r = R.run(s, pol)
        out.append((s, r['edge'], r['retail'], r['arb'], r['rv'], r['av'], pol.nmis, pol.ncls, pol.se, pol.nse, pol.rec))
    return out


def report(pool, path, opt, seeds, ref, label):
    t0 = time.time()
    res = pool.map(ev, [(path, opt, seeds[i:i + 8]) for i in range(0, len(seeds), 8)], chunksize=1)
    res = {r[0]: r for b in res for r in b}
    n = len(seeds); mean = lambda i: sum(res[s][i] for s in seeds) / n
    mh, seh = stats([res[s][1] - ref['HA'][s]['edge'] for s in seeds])
    mis = sum(res[s][6] for s in seeds) / max(1, sum(res[s][7] for s in seeds))
    rms = math.sqrt(sum(res[s][8] for s in seeds) / max(1, sum(res[s][9] for s in seeds)))
    print(f'{label:24s} edge {mean(1):8.3f} retail {mean(2):8.3f} arb {mean(3):7.3f} rv {mean(4):6.0f} av {mean(5):6.0f}  vs HA {mh:+.3f} se {seh:.3f}  misclass {mis:.2%}  p rms {rms:.1f} bps ({time.time()-t0:.0f}s)', flush=True)
    return res


if __name__ == '__main__':
    path = sys.argv[1]
    n = int(os.environ.get('NSEED', '2048'))
    seeds = list(range(10000, 10000 + n))
    rows = json.load(open(CALIB)); ref = {}
    for r in rows: ref.setdefault(r['kind'], {})[r['seed']] = r
    with Pool(14) as pool:
        for spec in (sys.argv[2:] or ['ce']):
            name, _, kv = spec.partition(':')
            opt = {}
            for item in kv.split(','):
                if item:
                    k, _, val = item.partition('=')
                    opt[k] = os.path.join(HERE, val) if k == 'qm' else float(val)
            report(pool, path, opt, seeds, ref, spec)
```