# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca """Stylised buyer–seller–lender model of `paper/sections/welfare.tex` (numpy only). All quantities are normalised by the market-consistent value V = 1. * Seller reservation S = 1 + u_S, buyer willingness to pay B = 1 + u_B; a match is mutually beneficial iff B >= S. * Appraisal A = 1 + b + s Z, Z ~ N(0, 1): bias b, standard deviation s. * Lender: loan L = ltv (1 - h(s)) A with haircut h(s) = h0 + h1 s; buyer cash W = (1 - ltv) kappa. * Transaction completes iff B >= S and S <= W + L; price P = min(S + theta (B - S), W + L). * Default with exogenous probability pd; recovery (1 - c) V1 with V1 = 1 + g, g ~ N(mu_g, sigma_g^2); loss given default max(L - (1 - c) V1, 0). * N lenders with AVM errors Z_n = sqrt(rho) C + sqrt(1 - rho) E_n (common-model correlation rho). Only ltv (80 %) and the haircut rule are transparent calibrations; every other number is ILLUSTRATIVE. """ from __future__ import annotations from dataclasses import dataclass, replace import numpy as np @dataclass(frozen=True) class WelfareParams: ltv: float = 0.80 # maximum loan-to-value (transparent: conventional uninsured mortgage cap in Canada) h0: float = 0.0 # haircut intercept h1: float = 0.0 # haircut slope in valuation uncertainty s kappa: float = 1.25 # buyer cash relative to the minimum down payment on V (illustrative) theta: float = 0.5 # seller's share of the bargaining surplus (illustrative) mu_s: float = -0.05 # seller reservation mean (illustrative) sd_s: float = 0.05 mu_b: float = 0.05 # buyer WTP mean (illustrative) sd_b: float = 0.08 pd: float = 0.02 # probability of default (illustrative) c_liq: float = 0.15 # liquidation cost share (illustrative) mu_g: float = -0.10 # value change at default, conditional on default (illustrative stress) sd_g: float = 0.10 def haircut(self, s: float) -> float: return min(max(self.h0 + self.h1 * s, 0.0), 0.95) def cash(self) -> float: return (1.0 - self.ltv) * self.kappa def _phi_cdf(x: np.ndarray) -> np.ndarray: """Standard normal CDF via the complementary error function (numpy only).""" from math import erf # scalar; vectorised below x = np.asarray(x, dtype=float) return 0.5 * (1.0 + np.vectorize(erf)(x / np.sqrt(2.0))) def failure_prob_given_seller(u_s: np.ndarray, b: float, s: float, p: WelfareParams) -> np.ndarray: """P(W + L < S | u_S): the mutually beneficial transaction is not financeable (Proposition W1).""" coef = p.ltv * (1.0 - p.haircut(s)) a = (1.0 + u_s) - p.cash() - coef * (1.0 + b) if s <= 0: return (a > 0).astype(float) return _phi_cdf(a / (coef * s)) def simulate_market(b: float, s: float, p: WelfareParams, n: int, rng: np.random.Generator, shocks: dict | None = None) -> dict[str, float]: """Monte Carlo of the transaction stage. Returns failure rate among beneficial matches, completion rate, mean price, mean seller and buyer surplus (conditional on completion) and the exact failure rate obtained by integrating Z analytically given the seller draw.""" if shocks is None: shocks = draw_market_shocks(n, rng) u_s, u_b, z = shocks["u_s"], shocks["u_b"], shocks["z"] S, B = 1.0 + u_s, 1.0 + u_b A = 1.0 + b + s * z L = p.ltv * (1.0 - p.haircut(s)) * A cap = p.cash() + L beneficial = B >= S completes = beneficial & (S <= cap) P = np.minimum(S + p.theta * (B - S), cap) fail_exact = failure_prob_given_seller(u_s[beneficial], b, s, p).mean() if beneficial.any() else np.nan out = { "beneficial_rate": float(beneficial.mean()), "failure_rate": float(1.0 - completes.sum() / max(beneficial.sum(), 1)), "failure_rate_exact": float(fail_exact), "completion_rate": float(completes.mean()), "mean_price": float(P[completes].mean()) if completes.any() else np.nan, "seller_surplus": float((P - S)[completes].mean()) if completes.any() else np.nan, "buyer_surplus": float((B - P)[completes].mean()) if completes.any() else np.nan, "total_surplus": float((B - S)[completes].mean()) if completes.any() else np.nan, "lost_surplus": float((B - S)[beneficial & ~completes].sum() / n), } return out def draw_market_shocks(n: int, rng: np.random.Generator, p: WelfareParams = WelfareParams()) -> dict: return { "u_s": p.mu_s + p.sd_s * rng.standard_normal(n), "u_b": p.mu_b + p.sd_b * rng.standard_normal(n), "z": rng.standard_normal(n), } # Gauss–Hermite nodes for integrating the future-value shock at default (numpy only) _GH_X, _GH_W = np.polynomial.hermite_e.hermegauss(40) _GH_W = _GH_W / _GH_W.sum() def expected_loss_given_z(z: np.ndarray, b: float, s: float, p: WelfareParams) -> np.ndarray: """g(z) = pd * E_{V1}[max(L(z) - (1 - c) V1, 0)], integrating V1 by Gauss–Hermite quadrature.""" z = np.asarray(z, dtype=float) L = p.ltv * (1.0 - p.haircut(s)) * (1.0 + b + s * z) v1 = 1.0 + p.mu_g + p.sd_g * _GH_X # quadrature nodes for V1 loss = np.maximum(L[:, None] - (1.0 - p.c_liq) * v1[None, :], 0.0) return p.pd * loss @ _GH_W def expected_loss(b: float, s: float, p: WelfareParams) -> float: """Unconditional expected loss E[g(Z)] by Gauss–Hermite quadrature over Z (Proposition W2).""" return float(expected_loss_given_z(_GH_X, b, s, p) @ _GH_W) def systemic_loss(rho: float, b: float, s: float, p: WelfareParams, n_lenders: int, n_draws: int, rng: np.random.Generator, shocks: dict | None = None) -> dict[str, float]: """System loss sum_n g(Z_n) with exchangeable AVM errors (Proposition W3). Common random numbers: pass `shocks` = {'C': (n_draws,), 'E': (n_draws, n_lenders)} to compare rho values smoothly.""" if shocks is None: shocks = {"C": rng.standard_normal(n_draws), "E": rng.standard_normal((n_draws, n_lenders))} Z = np.sqrt(rho) * shocks["C"][:, None] + np.sqrt(1.0 - rho) * shocks["E"] g = expected_loss_given_z(Z.ravel(), b, s, p).reshape(Z.shape) total = g.sum(axis=1) q95 = np.quantile(total, 0.95) return { "mean": float(total.mean()), "sd": float(total.std(ddof=1)), "es95": float(total[total >= q95].mean()), "cov_pair": float(np.cov(g[:, 0], g[:, 1])[0, 1]) if n_lenders > 1 else np.nan, } def with_params(p: WelfareParams, **kw) -> WelfareParams: return replace(p, **kw)