TeX 53.2%
Python 46.6%
1"""WP18 — model primitives for the Part B lifecycle/OLG housing economy.23Everything is in real terms (units of the consumption good), one period = 5 years.4This module holds: configuration loading, the parameter object (with the annual -> period5conversions), the grids, the income process (Rouwenhorst), utility functions, mortgage6arithmetic (constant nominal payment, real balance law of motion, LTV and DSTI caps) and the7construction / rental-arbitrage relations. No optimisation happens here.89Conventions (research/model_specification.md §0, §B):10 * p real price of one unit of owner-occupied housing; p_r real price of a rental unit;11 * rho real rent per unit of rental housing services per period;12 * b real outstanding mortgage balance at the START of the period (before interest);13 * n remaining number of payments (1..N); b == 0 means no mortgage whatever n is;14 * pay real payment made at the START of the period (annuity-due on the outstanding balance):15 pay = b * i_c / ((1+i_c) (1 - (1+i_c)^-n)), exact for a constant-nominal-payment loan;16 start-of-period timing makes borrowing and depositing present-value consistent17 (borrowing B and paying pay now leaves (B - pay) R real next period, like a deposit);18 * b' = (b - pay) (1+i_c) / (1 + pi') real balance next period.19"""20from __future__ import annotations2122import copy23import math24from dataclasses import dataclass, field25from pathlib import Path26from typing import Any2728import numpy as np29import yaml3031ROOT = Path(__file__).resolve().parents[1]32CONFIG_DIR = ROOT / "configs"33RESULTS_DIR = ROOT / "results"343536# --------------------------------------------------------------------------------------37# Configuration38# --------------------------------------------------------------------------------------39def load_config(path: str | Path | None = None, overrides: dict[str, Any] | None = None) -> dict:40 """Load configs/baseline.yaml (or `path`) and apply dotted-key overrides."""41 path = Path(path) if path is not None else CONFIG_DIR / "baseline.yaml"42 with open(path) as fh:43 cfg = yaml.safe_load(fh)44 cfg = copy.deepcopy(cfg)45 for key, val in (overrides or {}).items():46 node = cfg47 parts = key.split(".")48 for part in parts[:-1]:49 node = node[part]50 node[parts[-1]] = val51 return cfg525354def load_scenarios(path: str | Path | None = None) -> dict:55 path = Path(path) if path is not None else CONFIG_DIR / "scenarios.yaml"56 with open(path) as fh:57 return yaml.safe_load(fh)585960# --------------------------------------------------------------------------------------61# Income process62# --------------------------------------------------------------------------------------63def rouwenhorst(n: int, rho: float, sigma: float) -> tuple[np.ndarray, np.ndarray]:64 """Rouwenhorst (1995) discretisation of log z' = rho log z + sigma eps, eps ~ N(0,1).6566 Returns (grid of log z, transition matrix). The grid is symmetric around 0 with67 unconditional s.d. sigma / sqrt(1 - rho^2).68 """69 if n == 1:70 return np.zeros(1), np.ones((1, 1))71 q = (1.0 + rho) / 2.072 P = np.array([[q, 1 - q], [1 - q, q]])73 for m in range(3, n + 1):74 Pn = np.zeros((m, m))75 Pn[:-1, :-1] += q * P76 Pn[:-1, 1:] += (1 - q) * P77 Pn[1:, :-1] += (1 - q) * P78 Pn[1:, 1:] += q * P79 Pn[1:-1, :] /= 2.080 P = Pn81 sd = sigma / math.sqrt(1.0 - rho**2)82 psi = sd * math.sqrt(n - 1)83 grid = np.linspace(-psi, psi, n)84 return grid, P858687def stationary_dist(P: np.ndarray) -> np.ndarray:88 w, v = np.linalg.eig(P.T)89 k = int(np.argmin(np.abs(w - 1.0)))90 pi = np.real(v[:, k])91 pi = pi / pi.sum()92 return np.maximum(pi, 0.0) / np.maximum(pi, 0.0).sum()939495# --------------------------------------------------------------------------------------96# Parameters (period units)97# --------------------------------------------------------------------------------------98@dataclass99class Params:100 cfg: dict = field(repr=False)101 # demographics102 J: int = 12103 j_retire: int = 9104 j_gift: int = 7105 j_estate: int = 7106 # preferences107 beta: float = 0.8154108 gamma: float = 2.0109 alpha: float = 0.2110 chi: float = 0.0111 nu: float = 3.0112 x_shift: float = 0.3113 # income114 n_z: int = 5115 rho_z: float = 0.774116 sigma_z: float = 0.304117 z_grid: np.ndarray = field(default=None, repr=False)118 Pi_z: np.ndarray = field(default=None, repr=False)119 pi_z: np.ndarray = field(default=None, repr=False)120 Pi_g: np.ndarray = field(default=None, repr=False) # parent z -> child z_1121 ybar: np.ndarray = field(default=None, repr=False) # age profile (ages 1..J); retirement rows = pension/z122 rr: float = 0.5123 dep_premium: float = 0.0124 sigma_trans: float = 0.0125 # housing126 h_grid: np.ndarray = field(default=None, repr=False)127 s_min: float = 0.3128 s_max: float = 3.4129 m: float = 0.05130 kappa_b: float = 0.025131 kappa_s: float = 0.05132 kappa_f: float = 0.01133 delta_H: float = 0.049134 # mortgage135 N: int = 5136 theta: float = 0.2137 psi: float = 0.39138 stress: float = 0.0139 contract: str = "frm"140 spread: float = 0.0141 # rates142 r: float = 0.104143 R: float = 1.104144 pibar: float = 0.104145 i_nom: float = 0.219146 i_q: float = 0.0147 # rental / supply148 phi: float = 0.05149 sigma_seg: float = 1.0150 eta_rent: float = 1.0151 eta_wage: float = 1.0152 eps: float = 1.0153 eps_r: float = 1.0154 H_bar: float | None = None155 Hr_bar: float | None = None156 p_ref: float = 1.0157 tau_b: float = 0.0158 # grids159 a_grid: np.ndarray = field(default=None, repr=False)160 b_grid: np.ndarray = field(default=None, repr=False)161 x_grid: np.ndarray = field(default=None, repr=False)162 B_fracs: np.ndarray = field(default=None, repr=False)163 gift_fracs: np.ndarray = field(default=None, repr=False)164 K_estate: int = 5165 c_min: float = 0.01166 # simulation / equilibrium167 seed: int = 20261020168 n_dyn: int = 96000169 burn: int = 72170171 # ---- derived helpers -------------------------------------------------------------172 @property173 def n_h(self) -> int:174 return len(self.h_grid)175176 @property177 def n_a(self) -> int:178 return len(self.a_grid)179180 @property181 def n_b(self) -> int:182 return len(self.b_grid)183184 @property185 def n_x(self) -> int:186 return len(self.x_grid)187188 @property189 def period_years(self) -> int:190 return int(self.cfg["meta"]["period_years"])191192 def income(self, j: int, iz: int) -> float:193 """Real income at age j (1-based) in state iz (pension if retired; z frozen)."""194 return float(self.ybar[j - 1] * self.z_grid[iz])195196 def income_matrix(self) -> np.ndarray:197 """(J, n_z) matrix of real income by age and state."""198 return np.outer(self.ybar, self.z_grid)199200201def _to_period(x_annual: float, years: int) -> float:202 return (1.0 + x_annual) ** years - 1.0203204205def double_exp_grid(n: int, lo: float, hi: float, curv: float = 2.0) -> np.ndarray:206 """Double-exponential spacing: dense near `lo`."""207 u = np.linspace(0.0, 1.0, n)208 g = (np.exp(np.exp(curv * u) - 1.0) - 1.0)209 g = g / g[-1]210 return lo + (hi - lo) * g211212213def power_grid(n: int, lo: float, hi: float, power: float = 2.0) -> np.ndarray:214 """Polynomial spacing: dense near `lo`, spacing growing like u^(power-1)."""215 u = np.linspace(0.0, 1.0, n)216 return lo + (hi - lo) * u ** power217218219def build_params(cfg: dict) -> Params:220 T = int(cfg["meta"]["period_years"])221 d, pr, inc, hs, mg, rt, rl, sp, gv, gr, sm = (cfg[k] for k in222 ("demographics", "preferences", "income", "housing", "mortgage", "rates", "rental",223 "supply", "government", "grids", "simulation"))224 P = Params(cfg=cfg)225 P.J, P.j_retire, P.j_gift, P.j_estate = int(d["J"]), int(d["j_retire"]), int(d["j_gift"]), int(d["j_estate"])226 # preferences227 P.beta = float(pr["beta_annual"]) ** T228 P.gamma, P.alpha, P.chi, P.nu, P.x_shift = (float(pr[k]) for k in ("gamma", "alpha", "chi", "nu", "x_shift"))229 # income process: annual AR(1) -> 5-year AR(1) with the same unconditional variance230 rho_a, sig_a = float(inc["rho_annual"]), float(inc["sigma_annual"])231 P.n_z = int(inc["n_z"])232 P.rho_z = rho_a ** T233 var_u = sig_a**2 / (1.0 - rho_a**2)234 P.sigma_z = math.sqrt(var_u * (1.0 - P.rho_z**2))235 lz, Pi = rouwenhorst(P.n_z, P.rho_z, P.sigma_z)236 z = np.exp(lz)237 pi_z = stationary_dist(Pi)238 z = z / float(pi_z @ z) # normalise E[z] = 1 under the stationary distribution239 P.z_grid, P.Pi_z, P.pi_z = z, Pi, pi_z240 rho_g = float(inc.get("rho_intergenerational", 0.0))241 if rho_g > 0:242 # child's log z_1 = rho_g log z_parent + noise, with the same unconditional variance243 _, Pg = rouwenhorst(P.n_z, rho_g, math.sqrt(var_u * (1.0 - rho_g**2)))244 P.Pi_g = Pg245 else:246 P.Pi_g = np.tile(pi_z, (P.n_z, 1))247 prof = np.array(inc["age_profile"], dtype=float)248 assert len(prof) == P.j_retire - 1, "age_profile must cover ages 1..j_retire-1"249 prof = prof / prof.mean()250 P.rr = float(inc["replacement_rate"])251 ybar = np.zeros(P.J)252 ybar[: P.j_retire - 1] = prof253 ybar[P.j_retire - 1:] = P.rr * prof[-1]254 P.ybar = ybar255 P.dep_premium = _to_period(float(inc.get("dep_premium_annual", 0.0)), T)256 P.sigma_trans = float(inc.get("sigma_transitory", 0.0))257 # housing258 P.h_grid = np.array(hs["h_sizes"], dtype=float)259 P.s_min, P.s_max = float(hs["s_min"]), float(hs["s_max"])260 P.m = float(hs["m_annual"]) * T261 P.kappa_b, P.kappa_s, P.kappa_f = float(hs["kappa_b"]), float(hs["kappa_s"]), float(hs["kappa_f"])262 P.delta_H = 1.0 - (1.0 - float(hs["delta_H_annual"])) ** T263 # mortgage264 P.N = int(mg["N"])265 P.theta, P.psi = float(mg["theta"]), float(mg["psi"])266 P.contract = str(mg["contract"])267 P.spread = _to_period(float(mg.get("spread_annual", 0.0)), T)268 # rates (Fisher): 1 + i = R (1 + pibar)269 P.r = _to_period(float(rt["r_annual"]), T)270 P.R = 1.0 + P.r271 P.pibar = _to_period(float(rt["pi_annual"]), T)272 P.i_nom = P.R * (1.0 + P.pibar) - 1.0 + P.spread273 i_annual = (1.0 + P.i_nom) ** (1.0 / T) - 1.0274 P.i_q = (1.0 + i_annual + float(mg["stress_annual"])) ** T - 1.0275 P.stress = P.i_q - P.i_nom276 # rental, supply, government277 P.phi, P.sigma_seg = float(rl["phi"]), float(rl["sigma_seg"])278 P.eta_rent, P.eta_wage = float(rl["eta_rent"]), float(rl["eta_wage"])279 P.eps, P.eps_r = float(sp["eps"]), float(sp["eps_r"])280 P.H_bar = None if sp.get("H_bar") is None else float(sp["H_bar"])281 P.Hr_bar = None if sp.get("Hr_bar") is None else float(sp["Hr_bar"])282 P.p_ref = float(sp["p_ref"])283 P.tau_b = float(gv["tau_b"])284 # grids285 P.a_grid = power_grid(int(gr["n_a"]), 0.0, float(gr["a_max"]), power=2.2)286 p_cap = float(gr["p_cap_factor"]) * P.p_ref287 P.b_grid = np.linspace(0.0, (1.0 - min(P.theta, 0.95)) * p_cap * P.h_grid[-1], int(gr["n_b"]))288 P.x_grid = double_exp_grid(int(gr["n_x"]), 0.0, float(gr["x_max"]), curv=1.6)289 P.x_grid[0] = 0.0290 P.B_fracs = np.linspace(0.0, 1.0, int(gr["n_B"]))291 P.gift_fracs = np.linspace(0.0, float(gr["gift_frac_max"]), int(gr["n_gift"]))292 P.K_estate = int(gr["K_estate"])293 P.c_min = float(gr["c_min"])294 # simulation295 P.seed = int(cfg["meta"]["seed"])296 P.n_dyn = int(sm["n_dynasties"])297 P.burn = int(sm["burn_periods"])298 return P299300301# --------------------------------------------------------------------------------------302# Prices object (per period; stationary or along a transition path)303# --------------------------------------------------------------------------------------304@dataclass305class Prices:306 """Real prices faced by households in one period.307308 p, p_r house and rental-unit prices; rho rent; R_a realised gross real deposit return;309 R_next gross real rate for (t, t+1) (discounting / landlord arbitrage);310 i_new nominal rate on loans originated this period; i_q qualifying rate;311 pi_next inflation between this period and the next (deflates balances);312 pi_u realised surprise in this period (affects wages / rents of sitting tenants);313 i_c_by_n nominal rate applying to an existing FRM with n remaining payments (array, index n-1);314 reset True if existing loans are repriced to i_new this period (term5/arm);315 erosion multiplicative factor applied to the surprise-period balance (1 = full erosion kept,316 (1+pi_u) = ARM within-period repricing undoes it).317 """318 p: float319 p_r: float320 rho: float321 R_a: float322 R_next: float323 i_new: float324 i_q: float325 pi_next: float326 pi_u: float = 0.0327 i_c_by_n: np.ndarray | None = None328 reset: bool = False329 wage_factor: float = 1.0330 rent_rigid_factor: float = 1.0331 p_next: float | None = None # price at which the estate house is valued (default p)332 transfer_rebate: float = 0.0 # lump-sum rebate to age-1 households (transfer tax)333334 def i_c(self, n: int) -> float:335 if self.reset or self.i_c_by_n is None:336 return self.i_new337 return float(self.i_c_by_n[n - 1])338339340def stationary_prices(P: Params, p: float, p_r: float, rebate: float = 0.0) -> Prices:341 rho = rent_from_price(P, p_r, p_r, P.R)342 return Prices(p=p, p_r=p_r, rho=rho, R_a=P.R * (1.0 + P.dep_premium), R_next=P.R,343 i_new=P.i_nom, i_q=P.i_q, pi_next=P.pibar, pi_u=0.0,344 i_c_by_n=np.full(P.N, P.i_nom), reset=False, p_next=p, transfer_rebate=rebate)345346347# --------------------------------------------------------------------------------------348# Utility349# --------------------------------------------------------------------------------------350def utility(P: Params, c: np.ndarray | float, s: np.ndarray | float) -> np.ndarray | float:351 comp = np.power(c, 1.0 - P.alpha) * np.power(s, P.alpha)352 return np.power(comp, 1.0 - P.gamma) / (1.0 - P.gamma)353354355def warm_glow(P: Params, x: np.ndarray | float) -> np.ndarray | float:356 if P.nu == 0.0:357 return 0.0 * np.asarray(x, dtype=float)358 return P.nu * np.power(np.asarray(x, dtype=float) + P.x_shift, 1.0 - P.gamma) / (1.0 - P.gamma)359360361# --------------------------------------------------------------------------------------362# Mortgage arithmetic363# --------------------------------------------------------------------------------------364def annuity_factor(i: float, n: int) -> float:365 """Payment per unit of balance for n remaining payments at rate i (annuity-due: payment at366 the start of each period)."""367 if n <= 0:368 return 0.0369 if abs(i) < 1e-12:370 return 1.0 / n371 return i / ((1.0 + i) * (1.0 - (1.0 + i) ** (-n)))372373374def payment(b: float, n: int, i: float) -> float:375 return b * annuity_factor(i, n)376377378def next_balance(b: float, n: int, i: float, pi_next: float) -> float:379 """Real balance next period after this period's payment (annuity schedule)."""380 pay = payment(b, n, i)381 return max(b - pay, 0.0) * (1.0 + i) / (1.0 + pi_next)382383384def ltv_cap(P: Params, p: float, h: float) -> float:385 return (1.0 - P.theta) * p * h386387388def dsti_cap(P: Params, y: float, i_q: float) -> float:389 """Largest balance such that the qualifying payment / income <= psi."""390 af = annuity_factor(i_q, P.N)391 return P.psi * y / af if af > 0 else np.inf392393394def max_loan(P: Params, p: float, h: float, y: float, i_q: float) -> float:395 return max(min(ltv_cap(P, p, h), dsti_cap(P, y, i_q)), 0.0)396397398# --------------------------------------------------------------------------------------399# Rental arbitrage and construction400# --------------------------------------------------------------------------------------401def rent_from_price(P: Params, p_r: float, p_r_next: float, R_next: float) -> float:402 """Landlord zero-profit: rho = p_r (1 + m + phi) - E p_r' / R."""403 return p_r * (1.0 + P.m + P.phi) - p_r_next / R_next404405406def user_cost_owner(P: Params, p: float) -> float:407 """Stationary owner user cost per unit (interest + maintenance), for reporting only."""408 return p * (P.r + P.m)409410411def supply_stock(P: Params, p: float, scale: float, eps: float) -> float:412 """Stationary stock S(p) = scale * (p / p_ref)^eps."""413 return scale * (p / P.p_ref) ** eps414415416def construction(P: Params, p: float, scale: float, eps: float) -> float:417 """Flow I_t = delta_H * scale * (p/p_ref)^eps so that the stationary stock is S(p)."""418 return P.delta_H * supply_stock(P, p, scale, eps)419420421def stock_next(P: Params, H: float, p: float, scale: float, eps: float) -> float:422 return (1.0 - P.delta_H) * H + construction(P, p, scale, eps)423424425# --------------------------------------------------------------------------------------426# Contract-type helpers for transitions427# --------------------------------------------------------------------------------------428CONTRACTS = ("frm", "term5", "arm")429430431def contract_flags(contract: str) -> tuple[bool, bool]:432 """(reset_each_period, within_period_repricing)."""433 if contract == "frm":434 return False, False435 if contract == "term5":436 return True, False437 if contract == "arm":438 return True, True439 raise ValueError(contract)440441442def parameter_table(P: Params) -> list[dict]:443 """Flat list of parameters (period units) with sources, for results/data/parameters.csv."""444 cfg = P.cfg445 rows = [446 ("J", P.J, "periods", "spec B.1"),447 ("j_retire", P.j_retire, "age index", "spec B.1"),448 ("j_gift = j_estate", P.j_gift, "age index", "spec B.1"),449 ("beta (period)", P.beta, "-", f"{cfg['preferences']['beta_annual']} annual, illustrative"),450 ("gamma", P.gamma, "-", "illustrative"),451 ("alpha", P.alpha, "-", "illustrative"),452 ("chi", P.chi, "-", "spec B.2 baseline 0"),453 ("nu", P.nu, "-", "illustrative"),454 ("x_shift", P.x_shift, "income units", "illustrative"),455 ("rho_z (period)", P.rho_z, "-", f"{cfg['income']['rho_annual']}^5, illustrative"),456 ("sigma_z (period)", P.sigma_z, "-", "same unconditional variance as annual, illustrative"),457 ("n_z", P.n_z, "-", "Rouwenhorst"),458 ("replacement rate", P.rr, "-", "illustrative"),459 ("sigma_transitory (period, log)", P.sigma_trans, "-", "illustrative; unanticipated, breaks cash atoms"),460 ("h sizes", list(map(float, P.h_grid)), "housing units", "spec B.12"),461 ("s_min", P.s_min, "housing units", "illustrative"),462 ("m (period)", P.m, "share of value", f"{cfg['housing']['m_annual']} annual, illustrative"),463 ("kappa_b", P.kappa_b, "share of price", "illustrative"),464 ("kappa_s", P.kappa_s, "share of price", "illustrative"),465 ("kappa_f", P.kappa_f, "share of value", "illustrative"),466 ("delta_H (period)", P.delta_H, "-", f"{cfg['housing']['delta_H_annual']} annual, illustrative"),467 ("N", P.N, "periods", "25-year amortisation (literature_review §7)"),468 ("theta", P.theta, "-", "illustrative (uninsured threshold)"),469 ("psi", P.psi, "-", "illustrative"),470 ("stress add-on (period)", P.stress, "-", f"{cfg['mortgage']['stress_annual']} annual, illustrative"),471 ("r (period)", P.r, "-", f"{cfg['rates']['r_annual']} annual, illustrative"),472 ("pibar (period)", P.pibar, "-", f"{cfg['rates']['pi_annual']} annual (2 % target)"),473 ("i (period, Fisher)", P.i_nom, "-", "derived"),474 ("phi (landlord wedge)", P.phi, "share of value", "illustrative"),475 ("sigma_seg", P.sigma_seg, "-", "spec B.6"),476 ("eps", P.eps, "-", "ILLUSTRATIVE (no verified Canadian estimate)"),477 ("eps_r", P.eps_r, "-", "illustrative"),478 ("tau_b", P.tau_b, "-", "policy experiment"),479 ("n_a / a_max", f"{P.n_a} / {P.a_grid[-1]:.1f}", "-", "grid"),480 ("n_b", P.n_b, "-", "grid"),481 ("n_x", P.n_x, "-", "grid"),482 ("n_B", len(P.B_fracs), "-", "grid"),483 ("K_estate", P.K_estate, "-", "lottery"),484 ("n_dynasties", P.n_dyn, "-", "simulation"),485 ("seed", P.seed, "-", "simulation"),486 ]487 return [dict(parameter=a, value=b, unit=c, source=d) for a, b, c, d in rows]488