TeX 53.2%
Python 46.6%
1"""WP18 — cross-sectional distribution by Monte-Carlo simulation of linked dynasties.23Why Monte Carlo and not a Young (2010) histogram: the dynasty link pays the parent's gift to4the child at the child's age 1 and the parent's estate at the child's age j_estate = 7; the two5transfers are correlated through the parent's state, so the invariant distribution lives on the6joint parent-child state, which a histogram on the individual state grid does not carry. A panel7of `n_dyn` dynasty slots, each holding one "young" member (ages 1-6) and one "old" member8(ages 7-12) at all times, tracks the link exactly. Random numbers are common across price9iterations (seeded by (seed, period)) so that excess demand is a deterministic function of10prices; sample-size sensitivity is reported in experiment 10.1112Slot phase s in {0..5}: young age = 1 + s, old age = 7 + s. When the old member dies (after13age J) the young member becomes the old member (age 7) and receives the estate; a new young14member (age 1) is born in the same period and receives the gift chosen by the new old member.15Within a period, old members decide first (their gift is the newborn's initial asset).16"""17from __future__ import annotations1819from dataclasses import dataclass, field2021import numpy as np22from numba import njit, prange2324from .household_solver import NEG, AgePolicy, _lin, _locate, _af_due, _renter_cs, _u25from .model import Params, Prices2627CHOICE_NAMES = {-1: "forced_sale", 0: "stay", 1: "refinance", 2: "buy", 3: "rent"}282930# --------------------------------------------------------------------------------------31# Stacked policies (ages along axis 0) for the numba kernel32# --------------------------------------------------------------------------------------33@dataclass34class PolicyStack:35 W_S: np.ndarray; aS: np.ndarray; gS: np.ndarray36 W_F: np.ndarray; BF: np.ndarray; aF: np.ndarray; gF: np.ndarray37 W_B: np.ndarray; hB: np.ndarray; BB: np.ndarray; aB: np.ndarray; gB: np.ndarray38 W_R: np.ndarray; aR: np.ndarray; gR: np.ndarray39 W_Rs: np.ndarray; aRs: np.ndarray; gRs: np.ndarray40 W_Bh: np.ndarray = None; BBh: np.ndarray = None; aBh: np.ndarray = None; gBh: np.ndarray = None414243def stack_policies(pols: list[AgePolicy], dtype=np.float64) -> PolicyStack:44 def st(name, dt=dtype):45 return np.ascontiguousarray(np.stack([getattr(p, name) for p in pols]), dtype=dt)46 return PolicyStack(W_S=st("W_S"), aS=st("aS"), gS=st("gS"),47 W_F=st("W_F"), BF=st("BF"), aF=st("aF"), gF=st("gF"),48 W_B=st("W_B"), hB=st("hB", np.int64), BB=st("BB"), aB=st("aB"), gB=st("gB"),49 W_R=st("W_R"), aR=st("aR"), gR=st("gR"),50 W_Rs=st("W_Rs"), aRs=st("aRs"), gRs=st("gRs"),51 W_Bh=st("W_Bh"), BBh=st("BBh"), aBh=st("aBh"), gBh=st("gBh"))525354# --------------------------------------------------------------------------------------55# One-period decision kernel56# --------------------------------------------------------------------------------------57@njit(cache=True, inline="always")58def _ix(x_grid, x):59 j, w = _locate(x_grid, x)60 return j, w616263@njit(cache=True, inline="always")64def _v1(arr, j, w):65 return _lin(arr[j], arr[j + 1], w)666768@njit(cache=True, inline="always")69def _p1(arr, j, w):70 """Policy interpolation (linear, no infeasibility sentinel)."""71 return (1.0 - w) * arr[j] + w * arr[j + 1]727374@njit(cache=True, parallel=True)75def _sim_step(age, z, h, b, n, a, y_base, y_noise,76 W_S, aS, gS, W_F, BF, aF, gF, W_B, hB, BB, aB, gB, W_R, aR, gR, W_Rs, aRs, gRs, W_Bh, BBh, aBh, gBh,77 a_grid, b_grid, x_grid, h_grid, ic_by_n, reset,78 p, p_next, rho, rho_sit, R_a, i_new, pi_next, rebate, wage_factor,79 kappa_b, kappa_s, kappa_f, m, chi, alpha, gamma, s_min, s_max, c_min, tau_b, N, j_gift, J,80 out_choice, out_c, out_s, out_h, out_b, out_n, out_a, out_g, out_pay, out_Bnew, out_estate,81 out_flag, out_u, out_x):82 M = age.shape[0]83 af_new = _af_due(i_new, N)84 for i in prange(M):85 j = age[i]86 jj = j - 187 iz = z[i]88 y = y_base[jj, iz] * wage_factor * y_noise[i]89 ih = h[i]90 reb = rebate if j == 1 else 0.091 best = NEG92 ch = -193 jb = 0; wb = 0.0; ja = 0; wa = 0.094 hh = 0.0; bb = 0.0; nn = 1; ic = i_new; pay = 0.095 x_stay = 0.0; x_refi = 0.0; x_sell = 0.096 c = c_min; s = s_min; ap = 0.097 if ih >= 0:98 hh = h_grid[ih]99 bb = b[i]100 nn = n[i]101 ic = i_new if reset else ic_by_n[nn - 1]102 pay = bb * _af_due(ic, nn)103 base = y + a[i] * R_a104 x_stay = base - pay - m * p * hh105 x_refi = base - bb - (kappa_f + m) * p * hh106 x_sell = base + (1.0 - kappa_s) * p * hh - bb107 # stay (bilinear in b, a)108 jb, wb = _locate(b_grid, bb)109 ja, wa = _locate(a_grid, a[i])110 v00 = W_S[jj, iz, ih, jb, nn - 1, ja]; v01 = W_S[jj, iz, ih, jb, nn - 1, ja + 1]111 v10 = W_S[jj, iz, ih, jb + 1, nn - 1, ja]; v11 = W_S[jj, iz, ih, jb + 1, nn - 1, ja + 1]112 vS = _lin(_lin(v00, v01, wa), _lin(v10, v11, wa), wb)113 if vS > best:114 best = vS115 ch = 0116 if x_refi > 0.0:117 jx, wx = _locate(x_grid, x_refi)118 vF = _v1(W_F[jj, iz, ih], jx, wx)119 if vF > best:120 best = vF121 ch = 1122 if x_sell > 0.0:123 jx, wx = _locate(x_grid, x_sell)124 vB = _v1(W_B[jj, iz], jx, wx)125 if vB > best:126 best = vB127 ch = 2128 vR = _v1(W_R[jj, iz], jx, wx)129 if vR > best:130 best = vR131 ch = 3132 x_used = x_sell133 else:134 x = y + a[i] * R_a + reb135 x_used = x136 jx, wx = _locate(x_grid, x)137 vB = _v1(W_B[jj, iz], jx, wx)138 vR = _v1(W_Rs[jj, iz], jx, wx)139 if vR > best:140 best = vR141 ch = 3142 if vB > best:143 best = vB144 ch = 2145 flag = 0146 g = 0.0147 pay_out = 0.0148 Bnew = 0.0149 if ch == 0:150 # ---- stay151 g = (1.0 - wb) * ((1.0 - wa) * gS[jj, iz, ih, jb, nn - 1, ja] + wa * gS[jj, iz, ih, jb, nn - 1, ja + 1]) \152 + wb * ((1.0 - wa) * gS[jj, iz, ih, jb + 1, nn - 1, ja] + wa * gS[jj, iz, ih, jb + 1, nn - 1, ja + 1])153 ap = (1.0 - wb) * ((1.0 - wa) * aS[jj, iz, ih, jb, nn - 1, ja] + wa * aS[jj, iz, ih, jb, nn - 1, ja + 1]) \154 + wb * ((1.0 - wa) * aS[jj, iz, ih, jb + 1, nn - 1, ja] + wa * aS[jj, iz, ih, jb + 1, nn - 1, ja + 1])155 c = x_stay - g - ap156 if c < c_min:157 ap = max(x_stay - g - c_min, 0.0)158 c = x_stay - g - ap159 if c < c_min:160 g = max(x_stay - c_min, 0.0)161 c = x_stay - g162 flag = 2163 s = (1.0 + chi) * hh164 pay_out = pay165 bn = max(bb - pay, 0.0) * (1.0 + ic) / (1.0 + pi_next)166 nn_next = nn - 1 if nn > 1 else 1167 if bn <= 1e-12:168 bn = 0.0169 out_h[i] = ih; out_b[i] = bn; out_n[i] = nn_next; out_a[i] = ap170 elif ch == 1:171 # ---- refinance172 jx, wx = _locate(x_grid, x_refi)173 g = _p1(gF[jj, iz, ih], jx, wx)174 Bn = _p1(BF[jj, iz, ih], jx, wx)175 ap = _p1(aF[jj, iz, ih], jx, wx)176 payn = Bn * af_new177 c = x_refi + Bn - payn - g - ap178 if c < c_min:179 ap = max(x_refi + Bn - payn - g - c_min, 0.0)180 c = x_refi + Bn - payn - g - ap181 flag = 2182 s = (1.0 + chi) * hh183 pay_out = payn184 Bnew = Bn185 bn = max(Bn - payn, 0.0) * (1.0 + i_new) / (1.0 + pi_next)186 out_h[i] = ih; out_b[i] = bn; out_n[i] = N - 1; out_a[i] = ap187 elif ch == 2:188 # ---- buy: house size by interpolated per-size values, loan / saving / gift at the nearest grid point189 jx, wx = _locate(x_grid, x_used)190 kx = jx + 1 if wx > 0.5 else jx191 ihn = -1192 vbest = NEG193 for ihh in range(h_grid.shape[0]):194 vv = _lin(W_Bh[jj, iz, jx, ihh], W_Bh[jj, iz, jx + 1, ihh], wx)195 if vv > vbest:196 vbest = vv197 ihn = ihh198 if ihn < 0:199 ihn = 0200 if W_Bh[jj, iz, kx, ihn] <= NEG * 0.5:201 kx = jx + 1202 Bn = BBh[jj, iz, kx, ihn]203 ap = aBh[jj, iz, kx, ihn]204 g = gBh[jj, iz, kx, ihn]205 hh2 = h_grid[ihn]206 payn = Bn * af_new207 c = x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - ap208 if c < c_min:209 ap = max(x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - c_min, 0.0)210 c = x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - ap211 flag = 2212 if c < c_min:213 # fall back to renting with the cash available214 ch = 3215 if ch == 2:216 s = (1.0 + chi) * hh2217 pay_out = payn218 Bnew = Bn219 bn = max(Bn - payn, 0.0) * (1.0 + i_new) / (1.0 + pi_next)220 out_h[i] = ihn; out_b[i] = bn; out_n[i] = N - 1; out_a[i] = ap221 if ch == 3:222 # ---- rent (sitting tenant if renter already, new tenant if seller)223 jx, wx = _locate(x_grid, x_used)224 if ih >= 0:225 ap = _p1(aR[jj, iz], jx, wx)226 g = _p1(gR[jj, iz], jx, wx)227 rr = rho228 else:229 ap = _p1(aRs[jj, iz], jx, wx)230 g = _p1(gRs[jj, iz], jx, wx)231 rr = rho_sit232 e = x_used - g - ap233 c, s = _renter_cs(e, rr, alpha, s_min, s_max)234 if c < c_min:235 ap = 0.0236 e = x_used - g237 c, s = _renter_cs(e, rr, alpha, s_min, s_max)238 flag = 2239 if c < c_min:240 g = 0.0241 e = x_used242 c, s = _renter_cs(e, rr, alpha, s_min, s_max)243 out_h[i] = -1; out_b[i] = 0.0; out_n[i] = 1; out_a[i] = ap244 if ch == -1:245 # ---- no feasible branch: forced sale, consumption floor, lender absorbs the shortfall246 flag = 1247 s = s_min248 c = c_min249 g = 0.0250 out_h[i] = -1; out_b[i] = 0.0; out_n[i] = 1; out_a[i] = 0.0251 if g < 0.0:252 g = 0.0253 if j != j_gift:254 g = 0.0255 out_choice[i] = ch256 out_c[i] = c257 out_s[i] = s258 out_g[i] = g259 out_pay[i] = pay_out260 out_Bnew[i] = Bnew261 out_flag[i] = flag262 out_u[i] = _u(max(c, 1e-9), max(s, 1e-9), alpha, gamma)263 out_x[i] = x_used264 if j == J:265 hh3 = h_grid[out_h[i]] if out_h[i] >= 0 else 0.0266 gross = out_a[i] + (1.0 - kappa_s) * p_next * hh3 - out_b[i]267 if gross < 0.0:268 gross = 0.0269 out_estate[i] = gross270 else:271 out_estate[i] = 0.0272273274@njit(cache=True)275def _draw_z(z, age, u, cumPi, j_retire):276 """In-place Markov draw for ages 2..j_retire-1 (frozen afterwards; age 1 drawn at birth)."""277 M = z.shape[0]278 nz = cumPi.shape[1]279 for i in range(M):280 j = age[i]281 if j >= 2 and j < j_retire:282 r = u[i]283 zi = z[i]284 k = 0285 while k < nz - 1 and r > cumPi[zi, k]:286 k += 1287 z[i] = k288289290@njit(cache=True)291def _draw_z1(z_parent, u, cumPg):292 M = z_parent.shape[0]293 nz = cumPg.shape[1]294 out = np.empty(M, dtype=np.int64)295 for i in range(M):296 r = u[i]297 zi = z_parent[i]298 k = 0299 while k < nz - 1 and r > cumPg[zi, k]:300 k += 1301 out[i] = k302 return out303304305# --------------------------------------------------------------------------------------306# Panel state307# --------------------------------------------------------------------------------------308@dataclass309class Panel:310 """Dynasty slots: young (ages 1-6) and old (ages 7-12) members, start-of-period states."""311 phase: np.ndarray # 0..5312 age_y: np.ndarray; z_y: np.ndarray; h_y: np.ndarray; b_y: np.ndarray; n_y: np.ndarray; a_y: np.ndarray313 age_o: np.ndarray; z_o: np.ndarray; h_o: np.ndarray; b_o: np.ndarray; n_o: np.ndarray; a_o: np.ndarray314 estate_pending: np.ndarray # estate to be received by the young member at age 7315 gift_received: np.ndarray # gift (net) received at birth by the current young member316 parent_owner_at_gift: np.ndarray # parent tenure (1 owner) at the time of the gift (for mobility)317 parent_z_at_gift: np.ndarray318 parent_nw_at_gift: np.ndarray319 estate_received: np.ndarray # estate (net) received by the current old member at age 7320 parent_owner_old: np.ndarray = None # parent tenure at gift, carried after promotion321 parent_nw_old: np.ndarray = None322 parent_z_old: np.ndarray = None323 gift_received_old: np.ndarray = None324 t: int = 0325326 @classmethod327 def initial(cls, P: Params, rng: np.random.Generator) -> "Panel":328 M = P.n_dyn329 phase = np.arange(M) % 6330 z_y = rng.choice(P.n_z, size=M, p=P.pi_z)331 z_o = rng.choice(P.n_z, size=M, p=P.pi_z)332 zeros = np.zeros(M)333 return cls(phase=phase,334 age_y=1 + phase, z_y=z_y, h_y=np.full(M, -1), b_y=zeros.copy(), n_y=np.ones(M, dtype=np.int64), a_y=zeros.copy(),335 age_o=7 + phase, z_o=z_o, h_o=np.full(M, -1), b_o=zeros.copy(), n_o=np.ones(M, dtype=np.int64), a_o=zeros.copy(),336 estate_pending=zeros.copy(), gift_received=zeros.copy(), parent_owner_at_gift=np.zeros(M, dtype=np.int64),337 parent_z_at_gift=z_o.copy(), parent_nw_at_gift=zeros.copy(), estate_received=zeros.copy(),338 parent_owner_old=np.zeros(M, dtype=np.int64), parent_nw_old=zeros.copy(), parent_z_old=z_o.copy(),339 gift_received_old=zeros.copy(), t=0)340341 def copy(self) -> "Panel":342 kw = {k: (v.copy() if isinstance(v, np.ndarray) else v) for k, v in self.__dict__.items()}343 return Panel(**kw)344345346@dataclass347class PeriodRecord:348 """Decisions of one period for all households (young + old concatenated)."""349 age: np.ndarray; z: np.ndarray; h0: np.ndarray; b0: np.ndarray; a0: np.ndarray; n0: np.ndarray350 choice: np.ndarray; c: np.ndarray; s: np.ndarray; h1: np.ndarray; b1: np.ndarray; n1: np.ndarray; a1: np.ndarray351 gift: np.ndarray; pay: np.ndarray; Bnew: np.ndarray; estate: np.ndarray; flag: np.ndarray; u: np.ndarray; x: np.ndarray352 y_shock: np.ndarray353 gift_received: np.ndarray; parent_owner: np.ndarray; parent_nw: np.ndarray; estate_received: np.ndarray354 slot: np.ndarray; is_old: np.ndarray355 p: float; p_r: float; rho: float356357 def owner0(self) -> np.ndarray:358 return self.h0 >= 0359360 def owner1(self) -> np.ndarray:361 return self.h1 >= 0362363 def net_worth0(self, P: Params) -> np.ndarray:364 hv = np.where(self.h0 >= 0, P.h_grid[np.maximum(self.h0, 0)], 0.0)365 return self.a0 + (1.0 - P.kappa_s) * self.p * hv - self.b0366367 def net_worth1(self, P: Params) -> np.ndarray:368 hv = np.where(self.h1 >= 0, P.h_grid[np.maximum(self.h1, 0)], 0.0)369 return self.a1 + (1.0 - P.kappa_s) * self.p * hv - self.b1370371 def housing_demand(self, P: Params) -> float:372 """Owner-occupied units demanded per unit mass of households of one age."""373 hv = np.where(self.h1 >= 0, P.h_grid[np.maximum(self.h1, 0)], 0.0)374 return float(hv.sum()) / (len(self.age) / P.J)375376 def rental_demand(self, P: Params) -> float:377 return float(np.where(self.h1 < 0, self.s, 0.0).sum()) / (len(self.age) / P.J)378379380def _kernel_args(P: Params, S: PolicyStack, pr: Prices):381 ic_by_n = pr.i_c_by_n if pr.i_c_by_n is not None else np.full(P.N, pr.i_new)382 return (S.W_S, S.aS, S.gS, S.W_F, S.BF, S.aF, S.gF, S.W_B, S.hB, S.BB, S.aB, S.gB, S.W_R, S.aR, S.gR,383 S.W_Rs, S.aRs, S.gRs, S.W_Bh, S.BBh, S.aBh, S.gBh, P.a_grid, P.b_grid, P.x_grid, P.h_grid, np.ascontiguousarray(ic_by_n, dtype=float),384 bool(pr.reset), pr.p, (pr.p_next if pr.p_next is not None else pr.p), pr.rho,385 pr.rho * pr.rent_rigid_factor, pr.R_a, pr.i_new, pr.pi_next, pr.transfer_rebate, pr.wage_factor,386 P.kappa_b, P.kappa_s, P.kappa_f, P.m, P.chi, P.alpha, P.gamma, P.s_min, P.s_max, P.c_min, P.tau_b,387 P.N, P.j_gift, P.J)388389390def _run_group(P, S, pr, y_base, age, z, h, b, n, a, y_noise=None):391 M = len(age)392 if y_noise is None:393 y_noise = np.ones(M)394 out = dict(choice=np.empty(M, dtype=np.int64), c=np.empty(M), s=np.empty(M), h=np.empty(M, dtype=np.int64),395 b=np.empty(M), n=np.empty(M, dtype=np.int64), a=np.empty(M), g=np.empty(M), pay=np.empty(M),396 Bnew=np.empty(M), estate=np.empty(M), flag=np.empty(M, dtype=np.int64), u=np.empty(M), x=np.empty(M))397 _sim_step(age, z, h, b, n, a, y_base, np.ascontiguousarray(y_noise, dtype=float), *_kernel_args(P, S, pr),398 out["choice"], out["c"], out["s"], out["h"], out["b"], out["n"], out["a"], out["g"], out["pay"],399 out["Bnew"], out["estate"], out["flag"], out["u"], out["x"])400 return out401402403def simulate_period(P: Params, S: PolicyStack, pr: Prices, panel: Panel, rng_key: tuple,404 balance_factor: float = 1.0, record: bool = True) -> PeriodRecord | None:405 """Advance the panel by one period under policies S and prices pr; returns the record.406407 balance_factor multiplies all real mortgage balances at the start of the period (surprise408 inflation erosion 1/(1+pi_u), undone for the ARM contract); deposits are revalued through409 pr.R_a. rng_key seeds the period's draws (common random numbers).410 """411 M = P.n_dyn412 rng = np.random.default_rng(list(rng_key))413 U = rng.uniform(size=(3, M))414 if P.sigma_trans > 0:415 Nz = rng.standard_normal(size=(2, M))416 noise_o = np.exp(P.sigma_trans * Nz[0] - 0.5 * P.sigma_trans**2)417 noise_y = np.exp(P.sigma_trans * Nz[1] - 0.5 * P.sigma_trans**2)418 else:419 noise_o = np.ones(M); noise_y = np.ones(M)420 y_base = P.income_matrix()421 cumPi = np.cumsum(P.Pi_z, axis=1)422 cumPg = np.cumsum(P.Pi_g, axis=1)423 if balance_factor != 1.0:424 panel.b_o *= balance_factor425 panel.b_y *= balance_factor426 # income shocks (ages >= 2)427 _draw_z(panel.z_o, panel.age_o, U[0], cumPi, P.j_retire)428 _draw_z(panel.z_y, panel.age_y, U[1], cumPi, P.j_retire)429 # --- old members decide first430 o = _run_group(P, S, pr, y_base, panel.age_o, panel.z_o, panel.h_o, panel.b_o, panel.n_o, panel.a_o, noise_o)431 # newborns (young at age 1) receive the gift of the old member at age 7 (same slot)432 newborn = panel.age_y == 1433 if newborn.any():434 gift_net = (1.0 - P.tau_b) * o["g"]435 panel.gift_received[newborn] = gift_net[newborn]436 panel.a_y[newborn] = gift_net[newborn] + pr.transfer_rebate437 nw_parent = panel.a_o + np.where(panel.h_o >= 0, (1 - P.kappa_s) * pr.p * P.h_grid[np.maximum(panel.h_o, 0)], 0.0) - panel.b_o438 panel.parent_owner_at_gift[newborn] = (o["h"][newborn] >= 0).astype(np.int64)439 panel.parent_z_at_gift[newborn] = panel.z_o[newborn]440 panel.parent_nw_at_gift[newborn] = nw_parent[newborn]441 panel.z_y[newborn] = _draw_z1(panel.z_o, U[2], cumPg)[newborn]442 panel.h_y[newborn] = -1; panel.b_y[newborn] = 0.0; panel.n_y[newborn] = 1443 yg = _run_group(P, S, pr, y_base, panel.age_y, panel.z_y, panel.h_y, panel.b_y, panel.n_y, panel.a_y, noise_y)444 rec = None if not record else PeriodRecord(445 age=np.concatenate([panel.age_y, panel.age_o]), z=np.concatenate([panel.z_y, panel.z_o]),446 h0=np.concatenate([panel.h_y, panel.h_o]), b0=np.concatenate([panel.b_y, panel.b_o]),447 a0=np.concatenate([panel.a_y, panel.a_o]), n0=np.concatenate([panel.n_y, panel.n_o]),448 choice=np.concatenate([yg["choice"], o["choice"]]), c=np.concatenate([yg["c"], o["c"]]),449 s=np.concatenate([yg["s"], o["s"]]), h1=np.concatenate([yg["h"], o["h"]]), b1=np.concatenate([yg["b"], o["b"]]),450 n1=np.concatenate([yg["n"], o["n"]]), a1=np.concatenate([yg["a"], o["a"]]), gift=np.concatenate([yg["g"], o["g"]]),451 pay=np.concatenate([yg["pay"], o["pay"]]), Bnew=np.concatenate([yg["Bnew"], o["Bnew"]]),452 estate=np.concatenate([yg["estate"], o["estate"]]), flag=np.concatenate([yg["flag"], o["flag"]]),453 u=np.concatenate([yg["u"], o["u"]]), x=np.concatenate([yg["x"], o["x"]]),454 y_shock=np.concatenate([noise_y, noise_o]),455 gift_received=np.concatenate([panel.gift_received, panel.gift_received_old]),456 parent_owner=np.concatenate([panel.parent_owner_at_gift, panel.parent_owner_old]),457 parent_nw=np.concatenate([panel.parent_nw_at_gift, panel.parent_nw_old]),458 estate_received=np.concatenate([np.full(M, np.nan), panel.estate_received]),459 slot=np.concatenate([np.arange(M), np.arange(M)]), is_old=np.concatenate([np.zeros(M, bool), np.ones(M, bool)]),460 p=pr.p, p_r=pr.p_r, rho=pr.rho)461 # --- advance states462 panel.h_y, panel.b_y, panel.n_y, panel.a_y = yg["h"], yg["b"], yg["n"], yg["a"]463 panel.h_o, panel.b_o, panel.n_o, panel.a_o = o["h"], o["b"], o["n"], o["a"]464 dying = panel.age_o == P.J465 panel.estate_pending[dying] = (1.0 - P.tau_b) * o["estate"][dying]466 panel.age_y += 1467 panel.age_o += 1468 # young at 7 -> old; estate received now (start of next period); new young born next period469 promote = panel.age_y == P.j_estate470 if promote.any():471 panel.age_o[promote] = P.j_estate472 panel.z_o[promote] = panel.z_y[promote]473 panel.h_o[promote] = panel.h_y[promote]; panel.b_o[promote] = panel.b_y[promote]474 panel.n_o[promote] = panel.n_y[promote]475 panel.a_o[promote] = panel.a_y[promote] + panel.estate_pending[promote]476 panel.estate_received[promote] = panel.estate_pending[promote]477 panel.estate_pending[promote] = 0.0478 panel.parent_owner_old[promote] = panel.parent_owner_at_gift[promote]479 panel.parent_nw_old[promote] = panel.parent_nw_at_gift[promote]480 panel.parent_z_old[promote] = panel.parent_z_at_gift[promote]481 panel.gift_received_old[promote] = panel.gift_received[promote]482 panel.age_y[promote] = 1483 panel.h_y[promote] = -1; panel.b_y[promote] = 0.0; panel.n_y[promote] = 1; panel.a_y[promote] = 0.0484 panel.t += 1485 return rec486487488# --------------------------------------------------------------------------------------489# Stationary simulation and aggregation490# --------------------------------------------------------------------------------------491def simulate_stationary(P: Params, pols: list[AgePolicy], pr: Prices, n_record: int = 6,492 panel0: Panel | None = None, burn: int | None = None) -> tuple[list[PeriodRecord], Panel]:493 """Simulate `burn` periods from an arbitrary initial condition under stationary policies,494 then record `n_record` further periods (one full phase cycle = 6). Returns the records and495 the panel state at the end (start-of-period states for a subsequent transition)."""496 S = stack_policies(pols)497 burn = P.burn if burn is None else burn498 rng = np.random.default_rng(P.seed)499 panel = Panel.initial(P, rng) if panel0 is None else panel0.copy()500 recs: list[PeriodRecord] = []501 for t in range(burn + n_record):502 rec = simulate_period(P, S, pr, panel, (P.seed, t), record=(t >= burn))503 if t >= burn:504 recs.append(rec)505 return recs, panel506507508def pool_records(recs: list[PeriodRecord]) -> PeriodRecord:509 kw = {}510 for k, v in recs[0].__dict__.items():511 if isinstance(v, np.ndarray):512 kw[k] = np.concatenate([getattr(r, k) for r in recs])513 else:514 kw[k] = v515 return PeriodRecord(**kw)516517518def estate_lottery(P: Params, rec: PeriodRecord) -> tuple[np.ndarray, np.ndarray]:519 """K-point lottery: conditional means of net estates received within K equal-probability bins."""520 e = rec.estate_received[rec.is_old & (rec.age == P.j_estate)]521 e = e[np.isfinite(e)]522 if len(e) == 0:523 return np.zeros(P.K_estate), np.full(P.K_estate, 1.0 / P.K_estate)524 e = np.sort(e)525 bins = np.array_split(e, P.K_estate)526 pts = np.array([b.mean() if len(b) else 0.0 for b in bins])527 return pts, np.full(P.K_estate, 1.0 / P.K_estate)528529530def gini(x: np.ndarray) -> float:531 x = np.asarray(x, dtype=float)532 x = x[np.isfinite(x)]533 if len(x) == 0:534 return np.nan535 x = np.sort(x - min(x.min(), 0.0))536 n = len(x)537 if x.sum() <= 0:538 return np.nan539 return float((2.0 * np.sum(np.arange(1, n + 1) * x) / (n * x.sum())) - (n + 1.0) / n)540541542def rank_rank_slope(parent: np.ndarray, child: np.ndarray) -> float:543 ok = np.isfinite(parent) & np.isfinite(child)544 if ok.sum() < 10:545 return np.nan546 rp = np.argsort(np.argsort(parent[ok])) / ok.sum()547 rc = np.argsort(np.argsort(child[ok])) / ok.sum()548 return float(np.polyfit(rp, rc, 1)[0])549550551# --------------------------------------------------------------------------------------552# Cohort simulations (welfare decomposition) and deterministic paths (experiment 2)553# --------------------------------------------------------------------------------------554def simulate_cohort(P: Params, S: PolicyStack, pr: Prices, z1: np.ndarray, a1: np.ndarray, lottery,555 seed: int = 11, fixed_z: bool = False) -> dict:556 """Simulate a cohort from age 1 to J (no dynasty link; estate drawn from `lottery` at557 j_estate). Returns discounted flow utility U, warm-glow part Phi, and per-age arrays.558 fixed_z=True keeps z constant (equal-income deterministic design)."""559 n = len(a1)560 rng = np.random.default_rng(seed)561 y_base = P.income_matrix()562 cumPi = np.cumsum(P.Pi_z, axis=1)563 z = z1.astype(np.int64).copy(); a = a1.astype(float).copy()564 h = np.full(n, -1, dtype=np.int64); b = np.zeros(n); nn = np.ones(n, dtype=np.int64)565 U = np.zeros(n); Phi = np.zeros(n)566 out = dict(own=np.zeros((P.J, n)), nw=np.zeros((P.J, n)), a=np.zeros((P.J, n)), b=np.zeros((P.J, n)),567 c=np.zeros((P.J, n)), s=np.zeros((P.J, n)), gift=np.zeros(n), estate=np.zeros(n), h=np.zeros((P.J, n), dtype=np.int64))568 pts, w = lottery569 for j in range(1, P.J + 1):570 age = np.full(n, j, dtype=np.int64)571 if j >= 2 and not fixed_z:572 _draw_z(z, age, rng.uniform(size=n), cumPi, P.j_retire)573 if j == P.j_estate:574 k = rng.choice(len(pts), size=n, p=w)575 a = a + pts[k]576 out["estate"] = pts[k]577 o = _run_group(P, S, pr, y_base, age, z, h, b, nn, a)578 disc = P.beta ** (j - 1)579 U += disc * o["u"]580 if j == P.j_gift:581 Phi += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["g"]))582 out["gift"] = o["g"]583 if j == P.J:584 Phi += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["estate"]))585 hv = np.where(o["h"] >= 0, P.h_grid[np.maximum(o["h"], 0)], 0.0)586 out["own"][j - 1] = o["h"] >= 0587 out["nw"][j - 1] = o["a"] + (1 - P.kappa_s) * pr.p * hv - o["b"]588 out["a"][j - 1] = o["a"]; out["b"][j - 1] = o["b"]; out["c"][j - 1] = o["c"]; out["s"][j - 1] = o["s"]589 out["h"][j - 1] = o["h"]590 h, b, nn, a = o["h"], o["b"], o["n"], o["a"]591 out["U"] = U; out["Phi"] = Phi; out["V"] = U + Phi592 return out593594595def _wg(P: Params, x):596 from .model import warm_glow597 return warm_glow(P, np.asarray(x, dtype=float))598599600def remaining_lifetime(P: Params, S: PolicyStack, pr: Prices, rec: PeriodRecord, lottery, seed: int = 12,601 start_states: dict | None = None) -> tuple[np.ndarray, np.ndarray]:602 """Discounted flow utility U and warm-glow Phi over the remaining lifetime of every household603 in `rec` (start-of-period states, before this period's decision), under policies S and604 prices pr (stationary). Used for the welfare decomposition of households alive at a shock."""605 n = len(rec.age)606 rng = np.random.default_rng(seed)607 y_base = P.income_matrix()608 cumPi = np.cumsum(P.Pi_z, axis=1)609 st = start_states or {}610 age = rec.age.copy(); z = rec.z.copy(); h = st.get("h", rec.h0).copy(); b = st.get("b", rec.b0).copy()611 nn = rec.n0.copy(); a = st.get("a", rec.a0).copy()612 U = np.zeros(n); Phi = np.zeros(n)613 pts, w = lottery614 alive = np.ones(n, dtype=bool)615 step = 0616 while alive.any():617 idx = np.where(alive)[0]618 if step > 0:619 _draw_z(z, age, rng.uniform(size=n), cumPi, P.j_retire)620 est = age == P.j_estate621 if est.any():622 k = rng.choice(len(pts), size=int(est.sum()), p=w)623 a[est] += pts[k]624 o = _run_group(P, S, pr, y_base, age[idx], z[idx], h[idx], b[idx], nn[idx], a[idx])625 disc = P.beta ** step626 U[idx] += disc * o["u"]627 g_age = age[idx] == P.j_gift628 if g_age.any():629 Phi[idx[g_age]] += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["g"][g_age]))630 last = age[idx] == P.J631 if last.any():632 Phi[idx[last]] += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["estate"][last]))633 h[idx], b[idx], nn[idx], a[idx] = o["h"], o["b"], o["n"], o["a"]634 age[idx] += 1635 alive = age <= P.J636 step += 1637 return U, Phi638