1"""Run-level outcome metrics. One dict per run; every number is traceable to ticks/trades/agents/llm frames.23Primary outcomes (pre-registered in docs/identification.md):4 efficiency : msd_log (mean squared log deviation of mid from the latent value), mad_pct, time_to_absorb,5 post-shock response ratio, variance ratio VR(5) of log mid changes.6 liquidity : quoted spread (bp), effective spread (bp), depth (lots within 5 levels), Kyle lambda7 (bp per 100 lots signed volume), fill rate, resiliency (half-life of spread after a widening).8 instability : bubble amplitude/duration (P/V - 1 above `bubble_thr` for >= `bubble_min_len` ticks),9 max drawdown, crash count (drawdown >= crash_thr within crash_window), volatility clustering10 (sum of ACF(|r|) lags 1..10), excess kurtosis of tick returns, halts.11 contagion : (multi-asset) excess return correlation over fundamental correlation, Diebold-Yilmaz total12 spillover from a VAR(2) FEVD, lead-lag cross-correlation.13Secondary: volume, turnover, cancel ratio, P&L by policy, Gini, bankruptcies, MM inventory, LLM validity/cost.14"""15from __future__ import annotations1617import numpy as np18import pandas as pd192021def _acf(x: np.ndarray, lags: int) -> np.ndarray:22 x = np.asarray(x, dtype=float)23 x = x - x.mean()24 n = len(x)25 denom = (x ** 2).sum()26 if denom <= 0 or n <= lags + 1:27 return np.full(lags, np.nan)28 return np.array([(x[:-k] * x[k:]).sum() / denom for k in range(1, lags + 1)])293031def variance_ratio(logp: np.ndarray, q: int = 5) -> float:32 r = np.diff(logp)33 if len(r) < 5 * q or r.var() <= 0:34 return np.nan35 rq = logp[q:] - logp[:-q]36 return float(rq.var() / (q * r.var()))373839def gini(w: np.ndarray) -> float:40 w = np.sort(np.clip(np.asarray(w, dtype=float), 0, None))41 n = len(w)42 if n == 0 or w.sum() <= 0:43 return np.nan44 return float((2 * np.arange(1, n + 1) - n - 1).dot(w) / (n * w.sum()))454647def bubbles(ratio: np.ndarray, thr: float, min_len: int) -> tuple[int, float, float]:48 """Episodes where ratio = P/V - 1 stays above thr (or below -thr) for >= min_len consecutive ticks."""49 above = np.abs(ratio) > thr50 n, amp, dur = 0, 0.0, 0.051 i = 052 while i < len(above):53 if above[i]:54 j = i55 while j < len(above) and above[j]:56 j += 157 if j - i >= min_len:58 n += 159 amp = max(amp, float(np.abs(ratio[i:j]).max()))60 dur += j - i61 i = j62 else:63 i += 164 return n, amp, dur656667def drawdowns(p: np.ndarray, crash_thr: float, window: int) -> tuple[float, int]:68 p = np.asarray(p, dtype=float)69 run_max = np.maximum.accumulate(p)70 dd = 1.0 - p / run_max71 max_dd = float(dd.max()) if len(dd) else np.nan72 # crash: drop of >= crash_thr from the max of the trailing `window` ticks73 n_crash = 074 i = window75 while i < len(p):76 local_max = p[i - window:i].max()77 if p[i] <= (1 - crash_thr) * local_max:78 n_crash += 179 i += window # one crash per window80 else:81 i += 182 return max_dd, n_crash838485def kyle_lambda(ticks: pd.DataFrame) -> float:86 """OLS of mid change (bp) on signed volume (per 100 lots), tick by tick."""87 d = ticks.dropna(subset=["mid"])88 if len(d) < 30:89 return np.nan90 y = 1e4 * np.diff(np.log(d["mid"].to_numpy()))91 x = d["signed_volume"].to_numpy()[1:] / 100.092 if x.var() <= 0:93 return np.nan94 return float(np.cov(x, y)[0, 1] / x.var())959697def time_to_absorb(ticks: pd.DataFrame, shock_t: int, tol: float = 0.01, hold: int = 5, max_h: int = 100) -> float:98 """First horizon h after shock_t such that |P/V - 1| < tol for `hold` consecutive ticks."""99 d = ticks[(ticks.t >= shock_t) & (ticks.t < shock_t + max_h)]100 dev = np.abs(d["mid"].to_numpy() / d["full_value"].to_numpy() - 1.0)101 ok = dev < tol102 for h in range(len(ok) - hold + 1):103 if ok[h:h + hold].all():104 return float(h)105 return float(max_h)106107108def response_ratio(ticks: pd.DataFrame, shock_t: int, h: int = 20) -> float:109 """(P_{t+h} - P_{t-1}) / (V_{t+h} - V_{t-1}) : 1 = full absorption; <1 under-reaction; >1 over-reaction."""110 d = ticks.set_index("t")111 try:112 dp = d.loc[shock_t + h, "mid"] - d.loc[shock_t - 1, "mid"]113 dv = d.loc[shock_t + h, "full_value"] - d.loc[shock_t - 1, "full_value"]114 except KeyError:115 return np.nan116 return float(dp / dv) if abs(dv) > 1e-9 else np.nan117118119def resiliency(spread: np.ndarray) -> float:120 """Half-life (ticks) of spread deviations from their median: from AR(1) of log spread."""121 s = np.log(np.clip(spread[~np.isnan(spread)], 1e-6, None))122 if len(s) < 30:123 return np.nan124 s = s - np.median(s)125 phi = (s[:-1] * s[1:]).sum() / max((s[:-1] ** 2).sum(), 1e-12)126 phi = min(max(phi, 1e-6), 0.999)127 return float(np.log(0.5) / np.log(phi))128129130def series_outcomes(ticks: pd.DataFrame, warmup: int = 50, bubble_thr: float = 0.05, bubble_min_len: int = 10,131 crash_thr: float = 0.10, crash_window: int = 20, shocks: list[dict] | None = None) -> dict:132 out: dict = {}133 for a, d in ticks.groupby("asset"):134 d = d[d.t >= warmup].sort_values("t")135 P, V = d["mid"].to_numpy(), d["full_value"].to_numpy()136 logp = np.log(P)137 r = np.diff(logp)138 sfx = "" if a == 0 else f"_a{a}"139 dev = logp - np.log(V)140 out["msd_log" + sfx] = float(np.mean(dev ** 2))141 out["mad_pct" + sfx] = float(100 * np.mean(np.abs(dev)))142 out["bias_pct" + sfx] = float(100 * np.mean(dev))143 out["vr5" + sfx] = variance_ratio(logp, 5)144 out["vr10" + sfx] = variance_ratio(logp, 10)145 out["vol_bp" + sfx] = float(1e4 * r.std()) if len(r) > 2 else np.nan146 out["kurtosis" + sfx] = float(pd.Series(r).kurt()) if len(r) > 10 else np.nan147 out["acf_absr_1_10" + sfx] = float(np.nansum(_acf(np.abs(r), 10))) if len(r) > 30 else np.nan148 out["acf_r_1" + sfx] = float(_acf(r, 1)[0]) if len(r) > 30 else np.nan149 with np.errstate(all="ignore"):150 import warnings151 warnings.simplefilter("ignore", RuntimeWarning)152 out["spread_bp" + sfx] = float(np.nanmean(1e4 * d["spread"] / d["mid"])) if d["spread"].notna().any() else np.nan153 out["eff_spread_bp" + sfx] = float(1e4 * np.nanmean(d["eff_spread"] / d["mid"]))154 out["depth5" + sfx] = float(np.mean(d["depth_bid5"] + d["depth_ask5"]))155 out["lambda_bp_100" + sfx] = kyle_lambda(d)156 out["resiliency" + sfx] = resiliency(d["spread"].to_numpy())157 nb, amp, dur = bubbles(P / V - 1.0, bubble_thr, bubble_min_len)158 out["n_bubbles" + sfx], out["bubble_amp" + sfx], out["bubble_dur" + sfx] = nb, amp, dur159 out["max_dd" + sfx], out["n_crashes" + sfx] = drawdowns(P, crash_thr, crash_window)160 out["volume" + sfx] = float(d["volume"].sum())161 out["n_halts" + sfx] = int(d["n_halts"].max()) if "n_halts" in d else 0162 out["share_ticks_halted" + sfx] = float(d["halted"].mean()) if "halted" in d else 0.0163 if shocks:164 tta, rr = [], []165 for s in shocks:166 if s.get("shock") in ("fundamental_jump", "news_burst") and s["t"] >= warmup:167 tta.append(time_to_absorb(ticks[ticks.asset == a], s["t"]))168 rr.append(response_ratio(ticks[ticks.asset == a], s["t"]))169 if tta:170 out["time_to_absorb" + sfx] = float(np.mean(tta))171 out["response_ratio" + sfx] = float(np.nanmean(rr))172 assets = sorted(ticks["asset"].unique())173 if len(assets) >= 2:174 piv = ticks[ticks.t >= warmup].pivot(index="t", columns="asset", values="mid")175 fv = ticks[ticks.t >= warmup].pivot(index="t", columns="asset", values="full_value")176 rp, rv = np.log(piv).diff().dropna(), np.log(fv).diff().dropna()177 cp, cv = rp.corr().to_numpy(), rv.corr().to_numpy()178 iu = np.triu_indices(len(assets), 1)179 out["ret_corr"] = float(np.nanmean(cp[iu]))180 out["fund_corr"] = float(np.nanmean(cv[iu]))181 out["excess_comovement"] = out["ret_corr"] - out["fund_corr"]182 out["spillover_dy"] = diebold_yilmaz(rp.to_numpy())183 out["leadlag_01"] = float(np.corrcoef(rp.iloc[:-1, 0], rp.iloc[1:, 1])[0, 1]) if len(rp) > 10 else np.nan184 return out185186187def diebold_yilmaz(R: np.ndarray, p: int = 2, h: int = 10) -> float:188 """Total spillover index from the generalized FEVD of a VAR(p) (Diebold–Yilmaz 2012)."""189 T, K = R.shape190 if T < 10 * K * p + 10:191 return np.nan192 X = np.hstack([R[p - j - 1:T - j - 1] for j in range(p)])193 Y = R[p:]194 X = np.hstack([np.ones((len(Y), 1)), X])195 B, *_ = np.linalg.lstsq(X, Y, rcond=None)196 E = Y - X @ B197 Sigma = np.cov(E.T)198 A = [np.eye(K)]199 coefs = [B[1 + K * j:1 + K * (j + 1)].T for j in range(p)]200 for i in range(1, h):201 A.append(sum(coefs[j] @ A[i - j - 1] for j in range(p) if i - j - 1 >= 0))202 sig = np.diag(Sigma)203 theta = np.zeros((K, K))204 for i in range(K):205 for j in range(K):206 num = sum((A[l][i] @ Sigma[:, j]) ** 2 for l in range(h)) / sig[j]207 den = sum(A[l][i] @ Sigma @ A[l][i] for l in range(h))208 theta[i, j] = num / den209 theta = theta / theta.sum(axis=1, keepdims=True)210 return float(100 * (theta.sum() - np.trace(theta)) / K)211212213def agent_outcomes(agents: pd.DataFrame, trades: pd.DataFrame, llm: pd.DataFrame, n_ticks: int) -> dict:214 out: dict = {}215 g = agents.groupby("policy")["pnl_pct"].mean()216 for k, v in g.items():217 out[f"pnl_{k}"] = float(v)218 out["gini_final"] = gini(agents["final_wealth"].to_numpy())219 out["gini_initial"] = gini(agents["initial_wealth"].to_numpy())220 out["n_bankrupt"] = int(agents["frozen"].sum())221 out["cancel_ratio"] = float(agents["n_cancels"].sum() / max(agents["n_orders"].sum(), 1))222 out["reject_ratio"] = float(agents["n_rejects"].sum() / max(agents["n_orders"].sum(), 1))223 out["fill_ratio"] = float(agents["n_fills"].sum() / max(agents["n_orders"].sum(), 1))224 out["turnover"] = float(trades["qty"].sum() / max(agents["position"].abs().sum() + 1, 1)) if len(trades) else 0.0225 mm = agents[agents.policy == "market_maker"]226 if len(mm):227 out["mm_pnl"] = float(mm["pnl_pct"].mean())228 out["mm_inventory_abs"] = float(mm["position"].abs().mean())229 if len(llm):230 out["llm_valid_rate"] = float(llm["valid"].mean())231 out["llm_hold_rate"] = float((llm["action"] == "hold").mean())232 out["llm_market_share"] = float((llm["order_type"] == "market").mean())233 out["llm_usd"] = float(llm["usd"].sum())234 out["llm_tokens"] = float(llm["input_tokens"].sum() + llm["output_tokens"].sum())235 out["llm_latency_s"] = float(llm["latency_s"].mean())236 out["llm_cached_hits"] = float(llm["cached_hit"].mean())237 # synchrony and order-flow correlation among LLM agents: decisions grouped by decision round238 d = llm[llm.valid & (llm.action == "place")].copy()239 if len(d) > 5:240 d["sign"] = np.where(d["side"] == "buy", 1, -1)241 rounds = d.groupby("t")["sign"]242 sync = rounds.apply(lambda s: abs(s.mean()) if len(s) > 1 else np.nan).dropna()243 out["llm_synchrony"] = float(sync.mean()) if len(sync) else np.nan244 out["llm_buy_share"] = float((d["sign"] > 0).mean())245 out["llm_pairwise_corr"] = pairwise_sign_corr(d)246 return out247248249def pairwise_sign_corr(d: pd.DataFrame, window: int = 20) -> float:250 """Average pairwise correlation of agents' net order sign aggregated per `window` ticks."""251 d = d.copy()252 d["w"] = d["t"] // window253 piv = d.pivot_table(index="w", columns="agent_id", values="sign", aggfunc="sum").fillna(0.0)254 if piv.shape[1] < 2 or piv.shape[0] < 3:255 return np.nan256 c = piv.corr().to_numpy()257 iu = np.triu_indices(piv.shape[1], 1)258 vals = c[iu]259 return float(np.nanmean(vals)) if np.isfinite(vals).any() else np.nan260261262def flow_corr_by_policy(trades: pd.DataFrame, agents: pd.DataFrame, window: int = 20) -> dict:263 """Average pairwise correlation of agents' net signed traded volume per `window` ticks, by policy group264 (a classical analogue of llm_pairwise_corr that needs no decision log)."""265 if len(trades) < 50:266 return {}267 t = trades.copy()268 t["w"] = (t["time"] // window).astype(int)269 buy = t.groupby(["w", "buyer"])["qty"].sum().rename("q").reset_index().rename(columns={"buyer": "agent_id"})270 sell = t.groupby(["w", "seller"])["qty"].sum().rename("q").reset_index().rename(columns={"seller": "agent_id"})271 sell["q"] = -sell["q"]272 flow = pd.concat([buy, sell]).groupby(["w", "agent_id"])["q"].sum().reset_index()273 piv = flow.pivot(index="w", columns="agent_id", values="q").fillna(0.0)274 pol = agents.set_index("agent_id")["policy"]275 out = {}276 for name, grp in pol.groupby(pol):277 ids = [i for i in grp.index if i in piv.columns]278 if len(ids) < 3 or piv.shape[0] < 4:279 continue280 c = piv[ids].corr().to_numpy()281 iu = np.triu_indices(len(ids), 1)282 vals = c[iu]283 if np.isfinite(vals).any():284 out[f"flow_corr_{name}"] = float(np.nanmean(vals))285 return out286287288def run_outcomes(res, bubble_thr: float = 0.05, bubble_min_len: int = 10, crash_thr: float = 0.10, crash_window: int = 20) -> dict:289 cfg = res.config290 out = {"name": cfg.name, "seed": cfg.seed, "config_hash": cfg.config_hash(), "n_traders": cfg.population.n_traders,291 "llm_share": cfg.population.llm_share, **{f"tr_{k}": v for k, v in cfg.treatment.items()}}292 out.update({k: res.manifest[k] for k in ("llm_share_traders", "llm_share_orders", "llm_share_volume", "n_llm_decisions", "n_trades", "wall_s")})293 out.update(series_outcomes(res.ticks, cfg.warmup, bubble_thr, bubble_min_len, crash_thr, crash_window, res.manifest.get("shocks")))294 out.update(agent_outcomes(res.agents, res.trades, res.llm, cfg.n_ticks))295 out.update(flow_corr_by_policy(res.trades, res.agents))296 return out297