"""WP18 — cross-sectional distribution by Monte-Carlo simulation of linked dynasties. Why Monte Carlo and not a Young (2010) histogram: the dynasty link pays the parent's gift to the child at the child's age 1 and the parent's estate at the child's age j_estate = 7; the two transfers are correlated through the parent's state, so the invariant distribution lives on the joint parent-child state, which a histogram on the individual state grid does not carry. A panel of `n_dyn` dynasty slots, each holding one "young" member (ages 1-6) and one "old" member (ages 7-12) at all times, tracks the link exactly. Random numbers are common across price iterations (seeded by (seed, period)) so that excess demand is a deterministic function of prices; sample-size sensitivity is reported in experiment 10. Slot phase s in {0..5}: young age = 1 + s, old age = 7 + s. When the old member dies (after age J) the young member becomes the old member (age 7) and receives the estate; a new young member (age 1) is born in the same period and receives the gift chosen by the new old member. Within a period, old members decide first (their gift is the newborn's initial asset). """ from __future__ import annotations from dataclasses import dataclass, field import numpy as np from numba import njit, prange from .household_solver import NEG, AgePolicy, _lin, _locate, _af_due, _renter_cs, _u from .model import Params, Prices CHOICE_NAMES = {-1: "forced_sale", 0: "stay", 1: "refinance", 2: "buy", 3: "rent"} # -------------------------------------------------------------------------------------- # Stacked policies (ages along axis 0) for the numba kernel # -------------------------------------------------------------------------------------- @dataclass class PolicyStack: W_S: np.ndarray; aS: np.ndarray; gS: np.ndarray W_F: np.ndarray; BF: np.ndarray; aF: np.ndarray; gF: np.ndarray W_B: np.ndarray; hB: np.ndarray; BB: np.ndarray; aB: np.ndarray; gB: np.ndarray W_R: np.ndarray; aR: np.ndarray; gR: np.ndarray W_Rs: np.ndarray; aRs: np.ndarray; gRs: np.ndarray W_Bh: np.ndarray = None; BBh: np.ndarray = None; aBh: np.ndarray = None; gBh: np.ndarray = None def stack_policies(pols: list[AgePolicy], dtype=np.float64) -> PolicyStack: def st(name, dt=dtype): return np.ascontiguousarray(np.stack([getattr(p, name) for p in pols]), dtype=dt) return PolicyStack(W_S=st("W_S"), aS=st("aS"), gS=st("gS"), W_F=st("W_F"), BF=st("BF"), aF=st("aF"), gF=st("gF"), W_B=st("W_B"), hB=st("hB", np.int64), BB=st("BB"), aB=st("aB"), gB=st("gB"), W_R=st("W_R"), aR=st("aR"), gR=st("gR"), W_Rs=st("W_Rs"), aRs=st("aRs"), gRs=st("gRs"), W_Bh=st("W_Bh"), BBh=st("BBh"), aBh=st("aBh"), gBh=st("gBh")) # -------------------------------------------------------------------------------------- # One-period decision kernel # -------------------------------------------------------------------------------------- @njit(cache=True, inline="always") def _ix(x_grid, x): j, w = _locate(x_grid, x) return j, w @njit(cache=True, inline="always") def _v1(arr, j, w): return _lin(arr[j], arr[j + 1], w) @njit(cache=True, inline="always") def _p1(arr, j, w): """Policy interpolation (linear, no infeasibility sentinel).""" return (1.0 - w) * arr[j] + w * arr[j + 1] @njit(cache=True, parallel=True) def _sim_step(age, z, h, b, n, a, y_base, y_noise, 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, a_grid, b_grid, x_grid, h_grid, ic_by_n, reset, p, p_next, rho, rho_sit, R_a, i_new, pi_next, rebate, wage_factor, kappa_b, kappa_s, kappa_f, m, chi, alpha, gamma, s_min, s_max, c_min, tau_b, N, j_gift, J, out_choice, out_c, out_s, out_h, out_b, out_n, out_a, out_g, out_pay, out_Bnew, out_estate, out_flag, out_u, out_x): M = age.shape[0] af_new = _af_due(i_new, N) for i in prange(M): j = age[i] jj = j - 1 iz = z[i] y = y_base[jj, iz] * wage_factor * y_noise[i] ih = h[i] reb = rebate if j == 1 else 0.0 best = NEG ch = -1 jb = 0; wb = 0.0; ja = 0; wa = 0.0 hh = 0.0; bb = 0.0; nn = 1; ic = i_new; pay = 0.0 x_stay = 0.0; x_refi = 0.0; x_sell = 0.0 c = c_min; s = s_min; ap = 0.0 if ih >= 0: hh = h_grid[ih] bb = b[i] nn = n[i] ic = i_new if reset else ic_by_n[nn - 1] pay = bb * _af_due(ic, nn) base = y + a[i] * R_a x_stay = base - pay - m * p * hh x_refi = base - bb - (kappa_f + m) * p * hh x_sell = base + (1.0 - kappa_s) * p * hh - bb # stay (bilinear in b, a) jb, wb = _locate(b_grid, bb) ja, wa = _locate(a_grid, a[i]) v00 = W_S[jj, iz, ih, jb, nn - 1, ja]; v01 = W_S[jj, iz, ih, jb, nn - 1, ja + 1] v10 = W_S[jj, iz, ih, jb + 1, nn - 1, ja]; v11 = W_S[jj, iz, ih, jb + 1, nn - 1, ja + 1] vS = _lin(_lin(v00, v01, wa), _lin(v10, v11, wa), wb) if vS > best: best = vS ch = 0 if x_refi > 0.0: jx, wx = _locate(x_grid, x_refi) vF = _v1(W_F[jj, iz, ih], jx, wx) if vF > best: best = vF ch = 1 if x_sell > 0.0: jx, wx = _locate(x_grid, x_sell) vB = _v1(W_B[jj, iz], jx, wx) if vB > best: best = vB ch = 2 vR = _v1(W_R[jj, iz], jx, wx) if vR > best: best = vR ch = 3 x_used = x_sell else: x = y + a[i] * R_a + reb x_used = x jx, wx = _locate(x_grid, x) vB = _v1(W_B[jj, iz], jx, wx) vR = _v1(W_Rs[jj, iz], jx, wx) if vR > best: best = vR ch = 3 if vB > best: best = vB ch = 2 flag = 0 g = 0.0 pay_out = 0.0 Bnew = 0.0 if ch == 0: # ---- stay 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]) \ + wb * ((1.0 - wa) * gS[jj, iz, ih, jb + 1, nn - 1, ja] + wa * gS[jj, iz, ih, jb + 1, nn - 1, ja + 1]) 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]) \ + wb * ((1.0 - wa) * aS[jj, iz, ih, jb + 1, nn - 1, ja] + wa * aS[jj, iz, ih, jb + 1, nn - 1, ja + 1]) c = x_stay - g - ap if c < c_min: ap = max(x_stay - g - c_min, 0.0) c = x_stay - g - ap if c < c_min: g = max(x_stay - c_min, 0.0) c = x_stay - g flag = 2 s = (1.0 + chi) * hh pay_out = pay bn = max(bb - pay, 0.0) * (1.0 + ic) / (1.0 + pi_next) nn_next = nn - 1 if nn > 1 else 1 if bn <= 1e-12: bn = 0.0 out_h[i] = ih; out_b[i] = bn; out_n[i] = nn_next; out_a[i] = ap elif ch == 1: # ---- refinance jx, wx = _locate(x_grid, x_refi) g = _p1(gF[jj, iz, ih], jx, wx) Bn = _p1(BF[jj, iz, ih], jx, wx) ap = _p1(aF[jj, iz, ih], jx, wx) payn = Bn * af_new c = x_refi + Bn - payn - g - ap if c < c_min: ap = max(x_refi + Bn - payn - g - c_min, 0.0) c = x_refi + Bn - payn - g - ap flag = 2 s = (1.0 + chi) * hh pay_out = payn Bnew = Bn bn = max(Bn - payn, 0.0) * (1.0 + i_new) / (1.0 + pi_next) out_h[i] = ih; out_b[i] = bn; out_n[i] = N - 1; out_a[i] = ap elif ch == 2: # ---- buy: house size by interpolated per-size values, loan / saving / gift at the nearest grid point jx, wx = _locate(x_grid, x_used) kx = jx + 1 if wx > 0.5 else jx ihn = -1 vbest = NEG for ihh in range(h_grid.shape[0]): vv = _lin(W_Bh[jj, iz, jx, ihh], W_Bh[jj, iz, jx + 1, ihh], wx) if vv > vbest: vbest = vv ihn = ihh if ihn < 0: ihn = 0 if W_Bh[jj, iz, kx, ihn] <= NEG * 0.5: kx = jx + 1 Bn = BBh[jj, iz, kx, ihn] ap = aBh[jj, iz, kx, ihn] g = gBh[jj, iz, kx, ihn] hh2 = h_grid[ihn] payn = Bn * af_new c = x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - ap if c < c_min: ap = max(x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - c_min, 0.0) c = x_used - (1.0 + kappa_b + m) * p * hh2 + Bn - payn - g - ap flag = 2 if c < c_min: # fall back to renting with the cash available ch = 3 if ch == 2: s = (1.0 + chi) * hh2 pay_out = payn Bnew = Bn bn = max(Bn - payn, 0.0) * (1.0 + i_new) / (1.0 + pi_next) out_h[i] = ihn; out_b[i] = bn; out_n[i] = N - 1; out_a[i] = ap if ch == 3: # ---- rent (sitting tenant if renter already, new tenant if seller) jx, wx = _locate(x_grid, x_used) if ih >= 0: ap = _p1(aR[jj, iz], jx, wx) g = _p1(gR[jj, iz], jx, wx) rr = rho else: ap = _p1(aRs[jj, iz], jx, wx) g = _p1(gRs[jj, iz], jx, wx) rr = rho_sit e = x_used - g - ap c, s = _renter_cs(e, rr, alpha, s_min, s_max) if c < c_min: ap = 0.0 e = x_used - g c, s = _renter_cs(e, rr, alpha, s_min, s_max) flag = 2 if c < c_min: g = 0.0 e = x_used c, s = _renter_cs(e, rr, alpha, s_min, s_max) out_h[i] = -1; out_b[i] = 0.0; out_n[i] = 1; out_a[i] = ap if ch == -1: # ---- no feasible branch: forced sale, consumption floor, lender absorbs the shortfall flag = 1 s = s_min c = c_min g = 0.0 out_h[i] = -1; out_b[i] = 0.0; out_n[i] = 1; out_a[i] = 0.0 if g < 0.0: g = 0.0 if j != j_gift: g = 0.0 out_choice[i] = ch out_c[i] = c out_s[i] = s out_g[i] = g out_pay[i] = pay_out out_Bnew[i] = Bnew out_flag[i] = flag out_u[i] = _u(max(c, 1e-9), max(s, 1e-9), alpha, gamma) out_x[i] = x_used if j == J: hh3 = h_grid[out_h[i]] if out_h[i] >= 0 else 0.0 gross = out_a[i] + (1.0 - kappa_s) * p_next * hh3 - out_b[i] if gross < 0.0: gross = 0.0 out_estate[i] = gross else: out_estate[i] = 0.0 @njit(cache=True) def _draw_z(z, age, u, cumPi, j_retire): """In-place Markov draw for ages 2..j_retire-1 (frozen afterwards; age 1 drawn at birth).""" M = z.shape[0] nz = cumPi.shape[1] for i in range(M): j = age[i] if j >= 2 and j < j_retire: r = u[i] zi = z[i] k = 0 while k < nz - 1 and r > cumPi[zi, k]: k += 1 z[i] = k @njit(cache=True) def _draw_z1(z_parent, u, cumPg): M = z_parent.shape[0] nz = cumPg.shape[1] out = np.empty(M, dtype=np.int64) for i in range(M): r = u[i] zi = z_parent[i] k = 0 while k < nz - 1 and r > cumPg[zi, k]: k += 1 out[i] = k return out # -------------------------------------------------------------------------------------- # Panel state # -------------------------------------------------------------------------------------- @dataclass class Panel: """Dynasty slots: young (ages 1-6) and old (ages 7-12) members, start-of-period states.""" phase: np.ndarray # 0..5 age_y: np.ndarray; z_y: np.ndarray; h_y: np.ndarray; b_y: np.ndarray; n_y: np.ndarray; a_y: np.ndarray age_o: np.ndarray; z_o: np.ndarray; h_o: np.ndarray; b_o: np.ndarray; n_o: np.ndarray; a_o: np.ndarray estate_pending: np.ndarray # estate to be received by the young member at age 7 gift_received: np.ndarray # gift (net) received at birth by the current young member parent_owner_at_gift: np.ndarray # parent tenure (1 owner) at the time of the gift (for mobility) parent_z_at_gift: np.ndarray parent_nw_at_gift: np.ndarray estate_received: np.ndarray # estate (net) received by the current old member at age 7 parent_owner_old: np.ndarray = None # parent tenure at gift, carried after promotion parent_nw_old: np.ndarray = None parent_z_old: np.ndarray = None gift_received_old: np.ndarray = None t: int = 0 @classmethod def initial(cls, P: Params, rng: np.random.Generator) -> "Panel": M = P.n_dyn phase = np.arange(M) % 6 z_y = rng.choice(P.n_z, size=M, p=P.pi_z) z_o = rng.choice(P.n_z, size=M, p=P.pi_z) zeros = np.zeros(M) return cls(phase=phase, 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(), 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(), estate_pending=zeros.copy(), gift_received=zeros.copy(), parent_owner_at_gift=np.zeros(M, dtype=np.int64), parent_z_at_gift=z_o.copy(), parent_nw_at_gift=zeros.copy(), estate_received=zeros.copy(), parent_owner_old=np.zeros(M, dtype=np.int64), parent_nw_old=zeros.copy(), parent_z_old=z_o.copy(), gift_received_old=zeros.copy(), t=0) def copy(self) -> "Panel": kw = {k: (v.copy() if isinstance(v, np.ndarray) else v) for k, v in self.__dict__.items()} return Panel(**kw) @dataclass class PeriodRecord: """Decisions of one period for all households (young + old concatenated).""" age: np.ndarray; z: np.ndarray; h0: np.ndarray; b0: np.ndarray; a0: np.ndarray; n0: np.ndarray choice: np.ndarray; c: np.ndarray; s: np.ndarray; h1: np.ndarray; b1: np.ndarray; n1: np.ndarray; a1: np.ndarray gift: np.ndarray; pay: np.ndarray; Bnew: np.ndarray; estate: np.ndarray; flag: np.ndarray; u: np.ndarray; x: np.ndarray y_shock: np.ndarray gift_received: np.ndarray; parent_owner: np.ndarray; parent_nw: np.ndarray; estate_received: np.ndarray slot: np.ndarray; is_old: np.ndarray p: float; p_r: float; rho: float def owner0(self) -> np.ndarray: return self.h0 >= 0 def owner1(self) -> np.ndarray: return self.h1 >= 0 def net_worth0(self, P: Params) -> np.ndarray: hv = np.where(self.h0 >= 0, P.h_grid[np.maximum(self.h0, 0)], 0.0) return self.a0 + (1.0 - P.kappa_s) * self.p * hv - self.b0 def net_worth1(self, P: Params) -> np.ndarray: hv = np.where(self.h1 >= 0, P.h_grid[np.maximum(self.h1, 0)], 0.0) return self.a1 + (1.0 - P.kappa_s) * self.p * hv - self.b1 def housing_demand(self, P: Params) -> float: """Owner-occupied units demanded per unit mass of households of one age.""" hv = np.where(self.h1 >= 0, P.h_grid[np.maximum(self.h1, 0)], 0.0) return float(hv.sum()) / (len(self.age) / P.J) def rental_demand(self, P: Params) -> float: return float(np.where(self.h1 < 0, self.s, 0.0).sum()) / (len(self.age) / P.J) def _kernel_args(P: Params, S: PolicyStack, pr: Prices): 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) 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, 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), bool(pr.reset), pr.p, (pr.p_next if pr.p_next is not None else pr.p), pr.rho, pr.rho * pr.rent_rigid_factor, pr.R_a, pr.i_new, pr.pi_next, pr.transfer_rebate, pr.wage_factor, 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, P.N, P.j_gift, P.J) def _run_group(P, S, pr, y_base, age, z, h, b, n, a, y_noise=None): M = len(age) if y_noise is None: y_noise = np.ones(M) out = dict(choice=np.empty(M, dtype=np.int64), c=np.empty(M), s=np.empty(M), h=np.empty(M, dtype=np.int64), b=np.empty(M), n=np.empty(M, dtype=np.int64), a=np.empty(M), g=np.empty(M), pay=np.empty(M), Bnew=np.empty(M), estate=np.empty(M), flag=np.empty(M, dtype=np.int64), u=np.empty(M), x=np.empty(M)) _sim_step(age, z, h, b, n, a, y_base, np.ascontiguousarray(y_noise, dtype=float), *_kernel_args(P, S, pr), out["choice"], out["c"], out["s"], out["h"], out["b"], out["n"], out["a"], out["g"], out["pay"], out["Bnew"], out["estate"], out["flag"], out["u"], out["x"]) return out def simulate_period(P: Params, S: PolicyStack, pr: Prices, panel: Panel, rng_key: tuple, balance_factor: float = 1.0, record: bool = True) -> PeriodRecord | None: """Advance the panel by one period under policies S and prices pr; returns the record. balance_factor multiplies all real mortgage balances at the start of the period (surprise inflation erosion 1/(1+pi_u), undone for the ARM contract); deposits are revalued through pr.R_a. rng_key seeds the period's draws (common random numbers). """ M = P.n_dyn rng = np.random.default_rng(list(rng_key)) U = rng.uniform(size=(3, M)) if P.sigma_trans > 0: Nz = rng.standard_normal(size=(2, M)) noise_o = np.exp(P.sigma_trans * Nz[0] - 0.5 * P.sigma_trans**2) noise_y = np.exp(P.sigma_trans * Nz[1] - 0.5 * P.sigma_trans**2) else: noise_o = np.ones(M); noise_y = np.ones(M) y_base = P.income_matrix() cumPi = np.cumsum(P.Pi_z, axis=1) cumPg = np.cumsum(P.Pi_g, axis=1) if balance_factor != 1.0: panel.b_o *= balance_factor panel.b_y *= balance_factor # income shocks (ages >= 2) _draw_z(panel.z_o, panel.age_o, U[0], cumPi, P.j_retire) _draw_z(panel.z_y, panel.age_y, U[1], cumPi, P.j_retire) # --- old members decide first 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) # newborns (young at age 1) receive the gift of the old member at age 7 (same slot) newborn = panel.age_y == 1 if newborn.any(): gift_net = (1.0 - P.tau_b) * o["g"] panel.gift_received[newborn] = gift_net[newborn] panel.a_y[newborn] = gift_net[newborn] + pr.transfer_rebate 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_o panel.parent_owner_at_gift[newborn] = (o["h"][newborn] >= 0).astype(np.int64) panel.parent_z_at_gift[newborn] = panel.z_o[newborn] panel.parent_nw_at_gift[newborn] = nw_parent[newborn] panel.z_y[newborn] = _draw_z1(panel.z_o, U[2], cumPg)[newborn] panel.h_y[newborn] = -1; panel.b_y[newborn] = 0.0; panel.n_y[newborn] = 1 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) rec = None if not record else PeriodRecord( age=np.concatenate([panel.age_y, panel.age_o]), z=np.concatenate([panel.z_y, panel.z_o]), h0=np.concatenate([panel.h_y, panel.h_o]), b0=np.concatenate([panel.b_y, panel.b_o]), a0=np.concatenate([panel.a_y, panel.a_o]), n0=np.concatenate([panel.n_y, panel.n_o]), choice=np.concatenate([yg["choice"], o["choice"]]), c=np.concatenate([yg["c"], o["c"]]), s=np.concatenate([yg["s"], o["s"]]), h1=np.concatenate([yg["h"], o["h"]]), b1=np.concatenate([yg["b"], o["b"]]), n1=np.concatenate([yg["n"], o["n"]]), a1=np.concatenate([yg["a"], o["a"]]), gift=np.concatenate([yg["g"], o["g"]]), pay=np.concatenate([yg["pay"], o["pay"]]), Bnew=np.concatenate([yg["Bnew"], o["Bnew"]]), estate=np.concatenate([yg["estate"], o["estate"]]), flag=np.concatenate([yg["flag"], o["flag"]]), u=np.concatenate([yg["u"], o["u"]]), x=np.concatenate([yg["x"], o["x"]]), y_shock=np.concatenate([noise_y, noise_o]), gift_received=np.concatenate([panel.gift_received, panel.gift_received_old]), parent_owner=np.concatenate([panel.parent_owner_at_gift, panel.parent_owner_old]), parent_nw=np.concatenate([panel.parent_nw_at_gift, panel.parent_nw_old]), estate_received=np.concatenate([np.full(M, np.nan), panel.estate_received]), slot=np.concatenate([np.arange(M), np.arange(M)]), is_old=np.concatenate([np.zeros(M, bool), np.ones(M, bool)]), p=pr.p, p_r=pr.p_r, rho=pr.rho) # --- advance states panel.h_y, panel.b_y, panel.n_y, panel.a_y = yg["h"], yg["b"], yg["n"], yg["a"] panel.h_o, panel.b_o, panel.n_o, panel.a_o = o["h"], o["b"], o["n"], o["a"] dying = panel.age_o == P.J panel.estate_pending[dying] = (1.0 - P.tau_b) * o["estate"][dying] panel.age_y += 1 panel.age_o += 1 # young at 7 -> old; estate received now (start of next period); new young born next period promote = panel.age_y == P.j_estate if promote.any(): panel.age_o[promote] = P.j_estate panel.z_o[promote] = panel.z_y[promote] panel.h_o[promote] = panel.h_y[promote]; panel.b_o[promote] = panel.b_y[promote] panel.n_o[promote] = panel.n_y[promote] panel.a_o[promote] = panel.a_y[promote] + panel.estate_pending[promote] panel.estate_received[promote] = panel.estate_pending[promote] panel.estate_pending[promote] = 0.0 panel.parent_owner_old[promote] = panel.parent_owner_at_gift[promote] panel.parent_nw_old[promote] = panel.parent_nw_at_gift[promote] panel.parent_z_old[promote] = panel.parent_z_at_gift[promote] panel.gift_received_old[promote] = panel.gift_received[promote] panel.age_y[promote] = 1 panel.h_y[promote] = -1; panel.b_y[promote] = 0.0; panel.n_y[promote] = 1; panel.a_y[promote] = 0.0 panel.t += 1 return rec # -------------------------------------------------------------------------------------- # Stationary simulation and aggregation # -------------------------------------------------------------------------------------- def simulate_stationary(P: Params, pols: list[AgePolicy], pr: Prices, n_record: int = 6, panel0: Panel | None = None, burn: int | None = None) -> tuple[list[PeriodRecord], Panel]: """Simulate `burn` periods from an arbitrary initial condition under stationary policies, then record `n_record` further periods (one full phase cycle = 6). Returns the records and the panel state at the end (start-of-period states for a subsequent transition).""" S = stack_policies(pols) burn = P.burn if burn is None else burn rng = np.random.default_rng(P.seed) panel = Panel.initial(P, rng) if panel0 is None else panel0.copy() recs: list[PeriodRecord] = [] for t in range(burn + n_record): rec = simulate_period(P, S, pr, panel, (P.seed, t), record=(t >= burn)) if t >= burn: recs.append(rec) return recs, panel def pool_records(recs: list[PeriodRecord]) -> PeriodRecord: kw = {} for k, v in recs[0].__dict__.items(): if isinstance(v, np.ndarray): kw[k] = np.concatenate([getattr(r, k) for r in recs]) else: kw[k] = v return PeriodRecord(**kw) def estate_lottery(P: Params, rec: PeriodRecord) -> tuple[np.ndarray, np.ndarray]: """K-point lottery: conditional means of net estates received within K equal-probability bins.""" e = rec.estate_received[rec.is_old & (rec.age == P.j_estate)] e = e[np.isfinite(e)] if len(e) == 0: return np.zeros(P.K_estate), np.full(P.K_estate, 1.0 / P.K_estate) e = np.sort(e) bins = np.array_split(e, P.K_estate) pts = np.array([b.mean() if len(b) else 0.0 for b in bins]) return pts, np.full(P.K_estate, 1.0 / P.K_estate) def gini(x: np.ndarray) -> float: x = np.asarray(x, dtype=float) x = x[np.isfinite(x)] if len(x) == 0: return np.nan x = np.sort(x - min(x.min(), 0.0)) n = len(x) if x.sum() <= 0: return np.nan return float((2.0 * np.sum(np.arange(1, n + 1) * x) / (n * x.sum())) - (n + 1.0) / n) def rank_rank_slope(parent: np.ndarray, child: np.ndarray) -> float: ok = np.isfinite(parent) & np.isfinite(child) if ok.sum() < 10: return np.nan rp = np.argsort(np.argsort(parent[ok])) / ok.sum() rc = np.argsort(np.argsort(child[ok])) / ok.sum() return float(np.polyfit(rp, rc, 1)[0]) # -------------------------------------------------------------------------------------- # Cohort simulations (welfare decomposition) and deterministic paths (experiment 2) # -------------------------------------------------------------------------------------- def simulate_cohort(P: Params, S: PolicyStack, pr: Prices, z1: np.ndarray, a1: np.ndarray, lottery, seed: int = 11, fixed_z: bool = False) -> dict: """Simulate a cohort from age 1 to J (no dynasty link; estate drawn from `lottery` at j_estate). Returns discounted flow utility U, warm-glow part Phi, and per-age arrays. fixed_z=True keeps z constant (equal-income deterministic design).""" n = len(a1) rng = np.random.default_rng(seed) y_base = P.income_matrix() cumPi = np.cumsum(P.Pi_z, axis=1) z = z1.astype(np.int64).copy(); a = a1.astype(float).copy() h = np.full(n, -1, dtype=np.int64); b = np.zeros(n); nn = np.ones(n, dtype=np.int64) U = np.zeros(n); Phi = np.zeros(n) 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)), 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)) pts, w = lottery for j in range(1, P.J + 1): age = np.full(n, j, dtype=np.int64) if j >= 2 and not fixed_z: _draw_z(z, age, rng.uniform(size=n), cumPi, P.j_retire) if j == P.j_estate: k = rng.choice(len(pts), size=n, p=w) a = a + pts[k] out["estate"] = pts[k] o = _run_group(P, S, pr, y_base, age, z, h, b, nn, a) disc = P.beta ** (j - 1) U += disc * o["u"] if j == P.j_gift: Phi += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["g"])) out["gift"] = o["g"] if j == P.J: Phi += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["estate"])) hv = np.where(o["h"] >= 0, P.h_grid[np.maximum(o["h"], 0)], 0.0) out["own"][j - 1] = o["h"] >= 0 out["nw"][j - 1] = o["a"] + (1 - P.kappa_s) * pr.p * hv - o["b"] out["a"][j - 1] = o["a"]; out["b"][j - 1] = o["b"]; out["c"][j - 1] = o["c"]; out["s"][j - 1] = o["s"] out["h"][j - 1] = o["h"] h, b, nn, a = o["h"], o["b"], o["n"], o["a"] out["U"] = U; out["Phi"] = Phi; out["V"] = U + Phi return out def _wg(P: Params, x): from .model import warm_glow return warm_glow(P, np.asarray(x, dtype=float)) def remaining_lifetime(P: Params, S: PolicyStack, pr: Prices, rec: PeriodRecord, lottery, seed: int = 12, start_states: dict | None = None) -> tuple[np.ndarray, np.ndarray]: """Discounted flow utility U and warm-glow Phi over the remaining lifetime of every household in `rec` (start-of-period states, before this period's decision), under policies S and prices pr (stationary). Used for the welfare decomposition of households alive at a shock.""" n = len(rec.age) rng = np.random.default_rng(seed) y_base = P.income_matrix() cumPi = np.cumsum(P.Pi_z, axis=1) st = start_states or {} age = rec.age.copy(); z = rec.z.copy(); h = st.get("h", rec.h0).copy(); b = st.get("b", rec.b0).copy() nn = rec.n0.copy(); a = st.get("a", rec.a0).copy() U = np.zeros(n); Phi = np.zeros(n) pts, w = lottery alive = np.ones(n, dtype=bool) step = 0 while alive.any(): idx = np.where(alive)[0] if step > 0: _draw_z(z, age, rng.uniform(size=n), cumPi, P.j_retire) est = age == P.j_estate if est.any(): k = rng.choice(len(pts), size=int(est.sum()), p=w) a[est] += pts[k] o = _run_group(P, S, pr, y_base, age[idx], z[idx], h[idx], b[idx], nn[idx], a[idx]) disc = P.beta ** step U[idx] += disc * o["u"] g_age = age[idx] == P.j_gift if g_age.any(): Phi[idx[g_age]] += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["g"][g_age])) last = age[idx] == P.J if last.any(): Phi[idx[last]] += disc * np.asarray(_wg(P, (1 - P.tau_b) * o["estate"][last])) h[idx], b[idx], nn[idx], a[idx] = o["h"], o["b"], o["n"], o["a"] age[idx] += 1 alive = age <= P.J step += 1 return U, Phi