1"""Primitives of the WP19 search model: match-value distributions, seller-type grids, parameters.23All quantities are in units of the mean buyer valuation of the quality segment (q = 1), time in years.4"""5from __future__ import annotations67import hashlib8import json9from dataclasses import dataclass, field, asdict10from pathlib import Path11from typing import Optional1213import numpy as np14import yaml15from scipy import integrate, special, stats1617SQRT2PI = np.sqrt(2.0 * np.pi)181920# ----------------------------------------------------------------------------------------------21# Match-value distributions G (valuation w of a buyer for a given house, PV units)22# ----------------------------------------------------------------------------------------------23class MatchDist:24 """Interface: cdf, pdf, H(v)=E[(w-v)^+], M(v)=E[w-v|w>=v], cond_mean, cond_var, mean, support."""2526 lo: float27 hi: float28 name: str = "base"2930 def cdf(self, v):31 raise NotImplementedError3233 def pdf(self, v):34 raise NotImplementedError3536 def H(self, v):37 raise NotImplementedError3839 def surv(self, v):40 return 1.0 - self.cdf(v)4142 def hazard(self, v):43 s = self.surv(v)44 return np.where(s > 0, self.pdf(v) / np.maximum(s, 1e-300), np.inf)4546 def M(self, v):47 s = self.surv(v)48 return np.where(s > 0, self.H(v) / np.maximum(s, 1e-300), 0.0)4950 def cond_mean(self, v):51 return v + self.M(v)5253 def cond_var(self, v):54 """Var(w | w >= v) by quadrature (generic)."""55 v = float(v)56 s = float(self.surv(v))57 if s <= 0:58 return 0.059 hi = self.hi if np.isfinite(self.hi) else v + 50.0 * max(1.0, self.M(v))60 m2 = integrate.quad(lambda w: w * w * float(self.pdf(w)), v, hi, limit=200)[0] / s61 m1 = float(self.cond_mean(v))62 return max(m2 - m1 * m1, 0.0)6364 def mean(self):65 return float(self.H(self.lo) + self.lo) if np.isfinite(self.lo) else float(self.H(-1e6) - 1e6)6667 def sample(self, rng: np.random.Generator, n: int):68 raise NotImplementedError6970 def is_ifr(self) -> Optional[bool]:71 return None727374@dataclass75class TruncNormal(MatchDist):76 mu: float = 1.077 sigma: float = 0.1078 lo: float = 0.7079 hi: float = 1.3080 name: str = "truncnorm"8182 def __post_init__(self):83 self.a = (self.lo - self.mu) / self.sigma84 self.b = (self.hi - self.mu) / self.sigma85 self.Z = stats.norm.cdf(self.b) - stats.norm.cdf(self.a)8687 def _z(self, v):88 return (np.asarray(v, dtype=float) - self.mu) / self.sigma8990 def cdf(self, v):91 z = np.clip(self._z(v), self.a, self.b)92 return (stats.norm.cdf(z) - stats.norm.cdf(self.a)) / self.Z9394 def pdf(self, v):95 z = self._z(v)96 inside = (z >= self.a) & (z <= self.b)97 return np.where(inside, stats.norm.pdf(z) / (self.sigma * self.Z), 0.0)9899 def H(self, v):100 """E[(w-v)^+] closed form for the truncated normal."""101 v = np.asarray(v, dtype=float)102 zv = np.clip(self._z(v), self.a, self.b)103 term = (self.sigma * (stats.norm.pdf(zv) - stats.norm.pdf(self.b))104 + (self.mu - np.maximum(v, self.lo)) * (stats.norm.cdf(self.b) - stats.norm.cdf(zv))) / self.Z105 # below the support: E[w] - v106 below = v < self.lo107 mean = self.mu + self.sigma * (stats.norm.pdf(self.a) - stats.norm.pdf(self.b)) / self.Z108 return np.where(below, mean - v, np.maximum(term, 0.0))109110 def mean(self):111 return float(self.mu + self.sigma * (stats.norm.pdf(self.a) - stats.norm.pdf(self.b)) / self.Z)112113 def sample(self, rng, n):114 return stats.truncnorm(self.a, self.b, loc=self.mu, scale=self.sigma).rvs(size=n, random_state=rng)115116 def is_ifr(self):117 return True # log-concave density118119120@dataclass121class Exponential(MatchDist):122 rate: float = 10.0123 lo: float = 0.70124 hi: float = np.inf125 name: str = "exponential"126127 def cdf(self, v):128 v = np.asarray(v, dtype=float)129 return np.where(v < self.lo, 0.0, 1.0 - np.exp(-self.rate * (v - self.lo)))130131 def pdf(self, v):132 v = np.asarray(v, dtype=float)133 return np.where(v < self.lo, 0.0, self.rate * np.exp(-self.rate * (v - self.lo)))134135 def H(self, v):136 v = np.asarray(v, dtype=float)137 return np.where(v < self.lo, self.lo + 1.0 / self.rate - v, np.exp(-self.rate * (v - self.lo)) / self.rate)138139 def mean(self):140 return self.lo + 1.0 / self.rate141142 def sample(self, rng, n):143 return self.lo + rng.exponential(1.0 / self.rate, size=n)144145 def is_ifr(self):146 return True # constant failure rate (boundary case)147148149@dataclass150class Pareto(MatchDist):151 a: float = 3.0152 lo: float = 0.70153 hi: float = np.inf154 name: str = "pareto"155156 def cdf(self, v):157 v = np.asarray(v, dtype=float)158 return np.where(v < self.lo, 0.0, 1.0 - (self.lo / np.maximum(v, self.lo)) ** self.a)159160 def pdf(self, v):161 v = np.asarray(v, dtype=float)162 return np.where(v < self.lo, 0.0, self.a * self.lo ** self.a / np.maximum(v, self.lo) ** (self.a + 1))163164 def H(self, v):165 v = np.asarray(v, dtype=float)166 mean = self.lo * self.a / (self.a - 1.0)167 above = (self.lo / np.maximum(v, self.lo)) ** self.a * np.maximum(v, self.lo) / (self.a - 1.0)168 return np.where(v < self.lo, mean - v, above)169170 def mean(self):171 return self.lo * self.a / (self.a - 1.0)172173 def sample(self, rng, n):174 return self.lo * (1.0 - rng.random(n)) ** (-1.0 / self.a)175176 def is_ifr(self):177 return False # decreasing failure rate178179180@dataclass181class Weibull(MatchDist):182 """w = lo + scale * X, X ~ Weibull(shape). IFR iff shape >= 1."""183 shape: float = 2.0184 scale: float = 0.2185 lo: float = 0.70186 hi: float = np.inf187 name: str = "weibull"188189 def _x(self, v):190 return np.maximum((np.asarray(v, dtype=float) - self.lo) / self.scale, 0.0)191192 def cdf(self, v):193 return 1.0 - np.exp(-self._x(v) ** self.shape)194195 def pdf(self, v):196 x = self._x(v)197 return np.where(np.asarray(v) >= self.lo,198 (self.shape / self.scale) * x ** (self.shape - 1.0) * np.exp(-x ** self.shape), 0.0)199200 def H(self, v):201 v = np.asarray(v, dtype=float)202 x = self._x(v)203 # E[(X - x)^+] for Weibull: scale * Gamma(1+1/k) * Q(1+1/k, x^k) - x * exp(-x^k), Q = upper regularized gamma204 k = self.shape205 ex = self.scale * (special.gamma(1.0 + 1.0 / k) * special.gammaincc(1.0 + 1.0 / k, x ** k) - x * np.exp(-x ** k))206 mean = self.lo + self.scale * special.gamma(1.0 + 1.0 / k)207 return np.where(v < self.lo, mean - v, np.maximum(ex, 0.0))208209 def mean(self):210 return self.lo + self.scale * special.gamma(1.0 + 1.0 / self.shape)211212 def sample(self, rng, n):213 return self.lo + self.scale * rng.weibull(self.shape, size=n)214215 def is_ifr(self):216 return self.shape >= 1.0217218219def make_dist(cfg: dict) -> MatchDist:220 d = dict(cfg)221 kind = d.pop("dist")222 if kind == "truncnorm":223 return TruncNormal(**d)224 if kind == "exponential":225 return Exponential(**d)226 if kind == "pareto":227 return Pareto(**d)228 if kind == "weibull":229 return Weibull(**d)230 raise ValueError(f"unknown match-value distribution {kind}")231232233# ----------------------------------------------------------------------------------------------234# Seller-type grid (discrete approximation of F)235# ----------------------------------------------------------------------------------------------236def seller_type_grid(c_min: float, c_max: float, n_types: int, beta_a: float, beta_b: float):237 """Midpoint discretisation of Beta(a,b) rescaled to [c_min, c_max]: returns (c, f) with f summing to one."""238 if n_types == 1:239 mean = c_min + (c_max - c_min) * beta_a / (beta_a + beta_b)240 return np.array([mean]), np.array([1.0])241 edges = np.linspace(0.0, 1.0, n_types + 1)242 mass = np.diff(stats.beta.cdf(edges, beta_a, beta_b))243 # conditional means within each bin (exact for Beta via incomplete beta ratio)244 cond = np.diff(stats.beta.cdf(edges, beta_a + 1.0, beta_b)) * (beta_a / (beta_a + beta_b)) / mass245 c = c_min + (c_max - c_min) * cond246 return c, mass / mass.sum()247248249# ----------------------------------------------------------------------------------------------250# Parameters251# ----------------------------------------------------------------------------------------------252@dataclass253class Params:254 r: float = 0.04255 s: float = 1.0256 delta: float = 0.5257 w0: float = 0.8258 k: float = 0.02259 kappa_b: float = 0.05260 phi: float = 0.85261 A: float = 4.5262 alpha: float = 0.5263 beta: float = 0.5264 c: np.ndarray = field(default_factory=lambda: np.array([0.05]))265 f: np.ndarray = field(default_factory=lambda: np.array([1.0]))266 G: MatchDist = field(default_factory=TruncNormal)267 # financing extension268 pi: float = 0.3269 nu: float = 12.0270 phi_m: float = 0.85271 # numerics272 theta_bounds: tuple = (0.02, 60.0)273 root_xtol: float = 1e-12274 grid_check_points: int = 60275 # dynamics276 T: float = 8.0277 dt: float = 0.02278 damping: float = 0.3279 dyn_tol: float = 1e-7280 dyn_max_iter: int = 400281 name: str = "baseline"282 seed: int = 20261020283 raw: dict = field(default_factory=dict)284285 # --- derived -----------------------------------------------------------------------------286 def lam_s(self, theta):287 return self.A * np.asarray(theta, dtype=float) ** self.alpha288289 def lam_b(self, theta):290 return self.A * np.asarray(theta, dtype=float) ** (self.alpha - 1.0)291292 def copy(self, **changes) -> "Params":293 d = {k: v for k, v in self.__dict__.items()}294 d.update(changes)295 return Params(**d)296297 def with_types(self, c_min, c_max, n_types, beta_a, beta_b) -> "Params":298 c, f = seller_type_grid(c_min, c_max, n_types, beta_a, beta_b)299 return self.copy(c=c, f=f)300301 def config_hash(self) -> str:302 payload = {k: (v.tolist() if isinstance(v, np.ndarray) else (asdict(v) if hasattr(v, "__dataclass_fields__") else v))303 for k, v in self.__dict__.items() if k != "raw"}304 payload["G"] = {"name": self.G.name, **{k: float(v) for k, v in self.G.__dict__.items()305 if isinstance(v, (int, float)) and k not in ("Z", "a", "b")}}306 return hashlib.sha256(json.dumps(payload, sort_keys=True, default=str).encode()).hexdigest()[:16]307308309def load_params(path: str | Path) -> Params:310 cfg = yaml.safe_load(Path(path).read_text())311 c, f = seller_type_grid(cfg["sellers"]["c_min"], cfg["sellers"]["c_max"], cfg["sellers"]["n_types"],312 cfg["sellers"]["beta_a"], cfg["sellers"]["beta_b"])313 fin = cfg.get("financing_extension", {})314 num = cfg.get("numerics", {})315 dyn = cfg.get("dynamics", {})316 return Params(317 r=cfg["preferences"]["r"], s=cfg["sellers"]["s"], delta=cfg["sellers"]["delta"], w0=cfg["sellers"]["w0"],318 k=cfg["buyers"]["k"], kappa_b=cfg["buyers"]["kappa_b"], phi=cfg["buyers"]["phi"],319 A=cfg["matching"]["A"], alpha=cfg["matching"]["alpha"], beta=cfg["bargaining"]["beta"],320 c=c, f=f, G=make_dist(cfg["match_values"]),321 pi=fin.get("pi", 0.3), nu=fin.get("nu", 12.0), phi_m=fin.get("phi", cfg["buyers"]["phi"]),322 theta_bounds=tuple(num.get("theta_bounds", (0.02, 60.0))), root_xtol=float(num.get("root_xtol", 1e-12)),323 grid_check_points=int(num.get("grid_check_points", 60)),324 T=float(dyn.get("T", 8.0)), dt=float(dyn.get("dt", 0.02)), damping=float(dyn.get("damping", 0.3)),325 dyn_tol=float(dyn.get("tol", 1e-7)), dyn_max_iter=int(dyn.get("max_iter", 400)),326 name=cfg.get("name", "baseline"), seed=int(cfg.get("seed", 0)), raw=cfg,327 )328