1"""Experiment program for WP19 (Phase 4). Each experiment has a pre-registered hypothesis (stored with the outputs),2an exact configuration, CSV outputs in results/processed/, and a JSON record in results/raw/.34Run: python -m src.simulation --exp all (or --exp 1,3,9)5"""6from __future__ import annotations78import argparse9import json10import logging11import platform12import sys13import time14from pathlib import Path1516import numpy as np17import pandas as pd18import scipy1920from .model import load_params, Exponential, Pareto, Weibull, TruncNormal, Params21from .solver import (solve_equilibrium, partial_equilibrium, free_entry_residual, solve_equilibrium_posted,22 solve_equilibrium_financing, transition, frictionless_price, welfare_flow, planner_optimum)2324ROOT = Path(__file__).resolve().parents[1]25RAW = ROOT / "results" / "raw"26PRO = ROOT / "results" / "processed"27LOG = ROOT / "results" / "logs"28for d in (RAW, PRO, LOG):29 d.mkdir(parents=True, exist_ok=True)3031logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s",32 handlers=[logging.FileHandler(LOG / "simulation.log"), logging.StreamHandler(sys.stdout)])33log = logging.getLogger("wp19")3435HYPOTHESES = {36 1: "Baseline: unique stationary equilibrium; TOM of a few months; urgency discount of a few percent; price dispersion among identical houses with a positive between-type share; population hazard decreasing in duration (Lemma 3).",37 2: "Urgency gradient: V, pbar decrease and eta increases in c (Prop. 1); |pbar'| < |V'|; the gradient steepens (gap widens) when the arrival rate falls (Prop. 2(iii)).",38 3: "Arrival-rate collapse: lower A lowers prices and tightness, raises TOM and inventory, widens the urgency discount; GE response exceeds PE response (amplification); Gamma changes sign between thin and thick markets (Prop. 2(iv)).",39 4: "Financing: in the baseline screen, lower phi acts like lower A. In the contingency extension, mortgage deals need better matches; the cash-mortgage gap = delay-bargaining + outside-option + selection; the bargaining term rises with urgency (Prop. 7); tighter approval (lower phi_m, lower nu) lengthens TOM and lowers prices.",40 5: "Bargaining power: higher seller share beta raises prices and lowers tightness and hazards; total welfare is maximised near the Hosios value beta = 1-alpha with homogeneous sellers; with heterogeneous sellers the welfare-maximising beta deviates (numerical).",41 6: "Matching efficiency vs qualification: observationally equivalent in the baseline; in the extension a matching-efficiency improvement and an approval improvement calibrated to the same sales-hazard change have different price and TOM effects because approval shocks also change the pending wedge.",42 7: "Endogenous entry: the free-entry residual is monotone in theta (unique equilibrium); GE amplifies A and phi shocks relative to PE; lower entry/search costs raise theta and sellers' prices.",43 8: "Heterogeneous sellers: between-type dispersion rises with the spread of F; the transaction-price distribution is a mixture whose left tail is populated by urgent sellers; dispersion increases when liquidity falls.",44 9: "Dynamics: after a permanent negative A shock the observed index falls on impact, inventory builds up, volume drops then recovers; the composition term is transitory when delta=0; its sign depends on G (dampening for the truncated normal, amplifying for the exponential).",45 10: "Selection bias: impact decomposition of the index change into within-type and composition terms across G families and urgency spreads; the composition term is small for light tails, larger and amplifying for exponential/Pareto; steady-state composition effects vanish when delta=0.",46 11: "Posted prices: sellers hold out for better matches (lower hazard, higher price conditional on sale), all dispersion is between-type; the urgency gradient of prices is preserved.",47 12: "Sensitivity and convergence: qualitative results (signs) are robust to +-25% parameter changes and to the G family; equilibrium objects converge in the number of seller types and in tolerances; Monte Carlo simulation of the listing process matches analytic hazards and prices.",48}495051def env_record(p: Params) -> dict:52 return {"python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__,53 "pandas": pd.__version__, "config_hash": p.config_hash(), "seed": p.seed,54 "timestamp": time.strftime("%Y-%m-%dT%H:%M:%S")}555657def save(exp: int, name: str, df: pd.DataFrame, p: Params, meta: dict | None = None):58 path = PRO / f"exp{exp:02d}_{name}.csv"59 df.to_csv(path, index=False)60 rec = {"experiment": exp, "name": name, "hypothesis": HYPOTHESES[exp], "rows": len(df), "columns": list(df.columns),61 "env": env_record(p), **(meta or {})}62 (RAW / f"exp{exp:02d}_{name}.json").write_text(json.dumps(rec, indent=2, default=str))63 log.info("saved %s (%d rows)", path.name, len(df))646566def type_table(eq, p: Params) -> pd.DataFrame:67 sb = eq.sellers68 return pd.DataFrame({69 "c": p.c, "f": p.f, "V": sb.V, "vstar": sb.vstar, "acc": sb.acc, "eta": sb.eta, "pbar": sb.pbar,70 "D": 1 - sb.pbar / sb.pbar[0], "E_tau": 1 / (sb.eta + p.delta), "E_tau_months": 12 / (sb.eta + p.delta),71 "sold_prob": sb.eta / (sb.eta + p.delta), "L_holding_cost": p.c / (p.r + p.delta + sb.eta),72 "stock_share": eq.N, "sales_share": eq.omega, "M": sb.M,73 "price_sd_within": p.beta * np.sqrt([p.G.cond_var(v) for v in sb.vstar]),74 })757677# ------------------------------------------------------------------------------------------------78def exp01(p: Params):79 eq = solve_equilibrium(p)80 m = eq.moments(p)81 save(1, "baseline_moments", pd.DataFrame([m]), p, {"theta_roots": eq.theta_roots})82 save(1, "baseline_types", type_table(eq, p), p)83 tau = np.linspace(0, 2.0, 201)84 surv = np.exp(-np.outer(tau, eq.sellers.eta + p.delta))85 pop = surv @ p.f86 haz = -np.gradient(np.log(pop), tau)87 df = pd.DataFrame({"tau": tau, "tau_months": 12 * tau, "survival_population": pop, "hazard_population": haz})88 for i in (0, len(p.c) // 2, len(p.c) - 1):89 df[f"survival_c{i}"] = surv[:, i]90 save(1, "baseline_survival", df, p)91 # free-entry residual curve (uniqueness diagnostic)92 th = np.geomspace(p.theta_bounds[0], p.theta_bounds[1], 80)93 save(1, "free_entry_curve", pd.DataFrame({"theta": th, "residual": [free_entry_residual(t, p) for t in th]}), p)94 return eq959697def exp02(p: Params):98 eq = solve_equilibrium(p)99 cgrid = np.linspace(0.0, 0.5, 101)100 rows = []101 for mult in (0.5, 0.75, 1.0, 1.5, 2.0):102 pm = p.copy(c=cgrid, f=np.full_like(cgrid, 1 / len(cgrid)))103 pe = partial_equilibrium(eq.theta, pm.copy(A=p.A * mult)) # partial equilibrium: theta fixed at baseline104 sb = pe.sellers105 for i, c in enumerate(cgrid):106 rows.append({"A_mult": mult, "Lam_s": pe.Lam_s, "c": c, "V": sb.V[i], "vstar": sb.vstar[i], "eta": sb.eta[i],107 "pbar": sb.pbar[i], "D": 1 - sb.pbar[i] / sb.pbar[0], "E_tau_months": 12 / (sb.eta[i] + p.delta),108 "dV_dc_analytic": -1 / ((p.r + p.delta) + pe.Lam_s * p.beta * sb.acc[i]),109 "participation": sb.V[i] >= p.w0})110 df = pd.DataFrame(rows)111 df["dV_dc_fd"] = df.groupby("A_mult")["V"].transform(lambda s: np.gradient(s.values, cgrid))112 df["dpbar_dc_fd"] = df.groupby("A_mult")["pbar"].transform(lambda s: np.gradient(s.values, cgrid))113 save(2, "urgency_gradient", df, p, {"theta_fixed": eq.theta})114 # GE version: urgency discount of each baseline type as a function of A115 rows = []116 for mult in np.linspace(0.4, 1.6, 13):117 e = solve_equilibrium(p.copy(A=p.A * mult))118 for i, c in enumerate(p.c):119 rows.append({"A_mult": mult, "theta": e.theta, "c": c, "pbar": e.sellers.pbar[i], "V": e.sellers.V[i],120 "D": 1 - e.sellers.pbar[i] / e.sellers.pbar[0], "gapV": e.sellers.V[0] - e.sellers.V[i],121 "eta": e.sellers.eta[i]})122 save(2, "urgency_discount_vs_A_GE", pd.DataFrame(rows), p)123124125def exp03(p: Params):126 base = solve_equilibrium(p)127 rows = []128 for mult in np.concatenate([np.linspace(0.2, 1.0, 17), np.linspace(1.1, 3.0, 20)]):129 pm = p.copy(A=p.A * mult)130 ge = solve_equilibrium(pm); pe = partial_equilibrium(base.theta, pm)131 for kind, e in (("GE", ge), ("PE", pe)):132 m = e.moments(p)133 m.update({"A_mult": mult, "kind": kind, "Q_impact": float((e.sellers.eta * base.n).sum())})134 rows.append(m)135 save(3, "arrival_shock", pd.DataFrame(rows), p, {"theta_base": base.theta})136 # Gamma sign change for the patient type in partial equilibrium over a wide Lam_s range137 rows = []138 for L in np.geomspace(0.05, 2000, 60):139 pm = p.copy(c=np.array([0.0]), f=np.array([1.0]), delta=0.0)140 from .solver import solve_sellers141 sb = solve_sellers(L, p.k, pm)142 rows.append({"Lam_s": L, "pbar_patient": sb.pbar[0], "p_star": frictionless_price(p),143 "Gamma": 1 - sb.pbar[0] / frictionless_price(p), "vstar": sb.vstar[0]})144 save(3, "gamma_sign_change", pd.DataFrame(rows), p)145146147def exp04(p: Params):148 # (a) baseline screen149 rows = []150 for phi in np.linspace(0.4, 1.0, 13):151 e = solve_equilibrium(p.copy(phi=phi)); m = e.moments(p); m["phi"] = phi; rows.append(m)152 save(4, "qualification_screen", pd.DataFrame(rows), p)153 # (b) contingency extension154 def fin_rows(fe, label, val):155 tw_c = fe.Q_c / fe.Q_c.sum(); tw_m = fe.Q_m / fe.Q_m.sum()156 out = []157 for i, c in enumerate(p.c):158 out.append({label: val, "c": c, "f": p.f[i], "theta": fe.theta, "U_c": fe.U_c, "U_m": fe.U_m, "pi_tilde": fe.pi_tilde,159 "V": fe.V[i], "vstar_c": fe.vstar_c[i], "vstar_m": fe.vstar_m[i], "eta_c": fe.eta_c[i], "eta_m": fe.eta_m[i],160 "pbar_c": fe.pbar_c[i], "pbar_m": fe.pbar_m[i], "gap": fe.pbar_m[i] - fe.pbar_c[i],161 "Delta_b": fe.Delta_b[i], "outside_term": fe.outside_term, "sel_term": fe.sel_term[i],162 "share_cash_sales": fe.Q_c[i] / (fe.Q_c[i] + fe.Q_m[i]),163 "E_tau_listed_months": 12 / (fe.eta_c[i] + p.phi_m * fe.eta_m[i] + p.delta),164 "S": fe.S, "Q": fe.Q, "pending_stock": fe.m_pending.sum(), "converged": fe.converged,165 "E_total_duration_months": 12 * (fe.S + fe.m_pending.sum()) / p.s,166 "buyer_participation_ok": bool(min(fe.U_c, fe.U_m) >= 0), "seller_participation_ok": bool(np.all(fe.V >= p.w0)),167 "price_mean_cash": float((tw_c * fe.pbar_c).sum()), "price_mean_mortgage": float((tw_m * fe.pbar_m).sum()),168 "cash_share_of_sales": fe.Q_c.sum() / fe.Q})169 return out170 rows = []171 fe0 = solve_equilibrium_financing(p); rows += fin_rows(fe0, "phi_m", p.phi_m)172 for phi in (0.5, 0.6, 0.7, 0.75, 0.8, 0.9, 0.95, 1.0):173 rows += fin_rows(solve_equilibrium_financing(p.copy(phi_m=phi)), "phi_m", phi)174 save(4, "financing_phi", pd.DataFrame(rows), p)175 rows = []176 for nu in (3.0, 6.0, 12.0, 24.0, 52.0, 1e4):177 rows += fin_rows(solve_equilibrium_financing(p.copy(nu=nu)), "nu", nu)178 save(4, "financing_nu", pd.DataFrame(rows), p)179 rows = []180 for pi in (0.1, 0.3, 0.5, 0.7, 0.9):181 rows += fin_rows(solve_equilibrium_financing(p.copy(pi=pi)), "pi", pi)182 save(4, "financing_pi", pd.DataFrame(rows), p)183 rows = []184 for beta in (0.2, 0.35, 0.5, 0.65, 0.8):185 rows += fin_rows(solve_equilibrium_financing(p.copy(beta=beta)), "beta", beta)186 save(4, "financing_beta", pd.DataFrame(rows), p)187188189def exp05(p: Params):190 rows = []191 hom = p.copy(c=np.array([float((p.c * p.f).sum())]), f=np.array([1.0]))192 for beta in np.linspace(0.1, 0.9, 17):193 for label, pp in (("heterogeneous", p), ("homogeneous", hom)):194 e = solve_equilibrium(pp.copy(beta=beta)); m = e.moments(pp)195 m.update({"beta": beta, "sellers": label, "hosios_beta": 1 - p.alpha})196 rows.append(m)197 df = pd.DataFrame(rows)198 save(5, "bargaining_power", df, p)199 # planner benchmark for both economies200 rows = []201 for label, pp in (("heterogeneous", p), ("homogeneous", hom)):202 e = solve_equilibrium(pp); pl = planner_optimum(pp)203 rows.append({"sellers": label, "theta_decentralised": e.theta, "theta_planner": pl["theta"],204 "welfare_decentralised": welfare_flow(e, pp)["total"], "welfare_planner": pl["welfare"],205 "vstar_decentralised": json.dumps(np.round(e.sellers.vstar, 5).tolist()),206 "vstar_planner": json.dumps(np.round(pl["vstar"], 5).tolist()), "planner_success": pl["success"]})207 for beta in np.linspace(0.1, 0.9, 33):208 pass209 save(5, "planner_benchmark", pd.DataFrame(rows), p)210 # welfare surface over (beta, alpha)211 rows = []212 for alpha in (0.3, 0.4, 0.5, 0.6, 0.7):213 for beta in np.linspace(0.1, 0.9, 17):214 for label, pp in (("heterogeneous", p), ("homogeneous", hom)):215 e = solve_equilibrium(pp.copy(alpha=alpha, beta=beta))216 rows.append({"alpha": alpha, "beta": beta, "sellers": label, "welfare_total": welfare_flow(e, pp)["total"],217 "theta": e.theta, "converged": e.converged})218 save(5, "welfare_surface", pd.DataFrame(rows), p)219220221def exp06(p: Params):222 """Matching efficiency vs approval, matched on the change in the aggregate sales hazard (extension)."""223 fe0 = solve_equilibrium_financing(p)224 def agg(fe):225 n = fe.n; eta = fe.eta_c + p.phi_m * fe.eta_m226 return float((eta * n).sum() / n.sum())227 h0 = agg(fe0)228 rows = []229 for A_mult in (0.7, 0.85, 1.0, 1.15, 1.3):230 fe = solve_equilibrium_financing(p.copy(A=p.A * A_mult))231 rows.append({"shock": "A", "value": A_mult, "agg_hazard": agg(fe), "theta": fe.theta,232 "price_mean": float(((fe.Q_c * fe.pbar_c + fe.Q_m * fe.pbar_m).sum()) / fe.Q), "S": fe.S, "Q": fe.Q,233 "gap_mean": float(((fe.Q_m * (fe.pbar_m - fe.pbar_c)).sum()) / fe.Q_m.sum()),234 "mean_tom_months": 12 * float((p.f / (fe.eta_c + p.phi_m * fe.eta_m + p.delta)).sum()),235 "total_duration_months": 12 * (fe.S + fe.m_pending.sum()) / p.s, "U_m": fe.U_m, "U_c": fe.U_c,236 "buyer_participation_ok": bool(min(fe.U_c, fe.U_m) >= 0),237 "cash_share": fe.Q_c.sum() / fe.Q, "converged": fe.converged})238 for phi in (0.55, 0.7, 0.75, 0.8, 0.85, 0.95, 1.0):239 fe = solve_equilibrium_financing(p.copy(phi_m=phi))240 rows.append({"shock": "phi_m", "value": phi, "agg_hazard": agg(fe), "theta": fe.theta,241 "price_mean": float(((fe.Q_c * fe.pbar_c + fe.Q_m * fe.pbar_m).sum()) / fe.Q), "S": fe.S, "Q": fe.Q,242 "gap_mean": float(((fe.Q_m * (fe.pbar_m - fe.pbar_c)).sum()) / fe.Q_m.sum()),243 "mean_tom_months": 12 * float((p.f / (fe.eta_c + p.phi_m * fe.eta_m + p.delta)).sum()),244 "total_duration_months": 12 * (fe.S + fe.m_pending.sum()) / p.s, "U_m": fe.U_m, "U_c": fe.U_c,245 "buyer_participation_ok": bool(min(fe.U_c, fe.U_m) >= 0),246 "cash_share": fe.Q_c.sum() / fe.Q, "converged": fe.converged})247 df = pd.DataFrame(rows); df["agg_hazard_base"] = h0248 save(6, "matching_vs_approval", df, p)249250251def exp07(p: Params):252 base = solve_equilibrium(p)253 rows = []254 for k_mult in (0.5, 1.0, 2.0):255 for kb_mult in (0.5, 1.0, 2.0):256 e = solve_equilibrium(p.copy(k=p.k * k_mult, kappa_b=p.kappa_b * kb_mult)); m = e.moments(p)257 m.update({"k_mult": k_mult, "kappa_b_mult": kb_mult}); rows.append(m)258 save(7, "entry_costs", pd.DataFrame(rows), p)259 # amplification: elasticities of Lam_s, price, TOM to A and phi in GE vs PE260 rows = []261 for var in ("A", "phi"):262 for h in (0.02,):263 up = solve_equilibrium(p.copy(**{var: getattr(p, var) * (1 + h)})); dn = solve_equilibrium(p.copy(**{var: getattr(p, var) * (1 - h)}))264 upe = partial_equilibrium(base.theta, p.copy(**{var: getattr(p, var) * (1 + h)})); dne = partial_equilibrium(base.theta, p.copy(**{var: getattr(p, var) * (1 - h)}))265 for kind, (u, d) in (("GE", (up, dn)), ("PE", (upe, dne))):266 mu, md = u.moments(p), d.moments(p)267 el = lambda key: (np.log(mu[key]) - np.log(md[key])) / (np.log(1 + h) - np.log(1 - h))268 rows.append({"shock": var, "kind": kind, "el_Lam_s": el("Lam_s"), "el_theta": el("theta") if kind == "GE" else 0.0,269 "el_price": el("price_mean"), "el_tom": el("mean_tom_sold"), "el_Q": el("Q"), "el_S": el("S"),270 "el_D_urgent": el("D_urgent")})271 save(7, "amplification", pd.DataFrame(rows), p)272 # uniqueness check across parameter draws: count sign changes of R(theta)273 rng = np.random.default_rng(p.seed)274 rows = []275 for i in range(60):276 pp = p.copy(A=p.A * rng.uniform(0.5, 1.5), beta=rng.uniform(0.2, 0.8), k=p.k * rng.uniform(0.5, 2), kappa_b=p.kappa_b * rng.uniform(0.5, 2),277 delta=rng.uniform(0.1, 1.0), phi=rng.uniform(0.5, 1.0))278 e = solve_equilibrium(pp)279 rows.append({"draw": i, "A": pp.A, "beta": pp.beta, "k": pp.k, "kappa_b": pp.kappa_b, "delta": pp.delta, "phi": pp.phi,280 "n_roots": len(e.theta_roots), "theta": e.theta, "converged": e.converged})281 save(7, "uniqueness_draws", pd.DataFrame(rows), p)282283284def price_density(eq, p: Params, grid: np.ndarray) -> pd.DataFrame:285 """Mixture density of transaction prices: p = (1-beta)V(c) - beta U + beta w, w >= v*(c), weights omega(c)."""286 dens = np.zeros_like(grid); parts = {}287 for i in range(len(p.c)):288 V, vs = eq.sellers.V[i], eq.sellers.vstar[i]289 w = (grid - (1 - p.beta) * V + p.beta * eq.U) / p.beta290 di = np.where(w >= vs, np.asarray(p.G.pdf(w)) / max(float(p.G.surv(vs)), 1e-300) / p.beta, 0.0)291 parts[f"type{i}"] = eq.omega[i] * di; dens += eq.omega[i] * di292 df = pd.DataFrame({"price": grid, "density": dens}); 293 for k_, v in parts.items(): df[k_] = v294 return df295296297def exp08(p: Params):298 rows = []299 for c_max in (0.05, 0.1, 0.2, 0.3, 0.45, 0.6):300 for n_types in (9,):301 pp = p.with_types(0.02, c_max, n_types, 2.0, 4.0)302 for mult in (0.7, 1.0, 1.3):303 e = solve_equilibrium(pp.copy(A=p.A * mult)); m = e.moments(pp)304 m.update({"c_max": c_max, "A_mult": mult, "c_sd": float(np.sqrt((pp.f * (pp.c - (pp.f * pp.c).sum()) ** 2).sum())),305 "pbar_range": e.sellers.pbar[0] - e.sellers.pbar[-1], "V_range": e.sellers.V[0] - e.sellers.V[-1]})306 rows.append(m)307 save(8, "heterogeneity_spread", pd.DataFrame(rows), p)308 e = solve_equilibrium(p)309 grid = np.linspace(0.8, 1.25, 451)310 save(8, "price_density_baseline", price_density(e, p, grid), p)311 el = solve_equilibrium(p.copy(A=p.A * 0.6))312 save(8, "price_density_lowA", price_density(el, p, grid), p, {"A_mult": 0.6})313 # beta = 0 limit and beta = 1 limit decomposition314 rows = []315 for beta in (1e-6, 0.25, 0.5, 0.75, 1 - 1e-6):316 e = solve_equilibrium(p.copy(beta=beta)); m = e.moments(p); m["beta"] = beta; rows.append(m)317 save(8, "dispersion_vs_beta", pd.DataFrame(rows), p)318319320def exp09(p: Params):321 T, dt = p.T, p.dt322 def run(label, p0, path, extra=None):323 tr = transition(p0, path, T=T, dt=dt)324 df = pd.DataFrame({"t": tr.t, "theta": tr.theta, "P_obs": tr.P_obs, "P_fix": tr.P_fix, "Q": tr.Q, "S": tr.S,325 "share_urgent_sales": tr.share_urgent_sales, "mean_tom_sold_months": 12 * tr.mean_tom_sold})326 for i in (0, len(p0.c) - 1):327 df[f"pbar_c{i}"] = tr.pbar[:, i]; df[f"eta_c{i}"] = tr.eta[:, i]; df[f"n_c{i}"] = tr.n[:, i]328 df["scenario"] = label329 save(9, f"transition_{label}", df, p0, {"converged": tr.converged, "iterations": tr.iterations, "max_update": tr.max_update,330 **(extra or {})})331 return tr332 p_lo = p.copy(A=p.A * 0.8); p_hi = p.copy(A=p.A * 1.2); p_phi = p.copy(phi=p.phi * 0.8)333 run("perm_A_minus20", p, lambda t: p_lo)334 run("perm_A_plus20", p, lambda t: p_hi)335 run("temp_A_minus20_1y", p, lambda t: p_lo if t < 1.0 else p)336 run("perm_phi_minus20", p, lambda t: p_phi)337 pe = p.copy(G=Exponential(rate=8.0, lo=0.7)); run("perm_A_minus20_exponential", pe, lambda t: pe.copy(A=p.A * 0.8))338 pd0 = p.copy(delta=0.0); run("perm_A_minus20_delta0", pd0, lambda t: pd0.copy(A=p.A * 0.8))339340341def exp10(p: Params):342 dists = {"truncnorm": p.G, "exponential": Exponential(rate=8.0, lo=0.7), "weibull": Weibull(shape=2.0, scale=0.25, lo=0.75),343 "pareto": Pareto(a=4.0, lo=0.8)}344 rows = []345 for gname, G in dists.items():346 for c_max in (0.1, 0.3, 0.6):347 pp = p.with_types(0.02, c_max, 9, 2.0, 4.0).copy(G=G)348 base = solve_equilibrium(pp)349 if not base.converged:350 rows.append({"G": gname, "c_max": c_max, "note": "no equilibrium"}); continue351 for mult in (0.7, 0.8, 0.9, 1.1, 1.2, 1.3):352 pm = pp.copy(A=pp.A * mult)353 # impact (stock predetermined, theta re-solved instantly by free entry? -> use PE at base theta and GE both)354 for kind in ("PE", "GE"):355 e = partial_equilibrium(base.theta, pm) if kind == "PE" else solve_equilibrium(pm)356 eta_new = e.sellers.eta357 om_imp = eta_new * base.n / (eta_new * base.n).sum()358 P0 = float((base.omega * base.sellers.pbar).sum())359 P_fix = float((base.omega * e.sellers.pbar).sum()); P_imp = float((om_imp * e.sellers.pbar).sum())360 P_ss = float((e.omega * e.sellers.pbar).sum())361 eps = (np.log(eta_new) - np.log(base.sellers.eta)) / np.log(mult)362 rows.append({"G": gname, "c_max": c_max, "A_mult": mult, "kind": kind, "theta_base": base.theta, "theta_new": e.theta,363 "P0": P0, "P_fix": P_fix, "P_impact": P_imp, "P_ss": P_ss,364 "within_log": np.log(P_fix / P0), "composition_impact_log": np.log(P_imp / P_fix),365 "composition_ss_log": np.log(P_ss / P_fix), "total_impact_log": np.log(P_imp / P0),366 "amplification_ratio_impact": np.log(P_imp / P0) / np.log(P_fix / P0),367 "eps_patient": eps[0], "eps_urgent": eps[-1], "eps_slope_sign": np.sign(eps[-1] - eps[0]),368 "urgent_share_base": base.omega[p.c >= np.median(p.c)].sum() if len(p.c) == 9 else np.nan,369 "urgent_share_impact": om_imp[-4:].sum(), "urgent_share_base_": base.omega[-4:].sum(),370 "is_ifr": G.is_ifr(), "converged": e.converged})371 save(10, "selection_decomposition", pd.DataFrame(rows), p)372 # hazard-elasticity profiles and h_G*M curves373 rows = []374 for gname, G in dists.items():375 v = np.linspace(G.lo + 1e-3, (G.hi if np.isfinite(G.hi) else G.lo + 1.0) - 1e-3, 200)376 hm = np.asarray(G.hazard(v)) * np.asarray(G.M(v))377 for vi, hi_ in zip(v, hm):378 rows.append({"G": gname, "v": vi, "hG_times_M": hi_})379 save(10, "hazard_times_mrl", pd.DataFrame(rows), p)380381382def exp11(p: Params):383 rows = []384 for mult in (0.7, 1.0, 1.3):385 pm = p.copy(A=p.A * mult)386 en = solve_equilibrium(pm); ep = solve_equilibrium_posted(pm)387 for proto, e in (("nash", en), ("posted", ep)):388 m = e.moments(pm); m.update({"protocol": proto, "A_mult": mult}); rows.append(m)389 if mult == 1.0:390 save(11, f"types_{proto}", type_table(e, pm), pm)391 save(11, "protocol_comparison", pd.DataFrame(rows), p)392393394def exp12(p: Params):395 base = solve_equilibrium(p).moments(p)396 keys = ["theta", "price_mean", "mean_tom_sold_months", "D_urgent", "D_mean", "price_between_share", "price_sd_total", "Q", "S", "welfare_total"]397 rows = []398 for par in ("r", "delta", "w0", "k", "kappa_b", "phi", "A", "alpha", "beta"):399 for mult in (0.75, 1.25):400 val = getattr(p, par) * mult401 if par in ("phi", "alpha", "beta"):402 val = min(val, 0.98)403 e = solve_equilibrium(p.copy(**{par: val})); m = e.moments(p)404 row = {"parameter": par, "mult": mult, "value": val, "converged": e.converged}405 for k_ in keys:406 row[k_] = m[k_]; row[f"d_{k_}"] = m[k_] - base[k_]; row[f"dlog_{k_}"] = np.log(m[k_] / base[k_]) if m[k_] > 0 and base[k_] > 0 else np.nan407 rows.append(row)408 for sig in (0.05, 0.15, 0.2):409 G = TruncNormal(1.0, sig, 1 - 3 * sig, 1 + 3 * sig); e = solve_equilibrium(p.copy(G=G)); m = e.moments(p)410 row = {"parameter": "sigma_w", "mult": sig / 0.1, "value": sig, "converged": e.converged}411 for k_ in keys: row[k_] = m[k_]; row[f"d_{k_}"] = m[k_] - base[k_]; row[f"dlog_{k_}"] = np.log(m[k_] / base[k_]) if m[k_] > 0 and base[k_] > 0 else np.nan412 rows.append(row)413 for cm in (0.2, 0.45):414 pp = p.with_types(0.02, cm, 9, 2.0, 4.0); e = solve_equilibrium(pp); m = e.moments(pp)415 row = {"parameter": "c_max", "mult": cm / 0.3, "value": cm, "converged": e.converged}416 for k_ in keys: row[k_] = m[k_]; row[f"d_{k_}"] = m[k_] - base[k_]; row[f"dlog_{k_}"] = np.log(m[k_] / base[k_]) if m[k_] > 0 and base[k_] > 0 else np.nan417 rows.append(row)418 df = pd.DataFrame(rows); 419 for k_ in keys: df[f"base_{k_}"] = base[k_]420 save(12, "sensitivity_oat", df, p)421 # alternative G families at baseline422 rows = []423 for gname, G in (("truncnorm", p.G), ("exponential", Exponential(rate=8.0, lo=0.7)), ("weibull", Weibull(shape=2.0, scale=0.25, lo=0.75)),424 ("weibull_dfr", Weibull(shape=0.8, scale=0.25, lo=0.75)), ("pareto", Pareto(a=4.0, lo=0.8))):425 e = solve_equilibrium(p.copy(G=G)); m = e.moments(p); m.update({"G": gname, "G_mean": G.mean(), "is_ifr": G.is_ifr()}); rows.append(m)426 save(12, "alternative_distributions", pd.DataFrame(rows), p)427 # grid convergence in the number of seller types and in tolerance428 rows = []429 for n_types in (1, 3, 5, 9, 15, 25, 41, 81):430 pp = p.with_types(0.02, 0.3, n_types, 2.0, 4.0); t0 = time.time(); e = solve_equilibrium(pp); m = e.moments(pp)431 m.update({"n_types": n_types, "seconds": time.time() - t0}); rows.append(m)432 save(12, "convergence_types", pd.DataFrame(rows), p)433 rows = []434 for tol in (1e-6, 1e-8, 1e-10, 1e-12):435 for gp in (20, 60, 120):436 e = solve_equilibrium(p.copy(root_xtol=tol, grid_check_points=gp)); m = e.moments(p); m.update({"root_xtol": tol, "grid_points": gp}); rows.append(m)437 save(12, "convergence_tolerance", pd.DataFrame(rows), p)438 # dt convergence of the transition439 rows = []440 p_lo = p.copy(A=p.A * 0.8)441 for dt in (0.1, 0.05, 0.02, 0.01):442 t0 = time.time(); tr = transition(p, lambda t: p_lo, T=4.0, dt=dt)443 j = int(round(1.0 / dt))444 rows.append({"dt": dt, "theta_impact": tr.theta[0], "P_obs_impact": tr.P_obs[0], "P_obs_1y": tr.P_obs[j], "Q_1y": tr.Q[j],445 "S_1y": tr.S[j], "theta_end": tr.theta[-1], "converged": tr.converged, "iterations": tr.iterations, "seconds": time.time() - t0})446 save(12, "convergence_dt", pd.DataFrame(rows), p)447 # global sensitivity: Latin-hypercube-like random draws448 rng = np.random.default_rng(p.seed + 1)449 rows = []450 for i in range(150):451 pp = p.copy(r=rng.uniform(0.02, 0.08), delta=rng.uniform(0.1, 1.0), k=rng.uniform(0.005, 0.05), kappa_b=rng.uniform(0.02, 0.1),452 phi=rng.uniform(0.5, 1.0), A=rng.uniform(2.0, 8.0), alpha=rng.uniform(0.3, 0.7), beta=rng.uniform(0.2, 0.8))453 e = solve_equilibrium(pp); m = e.moments(pp)454 m.update({"draw": i, "r": pp.r, "delta": pp.delta, "k": pp.k, "kappa_b": pp.kappa_b, "phi": pp.phi, "A": pp.A, "alpha": pp.alpha, "beta": pp.beta,455 "sign_D_urgent_pos": m["D_urgent"] > 0, "sign_eta_increasing": bool(np.all(np.diff(e.sellers.eta) > 0)),456 "sign_pbar_decreasing": bool(np.all(np.diff(e.sellers.pbar) < 0))})457 rows.append(m)458 save(12, "global_sensitivity", pd.DataFrame(rows), p)459460461EXPERIMENTS = {1: exp01, 2: exp02, 3: exp03, 4: exp04, 5: exp05, 6: exp06, 7: exp07, 8: exp08, 9: exp09, 10: exp10, 11: exp11, 12: exp12}462463464def main(argv=None):465 ap = argparse.ArgumentParser()466 ap.add_argument("--config", default=str(ROOT / "config" / "baseline.yaml"))467 ap.add_argument("--exp", default="all")468 a = ap.parse_args(argv)469 p = load_params(a.config)470 np.random.seed(p.seed)471 todo = sorted(EXPERIMENTS) if a.exp == "all" else [int(x) for x in a.exp.split(",")]472 log.info("config %s hash %s; experiments %s", a.config, p.config_hash(), todo)473 for e in todo:474 t0 = time.time(); log.info("=== Experiment %d: %s", e, HYPOTHESES[e][:80])475 try:476 EXPERIMENTS[e](p)477 log.info("experiment %d done in %.1f s", e, time.time() - t0)478 except Exception as ex: # report, never hide479 log.exception("experiment %d FAILED: %s", e, ex)480 (RAW / f"exp{e:02d}_FAILED.json").write_text(json.dumps({"experiment": e, "error": repr(ex)}))481482483if __name__ == "__main__":484 main()485