1"""Validation diagnostics (Phase 3): Monte Carlo simulation of the listing process against the analytic solution,2residual checks, bounds, reproducibility across seeds and sample sizes. Writes results/logs/validation.json and3results/processed/val_*.csv.45Run: python -m src.diagnostics6"""7from __future__ import annotations89import json10import sys11import time12from pathlib import Path1314import numpy as np15import pandas as pd1617from .model import load_params, Params18from .solver import solve_equilibrium, solve_sellers1920ROOT = Path(__file__).resolve().parents[1]21PRO = ROOT / "results" / "processed"; LOG = ROOT / "results" / "logs"222324def simulate_listings(eq, p: Params, n_listings: int, seed: int) -> pd.DataFrame:25 """Simulate n_listings sellers: type c ~ f, meetings ~ Poisson(Lam_s), withdrawal ~ Exp(delta), match w ~ G,26 acceptance w >= v*(c), price = (1-beta)V - beta U + beta w. Returns one row per listing."""27 rng = np.random.default_rng(seed)28 types = rng.choice(len(p.c), size=n_listings, p=p.f)29 Lam = eq.Lam_s30 out = np.empty((n_listings, 5))31 for i, ti in enumerate(types):32 vstar, V = eq.sellers.vstar[ti], eq.sellers.V[ti]33 t = 0.0; sold = 0; price = np.nan; n_meet = 034 t_withdraw = rng.exponential(1.0 / p.delta) if p.delta > 0 else np.inf35 while True:36 t += rng.exponential(1.0 / Lam)37 if t > t_withdraw:38 t = t_withdraw; break39 n_meet += 140 w = float(p.G.sample(rng, 1)[0])41 if w >= vstar:42 sold = 1; price = (1 - p.beta) * V - p.beta * eq.U + p.beta * w; break43 out[i] = (ti, t, sold, price, n_meet)44 return pd.DataFrame(out, columns=["type", "duration", "sold", "price", "n_meetings"]).astype({"type": int, "sold": int})454647def main():48 p = load_params(ROOT / "config" / "baseline.yaml")49 eq = solve_equilibrium(p)50 report = {"config_hash": p.config_hash(), "timestamp": time.strftime("%Y-%m-%dT%H:%M:%S")}51 sb = eq.sellers52 report["residuals"] = {"seller_max_abs": float(np.abs(sb.resid).max()), "free_entry": float(eq.free_entry_resid),53 "stock_eq_max_abs": float(np.abs(p.s * p.f - (sb.eta + p.delta) * eq.n).max()),54 "omega_sum": float(eq.omega.sum()), "N_sum": float(eq.N.sum())}55 report["bounds"] = {"acc_in_[0,1]": bool(np.all((sb.acc >= 0) & (sb.acc <= 1))), "eta_nonneg": bool(np.all(sb.eta >= 0)),56 "prices_in_support": bool(np.all((sb.pbar >= p.G.lo) & (sb.pbar <= p.G.hi))),57 "V_monotone_dec": bool(np.all(np.diff(sb.V) < 0)), "eta_monotone_inc": bool(np.all(np.diff(sb.eta) > 0)),58 "participation": bool(np.all(sb.V >= p.w0))}59 # Monte Carlo cross-check60 rows = []61 for n in (2000, 20000, 100000):62 for seed in (1, 2, 3):63 t0 = time.time(); df = simulate_listings(eq, p, n, seed)64 sold = df[df.sold == 1]65 an_sold_share = float((p.f * sb.eta / (sb.eta + p.delta)).sum())66 an_dur = float((p.f / (sb.eta + p.delta)).sum())67 an_price = float((eq.omega * sb.pbar).sum())68 an_meet = float((p.f * eq.Lam_s / (sb.eta + p.delta)).sum())69 rows.append({"n": n, "seed": seed, "sold_share_mc": df.sold.mean(), "sold_share_an": an_sold_share,70 "mean_duration_mc": df.duration.mean(), "mean_duration_an": an_dur,71 "mean_price_mc": sold.price.mean(), "mean_price_an": an_price,72 "sd_price_mc": sold.price.std(), "sd_price_an": eq.moments(p)["price_sd_total"],73 "mean_meetings_mc": df.n_meetings.mean(), "mean_meetings_an": an_meet,74 "mean_tom_sold_mc": sold.duration.mean(), "mean_tom_sold_an": eq.moments(p)["mean_tom_sold"],75 "seconds": time.time() - t0})76 mc = pd.DataFrame(rows); mc.to_csv(PRO / "val_montecarlo.csv", index=False)77 big = mc[mc.n == 100000]78 report["montecarlo"] = {k: {"max_abs_rel_err": float((abs(big[f"{k}_mc"] - big[f"{k}_an"]) / big[f"{k}_an"]).max())}79 for k in ("sold_share", "mean_duration", "mean_price", "sd_price", "mean_meetings", "mean_tom_sold")}80 # per-type check with the largest sample81 df = simulate_listings(eq, p, 100000, 7)82 g = df.groupby("type").agg(sold_share=("sold", "mean"), mean_dur=("duration", "mean"), n=("sold", "size"),83 mean_price=("price", "mean"))84 g["sold_share_an"] = sb.eta / (sb.eta + p.delta); g["mean_dur_an"] = 1 / (sb.eta + p.delta); g["mean_price_an"] = sb.pbar85 g["c"] = p.c; g.to_csv(PRO / "val_montecarlo_types.csv")86 report["montecarlo_types_max_rel_err_price"] = float((abs(g.mean_price - g.mean_price_an) / g.mean_price_an).max())87 # reproducibility: same seed -> identical results88 d1 = simulate_listings(eq, p, 5000, 11); d2 = simulate_listings(eq, p, 5000, 11)89 report["reproducible_same_seed"] = bool(d1.equals(d2))90 # duration-dependence check (Lemma 3): exposure-based piecewise-constant hazard vs analytic mixture hazard91 bins = np.linspace(0, 0.6, 13); emp = []; ana = []92 mu = sb.eta + p.delta93 for a, b in zip(bins[:-1], bins[1:]):94 exposure = (np.clip(df.duration, a, b) - a).sum(); exits = ((df.duration >= a) & (df.duration < b)).sum()95 emp.append(exits / exposure if exposure > 0 else np.nan)96 mid = 0.5 * (a + b); wts = p.f * np.exp(-mu * mid); ana.append(float((wts * mu).sum() / wts.sum()))97 report["hazard_by_duration"] = {"bin_upper_years": bins[1:].tolist(), "empirical": [float(h) for h in emp], "analytic_mixture": ana,98 "max_abs_rel_err": float(np.nanmax(np.abs(np.array(emp) - np.array(ana)) / np.array(ana)))}99 report["hazard_decreasing_analytic"] = bool(np.all(np.diff(ana) < 0))100 report["hazard_decreasing_empirical_first3_vs_last3"] = bool(np.nanmean(emp[:3]) > np.nanmean(emp[-3:]))101 (LOG / "validation.json").write_text(json.dumps(report, indent=2))102 print(json.dumps(report, indent=2))103104105if __name__ == "__main__":106 main()107