Fact-checked by Grok 2 weeks ago

Heston model

The Heston model is a mathematical framework in that describes the joint dynamics of an asset's price and its instantaneous variance ( squared) as correlated processes, enabling the pricing of securities such as call options. Introduced by Steven L. Heston in , it assumes the asset price follows a driven by , while the variance process adheres to a mean-reverting square-root to ensure non-negativity, with a parameter capturing the leverage effect where negative asset returns coincide with increases. The model's core equations specify the asset price dynamics as dS_t = \mu S_t dt + \sqrt{v_t} S_t dZ_{1t} under the physical measure, where \mu is the drift, and the variance as dv_t = \kappa (\theta - v_t) dt + \sigma \sqrt{v_t} dZ_{2t}, with \kappa > 0 as the speed of mean reversion, \theta > 0 as the long-term variance mean, \sigma > 0 as the volatility of variance, and \rho as the correlation between the Wiener processes dZ_{1t} and dZ_{2t}. Under the risk-neutral measure, the drift \mu shifts to the risk-free rate r, facilitating a semi-closed-form solution for option prices via the characteristic function and Fourier inversion: the call price is C(S_0, K, T) = S_0 P_1 - K e^{-rT} P_2, where P_1 and P_2 are risk-neutral probabilities computed numerically. This formulation addresses limitations of the constant-volatility Black-Scholes model by generating volatility smiles, skewness, and excess kurtosis in return distributions, which align better with empirical option data. Key advantages include its ability to incorporate arbitrary correlations between spot returns and volatility changes, thus explaining strike-price biases and return asymmetries observed in markets, while remaining computationally tractable for to observed prices. Originally applied to options, the model extends to and options by integrating stochastic interest rates, and it has influenced subsequent extensions like jump-diffusion models for enhanced realism in volatile markets. Despite requiring numerical methods for the integrals in P_1 and P_2, its affine structure allows efficient Fourier-based implementations, making it a cornerstone for quantitative practitioners.

Introduction and Background

Overview

The Heston model is a mathematical framework in that models the dynamics of an asset price under , where the instantaneous variance of returns follows a mean-reverting square-root known as the Cox-Ingersoll-Ross (CIR) process. This approach allows the to fluctuate randomly over time, providing a more realistic representation of market behavior compared to models with fixed . Developed to improve upon the Black-Scholes model's assumption of constant volatility—which fails to replicate the empirical and skew in option prices across different strikes and maturities—the Heston model incorporates stochastic variance to generate these observed patterns through the interplay of mean reversion and . The leverage effect, captured by a negative between asset returns and variance changes, further enhances its ability to fit market-implied volatilities. At its core, the Heston model is governed by coupled equations describing the joint evolution of the asset price and its variance, enabling the derivation of option prices via methods. Key parameters include the initial variance v_0, which sets the starting level of ; the long-term variance \theta, representing the equilibrium mean; the mean reversion speed \kappa, controlling how quickly variance returns to its long-term level; and the volatility of \sigma, which measures the in the variance process itself. The model finds broad applications in pricing derivatives across asset classes, including equity options, foreign exchange (FX) contracts, and interest rate instruments, due to its flexibility in capturing volatility dynamics relevant to these markets.

Historical Development

The model was introduced by Steven L. Heston in through his seminal paper titled "A Closed-Form Solution for Options with with Applications to Bond and Currency Options," published in The Review of Financial Studies. In this work, Heston proposed a framework for pricing European call options on assets whose volatility follows a mean-reverting square-root process, deriving a closed-form solution via techniques that express option prices in terms of the of the log asset price. The model's development addressed key limitations of the Black-Scholes framework, which assumes constant and fails to capture observed volatility dynamics such as the and leverage effects in financial markets. Earlier stochastic volatility models, including the one by and in , had introduced time-varying volatility but often lacked tractable formulas, relying instead on numerical approximations. Heston's built on these foundations by ensuring the joint for asset price and variance admits an affine structure, enabling efficient computation while accommodating between asset returns and volatility shocks. Following its publication, the Heston model saw rapid adoption in quantitative finance practice during the late 1990s and beyond, becoming one of the most widely used models for option pricing and due to its balance of realism and computational feasibility. Its affine diffusion properties also exerted significant influence on the broader class of affine term structure models, providing a foundational example for pricing derivatives and extending to multi-factor settings as formalized in subsequent theoretical advances. While Heston himself did not publish major updates to the core model, it inspired numerous extensions, such as multifactor volatility variants that incorporate matrix Wishart processes for enhanced flexibility in capturing term structure dynamics.

Model Formulation

Asset Price Dynamics

The asset price dynamics in the Heston model are described by the stochastic differential equation (SDE) under the physical probability measure: dS_t = \mu S_t \, dt + \sqrt{v_t} \, S_t \, dW_t^S, where S_t denotes the price of the underlying asset at time t, \mu is the constant drift representing the expected return, v_t is the instantaneous variance process, and W_t^S is a standard Brownian motion driving the asset price innovations. This formulation extends the framework by incorporating through the term \sqrt{v_t}, which scales the diffusion component. Conditional on the realized path of v_t, the asset returns are log-normally distributed, reflecting the multiplicative nature of the shocks; the stochastic scaling by \sqrt{v_t} allows the model to capture time-varying volatility and empirical features such as in asset returns. The asset price Brownian motion W_t^S is correlated with the Brownian motion W_t^v of the variance process, characterized by the instantaneous correlation \rho = \frac{d\langle W^S, W^v \rangle_t}{dt}, where \rho \in (-1, 1). This correlation parameter introduces the leverage effect, typically with \rho < 0, whereby negative shocks to the asset price tend to coincide with increases in future volatility, consistent with empirical observations in equity markets. The model specifies initial conditions S_0 > 0 for the asset price and v_0 \geq 0 for the initial variance, with \mu calibrated to historical data under the physical measure. The derivation relies on and properties of stochastic integrals for multivariate diffusions. The variance process v_t exhibits mean reversion, linking the asset dynamics to a separate for (detailed in the Volatility Process section).

Volatility Process

The in the Heston model describes the evolution of the instantaneous variance v_t, which follows a –Ingersoll–Ross () to capture mean-reverting and non-negative values. The () for v_t is given by dv_t = \kappa (\theta - v_t) \, dt + \sigma \sqrt{v_t} \, dW_t^v, where \kappa > 0 denotes the speed of mean reversion, \theta > 0 is the long-term mean level of variance, \sigma > 0 represents the of , and W_t^v is a standard under the physical measure. This formulation, adapted from the context, allows to fluctuate realistically while reverting toward its equilibrium. The parameters carry clear economic interpretations in modeling financial volatility. The mean-reversion speed \kappa governs the persistence of volatility deviations, with larger \kappa implying quicker of shocks and less prolonged volatility regimes. The long-term variance \theta sets the baseline equilibrium level around which volatility oscillates, often calibrated to historical averages. The vol-of-vol parameter \sigma quantifies the intensity of volatility fluctuations, capturing phenomena like clustering where high-volatility periods beget further variability. These interpretations align with empirical observations of markets, where volatility exhibits persistence and bursts. Key properties of the CIR process ensure its suitability for variance modeling. It guarantees v_t > 0 almost surely when starting from a positive initial value, provided the Feller condition $2\kappa\theta > \sigma^2 holds; under this restriction, the process mean-reverts to \theta without reaching the origin, maintaining strict positivity. The process is ergodic under standard parameter assumptions, converging to a stationary gamma distribution with shape $2\kappa\theta / \sigma^2 and scale \sigma^2 / (2\kappa), which is equivalently a scaled non-central chi-squared distribution for transitional densities. This stationarity reflects long-run stability in volatility levels. Regarding boundary behavior, the origin acts as a : when the Feller is satisfied, zero is unattainable and entrance-regular, meaning the process cannot cross it and reflects instantaneously if approached. If violated ($2\kappa\theta < \sigma^2), zero becomes attainable but instantaneously reflecting, allowing brief touches without absorption, which preserves positivity while permitting extreme low-volatility states observed in calm markets. This behavior is analyzed via scale functions and speed measures in diffusion theory. The variance process operates independently of the asset price dynamics except for the correlation \rho between their driving Brownian motions, which introduces leverage effects without altering the marginal law of v_t. In simulations, the square-root diffusion term poses stability challenges, as naive Euler discretization can produce negative values, particularly when the Feller condition fails or time steps are large. To address this, exact simulation methods exploiting the non-central chi-squared transition density ensure positivity and accuracy, while approximate schemes like full truncation or Milstein variants enhance stability for Monte Carlo applications. The asset price process scales its diffusion by \sqrt{v_t} to incorporate this stochastic .

Pricing Framework

Risk-Neutral Measure

In derivative pricing within the Heston model, the risk-neutral measure \mathbb{Q} is employed to ensure no-arbitrage conditions, under which the discounted asset price process is a martingale. This measure adjusts the drifts of the stochastic processes while preserving their diffusion structures, allowing the fair price of a derivative to be the expected value of its discounted payoff under \mathbb{Q}. Under the physical measure \mathbb{P}, the asset price S_t evolves as dS_t = \mu S_t \, dt + \sqrt{v_t} S_t \, dW_t^{S,\mathbb{P}}, where \mu is the expected return, and the instantaneous variance v_t follows the dv_t = \kappa (\theta - v_t) \, dt + \sigma \sqrt{v_t} \, dW_t^{v,\mathbb{P}}, with \langle dW_t^{S,\mathbb{P}}, dW_t^{v,\mathbb{P}} \rangle = \rho \, dt. The transition to \mathbb{Q} is achieved via , which redefines the Brownian motions through a change-of-measure density that incorporates the market prices of risk. The stock price risk is priced by shifting the asset Brownian motion by (\mu - r)/\sqrt{v_t}, where r is the risk-free rate, while the volatility risk is priced via a premium \lambda(S_t, v_t, t) = \lambda v_t, leading to a shift of \lambda v_t / (\sigma \sqrt{v_t}) for the volatility Brownian motion. Due to the correlation \rho, the shifts interact, but the affine structure ensures the processes retain their form under \mathbb{Q}. Under \mathbb{Q}, the dynamics simplify to dS_t = r S_t \, dt + \sqrt{v_t} S_t \, dW_t^{S,\mathbb{Q}} and dv_t = \kappa^* (\theta^* - v_t) \, dt + \sigma \sqrt{v_t} \, dW_t^{v,\mathbb{Q}}, where \kappa^* = \kappa + \lambda, \theta^* = \kappa \theta / (\kappa + \lambda), and \langle dW_t^{S,\mathbb{Q}}, dW_t^{v,\mathbb{Q}} \rangle = \rho \, dt. The volatility process remains unchanged in its diffusion term and correlation, as the model is affine, facilitating equivalent risk-neutral pricing. This measure change implies that the forward price of the asset is S_0 e^{rT} = \mathbb{E}^{\mathbb{Q}}[S_T], enabling consistent no-arbitrage valuation of options and other derivatives by taking expectations under \mathbb{Q}.

Characteristic Function Derivation

The Heston model exhibits an affine structure, where the log-asset price \log S_T and the variance v_T jointly form an affine diffusion process under the risk-neutral measure. This property implies that the characteristic function of \log S_T, defined as \phi(u) = \mathbb{E}[\exp(i u \log S_T) \mid \log S_t, v_t], admits an exponential-affine form:
\phi(u) = \exp\left( C(\tau, u) + D(\tau, u) v_t + i u (\log S_t + (r - q)\tau) \right),
where \tau = T - t is the time to maturity, r is the risk-free rate, and q is the dividend yield.
The derivation of this characteristic function proceeds via the Feynman-Kac theorem, which links the expectation to the solution of a partial differential equation (PDE) satisfied by \phi(u). Under the risk-neutral dynamics, the PDE for \phi is obtained from the infinitesimal generator of the joint process for (\log S_t, v_t), incorporating the stochastic differential equations for the asset price and variance. Assuming an affine ansatz for \phi, the PDE reduces to a system of ordinary differential equations (ODEs) for the functions C(\tau, u) and D(\tau, u). The resulting Riccati ODE for D(\tau, u) is
\frac{d D}{d \tau} = \frac{1}{2} \sigma^2 D^2 + (\rho \sigma i u - \kappa^*) D + \frac{1}{2} (u^2 - i u),
with initial condition D(0, u) = 0, where \kappa^* is the risk-neutral mean-reversion speed, \theta^* is the long-term variance under the risk-neutral measure, \sigma > 0 as the volatility of variance, and \rho as the between asset and variance shocks. The ODE for C(\tau, u) is
\frac{d C}{d \tau} = \kappa \theta D - i u (r - q),
with C(0, u) = 0. These equations arise directly from substituting the affine form into the PDE and collecting coefficients of v_t and the constant terms.
The for D admits a closed-form :
D(\tau, u) = \frac{r_1 - \rho \sigma i u + d}{\sigma^2} \left( 1 - \frac{1 - g e^{d \tau}}{1 - g} \right),
where d = \sqrt{ (\kappa^* - \rho \sigma i u)^2 + \sigma^2 (u^2 - i u) }, r_1 = \kappa + \lambda (with \lambda the market price of risk, often set to zero), and g = \frac{r_1 - \rho \sigma i u + d}{r_1 - \rho \sigma i u - d}. Integrating the ODE for C yields
C(\tau, u) = (r - q) i u \tau + \frac{\kappa \theta}{\sigma^2} \left[ (r_1 - \rho \sigma i u + d) \tau - 2 \log \left( \frac{1 - g e^{d \tau}}{1 - g} \right) \right]. Due to the complex-valued arguments, the expressions involve branch cuts in the , particularly along the imaginary axis for u, which must be handled carefully in numerical implementations to ensure and avoid discontinuities.
In the special case where \sigma = 0, the variance process becomes deterministic, v_t = \theta^* for all t assuming v_0 = \theta^*, and the characteristic function reduces to the Black-Scholes form: \phi(u) = \exp\left( i u (\log S_t + (r - q - \theta^*/2) \tau) - \frac{1}{2} u^2 \theta^* \tau \right). This limit confirms the model's consistency with the constant-volatility case.

Numerical Methods

Fourier-Based Pricing

The Fourier-based pricing approach in the Heston model leverages the analytically derived characteristic function to compute European option prices through inverse Fourier transforms, enabling a semi-closed-form solution. The price of a European call option with strike K and time to maturity \tau is expressed as C(K) = S_0 \Pi_1 - K e^{-r\tau} \Pi_2, where S_0 is the current asset price, r is the risk-free rate, and \Pi_1 and \Pi_2 represent risk-neutral probabilities obtained via inverse Fourier transforms of adjusted versions of the characteristic function \phi(u). Specifically, \Pi_j = \frac{1}{2} + \frac{1}{\pi} \int_0^\infty \Re \left[ e^{-i u \ln K} \frac{\phi_j(u - i \beta_j)}{i u} \right] du for j = 1, 2, with \beta_1 = 1 and \beta_2 = 0, and \phi_j denoting the characteristic function under suitable measure adjustments to ensure integrability. Heston (1993) introduced this method using a cosine transform representation, where the real part of the integrand is damped by the denominator i u and the adjustment in the argument of \phi, avoiding singularities at the origin while integrating over [0, \infty) to capture the probability measures. This formulation exploits the symmetry of the Fourier transform for log-prices, transforming the pricing integral into a form suitable for numerical evaluation without direct simulation of paths. In practice, the integrals are evaluated using rules such as Gauss-Laguerre, which are well-suited to the semi-infinite domain and the exponential decay of the integrand under the Heston characteristic function. Care must be taken to handle branch cuts in \phi(u) arising from the in its expression, typically by selecting the principal branch and ensuring the integration path avoids discontinuities. For enhanced efficiency with multiple strikes, extensions like the (FFT) can be applied after damping the payoff with a factor e^{\alpha k} (where \alpha > 0 ensures square-integrability), though the original Heston integrals suffice for single-strike computations. This approach offers significant advantages for pricing options, delivering rapid and accurate results—often in milliseconds—while inherently capturing skew through the -of- \sigma, which influences the between asset returns and . Unlike constant- models, it generates realistic smiles observed in and FX markets. Error analysis focuses on truncation of the integration upper limit and the choice of damping or adjustment parameters, with convergence controlled by balancing oscillatory decay and numerical stability; for instance, truncating at u_{\max} = 100 typically yields errors below $10^{-4} for at-the-money options under standard Heston parameters, improvable via adaptive quadrature.

Monte Carlo Implementation

The Monte Carlo implementation of the Heston model relies on simulating discrete paths of the asset price S_t and variance v_t to estimate expectations for option pricing, particularly for payoffs beyond European calls where semi-closed forms are unavailable. This approach is versatile for handling path-dependent or early-exercise features, such as American options or barrier derivatives, by generating multiple sample paths under the risk-neutral measure and averaging the discounted payoffs. Discretization schemes approximate the continuous-time stochastic differential equations (SDEs) of the model. The Euler-Maruyama method is standard for the asset price dynamics in logarithmic form: \ln S_{t_{i+1}} = \ln S_{t_i} + \left( r - \frac{v_{t_i}}{2} \right) \Delta t + \sqrt{v_{t_i}} \sqrt{\Delta t} \, Z_i^S, where Z_i^S is a standard normal random variable, and dividends are omitted for simplicity. For the variance process, the Milstein scheme improves accuracy by including a second-order term to handle the \sqrt{v_t} diffusion coefficient: \begin{aligned} v_{t_{i+1}} &= v_{t_i} + \kappa (\theta - v_{t_i}) \Delta t + \sigma_v \sqrt{v_{t_i}} \sqrt{\Delta t} \, Z_i^v + \frac{\sigma_v^2}{4} ((\sqrt{\Delta t} \, Z_i^v)^2 - \Delta t), \\ Z_i^v &\sim \mathcal{N}(0,1). \end{aligned} These schemes can produce negative variance values, especially when the Feller condition $2\kappa\theta > \sigma_v^2 is violated, leading to bias in path simulations. To enforce positivity, the full truncation scheme modifies the Euler discretization by incorporating the positive part operator directly: v_{t_{i+1}} = v_{t_i} + \kappa (\theta - \max(v_{t_i}, 0)) \Delta t + \sigma_v \sqrt{\max(v_{t_i}, 0)} \sqrt{\Delta t} \, Z_i^v, and for the asset: \ln S_{t_{i+1}} = \ln S_{t_i} + \left( r - \frac{\max(v_{t_i}, 0)}{2} \right) \Delta t + \sqrt{\max(v_{t_i}, 0)} \sqrt{\Delta t} \, Z_i^S. This minimizes bias compared to reflection or absorption methods. This truncation introduces a small upward bias but outperforms basic Euler in root-mean-squared error for practical time steps like \Delta t = 1/252. Exact simulation avoids discretization bias entirely by sampling directly from the conditional distributions of the SDEs. For the CIR variance process, transitions v_{t + \Delta t} | v_t follow a scaled non-central chi-squared distribution: v_{t + \Delta t} = \frac{4 \kappa \theta}{\sigma_v^2 (1 - e^{-\kappa \Delta t})} \chi'^2_{d} \left( \frac{4 \kappa e^{-\kappa \Delta t}}{\sigma_v^2 (1 - e^{-\kappa \Delta t})} v_t \right), where d = 4 \kappa \theta / \sigma_v^2 is the degrees of freedom, and \chi'^2_d(\lambda) denotes the non-central chi-squared with non-centrality \lambda. Sampling involves generating the chi-squared via sums of squared normals or, for low d, a Poisson-augmented method to handle the distribution near zero, naturally respecting the Feller condition without truncation. Given the integrated variance \int_t^{t+\Delta t} v_s \, ds, the asset price increment is lognormal conditional on the variance path, allowing direct sampling of S_{t + \Delta t} using the model's correlation structure. This method enables "almost exact" schemes by applying exact steps over larger intervals, reducing computational steps while controlling bias to machine precision. To incorporate the leverage effect \rho, joint Brownian motions are generated using : Z_i^v = \rho Z_i^S + \sqrt{1 - \rho^2} Z_i^\perp, where Z_i^\perp \sim \mathcal{N}(0,1) is independent of Z_i^S. This ensures the correct in each time step for both and exact methods. techniques enhance efficiency by lowering the number of paths needed for a given . Antithetic variates pair simulations with negated random variables (-Z_i^S, -Z_i^v) to induce negative in payoffs, halving variance in symmetric cases. use the Black-Scholes formula as a baseline, since its expectation is known analytically; the estimator adjusts Heston payoffs by c (BS(S, K, T, r, \sigma_{imp}) - \widehat{BS}), where \sigma_{imp} is an implied volatility from and c is an optimal coefficient estimated from a pilot run, yielding variance reductions up to 90% for at-the-money options. These methods are particularly effective when combined with over the parameter. Computationally, requires balancing time steps \Delta t (e.g., T/100 to T/1000 for daily rebalancing) and count (typically $10^5 to $10^6 for 95% confidence intervals under 1% width). schemes tolerate larger steps (up to \Delta t = 0.1) without , while discretized methods demand finer grids when Feller is violated to limit negativity rates below 1%. Overall runtime scales as O(N M), where N is paths and M = T/\Delta t is steps, making GPU acceleration viable for high-dimensional extensions.

Parameter Estimation

Calibration Techniques

Calibration of the Heston model parameters—namely, the mean reversion rate κ, long-term variance θ, volatility of variance σ, ρ between asset returns and volatility, and initial variance v₀—typically aims to minimize the discrepancy between theoretical option prices or implied volatilities generated by the model and those observed in the market across a range of strikes and maturities. This process is essential for aligning the model's dynamics with empirical option pricing surfaces, enabling accurate forecasting and risk management. The most common approach employs least-squares optimization, where the objective function is the sum of squared errors between market and model implied volatilities, often weighted to emphasize relative errors or sensitivity via : \min_{\Theta} \sum_{t,k} w_{t,k} \left( \sigma_{t,k}^{\text{market}} - \sigma_{t,k}^{\text{model}}(\Theta) \right)^2, with Θ denoting the parameter vector, t the maturity, k the strike, and w_{t,k} optional weights such as squared for better . Alternative formulations include relative mean squared errors on prices or implied volatilities to account for the scale of option values, reducing bias in fitting deep out-of-the-money options. These objectives capture the model's ability to reproduce the and skew observed in equity or options markets. Due to the non-convex nature of the optimization landscape, a variety of algorithms are employed. Local methods such as or the Levenberg-Marquardt algorithm, which combines with Gauss-Newton updates for robustness near singularities, are popular for their efficiency when analytical gradients are available—derived from the model's . For global exploration to avoid local minima, evolutionary algorithms like are used, particularly in high-dimensional parameter spaces, with settings such as population sizes of 200–300 and mutation factors around 0.8 to balance convergence speed and solution quality. Parameter constraints are strictly enforced during optimization, including positivity for κ, θ, σ, and v₀ (typically bounded as 0 < κ ≤ 20, 0 < θ ≤ 1, 0 < σ ≤ 5, 0 < v₀ ≤ 1), |ρ| < 1 (often -0.9 ≤ ρ ≤ -0.1 for equity skew), and the Feller condition 2κθ ≥ σ² to ensure the CIR volatility process remains positive. The calibration process proceeds iteratively: initial parameter guesses are derived from historical data moments, such as setting v₀ to the squared realized volatility, κ and θ to match the sample mean reversion and long-term variance of returns, σ to the volatility of squared returns, and ρ to the correlation between returns and volatility changes. These starting values, or alternatives from asymptotic expansions of implied volatilities in short- or long-time regimes, provide a feasible point within constraints. Optimization then proceeds by evaluating the objective function repeatedly, adjusting parameters until convergence criteria are met, such as a residual norm below 10^{-10} or a maximum iteration limit, often requiring 10–20 evaluations per local minimum. This step-by-step fitting ensures the parameters reflect market-implied dynamics while maintaining model admissibility.

Empirical Challenges

One key empirical challenge in applying the Heston model arises during parameter calibration, where the four core parameters (κ, θ, σ, and ρ) along with the initial variance v₀ can effectively fit short-term implied volatility smiles observed in market data for . However, this fitting often leads to overfitting, resulting in poor extrapolation to longer maturities, as the mean-reverting structure of the variance process fails to capture the persistent skew and term structure dynamics beyond one to two years. Another practical issue stems from violations of the Feller condition (2κθ > σ²), which, while not strictly required for model validity, frequently occur in calibrations to real to match observed behaviors. Such violations allow the variance process to reach zero or negative values in simulations, necessitating schemes like full truncation or to ensure non-negativity during implementations, which can introduce bias and reduce simulation efficiency. The Heston model operates in an incomplete due to the unhedgeable , requiring of a via the λ, which adjusts the risk-neutral long-run variance θ* to differ from the physical measure θ (specifically, θ* = κθ / (κ + λ) when λ ≠ 0). This discrepancy complicates the separation of physical dynamics from premia in empirical , often leading to unstable identification when fitting to both historical returns and option prices simultaneously. Post-2008 analyses revealed the model's underperformance in capturing extreme market events, as its pure diffusion-based process could not adequately replicate the sudden spikes in implied volatilities and fat-tailed return distributions observed during turbulent periods like the subprime crisis. Comparative studies showed that jump-augmented extensions, such as the Bates model incorporating in the asset price, significantly outperformed the standard Heston in pricing and hedging short-term out-of-the-money options amid high regimes, highlighting the need for components to address crisis dynamics. In contemporary applications, the Heston model exhibits gaps in handling stochastic interest rates, assuming constant rates that limit its realism for long-dated derivatives in low-rate environments, prompting extensions like the Heston-Hull-White hybrid to incorporate rate volatility. Similarly, its single-factor structure struggles with multi-asset correlations, underestimating cross-asset volatility spillovers in portfolios, which multifactor Wishart-based variants address more effectively. Furthermore, post-2010 empirical evidence on realized volatility paths has favored rough volatility alternatives, such as the rough Heston model, which better captures the empirically observed roughness ( around 0.1) in volatility time series compared to the smoother paths generated by the classical Heston process.