"""WP18 — household dynamic programme (Part B), solved by backward induction over age. Method (justification in research/results_audit.md §algorithm): discrete-choice value function iteration with exhaustive search over the liquid-asset choice on the asset grid. The discrete choices (stay / refinance / move / sell-and-rent / rent / buy, house size, loan size, gift) make the problem non-concave; exhaustive search is robust to the kinks and, compiled with numba, is cheap enough (one age solves in ~50 ms on the baseline grids). EGM with upper envelopes (Iskhakov et al. 2017) would require taste shocks; we avoid them and report Bellman residuals at off-grid points instead (`bellman_residuals`). State: (j, z, h, b, n, a) with h the house index (-1 = renter), b the real balance at the start of the period, n the remaining payments, a the liquid real asset. All branches except "stay" depend on the state only through cash-on-hand x, so they are solved on a cash grid: W_R(z, x) rent W_B(z, x) buy (h', B', a') W_F(z, h, x) refinance W_S(z, h, b, n, a) stay (continuous a'). At the gift age the gift is chosen on top of each branch (gift transform). The estate at the estate age is anticipated as a K-point lottery (equilibrium object). Infeasible states carry the sentinel NEG = -1e10. """ from __future__ import annotations from dataclasses import dataclass import numpy as np from numba import njit, prange from .model import (Params, Prices, annuity_factor, next_balance, warm_glow) NEG = -1.0e10 # -------------------------------------------------------------------------------------- # Small numba helpers # -------------------------------------------------------------------------------------- @njit(cache=True, inline="always") def _locate(grid, x): """Index i and weight w such that x ~ (1-w) grid[i] + w grid[i+1]; clamped at the ends.""" n = grid.shape[0] if x <= grid[0]: return 0, 0.0 if x >= grid[n - 1]: return n - 2, 1.0 lo = 0 hi = n - 1 while hi - lo > 1: mid = (lo + hi) // 2 if grid[mid] <= x: lo = mid else: hi = mid w = (x - grid[lo]) / (grid[lo + 1] - grid[lo]) return lo, w @njit(cache=True, inline="always") def _lin(v0, v1, w): if v0 <= NEG * 0.5 or v1 <= NEG * 0.5: return NEG return (1.0 - w) * v0 + w * v1 @njit(cache=True, inline="always") def _u(c, s, alpha, gamma): comp = c ** (1.0 - alpha) * s ** alpha return comp ** (1.0 - gamma) / (1.0 - gamma) @njit(cache=True, inline="always") def _renter_cs(e, rho, alpha, s_min, s_max): """Intratemporal renter choice given expenditure e = c + rho s.""" s = alpha * e / rho if s < s_min: s = s_min elif s > s_max: s = s_max c = e - rho * s return c, s # -------------------------------------------------------------------------------------- # Branch solvers # -------------------------------------------------------------------------------------- GOLD = 0.6180339887498949 N_GOLD = 18 @njit(cache=True, inline="always") def _af_due(i, n): """Annuity-due factor: payment per unit of balance, n payments, first one now.""" if n <= 0: return 0.0 if abs(i) < 1e-12: return 1.0 / n return i / ((1.0 + i) * (1.0 - (1.0 + i) ** (-n))) @njit(cache=True, inline="always") def _ev_interp(EV2, jb, wb, a_grid, ap): """EV2 is a (nb, na) view; interpolate in b (jb, wb) and in a (ap).""" ja, wa = _locate(a_grid, ap) v0 = _lin(EV2[jb, ja], EV2[jb, ja + 1], wa) v1 = _lin(EV2[jb + 1, ja], EV2[jb + 1, ja + 1], wa) return _lin(v0, v1, wb) @njit(cache=True, inline="always") def _obj(ap, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max): """Period objective of saving ap. mode 0: owner (s given); mode 1: renter (rho given).""" e = cash - ap if mode == 0: c = e ss = s else: c, ss = _renter_cs(e, rho, alpha, s_min, s_max) if c <= c_min: return NEG ev = _ev_interp(EV2, jb, wb, a_grid, ap) if ev <= NEG * 0.5: return NEG return _u(c, ss, alpha, gamma) + beta * ev @njit(cache=True) def _best_a(cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max): """Grid search over a' on a_grid, then golden-section refinement between the neighbours.""" na = a_grid.shape[0] best = NEG kbest = -1 for ka in range(na): if cash - a_grid[ka] <= c_min: break v = _obj(a_grid[ka], cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max) if v > best: best = v kbest = ka if kbest < 0: return NEG, 0.0 lo = a_grid[kbest - 1] if kbest > 0 else a_grid[0] hi = a_grid[kbest + 1] if kbest < na - 1 else a_grid[na - 1] if cash - c_min < hi: hi = cash - c_min if hi <= lo: return best, a_grid[kbest] x1 = hi - GOLD * (hi - lo) x2 = lo + GOLD * (hi - lo) f1 = _obj(x1, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max) f2 = _obj(x2, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max) for _ in range(N_GOLD): if f1 < f2: lo = x1 x1 = x2 f1 = f2 x2 = lo + GOLD * (hi - lo) f2 = _obj(x2, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max) else: hi = x2 x2 = x1 f2 = f1 x1 = hi - GOLD * (hi - lo) f1 = _obj(x1, cash, s, rho, mode, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, s_min, s_max) abest = a_grid[kbest] if f1 > best: best = f1 abest = x1 if f2 > best: best = f2 abest = x2 return best, abest @njit(cache=True, parallel=True) def _solve_stay(y, a_grid, b_grid, h_grid, ic_by_n, pi_next, R_a, p, m, chi, alpha, gamma, beta, c_min, EV_own, W_S, aS, cash_S): """Owner stays: c + pay + m p h + a' = y + a R_a. EV_own indexed by (z, h, b', n', a').""" nz, nh, nb, nN, na = EV_own.shape for iz in prange(nz): for ih in range(nh): h = h_grid[ih] s = (1.0 + chi) * h for ib in range(nb): b = b_grid[ib] for n_idx in range(nN): n = n_idx + 1 ic = ic_by_n[n_idx] pay = b * _af_due(ic, n) bn = b - pay if bn < 0.0: bn = 0.0 bn = bn * (1.0 + ic) / (1.0 + pi_next) jb, wb = _locate(b_grid, bn) n_next = n_idx - 1 if n_idx > 0 else 0 EV2 = EV_own[iz, ih, :, n_next, :] for ia in range(na): x = y[iz] + a_grid[ia] * R_a - pay - m * p * h cash_S[iz, ih, ib, n_idx, ia] = x best, abest = _best_a(x, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0) W_S[iz, ih, ib, n_idx, ia] = best aS[iz, ih, ib, n_idx, ia] = abest @njit(cache=True, parallel=True) def _solve_rent(x_grid, a_grid, rho, alpha, gamma, beta, s_min, s_max, c_min, EV_rent, W_R, aR): nz = EV_rent.shape[0] nx = x_grid.shape[0] for iz in prange(nz): EV2 = EV_rent[iz:iz + 1, :] for ix in range(nx): 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) W_R[iz, ix] = best aR[iz, ix] = abest @njit(cache=True, parallel=True) def _solve_buy(x_grid, a_grid, b_grid, h_grid, B_fracs, maxloan, i_new, pi_next, p, kappa_b, m, chi, alpha, gamma, beta, N, c_min, EV_own, W_Bh, BBh, aBh): """Renter buys / owner moves: c + (1+kb+m) p h' - B' + pay(B') + a' = x. maxloan[z, h']. Outputs are per house size h' (W_Bh[z, x, h'] = best over B', a') so that the size choice can be interpolated in cash at simulation time (avoids synchronised jumps of the demand).""" nz = EV_own.shape[0] nh = h_grid.shape[0] nx = x_grid.shape[0] nB = B_fracs.shape[0] af = _af_due(i_new, N) for iz in prange(nz): for ix in range(nx): x = x_grid[ix] for ih in range(nh): h = h_grid[ih] s = (1.0 + chi) * h outlay = (1.0 + kappa_b + m) * p * h EV2 = EV_own[iz, ih, :, N - 2, :] best = NEG bestB = 0.0 besta = 0.0 for kB in range(nB): B = B_fracs[kB] * maxloan[iz, ih] pay = B * af bn = (B - pay) * (1.0 + i_new) / (1.0 + pi_next) jb, wb = _locate(b_grid, bn) cash = x - outlay + B - pay if cash <= c_min: continue v, ab = _best_a(cash, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0) if v > best: best = v bestB = B besta = ab W_Bh[iz, ix, ih] = best BBh[iz, ix, ih] = bestB aBh[iz, ix, ih] = besta @njit(cache=True, parallel=True) def _solve_refi(x_grid, a_grid, b_grid, h_grid, B_fracs, maxloan, i_new, pi_next, p, m, chi, alpha, gamma, beta, N, c_min, EV_own, W_F, BF, aF): """Owner refinances (keeps h): c + a' = x + B' - pay(B'), x = y + a R_a - b - (kf + m) p h.""" nz = EV_own.shape[0] nh = h_grid.shape[0] nx = x_grid.shape[0] nB = B_fracs.shape[0] af = _af_due(i_new, N) for iz in prange(nz): for ih in range(nh): h = h_grid[ih] s = (1.0 + chi) * h EV2 = EV_own[iz, ih, :, N - 2, :] for ix in range(nx): x = x_grid[ix] best = NEG bestB = 0.0 besta = 0.0 for kB in range(nB): B = B_fracs[kB] * maxloan[iz, ih] pay = B * af bn = (B - pay) * (1.0 + i_new) / (1.0 + pi_next) jb, wb = _locate(b_grid, bn) cash = x + B - pay if cash <= c_min: continue v, ab = _best_a(cash, s, 0.0, 0, EV2, jb, wb, a_grid, alpha, gamma, beta, c_min, 0.0, 0.0) if v > best: best = v bestB = B besta = ab W_F[iz, ih, ix] = best BF[iz, ih, ix] = bestB aF[iz, ih, ix] = besta @njit(cache=True) def _gift_transform_x(W, x_grid, gift_fracs, nu, x_shift, gamma, tau_b, Wg, Gg): """Wg(z[,h], x) = max_f phi((1-tau) f x) + W(z[,h], (1-f) x); W shaped (..., nx).""" lead = W.shape[0] nx = x_grid.shape[0] nf = gift_fracs.shape[0] for il in range(lead): for ix in range(nx): x = x_grid[ix] best = NEG bestg = 0.0 for kf in range(nf): g = gift_fracs[kf] * x net = (1.0 - tau_b) * g phi = nu * (net + x_shift) ** (1.0 - gamma) / (1.0 - gamma) if nu > 0.0 else 0.0 xr = x - g j, w = _locate(x_grid, xr) wv = _lin(W[il, j], W[il, j + 1], w) if wv <= NEG * 0.5: continue v = phi + wv if v > best: best = v bestg = g Wg[il, ix] = best Gg[il, ix] = bestg @njit(cache=True, parallel=True) def _gift_transform_stay(W_S, cash_S, a_grid, R_a, gift_fracs, nu, x_shift, gamma, tau_b, Wg, Gg): """Stay branch: giving g lowers cash by g, equivalent to a lower asset a - g/R_a.""" nz, nh, nb, nN, na = W_S.shape nf = gift_fracs.shape[0] for iz in prange(nz): for ih in range(nh): for ib in range(nb): for n_idx in range(nN): for ia in range(na): x = cash_S[iz, ih, ib, n_idx, ia] best = NEG bestg = 0.0 for kf in range(nf): g = gift_fracs[kf] * x if x > 0.0 else 0.0 if g < 0.0: g = 0.0 net = (1.0 - tau_b) * g phi = nu * (net + x_shift) ** (1.0 - gamma) / (1.0 - gamma) if nu > 0.0 else 0.0 ar = a_grid[ia] - g / R_a if ar < 0.0: continue j, w = _locate(a_grid, ar) wv = _lin(W_S[iz, ih, ib, n_idx, j], W_S[iz, ih, ib, n_idx, j + 1], w) if wv <= NEG * 0.5: continue v = phi + wv if v > best: best = v bestg = g if kf == 0 and nu == 0.0: break Wg[iz, ih, ib, n_idx, ia] = best Gg[iz, ih, ib, n_idx, ia] = bestg @njit(cache=True, parallel=True) def _combine_owner(y, a_grid, b_grid, h_grid, x_grid, R_a, p, kappa_s, kappa_f, m, W_S, W_F, W_B, W_R, V_own, choice_own): """V_own = max(stay, refi(x_refi), buy(x_sell), rent(x_sell)); choice: 0 S, 1 F, 2 B, 3 R, -1 none.""" nz, nh, nb, nN, na = W_S.shape for iz in prange(nz): for ih in range(nh): h = h_grid[ih] for ib in range(nb): b = b_grid[ib] for n_idx in range(nN): for ia in range(na): base = y[iz] + a_grid[ia] * R_a x_sell = base + (1.0 - kappa_s) * p * h - b x_refi = base - b - (kappa_f + m) * p * h best = W_S[iz, ih, ib, n_idx, ia] ch = 0 if best > NEG * 0.5 else -1 if x_refi > 0.0: j, w = _locate(x_grid, x_refi) v = _lin(W_F[iz, ih, j], W_F[iz, ih, j + 1], w) if v > best: best = v ch = 1 if x_sell > 0.0: j, w = _locate(x_grid, x_sell) v = _lin(W_B[iz, j], W_B[iz, j + 1], w) if v > best: best = v ch = 2 v = _lin(W_R[iz, j], W_R[iz, j + 1], w) if v > best: best = v ch = 3 V_own[iz, ih, ib, n_idx, ia] = best choice_own[iz, ih, ib, n_idx, ia] = ch @njit(cache=True) def _combine_renter(y, a_grid, x_grid, R_a, rebate, W_B, W_R_sit, V_rent, choice_rent): nz = W_B.shape[0] na = a_grid.shape[0] for iz in range(nz): for ia in range(na): x = y[iz] + a_grid[ia] * R_a + rebate j, w = _locate(x_grid, x) vb = _lin(W_B[iz, j], W_B[iz, j + 1], w) vr = _lin(W_R_sit[iz, j], W_R_sit[iz, j + 1], w) if vr >= vb: V_rent[iz, ia] = vr choice_rent[iz, ia] = 3 if vr > NEG * 0.5 else -1 else: V_rent[iz, ia] = vb choice_rent[iz, ia] = 2 @njit(cache=True) def _expect_lottery(V, a_grid, lot_pts, lot_w, out): """out[..., a] = sum_k w_k V[..., a + e_k] (linear interpolation / flat extrapolation in a).""" shp = V.shape na = shp[-1] V2 = V.reshape((-1, na)) o2 = out.reshape((-1, na)) K = lot_pts.shape[0] for r in range(V2.shape[0]): for ia in range(na): acc = 0.0 bad = False for k in range(K): j, w = _locate(a_grid, a_grid[ia] + lot_pts[k]) v = _lin(V2[r, j], V2[r, j + 1], w) if v <= NEG * 0.5: bad = True break acc += lot_w[k] * v o2[r, ia] = NEG if bad else acc # -------------------------------------------------------------------------------------- # Python-level age solve # -------------------------------------------------------------------------------------- @dataclass class AgePolicy: """Branch values and policies for one age (and one date along a transition).""" 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 # best over h' (reporting) W_Bh: np.ndarray; BBh: np.ndarray; aBh: np.ndarray; gBh: np.ndarray # per house size (simulation) W_R: np.ndarray; aR: np.ndarray; gR: np.ndarray # new-tenant rent branch W_Rs: np.ndarray; aRs: np.ndarray; gRs: np.ndarray # sitting-tenant rent branch V_own: np.ndarray; V_rent: np.ndarray choice_own: np.ndarray; choice_rent: np.ndarray maxloan: np.ndarray rho_sit: float def to_f32(self) -> "AgePolicy": kw = {} for k, v in self.__dict__.items(): kw[k] = v.astype(np.float32) if isinstance(v, np.ndarray) and v.dtype == np.float64 else v return AgePolicy(**kw) def expectation(P: Params, j: int, V_own_next: np.ndarray, V_rent_next: np.ndarray, lottery: tuple[np.ndarray, np.ndarray] | None) -> tuple[np.ndarray, np.ndarray]: """E[V_{j+1} | z] for a household of age j (1-based): Markov in z (frozen in retirement), estate lottery if j + 1 == j_estate.""" if lottery is not None and j + 1 == P.j_estate: pts, w = lottery Vo = np.empty_like(V_own_next); Vr = np.empty_like(V_rent_next) _expect_lottery(V_own_next, P.a_grid, pts, w, Vo) _expect_lottery(V_rent_next, P.a_grid, pts, w, Vr) else: Vo, Vr = V_own_next, V_rent_next if j + 1 >= P.j_retire: return Vo.copy(), Vr.copy() Pi = P.Pi_z bad_o = Vo <= NEG * 0.5 bad_r = Vr <= NEG * 0.5 EVo = np.tensordot(Pi, np.where(bad_o, 0.0, Vo), axes=(1, 0)) EVr = Pi @ np.where(bad_r, 0.0, Vr) # a state is infeasible in expectation if any reachable successor is infeasible EVo = np.where(np.tensordot(Pi, bad_o.astype(float), axes=(1, 0)) > 0, NEG, EVo) EVr = np.where((Pi @ bad_r.astype(float)) > 0, NEG, EVr) return EVo, EVr def terminal_continuation(P: Params, pr: Prices) -> tuple[np.ndarray, np.ndarray]: """Warm glow of the (net) estate as a function of the next-period state (h, b', a').""" p_next = pr.p_next if pr.p_next is not None else pr.p a = P.a_grid[None, None, None, None, :] b = P.b_grid[None, None, :, None, None] h = P.h_grid[None, :, None, None, None] estate_own = (1.0 - P.tau_b) * (a + (1.0 - P.kappa_s) * p_next * h - b) estate_own = np.broadcast_to(estate_own, (P.n_z, P.n_h, P.n_b, P.N, P.n_a)) EV_own = np.where(estate_own >= 0.0, warm_glow(P, np.maximum(estate_own, 0.0)), NEG) EV_rent = np.broadcast_to(warm_glow(P, (1.0 - P.tau_b) * P.a_grid)[None, :], (P.n_z, P.n_a)) return np.ascontiguousarray(EV_own, dtype=float), np.ascontiguousarray(EV_rent, dtype=float) def max_loan_matrix(P: Params, pr: Prices, y: np.ndarray) -> np.ndarray: """maxloan[z, h'] = min(LTV cap, DSTI cap at the qualifying rate).""" af_q = annuity_factor(pr.i_q, P.N) ltv = (1.0 - P.theta) * pr.p * P.h_grid[None, :] dsti = (P.psi * y / af_q)[:, None] if af_q > 0 else np.full((P.n_z, 1), np.inf) return np.maximum(np.minimum(ltv, dsti), 0.0) def solve_age(P: Params, j: int, pr: Prices, EV_own: np.ndarray, EV_rent: np.ndarray) -> AgePolicy: """Solve the age-j problem given continuation expectations EV_{j+1} (on the next-period grid).""" y = P.ybar[j - 1] * P.z_grid * pr.wage_factor nz, nh, nb, nN, na, nx = P.n_z, P.n_h, P.n_b, P.N, P.n_a, P.n_x 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) maxloan = max_loan_matrix(P, pr, y) # --- branches W_S = np.empty((nz, nh, nb, nN, na)); aS = np.empty_like(W_S); cash_S = np.empty_like(W_S) _solve_stay(y, P.a_grid, P.b_grid, P.h_grid, np.ascontiguousarray(ic_by_n, dtype=float), pr.pi_next, 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) W_R = np.empty((nz, nx)); aR = np.empty_like(W_R) _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, EV_rent, W_R, aR) rho_sit = pr.rho * pr.rent_rigid_factor if pr.rent_rigid_factor != 1.0: W_Rs = np.empty((nz, nx)); aRs = np.empty_like(W_Rs) _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, EV_rent, W_Rs, aRs) else: W_Rs, aRs = W_R, aR W_Bh = np.empty((nz, nx, nh)); BBh = np.empty_like(W_Bh); aBh = np.empty_like(W_Bh) _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, P.kappa_b, P.m, P.chi, P.alpha, P.gamma, P.beta, P.N, P.c_min, EV_own, W_Bh, BBh, aBh) W_F = np.empty((nz, nh, nx)); BF = np.empty_like(W_F); aF = np.empty_like(W_F) _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, P.m, P.chi, P.alpha, P.gamma, P.beta, P.N, P.c_min, EV_own, W_F, BF, aF) # --- gift transform at the gift age gS = np.zeros_like(W_S); gF = np.zeros_like(W_F); gBh = np.zeros_like(W_Bh); gR = np.zeros_like(W_R) gRs = np.zeros_like(W_Rs) if j == P.j_gift and P.nu > 0.0: Wg = np.empty_like(W_S) _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) W_S = Wg Wg = np.empty_like(W_F); gF = np.empty_like(W_F) _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, Wg.reshape(nz * nh, nx), gF.reshape(nz * nh, nx)) W_F = Wg Wg = np.empty_like(W_Bh); gB3 = np.empty_like(W_Bh) WB_t = np.ascontiguousarray(np.transpose(W_Bh, (0, 2, 1))).reshape(nz * nh, nx) Wg_t = np.empty_like(WB_t); g_t = np.empty_like(WB_t) _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) W_Bh = np.ascontiguousarray(np.transpose(Wg_t.reshape(nz, nh, nx), (0, 2, 1))) gBh = np.ascontiguousarray(np.transpose(g_t.reshape(nz, nh, nx), (0, 2, 1))) Wg = np.empty_like(W_R) _gift_transform_x(W_R, P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg, gR) W_R = Wg if pr.rent_rigid_factor != 1.0: Wg = np.empty_like(W_Rs) _gift_transform_x(W_Rs, P.x_grid, P.gift_fracs, P.nu, P.x_shift, P.gamma, P.tau_b, Wg, gRs) W_Rs = Wg else: W_Rs, aRs, gRs = W_R, aR, gR # --- best house size (reporting / combine) hB = np.argmax(W_Bh, axis=2) W_B = np.max(W_Bh, axis=2) idx = (np.arange(nz)[:, None], np.arange(nx)[None, :], hB) BB = BBh[idx]; aB = aBh[idx]; gB = gBh[idx] hB = np.where(W_B > NEG * 0.5, hB, -1) # --- combine V_own = np.empty_like(W_S); choice_own = np.empty(W_S.shape, dtype=np.int64) _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, W_S, W_F, W_B, W_R, V_own, choice_own) V_rent = np.empty((nz, na)); choice_rent = np.empty((nz, na), dtype=np.int64) rebate = pr.transfer_rebate if j == 1 else 0.0 _combine_renter(y, P.a_grid, P.x_grid, pr.R_a, rebate, W_B, W_Rs, V_rent, choice_rent) 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, W_Bh=W_Bh, BBh=BBh, aBh=aBh, gBh=gBh, W_R=W_R, aR=aR, gR=gR, W_Rs=W_Rs, aRs=aRs, gRs=gRs, V_own=V_own, V_rent=V_rent, choice_own=choice_own, choice_rent=choice_rent, maxloan=maxloan, rho_sit=rho_sit) def solve_lifecycle(P: Params, prices_by_age: list[Prices] | Prices, lottery: tuple[np.ndarray, np.ndarray] | None, V_next_terminal: tuple[np.ndarray, np.ndarray] | None = None, V_next_by_age: list[tuple[np.ndarray, np.ndarray]] | None = None) -> list[AgePolicy]: """Backward induction over ages j = J..1. prices_by_age: one Prices object (stationary) or a list of J objects (date-specific prices for each age, used along a transition where the household of age j at date t faces prices at t and continuation V_{j+1, t+1}). V_next_by_age[j-1] = (V_own, V_rent) of age j+1 at the next date (transition); None -> stationary (use this solve's own V_{j+1}). Returns policies indexed by age-1. """ stationary = not isinstance(prices_by_age, list) pols: list[AgePolicy | None] = [None] * P.J V_own_next = V_rent_next = None for j in range(P.J, 0, -1): pr = prices_by_age if stationary else prices_by_age[j - 1] if j == P.J: EV_own, EV_rent = terminal_continuation(P, pr) else: if V_next_by_age is not None and V_next_by_age[j - 1] is not None: Vo_n, Vr_n = V_next_by_age[j - 1] else: Vo_n, Vr_n = V_own_next, V_rent_next EV_own, EV_rent = expectation(P, j, Vo_n, Vr_n, lottery) pol = solve_age(P, j, pr, EV_own, EV_rent) pols[j - 1] = pol V_own_next, V_rent_next = pol.V_own, pol.V_rent return pols # type: ignore[return-value] # -------------------------------------------------------------------------------------- # Diagnostics: ownership threshold, HTM check, Bellman residuals # -------------------------------------------------------------------------------------- def ownership_threshold(P: Params, pol: AgePolicy, pr: Prices, j: int = 1) -> dict: """Smallest cash-on-hand / initial asset at which the renter's policy is "buy", by z, with a monotonicity flag (number of sign changes of W_B - W_R along the cash grid).""" out = {} y = P.ybar[j - 1] * P.z_grid * pr.wage_factor for iz in range(P.n_z): diff = pol.W_B[iz] - pol.W_Rs[iz] feas = (pol.W_B[iz] > NEG * 0.5) & (pol.W_Rs[iz] > NEG * 0.5) buy = (diff > 0) & feas idx = np.where(buy)[0] if len(idx) == 0: x_hat = np.nan else: k = idx[0] if k > 0 and feas[k - 1]: # linear root of diff between k-1 and k d0, d1 = diff[k - 1], diff[k] x_hat = P.x_grid[k - 1] + (P.x_grid[k] - P.x_grid[k - 1]) * (-d0 / (d1 - d0)) else: x_hat = P.x_grid[k] sgn = np.sign(diff[feas]) changes = int(np.sum(sgn[1:] != sgn[:-1])) if len(sgn) > 1 else 0 a_hat = (x_hat - y[iz]) / pr.R_a if np.isfinite(x_hat) else np.nan # hand-to-mouth check at the threshold: a' of buyer and renter at x_hat if np.isfinite(x_hat): k = min(int(np.searchsorted(P.x_grid, x_hat)), P.n_x - 1) htm_buyer = float(pol.aB[iz, k]) <= P.a_grid[1] htm_renter = float(pol.aRs[iz, max(k - 1, 0)]) <= P.a_grid[1] 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.0 else: htm_buyer = htm_renter = False feas_bound = np.nan out[iz] = dict(z=float(P.z_grid[iz]), x_hat=float(x_hat), a_hat=float(a_hat), sign_changes=changes, htm_buyer=bool(htm_buyer), htm_renter=bool(htm_renter), feasibility_cash=float((P.theta + P.kappa_b + P.m) * pr.p * P.h_grid[0] + annuity_factor(pr.i_new, P.N) * (1 - P.theta) * pr.p * P.h_grid[0])) return out def bellman_residuals(P: Params, pols: list[AgePolicy], pr: Prices, lottery, n_points: int = 400, seed: int = 1, n_fine: int = 400) -> dict: """Relative Bellman residuals at random off-grid owner and renter states, re-optimising the liquid-asset choice on a fine grid (n_fine points) for the branch chosen on the grid.""" rng = np.random.default_rng(seed) res_own, res_rent = [], [] for j in range(1, P.J): if j == P.j_gift and P.nu > 0: continue # gift age excluded: the gift grid is not re-optimised here (reported separately) pol = pols[j - 1] EV_own, EV_rent = expectation(P, j, pols[j].V_own, pols[j].V_rent, lottery) y = P.ybar[j - 1] * P.z_grid * pr.wage_factor # --- renters: random (z, a) for _ in range(n_points // (P.J - 1)): iz = int(rng.integers(P.n_z)) a = float(rng.uniform(0.0, P.a_grid[-1] * 0.3)) x = y[iz] + a * pr.R_a + (pr.transfer_rebate if j == 1 else 0.0) k, w = _locate_py(P.x_grid, x) v_grid = _lin_py(pol.W_Rs[iz, k], pol.W_Rs[iz, k + 1], w) vb = _lin_py(pol.W_B[iz, k], pol.W_B[iz, k + 1], w) v_grid = max(v_grid, vb) # fine re-optimisation of the rent branch (and buy branch with grid h', B') v_fine = _fine_rent(P, pr, x, EV_rent[iz], n_fine, j) v_fine = max(v_fine, _fine_buy(P, pr, x, EV_own[iz], pol.maxloan[iz], n_fine, j)) if v_grid > NEG * 0.5 and v_fine > NEG * 0.5: res_rent.append(abs(v_fine - v_grid) / abs(v_fine)) # --- owners: random (z, h, b, n, a) with stay re-optimised for _ in range(n_points // (P.J - 1)): iz = int(rng.integers(P.n_z)); ih = int(rng.integers(P.n_h)); n_idx = int(rng.integers(P.N)) 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)) jb, wb = _locate_py(P.b_grid, b); ja, wa = _locate_py(P.a_grid, a) v00 = pol.V_own[iz, ih, jb, n_idx, ja]; v01 = pol.V_own[iz, ih, jb, n_idx, ja + 1] v10 = pol.V_own[iz, ih, jb + 1, n_idx, ja]; v11 = pol.V_own[iz, ih, jb + 1, n_idx, ja + 1] if min(v00, v01, v10, v11) <= NEG * 0.5: continue v_grid = (1 - wb) * ((1 - wa) * v00 + wa * v01) + wb * ((1 - wa) * v10 + wa * v11) v_fine = _fine_stay(P, pr, iz, ih, b, n_idx + 1, a, EV_own, n_fine, j) x_sell = y[iz] + a * pr.R_a + (1 - P.kappa_s) * pr.p * P.h_grid[ih] - b v_fine = max(v_fine, _fine_rent(P, pr, x_sell, EV_rent[iz], n_fine, j), _fine_buy(P, pr, x_sell, EV_own[iz], pol.maxloan[iz], n_fine, j)) if v_fine > NEG * 0.5: res_own.append(abs(v_fine - v_grid) / abs(v_fine)) ro, rr = np.array(res_own), np.array(res_rent) return dict(n_owner=int(len(ro)), n_renter=int(len(rr)), owner_mean=float(ro.mean()) if len(ro) else np.nan, owner_max=float(ro.max()) if len(ro) else np.nan, renter_mean=float(rr.mean()) if len(rr) else np.nan, renter_max=float(rr.max()) if len(rr) else np.nan, owner_p99=float(np.quantile(ro, 0.99)) if len(ro) else np.nan, renter_p99=float(np.quantile(rr, 0.99)) if len(rr) else np.nan) def _locate_py(grid, x): n = len(grid) if x <= grid[0]: return 0, 0.0 if x >= grid[-1]: return n - 2, 1.0 k = int(np.searchsorted(grid, x, side="right")) - 1 return k, (x - grid[k]) / (grid[k + 1] - grid[k]) def _lin_py(v0, v1, w): if v0 <= NEG * 0.5 or v1 <= NEG * 0.5: return NEG return (1 - w) * v0 + w * v1 def _interp_a(P, EV_row, a): """EV_row shaped (..., na) -> interpolate in a (vector over leading dims).""" k, w = _locate_py(P.a_grid, a) v = (1 - w) * EV_row[..., k] + w * EV_row[..., k + 1] bad = (EV_row[..., k] <= NEG * 0.5) | (EV_row[..., k + 1] <= NEG * 0.5) return np.where(bad, NEG, v) def _fine_rent(P, pr, x, EV_rent_z, n_fine, j): if x <= 0: return NEG rho = pr.rho * pr.rent_rigid_factor aps = np.linspace(0.0, min(x, P.a_grid[-1]), n_fine) e = x - aps s = np.clip(P.alpha * e / rho, P.s_min, P.s_max) c = e - rho * s ok = c > P.c_min if not ok.any(): return NEG ev = np.array([_interp_a(P, EV_rent_z, ap) for ap in aps]) 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) return float(v.max()) def _fine_buy(P, pr, x, EV_own_z, maxloan_z, n_fine, j): best = NEG af = annuity_factor(pr.i_new, P.N) for ih, h in enumerate(P.h_grid): for frac in np.linspace(0, 1, 2 * len(P.B_fracs) - 1): B = frac * maxloan_z[ih] pay = B * af bn = next_balance(B, P.N, pr.i_new, pr.pi_next) jb, wb = _locate_py(P.b_grid, bn) cash = x - (1 + P.kappa_b + P.m) * pr.p * h + B - pay if cash <= P.c_min: continue aps = np.linspace(0.0, min(cash - P.c_min, P.a_grid[-1]), n_fine) c = cash - aps EVrow = (1 - wb) * EV_own_z[ih, jb, P.N - 2, :] + wb * EV_own_z[ih, jb + 1, P.N - 2, :] 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) ev = np.interp(aps, P.a_grid, EVrow) okk = (c > P.c_min) & (ev > NEG * 0.5) if okk.any(): v = (np.power(c[okk], 1 - P.alpha) * ((1 + P.chi) * h) ** P.alpha) ** (1 - P.gamma) / (1 - P.gamma) + P.beta * ev[okk] best = max(best, float(v.max())) return best def _fine_stay(P, pr, iz, ih, b, n, a, EV_own, n_fine, j): ic = pr.i_c(n) pay = b * annuity_factor(ic, n) bn = next_balance(b, n, ic, pr.pi_next) # annuity-due timing (model.next_balance) jb, wb = _locate_py(P.b_grid, bn) y = P.ybar[j - 1] * P.z_grid[iz] * pr.wage_factor x = y + a * pr.R_a - pay - P.m * pr.p * P.h_grid[ih] if x <= P.c_min: return NEG n_next = max(n - 2, 0) EVrow = (1 - wb) * EV_own[iz, ih, jb, n_next, :] + wb * EV_own[iz, ih, jb + 1, n_next, :] 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) aps = np.linspace(0.0, min(x - P.c_min, P.a_grid[-1]), n_fine) c = x - aps ev = np.interp(aps, P.a_grid, EVrow) okk = (c > P.c_min) & (ev > NEG * 0.5) if not okk.any(): return NEG s = (1 + P.chi) * P.h_grid[ih] v = (np.power(c[okk], 1 - P.alpha) * s ** P.alpha) ** (1 - P.gamma) / (1 - P.gamma) + P.beta * ev[okk] return float(v.max())