TeX 53.2%
Python 46.6%
1"""WP18 — experiment orchestration (CLAUDE.md Phase 3, experiments 1-10).23Usage: python -m src.simulation [--groups g1 g2 ...] [--workers N] [--quick]45Outputs (all under results/data/, every number used in a table or figure comes from here):6 parameters.csv period-unit parameters with sources7 baseline/summary.json prices, moments, convergence, Bellman residuals, thresholds8 baseline/profiles_by_age.csv ownership / wealth / debt / consumption by age9 baseline/by_gift_quintile.csv ownership and net worth by age x inheritance quintile10 baseline/mobility.json parent-child tenure transition, rank-rank slopes, lambda11 baseline/policy_age*.csv policy functions on the cash grid (figures)12 baseline/thresholds.csv ownership-entry thresholds by age and z (+ HTM check)13 baseline/cross_section.npz pooled cross-section (compressed)14 baseline/equal_income_paths.csv experiment 2 (deterministic equal-income paths)15 baseline/price_drop_accounting.csv experiment 5 (PE accounting of unanticipated price drops)16 <group>/<name>.json, <name>_profiles.csv stationary scenarios (GE or PE, labelled)17 inflation/<name>.json, _path.csv, _impact.csv transitions (GE / PE / static, labelled)18 welfare/welfare_<name>.csv consumption-equivalent variations (experiment 9)19 multiplicity/price_scan_<name>.csv experiment 3/7 multiplicity search20 diagnostics.json / diagnostics.csv convergence and residuals of every run21Logs: results/logs/<name>.jsonl (structured) and results/logs/simulate.log (summary).22"""23from __future__ import annotations2425import argparse26import json27import os28import pickle29import time30from concurrent.futures import ProcessPoolExecutor31from multiprocessing import get_context32from pathlib import Path3334import numpy as np35import pandas as pd3637from . import calibration as C38from . import distributions as D39from . import equilibrium as E40from . import model as M41from .household_solver import NEG, ownership_threshold4243DATA = M.RESULTS_DIR / "data"44LOGS = M.RESULTS_DIR / "logs"454647# --------------------------------------------------------------------------------------48# helpers49# --------------------------------------------------------------------------------------50def _js(o):51 if isinstance(o, dict):52 return {str(k): _js(v) for k, v in o.items()}53 if isinstance(o, (list, tuple)):54 return [_js(v) for v in o]55 if isinstance(o, np.ndarray):56 return o.tolist()57 if isinstance(o, (np.floating,)):58 return float(o)59 if isinstance(o, (np.integer,)):60 return int(o)61 if isinstance(o, (np.bool_,)):62 return bool(o)63 if isinstance(o, float) and not np.isfinite(o):64 return None65 return o666768def save_json(path: Path, obj: dict):69 path.parent.mkdir(parents=True, exist_ok=True)70 with open(path, "w") as fh:71 json.dump(_js(obj), fh, indent=1)727374def apply_overrides(cfg: dict, overrides: dict | None) -> dict:75 cfg = json.loads(json.dumps(cfg))76 for key, val in (overrides or {}).items():77 node = cfg78 parts = key.split(".")79 for part in parts[:-1]:80 node = node[part]81 node[parts[-1]] = val82 return cfg838485# --------------------------------------------------------------------------------------86# summaries of a stationary result87# --------------------------------------------------------------------------------------88def profiles_by_age(P: M.Params, rec: D.PeriodRecord) -> pd.DataFrame:89 nw = rec.net_worth1(P)90 own = rec.owner1()91 rows = []92 for j in range(1, P.J + 1):93 m = rec.age == j94 o = m & own; r = m & ~own95 hv = np.where(own, P.h_grid[np.maximum(rec.h1, 0)], 0.0)96 rows.append(dict(age=j, age_years=f"{25 + 5 * (j - 1)}-{29 + 5 * (j - 1)}", ownership=float(own[m].mean()),97 buy_rate=float((rec.choice[m] == 2).mean()), sell_to_rent_rate=float(((rec.choice[m] == 3) & (rec.h0[m] >= 0)).mean()),98 refinance_rate=float((rec.choice[m] == 1).mean()),99 nw_mean=float(nw[m].mean()), nw_median=float(np.median(nw[m])),100 nw_p25=float(np.quantile(nw[m], 0.25)), nw_p75=float(np.quantile(nw[m], 0.75)),101 nw_owner_median=float(np.median(nw[o])) if o.any() else np.nan,102 nw_renter_median=float(np.median(nw[r])) if r.any() else np.nan,103 nw_owner_mean=float(nw[o].mean()) if o.any() else np.nan,104 nw_renter_mean=float(nw[r].mean()) if r.any() else np.nan,105 liquid_mean=float(rec.a1[m].mean()), debt_mean=float(rec.b1[m].mean()),106 debt_owner_mean=float(rec.b1[o].mean()) if o.any() else np.nan,107 ltv_owner_mean=float((rec.b1[o] / (rec.p * hv[o])).mean()) if o.any() else np.nan,108 share_mortgage=float((rec.b1[m] > 1e-8).mean()),109 house_size_mean=float(hv[o].mean()) if o.any() else np.nan,110 rent_size_mean=float(rec.s[r].mean()) if r.any() else np.nan,111 c_mean=float(rec.c[m].mean()), c_owner=float(rec.c[o].mean()) if o.any() else np.nan,112 c_renter=float(rec.c[r].mean()) if r.any() else np.nan,113 income_mean=float((P.ybar[j - 1] * P.z_grid[rec.z[m]]).mean()),114 gini_nw=D.gini(nw[m]), n=int(m.sum())))115 return pd.DataFrame(rows)116117118def by_gift_quintile(P: M.Params, rec: D.PeriodRecord) -> pd.DataFrame:119 """Ownership and net worth by age and quintile of the gift received at age 1 (young and old120 members both carry it). Quintiles are defined on the age-1 gift distribution; ties at zero121 are broken by rank so that quintiles have equal mass (zero-gift households fall in the122 bottom groups)."""123 g = rec.gift_received124 ok = np.isfinite(g)125 g1 = np.sort(g[(rec.age == 1) & ok])126 # rank in the age-1 gift distribution with uniform tie-breaking (mass at zero is split evenly)127 rng = np.random.default_rng(3)128 lo = np.searchsorted(g1, g[ok], side="left"); hi = np.searchsorted(g1, g[ok], side="right")129 rank = (lo + rng.uniform(size=ok.sum()) * (hi - lo)) / len(g1)130 q = np.minimum((rank * 5).astype(int), 4) # 0..4131 nw = rec.net_worth1(P)[ok]; own = rec.owner1()[ok]; age = rec.age[ok]132 rows = []133 for j in range(1, P.J + 1):134 for k in range(5):135 m = (age == j) & (q == k)136 if m.sum() == 0:137 continue138 rows.append(dict(age=j, gift_quintile=k + 1, ownership=float(own[m].mean()), nw_mean=float(nw[m].mean()),139 nw_median=float(np.median(nw[m])), gift_mean=float(g[ok][m].mean()), n=int(m.sum())))140 return pd.DataFrame(rows)141142143def mobility_stats(P: M.Params, rec: D.PeriodRecord) -> dict:144 """Parent-child tenure transition at age j_gift (both measured at age 7, the gift age),145 persistence lambda = q_OO - q_RO (Part A analogue), rank-rank slope of net worth at age 7."""146 m = (rec.age == P.j_gift) & rec.is_old & (rec.parent_owner >= 0)147 child_own = rec.owner1()[m].astype(float)148 par_own = rec.parent_owner[m]149 q_OO = float(child_own[par_own == 1].mean()) if (par_own == 1).any() else np.nan150 q_RO = float(child_own[par_own == 0].mean()) if (par_own == 0).any() else np.nan151 nw = rec.net_worth1(P)[m]152 rr = D.rank_rank_slope(rec.parent_nw[m], nw)153 corr = float(np.corrcoef(par_own, child_own)[0, 1]) if par_own.std() > 0 and child_own.std() > 0 else np.nan154 lam = q_OO - q_RO155 # upward mobility of renter dynasties: child owns at age 3 given parent renter at the gift156 m3 = (rec.age == 3) & (~rec.is_old) & (rec.parent_owner >= 0)157 own3 = rec.owner1()[m3]; po3 = rec.parent_owner[m3]158 return dict(q_OO=q_OO, q_RO=q_RO, lambda_persistence=lam, half_life_generations=float(np.log(2) / -np.log(lam)) if 0 < lam < 1 else np.inf,159 tenure_correlation=corr, rank_rank_slope_nw=rr,160 q_RO_age3=float(own3[po3 == 0].mean()) if (po3 == 0).any() else np.nan,161 q_OO_age3=float(own3[po3 == 1].mean()) if (po3 == 1).any() else np.nan,162 share_parent_owner=float((par_own == 1).mean()), n=int(m.sum()))163164165def policy_tables(P: M.Params, res: E.SSResult) -> dict[str, pd.DataFrame]:166 out = {}167 pr = res.ev.prices168 for j in (1, 3, 7):169 pol = res.ev.pols[j - 1]170 rows = []171 for iz in range(P.n_z):172 y = P.ybar[j - 1] * P.z_grid[iz]173 for ix in range(P.n_x):174 x = P.x_grid[ix]175 wb, wr = pol.W_B[iz, ix], pol.W_Rs[iz, ix]176 rows.append(dict(age=j, z_index=iz, z=float(P.z_grid[iz]), cash=x, a_equiv=(x - y) / pr.R_a,177 W_buy=wb if wb > NEG / 2 else np.nan, W_rent=wr if wr > NEG / 2 else np.nan,178 buy=int((wb > wr) and wb > NEG / 2), h_buy=int(pol.hB[iz, ix]), B_buy=float(pol.BB[iz, ix]),179 a_buy=float(pol.aB[iz, ix]), a_rent=float(pol.aRs[iz, ix]), gift_buy=float(pol.gB[iz, ix]),180 gift_rent=float(pol.gRs[iz, ix]),181 c_rent=float(max(x - pol.aRs[iz, ix] - pol.gRs[iz, ix], 0) * (1 - P.alpha))))182 out[f"policy_age{j}"] = pd.DataFrame(rows)183 return out184185186def thresholds_table(P: M.Params, res: E.SSResult) -> pd.DataFrame:187 rows = []188 for j in (1, 2, 3):189 th = ownership_threshold(P, res.ev.pols[j - 1], res.ev.prices, j)190 for iz, v in th.items():191 rows.append(dict(age=j, z_index=iz, **v))192 return pd.DataFrame(rows)193194195def summarize_ss(P: M.Params, res: E.SSResult, name: str, group: str, scenario: dict, t_sec: float) -> dict:196 rec = res.ev.rec197 mom = C.model_moments(P, rec, res.p, res.rho)198 mob = mobility_stats(P, rec)199 nw = rec.net_worth1(P)200 own = rec.owner1()201 # owner-renter gaps at fixed ages (equal income by design: same income process for all)202 gaps = {}203 for j in (2, 4, 7, 10):204 m = rec.age == j205 gaps[f"nw_gap_owner_renter_age{j}"] = float(nw[m & own].mean() - nw[m & ~own].mean()) if (m & ~own).any() and (m & own).any() else np.nan206 z1 = rec.z[rec.age == 1]207 return dict(name=name, group=group, scenario=scenario, equilibrium=res.equilibrium, converged=bool(res.converged),208 iterations=int(res.iterations), time_sec=float(t_sec), p=res.p, p_r=res.p_r, ptil_r=res.ptil_r, rho=res.rho,209 rent_price_annual=res.rho / res.p_r / P.period_years, H=res.ev.H_d, Hr=res.ev.S_d,210 H_bar=P.H_bar, Hr_bar=P.Hr_bar, lottery=res.lottery[0], rebate=res.rebate,211 residuals=res.residuals, excess_owner=float(res.residuals[0]), excess_rental=float(res.residuals[1]),212 bellman=res.bellman, forced_share=res.ev.forced_share, clip_share=res.ev.clip_share,213 moments=mom, mobility=mob, gaps=gaps,214 thresholds={j: {iz: v["a_hat"] for iz, v in th.items()} for j, th in res.thresholds.items()},215 threshold_sign_changes={j: {iz: v["sign_changes"] for iz, v in th.items()} for j, th in res.thresholds.items()},216 threshold_htm={j: {iz: [v["htm_buyer"], v["htm_renter"]] for iz, v in th.items()} for j, th in res.thresholds.items()},217 mean_nw=float(nw.mean()), gini_nw=D.gini(nw), total_debt=float(rec.b1.sum() / (len(rec.age) / P.J)),218 total_deposits=float(rec.a1.sum() / (len(rec.age) / P.J)),219 params=dict(theta=P.theta, psi=P.psi, eps=P.eps, eps_r=P.eps_r, nu=P.nu, tau_b=P.tau_b, r=P.r, pibar=P.pibar,220 phi=P.phi, sigma_seg=P.sigma_seg, sigma_z=P.sigma_z, beta=P.beta, gamma=P.gamma, alpha=P.alpha,221 kappa_b=P.kappa_b, kappa_s=P.kappa_s, chi=P.chi, n_a=P.n_a, n_b=P.n_b, n_x=P.n_x, n_dyn=P.n_dyn,222 seed=P.seed, dep_premium=P.dep_premium))223224225def save_cross_section(P: M.Params, rec: D.PeriodRecord, path: Path, subsample: int = 1):226 path.parent.mkdir(parents=True, exist_ok=True)227 sl = slice(None, None, subsample)228 np.savez_compressed(path, age=rec.age[sl], z=rec.z[sl], h0=rec.h0[sl], b0=rec.b0[sl], a0=rec.a0[sl], choice=rec.choice[sl],229 c=rec.c[sl], s=rec.s[sl], h1=rec.h1[sl], b1=rec.b1[sl], a1=rec.a1[sl], gift=rec.gift[sl], estate=rec.estate[sl],230 gift_received=rec.gift_received[sl], parent_owner=rec.parent_owner[sl], parent_nw=rec.parent_nw[sl],231 estate_received=rec.estate_received[sl], is_old=rec.is_old[sl], nw0=rec.net_worth0(P)[sl],232 nw1=rec.net_worth1(P)[sl], pay=rec.pay[sl], Bnew=rec.Bnew[sl], flag=rec.flag[sl], u=rec.u[sl],233 p=rec.p, rho=rec.rho, h_grid=P.h_grid, z_grid=P.z_grid)234235236# --------------------------------------------------------------------------------------237# Experiment 2: equal-income deterministic paths238# --------------------------------------------------------------------------------------239def equal_income_paths(P: M.Params, res: E.SSResult) -> pd.DataFrame:240 """Households with the same (median) income path and different initial transfers; policies241 of the baseline equilibrium; estate at j_estate = the lottery's median point (deterministic).242 Also the two households straddling the age-1 ownership threshold (owner vs renter at equal243 income and almost equal endowment), and a renter who invests the down payment at R_a."""244 S = D.stack_policies(res.ev.pols)245 pr = res.ev.prices246 iz = P.n_z // 2247 g1 = res.ev.rec.gift_received[(res.ev.rec.age == 1) & ~res.ev.rec.is_old]248 qs = {"q10": 0.1, "q25": 0.25, "q50": 0.5, "q75": 0.75, "q90": 0.9}249 a1_vals = {k: float(np.quantile(g1, v)) for k, v in qs.items()}250 th = ownership_threshold(P, res.ev.pols[0], pr, 1)[iz]["a_hat"]251 if np.isfinite(th):252 a1_vals["threshold_minus"] = max(th - 0.01, 0.0)253 a1_vals["threshold_plus"] = th + 0.01254 names = list(a1_vals.keys())255 a1 = np.array([a1_vals[k] for k in names])256 z1 = np.full(len(a1), iz)257 med_lot = (np.array([float(np.median(res.lottery[0]))]), np.array([1.0]))258 out = D.simulate_cohort(P, S, pr, z1, a1, med_lot, fixed_z=True)259 rows = []260 for k, nm in enumerate(names):261 for j in range(1, P.J + 1):262 rows.append(dict(path=nm, a1=float(a1[k]), age=j, owner=int(out["own"][j - 1, k]), nw=float(out["nw"][j - 1, k]),263 liquid=float(out["a"][j - 1, k]), debt=float(out["b"][j - 1, k]), c=float(out["c"][j - 1, k]),264 s=float(out["s"][j - 1, k]), h=int(out["h"][j - 1, k]), gift_given=float(out["gift"][k]),265 lifetime_V=float(out["V"][k])))266 return pd.DataFrame(rows)267268269# --------------------------------------------------------------------------------------270# Experiment 5 (accounting part): unanticipated price drop on the cross-section271# --------------------------------------------------------------------------------------272def price_drop_accounting(P: M.Params, rec: D.PeriodRecord, drops=(0.1, 0.2, 0.3)) -> pd.DataFrame:273 own = rec.owner1()274 hv = np.where(own, P.h_grid[np.maximum(rec.h1, 0)], 0.0) * rec.p275 nw = rec.net_worth1(P)276 rows = []277 for x in drops:278 eq_after = (1 - P.kappa_s) * (1 - x) * hv - rec.b1279 loss = own * (1 - P.kappa_s) * x * hv280 for j in list(range(1, P.J + 1)) + ["all"]:281 m = own if j == "all" else (own & (rec.age == j))282 if m.sum() == 0:283 continue284 ltv = rec.b1[m] / hv[m]285 rows.append(dict(drop=x, age=j, n_owners=int(m.sum()), underwater_share=float((eq_after[m] < 0).mean()),286 equity_wiped_share=float((eq_after[m] <= 0.0).mean()),287 mean_loss_over_nw=float((loss[m] / np.maximum(nw[m], 1e-6)).mean()),288 mean_loss_over_income=float(loss[m].mean() / (P.ybar[min(j, P.J) - 1] if j != "all" else 1.0)),289 mean_ltv=float(ltv.mean()), share_ltv_above_x_star=float((ltv > (1 - P.kappa_s) * (1 - x)).mean())))290 # new borrowers evaluated at origination (before the first payment): LTV = Bnew / (p h)291 m = rec.Bnew > 0292 if m.any():293 ltv0 = rec.Bnew[m] / hv[m]294 rows.append(dict(drop=x, age="new_loans_at_origination", n_owners=int(m.sum()),295 underwater_share=float((((1 - P.kappa_s) * (1 - x) * hv[m]) < rec.Bnew[m]).mean()),296 equity_wiped_share=float((((1 - P.kappa_s) * (1 - x) * hv[m]) <= rec.Bnew[m]).mean()),297 mean_loss_over_nw=float(((1 - P.kappa_s) * x * hv[m] / np.maximum(nw[m], 1e-6)).mean()),298 mean_loss_over_income=float(((1 - P.kappa_s) * x * hv[m]).mean()),299 mean_ltv=float(ltv0.mean()), share_ltv_above_x_star=float((ltv0 > (1 - P.kappa_s) * (1 - x)).mean())))300 return pd.DataFrame(rows)301302303# --------------------------------------------------------------------------------------304# Experiment 9: consumption-equivalent welfare across stationary equilibria305# --------------------------------------------------------------------------------------306def welfare_vs_baseline(P_base: M.Params, base: E.SSResult, P_alt: M.Params, alt: E.SSResult, name: str) -> pd.DataFrame:307 """CEV at age 1 by (z_1, baseline gift quintile): lambda such that scaling non-housing308 consumption by (1 + lambda) in every period of the BASELINE makes the age-1 household as309 well off as under the alternative, at the same initial state (z_1, a_1).310 V_base(lambda) = (1+lambda)^((1-alpha)(1-gamma)) U_base + Phi_base, where U_base (flow311 utility) and Phi_base (warm glow) are obtained by cohort simulation under baseline policies;312 V_alt is the alternative's age-1 value function at (z_1, a_1).313 Note: the comparison is between stationary equilibria (not along a transition)."""314 S_b = D.stack_policies(base.ev.pols)315 g1 = base.ev.rec.gift_received[(base.ev.rec.age == 1) & ~base.ev.rec.is_old]316 edges = np.quantile(g1, [0.2, 0.4, 0.6, 0.8])317 rows = []318 expo = (1 - P_base.alpha) * (1 - P_base.gamma)319 n_per = 3000320 rng = np.random.default_rng(5)321 th = {iz: v["a_hat"] for iz, v in ownership_threshold(P_base, base.ev.pols[0], base.ev.prices, 1).items()}322 for k in range(6):323 if k < 5:324 lo = 0.0 if k == 0 else edges[k - 1]; hi = edges[k] if k < 4 else np.quantile(g1, 0.995)325 a1 = rng.uniform(lo, max(hi, lo + 1e-6), size=n_per) if hi > lo else np.full(n_per, lo)326 for iz in range(P_base.n_z):327 if k == 5:328 # marginal households: within +/- 10 % of the age-1 ownership threshold for this z329 if not np.isfinite(th[iz]) or th[iz] <= 0:330 continue331 a1 = rng.uniform(0.9 * th[iz], 1.1 * th[iz], size=n_per)332 z1 = np.full(n_per, iz)333 coh = D.simulate_cohort(P_base, S_b, base.ev.prices, z1, a1, base.lottery, seed=100 + k)334 U_sim, Phi = coh["U"].mean(), coh["Phi"].mean()335 Vb_sim = U_sim + Phi336 Vb = np.mean(np.interp(a1, P_base.a_grid, base.ev.pols[0].V_rent[iz]))337 Va = np.mean(np.interp(a1 + alt.rebate - base.rebate, P_alt.a_grid, alt.ev.pols[0].V_rent[iz]))338 # consistent decomposition: the flow part is V_base - Phi_base (Phi from the simulation),339 # so that the identity holds exactly at the baseline; U_sim is kept as a check340 U = Vb - Phi341 ratio = (Va - Phi) / U342 cev = ratio ** (1.0 / expo) - 1.0 if ratio > 0 else np.nan343 rows.append(dict(scenario=name, gift_quintile=(k + 1) if k < 5 else "marginal", z_index=iz, z=float(P_base.z_grid[iz]), a1_mean=float(a1.mean()),344 V_base=float(Vb), V_base_simulated=float(Vb_sim), V_alt=float(Va), U_base=float(U), U_base_simulated=float(U_sim), Phi_base=float(Phi),345 cev=float(cev), weight=float(P_base.pi_z[iz] / 5)))346 return pd.DataFrame(rows)347348349# --------------------------------------------------------------------------------------350# transition summaries351# --------------------------------------------------------------------------------------352def transition_path_table(P: M.Params, tr: E.TransitionResult, ss: E.SSResult) -> pd.DataFrame:353 rows = []354 for t, rec in enumerate(tr.recs):355 nw = rec.net_worth1(P); own = rec.owner1()356 young = rec.age <= 3357 rows.append(dict(t=t, years=5 * t, p=float(tr.p[t]), p_r=float(tr.p_r[t]), rho=float(tr.rho[t]),358 p_pct=100 * (tr.p[t] / ss.p - 1), rho_pct=100 * (tr.rho[t] / ss.rho - 1),359 H=float(tr.H[t]), Hr=float(tr.Hr[t]), ownership=float(own.mean()), ownership_young=float(own[young].mean()),360 buy_rate_young=float((rec.choice[young] == 2).mean()),361 excess_owner=float(tr.residuals[t]), excess_rental=float(tr.residuals[tr.T + t]),362 debt_total=float(rec.b1.sum() / (len(rec.age) / P.J)), deposits_total=float(rec.a1.sum() / (len(rec.age) / P.J)),363 nw_mean=float(nw.mean()), nw_owner_mean=float(nw[own].mean()), nw_renter_mean=float(nw[~own].mean()),364 gini_nw=D.gini(nw), c_mean=float(rec.c.mean()),365 i_new=float(tr.prices[t].i_new), R_a=float(tr.prices[t].R_a), pi_next=float(tr.prices[t].pi_next)))366 return pd.DataFrame(rows)367368369def impact_table(P: M.Params, tr: E.TransitionResult, ss: E.SSResult, pe: E.TransitionResult | None) -> pd.DataFrame:370 """Impact-period (t = 0) redistribution by age and start-of-period tenure: real revaluation371 of nominal positions (identity, PE), consumption and net-worth changes vs the stationary372 path, with and without price adjustment."""373 rec0 = tr.recs[0]374 pr0 = tr.prices[0]375 pi_u = tr.paths["pi_u"]376 ss_rec = ss.ev.rec377 # stationary comparison period: same households at t=0 under SS policies and prices (the378 # panel state is identical, only the shock differs) -> simulate one SS period from ss.ev.panel379 S_ss = D.stack_policies(ss.ev.pols)380 panel = ss.ev.panel.copy()381 rec_ss = D.simulate_period(P, S_ss, ss.ev.prices, panel, (P.seed, 100000))382 rows = []383 own0 = rec0.h0 >= 0384 b_pre = rec_ss.b0; a_pre = rec_ss.a0385 bf = float(tr.paths["balance_factor"])386 # revaluation gains: debt erosion (1 - bf) * b_pre ; deposit loss a_pre * (R - R_a0)387 debt_gain = (1 - bf) * b_pre388 dep_loss = a_pre * (P.R * (1 + P.dep_premium) - pr0.R_a)389 wage_loss = (P.ybar[rec0.age - 1] * P.z_grid[rec0.z]) * (1 - pr0.wage_factor)390 rent_gain = np.where(rec0.h0 < 0, rec_ss.s * ss.rho * (1 - pr0.rent_rigid_factor), 0.0)391 net_nominal = debt_gain - dep_loss392 for grp_name, grp in [("all", np.ones(len(rec0.age), bool)), ("owners", own0), ("renters", ~own0)]:393 for j in list(range(1, P.J + 1)) + ["all"]:394 m = grp if j == "all" else (grp & (rec0.age == j))395 if m.sum() == 0:396 continue397 rows.append(dict(scenario=tr.shock.contract, group=grp_name, age=j, n=int(m.sum()),398 debt_erosion_gain=float(debt_gain[m].mean()), deposit_loss=float(dep_loss[m].mean()),399 wage_loss=float(wage_loss[m].mean()), rent_rigidity_gain=float(rent_gain[m].mean()),400 net_nominal_gain=float(net_nominal[m].mean()),401 c_ss=float(rec_ss.c[m].mean()),402 d_consumption=float(rec0.c[m].mean() - rec_ss.c[m].mean()),403 d_consumption_pct=float(100 * (rec0.c[m].mean() / rec_ss.c[m].mean() - 1)),404 d_consumption_pe_pct=float(100 * (pe.recs[0].c[m].mean() / rec_ss.c[m].mean() - 1)) if pe is not None else np.nan,405 d_nw=float(rec0.net_worth1(P)[m].mean() - rec_ss.net_worth1(P)[m].mean()),406 d_ownership=float(rec0.owner1()[m].mean() - rec_ss.owner1()[m].mean()),407 d_consumption_pe=float(pe.recs[0].c[m].mean() - rec_ss.c[m].mean()) if pe is not None else np.nan,408 d_nw_pe=float(pe.recs[0].net_worth1(P)[m].mean() - rec_ss.net_worth1(P)[m].mean()) if pe is not None else np.nan,409 d_ownership_pe=float(pe.recs[0].owner1()[m].mean() - rec_ss.owner1()[m].mean()) if pe is not None else np.nan))410 # lender side (per unit mass of one age): mirror image of households' net nominal gain411 mass = len(rec0.age) / P.J412 lender = dict(scenario=tr.shock.contract, group="lender", age="all", n=int(len(rec0.age)),413 debt_erosion_gain=-float(debt_gain.sum() / mass), deposit_loss=-float(dep_loss.sum() / mass),414 wage_loss=0.0, rent_rigidity_gain=0.0, net_nominal_gain=-float(net_nominal.sum() / mass),415 c_ss=np.nan, d_consumption=np.nan, d_consumption_pct=np.nan, d_consumption_pe_pct=np.nan, d_nw=np.nan, d_ownership=np.nan,416 d_consumption_pe=np.nan, d_nw_pe=np.nan, d_ownership_pe=np.nan)417 rows.append(lender)418 return pd.DataFrame(rows)419420421def transition_welfare(P: M.Params, tr: E.TransitionResult, ss: E.SSResult, pols0: list) -> pd.DataFrame:422 """Group-level CEV for households alive at the shock: V_{j,0}(state after the shock) vs423 the stationary V_j(state); Phi_base (warm-glow part) by remaining-lifetime simulation under424 the stationary policies, U_base = V_ss - Phi_base (consistent flow part; the simulated U is425 reported as a check). lambda_group = ((mean V_alt - mean Phi_base)/mean U_base)^(1/expo) - 1."""426 rec0 = tr.recs[0]427 S_ss = D.stack_policies(ss.ev.pols)428 U, Phi = D.remaining_lifetime(P, S_ss, ss.ev.prices, rec0, ss.lottery)429 # stationary and post-shock values at the start-of-period state (interpolated)430 bf = float(tr.paths["balance_factor"])431 V_ss = _value_at_states(P, ss.ev.pols, rec0, b_scale=1.0 / bf if bf > 0 else 1.0) # pre-shock balances432 V_alt = _value_at_states(P, pols0, rec0, b_scale=1.0)433 expo = (1 - P.alpha) * (1 - P.gamma)434 rows = []435 own0 = rec0.h0 >= 0436 hasdebt = rec0.b0 > 1e-8437 for grp_name, grp in [("all", np.ones(len(rec0.age), bool)), ("owners", own0), ("owners_with_mortgage", own0 & hasdebt),438 ("owners_no_mortgage", own0 & ~hasdebt), ("renters", ~own0)]:439 for j in list(range(1, P.J + 1)) + ["all"]:440 m = grp if j == "all" else (grp & (rec0.age == j))441 ok = m & np.isfinite(V_ss) & np.isfinite(V_alt) & (V_ss > NEG / 2) & (V_alt > NEG / 2)442 if ok.sum() < 20:443 continue444 Ub_sim, Pb = U[ok].mean(), Phi[ok].mean()445 Vs = V_ss[ok].mean()446 Ub = Vs - Pb # consistent flow part (identity exact at the baseline)447 ratio = (V_alt[ok].mean() - Pb) / Ub448 cev = ratio ** (1 / expo) - 1 if ratio > 0 else np.nan449 rows.append(dict(scenario=tr.shock.contract, group=grp_name, age=j, n=int(ok.sum()), V_ss=float(Vs),450 V_shock=float(V_alt[ok].mean()), U_base=float(Ub), U_base_simulated=float(Ub_sim), Phi_base=float(Pb),451 V_check_sim_vs_fn=float((Ub_sim + Pb) / Vs - 1), cev=float(cev),452 share_better_off=float((V_alt[ok] > V_ss[ok]).mean())))453 return pd.DataFrame(rows)454455456def _value_at_states(P: M.Params, pols: list, rec: D.PeriodRecord, b_scale: float = 1.0) -> np.ndarray:457 V = np.full(len(rec.age), np.nan)458 for j in range(1, P.J + 1):459 m = rec.age == j460 pol = pols[j - 1]461 for iz in range(P.n_z):462 mm = m & (rec.z == iz)463 if not mm.any():464 continue465 r = mm & (rec.h0 < 0)466 if r.any():467 V[r] = np.interp(rec.a0[r], P.a_grid, pol.V_rent[iz])468 o = mm & (rec.h0 >= 0)469 if o.any():470 idx = np.where(o)[0]471 for i in idx:472 ih, n_idx = rec.h0[i], rec.n0[i] - 1473 bb = rec.b0[i] * b_scale474 jb, wb = _loc(P.b_grid, bb); ja, wa = _loc(P.a_grid, rec.a0[i])475 v = pol.V_own[iz, ih, jb:jb + 2, n_idx, ja:ja + 2]476 if v.min() <= NEG / 2:477 V[i] = NEG478 else:479 V[i] = (1 - wb) * ((1 - wa) * v[0, 0] + wa * v[0, 1]) + wb * ((1 - wa) * v[1, 0] + wa * v[1, 1])480 return V481482483def _loc(grid, x):484 n = len(grid)485 if x <= grid[0]:486 return 0, 0.0487 if x >= grid[-1]:488 return n - 2, 1.0489 k = int(np.searchsorted(grid, x, side="right")) - 1490 return k, (x - grid[k]) / (grid[k + 1] - grid[k])491492493def summarize_transition(P: M.Params, tr: E.TransitionResult, ss: E.SSResult, name: str, scen: dict, path_df: pd.DataFrame,494 impact_df: pd.DataFrame, welfare_df: pd.DataFrame | None, pe: E.TransitionResult | None) -> dict:495 rec0 = tr.recs[0]496 mass = len(rec0.age) / P.J497 bf = float(tr.paths["balance_factor"])498 S_ss = D.stack_policies(ss.ev.pols)499 panel = ss.ev.panel.copy()500 rec_ss = D.simulate_period(P, S_ss, ss.ev.prices, panel, (P.seed, 100000))501 debt_gain = (1 - bf) * rec_ss.b0; dep_loss = rec_ss.a0 * (P.R * (1 + P.dep_premium) - tr.prices[0].R_a)502 # mobility along the path: parent-child tenure at age 7 in the first generation after the shock503 mob = []504 for t, rec in enumerate(tr.recs):505 mob.append(mobility_stats(P, rec))506 return dict(name=name, group="inflation", scenario=scen, equilibrium=tr.equilibrium, converged=bool(tr.converged),507 iterations=int(tr.iterations), time_sec=float(tr.time_sec), max_residual=float(np.max(np.abs(tr.residuals))),508 T=tr.T, pi_u=float(tr.paths["pi_u"]), contract=tr.shock.contract, persistence=tr.shock.persistence,509 balance_factor=bf, p_impact_pct=float(100 * (tr.p[0] / ss.p - 1)), rho_impact_pct=float(100 * (tr.rho[0] / ss.rho - 1)),510 p_peak_pct=float(100 * (np.max(np.abs(tr.p / ss.p - 1)) * np.sign(tr.p[np.argmax(np.abs(tr.p / ss.p - 1))] - ss.p))),511 p_path_pct=(100 * (tr.p / ss.p - 1)).tolist(), rho_path_pct=(100 * (tr.rho / ss.rho - 1)).tolist(),512 ownership_impact=float(rec0.owner1().mean()), ownership_ss=float(rec_ss.owner1().mean()),513 d_ownership_impact_pp=float(100 * (rec0.owner1().mean() - rec_ss.owner1().mean())),514 d_ownership_impact_pe_pp=float(100 * (pe.recs[0].owner1().mean() - rec_ss.owner1().mean())) if pe is not None else None,515 aggregate_debt_erosion=float(debt_gain.sum() / mass), aggregate_deposit_loss=float(dep_loss.sum() / mass),516 lender_revaluation=-float((debt_gain - dep_loss).sum() / mass),517 # lender side computed from its balance sheet (M - D) x erosion: must equal minus the household sum518 zero_sum_check=float(((1 - bf) * rec_ss.b0.sum() - (P.R * (1 + P.dep_premium) - tr.prices[0].R_a) * rec_ss.a0.sum()) / mass519 + (-(debt_gain - dep_loss).sum() / mass)),520 p_path_sign_changes=int(np.sum(np.sign(tr.p[1:] / ss.p - 1) != np.sign(tr.p[:-1] / ss.p - 1))),521 residual_at_T=float(max(abs(tr.residuals[tr.T - 1]), abs(tr.residuals[2 * tr.T - 1]))),522 mobility_path=mob, mobility_ss=mobility_stats(P, ss.ev.rec),523 welfare_all=float(welfare_df.loc[(welfare_df.group == "all") & (welfare_df.age == "all"), "cev"].iloc[0]) if welfare_df is not None else None,524 path=path_df.to_dict(orient="list"))525526527# --------------------------------------------------------------------------------------528# runners529# --------------------------------------------------------------------------------------530def run_baseline(cfg: dict, quick: bool = False) -> tuple[M.Params, E.SSResult, dict]:531 t0 = time.time()532 P = M.build_params(cfg)533 log = E.RunLog("baseline")534 log("start", quick=quick, n_dyn=P.n_dyn)535 H_bar, Hr_bar, cal = E.calibrate_supply(P, log)536 log("calibrated_supply", H_bar=H_bar, Hr_bar=Hr_bar)537 res = E.solve_stationary(P, p0=P.p_ref, ptil0=P.p_ref, lottery0=cal.lottery, log=log, name="baseline")538 out = DATA / "baseline"539 out.mkdir(parents=True, exist_ok=True)540 summ = summarize_ss(P, res, "baseline", "baseline", {}, time.time() - t0)541 summ["supply_scales"] = dict(H_bar=H_bar, Hr_bar=Hr_bar, calibration_iterations=cal.iterations)542 save_json(out / "summary.json", summ)543 save_json(out / "supply_scales.json", dict(H_bar=H_bar, Hr_bar=Hr_bar, lottery=res.lottery[0].tolist(), p=res.p, ptil_r=res.ptil_r))544 rec = res.ev.rec545 profiles_by_age(P, rec).to_csv(out / "profiles_by_age.csv", index=False)546 by_gift_quintile(P, rec).to_csv(out / "by_gift_quintile.csv", index=False)547 save_json(out / "mobility.json", mobility_stats(P, rec))548 for k, df in policy_tables(P, res).items():549 df.to_csv(out / f"{k}.csv", index=False)550 thresholds_table(P, res).to_csv(out / "thresholds.csv", index=False)551 save_cross_section(P, rec, out / "cross_section.npz", subsample=1)552 equal_income_paths(P, res).to_csv(out / "equal_income_paths.csv", index=False)553 price_drop_accounting(P, rec).to_csv(out / "price_drop_accounting.csv", index=False)554 pd.DataFrame(M.parameter_table(P)).to_csv(DATA / "parameters.csv", index=False)555 C.targets_table(P, summ["moments"]).to_csv(DATA / "calibration_targets.csv", index=False)556 # income process557 pd.DataFrame(dict(z=P.z_grid, pi=P.pi_z)).to_csv(out / "income_grid.csv", index=False)558 pd.DataFrame(P.Pi_z).to_csv(out / "income_transition.csv", index=False)559 log("baseline_done", sec=time.time() - t0, p=res.p, converged=res.converged)560 with open(out / "ss_result.pkl", "wb") as fh:561 pickle.dump(dict(P=P, res=res), fh)562 return P, res, summ563564565def run_ss_scenario(args) -> dict:566 """Worker: one stationary scenario (GE or PE at baseline prices)."""567 cfg_base, scen, group, scales, base_lottery, base_prices, n_threads = args568 import numba569 numba.set_num_threads(max(1, n_threads))570 t0 = time.time()571 name = scen["name"]572 cfg = apply_overrides(cfg_base, scen.get("overrides"))573 P = M.build_params(cfg)574 P.H_bar, P.Hr_bar = scales["H_bar"], scales["Hr_bar"]575 log = E.RunLog(name, echo=False)576 kind = scen["kind"]577 K = P.K_estate578 lot0 = (np.array(base_lottery)[:K], np.full(K, 1.0 / K)) if len(base_lottery) >= K else None579 try:580 if kind == "ss_pe":581 res = E.solve_stationary(P, p0=base_prices["p"], ptil0=base_prices["ptil_r"], lottery0=lot0, fixed_prices=True,582 log=log, name=name)583 else:584 res = E.solve_stationary(P, p0=base_prices["p"], ptil0=base_prices["ptil_r"], lottery0=lot0, log=log, name=name)585 summ = summarize_ss(P, res, name, group, scen, time.time() - t0)586 out = DATA / group587 out.mkdir(parents=True, exist_ok=True)588 save_json(out / f"{name}.json", summ)589 profiles_by_age(P, res.ev.rec).to_csv(out / f"{name}_profiles.csv", index=False)590 by_gift_quintile(P, res.ev.rec).to_csv(out / f"{name}_by_gift_quintile.csv", index=False)591 thresholds_table(P, res).to_csv(out / f"{name}_thresholds.csv", index=False)592 if name.startswith("theta") or name.startswith("eps") or name in ("bequests_off", "noseg"):593 save_cross_section(P, res.ev.rec, out / f"{name}_cross_section.npz", subsample=4)594 # welfare vs baseline (experiment 9) for GE scenarios595 with open(out / f"{name}_ss.pkl", "wb") as fh:596 pickle.dump(dict(P=P, core=E.ss_core(res), pols0=res.ev.pols[0], lottery=res.lottery, rebate=res.rebate,597 p=res.p, ptil_r=res.ptil_r), fh)598 log("scenario_done", sec=time.time() - t0, converged=res.converged, p=res.p)599 return dict(name=name, group=group, ok=True, converged=bool(res.converged), p=res.p, sec=time.time() - t0,600 equilibrium=res.equilibrium, residuals=res.residuals.tolist(), bellman=res.bellman)601 except Exception as exc: # noqa: BLE001 — report computational failures, never hide them602 import traceback603 log("scenario_failed", error=repr(exc), trace=traceback.format_exc())604 return dict(name=name, group=group, ok=False, error=repr(exc), sec=time.time() - t0)605606607def run_welfare(cfg_base: dict, groups: dict, base_P: M.Params, base_res: E.SSResult):608 """Experiment 9: CEV by gift quintile and z_1 for every converged GE scenario with a saved pickle."""609 out = DATA / "welfare"610 out.mkdir(parents=True, exist_ok=True)611 frames = []612 for group, scens in groups.items():613 for scen in scens:614 if scen["kind"] != "ss" or scen["name"] == "baseline":615 continue616 pkl = DATA / group / f"{scen['name']}_ss.pkl"617 if not pkl.exists():618 continue619 with open(pkl, "rb") as fh:620 d = pickle.load(fh)621 P_alt = d["P"]622 # light-weight alternative result object with age-1 policies only623624 class _Alt:625 pass626 alt = _Alt(); alt.rebate = d["rebate"]627 alt.ev = _Alt(); alt.ev.pols = [d["pols0"]]628 if P_alt.n_a != base_P.n_a:629 continue # grid variants: value functions not comparable on the same grid630 df = welfare_vs_baseline(base_P, base_res, P_alt, alt, scen["name"])631 df.to_csv(out / f"welfare_{scen['name']}.csv", index=False)632 frames.append(df)633 if frames:634 pd.concat(frames).to_csv(out / "welfare_all.csv", index=False)635636637def run_transition_scenario(cfg_base: dict, scen: dict, scales: dict, base: dict, n_workers: int,638 warm_from: str | None = None) -> dict:639 """warm_from: name of a converged transition whose (p, p~_r) path initialises the shooting."""640 t0 = time.time()641 name = scen["name"]642 cfg = apply_overrides(cfg_base, scen.get("overrides"))643 shk = scen.get("shock", {})644 cfg["mortgage"]["contract"] = shk.get("contract", "frm")645 P = M.build_params(cfg)646 P.H_bar, P.Hr_bar = scales["H_bar"], scales["Hr_bar"]647 log = E.RunLog(name)648 T = int(P.cfg["transition"]["T"])649 # stationary equilibrium of this parameterisation (baseline unless overrides change the SS)650 ss_overrides = {k: v for k, v in (scen.get("overrides") or {}).items() if not k.startswith("rental.eta")}651 if ss_overrides:652 lot0 = (np.array(base["lottery"]), np.full(P.K_estate, 1.0 / P.K_estate))653 ss = E.solve_stationary(P, p0=base["p"], ptil0=base["ptil_r"], lottery0=lot0, log=log, name=name, diagnostics=False)654 else:655 with open(DATA / "baseline" / "ss_result.pkl", "rb") as fh:656 d = pickle.load(fh)657 ss = d["res"]658 ss.P = P # same parameters except the contract flag / indexation (SS-neutral)659 kind = shk.get("type", "inflation")660 sh = E.Shock(kind=kind, pi_u=float(P.cfg["transition"]["pi_surprise"]), persistence=int(shk.get("persistence", 1)),661 contract=shk.get("contract", "frm"), phi_r=float(shk.get("phi_r", 0.0)),662 dr=M._to_period(float(shk.get("dr_annual", 0.0)), P.period_years))663 core = E.ss_core(ss)664 mode = "static" if scen["kind"] == "static" else "GE"665 pe = E.solve_transition(P, core, sh, T=T, mode="PE", log=log, name=name)666 u0 = None667 if warm_from is not None and (DATA / "inflation" / f"{warm_from}_path.csv").exists():668 wp = pd.read_csv(DATA / "inflation" / f"{warm_from}_path.csv")669 if len(wp) == T:670 # scale the neighbour's path to this experiment's own stationary prices671 with open(DATA / "inflation" / f"{warm_from}.json") as fh:672 wj = json.load(fh)673 u0 = np.concatenate([np.log(ss.p * (1 + np.array(wj["p_path_pct"]) / 100)),674 np.log(ss.ptil_r * (wp.p_r.values / wp.p_r.values[-1]))])675 tr = E.solve_transition(P, core, sh, T=T, mode=mode, log=log, name=name, n_workers=n_workers, u0=u0)676 out = DATA / "inflation"677 out.mkdir(parents=True, exist_ok=True)678 path_df = transition_path_table(P, tr, ss)679 path_df.to_csv(out / f"{name}_path.csv", index=False)680 pe_df = transition_path_table(P, pe, ss)681 pe_df.to_csv(out / f"{name}_path_pe.csv", index=False)682 imp = impact_table(P, tr, ss, pe)683 imp.to_csv(out / f"{name}_impact.csv", index=False)684 # value functions at t = 0 come from the converged path: re-solve the household block once685 prices = tr.prices686 pols_by_t = E.solve_households_path(P, prices, core, T, static_expectations=(mode == "static"))687 wel = transition_welfare(P, tr, ss, pols_by_t[0])688 wel.to_csv(out / f"{name}_welfare.csv", index=False)689 # ownership-entry threshold at impact vs stationary (Proposition 9.3)690 thr0 = {j: {iz: v["a_hat"] for iz, v in ownership_threshold(P, pols_by_t[0][j - 1], prices[0], j).items()} for j in (1, 2, 3)}691 thr_ss = {j: {iz: v["a_hat"] for iz, v in ownership_threshold(P, ss.ev.pols[j - 1], ss.ev.prices, j).items()} for j in (1, 2, 3)}692 pe_pols = E.solve_households_path(P, pe.prices, core, T)693 wel_pe = transition_welfare(P, pe, ss, pe_pols[0])694 wel_pe.to_csv(out / f"{name}_welfare_pe.csv", index=False)695 summ = summarize_transition(P, tr, ss, name, scen, path_df, imp, wel, pe)696 summ["threshold_impact"] = thr0697 summ["threshold_ss"] = thr_ss698 summ["pe"] = dict(equilibrium=pe.equilibrium, max_residual=float(np.max(np.abs(pe.residuals))),699 d_ownership_impact_pp=summ["d_ownership_impact_pe_pp"],700 welfare_all=float(wel_pe.loc[(wel_pe.group == "all") & (wel_pe.age == "all"), "cev"].iloc[0]))701 summ["time_sec"] = time.time() - t0702 save_json(out / f"{name}.json", summ)703 log("transition_scenario_done", sec=time.time() - t0, converged=tr.converged, p_impact_pct=summ["p_impact_pct"])704 return dict(name=name, group="inflation", ok=True, converged=bool(tr.converged), equilibrium=tr.equilibrium,705 max_residual=float(np.max(np.abs(tr.residuals))), sec=time.time() - t0)706707708def _transition_job(args) -> dict:709 cfg, scen, scales, base, n_workers = args710 import numba711 numba.set_num_threads(max(2, n_workers))712 try:713 return run_transition_scenario(cfg, scen, scales, base, n_workers=n_workers)714 except Exception as exc: # noqa: BLE001 — computational failures are reported, never hidden715 import traceback716 return dict(name=scen["name"], group="inflation", ok=False, error=repr(exc), trace=traceback.format_exc())717718719def run_rent_growth_pe(cfg_base: dict, scen: dict, scales: dict, n_workers: int) -> dict:720 """Experiment 8 (rent growth): anticipated real rent growth g for T periods at fixed house721 prices (PE): households solve with the time-varying rent path; we report the tenure and722 wealth response by age. Labelled PE (prices fixed)."""723 t0 = time.time()724 name = scen["name"]725 cfg = apply_overrides(cfg_base, scen.get("overrides"))726 P = M.build_params(cfg)727 P.H_bar, P.Hr_bar = scales["H_bar"], scales["Hr_bar"]728 log = E.RunLog(name)729 with open(DATA / "baseline" / "ss_result.pkl", "rb") as fh:730 d = pickle.load(fh)731 ss = d["res"]; ss.P = P732 core = E.ss_core(ss)733 T = int(P.cfg["transition"]["T"])734 g = M._to_period(float(scen["shock"]["g_rent_annual"]), P.period_years)735 sh = E.Shock(kind="none")736 paths = E.shock_paths(P, sh, T)737 logp = np.full(T, np.log(ss.p)); logpt = np.full(T, np.log(ss.ptil_r))738 prices = E.build_price_path(P, paths, sh, T, logp, logpt, core)739 for t, pr in enumerate(prices):740 pr.rho = ss.rho * (1 + g) ** min(t, T)741 pols_by_t = E.solve_households_path(P, prices, core, T)742 recs = E.simulate_path(P, pols_by_t, prices, core.panel, T, 1.0)743 rows = []744 for t, rec in enumerate(recs):745 own = rec.owner1(); young = rec.age <= 3746 rows.append(dict(t=t, rho=prices[t].rho, ownership=float(own.mean()), ownership_young=float(own[young].mean()),747 buy_rate_young=float((rec.choice[young] == 2).mean()), nw_mean=float(rec.net_worth1(P).mean()),748 rental_demand=rec.rental_demand(P), housing_demand=rec.housing_demand(P)))749 out = DATA / "sensitivity"750 out.mkdir(parents=True, exist_ok=True)751 df = pd.DataFrame(rows)752 df.to_csv(out / f"{name}_path.csv", index=False)753 save_json(out / f"{name}.json", dict(name=name, group="sensitivity", scenario=scen, equilibrium="PE (prices fixed)",754 g_rent_period=g, path=df.to_dict(orient="list"), time_sec=time.time() - t0))755 return dict(name=name, group="sensitivity", ok=True, converged=True, equilibrium="PE (prices fixed)", sec=time.time() - t0)756757758def run_multiplicity(cfg_base: dict, scales: dict, base: dict, n_threads: int = 0) -> dict:759 """Experiment 3/7: scan the owner-market excess demand in p under near-fixed supply (and the760 baseline) with the rental market cleared, lottery fixed at the equilibrium one; count sign761 changes = number of candidate equilibria; also report d(ownership)/dp along the curve."""762 out = DATA / "multiplicity"763 out.mkdir(parents=True, exist_ok=True)764 results = {}765 for tag, ov in [("baseline", {}), ("eps_005", {"supply.eps": 0.05, "supply.eps_r": 0.05}),766 ("eps_005_theta005", {"supply.eps": 0.05, "supply.eps_r": 0.05, "mortgage.theta": 0.05})]:767 cfg = apply_overrides(cfg_base, ov)768 P = M.build_params(cfg)769 P.H_bar, P.Hr_bar = scales["H_bar"], scales["Hr_bar"]770 log = E.RunLog(f"multiplicity_{tag}", echo=False)771 lot = (np.array(base["lottery"]), np.full(P.K_estate, 1.0 / P.K_estate))772 p_values = np.round(np.arange(0.80, 1.251, 0.025), 3)773 scan = E.price_scan(P, p_values, lot, log)774 df = pd.DataFrame(scan)775 df.to_csv(out / f"price_scan_{tag}.csv", index=False)776 sg = np.sign(df.excess_owner.values)777 n_cross = int(np.sum(sg[1:] != sg[:-1]))778 slope = np.gradient(df.ownership.values, df.p.values)779 results[tag] = dict(n_sign_changes=n_cross, p_grid=p_values.tolist(), excess_owner=df.excess_owner.tolist(),780 ownership=df.ownership.tolist(), d_ownership_dp=slope.tolist(),781 monotone_decreasing=bool(np.all(np.diff(df.excess_owner.values) < 0)))782 log("scan_done", n_sign_changes=n_cross)783 save_json(out / "multiplicity_summary.json", results)784 return results785786787# --------------------------------------------------------------------------------------788# main789# --------------------------------------------------------------------------------------790def main(argv=None):791 ap = argparse.ArgumentParser()792 ap.add_argument("--groups", nargs="*", default=None, help="scenario groups to run (default: all)")793 ap.add_argument("--workers", type=int, default=max(1, min(9, (os.cpu_count() or 4) // 2)))794 ap.add_argument("--transition-workers", type=int, default=2, help="transitions solved concurrently")795 ap.add_argument("--quick", action="store_true", help="smaller panel / fewer scenarios (smoke test)")796 ap.add_argument("--skip-baseline", action="store_true")797 ap.add_argument("--rerun-transition", nargs="*", default=None, help="re-run only these transition scenarios")798 ap.add_argument("--warm-from", default=None, help="converged transition used as warm start for --rerun-transition")799 args = ap.parse_args(argv)800 T0 = time.time()801 DATA.mkdir(parents=True, exist_ok=True); LOGS.mkdir(parents=True, exist_ok=True)802 main_log = E.RunLog("simulate")803 cfg = M.load_config()804 if args.quick:805 cfg["simulation"]["n_dynasties"] = 24000806 cfg["transition"]["T"] = 8807 cfg["transition"]["max_iter"] = 3808 scenarios = M.load_scenarios()["groups"]809 groups = args.groups or list(scenarios.keys())810 diagnostics = []811 if args.rerun_transition:812 with open(DATA / "baseline" / "supply_scales.json") as fh:813 scales = json.load(fh)814 base = dict(p=scales["p"], ptil_r=scales["ptil_r"], lottery=scales["lottery"])815 all_scens = {sc["name"]: sc for g in scenarios.values() for sc in g}816 out = []817 for nm in args.rerun_transition:818 r = run_transition_scenario(cfg, all_scens[nm], scales, base, n_workers=min(16, (os.cpu_count() or 4) - 2),819 warm_from=args.warm_from)820 out.append(r)821 main_log("transition_rerun", **{k: v for k, v in r.items() if k != "trace"})822 # refresh diagnostics.json entries823 dpath = DATA / "diagnostics.json"824 if dpath.exists():825 with open(dpath) as fh:826 dj = json.load(fh)827 names = {r["name"] for r in out}828 dj["runs"] = [r for r in dj["runs"] if r.get("name") not in names] + out829 save_json(dpath, dj)830 pd.DataFrame([{k: v for k, v in r.items() if k not in ("bellman", "residuals", "trace")} for r in dj["runs"]]).to_csv(831 DATA / "diagnostics.csv", index=False)832 return833 # ---- baseline834 if not args.skip_baseline:835 P, res, summ = run_baseline(cfg, quick=args.quick)836 diagnostics.append(dict(name="baseline", group="baseline", ok=True, converged=summ["converged"], p=summ["p"],837 sec=summ["time_sec"], equilibrium=summ["equilibrium"], residuals=summ["residuals"], bellman=summ["bellman"]))838 main_log("baseline_done", sec=time.time() - T0, p=summ["p"], ownership=summ["moments"]["ownership_rate_all"])839 with open(DATA / "baseline" / "supply_scales.json") as fh:840 scales = json.load(fh)841 base = dict(p=scales["p"], ptil_r=scales["ptil_r"], lottery=scales["lottery"])842 # ---- stationary scenarios in parallel843 ss_jobs = []844 for g in groups:845 for scen in scenarios.get(g, []):846 if scen["kind"] in ("ss", "ss_pe") and scen["name"] != "baseline":847 if args.quick and g == "robustness" and scen["name"] not in ("grid_coarse", "sample_small"):848 continue849 ss_jobs.append((cfg, scen, g, scales, base["lottery"], base, max(1, (os.cpu_count() or 4) // args.workers)))850 if ss_jobs:851 main_log("ss_scenarios_start", n=len(ss_jobs), workers=args.workers)852 ctx = get_context("spawn")853 with ProcessPoolExecutor(max_workers=args.workers, mp_context=ctx) as ex:854 for r in ex.map(run_ss_scenario, ss_jobs):855 diagnostics.append(r)856 main_log("ss_scenario", **{k: v for k, v in r.items() if k not in ("bellman", "residuals")})857 # ---- multiplicity scan858 if any(g in groups for g in ("supply", "theta_sweep")):859 main_log("multiplicity_start")860 run_multiplicity(cfg, scales, base)861 # ---- transitions (a few concurrently; parallel Jacobian inside each)862 tr_jobs = []863 for g in groups:864 for scen in scenarios.get(g, []):865 if scen["kind"] in ("mit", "static"):866 if args.quick and scen["name"] not in ("frm_level", "term5_persistent", "frm_level_static"):867 continue868 tr_jobs.append((cfg, scen, scales, base, max(2, ((os.cpu_count() or 4) - 2) // max(1, args.transition_workers))))869 elif scen["kind"] == "rent_growth_pe":870 r = run_rent_growth_pe(cfg, scen, scales, args.workers)871 diagnostics.append(r)872 if tr_jobs:873 main_log("transitions_start", n=len(tr_jobs), concurrent=args.transition_workers)874 ctx = get_context("spawn")875 with ProcessPoolExecutor(max_workers=args.transition_workers, mp_context=ctx) as ex:876 for r in ex.map(_transition_job, tr_jobs):877 diagnostics.append(r)878 main_log("transition_done", **{k: v for k, v in r.items() if k != "trace"})879 # ---- welfare (experiment 9) across stationary scenarios880 with open(DATA / "baseline" / "ss_result.pkl", "rb") as fh:881 d = pickle.load(fh)882 main_log("welfare_start")883 run_welfare(cfg, {g: scenarios[g] for g in groups if g in scenarios}, d["P"], d["res"])884 # ---- diagnostics file885 save_json(DATA / "diagnostics.json", dict(runs=diagnostics, total_sec=time.time() - T0, quick=args.quick,886 cpu=os.cpu_count(), workers=args.workers))887 pd.DataFrame([{k: v for k, v in r.items() if k not in ("bellman", "residuals", "trace")} for r in diagnostics]).to_csv(888 DATA / "diagnostics.csv", index=False)889 main_log("all_done", sec=time.time() - T0)890891892if __name__ == "__main__":893 main()894