1#!/usr/bin/env python2"""Numerical verification of the closed forms of Propositions 1-5 (paper/sections/theory.tex).34Monte Carlo with N = 50 traders and R = 200,000 draws of the one-period model5 s_i = v + eta_{g(i)} + u_i, x_i = beta s_i, Y = sum_i x_i + z, P = lambda Y.6Every row of results/tables/theory_checks.csv compares a closed form with its simulated7counterpart and reports the relative error; the script exits with a non-zero code if any8relative error exceeds 2 %.910Run: .venv/bin/python scripts/90_theory_checks.py11"""12from __future__ import annotations1314import math15import sys16from pathlib import Path1718import numpy as np19import pandas as pd2021ROOT = Path(__file__).resolve().parents[1]22OUT = ROOT / "results" / "tables" / "theory_checks.csv"2324# ----------------------------------------------------------------------------- primitives25N = 50 # traders26R = 200_000 # Monte Carlo draws27SV, SE, SU, SZ = 1.0, 1.0, 0.8, 5.0 # sigma_v, sigma_eta, sigma_u, sigma_z28BETA = 1.029TOL = 0.0230rng = np.random.default_rng(20261027)3132SS2 = SV**2 + SE**2 + SU**2 # signal variance sigma_s^233RHO0 = SV**2 / SS2 # different / private models34RHO1 = (SV**2 + SE**2) / SS2 # same model3536rows: list[dict] = []373839def add(prop: str, case: str, quantity: str, theory: float, sim: float) -> None:40 rel = abs(sim - theory) / abs(theory) if theory != 0 else abs(sim - theory)41 rows.append(dict(proposition=prop, case=case, quantity=quantity, theory=theory,42 simulation=sim, rel_error=rel))434445# ----------------------------------------------------------------------------- closed forms46def omega(H: float, Phi: float) -> float:47 return SE**2 * (H + (1.0 - Phi) / N) + SU**2 / N + SZ**2 / (BETA**2 * N**2)484950def closed_forms(sizes: list[int]) -> dict:51 Phi = sum(sizes) / N52 H = sum((n / N) ** 2 for n in sizes)53 Om = omega(H, Phi)54 tau = SV**2 / (SV**2 + Om)55 lam = tau / (BETA * N)56 w = (N * H - Phi) / (N - 1) # share of same-model ordered pairs57 rhobar = RHO0 + w * SE**2 / SS258 sync = 0.5 + ((1 - w) * math.asin(RHO0) + w * math.asin(RHO1)) / math.pi59 return dict(H=H, Phi=Phi, Omega=Om, tau=tau, lam=lam, w=w, rhobar=rhobar, sync=sync,60 mse=SV**2 * (1 - tau))616263# ----------------------------------------------------------------------------- simulation64def simulate(sizes: list[int]) -> tuple[np.ndarray, np.ndarray, np.ndarray]:65 """Draw (v, x, Y) for a population with shared-model cohorts of the given sizes66 (first cohort occupies traders 0..n_1-1, etc.); remaining traders use private models."""67 v = rng.normal(0.0, SV, R)68 eps = rng.normal(0.0, SU, (R, N))69 start = 070 for n in sizes:71 eps[:, start:start + n] += rng.normal(0.0, SE, (R, 1))72 start += n73 if start < N:74 eps[:, start:] += rng.normal(0.0, SE, (R, N - start))75 x = BETA * (v[:, None] + eps)76 Y = x.sum(axis=1) + rng.normal(0.0, SZ, R)77 return v, x, Y787980def flow_stats(v: np.ndarray, x: np.ndarray, Y: np.ndarray) -> dict:81 cov = np.cov(v, Y, ddof=0)82 lam_hat = cov[0, 1] / cov[1, 1]83 tau_hat = cov[0, 1] ** 2 / (cov[0, 0] * cov[1, 1])84 P = lam_hat * Y85 mse = np.var(P - v)86 C = np.corrcoef(x, rowvar=False)87 off = ~np.eye(N, dtype=bool)88 rhobar = C[off].mean()89 sg = np.sign(x)90 A = (sg.T @ sg) / R # E[sign_i sign_j]91 sync = ((1.0 + A[off]) / 2.0).mean()92 return dict(tau=tau_hat, lam=lam_hat, mse=mse, rhobar=rhobar, sync=sync)939495def check_population(prop: str, label: str, sizes: list[int], what=("tau", "lam", "mse", "rhobar", "sync")):96 th = closed_forms(sizes)97 v, x, Y = simulate(sizes)98 sim = flow_stats(v, x, Y)99 for q in what:100 add(prop, label, q, th[q], sim[q])101 return th, sim, (v, x, Y)102103104# ============================================================================= Props 1-2105print(f"Monte Carlo: N={N}, R={R:,}, sigma_v={SV}, sigma_eta={SE}, sigma_u={SU}, sigma_z={SZ}, beta={BETA}")106single_cases = {"phi=0": [], "phi=0.1": [5], "phi=0.2": [10], "phi=0.5": [25], "phi=1": [50]}107taus = {}108keep = None109for label, sizes in single_cases.items():110 th, sim, draws = check_population("P1-P2", f"single model, {label}", sizes)111 taus[label] = (th["tau"], sim["tau"])112 if label == "phi=0.2":113 keep = draws114115# Prop 1(i): monotonicity across the grid (theory) and in simulation116tau_th = [taus[k][0] for k in single_cases]117tau_mc = [taus[k][1] for k in single_cases]118add("P1(i)", "tau strictly decreasing over phi grid (1=yes)", "monotone",119 1.0, float(all(np.diff(tau_th) < 0) and all(np.diff(tau_mc) < 0)))120# Prop 1(i): phi=0 and phi=1/N coincide in Omega121add("P1(i)", "Omega(phi=1/N) / Omega(phi=0)", "ratio", 1.0, omega((1 / N) ** 2, 1 / N) / omega(0.0, 0.0))122# Prop 1(ii): uniform bound tau(phi) <= tau_inf(phi) and N->infinity limit (theory, large N)123for phi in (0.2, 0.5):124 tau_inf = SV**2 / (SV**2 + phi**2 * SE**2)125 OmN = lambda n: SE**2 * (phi**2 + (1 - phi) / n) + SU**2 / n + SZ**2 / (BETA**2 * n**2)126 tau_big = SV**2 / (SV**2 + OmN(10**7))127 add("P1(ii)", f"N->inf limit, phi={phi}", "tau_inf", tau_inf, tau_big)128 add("P1(ii)", f"bound tau(N=50)<=tau_inf, phi={phi} (1=yes)", "bound", 1.0,129 float(SV**2 / (SV**2 + OmN(N)) <= tau_inf))130add("P1(iii)", "phi=0, N=1e7", "tau->1", 1.0,131 SV**2 / (SV**2 + (SE**2 + SU**2) / 1e7 + SZ**2 / (BETA**2 * 1e14)))132133# Prop 2(ii): same-model vs cross-model pair correlations in the phi=0.5 population134v, x, Y = simulate([25])135C = np.corrcoef(x, rowvar=False)136same = C[:25, :25][~np.eye(25, dtype=bool)].mean()137cross = C[:25, 25:].mean()138add("P2(ii)", "phi=0.5", "rho_same", RHO1, same)139add("P2(ii)", "phi=0.5", "rho_cross", RHO0, cross)140add("P2(iv)", "phi=0.5, same-model pair", "S(rho1)=1/2+asin(rho1)/pi", 0.5 + math.asin(RHO1) / math.pi,141 ((1 + (np.sign(x[:, :25]).T @ np.sign(x[:, :25]) / R)[~np.eye(25, dtype=bool)]) / 2).mean())142# Prop 2 placebo: sigma_eta = 0 -> no dependence on phi (theory)143add("P2 placebo", "sigma_eta=0: Omega(phi=1)/Omega(phi=0)", "ratio", 1.0,144 (SU**2 / N + SZ**2 / (BETA**2 * N**2)) / (SU**2 / N + SZ**2 / (BETA**2 * N**2)))145del v, x, Y, C146147# ============================================================================= Prop 3148div_cases = {149 "Phi=1, G=1 (50)": [50],150 "Phi=1, G=2 uniform (25,25)": [25, 25],151 "Phi=1, G=5 uniform (10x5)": [10] * 5,152 "Phi=1, G=3 skewed (30,10,10)": [30, 10, 10],153 "Phi=0.6, G=1 (30)": [30],154 "Phi=0.6, G=2 uniform (15,15)": [15, 15],155 "Phi=0.6, G=3 uniform (10,10,10)": [10, 10, 10],156}157div_tau = {}158for label, sizes in div_cases.items():159 th, sim, _ = check_population("P3", label, sizes, what=("tau", "rhobar", "sync"))160 div_tau[label] = (th["H"], th["tau"], sim["tau"])161# simulated tau is decreasing in H within each Phi group (hence single <= mixed <= uniform)162ok = True163for Phi_label in ("Phi=1", "Phi=0.6"):164 grp = sorted((H, tau_s) for lbl, (H, _, tau_s) in div_tau.items() if lbl.startswith(Phi_label + ","))165 ok = ok and all(grp[k][1] > grp[k + 1][1] for k in range(len(grp) - 1)) # H ascending -> tau descending166add("P3(iii)", "simulated tau decreasing in H within each Phi (1=yes)", "ordering", 1.0, float(ok))167# H bounds168for label, sizes in div_cases.items():169 H, Phi, G = sum((n / N) ** 2 for n in sizes), sum(sizes) / N, len(sizes)170 add("P3(ii)", f"Phi^2/G <= H <= Phi^2 for {label} (1=yes)", "bounds", 1.0,171 float(Phi**2 / G - 1e-12 <= H <= Phi**2 + 1e-12))172173# ============================================================================= Prop 4174v, x, Y = keep175th = closed_forms([10])176lam = th["lam"]177for lg in (0.3, 0.6):178 gamma = lg / lam179 dP = lam * Y.copy() # start of the iterated response180 for _ in range(400):181 dP = lam * (Y + gamma * dP)182 amp = np.var(dP) / np.var(lam * Y)183 add("P4(i)", f"lambda*gamma={lg}", "variance amplification 1/(1-lg)^2", 1 / (1 - lg) ** 2, amp)184 kappa = 1 / (1 - lg)185 add("P4(iii)", f"lambda*gamma={lg}", "Var(P1-v)=sv^2[(1-tau)+tau(kappa-1)^2]",186 SV**2 * ((1 - th["tau"]) + th["tau"] * (kappa - 1) ** 2), np.var(dP - v))187 # fixed point = kappa*lambda*Y188 add("P4(i)", f"lambda*gamma={lg}", "max|dP - kappa*lam*Y| (abs)", 0.0,189 float(np.max(np.abs(dP - kappa * lam * Y))))190# divergence when lambda*gamma >= 1191for lg in (1.0, 1.2):192 gamma = lg / lam193 dP = lam * Y[:1000].copy()194 for _ in range(60):195 dP = lam * (Y[:1000] + gamma * dP)196 # lambda*gamma = 1 diverges linearly (|dP| = 61|lam Y| after 60 rounds), > 1 geometrically197 add("P4(ii)", f"lambda*gamma={lg}: |dP| after 60 rounds >= 50 * |lam*Y| (1=yes)", "divergence",198 1.0, float(np.all(np.abs(dP) >= 50 * np.abs(lam * Y[:1000]))))199del v, x, Y200201# ============================================================================= Prop 5202# (ii)-(iii): cascade fixed point and multiplier203def cascade(a: float, dP0: float, b: float, X: float, lam_: float, m: float, rounds: int = 2000) -> float:204 dP = dP0205 for _ in range(rounds):206 F = max(0.0, (-X * dP - b) / m)207 dP = dP0 - lam_ * F208 if dP < -1e6:209 return dP210 return dP211212213lam_, m, X = 0.02, 1.0, 30.0 # a = 0.6214a = lam_ * X / m215for b in (0.0, 5.0):216 dP0 = -1.0217 fp = cascade(a, dP0, b, X, lam_, m)218 add("P5(ii)", f"a={a:.2f}, b={b}, dP0={dP0}", "fixed point (dP0+lam*b/m)/(1-a)",219 (dP0 + lam_ * b / m) / (1 - a), fp)220add("P5(ii)", f"a={a:.2f}, b=0", "multiplier 1/(1-a)", 1 / (1 - a), cascade(a, -1.0, 0.0, X, lam_, m) / -1.0)221add("P5(i)", f"a={a:.2f}, b=5, dP0=-0.1 (inside buffer)", "dP = dP0", -0.1, cascade(a, -0.1, 5.0, X, lam_, m))222for Xs in (50.0, 60.0): # a = 1.0, 1.2223 fp = cascade(lam_ * Xs / m, -1.0, 0.0, Xs, lam_, m, rounds=2000) # a=1: linear, -2001 after 2000 rounds224 add("P5(iii)", f"a={lam_ * Xs / m:.2f}: cascade below -1e3 (1=yes)", "no finite fixed point", 1.0,225 float(fp < -1e3))226227# (iv): distribution of the cohort net position X = |sum_{i in L} x_i|, N_L = 10 leveraged traders228NL = 10229c_thr = None230for phiL, sizes in (("phi_L=0", []), ("phi_L=0.5", [5]), ("phi_L=1", [10])):231 v = rng.normal(0.0, SV, R)232 eps = rng.normal(0.0, SU, (R, NL))233 if sizes:234 eps[:, :sizes[0]] += rng.normal(0.0, SE, (R, 1))235 eps[:, sizes[0]:] += rng.normal(0.0, SE, (R, NL - sizes[0]))236 else:237 eps += rng.normal(0.0, SE, (R, NL))238 xL = BETA * (v[:, None] + eps)239 Xnet = np.abs(xL.sum(axis=1))240 nL = sizes[0] if sizes else 0241 wL = nL * (nL - 1) / (NL * (NL - 1))242 rhoL = RHO0 + wL * SE**2 / SS2243 sigX = BETA * math.sqrt(SS2) * math.sqrt(NL * (1 + (NL - 1) * rhoL))244 if c_thr is None:245 c_thr = 1.0 * sigX # fixed spiral threshold m/lambda across cases246 add("P5(iv)", phiL, "sd(X_net)", sigX, float(np.std(xL.sum(axis=1))))247 add("P5(iv)", phiL, "E|X| = sigX*sqrt(2/pi)", sigX * math.sqrt(2 / math.pi), float(Xnet.mean()))248 p_th = 2 * (1 - 0.5 * (1 + math.erf(c_thr / sigX / math.sqrt(2))))249 add("P5(iv)", f"{phiL}, threshold m/lambda={c_thr:.3f}", "P(X>=m/lambda)=2*Phibar", p_th,250 float((Xnet >= c_thr).mean()))251 del v, eps, xL, Xnet252253# ============================================================================= report254df = pd.DataFrame(rows)255OUT.parent.mkdir(parents=True, exist_ok=True)256df.to_csv(OUT, index=False)257pd.set_option("display.width", 200, "display.max_rows", 500, "display.max_colwidth", 60)258print(df.to_string(index=False, float_format=lambda f: f"{f:.5f}"))259worst = df["rel_error"].max()260print(f"\n{len(df)} checks written to {OUT.relative_to(ROOT)}; max relative error = {worst:.4%}")261if worst >= TOL:262 print(f"FAIL: at least one relative error >= {TOL:.0%}")263 sys.exit(1)264print(f"OK: all relative errors < {TOL:.0%}")265