# r6-arb / ar3: quote against the whole posterior instead of its collapsed mean

Credit: the solved policy, replica, tables, Gaussian filter (`real.py`) and the list of open ideas all come from Taylor (gtmcknight), github.com/gtmcknight/amm-sim at commit 4812a6382311624e1506aca131ca35c378f930d9, folder `research/solved`. Everything below was run on that replica under plain `python3 -I`.

## Gate P: PASS

Unmodified `python3 -I gates.py` printed (3.6 s):

```
seed   perfect-info      gaussian+plain    gaussian+qm       lattice+qm
10000  448.261189  443.455598  443.522608  443.663520
10001  576.779199  571.463409  571.623551  571.977293
10002  531.886497  525.168028  525.865859  526.411631
```

This matches the README's gate block digit for digit, gaussian+plain 443.455598 / 571.463409 / 525.168028 included. My BASE arm, run through my harness, gives 443.4555976 / 571.4634094 / 525.1680278 on those seeds. Across all of 10000..12047 it scores +2.828 se 0.041 vs HA, which is the README's +2.828.

## Which idea, and why

README Open item 3 lists three untried ideas: (a) "choosing quotes against the whole lattice posterior instead of its mean and sd", (b) per-run sigma and order rate inside the classifier, and (c) the router likelihood per price hypothesis. Earlier items already covered (b) (+0.024) and (c) (+0.049). Size-to-gap and arrival-pattern are excluded. That leaves **(a)**, and it is also the one I would have picked on merit. The README says the classifier error "is near what those features allow" (table classifier 18.8% vs the lattice's 16.2%), and the two earlier items confirm it: both moved the result by less than 0.05. So a better *classifier* has little room left. What can still be recovered is the *cost* of each misclassification. Right now `real.py` collapses an ambiguous first trade into one Gaussian (m, v) by moment matching. With plain tables v is then ignored, so the quote is set at one price that may be far from both hypotheses. Quoting against the whole posterior hedges both cases.

Adapting (a) to the Gaussian filter (the brief requires a copy of `real.py`, so I used the Gaussian version rather than the lattice): the belief update is unchanged. The filter state m, v, the classification `pi > 0.5`, the rival tracking and the phase logic are all exactly as in BASE. The only change is how the plain-table quote is chosen.

## What the TEST arm computes

At the first own trade of a step, the filter already builds a two-part posterior: an arb point mass at `lp_arb` with weight pi, and a retail part N(m, v) truncated to the no-arb interval, with mass 1 - pi, mean mt and variance vt. TEST stores this as a list of price nodes `mix`:
- arb branch (pi >= 0.995): one node at `lp_arb`;
- retail branch (pi <= 0.005): 3-point Gauss-Hermite nodes mt - h, mt, mt + h with weights 1/6, 2/3, 1/6, where h = qs * sqrt(3 vt);
- soft branch: weight qa = pi^k / (pi^k + (1-pi)^k) on `lp_arb`, plus the three retail nodes sharing weight 1 - qa.

At every callback in that step (the filter does not update m within a step, so the posterior stays valid), each node x gives its own table lookup: own displacement w_x = ln s - x, rival displacement ln s2 - x, same phase. That gives offsets (a_x, b_x). The quoted log ask is the weighted mean of x + a_x and the quoted log bid is the weighted mean of x - b_x. The fees follow from the spot exactly as in BASE (ask fee = 1 - exp(-(LA - ln s)) if LA > ln s, else 0; bid likewise). With a single node this reduces exactly to BASE's `_fees`. Before the first step, or with `qm` tables, behavior is unchanged.

Parameters: qs (quadrature width scale) and qk (power sharpening of pi). I tuned them on 60000..61023 only:

| setting (choose set, TEST - BASE, n=1024) | mean_diff | se |
|---|---|---|
| qs=0 (arb point + retail mean only, 2 lookups) | +0.4907 | 0.0164 |
| qs=0.5 | +0.5430 | 0.0175 |
| **qs=1, qk=1 (chosen = file defaults)** | **+0.5667** | 0.0182 |
| qs=1.5 | +0.5405 | 0.0191 |
| qs=2 | +0.4286 | 0.0196 |
| qk=0.5 | +0.2278 | 0.0197 |
| qk=2 | +0.4286 | 0.0195 |

## Results, TEST minus BASE (paired, same seeds; se = sample sd / sqrt(n))

| set | n | mean_diff | se | wins | worst5 (seed, diff) |
|---|---|---|---|---|---|
| choose 60000..61023 | 1024 | +0.5667 | 0.0182 | 876 | (60773, -2.1817), (60225, -1.7215), (60210, -1.5763), (60204, -1.3188), (60084, -1.2056) |
| **report 10000..12047** | 2048 | **+0.5662** | **0.0125** | 1746 | (10212, -1.5005), (11022, -1.2620), (11179, -1.2197), (11074, -1.1872), (12034, -1.1469) |

The report seeds were scored once per arm, with no tuning on them. Means on the report set: BASE 532.3857, TEST 532.9519, true flag 536.1083. Against Taylor's stored HA rows: BASE +2.828 se 0.041, **TEST +3.394 se 0.040**, true flag +6.551 se 0.046.

**True-flag gap.** BASE with the true arb/retail flag (`oraclearb=1`, plain tables, soft=1) minus BASE: n=2048, +3.7226, se 0.0237, wins 2044. TEST closes 0.566 of that 3.72, about 15%. (The README's "worth 2.7" and +6.985 refer to the qm tables. With plain tables the flag is worth 3.72 on this replica.)

## Confusion matrix (report seeds, first own trade of each step, all 2048 seeds summed)

| arm | arb called arb | arb called retail | retail called arb | retail called retail | misclassified |
|---|---|---|---|---|---|
| BASE | 6,301,005 | 1,309,754 | 1,048,181 | 5,566,964 | 16.57% |
| TEST | 6,394,351 | 1,244,969 | 1,128,646 | 5,470,053 | 16.67% |
| true flag | 7,385,296 | 0 | 0 | 6,708,520 | 0.00% |

The classifier itself is unchanged, so the error rate stays at about 16.6%. The gain comes from making each error cheaper, not from fewer errors. That is consistent with the README's finding that these features cap the classifier near 16%. (A trade counts as called arb when the policy ended the callback in phase 0. The truth is the replica's `arb` argument.)

## Cost

Seconds per seed (one core, replica): BASE 0.098, TEST 0.157, flag 0.088. The choose set ran in 12 s (BASE) and 21 s (TEST) on 8 processes. Each callback needs up to 4 table lookups instead of 1 (2 with qs=0, which keeps +0.49 of the +0.57). That matters for the 245k gas budget in the Solidity port.

## Verdict: worth porting

+0.566 with se 0.0125 on the untouched report set clears the bar (> +0.3, se < 0.1). The choose and report sets agree (+0.567 / +0.566), TEST wins 85% of report seeds, and the worst seed loses only 1.5. It is about 10x the gain of either earlier item. Caveats:
- The extra lookups cost gas. qs=0 is a cheaper fallback (2 lookups, +0.49 on the choose set; not scored on the report set).
- I tested it only with the plain tables. Whether it adds on top of the `qm` tables or the lattice filter (which already use sd) is untested.
- Replica gaps have so far tracked the engine, per the README, but nothing here ran on the Rust engine.

## The complete TEST file (`realpq.py`, a copy of `real.py` with only this idea added)

sha256 `488ad8a89884ed22713ebc38031cbc71ea98dc78e759e65a790bc478247c32e1`. It goes in `research/solved/` next to `real.py`. Usage in my harness: `realpq.RealPolicy(T, {'soft': 1.0})` (defaults qs=1, qk=1). Its `__main__` works like `real.py`'s, e.g. `python3 -I realpq.py pol_mid.json plain:soft=1`.

```python
"""[posterior-quote variant: quotes are chosen against the step's whole posterior
(arb point mass + retail truncated-Gaussian quadrature) instead of its collapsed mean.]
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
        self.mix = None               # posterior of this step as [(weight, ln p)]
        # diagnostics
        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 _ab(self, w, w2, ph):
        A, B = self.tabs[ph]
        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 ph == 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
        a = (A[i][j] * (1 - f) + A[i + 1][j] * f) * (1 - g) + (A[i][j + 1] * (1 - f) + A[i + 1][j + 1] * f) * g
        b = (B[i][j] * (1 - f) + B[i + 1][j] * f) * (1 - g) + (B[i][j + 1] * (1 - f) + B[i + 1][j + 1] * f) * g
        return a, b
    def _retail_nodes(self, wt, mt, vt):
        qs = self.o.get('qs', 1.0)
        if vt <= 0.0 or qs <= 0.0 or wt <= 0.0:
            return [(wt, mt)]
        h = qs * math.sqrt(3.0 * vt)          # 3-point Gauss-Hermite
        return [(wt / 6.0, mt - h), (wt * 2.0 / 3.0, mt), (wt / 6.0, mt + h)]
    def lookup(self, w, w2, phs, sd):
        Q = self.Q
        if Q is None:
            mix = self.mix
            if not mix:
                return self._fees(w, w2, phs)
            ls = math.log(self.sp); l2 = math.log(self.s2h)
            LA = LB = 0.0
            for wt, x in mix:
                a, b = self._ab((ls - x) * 1e4, (l2 - x) * 1e4, phs)
                LA += wt * (x + a * 1e-4); LB += wt * (x - b * 1e-4)
            fa = 1.0 - math.exp(-(LA - ls)) if LA > ls else 0.0
            fb = 1.0 - math.exp(-(ls - LB)) if LB < ls else 0.0
            return fb, fa
        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
            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
                k = o.get('qk', 1.0)
                qa = pi ** k / (pi ** k + (1.0 - pi) ** k)
                self.mix = [(qa, lp_arb)] + self._retail_nodes(1.0 - qa, mt, vt)
            elif is_arb:
                self.m, self.v = lp_arb, 0.0
                self.mix = [(1.0, lp_arb)]
            else:
                self.m, self.v = mt, vt
                self.mix = self._retail_nodes(1.0, 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))
    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)
```

Harness, not part of the policy: each arm ran `replica2.run(seed, policy)` over a pool of 8 processes. A subclass wrapper counted the confusion cells on callbacks with `t - tlast >= 1`, and pairing was done per seed.