From 3aec5cb5f2008464d1dc3cbc4b4ebf063614591b Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 29 Jul 2026 08:07:51 +0000 Subject: [PATCH] Add rough Heston stochastic volatility model Introduce a preliminary first-slice RoughHeston model that prices through the existing characteristic-function pricing pipeline. The variance is a rough (fractional) square-root process with Hurst exponent H, so the roughness index is alpha = H + 1/2. Unlike the classical Heston model, the characteristic function has no closed form: it is obtained by solving the fractional Riccati equation with a fractional Adams (Diethelm-Ford-Freed) predictor-corrector scheme, vectorized over the frequency grid. The model reuses the CIR variance parameters and follows the same time-changed convention as Heston, so the martingale convexity correction is applied by the pricing layer. At H = 1/2 it recovers the classical Heston model. - quantflow/sp/rough_heston.py: RoughHeston model, fractional Adams solver, and an approximate Euler sampler - quantflow_tests/test_rough_heston.py: characteristic-function properties, Heston-limit recovery, Adams convergence, and steeper-skew checks - docs: API stub, COS pricing example, nav entry, index table row, and El Euch-Rosenbaum / Gatheral-Jaisson-Rosenbaum bibliography entries Calibration, a faster Pade characteristic-function solver, a refined Monte Carlo scheme, and a full tutorial are deferred to a later phase. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01E2vqE1bifN1kmRboBroqLR --- docs/api/sp/index.md | 1 + docs/api/sp/rough_heston.md | 3 + docs/bibliography.md | 16 ++ docs/examples/rough_heston_pricer.py | 37 +++ docs/references.bib | 22 ++ mkdocs.yml | 1 + quantflow/sp/rough_heston.py | 357 +++++++++++++++++++++++++++ quantflow_tests/test_rough_heston.py | 115 +++++++++ 8 files changed, 552 insertions(+) create mode 100644 docs/api/sp/rough_heston.md create mode 100644 docs/examples/rough_heston_pricer.py create mode 100644 quantflow/sp/rough_heston.py create mode 100644 quantflow_tests/test_rough_heston.py diff --git a/docs/api/sp/index.md b/docs/api/sp/index.md index a4e5223c..966b097a 100644 --- a/docs/api/sp/index.md +++ b/docs/api/sp/index.md @@ -32,6 +32,7 @@ This page gives an overview of all stochastic processes available in the library | Process | Description | |---|---| | [Heston][quantflow.sp.heston.Heston] | Classical square-root stochastic volatility model | +| [RoughHeston][quantflow.sp.rough_heston.RoughHeston] | Rough (fractional) Heston model with Hurst exponent H < 1/2 | | [HestonJ][quantflow.sp.heston.HestonJ] | Heston model with compound Poisson jumps | | [DoubleHeston][quantflow.sp.heston.DoubleHeston] | Two independent Heston variance processes | | [DoubleHestonJ][quantflow.sp.heston.DoubleHestonJ] | Double Heston with compound Poisson jumps on the first component | diff --git a/docs/api/sp/rough_heston.md b/docs/api/sp/rough_heston.md new file mode 100644 index 00000000..f0262060 --- /dev/null +++ b/docs/api/sp/rough_heston.md @@ -0,0 +1,3 @@ +# Rough Heston process + +::: quantflow.sp.rough_heston.RoughHeston diff --git a/docs/bibliography.md b/docs/bibliography.md index 462a0f4d..470cadf9 100644 --- a/docs/bibliography.md +++ b/docs/bibliography.md @@ -86,6 +86,14 @@ Wang X., He X., Zhao Y., Zuo Z. (2017) [Parameter Estimations of Heston Model Ba
+#### eleuch_rosenbaum + +Omar El Euch, Mathieu Rosenbaum. (2019) [The characteristic function of rough Heston models](https://doi.org/10.1111/mafi.12173){target="_blank" rel="noopener"}, Mathematical Finance, 29(1):3-38 + +
+ +
+ #### gamma-ou P. Sabino, C. Petroni. (2021) [Gamma Related Ornstein-Uhlenbeck Processes and their Simulation](https://doi.org/10.1080/00949655.2020.1842408){target="_blank" rel="noopener"}, Journal of Statistical Computation and Simulation, 91(6) @@ -102,6 +110,14 @@ Jim Gatheral, Antoine Jacquier. (2014) [Arbitrage-free SVI volatility surfaces](
+#### gatheral_jaisson_rosenbaum + +Jim Gatheral, Thibault Jaisson, Mathieu Rosenbaum. (2018) [Volatility is rough](https://doi.org/10.1080/14697688.2017.1393551){target="_blank" rel="noopener"}, Quantitative Finance, 18(6):933-949 + +
+ +
+ #### gatheral_svi Jim Gatheral. (2004) [A parsimonious arbitrage-free implied volatility parameterization with application to the valuation of volatility derivatives](https://faculty.baruch.cuny.edu/jgatheral/madrid2004.pdf){target="_blank" rel="noopener"} diff --git a/docs/examples/rough_heston_pricer.py b/docs/examples/rough_heston_pricer.py new file mode 100644 index 00000000..83bd1e6a --- /dev/null +++ b/docs/examples/rough_heston_pricer.py @@ -0,0 +1,37 @@ +"""Rough Heston: price a short-maturity smile and compare with the Heston model. + +The rough model (Hurst H < 1/2) produces a steeper short-maturity skew than the +classical Heston model with the same parameters, without adding jumps. +""" + +from quantflow.options.inputs import OptionType +from quantflow.options.pricer import OptionPricer, OptionPricingMethod +from quantflow.sp.heston import Heston +from quantflow.sp.rough_heston import RoughHeston + +TTM = 0.1 + +rough = OptionPricer( + model=RoughHeston.create(vol=0.5, kappa=1.5, sigma=0.7, rho=-0.6, hurst=0.1), + method=OptionPricingMethod.COS, +) +heston = OptionPricer( + model=Heston.create(vol=0.5, kappa=1.5, sigma=0.7, rho=-0.6), + method=OptionPricingMethod.COS, +) + + +def implied_vol(pricer: OptionPricer, strike: float) -> float: + option_type = OptionType.PUT if strike < 100.0 else OptionType.CALL + price = pricer.price(option_type=option_type, strike=strike, forward=100.0, ttm=TTM) + return float(price.black.iv) + + +print(f"Short-maturity smile at ttm={TTM} (forward=100)") +print(f"{'strike':>8}{'rough IV':>12}{'heston IV':>12}") +for strike in (90.0, 95.0, 100.0, 105.0, 110.0): + print( + f"{strike:>8.1f}" + f"{implied_vol(rough, strike):>12.4f}" + f"{implied_vol(heston, strike):>12.4f}" + ) diff --git a/docs/references.bib b/docs/references.bib index 1afe76f1..3e5c376f 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -90,6 +90,17 @@ @article{ekf url={https://www.sciencedirect.com/science/article/pii/S2405896317324758}, } +@article{eleuch_rosenbaum, + title={The characteristic function of rough Heston models}, + author={Omar El Euch and Mathieu Rosenbaum}, + journal={Mathematical Finance}, + year={2019}, + volume={29}, + number={1}, + pages={3--38}, + url={https://doi.org/10.1111/mafi.12173}, +} + @article{gamma-ou, title={Gamma Related Ornstein-Uhlenbeck Processes and their Simulation}, author={Sabino, P. & Petroni, C.}, @@ -111,6 +122,17 @@ @article{gatheral_jacquier url={https://doi.org/10.1080/14697688.2013.819986}, } +@article{gatheral_jaisson_rosenbaum, + title={Volatility is rough}, + author={Jim Gatheral and Thibault Jaisson and Mathieu Rosenbaum}, + journal={Quantitative Finance}, + year={2018}, + volume={18}, + number={6}, + pages={933--949}, + url={https://doi.org/10.1080/14697688.2017.1393551}, +} + @misc{gatheral_svi, title={A parsimonious arbitrage-free implied volatility parameterization with application to the valuation of volatility derivatives}, diff --git a/mkdocs.yml b/mkdocs.yml index 01727a3e..743d1c9e 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -96,6 +96,7 @@ nav: - Jump Diffusion: api/sp/jump_diffusion.md - Ornstein-Uhlenbeck: api/sp/ou.md - Poisson Process: api/sp/poisson.md + - Rough Heston Model: api/sp/rough_heston.md - Wiener Process: api/sp/wiener.md - Copulas: api/sp/copula.md - Timeseries Analysis: diff --git a/quantflow/sp/rough_heston.py b/quantflow/sp/rough_heston.py new file mode 100644 index 00000000..d4395823 --- /dev/null +++ b/quantflow/sp/rough_heston.py @@ -0,0 +1,357 @@ +from __future__ import annotations + +from typing import Self + +import numpy as np +from pydantic import Field +from scipy.special import gamma +from typing_extensions import Annotated, Doc + +from quantflow.ta.paths import Paths +from quantflow.utils.types import FloatArrayLike, Vector + +from .base import StochasticProcess1D +from .cir import CIR + + +class RoughHeston(StochasticProcess1D): + r"""The rough Heston stochastic volatility model. + + The rough Heston model of + [El Euch and Rosenbaum](../../bibliography.md#eleuch_rosenbaum) replaces the + Markovian square-root variance of the classical + [Heston][quantflow.sp.heston.Heston] model with a rough (fractional) + square-root process. The variance has long memory and is not diffusive at + short time scales (see + [Volatility is rough](../../bibliography.md#gatheral_jaisson_rosenbaum)), + which produces the steep short-maturity skew observed in equity and crypto + option markets without resorting to jumps. + + Following the same time-changed convention as the classical + [Heston][quantflow.sp.heston.Heston] model, the driving log-price $x_t$ is a + Brownian motion time changed by the rough variance $v_t$, + + \begin{equation} + \begin{aligned} + d x_t &= \sqrt{v_t}\, d w^1_t \\ + v_t &= v_0 + \frac{1}{\Gamma(\alpha)} \int_0^t (t-s)^{\alpha-1} + \kappa(\theta - v_s)\, ds + + \frac{1}{\Gamma(\alpha)} \int_0^t (t-s)^{\alpha-1} + \nu \sqrt{v_s}\, d w^2_s + \end{aligned} + \end{equation} + + with correlation $\rho\, dt = {\tt E}[d w^1 d w^2]$ and roughness index + $\alpha = H + \tfrac{1}{2}$, where $H$ is the + [Hurst exponent](../../glossary.md#hurst-exponent). The choice $H = 1/2$ + (so $\alpha = 1$) recovers the classical Heston model. The deterministic + [convexity correction](../../theory/convexity_correction.md) that enforces + the martingale property is applied by the pricing layer, exactly as for the + other stochastic-volatility models. + + The characteristic function has no closed form. It is computed by solving a + fractional Riccati equation numerically with a fractional Adams + predictor-corrector scheme (see + [characteristic_exponent][.characteristic_exponent]). + """ + + variance_process: CIR = Field( + default_factory=CIR, + description=( + "The rough variance is a fractional extension of the " + "[Cox-Ingersoll-Ross][quantflow.sp.cir.CIR] process, from which the " + "mean-reversion speed, long-term variance, vol of vol and initial " + "variance are taken" + ), + ) + rho: float = Field( + default=0, + ge=-1, + le=1, + description=( + "Correlation between the Brownian motions, it provides " + "the leverage effect and therefore the skewness of the distribution" + ), + ) + hurst: float = Field( + default=0.1, + gt=0, + le=0.5, + description=( + "Hurst exponent H of the variance process, the roughness index is " + "alpha = H + 1/2 in (1/2, 1]. Lower H means rougher paths and a " + "steeper short-maturity skew. H = 1/2 recovers the classical Heston " + "model" + ), + ) + adams_steps: int = Field( + default=200, + ge=10, + description=( + "Number of time steps of the fixed-step fractional Adams solver used " + "for the Riccati equation on the interval [0, t]. The solver cost is " + "quadratic in this value for each frequency, so increase it only when " + "short-maturity accuracy requires it" + ), + ) + + @property + def alpha(self) -> float: + """Roughness index of the variance process, $\\alpha = H + 1/2$.""" + return self.hurst + 0.5 + + @classmethod + def create( + cls, + *, + rate: Annotated[float, Doc("Initial rate of the variance process")] = 1.0, + vol: Annotated[ + float, + Doc( + "Volatility of the price process, normalized by the " + "square root of time, as time tends to infinity " + "(the long term standard deviation)" + ), + ] = 0.5, + kappa: Annotated[ + float, + Doc("Mean reversion speed for the variance process"), + ] = 1, + sigma: Annotated[ + float, Doc("Volatility of the variance process (a.k.a. the vol of vol)") + ] = 0.8, + rho: Annotated[ + float, + Doc( + "Correlation between the Brownian motions of the " + "variance and price processes" + ), + ] = 0, + theta: Annotated[ + float | None, + Doc( + "Long-term mean of the variance process. " + r"If `None`, it defaults to the variance given by ${\tt vol}^2$." + ), + ] = None, + hurst: Annotated[ + float, + Doc("Hurst exponent H in (0, 1/2]; H = 1/2 recovers the Heston model"), + ] = 0.1, + adams_steps: Annotated[ + int, Doc("Number of steps of the fractional Adams Riccati solver") + ] = 200, + ) -> Self: + r"""Create a rough Heston model. + + To understand the parameters lets introduce the following notation: + + \begin{align} + {\tt var} &= {\tt vol}^2 \\ + v_0 &= {\tt rate}\cdot{\tt var} + \end{align} + """ + variance = vol * vol + return cls( + variance_process=CIR( + rate=rate * variance, + kappa=kappa, + sigma=sigma, + theta=theta if theta is not None else variance, + ), + rho=rho, + hurst=hurst, + adams_steps=adams_steps, + ) + + def characteristic_exponent( + self, + t: Annotated[FloatArrayLike, Doc("Time horizon (must be a scalar)")], + u: Annotated[Vector, Doc("Characteristic exponent argument")], + ) -> Vector: + r"""Characteristic exponent of the rough Heston model. + + The characteristic function is semi-closed and given by + + \begin{equation} + \Phi_{x_t, u} = \exp\left( + \kappa\theta\, I^1 \psi(t) + + v_0\, I^{1-\alpha} \psi(t) + \right) + \end{equation} + + where $I^1$ is the ordinary integral, $I^{1-\alpha}$ is the + Riemann-Liouville fractional integral of order $1-\alpha$, and + $\psi(\cdot) = \psi(u, \cdot)$ solves the fractional Riccati equation + $\psi(u, 0) = 0$, + + \begin{equation} + D^\alpha \psi(t) = -\tfrac{1}{2}u^2 + + (i u \rho \nu - \kappa)\, \psi(t) + + \tfrac{1}{2}\nu^2\, \psi(t)^2 + \end{equation} + + with $\kappa$ the mean-reversion speed, $\theta$ the long-term variance, + $\nu$ the vol of vol, $\rho$ the correlation and $v_0$ the initial + variance. The equation is integrated with a fractional Adams + predictor-corrector scheme on a grid of + [adams_steps][..adams_steps] points. + + The characteristic exponent returned is $\phi = -\log \Phi$, consistent + with the library convention $\Phi = e^{-\phi}$. + """ + if np.ndim(t) != 0: + raise ValueError( + "RoughHeston.characteristic_exponent requires a scalar time " + "horizon; the fractional Riccati equation is solved on a single " + "interval [0, t]" + ) + u_arr = np.atleast_1d(np.asarray(u, dtype=complex)) + psi = self._solve_fractional_riccati(float(t), u_arr) + vp = self.variance_process + dt = float(t) / self.adams_steps + i1 = dt * (0.5 * psi[0] + psi[1:-1].sum(axis=0) + 0.5 * psi[-1]) + i_frac = self._fractional_integral_weights(1.0 - self.alpha) @ psi + log_phi = vp.kappa * vp.theta * i1 + vp.rate * i_frac + phi = -log_phi + return phi[0] if np.ndim(u) == 0 else phi + + def _riccati_f(self, u: np.ndarray, x: np.ndarray) -> np.ndarray: + """Right-hand side of the fractional Riccati equation, vectorized over u.""" + vp = self.variance_process + kappa = vp.kappa + nu = vp.sigma + return ( + -0.5 * u * u + (1j * u * self.rho * nu - kappa) * x + 0.5 * nu * nu * x * x + ) + + def _solve_fractional_riccati(self, t: float, u: np.ndarray) -> np.ndarray: + """Fractional Adams (Diethelm-Ford-Freed) predictor-corrector solver. + + Returns the Riccati trajectory psi on the uniform grid of + ``adams_steps + 1`` points, with shape ``(adams_steps + 1, u.size)``. + Every operation is vectorized over the frequency grid ``u``. + """ + alpha = self.alpha + n = self.adams_steps + h = t / n + ha = h**alpha + inv_ga1 = ha / gamma(alpha + 1.0) + inv_ga2 = ha / gamma(alpha + 2.0) + a1 = alpha + 1.0 + + m = u.size + psi = np.zeros((n + 1, m), dtype=complex) + fh = np.zeros((n + 1, m), dtype=complex) + fh[0] = self._riccati_f(u, psi[0]) + + for k in range(n): + j = np.arange(k + 1) + # predictor (fractional Adams-Bashforth) weights + b = (k + 1 - j) ** alpha - (k - j) ** alpha + psi_pred = inv_ga1 * (b @ fh[: k + 1]) + f_pred = self._riccati_f(u, psi_pred) + # corrector (fractional Adams-Moulton) weights + p = k - j + a = (p + 2) ** a1 + p**a1 - 2 * (p + 1) ** a1 + a[0] = k**a1 - (k - alpha) * (k + 1) ** alpha + psi[k + 1] = inv_ga2 * (a @ fh[: k + 1] + f_pred) + fh[k + 1] = self._riccati_f(u, psi[k + 1]) + return psi + + def _fractional_integral_weights(self, beta: float) -> np.ndarray: + r"""Product-trapezoidal weights for the Riemann-Liouville integral. + + Returns the vector $c$ of length ``adams_steps + 1`` such that + $I^\beta \psi(t) \approx c \cdot \psi$ on the solver grid. For + $\beta = 0$ (i.e. $\alpha = 1$) the weights collapse to picking + $\psi(t)$, recovering the Heston limit. + """ + n = self.adams_steps + b1 = beta + 1.0 + j = np.arange(n + 1) + # middle-node weights (valid for 1 <= j <= n - 1); the (n - 1 - j) base + # is negative only at j = n, whose weight is overwritten below, so clip + # it to avoid a fractional power of a negative number + low = np.clip(n - 1 - j, 0, None) + c = (n + 1 - j) ** b1 + low**b1 - 2 * (n - j) ** b1 + # first node + c[0] = (n - 1) ** b1 - (n - 1 - beta) * n**beta + # last node (integrand endpoint) + c[n] = 1.0 + h = 1.0 / n # dt = t / n; the t**beta factor cancels with the psi grid step + return (h**beta / gamma(beta + 2.0)) * c + + def sample( + self, + n: Annotated[int, Doc("Number of sample paths")], + time_horizon: Annotated[float, Doc("Time horizon")] = 1, + time_steps: Annotated[int, Doc("Number of discrete time steps")] = 100, + ) -> Paths: + dw1 = Paths.normal_draws(n, time_horizon, time_steps) + dw2 = Paths.normal_draws(n, time_horizon, time_steps) + return self.sample_from_draws(dw1, dw2) + + def sample_from_draws( + self, + draws: Annotated[ + Paths, + Doc( + "Pre-drawn standard normal increments for the variance Brownian motion" + ), + ], + *args: Annotated[ + Paths, + Doc( + "Optional pre-drawn increments for the independent Brownian motion; " + "new draws are generated if omitted" + ), + ], + ) -> Paths: + r"""Sample paths with an explicit Euler scheme on the Volterra equation. + + This is a simple approximate scheme: the fractional variance is + discretized as a lower-triangular convolution of the power-law kernel + $g(\tau) = \tau^{\alpha-1}/\Gamma(\alpha)$ with the Euler increments of + the square-root process, and the driftless log-price integrates + $d x = \sqrt{v}\, d w^1$ (the martingale convexity correction is applied + by the pricing layer). It is intended for illustration only; the + characteristic-function pricing path does not use it. + """ + if args: + path2 = args[0] + else: + path2 = Paths.normal_draws(draws.samples, draws.t, draws.time_steps) + vp = self.variance_process + alpha = self.alpha + kappa = vp.kappa + nu = vp.sigma + theta = vp.theta + v0 = vp.rate + dt = draws.dt + steps = draws.time_steps + sqrt_dt = np.sqrt(dt) + + # power-law kernel weights g(k*dt) for k = 1, 2, ... + lag = np.arange(1, steps + 1) + g = (lag * dt) ** (alpha - 1.0) / gamma(alpha) + + db = draws.data * sqrt_dt # variance Brownian increments + dw_price = ( + self.rho * draws.data + np.sqrt(1.0 - self.rho * self.rho) * path2.data + ) * sqrt_dt + + v = np.zeros(draws.data.shape) + x = np.zeros(draws.data.shape) + v[0] = v0 + # forcing increment at each past step: kappa(theta - v_j) dt + nu sqrt(v_j) dB_j + forcing = np.zeros(draws.data.shape) + for i in range(steps): + vplus = np.clip(v[i], 0.0, None) + forcing[i] = kappa * (theta - vplus) * dt + nu * np.sqrt(vplus) * db[i] + # convolve the kernel with all forcing increments up to i: + # v[i+1] = v0 + sum_{j=0}^{i} g_{i-j} * forcing[j] + weights = g[i::-1] + v[i + 1] = v0 + np.tensordot(weights, forcing[: i + 1], axes=(0, 0)) + x[i + 1] = x[i] + np.sqrt(vplus) * dw_price[i] + return Paths(t=draws.t, data=x) diff --git a/quantflow_tests/test_rough_heston.py b/quantflow_tests/test_rough_heston.py new file mode 100644 index 00000000..8742e36d --- /dev/null +++ b/quantflow_tests/test_rough_heston.py @@ -0,0 +1,115 @@ +import numpy as np +import pytest + +from quantflow.options.inputs import OptionType +from quantflow.options.pricer import OptionPricer, OptionPricingMethod +from quantflow.sp.heston import Heston +from quantflow.sp.rough_heston import RoughHeston +from quantflow_tests.utils import characteristic_tests + + +@pytest.fixture +def rough_heston() -> RoughHeston: + return RoughHeston.create( + vol=0.5, kappa=1.5, sigma=0.7, rho=-0.5, hurst=0.1, adams_steps=200 + ) + + +def test_alpha(rough_heston: RoughHeston) -> None: + assert rough_heston.alpha == pytest.approx(0.6) + assert rough_heston.hurst == 0.1 + + +def test_characteristic(rough_heston: RoughHeston) -> None: + assert rough_heston.variance_process.is_positive is True + assert rough_heston.characteristic(1, 0) == 1 + m = rough_heston.marginal(0.5) + characteristic_tests(m) + assert m.mean() == pytest.approx(0.0, abs=1e-6) + + +def test_scalar_time_required(rough_heston: RoughHeston) -> None: + with pytest.raises(ValueError): + rough_heston.characteristic_exponent(np.array([0.5, 1.0]), 1.0) + + +def test_heston_limit_characteristic() -> None: + """At H=0.5 the rough variance is Markovian and the model is Heston.""" + kwargs = dict(vol=0.5, kappa=2.0, sigma=0.6, rho=-0.4) + rough = RoughHeston.create(hurst=0.5, adams_steps=300, **kwargs) + heston = Heston.create(**kwargs) + mr = rough.marginal(1.0) + mh = heston.marginal(1.0) + u = np.linspace(0.0, 8.0, 40) + # the martingale-corrected characteristic functions must agree + np.testing.assert_allclose( + mr.characteristic_corrected(u), mh.characteristic_corrected(u), atol=1e-3 + ) + assert mr.std() == pytest.approx(mh.std(), rel=1e-3) + + +def test_heston_limit_prices() -> None: + """At H=0.5 the COS option prices must match the classical Heston model.""" + kwargs = dict(vol=0.5, kappa=2.0, sigma=0.6, rho=-0.4) + rough = OptionPricer( + model=RoughHeston.create(hurst=0.5, adams_steps=300, **kwargs), + method=OptionPricingMethod.COS, + ) + heston = OptionPricer(model=Heston.create(**kwargs), method=OptionPricingMethod.COS) + for strike in (70.0, 85.0, 100.0, 115.0, 130.0): + pr = rough.price( + option_type=OptionType.CALL, strike=strike, forward=100.0, ttm=1.0 + ) + ph = heston.price( + option_type=OptionType.CALL, strike=strike, forward=100.0, ttm=1.0 + ) + assert pr.price == pytest.approx(ph.price, abs=1e-4) + + +def test_adams_convergence() -> None: + """Prices are stable once the fractional Adams grid is fine enough.""" + kwargs = dict(vol=0.5, kappa=1.5, sigma=0.7, rho=-0.5, hurst=0.2) + prices = [] + for steps in (100, 200, 400): + pricer = OptionPricer( + model=RoughHeston.create(adams_steps=steps, **kwargs), + method=OptionPricingMethod.COS, + ) + prices.append( + pricer.price( + option_type=OptionType.CALL, strike=100.0, forward=100.0, ttm=0.5 + ).price + ) + assert prices[2] == pytest.approx(prices[1], abs=1e-4) + assert prices[1] == pytest.approx(prices[0], abs=1e-3) + + +def test_rough_skew_steeper_than_heston() -> None: + """A rough model (H<1/2) produces a steeper short-maturity skew than Heston.""" + kwargs = dict(vol=0.5, kappa=1.5, sigma=0.7, rho=-0.6) + rough = OptionPricer( + model=RoughHeston.create(hurst=0.1, adams_steps=200, **kwargs), + method=OptionPricingMethod.COS, + ) + heston = OptionPricer(model=Heston.create(**kwargs), method=OptionPricingMethod.COS) + ttm = 0.1 + + def skew(pricer: OptionPricer) -> float: + lo = pricer.price( + option_type=OptionType.PUT, strike=90.0, forward=100.0, ttm=ttm + ) + hi = pricer.price( + option_type=OptionType.CALL, strike=110.0, forward=100.0, ttm=ttm + ) + return float(lo.black.iv) - float(hi.black.iv) + + assert skew(rough) > skew(heston) + + +def test_sample(rough_heston: RoughHeston) -> None: + np.random.seed(42) + paths = rough_heston.sample(2000, time_horizon=1.0, time_steps=100) + assert paths.samples == 2000 + assert paths.time_steps == 100 + assert np.all(np.isfinite(paths.data)) + assert paths.data[0].std() == 0.0