1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Stylised buyer–seller–lender model of `paper/sections/welfare.tex` (numpy only).34All quantities are normalised by the market-consistent value V = 1.5 * Seller reservation S = 1 + u_S, buyer willingness to pay B = 1 + u_B; a match is mutually6 beneficial iff B >= S.7 * Appraisal A = 1 + b + s Z, Z ~ N(0, 1): bias b, standard deviation s.8 * Lender: loan L = ltv (1 - h(s)) A with haircut h(s) = h0 + h1 s; buyer cash W = (1 - ltv) kappa.9 * Transaction completes iff B >= S and S <= W + L; price P = min(S + theta (B - S), W + L).10 * Default with exogenous probability pd; recovery (1 - c) V1 with V1 = 1 + g, g ~ N(mu_g, sigma_g^2);11 loss given default max(L - (1 - c) V1, 0).12 * N lenders with AVM errors Z_n = sqrt(rho) C + sqrt(1 - rho) E_n (common-model correlation rho).13Only ltv (80 %) and the haircut rule are transparent calibrations; every other number is ILLUSTRATIVE.14"""15from __future__ import annotations1617from dataclasses import dataclass, replace1819import numpy as np202122@dataclass(frozen=True)23class WelfareParams:24 ltv: float = 0.80 # maximum loan-to-value (transparent: conventional uninsured mortgage cap in Canada)25 h0: float = 0.0 # haircut intercept26 h1: float = 0.0 # haircut slope in valuation uncertainty s27 kappa: float = 1.25 # buyer cash relative to the minimum down payment on V (illustrative)28 theta: float = 0.5 # seller's share of the bargaining surplus (illustrative)29 mu_s: float = -0.05 # seller reservation mean (illustrative)30 sd_s: float = 0.0531 mu_b: float = 0.05 # buyer WTP mean (illustrative)32 sd_b: float = 0.0833 pd: float = 0.02 # probability of default (illustrative)34 c_liq: float = 0.15 # liquidation cost share (illustrative)35 mu_g: float = -0.10 # value change at default, conditional on default (illustrative stress)36 sd_g: float = 0.103738 def haircut(self, s: float) -> float:39 return min(max(self.h0 + self.h1 * s, 0.0), 0.95)4041 def cash(self) -> float:42 return (1.0 - self.ltv) * self.kappa434445def _phi_cdf(x: np.ndarray) -> np.ndarray:46 """Standard normal CDF via the complementary error function (numpy only)."""47 from math import erf # scalar; vectorised below4849 x = np.asarray(x, dtype=float)50 return 0.5 * (1.0 + np.vectorize(erf)(x / np.sqrt(2.0)))515253def failure_prob_given_seller(u_s: np.ndarray, b: float, s: float, p: WelfareParams) -> np.ndarray:54 """P(W + L < S | u_S): the mutually beneficial transaction is not financeable (Proposition W1)."""55 coef = p.ltv * (1.0 - p.haircut(s))56 a = (1.0 + u_s) - p.cash() - coef * (1.0 + b)57 if s <= 0:58 return (a > 0).astype(float)59 return _phi_cdf(a / (coef * s))606162def simulate_market(b: float, s: float, p: WelfareParams, n: int, rng: np.random.Generator,63 shocks: dict | None = None) -> dict[str, float]:64 """Monte Carlo of the transaction stage. Returns failure rate among beneficial matches, completion65 rate, mean price, mean seller and buyer surplus (conditional on completion) and the exact failure66 rate obtained by integrating Z analytically given the seller draw."""67 if shocks is None:68 shocks = draw_market_shocks(n, rng)69 u_s, u_b, z = shocks["u_s"], shocks["u_b"], shocks["z"]70 S, B = 1.0 + u_s, 1.0 + u_b71 A = 1.0 + b + s * z72 L = p.ltv * (1.0 - p.haircut(s)) * A73 cap = p.cash() + L74 beneficial = B >= S75 completes = beneficial & (S <= cap)76 P = np.minimum(S + p.theta * (B - S), cap)77 fail_exact = failure_prob_given_seller(u_s[beneficial], b, s, p).mean() if beneficial.any() else np.nan78 out = {79 "beneficial_rate": float(beneficial.mean()),80 "failure_rate": float(1.0 - completes.sum() / max(beneficial.sum(), 1)),81 "failure_rate_exact": float(fail_exact),82 "completion_rate": float(completes.mean()),83 "mean_price": float(P[completes].mean()) if completes.any() else np.nan,84 "seller_surplus": float((P - S)[completes].mean()) if completes.any() else np.nan,85 "buyer_surplus": float((B - P)[completes].mean()) if completes.any() else np.nan,86 "total_surplus": float((B - S)[completes].mean()) if completes.any() else np.nan,87 "lost_surplus": float((B - S)[beneficial & ~completes].sum() / n),88 }89 return out909192def draw_market_shocks(n: int, rng: np.random.Generator, p: WelfareParams = WelfareParams()) -> dict:93 return {94 "u_s": p.mu_s + p.sd_s * rng.standard_normal(n),95 "u_b": p.mu_b + p.sd_b * rng.standard_normal(n),96 "z": rng.standard_normal(n),97 }9899100# Gauss–Hermite nodes for integrating the future-value shock at default (numpy only)101_GH_X, _GH_W = np.polynomial.hermite_e.hermegauss(40)102_GH_W = _GH_W / _GH_W.sum()103104105def expected_loss_given_z(z: np.ndarray, b: float, s: float, p: WelfareParams) -> np.ndarray:106 """g(z) = pd * E_{V1}[max(L(z) - (1 - c) V1, 0)], integrating V1 by Gauss–Hermite quadrature."""107 z = np.asarray(z, dtype=float)108 L = p.ltv * (1.0 - p.haircut(s)) * (1.0 + b + s * z)109 v1 = 1.0 + p.mu_g + p.sd_g * _GH_X # quadrature nodes for V1110 loss = np.maximum(L[:, None] - (1.0 - p.c_liq) * v1[None, :], 0.0)111 return p.pd * loss @ _GH_W112113114def expected_loss(b: float, s: float, p: WelfareParams) -> float:115 """Unconditional expected loss E[g(Z)] by Gauss–Hermite quadrature over Z (Proposition W2)."""116 return float(expected_loss_given_z(_GH_X, b, s, p) @ _GH_W)117118119def systemic_loss(rho: float, b: float, s: float, p: WelfareParams, n_lenders: int, n_draws: int,120 rng: np.random.Generator, shocks: dict | None = None) -> dict[str, float]:121 """System loss sum_n g(Z_n) with exchangeable AVM errors (Proposition W3). Common random numbers:122 pass `shocks` = {'C': (n_draws,), 'E': (n_draws, n_lenders)} to compare rho values smoothly."""123 if shocks is None:124 shocks = {"C": rng.standard_normal(n_draws), "E": rng.standard_normal((n_draws, n_lenders))}125 Z = np.sqrt(rho) * shocks["C"][:, None] + np.sqrt(1.0 - rho) * shocks["E"]126 g = expected_loss_given_z(Z.ravel(), b, s, p).reshape(Z.shape)127 total = g.sum(axis=1)128 q95 = np.quantile(total, 0.95)129 return {130 "mean": float(total.mean()),131 "sd": float(total.std(ddof=1)),132 "es95": float(total[total >= q95].mean()),133 "cov_pair": float(np.cov(g[:, 0], g[:, 1])[0, 1]) if n_lenders > 1 else np.nan,134 }135136137def with_params(p: WelfareParams, **kw) -> WelfareParams:138 return replace(p, **kw)139