Double Heston Model Explained

• Double Heston model is a two-factor stochastic volatility framework that extends the classic Heston model to capture complex spot/volatility correlation and volatility-of-volatility clustering.
• It employs advanced methods such as Fourier inversion and Almost-Exact Simulation to enhance accuracy and efficiency in pricing exotic options.
• Calibration techniques using joint characteristic functions and maximum likelihood estimation ensure robust parameter inference and risk management.

The Double Heston model is a two-factor extension of the well-known Heston stochastic volatility model, designed to capture more complex dynamics observed in equity and foreign exchange option markets, particularly those arising from stochastic spot/volatility correlation and volatility-of-volatility clustering. The model prescribes asset price dynamics governed by two independent Cox–Ingersoll–Ross (CIR) variance factors, each with distinct means, mean-reversion speeds, volatilities-of-volatility, and asset-variance correlation structures, producing richer implied volatility surfaces and more realistic exotic option pricing (Higgins, 1 Feb 2026, Dimitrov et al., 22 Dec 2025, Alaya et al., 28 Jan 2025).

1. Model Specification and Mathematical Framework

The standard risk-neutral form of the Double Heston model for a non-dividend-paying asset $S_t$ with two stochastic variances $\nu_t^{(1)}, \nu_t^{(2)}$ is:

$$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$

with

$$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$

Key parameters include the risk-free rate $r$, mean-reversions $\kappa_i > 0$, long-run variances $\theta_i > 0$, volatilities-of-volatility $\sigma_i > 0$, and asset-variance correlations $\rho_i \in [-1,1]$. Positivity is enforced by the Feller conditions $2\kappa_i\theta_i \geq \sigma_i^2$ (Dimitrov et al., 22 Dec 2025).

A prominent variant focuses on the case where both variance factors share $\nu_t^{(1)}, \nu_t^{(2)}$0, $\nu_t^{(1)}, \nu_t^{(2)}$1, but have different $\nu_t^{(1)}, \nu_t^{(2)}$2 and crucially, different spot–vol correlations $\nu_t^{(1)}, \nu_t^{(2)}$3. The total variance is $\nu_t^{(1)}, \nu_t^{(2)}$4, and correlation parameters are often expressed as $\nu_t^{(1)}, \nu_t^{(2)}$5, thus enabling stochastic spot/volatility correlation effects (Higgins, 1 Feb 2026).

2. Calibration and Statistical Estimation

Calibration of the Double Heston model to vanilla options employs joint characteristic function methods and least-squares fits to implied volatility surfaces. The parameter vector $\nu_t^{(1)}, \nu_t^{(2)}$6 is typically bounded as $\nu_t^{(1)}, \nu_t^{(2)}$7, $\nu_t^{(1)}, \nu_t^{(2)}$8, $\nu_t^{(1)}, \nu_t^{(2)}$9, $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$0, $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$1.

Fourier inversion approaches, such as Carr–Madan single-integral or Lewis’ $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$2 representations, are employed. The affine structure yields a joint characteristic function: $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$3 with explicit formulas for $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$4 via Riccati ODEs (Higgins, 1 Feb 2026).

For statistical inference, global parameter estimation is performed via maximum likelihood (MLE) and conditional least squares (CLS) based on continuous-time observations. In the ergodic regime, the drift parameter vector $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$5 admits closed-form estimators, with strong consistency and asymptotic normality following from martingale central limit theory: $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$6 where $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$7 is a design matrix, and $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$8 is the diffusion covariance (Alaya et al., 28 Jan 2025).

3. Simulation Schemes and Numerical Methods

Exact simulation of the CIR variance factors is achieved via non-central chi-square sampling: $$\begin{aligned} dS_t &= r S_t dt + \sqrt{\nu_t^{(1)}} S_t dW_t^{(1)} + \sqrt{\nu_t^{(2)}} S_t dW_t^{(2)}, \qquad S_0 = s>0, \ d\nu_t^{(i)} &= \kappa_i (\theta_i - \nu_t^{(i)}) dt + \sigma_i \sqrt{\nu_t^{(i)}} dW_t^{(2+i)}, \qquad \nu_0^{(i)} = v_i>0, \quad i=1,2, \end{aligned}$$9 where $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$0, $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$1, and $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$2.

The Almost-Exact Simulation (AES) scheme couples this sampling with a partially corrected update for $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$3, replacing two of the four stochastic integrals in the discretization by exact CIR-integrals. The AES method substantially outperforms Euler–Maruyama in accuracy and time-step economy, especially when the number of simulation steps matches exercise opportunities in Bermudan and American path-dependent options. For typical scenarios, AES achieves pricing errors $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$4 for American puts ($$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$5, $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$6) where Euler errors exceed $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$7, with AES CPU time $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$8 s versus $$\langle W^{(1)}, W^{(2+i)}\rangle_t = \rho_i t, \quad \langle W^{(1)}, W^{(2)}\rangle_t = \langle W^{(2)}, W^{(3)}\rangle_t = \langle W^{(2)}, W^{(4)}\rangle_t = 0.$$9 s for Euler (Dimitrov et al., 22 Dec 2025).

4. Stochastic Spot/Volatility Correlation and Model Dynamics

The Double Heston construction enables stochastic spot/volatility correlation: $r$0 so that $r$1 fluctuates within $r$2 as the relative weights of the two variance factors evolve (Higgins, 1 Feb 2026).

Instantaneous correlation between spot and the spot–volatility correlation,

$r$3

translates into a positive “skew–spot beta,” reproducing phenomena such as the risk-reversal beta observed in FX markets.

Exotic derivative pricing is affected: the additivity and stochasticity of correlation generate path-dependent implied skew and realized correlations, which modify one-touch and barrier option pricing as well as volatility swap fair strikes. Numerical experiments show price uplifts for one-touch options of order $r$4–$r$5 for realistic parameters; volatility-swap strikes experience upward shifts of $r$6 bps as $r$7 increases from $r$8 to $r$9 (Higgins, 1 Feb 2026).

5. Stationarity and Ergodicity

Under “subcritical” conditions ($\kappa_i > 0$0, $\kappa_i > 0$1, $\kappa_i > 0$2, $\kappa_i > 0$3), the Double Heston process admits a unique invariant distribution. The system exhibits exponential ergodicity; averages of functionals converge rapidly to their stationary means. The stationary joint law can be expressed in terms of affine transforms and is characterized via solutions to two-dimensional Riccati ODEs.

A direct implication is the strong law of large numbers for long-time empirical averages, enabling robust risk management (e.g., long-run realized variance targeting) and statistical inference for model parameters (Alaya et al., 28 Jan 2025).

6. Applications, Performance, and Practical Guidance

The Double Heston model is calibrated to market data using least-squares fits to implied volatilities and global optimizers (e.g., Differential Evolution followed by Levenberg-Marquardt). In single-expiry calibration, $\kappa_i > 0$4 is often fixed for robustness.

AES simulation excels for Bermudan options when the number of time steps matches the number of exercise dates, rendering Monte Carlo option prices virtually bias-free. For American options, adopting a reasonable exercise grid (e.g., 12 or 50 per year) preserves high accuracy. Typical runs use $\kappa_i > 0$5–$\kappa_i > 0$6 paths.

The model is particularly effective for capturing risk-reversal beta term structure, stochastic skew effects, and providing realistic pricing for barrier products and volatility swaps in foreign exchange and equity option markets (Higgins, 1 Feb 2026, Dimitrov et al., 22 Dec 2025).

7. Limitations and Best-Practice Recommendations

The AES scheme retains Euler-type bias for two log-price integrals, so in the presence of extreme “vol-of-vol” or correlations, smaller time steps ($\kappa_i > 0$7) are advised. The scheme generalizes to jump processes or further variance factors, with complexity scaling with the number of CIR draws. For out-of-the-money options with few time steps, relative error may be larger and can be mitigated by moderately increasing $\kappa_i > 0$8.

In summary, the Double Heston model, when paired with advanced simulation and calibration methods, enables a tractable, expressive framework for modeling and pricing a broad class of derivatives with features closely matching observed market phenomena in volatility surfaces and exotics pricing (Higgins, 1 Feb 2026, Dimitrov et al., 22 Dec 2025, Alaya et al., 28 Jan 2025).