Multilevel Monte Carlo simulation for the Heston stochastic volatility model

Abstract

We combine the multilevel Monte Carlo (MLMC) method with a numerical scheme for the Heston model that simulates the variance process exactly or almost exactly and applies the stochastic trapezoidal rule to approximate the time-integrated variance process within the SDE of the logarithmic asset process. We conduct separate simulations for path-independent options and path-dependent options. In both situations, novel MLMC estimators are established, and the theoretical convergence rates are derived for the full parameter regime. We present numerical results to demonstrate the efficiency of our MLMC estimators.

Appendix. The NCI method

We revisit the NCI method by Van Haastrecht and Pelsser [40]. Let $Q_{max}$ be a positive integer, $\mathcal {Q}:=\{0,1,2,...,Q_{max}\}$ be a set of some Poisson-values, and $\mathcal {U}:=\{0,...,1-\delta \}\in [0,1)$, in which $\delta$ is a small number such as $10^{-15}$. The NCI method starts with a precalculation of the inverse of the chi-square distributions, i.e., it calculates

$$\begin{aligned} H_{Q}^{-1}(U):=G_{\mathcal {X}_{d+2Q}^{2}}^{-1}(U),\quad \forall Q\in \mathcal {Q}, \forall U\in \mathcal {U} \end{aligned}$$

(28)

where $G_{\mathcal {X}_{d+2Q}^{2}}^{-1}$ is the inverse of the chi-square distribution with $d+2Q$ degrees of freedom. Here, $d=\frac{4\theta k}{\sigma ^{2}}$. Recall that for the variance process, $V_{t}$ at time t follows the scaled noncentral chi-square distribution

$$ C_{0}\chi _{d}^{2}(\lambda ) $$

given $V_{u}$ for any $u<t$, with $C_{0}=\frac{\sigma ^{2}(1-\textrm{e}^{-k(t-u)})}{4k}$ and $\lambda =\frac{4k\textrm{e}^{-k(t-u)}}{\sigma ^{2}(1-\textrm{e}^{-k(t-u)})}V_{u}$. The NCI method then generates one sample Q from the Poisson distribution with the mean $\frac{\lambda }{2}$ and another sample U from the uniform distribution on (0, 1). A sample of $V_{t}$ is generated by

$$\begin{aligned} F_{Q}^{-1}(U):={\left\{ \begin{array}{ll} C_{0}J(U), &{} Q\le Q_{max}\\ C_{0}F_{\mathcal {X}_{d+2Q}^{2}}^{-1}(U), &{} Q> Q_{max} \end{array}\right. } \end{aligned}$$

(29)

where J(.) represents an interpolation rule based on $H_{Q}^{-1}(.)$, and $F_{\mathcal {X}_{d+2Q}^{2}}^{-1}(.)$ is the inverse of the noncentral chi-square cumulative distribution function.

The authors suggested two interpolations: the linear interpolation and the monotone cubic Hermite spline interpolation. Here, we prefer the latter, since it is much more accurate with the same number of points in the cache. We summarize the algorithm of the NCI method as Algorithm 2.

Algorithm 2

NCI method

As we can see, the accuracy of the NCI method relies on the parameter $Q_{max}$ and the size of the uniform grid $\mathcal {U}$. In our numerical test, for very high precision, we set $Q_{max}=20$ and $m=10000$, where m is the size of $\mathcal {U}$, such that $\mathcal {U}=\{0,1/m,2/m.,..,(m-1)/m,1-\delta \}$, with $\delta =10^{-15}$.