TeX 53.2%
Python 46.6%
1"""WP18 — Part A benchmark model (three-stage lifecycle, two-period generations).23Implements the closed-form hand-to-mouth (HTM) solution of `research/model_specification.md`4(Part A) and the objects used by the propositions of `paper/appendix_proofs.tex`:56* mortgage schedule (FRM / ARM, level or persistent surprise), real payments D2, D3;7* ownership premium Delta, Lambda = exp(-Delta), entry threshold b_hat (Proposition 1);8* HTM value functions and a general numerical household problem used to *verify*9 that the HTM corner is optimal (Assumption A4) rather than to assume it;10* dynasty transition probabilities, persistence lambda, stationary owner share (Prop. 5);11* stationary general equilibrium with segmented owner / rental markets (Prop. 6);12* leverage statistics (Prop. 4);13* surprise-inflation accounting, FRM vs ARM, level vs persistent (Prop. 8);14* impact general-equilibrium response to a surprise (Prop. 9) and a myopic-expectations15 transition simulation (diagnostic only — not covered by a proposition).1617Units: one model period = one stage of life; all flows real (consumption units); the house18price p and rent rho are per unit of housing; nominal principal B_nom = P * (1-theta) * p.19Parameters are illustrative (no calibration claim).20"""21from __future__ import annotations2223from dataclasses import dataclass, replace24from typing import Literal2526import numpy as np27from scipy import optimize, stats2829Contract = Literal["FRM", "ARM"]30Experiment = Literal["level", "persistent"]313233# --------------------------------------------------------------------------------------34# Parameters35# --------------------------------------------------------------------------------------36@dataclass(frozen=True)37class Params:38 # preferences39 beta: float = 0.80 # discount factor per stage40 nu: float = 0.30 # warm-glow weight: transfer share nu/(1+nu) of end-of-life wealth41 chi: float = 0.02 # ownership utility premium (utils, stage 1 + 2 combined)42 # incomes (real, per stage); y3 ~ y_lo + (y_hi - y_lo) * Beta(a_y, b_y)43 y1: float = 1.044 y2: float = 3.045 y3_lo: float = 2.046 y3_hi: float = 7.047 y3_a: float = 2.548 y3_b: float = 2.549 # housing finance and costs50 theta: float = 0.30 # minimum down-payment share (max LTV = 0.70)51 kappa_b: float = 0.03 # buying cost (share of price)52 kappa_s: float = 0.05 # selling cost (share of price)53 m: float = 0.05 # maintenance per stage (share of price)54 amort: float = 0.50 # share of principal amortised at stage 255 # rates and inflation56 r: float = 0.10 # real rate per stage57 pibar: float = 0.10 # steady-state inflation per stage58 # supply curves (per cohort): S(p) = Hbar p^eps ; S_r(rho) = Hbar_r rho^eps_r59 Hbar: float = 0.6060 eps: float = 1.061 Hbar_r: float = 0.8062 eps_r: float = 1.063 # indexation of nominal wages / rents to a surprise (1 = full)64 eta_w: float = 1.065 eta_rho: float = 1.066 # quadrature nodes for expectations over y367 n_quad: int = 646869 @property70 def R(self) -> float:71 return 1.0 + self.r7273 @property74 def d(self) -> float:75 """Stage-1 cash outlay per unit of price: down payment + buying cost + maintenance."""76 return self.theta + self.kappa_b + self.m777879# --------------------------------------------------------------------------------------80# Income distribution G (stage-3 income) — bounded support, smooth density81# --------------------------------------------------------------------------------------82def y3_dist(par: Params):83 return stats.beta(par.y3_a, par.y3_b, loc=par.y3_lo, scale=par.y3_hi - par.y3_lo)848586def y3_quad(par: Params):87 """Gauss–Legendre nodes/weights for E[f(y3)] with y3 ~ G on [y_lo, y_hi]."""88 x, w = np.polynomial.legendre.leggauss(par.n_quad)89 lo, hi = par.y3_lo, par.y3_hi90 y = lo + (hi - lo) * (x + 1) / 291 pdf = y3_dist(par).pdf(y)92 wts = w * (hi - lo) / 2 * pdf93 wts = wts / wts.sum() # renormalise (quadrature error of the density integral)94 return y, wts959697def E_y3(par: Params, f) -> float:98 y, w = y3_quad(par)99 return float(np.sum(w * f(y)))100101102def G(par: Params, y) -> np.ndarray | float:103 return y3_dist(par).cdf(y)104105106def g(par: Params, y) -> np.ndarray | float:107 return y3_dist(par).pdf(y)108109110# --------------------------------------------------------------------------------------111# Mortgage schedule112# --------------------------------------------------------------------------------------113def nominal_rate(par: Params, pi_e: float, r: float | None = None) -> float:114 """Fisher: 1 + i = R (1 + pi_e)."""115 R = par.R if r is None else 1.0 + r116 return R * (1.0 + pi_e) - 1.0117118119def real_payments(par: Params, p: float, *, pi1: float, pi2: float, pi_e1: float, pi_e2: float,120 contract: Contract = "FRM", r2: float | None = None,121 principal_share: float | None = None) -> tuple[float, float]:122 """Real stage-2 and stage-3 payments of a mortgage originated at price p.123124 pi1, pi2 : realised inflation over (t, t+1) and (t+1, t+2)125 pi_e1 : inflation expected at origination for (t, t+1) -> fixes i_1 for all contracts126 pi_e2 : inflation expected at t+1 for (t+1, t+2) -> ARM reset rate127 r2 : real rate used at the ARM reset (None = par.r)128 """129 B = ((1.0 - par.theta) if principal_share is None else principal_share) * p130 i1 = nominal_rate(par, pi_e1)131 if contract == "FRM":132 i2 = i1133 elif contract == "ARM":134 i2 = nominal_rate(par, pi_e2, r2)135 else:136 raise ValueError(contract)137 D2 = (i1 + par.amort) * B / (1.0 + pi1)138 D3 = (1.0 + i2) * (1.0 - par.amort) * B / ((1.0 + pi1) * (1.0 + pi2))139 return D2, D3140141142def stationary_schedule(par: Params) -> tuple[float, float, float]:143 """(delta2, delta3, e3): payments per unit of price and net equity share at sale."""144 D2, D3 = real_payments(par, 1.0, pi1=par.pibar, pi2=par.pibar, pi_e1=par.pibar, pi_e2=par.pibar)145 e3 = (1.0 - par.kappa_s) - D3146 return D2, D3, e3147148149# --------------------------------------------------------------------------------------150# HTM closed forms (Part A, §A.5–A.6)151# --------------------------------------------------------------------------------------152def ownership_premium(par: Params, p: float, rho: float, *, p_sale: float | None = None,153 rho2: float | None = None, D2: float | None = None, D3: float | None = None) -> float:154 """Delta = chi + beta E[ln c2^O - ln c2^R] + beta^2 (1+nu) E[ln W3^O - ln W3^R] under HTM."""155 p_sale = p if p_sale is None else p_sale156 rho2 = rho if rho2 is None else rho2157 if D2 is None or D3 is None:158 d2, d3, _ = stationary_schedule(par)159 D2 = d2 * p if D2 is None else D2160 D3 = d3 * p if D3 is None else D3161 c2O = par.y2 - D2 - par.m * p162 c2R = par.y2 - rho2163 if c2O <= 0 or c2R <= 0:164 raise ValueError("Assumption A3 violated: non-positive stage-2 consumption")165 eq3 = (1.0 - par.kappa_s) * p_sale - D3166 if par.y3_lo + eq3 <= 0:167 raise ValueError("Assumption A3 violated: owner end-of-life wealth can be negative")168 term3 = E_y3(par, lambda y: np.log(y + eq3) - np.log(y))169 return par.chi + par.beta * (np.log(c2O) - np.log(c2R)) + par.beta**2 * (1 + par.nu) * term3170171172def threshold(par: Params, p: float, rho: float, **kw) -> float:173 """Ownership-entry threshold b_hat (Proposition 1). Returns -inf-like value if nobody is excluded."""174 Delta = ownership_premium(par, p, rho, **kw)175 d = par.d176 b_min = d * p - par.y1177 if Delta <= 0:178 return np.inf # nobody owns (A5 fails): owning is dominated at every wealth level179 Lam = np.exp(-Delta)180 if d * p <= rho:181 return b_min182 return (d * p - Lam * rho) / (1.0 - Lam) - par.y1183184185def value_htm(par: Params, b: float, p: float, rho: float, tenure: Literal["O", "R"], **kw) -> float:186 """Lifetime value of a hand-to-mouth household (n1 = n2 = 0, max mortgage if owner)."""187 K = par.beta**2 * (par.nu * np.log(par.nu) - (1 + par.nu) * np.log(1 + par.nu))188 d2, d3, _ = stationary_schedule(par)189 D2 = kw.get("D2", d2 * p)190 D3 = kw.get("D3", d3 * p)191 p_sale = kw.get("p_sale", p)192 rho2 = kw.get("rho2", rho)193 if tenure == "O":194 c1 = b + par.y1 - par.d * p195 if c1 <= 0:196 return -np.inf197 c2 = par.y2 - D2 - par.m * p198 eq3 = (1 - par.kappa_s) * p_sale - D3199 return (np.log(c1) + par.chi + par.beta * np.log(c2)200 + par.beta**2 * (1 + par.nu) * E_y3(par, lambda y: np.log(y + eq3)) + K)201 c1 = b + par.y1 - rho202 if c1 <= 0:203 return -np.inf204 c2 = par.y2 - rho2205 return (np.log(c1) + par.beta * np.log(c2)206 + par.beta**2 * (1 + par.nu) * E_y3(par, lambda y: np.log(y)) + K)207208209def htm_slack(par: Params, b: float, p: float, rho: float, tenure: Literal["O", "R"]) -> tuple[float, float]:210 """Euler-inequality slacks (Assumption A4): positive values mean the household would like211 to borrow more at stage 1 and at stage 2 (so n1 = n2 = 0 is the optimum)."""212 d2, d3, e3 = stationary_schedule(par)213 R = par.R214 if tenure == "O":215 c1 = b + par.y1 - par.d * p216 c2 = par.y2 - d2 * p - par.m * p217 EW3inv = E_y3(par, lambda y: 1.0 / (y + e3 * p))218 else:219 c1 = b + par.y1 - rho220 c2 = par.y2 - rho221 EW3inv = E_y3(par, lambda y: 1.0 / y)222 s1 = 1.0 / c1 - par.beta * R * (1.0 / c2) # stage 1 vs 2 (c2 deterministic)223 s2 = 1.0 / c2 - par.beta * R * (1 + par.nu) * EW3inv # stage 2 vs 3224 return s1, s2225226227# --------------------------------------------------------------------------------------228# General (non-HTM) household problem — numerical, used to verify the corner solution229# --------------------------------------------------------------------------------------230def value_general(par: Params, b: float, p: float, rho: float, tenure: Literal["O", "R"],231 return_policy: bool = False):232 """Owner chooses mortgage share phi in [0, 1-theta], savings n1 >= 0, n2 >= 0.233 Renter chooses n1 >= 0, n2 >= 0. Stage-3 split by Lemma 1. Returns the max value."""234 y, w = y3_quad(par)235 K = par.beta**2 * (par.nu * np.log(par.nu) - (1 + par.nu) * np.log(1 + par.nu))236 R = par.R237238 def obj_owner(x):239 phi, n1, n2 = x240 D2, D3 = real_payments(par, p, pi1=par.pibar, pi2=par.pibar, pi_e1=par.pibar, pi_e2=par.pibar,241 principal_share=phi)242 c1 = b + par.y1 - (par.theta + par.kappa_b + par.m) * p - ((1 - par.theta) - phi) * p - n1243 c2 = par.y2 - D2 - par.m * p + n1 * R - n2244 W3 = y + (1 - par.kappa_s) * p - D3 + n2 * R245 if c1 <= 1e-9 or c2 <= 1e-9 or np.any(W3 <= 1e-9):246 return 1e6247 return -(np.log(c1) + par.chi + par.beta * np.log(c2)248 + par.beta**2 * (1 + par.nu) * np.sum(w * np.log(W3)) + K)249250 def obj_renter(x):251 n1, n2 = x252 c1 = b + par.y1 - rho - n1253 c2 = par.y2 - rho + n1 * R - n2254 W3 = y + n2 * R255 if c1 <= 1e-9 or c2 <= 1e-9 or np.any(W3 <= 1e-9):256 return 1e6257 return -(np.log(c1) + par.beta * np.log(c2)258 + par.beta**2 * (1 + par.nu) * np.sum(w * np.log(W3)) + K)259260 best = None261 if tenure == "O":262 bounds = [(0.0, 1 - par.theta), (0.0, 5.0), (0.0, 5.0)]263 starts = [(1 - par.theta, 0.0, 0.0), (0.5 * (1 - par.theta), 0.1, 0.1), (0.0, 0.0, 0.5)]264 f = obj_owner265 else:266 bounds = [(0.0, 5.0), (0.0, 5.0)]267 starts = [(0.0, 0.0), (0.2, 0.2), (0.0, 1.0)]268 f = obj_renter269 for x0 in starts:270 res = optimize.minimize(f, np.array(x0), bounds=bounds, method="L-BFGS-B")271 if res.fun < 1e5 and (best is None or res.fun < best.fun):272 best = res273 if best is None:274 return (-np.inf, None) if return_policy else -np.inf275 return (-best.fun, best.x) if return_policy else -best.fun276277278def threshold_general(par: Params, p: float, rho: float) -> float:279 """Threshold from the general (numerically optimised) values; verifies Proposition 1."""280 lo = par.d * p - par.y1 + 1e-6281 hi = lo + 10.0282 f = lambda b: value_general(par, b, p, rho, "O") - value_general(par, b, p, rho, "R")283 if f(hi) < 0:284 return np.inf285 if f(lo) > 0:286 return lo287 return optimize.brentq(f, lo, hi, xtol=1e-8)288289290# --------------------------------------------------------------------------------------291# Dynasties and persistence (Proposition 5)292# --------------------------------------------------------------------------------------293def transitions(par: Params, p: float, rho: float, b_hat: float | None = None, *,294 eq3: float | None = None) -> dict:295 """q_RO, q_OO, lambda, stationary owner share mu* for given prices.296 eq3 = owner's net housing equity at sale (default e3 p)."""297 if b_hat is None:298 b_hat = threshold(par, p, rho)299 _, _, e3 = stationary_schedule(par)300 eq3 = e3 * p if eq3 is None else eq3301 yR = (1 + par.nu) / par.nu * b_hat302 yO = yR - eq3303 qRO = 1.0 - float(G(par, yR)) if np.isfinite(b_hat) else 0.0304 qOO = 1.0 - float(G(par, yO)) if np.isfinite(b_hat) else 0.0305 lam = qOO - qRO306 mu_star = qRO / (1.0 - lam) if (1.0 - lam) > 0 else 1.0307 half_life = np.log(2) / (-np.log(lam)) if 0 < lam < 1 else (np.inf if lam >= 1 else 0.0)308 return dict(b_hat=b_hat, yR=yR, yO=yO, qRO=qRO, qOO=qOO, lam=lam, mu_star=mu_star,309 half_life=half_life, transfer_gap=par.nu / (1 + par.nu) * eq3)310311312# --------------------------------------------------------------------------------------313# Supply and stationary equilibrium (Proposition 6)314# --------------------------------------------------------------------------------------315def S(par: Params, p: float) -> float:316 return par.Hbar * p ** par.eps317318319def S_inv(par: Params, mu: float) -> float:320 return (mu / par.Hbar) ** (1.0 / par.eps)321322323def Sr(par: Params, rho: float) -> float:324 return par.Hbar_r * rho ** par.eps_r325326327def Sr_inv(par: Params, share: float) -> float:328 return (share / par.Hbar_r) ** (1.0 / par.eps_r)329330331def p_upper(par: Params) -> float:332 """Largest price compatible with Assumption A3 (positive stage-2 owner consumption)."""333 d2, _, _ = stationary_schedule(par)334 return 0.999 * par.y2 / (d2 + par.m)335336337def excess_map(par: Params, mu: float) -> float:338 """f(mu) = mu*(S^-1(mu), S_r^-1(1-mu)) - mu."""339 p = S_inv(par, mu)340 rho = Sr_inv(par, 1.0 - mu)341 return transitions(par, p, rho)["mu_star"] - mu342343344def stationary_equilibrium(par: Params, mu_lo: float = 1e-4, mu_hi: float = 1 - 1e-4) -> dict:345 """Solve f(mu) = 0 on [mu_lo, mu_hi] (IVT); returns prices, shares, thresholds, diagnostics.346 With a fixed owner stock (eps = 0) the owner share is mu = Hbar and the price solves347 mu*(p, rho) = Hbar (Corollary to Proposition 6)."""348 if par.eps == 0:349 mu = par.Hbar350 rho = Sr_inv(par, 1.0 - mu)351 h = lambda p: transitions(par, p, rho)["mu_star"] - mu352 p = optimize.brentq(h, 1e-3, p_upper(par), xtol=1e-12)353 else:354 mu_hi = min(mu_hi, S(par, p_upper(par)) * (1 - 1e-9)) # keep A3 along the bracket355 f_lo, f_hi = excess_map(par, mu_lo), excess_map(par, mu_hi)356 if f_lo <= 0:357 mu = mu_lo358 elif f_hi >= 0:359 mu = mu_hi360 else:361 mu = optimize.brentq(lambda x: excess_map(par, x), mu_lo, mu_hi, xtol=1e-12)362 p = S_inv(par, mu)363 rho = Sr_inv(par, 1.0 - mu)364 tr = transitions(par, p, rho)365 out = dict(mu=mu, p=p, rho=rho, **tr)366 out["Delta"] = ownership_premium(par, p, rho)367 out["Lambda"] = np.exp(-out["Delta"])368 out["owner_clearing"] = S(par, p) - mu369 out["rental_clearing"] = Sr(par, rho) - (1 - mu)370 out["A5_cash"] = par.d * p - rho # > 0 : owning costs more cash up front371 return out372373374def uniqueness_condition(par: Params, p: float, rho: float, h: float = 1e-5) -> float:375 """Condition (C): (1+nu)/nu * d b_hat/dp - e3 >= 0 (b_hat channel dominates ladder channel)."""376 _, _, e3 = stationary_schedule(par)377 db = (threshold(par, p + h, rho) - threshold(par, p - h, rho)) / (2 * h)378 return (1 + par.nu) / par.nu * db - e3379380381def no_segmentation_equilibrium(par: Params, phi: float = 0.0) -> dict:382 """Proposition 7 benchmark: a single stock with costless conversion and unit demand.383 Total demand per cohort is 1, so p solves S(p) = 1 and does not depend on theta, G, transfers384 or surprises; the rent follows landlord arbitrage rho = p (1 + m + phi - 1/R)."""385 p = S_inv(par, 1.0)386 rho = p * (1.0 + par.m + phi - 1.0 / par.R)387 bh = threshold(par, p, rho)388 tr = transitions(par, p, rho, bh)389 return dict(p=p, rho=rho, b_hat=bh, **{k: v for k, v in tr.items() if k != "b_hat"})390391392# --------------------------------------------------------------------------------------393# Leverage amplification with costs (Proposition 4)394# --------------------------------------------------------------------------------------395def leverage_stats(par: Params, p: float, rho: float) -> dict:396 d2, d3, e3 = stationary_schedule(par)397 R = par.R398 multiplier = (1 - par.kappa_s) / par.d # dW3 per unit of initial cash, per unit price growth399 rhs = R**2 * (1 + par.kappa_b + par.m) + par.m * R - (rho / p) * (R + R**2)400 g_star = rhs / (1 - par.kappa_s) - 1.0 # break-even real price growth401 x_star = 1.0 - d3 / (1 - par.kappa_s) # underwater threshold (price drop)402 return dict(multiplier=multiplier, g_star=g_star, x_star=x_star, e3=e3,403 wealth_gap=e3 * p, transfer_gap=par.nu / (1 + par.nu) * e3 * p)404405406def owner_vs_renter_investor(par: Params, p: float, rho: float, growth: float) -> float:407 """Terminal wealth of an owner minus that of a renter who invests the cash-flow differences408 at R (same consumption path). Sign flips at g_star (Proposition 4)."""409 d2, d3, _ = stationary_schedule(par)410 R = par.R411 owner = (1 - par.kappa_s) * p * (1 + growth) - d3 * p412 renter = (par.d * p - rho) * R**2 + (d2 * p + par.m * p - rho) * R413 return owner - renter414415416# --------------------------------------------------------------------------------------417# Surprise inflation accounting (Proposition 8)418# --------------------------------------------------------------------------------------419def surprise_accounting(par: Params, p: float, rho: float, pi_u: float, *, contract: Contract = "FRM",420 experiment: Experiment = "level", phi_r: float = 0.0,421 n_savings_age2: float = 0.0, n_savings_age3: float = 0.0) -> dict:422 """Real gains (+) / losses (-) relative to the no-surprise path at fixed real prices,423 for the cohorts aged 2 and 3 when the surprise hits, the lender, landlords and wages.424425 n_savings_age2 / age3: nominal deposits held by (non-HTM) households of those ages, to show426 that nominal savers lose. Returns a dict with a zero-sum check on nominal claims."""427 pb = par.pibar428 pi1 = (1 + pb) * (1 + pi_u) - 1429 if experiment == "level":430 pi2 = pb431 pi_e2 = pb432 elif experiment == "persistent":433 pi2 = pi1434 pi_e2 = pi1435 else:436 raise ValueError(experiment)437 r2 = par.r + phi_r * pi_u438 # baseline (no surprise)439 D2_0, D3_0 = real_payments(par, p, pi1=pb, pi2=pb, pi_e1=pb, pi_e2=pb, contract=contract)440 # cohort aged 2 at the shock (bought at t)441 D2_s, D3_s = real_payments(par, p, pi1=pi1, pi2=pi2, pi_e1=pb, pi_e2=pi_e2, contract=contract, r2=r2)442 gain_age2_D2 = D2_0 - D2_s443 gain_age2_D3 = D3_0 - D3_s444 # cohort aged 3 at the shock (bought at t-1): stage-3 nominal payment fixed before the shock445 gain_age3_D3 = D3_0 - D3_0 / (1 + pi_u)446 # nominal savers447 R = par.R448 loss_savers = (n_savings_age2 + n_savings_age3) * R * (1 - 1 / (1 + pi_u))449 # wages (imperfect indexation) — cohort aged 2 and 3 real income at the shock date450 wage_factor = (1 + pi_u) ** (-(1 - par.eta_w))451 wage_loss_age2 = par.y2 * (1 - wage_factor)452 # rents (rigidity): sitting tenants aged 2 pay less; landlords absorb453 rent_factor = (1 + pi_u) ** (-(1 - par.eta_rho))454 renter_gain_age2 = rho * (1 - rent_factor)455 # lender: mirror image of all nominal revaluations on the shock date and (for the age-2 cohort)456 # on the continuation, both contracts457 lender_loss = (gain_age2_D2 + gain_age2_D3 + gain_age3_D3) - loss_savers458 zero_sum = (gain_age2_D2 + gain_age2_D3 + gain_age3_D3 - loss_savers) - lender_loss459 net_owner_age2 = gain_age2_D2 + gain_age2_D3 - wage_loss_age2 - n_savings_age2 * R * (1 - 1 / (1 + pi_u))460 return dict(D2_base=D2_0, D3_base=D3_0, D2_shock=D2_s, D3_shock=D3_s,461 gain_age2_D2=gain_age2_D2, gain_age2_D3=gain_age2_D3, gain_age3_D3=gain_age3_D3,462 loss_savers=loss_savers, wage_loss_age2=wage_loss_age2,463 renter_gain_age2=renter_gain_age2, landlord_loss=renter_gain_age2,464 lender_loss=lender_loss, zero_sum_residual=zero_sum, net_owner_age2=net_owner_age2)465466467# --------------------------------------------------------------------------------------468# Impact general-equilibrium response (Proposition 9) under static expectations (A6)469# --------------------------------------------------------------------------------------470def impact_response(par: Params, eq: dict, pi_u: float, *, contract: Contract = "FRM") -> dict:471 """Date t+1 equilibrium after a level surprise pi_u, with (i) stage-2 owners not transacting,472 (ii) stage-3 owner-parents' D3 eroded by 1/(1+pi_u), (iii) the stage-1 cohort's threshold473 evaluated at current (p, rho) and stationary future prices (A6).474475 Owner market: S(p) = (mu* + mu_{t+1})/2 ; rental: S_r(rho) = 1 - (mu* + mu_{t+1})/2.476 Solved by a root search on mu_{t+1}."""477 mu_s, p_s, rho_s = eq["mu"], eq["p"], eq["rho"]478 d2, d3, _ = stationary_schedule(par)479 D3_parent = d3 * p_s / (1 + pi_u) # stage-3 parents bought at p_s; payment eroded480481 def mu_next_given(mu1, p_fixed=None):482 avg = (mu_s + mu1) / 2483 p = S_inv(par, avg) if p_fixed is None else p_fixed484 rho = Sr_inv(par, 1 - avg)485 # child's threshold: buys at p now, pays D2, D3 on (1-theta) p, sells at p_s (static exp.), rents at rho now, rho_s later486 bh = threshold(par, p, rho, p_sale=p_s, rho2=rho_s, D2=d2 * p, D3=d3 * p)487 eq3_parent = (1 - par.kappa_s) * p - D3_parent # parent's equity at sale (sells at p)488 tr = transitions(par, p, rho, bh, eq3=eq3_parent)489 mu_new = mu_s * tr["qOO"] + (1 - mu_s) * tr["qRO"]490 return mu_new, p, rho, bh, tr491492 if par.eps == 0:493 # fixed owner stock: mu_{t+1} = mu* and the price clears the owner market494 mu1 = mu_s495 hfun = lambda p: mu_next_given(mu1, p_fixed=p)[0] - mu1496 p_sol = optimize.brentq(hfun, 1e-3, p_upper(par), xtol=1e-12)497 mu_new, p, rho, bh, tr = mu_next_given(mu1, p_fixed=p_sol)498 else:499 f = lambda mu1: mu_next_given(mu1)[0] - mu1500 lo, hi = 1e-6, min(1 - 1e-6, 2 * S(par, p_upper(par)) * (1 - 1e-9) - mu_s)501 if f(lo) <= 0:502 mu1 = lo503 elif f(hi) >= 0:504 mu1 = hi505 else:506 mu1 = optimize.brentq(f, lo, hi, xtol=1e-12)507 mu_new, p, rho, bh, tr = mu_next_given(mu1)508 return dict(mu_next=mu1, p=p, rho=rho, b_hat=bh, qRO=tr["qRO"], qOO=tr["qOO"], lam=tr["lam"],509 dp=p - p_s, drho=rho - rho_s, dmu=mu1 - mu_s, db_hat=bh - eq["b_hat"],510 rent_price=rho / p, rent_price_base=rho_s / p_s,511 seller_gain=d3 * p_s * (1 - 1 / (1 + pi_u)) + (1 - par.kappa_s) * (p - p_s),512 new_owner_cash_change=par.nu / (1 + par.nu) * (d3 * p_s * (1 - 1 / (1 + pi_u))513 + (1 - par.kappa_s) * (p - p_s))514 - par.d * (p - p_s))515516517def simulate_myopic_transition(par: Params, eq: dict, pi_u: float, T: int = 12,518 experiment: Experiment = "level", contract: Contract = "FRM") -> dict:519 """Diagnostic simulation of the dynasty composition after a surprise at date 1 under the520 Part A structure with per-period supply and static expectations. Tracks, for each cohort,521 the tenure share and the purchase price, so that parents' equity at sale uses the right522 purchase price and payment erosion. Not a proposition: a check that the economy returns to523 the stationary state and that the impact signs match `impact_response`."""524 d2, d3, _ = stationary_schedule(par)525 pb = par.pibar526 # histories are indexed by cohort c = -1, 0, 1, ... at position c + 1527 mu_hist = [eq["mu"], eq["mu"]] # cohorts born at t=-1, 0528 p_hist = [eq["p"], eq["p"]]529 rho_hist = [eq["rho"], eq["rho"]]530 rows = []531 for t in range(1, T + 1):532 # parents of this cohort: born at t-2 (stage 3 now), position t-1; previous cohort: position t533 p_buy = p_hist[t - 1]534 mu_parents = mu_hist[t - 1]535 # erosion of the parents' D3 and of D2/D3 of the cohort born at t-1 (not needed for mu)536 if t == 1:537 pi1, pi2, pi_e2 = (1 + pb) * (1 + pi_u) - 1, pb, pb538 _, D3_par = real_payments(par, p_buy, pi1=pb, pi2=pi1, pi_e1=pb, pi_e2=pb, contract=contract)539 elif t == 2 and experiment == "persistent":540 pi1 = (1 + pb) * (1 + pi_u) - 1541 _, D3_par = real_payments(par, p_buy, pi1=pi1, pi2=pi1, pi_e1=pb, pi_e2=pi1, contract=contract)542 elif t == 2:543 pi1 = (1 + pb) * (1 + pi_u) - 1544 _, D3_par = real_payments(par, p_buy, pi1=pi1, pi2=pb, pi_e1=pb, pi_e2=pb, contract=contract)545 else:546 D3_par = d3 * p_buy547 mu_prev = mu_hist[t]548549 def mu_given(mu1):550 avg = (mu_prev + mu1) / 2551 p = S_inv(par, avg)552 rho = Sr_inv(par, 1 - avg)553 bh = threshold(par, p, rho, p_sale=eq["p"], rho2=eq["rho"], D2=d2 * p, D3=d3 * p)554 eq3_par = (1 - par.kappa_s) * p - D3_par555 tr = transitions(par, p, rho, bh, eq3=eq3_par)556 return mu_parents * tr["qOO"] + (1 - mu_parents) * tr["qRO"], p, rho, bh557558 f = lambda x: mu_given(x)[0] - x559 hi = min(1 - 1e-6, 2 * S(par, p_upper(par)) * (1 - 1e-9) - mu_prev)560 mu1 = optimize.brentq(f, 1e-6, hi, xtol=1e-12) if f(1e-6) > 0 > f(hi) else (1e-6 if f(1e-6) <= 0 else hi)561 _, p, rho, bh = mu_given(mu1)562 mu_hist.append(mu1); p_hist.append(p); rho_hist.append(rho)563 rows.append(dict(t=t, mu=mu1, p=p, rho=rho, b_hat=bh))564 return dict(rows=rows, mu=np.array(mu_hist[2:]), p=np.array(p_hist[2:]), rho=np.array(rho_hist[2:]))565566567# --------------------------------------------------------------------------------------568# Budget identities (used by the tests)569# --------------------------------------------------------------------------------------570def htm_allocation(par: Params, b: float, p: float, rho: float, tenure: Literal["O", "R"], y3: float,571 *, pi_u1: float = 0.0, pi_u2: float = 0.0, contract: Contract = "FRM",572 experiment: Experiment = "level", p_sale: float | None = None) -> dict:573 """Full real flows of an HTM household for a given y3 draw; used to check budget identities574 state by state. Returns sources and uses at each stage, the nominal principal path and the575 real debt path (B_nom / P)."""576 pb = par.pibar577 P0 = 1.0578 pi1 = (1 + pb) * (1 + pi_u1) - 1579 if experiment == "persistent" and pi_u1 != 0:580 pi2, pi_e2 = pi1, pi1581 else:582 pi2 = (1 + pb) * (1 + pi_u2) - 1583 pi_e2 = pb584 P1, P2 = P0 * (1 + pi1), P0 * (1 + pi1) * (1 + pi2)585 p_sale = p if p_sale is None else p_sale586 out = dict(P=(P0, P1, P2))587 if tenure == "O":588 B_nom0 = P0 * (1 - par.theta) * p589 i1 = nominal_rate(par, pb)590 i2 = i1 if contract == "FRM" else nominal_rate(par, pi_e2)591 pay2_nom = (i1 + par.amort) * B_nom0592 B_nom1 = B_nom0 * (1 + i1) - pay2_nom # balance after stage-2 payment = (1-a) B_nom0593 pay3_nom = (1 + i2) * B_nom1594 D2, D3 = pay2_nom / P1, pay3_nom / P2595 c1 = b + par.y1 - par.d * p596 c2 = par.y2 - D2 - par.m * p597 W3 = y3 + (1 - par.kappa_s) * p_sale - D3598 out.update(B_nom=(B_nom0, B_nom1, 0.0), real_debt=(B_nom0 / P0, B_nom1 / P1, 0.0),599 D2=D2, D3=D3, c1=c1, c2=c2,600 sources=(b + par.y1 + B_nom0 / P0, par.y2, y3 + (1 - par.kappa_s) * p_sale),601 uses=(c1 + (1 + par.kappa_b + par.m) * p, c2 + D2 + par.m * p, None))602 # check D2, D3 against real_payments603 D2c, D3c = real_payments(par, p, pi1=pi1, pi2=pi2, pi_e1=pb, pi_e2=pi_e2, contract=contract)604 out["schedule_check"] = (D2 - D2c, D3 - D3c)605 else:606 c1 = b + par.y1 - rho607 c2 = par.y2 - rho608 W3 = y3609 out.update(c1=c1, c2=c2, sources=(b + par.y1, par.y2, y3), uses=(c1 + rho, c2 + rho, None))610 c3 = W3 / (1 + par.nu)611 bq = par.nu * W3 / (1 + par.nu)612 out.update(W3=W3, c3=c3, bequest=bq)613 out["uses"] = (out["uses"][0], out["uses"][1], c3 + bq + (D3 if tenure == "O" else 0.0))614 return out615616617if __name__ == "__main__": # quick report618 par = Params()619 eq = stationary_equilibrium(par)620 print("Stationary equilibrium:")621 for k, v in eq.items():622 print(f" {k:>18s} = {v:.6f}" if isinstance(v, float) else f" {k:>18s} = {v}")623 print("Uniqueness condition (C) at eq:", uniqueness_condition(par, eq["p"], eq["rho"]))624 print("Leverage:", leverage_stats(par, eq["p"], eq["rho"]))625 for contract in ("FRM", "ARM"):626 for exp_ in ("level", "persistent"):627 acc = surprise_accounting(par, eq["p"], eq["rho"], 0.10, contract=contract, experiment=exp_)628 print(contract, exp_, {k: round(v, 5) for k, v in acc.items() if k.startswith(("gain", "lender"))})629 imp = impact_response(par, eq, 0.10)630 print("Impact:", {k: round(v, 5) for k, v in imp.items()})631