1"""Linear-quadratic delegation game (Model A, §1.8 of docs/model_specification.md).23Closed forms derived in scripts/00_lq_derivations.py (sympy) and re-implemented here in numpy so that every formula is4cross-checked by Monte Carlo in tests/unit/test_lq_game.py.56Conventions7-----------8theta ~ N(0, 1/tau_theta); human signal s_i = theta + e_i (precision tau_H); model-m signal s_m = theta + eta_m9(precision tau_A), common to all users of model m; M models, delegation share D split evenly.10Best response x_i = E_i[kappa0 * theta + alpha * X]; alpha < 0 strategic substitutes.11Payoff normalised by |U_xx| = 1; the social curvature on the aggregate action is (1 - alpha) under the CRS restriction12U_xX + U_XX = 0 (Lemma 2); the general case U_XX = -upsilon * U_xX is kept for Proposition 4(a).13"""14from __future__ import annotations1516from dataclasses import dataclass, replace1718import numpy as np192021@dataclass(frozen=True)22class LQParams:23 kappa0: float = 1.024 alpha: float = -0.325 tau_theta: float = 1.026 tau_H: float = 1.027 tau_A: float = 4.028 M: int = 129 upsilon: float = 1.0 # U_XX = -upsilon U_xX ; 1 = CRS/pecuniary restriction30 U_Xtheta: float = 0.0 # direct aggregate-state interaction (0 under CRS)3132 @property33 def beta_H(self):34 return self.tau_H / (self.tau_theta + self.tau_H)3536 @property37 def beta_A(self):38 return self.tau_A / (self.tau_theta + self.tau_A)394041# --------------------------------------------------------------------------- equilibrium strategies (A.8)42def strategies(p: LQParams, D: float) -> dict:43 bH, bA = p.beta_H, p.beta_A44 xi = 1.0 / (1.0 - p.alpha * (D / p.M) * (1.0 - bA))45 c = p.kappa0 / (1.0 - p.alpha * ((1.0 - D) * bH + D * bA * xi))46 aH, aA = bH * c, bA * xi * c47 Abar = (1.0 - D) * aH + D * aA48 return dict(aH=aH, aA=aA, xi=xi, c=c, Abar=Abar, kappa_full=p.kappa0 / (1.0 - p.alpha))495051# --------------------------------------------------------------------------- second moments (A.9)52def moments(p: LQParams, D: float) -> dict:53 s = strategies(p, D)54 st2, sH2, sA2 = 1.0 / p.tau_theta, 1.0 / p.tau_H, 1.0 / p.tau_A55 aH, aA, Abar, kf = s['aH'], s['aA'], s['Abar'], s['kappa_full']56 V_noise = D ** 2 * aA ** 2 * sA2 / p.M57 V_track = (Abar - kf) ** 2 * st258 V_agg = V_track + V_noise59 # dispersion E (x_i - X)^260 V_disp_H = aH ** 2 * sH2 + (aH - Abar) ** 2 * st2 + V_noise61 V_disp_A = (aA - Abar) ** 2 * st2 + aA ** 2 * sA2 - 2 * aA * (D * aA / p.M) * sA2 + V_noise62 V_disp = (1.0 - D) * V_disp_H + D * V_disp_A63 # raw moments for the general quadratic welfare64 Ex2 = (1.0 - D) * aH ** 2 * (st2 + sH2) + D * aA ** 2 * (st2 + sA2)65 EX2 = Abar ** 2 * st2 + V_noise66 ExX = (1.0 - D) * aH * Abar * st2 + D * aA * (Abar * st2 + D * aA / p.M * sA2)67 Ext = ((1.0 - D) * aH + D * aA) * st268 EXt = Abar * st269 return dict(V_noise=V_noise, V_track=V_track, V_agg=V_agg, V_disp=V_disp,70 Ex2=Ex2, EX2=EX2, ExX=ExX, Ext=Ext, EXt=EXt, **s)717273# --------------------------------------------------------------------------- welfare (A.11) and values74def welfare_gross(p: LQParams, D: float) -> float:75 """Per-capita expected payoff E U(x_i, X, theta) (constant terms dropped), general quadratic form."""76 m = moments(p, D)77 a = p.alpha78 return (-m['Ex2'] / 2.0 + a * m['ExX'] - p.upsilon * a * m['EX2'] / 2.079 + p.kappa0 * m['Ext'] + p.U_Xtheta * m['EXt'])808182def welfare_decomposed(p: LQParams, D: float) -> dict:83 """Angeletos–Pavan decomposition under the CRS restriction: W = W* - 1/2[(1-alpha) V_agg + V_disp]."""84 m = moments(p, D)85 return dict(loss_agg=0.5 * (1.0 - p.alpha) * m['V_agg'], loss_disp=0.5 * m['V_disp'],86 W_rel=-0.5 * ((1.0 - p.alpha) * m['V_agg'] + m['V_disp']))878889def private_value(p: LQParams, D: float) -> float:90 """Gain in expected payoff of a measure-zero agent switching from own signal to model signal (units |U_xx|)."""91 m = moments(p, D)92 st2, sH2, sA2 = 1.0 / p.tau_theta, 1.0 / p.tau_H, 1.0 / p.tau_A93 g = p.kappa0 + p.alpha * m['Abar']94 Vx = g ** 2 * st2 + p.alpha ** 2 * (D * m['aA'] / p.M) ** 2 * p.M * sA295 loss_H = Vx - (g * st2) ** 2 / (st2 + sH2)96 loss_A = Vx - (g * st2 + p.alpha * (D * m['aA'] / p.M) * sA2) ** 2 / (st2 + sA2)97 return 0.5 * (loss_H - loss_A)9899100def social_marginal_value(p: LQParams, D: float, h: float = 1e-6) -> float:101 """d E U / dD (gross of fees), central finite difference of the exact quadratic welfare."""102 D0, D1 = max(0.0, D - h), min(1.0, D + h)103 return (welfare_gross(p, D1) - welfare_gross(p, D0)) / (D1 - D0)104105106def wedge(p: LQParams, D: float) -> float:107 """Private minus social marginal value; identically 0 under Lemma 2 (Proposition 3)."""108 return private_value(p, D) - social_marginal_value(p, D)109110111# --------------------------------------------------------------------------- delegation equilibrium (A.10)112def delegation_equilibrium(p: LQParams, cost_cdf, scale: float = 1.0, tol: float = 1e-12) -> float:113 """Fixed point D = F(scale * PV(D) - 0) where cost_cdf(c) = Pr(kappa_i + fee <= c).114115 PV is evaluated in units |U_xx|; `scale` converts to the units of the cost distribution.116 Solved by bisection on D in [0,1]; the map D -> F(scale*PV(D)) is continuous, so a fixed point exists (Brouwer).117 Uniqueness is checked numerically by scanning for sign changes of D - F(.).118 """119 f = lambda D: D - cost_cdf(scale * private_value(p, D))120 grid = np.linspace(0.0, 1.0, 201)121 vals = np.array([f(D) for D in grid])122 if vals[0] >= 0:123 return 0.0124 if vals[-1] <= 0:125 return 1.0126 idx = np.where(np.diff(np.sign(vals)) != 0)[0]127 lo, hi = grid[idx[0]], grid[idx[0] + 1]128 for _ in range(200):129 mid = 0.5 * (lo + hi)130 if f(mid) > 0:131 hi = mid132 else:133 lo = mid134 if hi - lo < tol:135 break136 return 0.5 * (lo + hi)137138139def planner_delegation(p: LQParams, cost_cdf, cost_quantile, scale: float = 1.0) -> float:140 """Planner's D: maximises E U(D) - ∫_0^D cost_quantile(q) dq / scale over a grid + golden refinement."""141 grid = np.linspace(0.0, 1.0, 401)142143 def obj(D):144 if D <= 0:145 return welfare_gross(p, 0.0)146 qs = np.linspace(0.0, D, 400)147 fees = np.trapezoid(cost_quantile(qs), qs) / scale148 return welfare_gross(p, D) - fees149150 vals = np.array([obj(D) for D in grid])151 return float(grid[int(np.argmax(vals))])152153154# --------------------------------------------------------------------------- Monte Carlo verifier155def simulate(p: LQParams, D: float, n_agents: int = 20000, n_draws: int = 4000, seed: int = 0) -> dict:156 """Simulate the signal game with the closed-form strategies; returns sample moments to compare with moments()."""157 rng = np.random.default_rng(seed)158 s = strategies(p, D)159 nD = int(round(D * n_agents))160 nH = n_agents - nD161 model_of = np.arange(nD) % p.M if nD > 0 else np.zeros(0, int)162 theta = rng.normal(0.0, p.tau_theta ** -0.5, n_draws)163 X = np.empty(n_draws)164 Ex2 = 0.0165 disp = 0.0166 for t in range(n_draws):167 e = rng.normal(0.0, p.tau_H ** -0.5, nH)168 eta = rng.normal(0.0, p.tau_A ** -0.5, p.M)169 xH = s['aH'] * (theta[t] + e)170 xA = s['aA'] * (theta[t] + eta[model_of]) if nD > 0 else np.zeros(0)171 x = np.concatenate([xH, xA])172 X[t] = x.mean()173 Ex2 += (x ** 2).mean()174 disp += ((x - X[t]) ** 2).mean()175 return dict(EX2=float(np.mean(X ** 2)), Ex2=Ex2 / n_draws, V_disp=disp / n_draws,176 EXt=float(np.mean(X * theta)), Abar_hat=float(np.sum(X * theta) / np.sum(theta ** 2)))177