Simulating the Heston Process

Fourier inversion prices European options under Heston in milliseconds, but exotic payoffs — barriers, cliquets, path-dependent structures — require Monte Carlo simulation of the full $(S_t, v_t)$ path. Simulating the Heston variance process is non-trivial: the CIR dynamics $dv_t = \kappa(\theta - v_t)\,dt + \eta\sqrt{v_t}\,dW_t^v$ have non-Lipschitz diffusion coefficient $\sqrt{v_t}$, which causes the naive Euler scheme to produce negative variance values even when the Feller condition is satisfied (because the discrete step can overshoot zero).

This note covers three schemes of increasing sophistication: the naive Euler (for context only), the full-truncation Euler (the standard practical fix), the Milstein scheme (one order of convergence higher), and the exact Broadie-Kaya scheme for the joint distribution of $(v_t, v_T, \int_t^T v_s\,ds)$ when an exact sample is feasible.

Learning objectives. By the end of this note you should be able to: