feat(quantlib): add Heston (1993) stochastic volatility option pricing

README.md

@@ -614,7 +614,7 @@ Intraday bars: 1m / 5m / 15m / 30m / 1H / 4H / 1D. 15 metrics + benchmark compar
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>286 tested functions across 19 modules, callable from every transport</sub></summary>
+<summary><b>Quant Library</b> <sub>289 tested functions across 19 modules, callable from every transport</sub></summary>
 
 `src/quantlib` holds one tested implementation of each piece of finance math the
 agent needs. Skills **import** these rather than carrying formulas inside

README_ar.md

@@ -610,7 +610,7 @@ LONGBRIDGE_ACCESS_TOKEN=...
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>286 دالة مختبَرة عبر 19 وحدة، قابلة للاستدعاء من كل المسارات</sub></summary>
+<summary><b>Quant Library</b> <sub>289 دالة مختبَرة عبر 19 وحدة، قابلة للاستدعاء من كل المسارات</sub></summary>
 
 يحتفظ `src/quantlib` بتنفيذ مختبَر **واحد فقط** لكل قطعة من الرياضيات المالية التي
 يحتاجها الـ agent. صارت الـ skills **تستورد** هذه الدوال بدلاً من حمل الصيغ داخل كتل

README_es.md

@@ -615,7 +615,7 @@ Barras intradía: 1m / 5m / 15m / 30m / 1H / 4H / 1D. 15 métricas + comparació
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>286 funciones probadas en 19 módulos, invocables desde cualquier transporte</sub></summary>
+<summary><b>Quant Library</b> <sub>289 funciones probadas en 19 módulos, invocables desde cualquier transporte</sub></summary>
 
 `src/quantlib` contiene una implementación probada de cada pieza de matemática
 financiera que el agente necesita. Las skills **importan** estas funciones en

README_ja.md

@@ -610,7 +610,7 @@ Intraday bars: 1m / 5m / 15m / 30m / 1H / 4H / 1D. 15 metrics + benchmark compar
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>19 モジュール・286 個のテスト済み関数、すべての経路から呼び出し可能</sub></summary>
+<summary><b>Quant Library</b> <sub>19 モジュール・289 個のテスト済み関数、すべての経路から呼び出し可能</sub></summary>
 
 `src/quantlib` は、agent が必要とする金融数学のそれぞれについて、テスト済みの実装を
 **1 つだけ**保持します。skill はこれらを **import** するようになり、markdown コード

README_ko.md

@@ -610,7 +610,7 @@ Intraday bars: 1m / 5m / 15m / 30m / 1H / 4H / 1D. 15 metrics + benchmark compar
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>19개 모듈 286개의 테스트된 함수, 모든 경로에서 호출 가능</sub></summary>
+<summary><b>Quant Library</b> <sub>19개 모듈 289개의 테스트된 함수, 모든 경로에서 호출 가능</sub></summary>
 
 `src/quantlib`는 agent가 필요로 하는 각 금융 수학에 대해 테스트된 구현을 **하나씩만**
 보유합니다. skill은 이제 이 함수들을 **import**하며, markdown 코드 블록 안에 수식을

README_zh.md

@@ -607,7 +607,7 @@ LONGBRIDGE_ACCESS_TOKEN=...
 </details>
 
 <details>
-<summary><b>Quant Library</b> <sub>19 个模块 286 个经测试的函数，四条通路皆可调用</sub></summary>
+<summary><b>Quant Library</b> <sub>19 个模块 289 个经测试的函数，四条通路皆可调用</sub></summary>
 
 `src/quantlib` 为 agent 需要的每一块金融数学各提供**一份**经测试的实现。skill 现在是
 **import** 这些函数，而不再把公式抄在 markdown 代码块里——如果你在某个 `SKILL.md`

agent/src/quantlib/volatility.py

@@ -0,0 +1,181 @@
+"""Heston (1993) stochastic volatility option pricing model.
+
+Implements semi-analytical European option pricing under the Heston model
+via numerical quadrature of the characteristic function (Lewis 2001 / Albrecher et al. 2007 formulation).
+
+Model dynamics:
+    dS_t = (r - q) * S_t * dt + sqrt(V_t) * S_t * dW_1,t
+    dV_t = kappa * (theta - V_t) * dt + sigma_v * sqrt(V_t) * dW_2,t
+    d<W_1, W_2>_t = rho * dt
+
+Parameters:
+    S0: Initial spot price (> 0)
+    K: Strike price (> 0)
+    T: Time to expiration in years (> 0)
+    r: Continuously compounded risk-free rate
+    q: Continuously compounded dividend yield
+    v0: Initial variance (> 0)
+    kappa: Mean-reversion rate (> 0)
+    theta: Long-term variance (> 0)
+    sigma_v: Volatility of variance (> 0)
+    rho: Correlation between spot and variance Brownian motions in [-1, 1]
+"""
+
+from __future__ import annotations
+
+import cmath
+import math
+
+from scipy.integrate import quad
+from src.quantlib.options import normalise_option_type
+
+__all__ = [
+    "heston_characteristic_function",
+    "heston_price",
+    "heston_feller_condition",
+]
+
+
+def heston_feller_condition(kappa: float, theta: float, sigma_v: float) -> dict[str, float | bool]:
+    """Check Feller condition (2 * kappa * theta > sigma_v**2).
+
+    When satisfied, variance process V_t strictly stays positive (never reaches zero).
+
+    Args:
+        kappa: Mean reversion rate.
+        theta: Long-term variance.
+        sigma_v: Volatility of variance.
+
+    Returns:
+        dict with keys:
+            * ``feller_ratio`` (float): 2 * kappa * theta / (sigma_v**2)
+            * ``is_satisfied`` (bool): True if feller_ratio > 1.0.
+    """
+    if sigma_v <= 0.0:
+        raise ValueError(f"sigma_v must be strictly positive, got {sigma_v}")
+    if kappa <= 0.0 or theta <= 0.0:
+        raise ValueError(f"kappa and theta must be strictly positive, got kappa={kappa}, theta={theta}")
+
+    ratio = float((2.0 * kappa * theta) / (sigma_v**2))
+    return {
+        "feller_ratio": ratio,
+        "is_satisfied": bool(ratio > 1.0),
+    }
+
+
+def heston_characteristic_function(
+    u: complex | float,
+    S0: float,
+    T: float,
+    r: float,
+    q: float,
+    v0: float,
+    kappa: float,
+    theta: float,
+    sigma_v: float,
+    rho: float,
+) -> complex:
+    """Evaluate the Heston characteristic function phi(u) for log(S_T).
+
+    Uses the Albrecher et al. (2007) formulation (Little Heston Trap stable branch).
+    """
+    i = 1j
+    x0 = math.log(S0)
+    sigma_sq = sigma_v**2
+
+    # Continuous branch under risk-neutral measure
+    term = kappa - i * rho * sigma_v * u
+    d = cmath.sqrt(term**2 + sigma_sq * (i * u + u**2))
+    g = (term - d) / (term + d)
+
+    exp_neg_dt = cmath.exp(-d * T)
+    one_minus_g_exp = 1.0 - g * exp_neg_dt
+    one_minus_g = 1.0 - g
+
+    C = (r - q) * i * u * T + (kappa * theta / sigma_sq) * (
+        (term - d) * T - 2.0 * cmath.log(one_minus_g_exp / one_minus_g)
+    )
+    D = ((term - d) / sigma_sq) * ((1.0 - exp_neg_dt) / one_minus_g_exp)
+
+    return cmath.exp(C + D * v0 + i * u * x0)
+
+
+def heston_price(
+    S0: float,
+    K: float,
+    T: float,
+    r: float,
+    v0: float,
+    kappa: float,
+    theta: float,
+    sigma_v: float,
+    rho: float,
+    option_type: str = "call",
+    q: float = 0.0,
+    integration_limit: float = 200.0,
+) -> float:
+    """Price a European option under the Heston stochastic volatility model.
+
+    Uses Lewis (2001) / Carr-Madan single-integral formulation with stable characteristic function.
+
+    Args:
+        S0: Current spot price (> 0).
+        K: Strike price (> 0).
+        T: Time to expiry in years (> 0).
+        r: Risk-free rate.
+        v0: Initial variance (> 0).
+        kappa: Mean reversion speed (> 0).
+        theta: Long-term variance (> 0).
+        sigma_v: Volatility of variance (> 0).
+        rho: Spot-variance correlation in [-1.0, 1.0].
+        option_type: 'call' or 'put'.
+        q: Continuous dividend yield.
+        integration_limit: Upper limit for numerical quadrature.
+
+    Returns:
+        Option price as a non-negative float.
+    """
+    if S0 <= 0.0 or K <= 0.0:
+        raise ValueError(f"Spot S0 and strike K must be positive, got S0={S0}, K={K}")
+    if T <= 0.0:
+        opt_type = normalise_option_type(option_type)
+        return float(max(0.0, S0 - K) if opt_type == "call" else max(0.0, K - S0))
+    if v0 < 0.0 or kappa <= 0.0 or theta <= 0.0 or sigma_v <= 0.0:
+        raise ValueError("v0 must be >= 0, and kappa, theta, sigma_v must be strictly positive")
+    if not (-1.0 <= rho <= 1.0):
+        raise ValueError(f"rho must be in [-1.0, 1.0], got {rho}")
+
+    opt_type = normalise_option_type(option_type)
+    k = math.log(S0 / K) + (r - q) * T
+
+    def integrand(u: float) -> float:
+        if u == 0.0:
+            return 0.0
+        # Centered characteristic function for u - i/2
+        phi = heston_characteristic_function(
+            u=u - 0.5j,
+            S0=1.0,
+            T=T,
+            r=0.0,
+            q=0.0,
+            v0=v0,
+            kappa=kappa,
+            theta=theta,
+            sigma_v=sigma_v,
+            rho=rho,
+        )
+        val = cmath.exp(-1j * u * k) * phi / (u**2 + 0.25)
+        return float(val.real)
+
+    integ, _ = quad(integrand, 0.0, integration_limit, limit=2000)
+    prefactor = (1.0 / math.pi) * math.sqrt(S0 * K) * math.exp(-0.5 * (r + q) * T)
+    call = float(S0 * math.exp(-q * T) - prefactor * integ)
+    call = float(max(0.0, call))
+
+    if opt_type == "call":
+        return call
+    else:
+        discounted_F = S0 * math.exp(-q * T)
+        discounted_K = K * math.exp(-r * T)
+        put = call - discounted_F + discounted_K
+        return float(max(0.0, put))

agent/src/tools/quantlib_tool.py

@@ -82,6 +82,7 @@
     "valuation.threestatement": "src.quantlib.valuation.threestatement",
     "valuation.artifact": "src.quantlib.valuation.artifact",
     "valuation.contracts": "src.quantlib.valuation.contracts",
+    "volatility": "src.quantlib.volatility",
 }
 
 #: Exported names refused because they write to a caller-supplied path. This

agent/tests/quantlib/test_volatility.py

@@ -0,0 +1,95 @@
+"""Unit and property tests for Heston (1993) stochastic volatility model."""
+
+import math
+import pytest
+from src.quantlib.volatility import (
+    heston_characteristic_function,
+    heston_price,
+    heston_feller_condition,
+)
+from src.quantlib.options import bs_price
+
+
+class TestHestonModel:
+    def test_feller_condition(self):
+        # 2 * 1.5 * 0.04 = 0.12, sigma_v**2 = 0.575**2 = 0.330625 -> ratio ~ 0.363 (< 1, violated)
+        res = heston_feller_condition(kappa=1.5, theta=0.04, sigma_v=0.575)
+        assert not res["is_satisfied"]
+        assert pytest.approx(res["feller_ratio"], abs=1e-3) == 0.363
+
+        # Satisfied case: kappa=2.0, theta=0.1, sigma_v=0.2 -> 2*2*0.1 / 0.04 = 10.0 (> 1)
+        res_sat = heston_feller_condition(kappa=2.0, theta=0.1, sigma_v=0.2)
+        assert res_sat["is_satisfied"]
+        assert pytest.approx(res_sat["feller_ratio"], rel=1e-7) == 10.0
+
+    def test_characteristic_function_at_zero(self):
+        # phi(0) = E[e^0] = 1.0
+        cf = heston_characteristic_function(
+            u=0.0,
+            S0=100.0,
+            T=1.0,
+            r=0.05,
+            q=0.02,
+            v0=0.04,
+            kappa=1.5,
+            theta=0.04,
+            sigma_v=0.3,
+            rho=-0.5,
+        )
+        assert pytest.approx(cf.real, abs=1e-7) == 1.0
+        assert pytest.approx(cf.imag, abs=1e-7) == 0.0
+
+    def test_heston_moodley_benchmark(self):
+        # Moodley (2005) Table 2: S0=100, K=100, T=0.5, r=0.0, q=0.0, v0=0.04, kappa=1.5, theta=0.04, sigma_v=0.575, rho=-0.5711
+        # Call price = 5.0272...
+        price = heston_price(
+            S0=100.0,
+            K=100.0,
+            T=0.5,
+            r=0.0,
+            v0=0.04,
+            kappa=1.5,
+            theta=0.04,
+            sigma_v=0.575,
+            rho=-0.5711,
+            option_type="call",
+        )
+        assert pytest.approx(price, abs=1e-3) == 5.027
+
+    def test_put_call_parity(self):
+        S0, K, T, r, q = 100.0, 95.0, 1.0, 0.04, 0.01
+        v0, kappa, theta, sigma_v, rho = 0.04, 2.0, 0.04, 0.3, -0.7
+        call = heston_price(S0, K, T, r, v0, kappa, theta, sigma_v, rho, option_type="call", q=q)
+        put = heston_price(S0, K, T, r, v0, kappa, theta, sigma_v, rho, option_type="put", q=q)
+        # Parity: Call - Put = S0*exp(-qT) - K*exp(-rT)
+        parity = S0 * math.exp(-q * T) - K * math.exp(-r * T)
+        assert pytest.approx(call - put, abs=1e-4) == parity
+
+    def test_convergence_to_black_scholes_zero_vol_of_vol(self):
+        # When sigma_v -> 0, Heston approaches Black-Scholes with constant variance v0=theta
+        S0, K, T, r, q = 100.0, 100.0, 1.0, 0.05, 0.0
+        vol = 0.2
+        v0 = theta = vol**2
+        h_price = heston_price(
+            S0=S0,
+            K=K,
+            T=T,
+            r=r,
+            v0=v0,
+            kappa=1.0,
+            theta=theta,
+            sigma_v=1e-4,
+            rho=0.0,
+            option_type="call",
+            q=q,
+        )
+        bs = bs_price(S=S0, K=K, T=T, r=r, sigma=vol, option_type="call", q=q)
+        assert pytest.approx(h_price, abs=1e-3) == bs
+
+    def test_invalid_parameters_raise(self):
+        with pytest.raises(ValueError):
+            heston_price(S0=-100.0, K=100.0, T=1.0, r=0.05, v0=0.04, kappa=1.0, theta=0.04, sigma_v=0.3, rho=0.0)
+        with pytest.raises(ValueError):
+            heston_price(S0=100.0, K=100.0, T=1.0, r=0.05, v0=0.04, kappa=1.0, theta=0.04, sigma_v=0.3, rho=1.5)
+        with pytest.raises(ValueError):
+            heston_feller_condition(kappa=-1.0, theta=0.04, sigma_v=0.3)