Explain calibrated Heston parameters in plain English

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)

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)))])