diff --git a/osl/models/heston.py b/osl/models/heston.py index 7cbbaf3..61c7bfd 100644 --- a/osl/models/heston.py +++ b/osl/models/heston.py @@ -9,6 +9,8 @@ from __future__ import annotations +import math +import statistics from collections.abc import Sequence from dataclasses import dataclass @@ -18,6 +20,12 @@ _EVAL_MONTH = 5 _EVAL_DAY = 15 +# Calibration parameter bounds, in QuantLib's HestonModel parameter order +# (theta, kappa, sigma, rho, v0). Bounding theta and v0 away from zero is what +# stops the optimiser collapsing to the degenerate theta=0 corner. +_BOUNDS_LO = (1e-4, 1e-2, 1e-3, -0.999, 1e-4) +_BOUNDS_HI = (1.0, 50.0, 5.0, 0.999, 1.0) + @dataclass(frozen=True) class HestonParams: @@ -70,21 +78,46 @@ def heston_price( return float(option.NPV()) -def calibrate_heston( +def _seed_guesses( quotes: Sequence[tuple[float, float, float]], - *, spot: float, rate: float, dividend_yield: float, - initial: HestonParams | None = None, -) -> HestonCalibration: - """Calibrate Heston to ``(T_years, strike, implied_vol)`` quotes. +) -> list[HestonParams]: + """Data-driven Heston starting points for a multi-start calibration. - Uses Levenberg-Marquardt over QuantLib Heston helpers and reports the - root-mean-square implied-vol error of the fit. + Levenberg-Marquardt is local and Heston's objective is riddled with poor + local minima, so the starting point matters more than the optimiser. Seed + ``v0`` from the at-the-money variance and ``theta`` from the median variance, + then fan out over a few (kappa, sigma, rho) regimes and keep the best fit. """ - ql, _today, _dc, spot_h, r_ts, q_ts = _ql_context(spot, rate, dividend_yield) - guess = initial or HestonParams(v0=0.04, kappa=1.0, theta=0.04, sigma=0.5, rho=-0.5) + atm_var = quotes[0][2] ** 2 + best_moneyness = float("inf") + variances: list[float] = [] + for t_years, strike, vol in quotes: + forward = spot * math.exp((rate - dividend_yield) * t_years) + moneyness = abs(math.log(strike / forward)) + variances.append(vol**2) + if moneyness < best_moneyness: + best_moneyness, atm_var = moneyness, vol**2 + theta0 = statistics.median(variances) + return [ + HestonParams(v0=atm_var, kappa=2.0, theta=theta0, sigma=0.5, rho=-0.5), + HestonParams(v0=atm_var, kappa=5.0, theta=theta0, sigma=1.0, rho=-0.7), + HestonParams(v0=atm_var, kappa=1.0, theta=theta0, sigma=0.3, rho=-0.3), + ] + + +def _calibrate_once( # type: ignore[no-untyped-def] + ql, + spot_h, + r_ts, + q_ts, + quotes: Sequence[tuple[float, float, float]], + spot: float, + guess: HestonParams, +) -> tuple[HestonParams, float]: + """Run one bounded Levenberg-Marquardt fit from ``guess``; return (params, vol RMSE).""" process = ql.HestonProcess( r_ts, q_ts, spot_h, guess.v0, guess.kappa, guess.theta, guess.sigma, guess.rho ) @@ -102,24 +135,62 @@ def calibrate_heston( ql.QuoteHandle(ql.SimpleQuote(vol)), r_ts, q_ts, + # Minimise implied-vol error (not price error) so the wings — where + # prices are tiny but vols large — are not effectively ignored. + ql.BlackCalibrationHelper.ImpliedVolError, ) helper.setPricingEngine(engine) helpers.append(helper) lm = ql.LevenbergMarquardt(1e-8, 1e-8, 1e-8) - model.calibrate(helpers, lm, ql.EndCriteria(500, 50, 1e-8, 1e-8, 1e-8)) + constraint = ql.NonhomogeneousBoundaryConstraint(ql.Array(_BOUNDS_LO), ql.Array(_BOUNDS_HI)) + model.calibrate(helpers, lm, ql.EndCriteria(1000, 100, 1e-8, 1e-8, 1e-8), constraint) theta, kappa, sigma, rho, v0 = model.params() params = HestonParams(v0=v0, kappa=kappa, theta=theta, sigma=sigma, rho=rho) - # RMSE in vol via each helper's implied-vol error. sq = 0.0 for helper, (_t, _k, vol) in zip(helpers, quotes, strict=True): - model_price = helper.modelValue() try: - iv = helper.impliedVolatility(model_price, 1e-6, 200, 1e-4, 5.0) + iv = helper.impliedVolatility(helper.modelValue(), 1e-6, 500, 1e-4, 5.0) except RuntimeError: iv = vol sq += (iv - vol) ** 2 - rmse = (sq / len(quotes)) ** 0.5 if quotes else 0.0 - return HestonCalibration(params=params, rmse_vol=rmse) + rmse = (sq / len(quotes)) ** 0.5 + return params, rmse + + +def calibrate_heston( + quotes: Sequence[tuple[float, float, float]], + *, + spot: float, + rate: float, + dividend_yield: float, + initial: HestonParams | None = None, +) -> HestonCalibration: + """Calibrate Heston to ``(T_years, strike, implied_vol)`` quotes. + + Multi-start bounded Levenberg-Marquardt over QuantLib Heston helpers, + minimising implied-vol error. Reports the root-mean-square implied-vol error + of the best fit. Pass ``initial`` to pin a single starting point. + """ + if not quotes: + raise ValueError("calibrate_heston requires at least one quote") + ql, _today, _dc, spot_h, r_ts, q_ts = _ql_context(spot, rate, dividend_yield) + guesses = ( + [initial] if initial is not None else _seed_guesses(quotes, spot, rate, dividend_yield) + ) + + best_params: HestonParams | None = None + best_rmse = float("inf") + for guess in guesses: + try: + params, rmse = _calibrate_once(ql, spot_h, r_ts, q_ts, quotes, spot, guess) + except RuntimeError: + continue # this start failed to converge; try the next + if rmse < best_rmse: + best_params, best_rmse = params, rmse + + if best_params is None: + raise RuntimeError("Heston calibration failed from every starting point") + return HestonCalibration(params=best_params, rmse_vol=best_rmse) diff --git a/pages/11_Advanced_Models.py b/pages/11_Advanced_Models.py index e73f13b..4aa32c9 100644 --- a/pages/11_Advanced_Models.py +++ b/pages/11_Advanced_Models.py @@ -78,16 +78,56 @@ st.error(f"Heston calibration failed: {exc}") else: p = cal.params - st.write( - { - "v0": round(p.v0, 4), - "kappa": round(p.kappa, 4), - "theta": round(p.theta, 4), - "sigma(vol-of-vol)": round(p.sigma, 4), - "rho": round(p.rho, 4), - "RMSE (vol pts)": round(cal.rmse_vol * 100, 3), - } + feller = 2 * p.kappa * p.theta - p.sigma**2 + inst_vol = float(np.sqrt(max(p.v0, 0.0))) + lr_vol = float(np.sqrt(max(p.theta, 0.0))) + half_life_days = float(np.log(2) / p.kappa * 365) if p.kappa > 0 else float("inf") + rmse_pts = cal.rmse_vol * 100 + fit_quality = "excellent" if rmse_pts < 0.5 else "acceptable" if rmse_pts < 2 else "poor" + trend = ( + "rise toward it" + if p.theta > p.v0 + else "fall toward it" if p.theta < p.v0 else "stay near it" ) + if p.rho < -0.05: + rho_txt = "spot falls → vol rises (leverage effect), producing a **downside put skew**." + elif p.rho > 0.05: + rho_txt = "spot and vol rise together, an **upside skew** (atypical for equities)." + else: + rho_txt = "near zero, so a roughly **symmetric smile** with little skew." + + c1, c2, c3 = st.columns(3) + c1.metric("Instantaneous vol √v0", f"{inst_vol:.1%}") + c2.metric("Long-run vol √θ", f"{lr_vol:.1%}") + c3.metric("Variance shock half-life", f"{half_life_days:.0f} days") + + st.markdown( + f"- **v0 = {p.v0:.4f}** — instantaneous variance *now*; the market's current " + f"spot vol is √v0 ≈ **{inst_vol:.1%}**.\n" + f"- **κ (kappa) = {p.kappa:.4f}** — how fast variance mean-reverts; shocks decay " + f"with a half-life of ln2/κ ≈ **{half_life_days:.0f} days** (higher κ = settles " + f"faster).\n" + f"- **θ (theta) = {p.theta:.4f}** — the long-run variance it reverts to, i.e. a " + f"vol of √θ ≈ **{lr_vol:.1%}**; θ is {'above' if p.theta > p.v0 else 'below' if p.theta < p.v0 else 'at'} " + f"today's level, so the model expects vol to {trend}.\n" + f"- **σ (sigma) = {p.sigma:.4f}** — vol-of-vol; sets smile **curvature** / tail " + f"fatness — higher σ lifts both wings relative to ATM.\n" + f"- **ρ (rho) = {p.rho:.4f}** — spot/variance correlation; {rho_txt}\n" + f"- **RMSE = {rmse_pts:.3f} vol pts** — average gap between the model and market " + f"IV across the fitted quotes; this fit is **{fit_quality}**.\n" + f"- **Feller 2κθ−σ² = {feller:.4f}** — " + + ( + "≥ 0, so variance stays strictly positive." + if feller >= 0 + else "< 0, so variance can reach zero (stressed fit — see note)." + ) + ) + if feller < 0: + st.caption( + "Feller condition violated (2κθ < σ²): variance can reach zero. Common " + "when a steep short-dated equity skew forces a high vol-of-vol — read the " + "fitted parameters as a stressed best-fit, not a structural estimate." + ) near = smiles[0] model_iv = [] for k in near.strikes: @@ -107,6 +147,10 @@ lam = c1.slider("Jump intensity λ", 0.0, 3.0, 1.0) mu_j = c2.slider("Mean log-jump μ", -0.3, 0.1, -0.1) delta = c3.slider("Jump vol δ", 0.0, 0.4, 0.15) + st.caption( + "λ = expected number of jumps per year · μ = average jump size in log-return " + "(negative = downward crashes) · δ = how dispersed the jump sizes are." + ) params = MertonParams(jump_intensity=lam, jump_mean=mu_j, jump_std=delta) near = smiles[0] atm_iv = float(near.iv[int(np.argmin(np.abs(near.k)))])