Certified Joint Pricing Errors in Rough Heston

I. Introduction

The characteristic function of the Heston model is known in closed form, and for more than a decade one of its two standard expressions carried a flaw that went largely unnoticed. Written as in Heston (1993), the formula crosses a branch cut of the complex logarithm at long maturities, and an implementation on the principal branch returns a discontinuous transform (Kahl and Jäckel 2005; Albrecher, Mayer, Schoutens and Tistaert 2007). The error entered every strike of a maturity together, because every strike is an integral of the same function. A second episode shows the other side of the same fact. The expansion of Hagan, Kumar, Lesniewski and Woodward (2002) for the SABR model tracks the exact implied volatility closely, strike by strike, yet at low strikes and long expiries the density it implies turns negative: a butterfly, which is a portfolio of three options, receives a negative price (Hagan, Kumar, Lesniewski and Woodward 2014). Accuracy in each price had not delivered accuracy in a direction.

Both episodes say that the accuracy of an option pricer is a property of directions in strike space. This paper makes that statement quantitative. Its portfolio results hold for every pricer built on a common transform. Its bound on the transform itself is proved for the rough Heston model of El Euch and Rosenbaum (2019), whose transform solves a fractional Riccati equation with no closed form, so that the prices used in practice come from approximations: the rational approximation of Gatheral and Radoičić (2019), fractional Adams schemes (Diethelm, Ford and Freed 2004), Markovian lifts (Abi Jaber and El Euch 2019b; Abi Jaber 2019) and neural-network surrogates (Horvath, Muguruza and Tomas 2021). Their accuracy is usually reported by comparison with a slower solver on a test set. Agreement between two solvers bounds neither.

We ask what can be proved after the fact. A vector $c^{\rm fast}$ of rough Heston prices has been computed and stored. Let $c^*$ be the exact model prices of the same contracts, normalised by the forward $F_T$ and the discount factor $D$. For a portfolio $w$ and a monetary tolerance $\tau$, the task is to prove

The stored output is taken as given. An independent reference $\bar c$ is built with its own error bound, and the stored output is located relative to it:

Here $z_n$ is the error of the reference transform at the frequency $u_n$, one complex number shared by all strikes; $A$ is the matrix of known pricing coefficients; and $r$ collects quadrature, the infinite tail and arithmetic. Two questions follow. How large can $z_n$ be, and how should the bounds on the $z_n$ be combined for a portfolio?

A. What a Riccati residual costs

A reference trajectory $\widehat h$ fails to solve the Riccati equation by a residual $r(t)$, which can be bounded on every time interval by exact arithmetic. The first result converts that residual into an error of the log-characteristic function $L_T$. For a reference in the left half-plane and an admissible forward variance curve,

The weight $\xi^\lambda$ is a forward variance curve: it is the expected variance of the same model under the measure whose Girsanov kernel is half that of the share measure, which is the measure the Lewis contour sees (Theorem 3.1 and Proposition 3.2). Three features of (1.3) deserve emphasis. The constant is one; there is no Grönwall exponential, although the equation is nonlinear and carries its whole history. The damping rate $\lambda$ is the amount by which the leverage effect speeds up mean reversion under that measure, so the feature of the model that produces the equity skew is also the one that makes numerical errors contract. And the single hypothesis on the curve, nonnegativity of its Riemann–Liouville derivative, is the condition under which El Euch and Rosenbaum (2018) and Abi Jaber and El Euch (2019a) construct the model with nonnegative variance. The certificate asks nothing of the curve that the model does not already ask.

The bound is the rough counterpart of variation of constants for the Heston Riccati equation. A fractional equation has no semigroup to restart, and two classical tools take its place: convexity of the Caputo derivative on a modulus and positivity of the scalar fractional resolvent. On the six-month example of Section 7, replacing a uniform bound on the state by (1.3) lowers the complete error bound of a call spread from 1.398 to 0.379 index points, a reduction of 72.91 per cent, with the stored prices and the reference unchanged. The same proof covers every damped contour and Volterra Heston models with completely monotone kernels, the classical and the lifted Heston models among them (Remarks 3.5 and 3.6).

B. Transform error as a common factor

Once each $z_n$ lies in a disk of radius $\epsilon_n$, the set of price errors compatible with (1.2) is a sum of ellipses, one for each frequency, and its support in the direction $w$ can be written down exactly (Theorem 5.1):

where $k_i=\log m_i$ is log-moneyness, $h_*$ is the grid step and $\omega_n$ are the trapezoidal weights. The function $\widehat w$ is the Fourier transform of the payoff of the portfolio. A portfolio of options with a common expiry is one payoff, and (1.4) prices its error as one payoff. Bounding each leg separately replaces $|\widehat w(u)|$ by $\sum_i|w_i|\sqrt{m_i}$, and the ratio of the two, which we call the netting ratio, measures what is lost.

The netting ratio has a simple anatomy. For a vertical spread with log-strike gap $\Delta$ it is $|\sin(u\Delta/2)|$ up to terms of order $\Delta$: the spread is hedged against transform error of wavelength longer than twice its strike spacing and fully exposed at the Nyquist frequency $\pi/\Delta$. For a butterfly the ratio is squared. For any portfolio its mean square over a long band of frequencies is the Herfindahl index of the position (Propositions 5.2 and 5.3). Transform error therefore behaves like a factor model in frequency. Low-frequency error is systematic across strikes and is removed by hedging. High-frequency error is idiosyncratic across strikes and is removed by diversification. On the S&P 500 example the eleven adjacent call spreads keep on average 0.708 of their leg-by-leg node error, against a Herfindahl benchmark of $1/\sqrt2\approx0.707$; an equal-weight basket of all twelve calls keeps 0.076.

This part of the paper is model-free. It applies to any pricer that obtains a strip of options from one transform, whatever the model and whatever the method that produced the transform.

C. Evidence and economic scale

Section 7 applies the certificate to S&P 500 index options quoted on 18 June 2021 (Jeng and Kiliçman 2021), with the index at 4221.86. An SPX option has a multiplier of 100 dollars per index point and a minimum tick of 0.05 points below a premium of 3.00 and 0.10 above (Cboe Global Markets 2024); in 2025 an average of 3.9 million contracts changed hands each day (Cboe Global Markets 2026). On the sample date the six-month calls between 3700 and 4800 were quoted between 0.8 and 2.5 points wide. A tolerance of a quarter point is 25 dollars a contract and about one fifth of the quoted width of the 4400 call. Our tolerances of 1, 1/2 and 1/4 point thus ask whether numerical error is safely inside the bid–ask spread.

The answers are as follows. The stored rational-approximation output for the six-month 4400–4500 call spread is certified to 0.115 points, against a quoted width of about 2.5. At three months the netted bound is 0.233 points and the leg-by-leg bound 0.252, from the same stored output, reference centre, 1,025 node radii and remainders; only netting certifies a quarter point. Across 28 spreads, butterflies and baskets at six months, a half-point tolerance is certified for all 28 with netting and for 13 without. Five neighbouring roughness candidates, from $H=0.02$ to $H=0.05$, are ranked by their least-squares distance to the quotes with certified separation in all ten pairs, against seven of ten without netting, so the ranking belongs to the model and the quotes. And at the best candidate the model price of the 4400 call lies at least 0.49 points above the quoted ask: a discrepancy between model and market that the certificate attributes to the model.

The setting is one that supervisors already describe in portfolio terms. European prudent-valuation rules require a model-risk valuation adjustment for each valuation model and aggregate the individual adjustments with a flat factor of 50 per cent (Commission Delegated Regulation (EU) 2016/101, Article 11 and Annex). The Basel market-risk framework tests, desk by desk, whether the prices of the risk engine track those of the front office, with the stated aim to "measure the materiality of simplifications" in the former (Basel Committee on Banking Supervision 2019, MAR32.20). United States guidance asks that model risk be assessed "both individually and in aggregate", where the aggregate reflects reliance on common methodologies (Board of Governors of the Federal Reserve System, Federal Deposit Insurance Corporation and Office of the Comptroller of the Currency 2026). A shared transform is a common methodology in the most literal sense. For the numerical component of valuation uncertainty, (1.4) replaces a flat diversification factor by the exact one.

D. Relation to the literature

The rough Heston model arises as the scaling limit of Hawkes-based microstructure models (El Euch, Fukasawa and Rosenbaum 2018) and belongs to the affine Volterra class (Abi Jaber, Larsson and Pulido 2019); its forward variance form is due to El Euch and Rosenbaum (2018) and Gatheral and Keller-Ressel (2019). The empirical case for rough volatility, and the debate around it, rests on precise pricing at the short end, where numerical methods are under most strain (Gatheral, Jaisson and Rosenbaum 2018; Bayer, Friz and Gatheral 2016; Rømer 2022; Cont and Das 2024).

Gatheral and Radoičić (2019, 2024) derive rational approximations by matching the short- and long-time behaviour of the fractional Riccati solution, and Jeng and Kiliçman (2020, 2021) study global Padé approximations and their calibration to SPX. The stored output we audit is the third-order approximation of Gatheral and Radoičić; Section 6 records a structural property of it on the parameter set we use. Bayer and Breneis (2023a, 2023b) study Markovian kernel approximations and weak simulation schemes. Boyarchenko, de Innocentis and Levendorskii (2025) study reliable rough Heston pricing and calibration, and their modified Adams method is one of our controls. Ben Hammouda, Ben Romdhane, Samet and Tempone (2026) allocate time-discretisation and quadrature errors in single- and multilevel methods, and Hager and Kreher (2026) expand the model in the Hurst parameter. These works bound or estimate the error of a method. Our object is the error of a given output, with every remainder charged, in a given portfolio direction.

The analytical tools are established. Caputo convexity is due to Li and Liu (2018), pointwise a posteriori control of fractional equations by positive resolvents to Kopteva (2022), complete monotonicity of Mittag-Leffler functions to Simon (2015), and the strip analysis of the trapezoidal rule to Trefethen and Weideman (2014). What we add is their composition along the valuation map, from the Riccati residual through the forward variance curve to the payoff transform, the identification of the resulting kernel as a forward variance curve under a tilted measure, and the frequency anatomy of portfolio error.

Section 2 states the pricing problem. Section 3 proves the residual bound and Section 4 turns it into transform disks and price bounds. Section 5 treats portfolios. Section 6 describes the stored output and Section 7 the evidence. Appendices A–E contain the model verification, the complete proofs and the experimental record.

II. The pricing problem

A. Model, transform and price

Fix a maturity $T$ and let $S$ be the forward price for delivery at $T$, so that $S_0=F_T$. Under the pricing measure $\mathbb Q$,

with $g_\beta(t)=t^{\beta-1}/\Gamma(\beta)$, roughness index $\alpha=H+\tfrac12\in(\tfrac12,1)$, volatility of variance $\nu>0$ and correlation $-1\le\rho\le0$. The function $\xi(t)=\mathbb E V_t$ is the forward variance curve. In this form the model has no separate mean-reversion parameters: whatever drift the variance has is recorded in $\xi$, which is in principle observable from variance swaps (El Euch and Rosenbaum 2018; Gatheral and Keller-Ressel 2019).

Throughout, convolutions start at zero, $(f*g)(t)=\int_0^tf(t-s)g(s)\,ds$, and

are the fractional integral and the Caputo derivative. We work on the contour of Lewis (2001), on which the transform is $\phi_T(u-\tfrac i2)=\mathbb E[(S_T/F_T)^{1/2+iu}]$ for real $u$. It is given by a fractional Riccati equation (El Euch and Rosenbaum 2019):

With $m=K/F_T$, $k=\log m$ and $c=C/(DF_T)$, the normalised price of a call struck at $K$ is

Proposition 2.1 (the exact Riccati trajectory). Let $\alpha\in(\tfrac12,1)$, $-1\le\rho\le0$ and $u\in\mathbb R$. Equation (2.3) has a unique global solution. It is absolutely continuous on every finite interval, continuously differentiable at positive times, and satisfies $\operatorname{Re}h(t)<0$ for $t>0$. With $b_u=\tfrac12(u^2+\tfrac14)$ and $\lambda_*=\tfrac12|\rho|\nu$,

where $E_\alpha$ is the Mittag-Leffler function and the last bound remains valid for $\rho=0$.

Proof. Local Volterra contraction gives existence and uniqueness. On the imaginary axis the real part of $F_u$ is at most $-\tfrac18$ after completing a square, and the fractional maximum principle then prevents the trajectory from leaving the left half-plane. Convexity of a regularised modulus and comparison with the scalar linear equation give (2.6), which excludes blow-up, and the continuation keeps the preceding history. Appendix C.1 gives the details. Q.E.D.

The proposition identifies the object to be approximated without reference to any numerical scheme. Its two conclusions, the invariant half-plane and the scale $\lambda_*$, are the source of the stability used below.

B. Two faces of the exponent

The exponent $L_T$ has a second expression, and the pair is what makes residual control possible.

Proposition 2.2 (two faces). Let $\xi\in AC[0,T]$ and define

Then $q\in L^1(0,T)$, $g_\alpha*q=\xi$, and for every $v\in AC[0,T]$ with $v(0)=0$,

Proof. By Fubini's theorem the left side is $((g_{1-\alpha}*\xi)*v')(T)$. The function $g_{1-\alpha}*\xi$ is absolutely continuous and vanishes at zero, as does $v$, so integration by parts moves the derivative and gives $(q*v)(T)$. The identity $g_\alpha*g_{1-\alpha}=1$ gives $g_\alpha*q=\xi$. Appendix C.2 records the integrability bounds. Q.E.D.

The function $q$ is the Riemann–Liouville derivative of the curve, and it has a direct meaning. If the variance is written with an explicit drift, $V_t=V_0+\int_0^tg_\alpha(t-s)\theta(s)\,ds+\nu\int_0^tg_\alpha(t-s)\sqrt{V_s}\,dW_s$, then $\xi=V_0+g_\alpha*\theta$ and

Thus $q$ is the rate at which variance is supplied: the drift $\theta$, together with the term a power-law memory requires in order to hold the initial level $V_0$ in place. Applied to the exact solution, (2.8) reads $L_T=\int_0^T\xi(T-t)F_u(h(t))\,dt=\int_0^Tq(T-t)h(t)\,dt$, and both forms appear in El Euch and Rosenbaum (2018, Corollaries 3.2 and 3.3). The first face says what the exponent is, the Riccati polynomial accumulated against forward variance. The second says how the exponent responds to an error in the state: linearly, with weight $q$.

Assumption 2.3 (admissible curve). $\xi\in AC[0,T]$, $\xi(0)\ge0$, and $q\ge0$ almost everywhere on $(0,T)$.

Under this assumption $\xi=g_\alpha*q\ge0$. The assumption is the model's own. In the notation of (2.9) it says that $\theta(t)\,dt+V_0\,t^{-\alpha}dt/\Gamma(1-\alpha)$ is a nonnegative measure, which is the condition imposed by El Euch and Rosenbaum (2018, condition (4)) and by Abi Jaber and El Euch (2019a, Example 2.2) to construct the model with nonnegative variance. Every nondecreasing curve with $\xi(0)\ge0$ is admissible. A decreasing curve is admissible as long as it does not fall faster than the memory of the kernel allows: $\xi(t)=V_0(1-ct)$ is admissible on $[0,T]$ exactly when $cT\le1-\alpha$ (Corollary C.4). Appendix A verifies the hypotheses of the existence and transform theorems for the curve used in Section 7.

III. What a Riccati residual costs

A. A martingale that fails to be one

Let $\widehat h\in AC[0,T]$ with $\widehat h(0)=0$ be an approximate solution of (2.3), and let

be its residual. A formal computation shows what $r$ should cost. Write $\xi_t(s)=\mathbb E[V_s\mid\mathcal F_t]$ for the forward variance curve seen at time $t$, so that $d\xi_t(s)=\nu g_\alpha(s-t)\sqrt{V_t}\,dW_t$, put $X_t=\log(S_t/S_0)$, and consider

At $t=T$ this is $(S_T/S_0)^{1/2+iu}$ and at $t=0$ it is $e^{\widehat L_T}$, where $\widehat L_T=\int_0^T\xi(T-t)D^\alpha\widehat h(t)\,dt$. Itô's formula and the identity $\int_t^TD^\alpha\widehat h(T-s)\,g_\alpha(s-t)\,ds=\widehat h(T-t)$ give

with $N$ a local martingale, and taking expectations,

For the exact solution $r=0$, and this is the usual verification of (2.4). For an approximation the would-be martingale drifts, and it drifts in proportion to variance, because per unit of variance the Riccati polynomial is the generator of the pair $(X,\xi)$ acting on exponentials. A residual at time to maturity $T-t$ therefore costs the variance expected at calendar time $t$, weighted by $\widehat M_t$. Since $|\widehat M_t|$ carries the factor $(S_t/S_0)^{1/2}$, the weighting is that of a measure tilted halfway towards the share measure.

Identity (3.2) is the heuristic behind the theorem that follows, which proves the same conclusion for the exponent by deterministic means and with a weight that depends only on the model.

B. The residual bound

For $\lambda\ge0$ let $k_\lambda(t)=t^{\alpha-1}E_{\alpha,\alpha}(-\lambda t^\alpha)$ be the resolvent kernel of the scalar fractional equation, the solution of $k_\lambda+\lambda\,g_\alpha*k_\lambda=g_\alpha$, with $k_0=g_\alpha$. For an admissible curve define

Theorem 3.1 (a residual is priced by forward variance). Let $\alpha\in(\tfrac12,1)$, $-1\le\rho\le0$, $u\in\mathbb R$, and let $\xi$ satisfy Assumption 2.3. Let $h$ be the solution of (2.3) and let $\widehat h\in AC[0,T]$ satisfy $\widehat h(0)=0$ and $\operatorname{Re}\widehat h\le\varepsilon$ on $[0,T]$, with $\lambda=\tfrac12|\rho|\nu-\tfrac12\nu^2\varepsilon\ge0$. Suppose the residual (3.1) satisfies $|r|\le R$ almost everywhere for some $R\in L^\infty(0,T)$. Then, for $0\le t\le T$,

and the exponents $L_T=\int_0^T\xi(T-t)D^\alpha h(t)\,dt$ and $\widehat L_T=\int_0^T\xi(T-t)D^\alpha\widehat h(t)\,dt$ satisfy

In particular, if $R\le\delta$ then $|L_T-\widehat L_T|\le\delta\int_0^T\xi(s)\,ds$.

Proof. Put $e=h-\widehat h$. Because $F_u$ is quadratic, subtraction gives an exact linear equation,

and by Proposition 2.1 the coefficient has real part at most $-\lambda$. The Caputo derivative of a convex function of a path is at most the gradient applied to the Caputo derivative of the path. Applied to a regularised modulus this gives $D^\alpha|e|+\lambda|e|\le|r|\le R$, and positivity of the resolvent $k_\lambda$ gives (3.4). By Proposition 2.2, $L_T-\widehat L_T=\int_0^Tq(T-t)e(t)\,dt$. Since $q\ge0$, (3.4) gives $|L_T-\widehat L_T|\le(q*k_\lambda*R)(T)=(\xi^\lambda*R)(T)$. Finally $g_\alpha-k_\lambda=\lambda\,g_\alpha*k_\lambda\ge0$, so $\xi^\lambda\le g_\alpha*q=\xi$. Appendix C.2 supplies the regularity, including the passage to the terminal time. Q.E.D.