TeX 53.2%
Python 46.6%
1"""WP18 — household dynamic programme (Part B), solved by backward induction over age.23Method (justification in research/results_audit.md §algorithm): discrete-choice value function4iteration with exhaustive search over the liquid-asset choice on the asset grid. The discrete5choices (stay / refinance / move / sell-and-rent / rent / buy, house size, loan size, gift) make6the problem non-concave; exhaustive search is robust to the kinks and, compiled with numba, is7cheap enough (one age solves in ~50 ms on the baseline grids). EGM with upper envelopes8(Iskhakov et al. 2017) would require taste shocks; we avoid them and report Bellman residuals9at off-grid points instead (`bellman_residuals`).1011State: (j, z, h, b, n, a) with h the house index (-1 = renter), b the real balance at the start12of the period, n the remaining payments, a the liquid real asset. All branches except "stay"13depend on the state only through cash-on-hand x, so they are solved on a cash grid:14 W_R(z, x) rent W_B(z, x) buy (h', B', a') W_F(z, h, x) refinance15 W_S(z, h, b, n, a) stay (continuous a').16At the gift age the gift is chosen on top of each branch (gift transform). The estate at the17estate age is anticipated as a K-point lottery (equilibrium object).1819Infeasible states carry the sentinel NEG = -1e10.20"""21from __future__ import annotations2223from dataclasses import dataclass2425import numpy as np26from numba import njit, prange2728from .model import (Params, Prices, annuity_factor, next_balance, warm_glow)2930NEG = -1.0e10313233# --------------------------------------------------------------------------------------34# Small numba helpers35# --------------------------------------------------------------------------------------36@njit(cache=True, inline="always")37def _locate(grid, x):38 """Index i and weight w such that x ~ (1-w) grid[i] + w grid[i+1]; clamped at the ends."""39 n = grid.shape[0]40 if x <= grid[0]:41 return 0, 0.042 if x >= grid[n - 1]:43 return n - 2, 1.044 lo = 045 hi = n - 146 while hi - lo > 1:47 mid = (lo + hi) // 248 if grid[mid] <= x:49 lo = mid50 else:51 hi = mid52 w = (x - grid[lo]) / (grid[lo + 1] - grid[lo])53 return lo, w545556@njit(cache=True, inline="always")57def _lin(v0, v1, w):58 if v0 <= NEG * 0.5 or v1 <= NEG * 0.5:59 return NEG60 return (1.0 - w) * v0 + w * v1616263@njit(cache=True, inline="always")64def _u(c, s, alpha, gamma):65 comp = c ** (1.0 - alpha) * s ** alpha66 return comp ** (1.0 - gamma) / (1.0 - gamma)676869@njit(cache=True, inline="always")70def _renter_cs(e, rho, alpha, s_min, s_max):71 """Intratemporal renter choice given expenditure e = c + rho s."""72 s = alpha * e / rho73 if s < s_min:74 s = s_min75 elif s > s_max:76 s = s_max77 c = e - rho * s78 return c, s798081# --------------------------------------------------------------------------------------82# Branch solvers83# --------------------------------------------------------------------------------------84GOLD = 0.618033988749894985N_GOLD = 18868788@njit(cache=True, inline="always")89def _af_due(i, n):90 """Annuity-due factor: payment per unit of balance, n payments, first one now."""91 if n <= 0:92 return 0.093 if abs(i) < 1e-12:94 return 1.0 / n95 return i / ((1.0 + i) * (1.0 - (1.0 + i) ** (-n)))969798@njit(cache=True, inline="always")99def _ev_interp(EV2, jb, wb, a_grid, ap):100 """EV2 is a (nb, na) view; interpolate in b (jb, wb) and in a (ap)."""101 ja, wa = _locate(a_grid, ap)102 v0 = _lin(EV2[jb, ja], EV2[jb, ja + 1], wa)103 v1 = _lin(EV2[jb + 1, ja], EV2[jb + 1, ja + 1], wa)104 return _lin(v0, v1, wb)105106107@njit(cache=True, inline="always")108def _obj(ap, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max):109 """Period objective of saving ap. mode 0: owner (s given); mode 1: renter (rho given)."""110 e = cash - ap111 if mode == 0:112 c = e113 ss = s114 else:115 c, ss = _renter_cs(e, rho, alpha, s_min, s_max)116 if c <= c_min:117 return NEG118 ev = _ev_interp(EV2, jb, wb, a_grid, ap)119 if ev <= NEG * 0.5:120 return NEG121 return _u(c, ss, alpha, gamma) + beta * ev122123124@njit(cache=True)125def _best_a(cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max):126 """Grid search over a' on a_grid, then golden-section refinement between the neighbours."""127 na = a_grid.shape[0]128 best = NEG129 kbest = -1130 for ka in range(na):131 if cash - a_grid[ka] <= c_min:132 break133 v = _obj(a_grid[ka], cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max)134 if v > best:135 best = v136 kbest = ka137 if kbest < 0:138 return NEG, 0.0139 lo = a_grid[kbest - 1] if kbest > 0 else a_grid[0]140 hi = a_grid[kbest + 1] if kbest < na - 1 else a_grid[na - 1]141 if cash - c_min < hi:142 hi = cash - c_min143 if hi <= lo:144 return best, a_grid[kbest]145 x1 = hi - GOLD * (hi - lo)146 x2 = lo + GOLD * (hi - lo)147 f1 = _obj(x1, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max)148 f2 = _obj(x2, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max)149 for _ in range(N_GOLD):150 if f1 < f2:151 lo = x1152 x1 = x2153 f1 = f2154 x2 = lo + GOLD * (hi - lo)155 f2 = _obj(x2, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max)156 else:157 hi = x2158 x2 = x1159 f2 = f1160 x1 = hi - GOLD * (hi - lo)161 f1 = _obj(x1, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max)162 abest = a_grid[kbest]163 if f1 > best:164 best = f1165 abest = x1166 if f2 > best:167 best = f2168 abest = x2169 return best, abest170171172@njit(cache=True, parallel=True)173def _solve_stay(y, a_grid, b_grid, h_grid, ic_by_n, pi_next, R_a, p, m, chi, alpha, gamma, beta,174 c_min, EV_own, W_S, aS, cash_S):175 """Owner stays: c + pay + m p h + a' = y + a R_a. EV_own indexed by (z, h, b', n', a')."""176 nz, nh, nb, nN, na = EV_own.shape177 for iz in prange(nz):178 for ih in range(nh):179 h = h_grid[ih]180 s = (1.0 + chi) * h181 for ib in range(nb):182 b = b_grid[ib]183 for n_idx in range(nN):184 n = n_idx + 1185 ic = ic_by_n[n_idx]186 pay = b * _af_due(ic, n)187 bn = b - pay188 if bn < 0.0:189 bn = 0.0190 bn = bn * (1.0 + ic) / (1.0 + pi_next)191 jb, wb = _locate(b_grid, bn)192 n_next = n_idx - 1 if n_idx > 0 else 0193 EV2 = EV_own[iz, ih, :, n_next, :]194 for ia in range(na):195 x = y[iz] + a_grid[ia] * R_a - pay - m * p * h196 cash_S[iz, ih, ib, n_idx, ia] = x197 best, abest = _best_a(x, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0)198 W_S[iz, ih, ib, n_idx, ia] = best199 aS[iz, ih, ib, n_idx, ia] = abest200201202@njit(cache=True, parallel=True)203def _solve_rent(x_grid, a_grid, rho, alpha, gamma, beta, s_min, s_max, c_min, EV_rent, W_R, aR):204 nz = EV_rent.shape[0]205 nx = x_grid.shape[0]206 for iz in prange(nz):207 EV2 = EV_rent[iz:iz + 1, :]208 for ix in range(nx):209 best, abest = _best_a(x_grid[ix], 0.0, rho, 1, EV2, 0, 0.0, a_grid, alpha, gamma, beta, c_min, s_min, s_max)210 W_R[iz, ix] = best211 aR[iz, ix] = abest212213214@njit(cache=True, parallel=True)215def _solve_buy(x_grid, a_grid, b_grid, h_grid, B_fracs, maxloan, i_new, pi_next, p, kappa_b, m,216 chi, alpha, gamma, beta, N, c_min, EV_own, W_Bh, BBh, aBh):217 """Renter buys / owner moves: c + (1+kb+m) p h' - B' + pay(B') + a' = x. maxloan[z, h'].218 Outputs are per house size h' (W_Bh[z, x, h'] = best over B', a') so that the size choice can219 be interpolated in cash at simulation time (avoids synchronised jumps of the demand)."""220 nz = EV_own.shape[0]221 nh = h_grid.shape[0]222 nx = x_grid.shape[0]223 nB = B_fracs.shape[0]224 af = _af_due(i_new, N)225 for iz in prange(nz):226 for ix in range(nx):227 x = x_grid[ix]228 for ih in range(nh):229 h = h_grid[ih]230 s = (1.0 + chi) * h231 outlay = (1.0 + kappa_b + m) * p * h232 EV2 = EV_own[iz, ih, :, N - 2, :]233 best = NEG234 bestB = 0.0235 besta = 0.0236 for kB in range(nB):237 B = B_fracs[kB] * maxloan[iz, ih]238 pay = B * af239 bn = (B - pay) * (1.0 + i_new) / (1.0 + pi_next)240 jb, wb = _locate(b_grid, bn)241 cash = x - outlay + B - pay242 if cash <= c_min:243 continue244 v, ab = _best_a(cash, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0)245 if v > best:246 best = v247 bestB = B248 besta = ab249 W_Bh[iz, ix, ih] = best250 BBh[iz, ix, ih] = bestB251 aBh[iz, ix, ih] = besta252253254@njit(cache=True, parallel=True)255def _solve_refi(x_grid, a_grid, b_grid, h_grid, B_fracs, maxloan, i_new, pi_next, p, m, chi,256 alpha, gamma, beta, N, c_min, EV_own, W_F, BF, aF):257 """Owner refinances (keeps h): c + a' = x + B' - pay(B'), x = y + a R_a - b - (kf + m) p h."""258 nz = EV_own.shape[0]259 nh = h_grid.shape[0]260 nx = x_grid.shape[0]261 nB = B_fracs.shape[0]262 af = _af_due(i_new, N)263 for iz in prange(nz):264 for ih in range(nh):265 h = h_grid[ih]266 s = (1.0 + chi) * h267 EV2 = EV_own[iz, ih, :, N - 2, :]268 for ix in range(nx):269 x = x_grid[ix]270 best = NEG271 bestB = 0.0272 besta = 0.0273 for kB in range(nB):274 B = B_fracs[kB] * maxloan[iz, ih]275 pay = B * af276 bn = (B - pay) * (1.0 + i_new) / (1.0 + pi_next)277 jb, wb = _locate(b_grid, bn)278 cash = x + B - pay279 if cash <= c_min:280 continue281 v, ab = _best_a(cash, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0)282 if v > best:283 best = v284 bestB = B285 besta = ab286 W_F[iz, ih, ix] = best287 BF[iz, ih, ix] = bestB288 aF[iz, ih, ix] = besta289290291@njit(cache=True)292def _gift_transform_x(W, x_grid, gift_fracs, nu, x_shift, gamma, tau_b, Wg, Gg):293 """Wg(z[,h], x) = max_f phi((1-tau) f x) + W(z[,h], (1-f) x); W shaped (..., nx)."""294 lead = W.shape[0]295 nx = x_grid.shape[0]296 nf = gift_fracs.shape[0]297 for il in range(lead):298 for ix in range(nx):299 x = x_grid[ix]300 best = NEG301 bestg = 0.0302 for kf in range(nf):303 g = gift_fracs[kf] * x304 net = (1.0 - tau_b) * g305 phi = nu * (net + x_shift) ** (1.0 - gamma) / (1.0 - gamma) if nu > 0.0 else 0.0306 xr = x - g307 j, w = _locate(x_grid, xr)308 wv = _lin(W[il, j], W[il, j + 1], w)309 if wv <= NEG * 0.5:310 continue311 v = phi + wv312 if v > best:313 best = v314 bestg = g315 Wg[il, ix] = best316 Gg[il, ix] = bestg317318319@njit(cache=True, parallel=True)320def _gift_transform_stay(W_S, cash_S, a_grid, R_a, gift_fracs, nu, x_shift, gamma, tau_b, Wg, Gg):321 """Stay branch: giving g lowers cash by g, equivalent to a lower asset a - g/R_a."""322 nz, nh, nb, nN, na = W_S.shape323 nf = gift_fracs.shape[0]324 for iz in prange(nz):325 for ih in range(nh):326 for ib in range(nb):327 for n_idx in range(nN):328 for ia in range(na):329 x = cash_S[iz, ih, ib, n_idx, ia]330 best = NEG331 bestg = 0.0332 for kf in range(nf):333 g = gift_fracs[kf] * x if x > 0.0 else 0.0334 if g < 0.0:335 g = 0.0336 net = (1.0 - tau_b) * g337 phi = nu * (net + x_shift) ** (1.0 - gamma) / (1.0 - gamma) if nu > 0.0 else 0.0338 ar = a_grid[ia] - g / R_a339 if ar < 0.0:340 continue341 j, w = _locate(a_grid, ar)342 wv = _lin(W_S[iz, ih, ib, n_idx, j], W_S[iz, ih, ib, n_idx, j + 1], w)343 if wv <= NEG * 0.5:344 continue345 v = phi + wv346 if v > best:347 best = v348 bestg = g349 if kf == 0 and nu == 0.0:350 break351 Wg[iz, ih, ib, n_idx, ia] = best352 Gg[iz, ih, ib, n_idx, ia] = bestg353354355@njit(cache=True, parallel=True)356def _combine_owner(y, a_grid, b_grid, h_grid, x_grid, R_a, p, kappa_s, kappa_f, m,357 W_S, W_F, W_B, W_R, V_own, choice_own):358 """V_own = max(stay, refi(x_refi), buy(x_sell), rent(x_sell)); choice: 0 S, 1 F, 2 B, 3 R, -1 none."""359 nz, nh, nb, nN, na = W_S.shape360 for iz in prange(nz):361 for ih in range(nh):362 h = h_grid[ih]363 for ib in range(nb):364 b = b_grid[ib]365 for n_idx in range(nN):366 for ia in range(na):367 base = y[iz] + a_grid[ia] * R_a368 x_sell = base + (1.0 - kappa_s) * p * h - b369 x_refi = base - b - (kappa_f + m) * p * h370 best = W_S[iz, ih, ib, n_idx, ia]371 ch = 0 if best > NEG * 0.5 else -1372 if x_refi > 0.0:373 j, w = _locate(x_grid, x_refi)374 v = _lin(W_F[iz, ih, j], W_F[iz, ih, j + 1], w)375 if v > best:376 best = v377 ch = 1378 if x_sell > 0.0:379 j, w = _locate(x_grid, x_sell)380 v = _lin(W_B[iz, j], W_B[iz, j + 1], w)381 if v > best:382 best = v383 ch = 2384 v = _lin(W_R[iz, j], W_R[iz, j + 1], w)385 if v > best:386 best = v387 ch = 3388 V_own[iz, ih, ib, n_idx, ia] = best389 choice_own[iz, ih, ib, n_idx, ia] = ch390391392@njit(cache=True)393def _combine_renter(y, a_grid, x_grid, R_a, rebate, W_B, W_R_sit, V_rent, choice_rent):394 nz = W_B.shape[0]395 na = a_grid.shape[0]396 for iz in range(nz):397 for ia in range(na):398 x = y[iz] + a_grid[ia] * R_a + rebate399 j, w = _locate(x_grid, x)400 vb = _lin(W_B[iz, j], W_B[iz, j + 1], w)401 vr = _lin(W_R_sit[iz, j], W_R_sit[iz, j + 1], w)402 if vr >= vb:403 V_rent[iz, ia] = vr404 choice_rent[iz, ia] = 3 if vr > NEG * 0.5 else -1405 else:406 V_rent[iz, ia] = vb407 choice_rent[iz, ia] = 2408409410@njit(cache=True)411def _expect_lottery(V, a_grid, lot_pts, lot_w, out):412 """out[..., a] = sum_k w_k V[..., a + e_k] (linear interpolation / flat extrapolation in a)."""413 shp = V.shape414 na = shp[-1]415 V2 = V.reshape((-1, na))416 o2 = out.reshape((-1, na))417 K = lot_pts.shape[0]418 for r in range(V2.shape[0]):419 for ia in range(na):420 acc = 0.0421 bad = False422 for k in range(K):423 j, w = _locate(a_grid, a_grid[ia] + lot_pts[k])424 v = _lin(V2[r, j], V2[r, j + 1], w)425 if v <= NEG * 0.5:426 bad = True427 break428 acc += lot_w[k] * v429 o2[r, ia] = NEG if bad else acc430431432# --------------------------------------------------------------------------------------433# Python-level age solve434# --------------------------------------------------------------------------------------435@dataclass436class AgePolicy:437 """Branch values and policies for one age (and one date along a transition)."""438 W_S: np.ndarray; aS: np.ndarray; gS: np.ndarray439 W_F: np.ndarray; BF: np.ndarray; aF: np.ndarray; gF: np.ndarray440 W_B: np.ndarray; hB: np.ndarray; BB: np.ndarray; aB: np.ndarray; gB: np.ndarray # best over h' (reporting)441 W_Bh: np.ndarray; BBh: np.ndarray; aBh: np.ndarray; gBh: np.ndarray # per house size (simulation)442 W_R: np.ndarray; aR: np.ndarray; gR: np.ndarray # new-tenant rent branch443 W_Rs: np.ndarray; aRs: np.ndarray; gRs: np.ndarray # sitting-tenant rent branch444 V_own: np.ndarray; V_rent: np.ndarray445 choice_own: np.ndarray; choice_rent: np.ndarray446 maxloan: np.ndarray447 rho_sit: float448449 def to_f32(self) -> "AgePolicy":450 kw = {}451 for k, v in self.__dict__.items():452 kw[k] = v.astype(np.float32) if isinstance(v, np.ndarray) and v.dtype == np.float64 else v453 return AgePolicy(**kw)454455456def expectation(P: Params, j: int, V_own_next: np.ndarray, V_rent_next: np.ndarray,457 lottery: tuple[np.ndarray, np.ndarray] | None) -> tuple[np.ndarray, np.ndarray]:458 """E[V_{j+1} | z] for a household of age j (1-based): Markov in z (frozen in retirement),459 estate lottery if j + 1 == j_estate."""460 if lottery is not None and j + 1 == P.j_estate:461 pts, w = lottery462 Vo = np.empty_like(V_own_next); Vr = np.empty_like(V_rent_next)463 _expect_lottery(V_own_next, P.a_grid, pts, w, Vo)464 _expect_lottery(V_rent_next, P.a_grid, pts, w, Vr)465 else:466 Vo, Vr = V_own_next, V_rent_next467 if j + 1 >= P.j_retire:468 return Vo.copy(), Vr.copy()469 Pi = P.Pi_z470 bad_o = Vo <= NEG * 0.5471 bad_r = Vr <= NEG * 0.5472 EVo = np.tensordot(Pi, np.where(bad_o, 0.0, Vo), axes=(1, 0))473 EVr = Pi @ np.where(bad_r, 0.0, Vr)474 # a state is infeasible in expectation if any reachable successor is infeasible475 EVo = np.where(np.tensordot(Pi, bad_o.astype(float), axes=(1, 0)) > 0, NEG, EVo)476 EVr = np.where((Pi @ bad_r.astype(float)) > 0, NEG, EVr)477 return EVo, EVr478479480def terminal_continuation(P: Params, pr: Prices) -> tuple[np.ndarray, np.ndarray]:481 """Warm glow of the (net) estate as a function of the next-period state (h, b', a')."""482 p_next = pr.p_next if pr.p_next is not None else pr.p483 a = P.a_grid[None, None, None, None, :]484 b = P.b_grid[None, None, :, None, None]485 h = P.h_grid[None, :, None, None, None]486 estate_own = (1.0 - P.tau_b) * (a + (1.0 - P.kappa_s) * p_next * h - b)487 estate_own = np.broadcast_to(estate_own, (P.n_z, P.n_h, P.n_b, P.N, P.n_a))488 EV_own = np.where(estate_own >= 0.0, warm_glow(P, np.maximum(estate_own, 0.0)), NEG)489 EV_rent = np.broadcast_to(warm_glow(P, (1.0 - P.tau_b) * P.a_grid)[None, :], (P.n_z, P.n_a))490 return np.ascontiguousarray(EV_own, dtype=float), np.ascontiguousarray(EV_rent, dtype=float)491492493def max_loan_matrix(P: Params, pr: Prices, y: np.ndarray) -> np.ndarray:494 """maxloan[z, h'] = min(LTV cap, DSTI cap at the qualifying rate)."""495 af_q = annuity_factor(pr.i_q, P.N)496 ltv = (1.0 - P.theta) * pr.p * P.h_grid[None, :]497 dsti = (P.psi * y / af_q)[:, None] if af_q > 0 else np.full((P.n_z, 1), np.inf)498 return np.maximum(np.minimum(ltv, dsti), 0.0)499500501def solve_age(P: Params, j: int, pr: Prices, EV_own: np.ndarray, EV_rent: np.ndarray) -> AgePolicy:502 """Solve the age-j problem given continuation expectations EV_{j+1} (on the next-period grid)."""503 y = P.ybar[j - 1] * P.z_grid * pr.wage_factor504 nz, nh, nb, nN, na, nx = P.n_z, P.n_h, P.n_b, P.N, P.n_a, P.n_x505 ic_by_n = pr.i_c_by_n if (pr.i_c_by_n is not None and not pr.reset) else np.full(P.N, pr.i_new)506 maxloan = max_loan_matrix(P, pr, y)507 # --- branches508 W_S = np.empty((nz, nh, nb, nN, na)); aS = np.empty_like(W_S); cash_S = np.empty_like(W_S)509 _solve_stay(y, P.a_grid, P.b_grid, P.h_grid, np.ascontiguousarray(ic_by_n, dtype=float), pr.pi_next,510 pr.R_a, pr.p, P.m, P.chi, P.alpha, P.gamma, P.beta, P.c_min, EV_own, W_S, aS, cash_S)511 W_R = np.empty((nz, nx)); aR = np.empty_like(W_R)512 _solve_rent(P.x_grid, P.a_grid, pr.rho, P.alpha, P.gamma, P.beta, P.s_min, P.s_max, P.c_min,513 EV_rent, W_R, aR)514 rho_sit = pr.rho * pr.rent_rigid_factor515 if pr.rent_rigid_factor != 1.0:516 W_Rs = np.empty((nz, nx)); aRs = np.empty_like(W_Rs)517 _solve_rent(P.x_grid, P.a_grid, rho_sit, P.alpha, P.gamma, P.beta, P.s_min, P.s_max, P.c_min,518 EV_rent, W_Rs, aRs)519 else:520 W_Rs, aRs = W_R, aR521 W_Bh = np.empty((nz, nx, nh)); BBh = np.empty_like(W_Bh); aBh = np.empty_like(W_Bh)522 _solve_buy(P.x_grid, P.a_grid, P.b_grid, P.h_grid, P.B_fracs, maxloan, pr.i_new, pr.pi_next, pr.p,523 P.kappa_b, P.m, P.chi, P.alpha, P.gamma, P.beta, P.N, P.c_min, EV_own, W_Bh, BBh, aBh)524 W_F = np.empty((nz, nh, nx)); BF = np.empty_like(W_F); aF = np.empty_like(W_F)525 _solve_refi(P.x_grid, P.a_grid, P.b_grid, P.h_grid, P.B_fracs, maxloan, pr.i_new, pr.pi_next, pr.p,526 P.m, P.chi, P.alpha, P.gamma, P.beta, P.N, P.c_min, EV_own, W_F, BF, aF)527 # --- gift transform at the gift age528 gS = np.zeros_like(W_S); gF = np.zeros_like(W_F); gBh = np.zeros_like(W_Bh); gR = np.zeros_like(W_R)529 gRs = np.zeros_like(W_Rs)530 if j == P.j_gift and P.nu > 0.0:531 Wg = np.empty_like(W_S)532 _gift_transform_stay(W_S, cash_S, P.a_grid, pr.R_a, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg, gS)533 W_S = Wg534 Wg = np.empty_like(W_F); gF = np.empty_like(W_F)535 _gift_transform_x(W_F.reshape(nz * nh, nx), P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b,536 Wg.reshape(nz * nh, nx), gF.reshape(nz * nh, nx))537 W_F = Wg538 Wg = np.empty_like(W_Bh); gB3 = np.empty_like(W_Bh)539 WB_t = np.ascontiguousarray(np.transpose(W_Bh, (0, 2, 1))).reshape(nz * nh, nx)540 Wg_t = np.empty_like(WB_t); g_t = np.empty_like(WB_t)541 _gift_transform_x(WB_t, P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg_t, g_t)542 W_Bh = np.ascontiguousarray(np.transpose(Wg_t.reshape(nz, nh, nx), (0, 2, 1)))543 gBh = np.ascontiguousarray(np.transpose(g_t.reshape(nz, nh, nx), (0, 2, 1)))544 Wg = np.empty_like(W_R)545 _gift_transform_x(W_R, P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg, gR)546 W_R = Wg547 if pr.rent_rigid_factor != 1.0:548 Wg = np.empty_like(W_Rs)549 _gift_transform_x(W_Rs, P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg, gRs)550 W_Rs = Wg551 else:552 W_Rs, aRs, gRs = W_R, aR, gR553 # --- best house size (reporting / combine)554 hB = np.argmax(W_Bh, axis=2)555 W_B = np.max(W_Bh, axis=2)556 idx = (np.arange(nz)[:, None], np.arange(nx)[None, :], hB)557 BB = BBh[idx]; aB = aBh[idx]; gB = gBh[idx]558 hB = np.where(W_B > NEG * 0.5, hB, -1)559 # --- combine560 V_own = np.empty_like(W_S); choice_own = np.empty(W_S.shape, dtype=np.int64)561 _combine_owner(y, P.a_grid, P.b_grid, P.h_grid, P.x_grid, pr.R_a, pr.p, P.kappa_s, P.kappa_f, P.m,562 W_S, W_F, W_B, W_R, V_own, choice_own)563 V_rent = np.empty((nz, na)); choice_rent = np.empty((nz, na), dtype=np.int64)564 rebate = pr.transfer_rebate if j == 1 else 0.0565 _combine_renter(y, P.a_grid, P.x_grid, pr.R_a, rebate, W_B, W_Rs, V_rent, choice_rent)566 return AgePolicy(W_S=W_S, aS=aS, gS=gS, W_F=W_F, BF=BF, aF=aF, gF=gF, W_B=W_B, hB=hB, BB=BB, aB=aB, gB=gB,567 W_Bh=W_Bh, BBh=BBh, aBh=aBh, gBh=gBh,568 W_R=W_R, aR=aR, gR=gR, W_Rs=W_Rs, aRs=aRs, gRs=gRs, V_own=V_own, V_rent=V_rent,569 choice_own=choice_own, choice_rent=choice_rent, maxloan=maxloan, rho_sit=rho_sit)570571572def solve_lifecycle(P: Params, prices_by_age: list[Prices] | Prices,573 lottery: tuple[np.ndarray, np.ndarray] | None,574 V_next_terminal: tuple[np.ndarray, np.ndarray] | None = None,575 V_next_by_age: list[tuple[np.ndarray, np.ndarray]] | None = None) -> list[AgePolicy]:576 """Backward induction over ages j = J..1.577578 prices_by_age: one Prices object (stationary) or a list of J objects (date-specific prices579 for each age, used along a transition where the household of age j at date t faces prices580 at t and continuation V_{j+1, t+1}). V_next_by_age[j-1] = (V_own, V_rent) of age j+1 at the581 next date (transition); None -> stationary (use this solve's own V_{j+1}).582 Returns policies indexed by age-1.583 """584 stationary = not isinstance(prices_by_age, list)585 pols: list[AgePolicy | None] = [None] * P.J586 V_own_next = V_rent_next = None587 for j in range(P.J, 0, -1):588 pr = prices_by_age if stationary else prices_by_age[j - 1]589 if j == P.J:590 EV_own, EV_rent = terminal_continuation(P, pr)591 else:592 if V_next_by_age is not None and V_next_by_age[j - 1] is not None:593 Vo_n, Vr_n = V_next_by_age[j - 1]594 else:595 Vo_n, Vr_n = V_own_next, V_rent_next596 EV_own, EV_rent = expectation(P, j, Vo_n, Vr_n, lottery)597 pol = solve_age(P, j, pr, EV_own, EV_rent)598 pols[j - 1] = pol599 V_own_next, V_rent_next = pol.V_own, pol.V_rent600 return pols # type: ignore[return-value]601602603# --------------------------------------------------------------------------------------604# Diagnostics: ownership threshold, HTM check, Bellman residuals605# --------------------------------------------------------------------------------------606def ownership_threshold(P: Params, pol: AgePolicy, pr: Prices, j: int = 1) -> dict:607 """Smallest cash-on-hand / initial asset at which the renter's policy is "buy", by z, with a608 monotonicity flag (number of sign changes of W_B - W_R along the cash grid)."""609 out = {}610 y = P.ybar[j - 1] * P.z_grid * pr.wage_factor611 for iz in range(P.n_z):612 diff = pol.W_B[iz] - pol.W_Rs[iz]613 feas = (pol.W_B[iz] > NEG * 0.5) & (pol.W_Rs[iz] > NEG * 0.5)614 buy = (diff > 0) & feas615 idx = np.where(buy)[0]616 if len(idx) == 0:617 x_hat = np.nan618 else:619 k = idx[0]620 if k > 0 and feas[k - 1]:621 # linear root of diff between k-1 and k622 d0, d1 = diff[k - 1], diff[k]623 x_hat = P.x_grid[k - 1] + (P.x_grid[k] - P.x_grid[k - 1]) * (-d0 / (d1 - d0))624 else:625 x_hat = P.x_grid[k]626 sgn = np.sign(diff[feas])627 changes = int(np.sum(sgn[1:] != sgn[:-1])) if len(sgn) > 1 else 0628 a_hat = (x_hat - y[iz]) / pr.R_a if np.isfinite(x_hat) else np.nan629 # hand-to-mouth check at the threshold: a' of buyer and renter at x_hat630 if np.isfinite(x_hat):631 k = min(int(np.searchsorted(P.x_grid, x_hat)), P.n_x - 1)632 htm_buyer = float(pol.aB[iz, k]) <= P.a_grid[1]633 htm_renter = float(pol.aRs[iz, max(k - 1, 0)]) <= P.a_grid[1]634 feas_bound = (1.0 + P.kappa_b + P.m + P.theta) * pr.p * P.h_grid[0] - (1.0 - P.theta) * pr.p * P.h_grid[0] * 0.0635 else:636 htm_buyer = htm_renter = False637 feas_bound = np.nan638 out[iz] = dict(z=float(P.z_grid[iz]), x_hat=float(x_hat), a_hat=float(a_hat),639 sign_changes=changes, htm_buyer=bool(htm_buyer), htm_renter=bool(htm_renter),640 feasibility_cash=float((P.theta + P.kappa_b + P.m) * pr.p * P.h_grid[0]641 + annuity_factor(pr.i_new, P.N) * (1 - P.theta) * pr.p * P.h_grid[0]))642 return out643644645def bellman_residuals(P: Params, pols: list[AgePolicy], pr: Prices, lottery, n_points: int = 400,646 seed: int = 1, n_fine: int = 400) -> dict:647 """Relative Bellman residuals at random off-grid owner and renter states, re-optimising the648 liquid-asset choice on a fine grid (n_fine points) for the branch chosen on the grid."""649 rng = np.random.default_rng(seed)650 res_own, res_rent = [], []651 for j in range(1, P.J):652 if j == P.j_gift and P.nu > 0:653 continue # gift age excluded: the gift grid is not re-optimised here (reported separately)654 pol = pols[j - 1]655 EV_own, EV_rent = expectation(P, j, pols[j].V_own, pols[j].V_rent, lottery)656 y = P.ybar[j - 1] * P.z_grid * pr.wage_factor657 # --- renters: random (z, a)658 for _ in range(n_points // (P.J - 1)):659 iz = int(rng.integers(P.n_z))660 a = float(rng.uniform(0.0, P.a_grid[-1] * 0.3))661 x = y[iz] + a * pr.R_a + (pr.transfer_rebate if j == 1 else 0.0)662 k, w = _locate_py(P.x_grid, x)663 v_grid = _lin_py(pol.W_Rs[iz, k], pol.W_Rs[iz, k + 1], w)664 vb = _lin_py(pol.W_B[iz, k], pol.W_B[iz, k + 1], w)665 v_grid = max(v_grid, vb)666 # fine re-optimisation of the rent branch (and buy branch with grid h', B')667 v_fine = _fine_rent(P, pr, x, EV_rent[iz], n_fine, j)668 v_fine = max(v_fine, _fine_buy(P, pr, x, EV_own[iz], pol.maxloan[iz], n_fine, j))669 if v_grid > NEG * 0.5 and v_fine > NEG * 0.5:670 res_rent.append(abs(v_fine - v_grid) / abs(v_fine))671 # --- owners: random (z, h, b, n, a) with stay re-optimised672 for _ in range(n_points // (P.J - 1)):673 iz = int(rng.integers(P.n_z)); ih = int(rng.integers(P.n_h)); n_idx = int(rng.integers(P.N))674 b = float(rng.uniform(0.0, (1 - P.theta) * pr.p * P.h_grid[ih])); a = float(rng.uniform(0.0, P.a_grid[-1] * 0.3))675 jb, wb = _locate_py(P.b_grid, b); ja, wa = _locate_py(P.a_grid, a)676 v00 = pol.V_own[iz, ih, jb, n_idx, ja]; v01 = pol.V_own[iz, ih, jb, n_idx, ja + 1]677 v10 = pol.V_own[iz, ih, jb + 1, n_idx, ja]; v11 = pol.V_own[iz, ih, jb + 1, n_idx, ja + 1]678 if min(v00, v01, v10, v11) <= NEG * 0.5:679 continue680 v_grid = (1 - wb) * ((1 - wa) * v00 + wa * v01) + wb * ((1 - wa) * v10 + wa * v11)681 v_fine = _fine_stay(P, pr, iz, ih, b, n_idx + 1, a, EV_own, n_fine, j)682 x_sell = y[iz] + a * pr.R_a + (1 - P.kappa_s) * pr.p * P.h_grid[ih] - b683 v_fine = max(v_fine, _fine_rent(P, pr, x_sell, EV_rent[iz], n_fine, j),684 _fine_buy(P, pr, x_sell, EV_own[iz], pol.maxloan[iz], n_fine, j))685 if v_fine > NEG * 0.5:686 res_own.append(abs(v_fine - v_grid) / abs(v_fine))687 ro, rr = np.array(res_own), np.array(res_rent)688 return dict(n_owner=int(len(ro)), n_renter=int(len(rr)),689 owner_mean=float(ro.mean()) if len(ro) else np.nan, owner_max=float(ro.max()) if len(ro) else np.nan,690 renter_mean=float(rr.mean()) if len(rr) else np.nan, renter_max=float(rr.max()) if len(rr) else np.nan,691 owner_p99=float(np.quantile(ro, 0.99)) if len(ro) else np.nan,692 renter_p99=float(np.quantile(rr, 0.99)) if len(rr) else np.nan)693694695def _locate_py(grid, x):696 n = len(grid)697 if x <= grid[0]:698 return 0, 0.0699 if x >= grid[-1]:700 return n - 2, 1.0701 k = int(np.searchsorted(grid, x, side="right")) - 1702 return k, (x - grid[k]) / (grid[k + 1] - grid[k])703704705def _lin_py(v0, v1, w):706 if v0 <= NEG * 0.5 or v1 <= NEG * 0.5:707 return NEG708 return (1 - w) * v0 + w * v1709710711def _interp_a(P, EV_row, a):712 """EV_row shaped (..., na) -> interpolate in a (vector over leading dims)."""713 k, w = _locate_py(P.a_grid, a)714 v = (1 - w) * EV_row[..., k] + w * EV_row[..., k + 1]715 bad = (EV_row[..., k] <= NEG * 0.5) | (EV_row[..., k + 1] <= NEG * 0.5)716 return np.where(bad, NEG, v)717718719def _fine_rent(P, pr, x, EV_rent_z, n_fine, j):720 if x <= 0:721 return NEG722 rho = pr.rho * pr.rent_rigid_factor723 aps = np.linspace(0.0, min(x, P.a_grid[-1]), n_fine)724 e = x - aps725 s = np.clip(P.alpha * e / rho, P.s_min, P.s_max)726 c = e - rho * s727 ok = c > P.c_min728 if not ok.any():729 return NEG730 ev = np.array([_interp_a(P, EV_rent_z, ap) for ap in aps])731 v = np.where(ok & (ev > NEG * 0.5), (np.power(np.maximum(c, 1e-12), 1 - P.alpha) * s ** P.alpha) ** (1 - P.gamma) / (1 - P.gamma) + P.beta * ev, NEG)732 return float(v.max())733734735def _fine_buy(P, pr, x, EV_own_z, maxloan_z, n_fine, j):736 best = NEG737 af = annuity_factor(pr.i_new, P.N)738 for ih, h in enumerate(P.h_grid):739 for frac in np.linspace(0, 1, 2 * len(P.B_fracs) - 1):740 B = frac * maxloan_z[ih]741 pay = B * af742 bn = next_balance(B, P.N, pr.i_new, pr.pi_next)743 jb, wb = _locate_py(P.b_grid, bn)744 cash = x - (1 + P.kappa_b + P.m) * pr.p * h + B - pay745 if cash <= P.c_min:746 continue747 aps = np.linspace(0.0, min(cash - P.c_min, P.a_grid[-1]), n_fine)748 c = cash - aps749 EVrow = (1 - wb) * EV_own_z[ih, jb, P.N - 2, :] + wb * EV_own_z[ih, jb + 1, P.N - 2, :]750 EVrow = np.where((EV_own_z[ih, jb, P.N - 2, :] <= NEG * 0.5) | (EV_own_z[ih, jb + 1, P.N - 2, :] <= NEG * 0.5), NEG, EVrow)751 ev = np.interp(aps, P.a_grid, EVrow)752 okk = (c > P.c_min) & (ev > NEG * 0.5)753 if okk.any():754 v = (np.power(c[okk], 1 - P.alpha) * ((1 + P.chi) * h) ** P.alpha) ** (1 - P.gamma) / (1 - P.gamma) + P.beta * ev[okk]755 best = max(best, float(v.max()))756 return best757758759def _fine_stay(P, pr, iz, ih, b, n, a, EV_own, n_fine, j):760 ic = pr.i_c(n)761 pay = b * annuity_factor(ic, n)762 bn = next_balance(b, n, ic, pr.pi_next) # annuity-due timing (model.next_balance)763 jb, wb = _locate_py(P.b_grid, bn)764 y = P.ybar[j - 1] * P.z_grid[iz] * pr.wage_factor765 x = y + a * pr.R_a - pay - P.m * pr.p * P.h_grid[ih]766 if x <= P.c_min:767 return NEG768 n_next = max(n - 2, 0)769 EVrow = (1 - wb) * EV_own[iz, ih, jb, n_next, :] + wb * EV_own[iz, ih, jb + 1, n_next, :]770 EVrow = np.where((EV_own[iz, ih, jb, n_next, :] <= NEG * 0.5) | (EV_own[iz, ih, jb + 1, n_next, :] <= NEG * 0.5), NEG, EVrow)771 aps = np.linspace(0.0, min(x - P.c_min, P.a_grid[-1]), n_fine)772 c = x - aps773 ev = np.interp(aps, P.a_grid, EVrow)774 okk = (c > P.c_min) & (ev > NEG * 0.5)775 if not okk.any():776 return NEG777 s = (1 + P.chi) * P.h_grid[ih]778 v = (np.power(c[okk], 1 - P.alpha) * s ** P.alpha) ** (1 - P.gamma) / (1 - P.gamma) + P.beta * ev[okk]779 return float(v.max())780