"""WP18 — model primitives for the Part B lifecycle/OLG housing economy. Everything is in real terms (units of the consumption good), one period = 5 years. This module holds: configuration loading, the parameter object (with the annual -> period conversions), the grids, the income process (Rouwenhorst), utility functions, mortgage arithmetic (constant nominal payment, real balance law of motion, LTV and DSTI caps) and the construction / rental-arbitrage relations. No optimisation happens here. Conventions (research/model_specification.md §0, §B): * p real price of one unit of owner-occupied housing; p_r real price of a rental unit; * rho real rent per unit of rental housing services per period; * b real outstanding mortgage balance at the START of the period (before interest); * n remaining number of payments (1..N); b == 0 means no mortgage whatever n is; * pay real payment made at the START of the period (annuity-due on the outstanding balance): pay = b * i_c / ((1+i_c) (1 - (1+i_c)^-n)), exact for a constant-nominal-payment loan; start-of-period timing makes borrowing and depositing present-value consistent (borrowing B and paying pay now leaves (B - pay) R real next period, like a deposit); * b' = (b - pay) (1+i_c) / (1 + pi') real balance next period. """ from __future__ import annotations import copy import math from dataclasses import dataclass, field from pathlib import Path from typing import Any import numpy as np import yaml ROOT = Path(__file__).resolve().parents[1] CONFIG_DIR = ROOT / "configs" RESULTS_DIR = ROOT / "results" # -------------------------------------------------------------------------------------- # Configuration # -------------------------------------------------------------------------------------- def load_config(path: str | Path | None = None, overrides: dict[str, Any] | None = None) -> dict: """Load configs/baseline.yaml (or `path`) and apply dotted-key overrides.""" path = Path(path) if path is not None else CONFIG_DIR / "baseline.yaml" with open(path) as fh: cfg = yaml.safe_load(fh) cfg = copy.deepcopy(cfg) for key, val in (overrides or {}).items(): node = cfg parts = key.split(".") for part in parts[:-1]: node = node[part] node[parts[-1]] = val return cfg def load_scenarios(path: str | Path | None = None) -> dict: path = Path(path) if path is not None else CONFIG_DIR / "scenarios.yaml" with open(path) as fh: return yaml.safe_load(fh) # -------------------------------------------------------------------------------------- # Income process # -------------------------------------------------------------------------------------- def rouwenhorst(n: int, rho: float, sigma: float) -> tuple[np.ndarray, np.ndarray]: """Rouwenhorst (1995) discretisation of log z' = rho log z + sigma eps, eps ~ N(0,1). Returns (grid of log z, transition matrix). The grid is symmetric around 0 with unconditional s.d. sigma / sqrt(1 - rho^2). """ if n == 1: return np.zeros(1), np.ones((1, 1)) q = (1.0 + rho) / 2.0 P = np.array([[q, 1 - q], [1 - q, q]]) for m in range(3, n + 1): Pn = np.zeros((m, m)) Pn[:-1, :-1] += q * P Pn[:-1, 1:] += (1 - q) * P Pn[1:, :-1] += (1 - q) * P Pn[1:, 1:] += q * P Pn[1:-1, :] /= 2.0 P = Pn sd = sigma / math.sqrt(1.0 - rho**2) psi = sd * math.sqrt(n - 1) grid = np.linspace(-psi, psi, n) return grid, P def stationary_dist(P: np.ndarray) -> np.ndarray: w, v = np.linalg.eig(P.T) k = int(np.argmin(np.abs(w - 1.0))) pi = np.real(v[:, k]) pi = pi / pi.sum() return np.maximum(pi, 0.0) / np.maximum(pi, 0.0).sum() # -------------------------------------------------------------------------------------- # Parameters (period units) # -------------------------------------------------------------------------------------- @dataclass class Params: cfg: dict = field(repr=False) # demographics J: int = 12 j_retire: int = 9 j_gift: int = 7 j_estate: int = 7 # preferences beta: float = 0.8154 gamma: float = 2.0 alpha: float = 0.2 chi: float = 0.0 nu: float = 3.0 x_shift: float = 0.3 # income n_z: int = 5 rho_z: float = 0.774 sigma_z: float = 0.304 z_grid: np.ndarray = field(default=None, repr=False) Pi_z: np.ndarray = field(default=None, repr=False) pi_z: np.ndarray = field(default=None, repr=False) Pi_g: np.ndarray = field(default=None, repr=False) # parent z -> child z_1 ybar: np.ndarray = field(default=None, repr=False) # age profile (ages 1..J); retirement rows = pension/z rr: float = 0.5 dep_premium: float = 0.0 sigma_trans: float = 0.0 # housing h_grid: np.ndarray = field(default=None, repr=False) s_min: float = 0.3 s_max: float = 3.4 m: float = 0.05 kappa_b: float = 0.025 kappa_s: float = 0.05 kappa_f: float = 0.01 delta_H: float = 0.049 # mortgage N: int = 5 theta: float = 0.2 psi: float = 0.39 stress: float = 0.0 contract: str = "frm" spread: float = 0.0 # rates r: float = 0.104 R: float = 1.104 pibar: float = 0.104 i_nom: float = 0.219 i_q: float = 0.0 # rental / supply phi: float = 0.05 sigma_seg: float = 1.0 eta_rent: float = 1.0 eta_wage: float = 1.0 eps: float = 1.0 eps_r: float = 1.0 H_bar: float | None = None Hr_bar: float | None = None p_ref: float = 1.0 tau_b: float = 0.0 # grids a_grid: np.ndarray = field(default=None, repr=False) b_grid: np.ndarray = field(default=None, repr=False) x_grid: np.ndarray = field(default=None, repr=False) B_fracs: np.ndarray = field(default=None, repr=False) gift_fracs: np.ndarray = field(default=None, repr=False) K_estate: int = 5 c_min: float = 0.01 # simulation / equilibrium seed: int = 20261020 n_dyn: int = 96000 burn: int = 72 # ---- derived helpers ------------------------------------------------------------- @property def n_h(self) -> int: return len(self.h_grid) @property def n_a(self) -> int: return len(self.a_grid) @property def n_b(self) -> int: return len(self.b_grid) @property def n_x(self) -> int: return len(self.x_grid) @property def period_years(self) -> int: return int(self.cfg["meta"]["period_years"]) def income(self, j: int, iz: int) -> float: """Real income at age j (1-based) in state iz (pension if retired; z frozen).""" return float(self.ybar[j - 1] * self.z_grid[iz]) def income_matrix(self) -> np.ndarray: """(J, n_z) matrix of real income by age and state.""" return np.outer(self.ybar, self.z_grid) def _to_period(x_annual: float, years: int) -> float: return (1.0 + x_annual) ** years - 1.0 def double_exp_grid(n: int, lo: float, hi: float, curv: float = 2.0) -> np.ndarray: """Double-exponential spacing: dense near `lo`.""" u = np.linspace(0.0, 1.0, n) g = (np.exp(np.exp(curv * u) - 1.0) - 1.0) g = g / g[-1] return lo + (hi - lo) * g def power_grid(n: int, lo: float, hi: float, power: float = 2.0) -> np.ndarray: """Polynomial spacing: dense near `lo`, spacing growing like u^(power-1).""" u = np.linspace(0.0, 1.0, n) return lo + (hi - lo) * u ** power def build_params(cfg: dict) -> Params: T = int(cfg["meta"]["period_years"]) d, pr, inc, hs, mg, rt, rl, sp, gv, gr, sm = (cfg[k] for k in ("demographics", "preferences", "income", "housing", "mortgage", "rates", "rental", "supply", "government", "grids", "simulation")) P = Params(cfg=cfg) 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"]) # preferences P.beta = float(pr["beta_annual"]) ** T P.gamma, P.alpha, P.chi, P.nu, P.x_shift = (float(pr[k]) for k in ("gamma", "alpha", "chi", "nu", "x_shift")) # income process: annual AR(1) -> 5-year AR(1) with the same unconditional variance rho_a, sig_a = float(inc["rho_annual"]), float(inc["sigma_annual"]) P.n_z = int(inc["n_z"]) P.rho_z = rho_a ** T var_u = sig_a**2 / (1.0 - rho_a**2) P.sigma_z = math.sqrt(var_u * (1.0 - P.rho_z**2)) lz, Pi = rouwenhorst(P.n_z, P.rho_z, P.sigma_z) z = np.exp(lz) pi_z = stationary_dist(Pi) z = z / float(pi_z @ z) # normalise E[z] = 1 under the stationary distribution P.z_grid, P.Pi_z, P.pi_z = z, Pi, pi_z rho_g = float(inc.get("rho_intergenerational", 0.0)) if rho_g > 0: # child's log z_1 = rho_g log z_parent + noise, with the same unconditional variance _, Pg = rouwenhorst(P.n_z, rho_g, math.sqrt(var_u * (1.0 - rho_g**2))) P.Pi_g = Pg else: P.Pi_g = np.tile(pi_z, (P.n_z, 1)) prof = np.array(inc["age_profile"], dtype=float) assert len(prof) == P.j_retire - 1, "age_profile must cover ages 1..j_retire-1" prof = prof / prof.mean() P.rr = float(inc["replacement_rate"]) ybar = np.zeros(P.J) ybar[: P.j_retire - 1] = prof ybar[P.j_retire - 1:] = P.rr * prof[-1] P.ybar = ybar P.dep_premium = _to_period(float(inc.get("dep_premium_annual", 0.0)), T) P.sigma_trans = float(inc.get("sigma_transitory", 0.0)) # housing P.h_grid = np.array(hs["h_sizes"], dtype=float) P.s_min, P.s_max = float(hs["s_min"]), float(hs["s_max"]) P.m = float(hs["m_annual"]) * T P.kappa_b, P.kappa_s, P.kappa_f = float(hs["kappa_b"]), float(hs["kappa_s"]), float(hs["kappa_f"]) P.delta_H = 1.0 - (1.0 - float(hs["delta_H_annual"])) ** T # mortgage P.N = int(mg["N"]) P.theta, P.psi = float(mg["theta"]), float(mg["psi"]) P.contract = str(mg["contract"]) P.spread = _to_period(float(mg.get("spread_annual", 0.0)), T) # rates (Fisher): 1 + i = R (1 + pibar) P.r = _to_period(float(rt["r_annual"]), T) P.R = 1.0 + P.r P.pibar = _to_period(float(rt["pi_annual"]), T) P.i_nom = P.R * (1.0 + P.pibar) - 1.0 + P.spread i_annual = (1.0 + P.i_nom) ** (1.0 / T) - 1.0 P.i_q = (1.0 + i_annual + float(mg["stress_annual"])) ** T - 1.0 P.stress = P.i_q - P.i_nom # rental, supply, government P.phi, P.sigma_seg = float(rl["phi"]), float(rl["sigma_seg"]) P.eta_rent, P.eta_wage = float(rl["eta_rent"]), float(rl["eta_wage"]) P.eps, P.eps_r = float(sp["eps"]), float(sp["eps_r"]) P.H_bar = None if sp.get("H_bar") is None else float(sp["H_bar"]) P.Hr_bar = None if sp.get("Hr_bar") is None else float(sp["Hr_bar"]) P.p_ref = float(sp["p_ref"]) P.tau_b = float(gv["tau_b"]) # grids P.a_grid = power_grid(int(gr["n_a"]), 0.0, float(gr["a_max"]), power=2.2) p_cap = float(gr["p_cap_factor"]) * P.p_ref P.b_grid = np.linspace(0.0, (1.0 - min(P.theta, 0.95)) * p_cap * P.h_grid[-1], int(gr["n_b"])) P.x_grid = double_exp_grid(int(gr["n_x"]), 0.0, float(gr["x_max"]), curv=1.6) P.x_grid[0] = 0.0 P.B_fracs = np.linspace(0.0, 1.0, int(gr["n_B"])) P.gift_fracs = np.linspace(0.0, float(gr["gift_frac_max"]), int(gr["n_gift"])) P.K_estate = int(gr["K_estate"]) P.c_min = float(gr["c_min"]) # simulation P.seed = int(cfg["meta"]["seed"]) P.n_dyn = int(sm["n_dynasties"]) P.burn = int(sm["burn_periods"]) return P # -------------------------------------------------------------------------------------- # Prices object (per period; stationary or along a transition path) # -------------------------------------------------------------------------------------- @dataclass class Prices: """Real prices faced by households in one period. p, p_r house and rental-unit prices; rho rent; R_a realised gross real deposit return; R_next gross real rate for (t, t+1) (discounting / landlord arbitrage); i_new nominal rate on loans originated this period; i_q qualifying rate; pi_next inflation between this period and the next (deflates balances); pi_u realised surprise in this period (affects wages / rents of sitting tenants); i_c_by_n nominal rate applying to an existing FRM with n remaining payments (array, index n-1); reset True if existing loans are repriced to i_new this period (term5/arm); erosion multiplicative factor applied to the surprise-period balance (1 = full erosion kept, (1+pi_u) = ARM within-period repricing undoes it). """ p: float p_r: float rho: float R_a: float R_next: float i_new: float i_q: float pi_next: float pi_u: float = 0.0 i_c_by_n: np.ndarray | None = None reset: bool = False wage_factor: float = 1.0 rent_rigid_factor: float = 1.0 p_next: float | None = None # price at which the estate house is valued (default p) transfer_rebate: float = 0.0 # lump-sum rebate to age-1 households (transfer tax) def i_c(self, n: int) -> float: if self.reset or self.i_c_by_n is None: return self.i_new return float(self.i_c_by_n[n - 1]) def stationary_prices(P: Params, p: float, p_r: float, rebate: float = 0.0) -> Prices: rho = rent_from_price(P, p_r, p_r, P.R) return Prices(p=p, p_r=p_r, rho=rho, R_a=P.R * (1.0 + P.dep_premium), R_next=P.R, i_new=P.i_nom, i_q=P.i_q, pi_next=P.pibar, pi_u=0.0, i_c_by_n=np.full(P.N, P.i_nom), reset=False, p_next=p, transfer_rebate=rebate) # -------------------------------------------------------------------------------------- # Utility # -------------------------------------------------------------------------------------- def utility(P: Params, c: np.ndarray | float, s: np.ndarray | float) -> np.ndarray | float: comp = np.power(c, 1.0 - P.alpha) * np.power(s, P.alpha) return np.power(comp, 1.0 - P.gamma) / (1.0 - P.gamma) def warm_glow(P: Params, x: np.ndarray | float) -> np.ndarray | float: if P.nu == 0.0: return 0.0 * np.asarray(x, dtype=float) return P.nu * np.power(np.asarray(x, dtype=float) + P.x_shift, 1.0 - P.gamma) / (1.0 - P.gamma) # -------------------------------------------------------------------------------------- # Mortgage arithmetic # -------------------------------------------------------------------------------------- def annuity_factor(i: float, n: int) -> float: """Payment per unit of balance for n remaining payments at rate i (annuity-due: payment at the start of each period).""" 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))) def payment(b: float, n: int, i: float) -> float: return b * annuity_factor(i, n) def next_balance(b: float, n: int, i: float, pi_next: float) -> float: """Real balance next period after this period's payment (annuity schedule).""" pay = payment(b, n, i) return max(b - pay, 0.0) * (1.0 + i) / (1.0 + pi_next) def ltv_cap(P: Params, p: float, h: float) -> float: return (1.0 - P.theta) * p * h def dsti_cap(P: Params, y: float, i_q: float) -> float: """Largest balance such that the qualifying payment / income <= psi.""" af = annuity_factor(i_q, P.N) return P.psi * y / af if af > 0 else np.inf def max_loan(P: Params, p: float, h: float, y: float, i_q: float) -> float: return max(min(ltv_cap(P, p, h), dsti_cap(P, y, i_q)), 0.0) # -------------------------------------------------------------------------------------- # Rental arbitrage and construction # -------------------------------------------------------------------------------------- def rent_from_price(P: Params, p_r: float, p_r_next: float, R_next: float) -> float: """Landlord zero-profit: rho = p_r (1 + m + phi) - E p_r' / R.""" return p_r * (1.0 + P.m + P.phi) - p_r_next / R_next def user_cost_owner(P: Params, p: float) -> float: """Stationary owner user cost per unit (interest + maintenance), for reporting only.""" return p * (P.r + P.m) def supply_stock(P: Params, p: float, scale: float, eps: float) -> float: """Stationary stock S(p) = scale * (p / p_ref)^eps.""" return scale * (p / P.p_ref) ** eps def construction(P: Params, p: float, scale: float, eps: float) -> float: """Flow I_t = delta_H * scale * (p/p_ref)^eps so that the stationary stock is S(p).""" return P.delta_H * supply_stock(P, p, scale, eps) def stock_next(P: Params, H: float, p: float, scale: float, eps: float) -> float: return (1.0 - P.delta_H) * H + construction(P, p, scale, eps) # -------------------------------------------------------------------------------------- # Contract-type helpers for transitions # -------------------------------------------------------------------------------------- CONTRACTS = ("frm", "term5", "arm") def contract_flags(contract: str) -> tuple[bool, bool]: """(reset_each_period, within_period_repricing).""" if contract == "frm": return False, False if contract == "term5": return True, False if contract == "arm": return True, True raise ValueError(contract) def parameter_table(P: Params) -> list[dict]: """Flat list of parameters (period units) with sources, for results/data/parameters.csv.""" cfg = P.cfg rows = [ ("J", P.J, "periods", "spec B.1"), ("j_retire", P.j_retire, "age index", "spec B.1"), ("j_gift = j_estate", P.j_gift, "age index", "spec B.1"), ("beta (period)", P.beta, "-", f"{cfg['preferences']['beta_annual']} annual, illustrative"), ("gamma", P.gamma, "-", "illustrative"), ("alpha", P.alpha, "-", "illustrative"), ("chi", P.chi, "-", "spec B.2 baseline 0"), ("nu", P.nu, "-", "illustrative"), ("x_shift", P.x_shift, "income units", "illustrative"), ("rho_z (period)", P.rho_z, "-", f"{cfg['income']['rho_annual']}^5, illustrative"), ("sigma_z (period)", P.sigma_z, "-", "same unconditional variance as annual, illustrative"), ("n_z", P.n_z, "-", "Rouwenhorst"), ("replacement rate", P.rr, "-", "illustrative"), ("sigma_transitory (period, log)", P.sigma_trans, "-", "illustrative; unanticipated, breaks cash atoms"), ("h sizes", list(map(float, P.h_grid)), "housing units", "spec B.12"), ("s_min", P.s_min, "housing units", "illustrative"), ("m (period)", P.m, "share of value", f"{cfg['housing']['m_annual']} annual, illustrative"), ("kappa_b", P.kappa_b, "share of price", "illustrative"), ("kappa_s", P.kappa_s, "share of price", "illustrative"), ("kappa_f", P.kappa_f, "share of value", "illustrative"), ("delta_H (period)", P.delta_H, "-", f"{cfg['housing']['delta_H_annual']} annual, illustrative"), ("N", P.N, "periods", "25-year amortisation (literature_review §7)"), ("theta", P.theta, "-", "illustrative (uninsured threshold)"), ("psi", P.psi, "-", "illustrative"), ("stress add-on (period)", P.stress, "-", f"{cfg['mortgage']['stress_annual']} annual, illustrative"), ("r (period)", P.r, "-", f"{cfg['rates']['r_annual']} annual, illustrative"), ("pibar (period)", P.pibar, "-", f"{cfg['rates']['pi_annual']} annual (2 % target)"), ("i (period, Fisher)", P.i_nom, "-", "derived"), ("phi (landlord wedge)", P.phi, "share of value", "illustrative"), ("sigma_seg", P.sigma_seg, "-", "spec B.6"), ("eps", P.eps, "-", "ILLUSTRATIVE (no verified Canadian estimate)"), ("eps_r", P.eps_r, "-", "illustrative"), ("tau_b", P.tau_b, "-", "policy experiment"), ("n_a / a_max", f"{P.n_a} / {P.a_grid[-1]:.1f}", "-", "grid"), ("n_b", P.n_b, "-", "grid"), ("n_x", P.n_x, "-", "grid"), ("n_B", len(P.B_fracs), "-", "grid"), ("K_estate", P.K_estate, "-", "lottery"), ("n_dynasties", P.n_dyn, "-", "simulation"), ("seed", P.seed, "-", "simulation"), ] return [dict(parameter=a, value=b, unit=c, source=d) for a, b, c, d in rows]