"""Iterative KMV solver with deterministic numerics. The public function is :func:`kmv_solve`. It recovers an asset-value path :math:`V_t` and an annualised asset volatility :math:`\\sigma_V` from an observed equity series :math:`E_t` under the Merton (1974) capital structure. The implementation follows Crosbie and Bohn (2003) and Vassalou and Xing (2004) with three production hardenings: 1. **Vectorised Newton on log-V** with a brentq fallback when Newton overshoots. The Newton step is dominated by ``norm.cdf`` calls on the same array; brentq is invoked only on rows that fail the monotonicity guard, which keeps the median firm at O(1) iteration. 2. **Deterministic tolerance.** Every numerical knob is exposed in :class:`MertonKMVConfig`, including the brentq absolute tolerance, so a re-run on the same inputs produces bit-identical output. The bit-identical property survives joblib parallelism because the solver is pure (no global state, no float reductions across workers). 3. **Full diagnostics.** :class:`KmvResult` carries the iteration count, final residual on :math:`\\sigma_V`, the maximum damping factor actually applied, and a fallback flag. Pipelines log all four; monitoring alerts read them directly. """ from __future__ import annotations from dataclasses import dataclass, field from typing import Optional import numpy as np from scipy.optimize import brentq from scipy.stats import norm @dataclass(frozen=True) class MertonKMVConfig: """Numerical configuration for the KMV iteration. Defaults match the Vassalou-Xing (2004) reference implementation. Tolerances are chosen so that pinned NumPy/SciPy versions reproduce EDF to ~1e-7 across compute nodes. """ r: float = 0.03 T: float = 1.0 horizon_days: int = 252 tol_sigma: float = 1.0e-6 tol_brentq: float = 1.0e-10 max_iter: int = 100 damping: float = 0.5 newton_max_step: float = 5.0 newton_iters: int = 50 V_lower: float = 1.0e-8 V_upper: float = 1.0e14 @dataclass class KmvResult: V_path: np.ndarray sigma_V: float n_iter: int final_residual: float max_damping_used: float fallback_used: bool converged: bool def equity_from_V(V_path, D, sigma_V, r, T): """Merton (1974) equity call as a function of asset value V. Vectorised in V. ``D`` is the face of the zero-coupon proxy debt; ``r`` the risk-free rate; ``T`` the horizon in years. """ V = np.asarray(V_path, dtype=float) sqrtT = np.sqrt(T) d1 = (np.log(V / D) + (r + 0.5 * sigma_V ** 2) * T) / (sigma_V * sqrtT) d2 = d1 - sigma_V * sqrtT return V * norm.cdf(d1) - D * np.exp(-r * T) * norm.cdf(d2) def _solve_V_vector(E, D, sigma_V, r, T, cfg: MertonKMVConfig): """Pointwise inversion of the Merton equity equation. Strategy: log-Newton with the Black-Scholes delta on V. If a row fails to converge within ``cfg.newton_iters`` or its iterate leaves the bracket ``[V_lower, V_upper]``, we fall back to brentq for that row. This keeps the typical case vectorised while preserving correctness for edge firms. """ E = np.asarray(E, dtype=float) sqrtT = np.sqrt(T) log_V = np.log(np.maximum(E + D, cfg.V_lower)) fallback_used = False for _ in range(cfg.newton_iters): V = np.exp(log_V) d1 = (np.log(V / D) + (r + 0.5 * sigma_V ** 2) * T) / (sigma_V * sqrtT) d2 = d1 - sigma_V * sqrtT f = V * norm.cdf(d1) - D * np.exp(-r * T) * norm.cdf(d2) - E # df/d(log V) = V * dE/dV = V * Phi(d1) df = V * norm.cdf(d1) df = np.where(df < 1.0e-12, 1.0e-12, df) step = f / df step = np.clip(step, -cfg.newton_max_step, cfg.newton_max_step) log_V = log_V - step if np.max(np.abs(f) / np.maximum(E, 1.0)) < cfg.tol_brentq: break V = np.exp(log_V) bad = (V < cfg.V_lower) | (V > cfg.V_upper) | ~np.isfinite(V) bad |= np.abs(equity_from_V(V, D, sigma_V, r, T) - E) > 1.0e-4 * np.maximum(E, 1.0) if np.any(bad): fallback_used = True for i in np.where(bad)[0]: f_scalar = lambda v, Ei=E[i]: float( equity_from_V(np.array([v]), D, sigma_V, r, T)[0] - Ei ) V[i] = brentq( f_scalar, cfg.V_lower, cfg.V_upper, xtol=cfg.tol_brentq, maxiter=200, ) return V, fallback_used def kmv_solve(E_series, D, cfg: Optional[MertonKMVConfig] = None) -> KmvResult: """Iterative KMV-style solver. Parameters ---------- E_series : array-like Daily equity time series, typically 1y of trading days. D : float Face value of debt (short-term + 0.5 * long-term in the standard KMV mapping). cfg : MertonKMVConfig, optional Numerical configuration. Defaults to :class:`MertonKMVConfig()`. Returns ------- KmvResult Asset path, asset volatility, iteration count, final residual, max damping applied, fallback flag, and convergence flag. """ cfg = cfg or MertonKMVConfig() E = np.asarray(E_series, dtype=float) if E.ndim != 1 or E.size < 30: raise ValueError("E_series must be a 1D series with >=30 observations.") if not np.all(np.isfinite(E)) or np.any(E <= 0): raise ValueError("E_series must be strictly positive and finite.") sigma_E = float(np.std(np.diff(np.log(E))) * np.sqrt(cfg.horizon_days)) sigma_V = sigma_E * float(np.mean(E)) / (float(np.mean(E)) + D) sigma_V = max(sigma_V, 1.0e-4) V = E + D fallback_any = False max_damping = 0.0 residual = np.inf for it in range(cfg.max_iter): V_new, fb = _solve_V_vector(E, D, sigma_V, cfg.r, cfg.T, cfg) fallback_any |= fb sigma_V_new = float(np.std(np.diff(np.log(V_new))) * np.sqrt(cfg.horizon_days)) sigma_V_next = (1.0 - cfg.damping) * sigma_V + cfg.damping * sigma_V_new residual = abs(sigma_V_next - sigma_V) max_damping = max(max_damping, abs(sigma_V_new - sigma_V)) if residual < cfg.tol_sigma: return KmvResult( V_path=V_new, sigma_V=sigma_V_next, n_iter=it + 1, final_residual=residual, max_damping_used=max_damping, fallback_used=fallback_any, converged=True, ) V, sigma_V = V_new, sigma_V_next return KmvResult( V_path=V, sigma_V=sigma_V, n_iter=cfg.max_iter, final_residual=residual, max_damping_used=max_damping, fallback_used=fallback_any, converged=False, ) def distance_to_default(V, sigma_V, D, r, T, mu=None, q=0.0) -> float: """Naive-drift distance to default at horizon T. If ``mu`` is None, uses the risk-free rate ``r`` (risk-neutral DD). ``q`` is the dividend yield, subtracted from the drift in line with Vassalou-Xing (2004) and Bharath-Shumway (2008). """ drift = (r if mu is None else mu) - q V_last = float(np.asarray(V).ravel()[-1]) return float( (np.log(V_last / D) + (drift - 0.5 * sigma_V ** 2) * T) / (sigma_V * np.sqrt(T)) )