1"""Solvers for the WP19 model: seller values, stationary equilibrium with free entry, posted prices,2financing contingency extension, planner benchmark and transition dynamics.34Equations (model_specification.md):5 (S) (r+delta) V(c) = -c + delta*w0 + Lam_s * beta * H(V(c)+U)6 (B) r U = -kappa_b + Lam_b (1-beta) sum_c H(V(c)+U) N(c)7 free entry U = k; stocks n(c) = s f(c) / (eta(c)+delta), eta(c) = Lam_s (1 - G(V(c)+U)).8"""9from __future__ import annotations1011from dataclasses import dataclass, field12from typing import Optional1314import numpy as np15from scipy import optimize1617from .model import Params181920# ----------------------------------------------------------------------------------------------21# Seller block (partial equilibrium: Lam_s and U given)22# ----------------------------------------------------------------------------------------------23def seller_value(c: float, Lam_s: float, U: float, p: Params, xtol: float = 1e-12) -> float:24 """Unique root of Phi(V) = -c + delta*w0 + Lam_s*beta*H(V+U) - (r+delta) V (Lemma 2)."""25 G, r, d = p.G, p.r, p.delta2627 def Phi(V):28 return -c + d * p.w0 + Lam_s * p.beta * float(G.H(V + U)) - (r + d) * V2930 lo = (-c + d * p.w0) / (r + d) - 1e-9 # Phi(lo) >= 0 since H >= 031 hi_support = G.hi if np.isfinite(G.hi) else G.mean() + 50.032 hi = hi_support - U + 1e-9 # H = 0 beyond the support => Phi(hi) < 033 hi = max(hi, lo + 1.0)34 guard = 035 while Phi(hi) >= 0 and guard < 60: # unbounded (heavy-tailed) support safeguard36 hi *= 2.0; guard += 137 return optimize.brentq(Phi, lo, hi, xtol=xtol, rtol=4 * np.finfo(float).eps, maxiter=500)383940@dataclass41class SellerBlock:42 V: np.ndarray43 vstar: np.ndarray44 acc: np.ndarray # conditional sale probability 1-G(v*)45 eta: np.ndarray # sale hazard46 pbar: np.ndarray # expected transaction price by type47 M: np.ndarray # mean residual valuation at v*48 H: np.ndarray49 resid: np.ndarray # residual of (S)505152def solve_sellers(Lam_s: float, U: float, p: Params) -> SellerBlock:53 V = np.array([seller_value(float(ci), Lam_s, U, p, p.root_xtol) for ci in p.c])54 vstar = V + U55 acc = np.asarray(p.G.surv(vstar), dtype=float)56 H = np.asarray(p.G.H(vstar), dtype=float)57 M = np.where(acc > 0, H / np.maximum(acc, 1e-300), 0.0)58 eta = Lam_s * acc59 pbar = V + p.beta * M60 resid = -p.c + p.delta * p.w0 + Lam_s * p.beta * H - (p.r + p.delta) * V61 return SellerBlock(V, vstar, acc, eta, pbar, M, H, resid)626364def stationary_stocks(eta: np.ndarray, p: Params):65 n = p.s * p.f / (eta + p.delta)66 S = n.sum()67 N = n / S68 Q = float((eta * n).sum())69 omega = eta * n / Q if Q > 0 else N70 return n, S, N, Q, omega717273# ----------------------------------------------------------------------------------------------74# Stationary equilibrium with free entry (Nash bargaining, baseline screen phi)75# ----------------------------------------------------------------------------------------------76@dataclass77class Equilibrium:78 theta: float79 Lam_s: float80 Lam_b: float81 U: float82 sellers: SellerBlock83 n: np.ndarray84 S: float85 B: float86 Q: float87 omega: np.ndarray88 N: np.ndarray89 free_entry_resid: float90 converged: bool91 theta_roots: list = field(default_factory=list)92 protocol: str = "nash"93 extra: dict = field(default_factory=dict)9495 # --- moments ----------------------------------------------------------------------------96 def moments(self, p: Params) -> dict:97 sb = self.sellers98 G = p.G99 pT = float((self.omega * sb.pbar).sum())100 # within-type variance of prices: beta^2 Var(w | w>=v*)101 cvar = np.array([G.cond_var(v) for v in sb.vstar])102 within = float((self.omega * (p.beta ** 2) * cvar).sum())103 between = float((self.omega * (sb.pbar - pT) ** 2).sum())104 mean_dur_listing = float((p.f * 1.0 / (sb.eta + p.delta)).sum()) # expected duration of a new listing105 sold_share = float((p.f * sb.eta / (sb.eta + p.delta)).sum())106 # mean duration of completed sales: E[tau | sale, c] = 1/(eta+delta) (exponential competing risks)107 mean_tom_sold = float((self.omega / (sb.eta + p.delta)).sum())108 pstar = frictionless_price(p)109 D = 1.0 - sb.pbar / sb.pbar[0]110 Dbar = float((self.omega * D).sum())111 welfare = welfare_flow(self, p)112 return {113 "theta": self.theta, "Lam_s": self.Lam_s, "Lam_b": self.Lam_b,114 "S": self.S, "B": self.B, "Q": self.Q, "sold_share": sold_share,115 "mean_listing_duration": mean_dur_listing, "mean_tom_sold": mean_tom_sold,116 "mean_tom_sold_months": 12.0 * mean_tom_sold,117 "buyer_mean_duration": 1.0 / (self.Lam_b * float((self.N * sb.acc).sum())),118 "price_mean": pT, "price_var_within": within, "price_var_between": between,119 "price_sd_total": float(np.sqrt(within + between)),120 "price_between_share": between / (within + between) if within + between > 0 else np.nan,121 "V_min": float(sb.V.min()), "V_max": float(sb.V.max()),122 "pbar_patient": float(sb.pbar[0]), "pbar_urgent": float(sb.pbar[-1]),123 "D_urgent": float(D[-1]), "D_mean": Dbar,124 "eta_patient": float(sb.eta[0]), "eta_urgent": float(sb.eta[-1]),125 "acc_patient": float(sb.acc[0]), "acc_urgent": float(sb.acc[-1]),126 "p_star": pstar, "Gamma": (pstar - pT) / pstar,127 "participation_ok": bool((sb.V >= p.w0 - 1e-9).all()),128 "welfare_total": welfare["total"], "welfare_sellers": welfare["sellers"],129 "welfare_buyers_net": welfare["buyers_net"], "welfare_holding_costs": welfare["holding_costs"],130 "free_entry_resid": self.free_entry_resid, "seller_resid_max": float(np.abs(sb.resid).max()),131 "converged": self.converged, "n_theta_roots": len(self.theta_roots),132 }133134135def frictionless_price(p: Params) -> float:136 """Walrasian free-entry benchmark: H(p*) = r k."""137 G = p.G138 target = p.r * p.k139 hi = G.hi if np.isfinite(G.hi) else G.mean() + 50.0140 return optimize.brentq(lambda v: float(G.H(v)) - target, G.lo - 10.0, hi, xtol=1e-13)141142143def free_entry_residual(theta: float, p: Params, U: Optional[float] = None) -> float:144 U = p.k if U is None else U145 Lam_s = p.phi * float(p.lam_s(theta))146 Lam_b = p.phi * float(p.lam_b(theta))147 sb = solve_sellers(Lam_s, U, p)148 n, S, N, Q, omega = stationary_stocks(sb.eta, p)149 return -p.kappa_b + Lam_b * (1.0 - p.beta) * float((N * sb.H).sum()) - p.r * U150151152def _build_equilibrium(theta: float, p: Params, roots, converged: bool) -> Equilibrium:153 Lam_s = p.phi * float(p.lam_s(theta))154 Lam_b = p.phi * float(p.lam_b(theta))155 sb = solve_sellers(Lam_s, p.k, p)156 n, S, N, Q, omega = stationary_stocks(sb.eta, p)157 resid = -p.kappa_b + Lam_b * (1.0 - p.beta) * float((N * sb.H).sum()) - p.r * p.k158 return Equilibrium(theta, Lam_s, Lam_b, p.k, sb, n, S, theta * S, Q, omega, N, resid, converged, list(roots))159160161def solve_equilibrium(p: Params, check_grid: bool = True) -> Equilibrium:162 """Find theta such that the free-entry residual is zero. Scans a log grid for sign changes (reports multiplicity)."""163 lo, hi = p.theta_bounds164 roots = []165 if check_grid:166 grid = np.geomspace(lo, hi, p.grid_check_points)167 vals = np.array([free_entry_residual(t, p) for t in grid])168 for i in range(len(grid) - 1):169 if vals[i] == 0.0:170 roots.append(float(grid[i]))171 elif vals[i] * vals[i + 1] < 0:172 roots.append(optimize.brentq(free_entry_residual, grid[i], grid[i + 1], args=(p,), xtol=1e-12))173 if not roots:174 # no sign change: report the end with smallest |residual|175 j = int(np.argmin(np.abs(vals)))176 return _build_equilibrium(float(grid[j]), p, [], converged=False)177 theta = roots[0] if len(roots) == 1 else roots[int(np.argmin(np.abs(np.log(np.array(roots)))))]178 return _build_equilibrium(theta, p, roots, converged=True)179 f_lo, f_hi = free_entry_residual(lo, p), free_entry_residual(hi, p)180 if f_lo * f_hi > 0:181 return _build_equilibrium(lo if abs(f_lo) < abs(f_hi) else hi, p, [], converged=False)182 theta = optimize.brentq(free_entry_residual, lo, hi, args=(p,), xtol=1e-12)183 return _build_equilibrium(theta, p, [theta], converged=True)184185186def partial_equilibrium(theta: float, p: Params) -> Equilibrium:187 """Hold tightness fixed (U = k), re-solve the seller block and stocks."""188 return _build_equilibrium(theta, p, [theta], converged=True)189190191# ----------------------------------------------------------------------------------------------192# Welfare193# ----------------------------------------------------------------------------------------------194def welfare_flow(eq: Equilibrium, p: Params) -> dict:195 """Steady-state flow surplus: sum_c eta n (E[w|w>=v*] - w0) - sum_c c n - kappa_b B - k Q."""196 sb = eq.sellers197 gains = float((sb.eta * eq.n * (sb.vstar + sb.M - p.w0)).sum())198 holding = float((p.c * eq.n).sum())199 search = p.kappa_b * eq.B200 entry = p.k * eq.Q201 total = gains - holding - search - entry202 sellers = p.s * float((p.f * (sb.V - p.w0)).sum()) # ex-ante value of a listing cohort net of own-use value203 buyers_net = 0.0 # free entry204 return {"total": total, "gains": gains, "holding_costs": holding, "search_costs": search,205 "entry_costs": entry, "sellers": sellers, "buyers_net": buyers_net}206207208def planner_optimum(p: Params, x0: Optional[np.ndarray] = None) -> dict:209 """Constrained-efficient steady state: choose theta and acceptance thresholds v*(c) to maximise flow surplus210 subject to the stock equations. Deterministic Nelder-Mead from the decentralised allocation."""211 G = p.G212213 def objective(x):214 theta = np.exp(x[0]); vstar = x[1:]215 Lam_s = p.phi * float(p.lam_s(theta))216 acc = np.asarray(G.surv(vstar)); eta = Lam_s * acc217 n = p.s * p.f / (eta + p.delta); S = n.sum(); B = theta * S; Q = (eta * n).sum()218 gains = (eta * n * (np.asarray(G.cond_mean(vstar)) - p.w0)).sum()219 return -(gains - (p.c * n).sum() - p.kappa_b * B - p.k * Q)220221 if x0 is None:222 eq = solve_equilibrium(p)223 x0 = np.concatenate([[np.log(eq.theta)], eq.sellers.vstar])224 res = optimize.minimize(objective, x0, method="Nelder-Mead",225 options={"xatol": 1e-9, "fatol": 1e-12, "maxiter": 40000, "maxfev": 80000})226 res2 = optimize.minimize(objective, res.x, method="Powell", options={"xtol": 1e-10, "ftol": 1e-13, "maxiter": 20000})227 best = res2 if res2.fun <= res.fun else res228 return {"theta": float(np.exp(best.x[0])), "vstar": best.x[1:], "welfare": -float(best.fun),229 "success": bool(best.success), "nfev": int(best.nfev) + int(res.nfev)}230231232# ----------------------------------------------------------------------------------------------233# Posted-price protocol (Proposition 8)234# ----------------------------------------------------------------------------------------------235def posted_threshold(a: float, p: Params) -> float:236 """y solving y - a = (1-G(y))/g(y) on (a, hi) (unique under log-concave g). Falls back to direct maximisation."""237 G = p.G238 hi = G.hi if np.isfinite(G.hi) else G.mean() + 50.0239240 def foc(y):241 s = float(G.surv(y)); g = float(G.pdf(y))242 return (y - a) * g - s243244 lo_y = max(a, G.lo) + 1e-10245 try:246 if foc(lo_y) < 0 and foc(hi - 1e-10) > 0:247 return optimize.brentq(foc, lo_y, hi - 1e-10, xtol=1e-12)248 except ValueError:249 pass250 res = optimize.minimize_scalar(lambda y: -float(G.surv(y)) * (y - a), bounds=(lo_y, hi), method="bounded",251 options={"xatol": 1e-12})252 return float(res.x)253254255def seller_value_posted(c: float, Lam_s: float, U: float, p: Params) -> tuple[float, float]:256 G, r, d = p.G, p.r, p.delta257258 def J(V):259 y = posted_threshold(U + V, p)260 return float(G.surv(y)) * (y - U - V), y261262 def Phi(V):263 return -c + d * p.w0 + Lam_s * J(V)[0] - (r + d) * V264265 lo = (-c + d * p.w0) / (r + d) - 1e-9266 hi_support = G.hi if np.isfinite(G.hi) else G.mean() + 50.0267 hi = hi_support - U - 1e-9268 V = optimize.brentq(Phi, lo, hi, xtol=p.root_xtol, maxiter=500)269 return V, J(V)[1]270271272def solve_sellers_posted(Lam_s: float, U: float, p: Params) -> SellerBlock:273 out = [seller_value_posted(float(ci), Lam_s, U, p) for ci in p.c]274 V = np.array([o[0] for o in out]); y = np.array([o[1] for o in out])275 acc = np.asarray(p.G.surv(y), dtype=float); H = np.asarray(p.G.H(y), dtype=float)276 M = np.where(acc > 0, H / np.maximum(acc, 1e-300), 0.0)277 eta = Lam_s * acc278 pbar = y - U # deterministic posted price279 resid = -p.c + p.delta * p.w0 + Lam_s * acc * (y - U - V) - (p.r + p.delta) * V280 return SellerBlock(V, y, acc, eta, pbar, M, H, resid)281282283def free_entry_residual_posted(theta: float, p: Params) -> float:284 Lam_s = p.phi * float(p.lam_s(theta)); Lam_b = p.phi * float(p.lam_b(theta))285 sb = solve_sellers_posted(Lam_s, p.k, p)286 n, S, N, Q, omega = stationary_stocks(sb.eta, p)287 # buyer surplus per meeting with a type-c seller: E[(w - y(c))^+] = H(y)288 return -p.kappa_b + Lam_b * float((N * sb.H).sum()) - p.r * p.k289290291def solve_equilibrium_posted(p: Params) -> Equilibrium:292 lo, hi = p.theta_bounds293 grid = np.geomspace(lo, hi, p.grid_check_points)294 vals = np.array([free_entry_residual_posted(t, p) for t in grid])295 roots = [optimize.brentq(free_entry_residual_posted, grid[i], grid[i + 1], args=(p,), xtol=1e-12)296 for i in range(len(grid) - 1) if vals[i] * vals[i + 1] < 0]297 converged = bool(roots)298 theta = roots[0] if roots else float(grid[int(np.argmin(np.abs(vals)))])299 Lam_s = p.phi * float(p.lam_s(theta)); Lam_b = p.phi * float(p.lam_b(theta))300 sb = solve_sellers_posted(Lam_s, p.k, p)301 n, S, N, Q, omega = stationary_stocks(sb.eta, p)302 resid = -p.kappa_b + Lam_b * float((N * sb.H).sum()) - p.r * p.k303 eq = Equilibrium(theta, Lam_s, Lam_b, p.k, sb, n, S, theta * S, Q, omega, N, resid, converged, roots, "posted")304 return eq305306307# ----------------------------------------------------------------------------------------------308# Financing-contingency extension (cash vs mortgage buyers, Proposition 7)309# ----------------------------------------------------------------------------------------------310@dataclass311class FinancingEquilibrium:312 theta: float313 lam_s: float314 lam_b: float315 U_c: float316 U_m: float317 pi_tilde: float318 V: np.ndarray319 vstar_c: np.ndarray320 vstar_m: np.ndarray321 eta_c: np.ndarray322 eta_m: np.ndarray323 n: np.ndarray324 m_pending: np.ndarray325 S: float326 B_c: float327 B_m: float328 Q: float329 Q_c: np.ndarray330 Q_m: np.ndarray331 pbar_c: np.ndarray332 pbar_m: np.ndarray333 Delta_b: np.ndarray334 sel_term: np.ndarray335 outside_term: float336 resid: dict337 converged: bool338 iterations: int339340341def _fin_state(theta: float, U_c: float, U_m: float, pit: float, p: Params) -> dict:342 """Seller block, stocks and buyer-side objects for given (theta, U_c, U_m, pi_tilde)."""343 G = p.G344 lam_s = float(p.lam_s(theta)); lam_b = float(p.lam_b(theta))345 nu, phi, r, d = p.nu, p.phi_m, p.r, p.delta346 kappa = nu * phi / (r + nu)347348 def vstar_m_of(Vc, c):349 return Vc + U_m + (c + p.kappa_b + r * (Vc + U_m)) / (nu * phi)350351 V = np.empty_like(p.c)352 for i, c in enumerate(p.c):353 def Phi(Vc):354 return (-c + d * p.w0 + lam_s * p.beta * (pit * float(G.H(Vc + U_c)) + (1 - pit) * kappa * float(G.H(vstar_m_of(Vc, c))))355 - (r + d) * Vc)356 lo = (-c + d * p.w0) / (r + d) - 1e-9357 hi = max((G.hi if np.isfinite(G.hi) else G.mean() + 50.0) - min(U_c, U_m) + 1e-9, lo + 1.0)358 guard = 0359 while Phi(hi) >= 0 and guard < 50:360 hi += 1.0; guard += 1361 V[i] = optimize.brentq(Phi, lo, hi, xtol=p.root_xtol, maxiter=500)362 vs_c = V + U_c363 vs_m = np.array([vstar_m_of(V[i], c) for i, c in enumerate(p.c)])364 acc_c = np.asarray(G.surv(vs_c)); acc_m = np.asarray(G.surv(vs_m))365 eta_c = lam_s * pit * acc_c; eta_m = lam_s * (1 - pit) * acc_m366 n = p.s * p.f / (eta_c + phi * eta_m + d); S = n.sum(); N = n / S367 R_c = r * U_c - (-p.kappa_b + lam_b * (1 - p.beta) * float((N * np.asarray(G.H(vs_c))).sum()))368 R_m = r * U_m - (-p.kappa_b + lam_b * (1 - p.beta) * kappa * float((N * np.asarray(G.H(vs_m))).sum()))369 Q = float(((eta_c + phi * eta_m) * n).sum()); b = Q370 xi_c = lam_b * float((N * acc_c).sum()); xi_m = lam_b * float((N * acc_m).sum())371 B_c = p.pi * b / max(xi_c, 1e-300); B_m = (1 - p.pi) * b / max(phi * xi_m, 1e-300)372 pit_cons = B_c / (B_c + B_m) if (B_c + B_m) > 0 else p.pi373 R_theta = theta - (B_c + B_m) / S374 return dict(lam_s=lam_s, lam_b=lam_b, V=V, vs_c=vs_c, vs_m=vs_m, eta_c=eta_c, eta_m=eta_m, n=n, S=S, N=N, Q=Q, B_c=B_c, B_m=B_m,375 acc_c=acc_c, acc_m=acc_m, xi_c=xi_c, xi_m=xi_m, R_c=R_c, R_m=R_m, R_pit=pit_cons - pit, R_theta=R_theta, pit_cons=pit_cons)376377378def solve_equilibrium_financing(p: Params, x0: Optional[np.ndarray] = None) -> FinancingEquilibrium:379 """Stationary equilibrium of the financing extension as a 3-equation root problem in (log theta, U_c, logit pi_tilde),380 with U_m = (k - pi U_c)/(1 - pi) imposed by ex-ante free entry. Equations: buyer HJBs for cash and mortgage buyers and381 stock-consistency of the cash share of meetings. (The theta = B/S identity is implied by the flow accounting and is382 reported as a diagnostic.)"""383 pi = p.pi384385 def unpack(x):386 theta = float(np.exp(x[0])); U_c = float(x[1]); pit = float(1 / (1 + np.exp(-x[2])))387 U_m = (p.k - pi * U_c) / (1 - pi)388 return theta, U_c, U_m, pit389390 def F(x):391 theta, U_c, U_m, pit = unpack(x)392 st = _fin_state(theta, U_c, U_m, pit, p)393 return np.array([st["R_c"] / p.r, st["R_m"] / p.r, st["R_pit"]])394395 if x0 is None:396 eq0 = solve_equilibrium(p.copy(phi=pi + (1 - pi) * p.phi_m))397 x0 = np.array([np.log(eq0.theta), p.k, np.log(pi / (1 - pi))])398 best = None399 for method in ("hybr", "lm"):400 try:401 sol = optimize.root(F, x0, method=method, tol=1e-12, options={"maxfev": 4000} if method == "hybr" else {"maxiter": 4000})402 except (ValueError, FloatingPointError, ZeroDivisionError, OverflowError):403 continue404 if best is None or np.abs(sol.fun).max() < np.abs(best.fun).max():405 best = sol406 if np.abs(sol.fun).max() < 1e-9:407 break408 sol = best409 theta, U_c, U_m, pit = unpack(sol.x)410 st = _fin_state(theta, U_c, U_m, pit, p)411 G = p.G; nu, phi, r = p.nu, p.phi_m, p.r412 V = st["V"]413 Delta_b = ((1 - p.beta) * (p.c + r * V) - p.beta * (p.kappa_b + r * U_m)) / (nu * phi)414 outside = p.beta * (U_c - U_m)415 Mc = np.asarray(G.M(st["vs_c"])); Mm = np.asarray(G.M(st["vs_m"]))416 pbar_c = V + p.beta * Mc417 pbar_m = V + p.beta * (st["vs_m"] + Mm - V - U_m) + Delta_b418 sel = p.beta * ((st["vs_m"] + Mm) - (st["vs_c"] + Mc))419 m_pending = st["eta_m"] * st["n"] / nu420 Q_c = st["eta_c"] * st["n"]; Q_m = phi * st["eta_m"] * st["n"]421 resid = {"R_c": st["R_c"], "R_m": st["R_m"], "R_pit": st["R_pit"], "free_entry": pi * U_c + (1 - pi) * U_m - p.k,422 "theta_identity": st["R_theta"], "max_abs": float(np.abs(sol.fun).max())}423 converged = bool(np.abs(sol.fun).max() < 1e-8)424 return FinancingEquilibrium(theta, st["lam_s"], st["lam_b"], U_c, U_m, pit, V, st["vs_c"], st["vs_m"], st["eta_c"], st["eta_m"],425 st["n"], m_pending, st["S"], st["B_c"], st["B_m"], st["Q"], Q_c, Q_m, pbar_c, pbar_m, Delta_b, sel,426 outside, resid, converged, int(sol.nfev))427428429# ----------------------------------------------------------------------------------------------430# Transition dynamics (Experiment 9): unanticipated shock to A or phi at t = 0431# ----------------------------------------------------------------------------------------------432@dataclass433class Transition:434 t: np.ndarray435 theta: np.ndarray436 V: np.ndarray # (T, n_c)437 eta: np.ndarray438 n: np.ndarray439 pbar: np.ndarray440 Q: np.ndarray441 S: np.ndarray442 P_obs: np.ndarray # transaction-weighted mean price443 P_fix: np.ndarray # fixed (initial) composition index444 omega: np.ndarray445 share_urgent_sales: np.ndarray446 mean_tom_sold: np.ndarray447 converged: bool448 iterations: int449 max_update: float450451452def transition(p_old: Params, p_path, T: Optional[float] = None, dt: Optional[float] = None) -> Transition:453 """p_path(t) -> Params in force at time t (piecewise constant shocks). Terminal condition = steady state of p_path(T)."""454 T = p_old.T if T is None else T455 dt = p_old.dt if dt is None else dt456 nt = int(round(T / dt)) + 1457 t = np.linspace(0.0, T, nt)458 eq0 = solve_equilibrium(p_old)459 pT = p_path(T)460 eqT = solve_equilibrium(pT)461 n_c = len(p_old.c)462 theta = np.full(nt, eqT.theta)463 # parameter paths464 A_t = np.array([p_path(ti).A for ti in t]); phi_t = np.array([p_path(ti).phi for ti in t])465 pp = pT # structural constants other than A, phi are assumed unchanged along the path466 G = pp.G; r, d, beta, k = pp.r, pp.delta, pp.beta, pp.k467 V = np.tile(eqT.sellers.V, (nt, 1)); n = np.tile(eq0.n, (nt, 1)); eta = np.zeros((nt, n_c))468 converged = False; it = 0; upd = np.inf469 for it in range(1, pp.dyn_max_iter + 1):470 # backward HJB (semi-implicit): V_j = [V_{j+1} + dt(-c + d w0 + Lam_s beta H(V_{j+1}+k))]/(1+(r+d)dt)471 V[-1] = eqT.sellers.V472 for j in range(nt - 2, -1, -1):473 Lam_s = phi_t[j] * A_t[j] * theta[j] ** pp.alpha474 V[j] = (V[j + 1] + dt * (-pp.c + d * pp.w0 + Lam_s * beta * np.asarray(G.H(V[j + 1] + k)))) / (1 + (r + d) * dt)475 # forward KFE with hazards from V_j476 n[0] = eq0.n477 for j in range(nt):478 Lam_s = phi_t[j] * A_t[j] * theta[j] ** pp.alpha479 eta[j] = Lam_s * np.asarray(G.surv(V[j] + k))480 if j < nt - 1:481 n[j + 1] = (n[j] + dt * pp.s * pp.f) / (1 + dt * (eta[j] + d)) # implicit Euler482 # free entry pointwise: r k + kappa_b = phi lam_b(theta) (1-beta) sum H(V_j + k) N_j483 theta_new = np.empty(nt)484 for j in range(nt):485 N = n[j] / n[j].sum(); Hj = float((N * np.asarray(G.H(V[j] + k))).sum())486 rhs = (r * k + pp.kappa_b) / (phi_t[j] * A_t[j] * (1 - beta) * max(Hj, 1e-300))487 # lam_b = theta^(alpha-1) => theta = rhs^(1/(alpha-1))488 theta_new[j] = rhs ** (1.0 / (pp.alpha - 1.0))489 upd = float(np.abs(theta_new - theta).max())490 theta = (1 - pp.damping) * theta + pp.damping * theta_new491 if upd < pp.dyn_tol:492 converged = True493 break494 # outputs495 vstar = V + k496 M = np.asarray(G.M(vstar)); pbar = V + beta * M497 Q = (eta * n).sum(axis=1); S = n.sum(axis=1)498 omega = eta * n / Q[:, None]499 P_obs = (omega * pbar).sum(axis=1)500 P_fix = (eq0.omega[None, :] * pbar).sum(axis=1) # pre-shock (predetermined) transaction composition501 urgent = p_old.c >= np.median(p_old.c)502 share_urgent = omega[:, urgent].sum(axis=1)503 mean_tom = (omega / (eta + d)).sum(axis=1)504 return Transition(t, theta, V, eta, n, pbar, Q, S, P_obs, P_fix, omega, share_urgent, mean_tom, converged, it, upd)505