I am currently experimenting with various implementations for simulating the standard Heston model. \begin{eqnarray*} dS_t &=& \mu S_t \, dt + \sqrt{v_t} \cdot S_t \, dW_t^S \\ dv_t &=& \kappa(\theta - v_t) \, dt + \xi \cdot \sqrt{v_t} \, dW_t^v, \end{eqnarray*} where the correlation between the Brownian motions is $\rho$.
I am however struggling to find a decent reference article with an implementation which is accurate for all choices of parameter values.

I have, for example, implemented the method described in the article "A Simple and Exact Simulation Approach to Heston Model" by J. Zhu. This has the advantage of being very easy to implement and understand. It also give good results even for higher values of the correlation parameter. It is also very fast.

However, when the "vol-vol", $\xi$, is large and the Feller condition $2 \kappa \theta > \xi^2$ is violated by a large margin, the method fails. The option prices seem to become too large in general. The reason this is happening is not too hard to understand. The method of Zhu is based on a moment matching procedure for the volatility process. When $\xi$ is too large the equations you need to solve to make the moments match lack solution. The authors "solved" this by flooring a negative value to zero. If values are only slightly negative the effect of this should not be too bad, but for larger negative values the error should be significant, which is exactly what is seen for larger $\xi$.

What is the current state of the art regarding to simulation of the Heston method? Are there any good references to point to? The most important thing for me is that the method produces at least decently accurate results. After that, a faster method is of course preferable. Simplicity of implementation comes third.


2 Answers 2


This link presents several discretization schemes for Heston: https://www.degruyter.com/view/journals/math/15/1/article-p679.xml

For example, Milstein is a popular one.

An alternative to the Euler discretization scheme for the Heston model is the second-order discretization method. The system of SDE under the risk-neutral measure \begin{eqnarray*} dS_t &=& r S_t \, dt + \sqrt{v_t} S_t \, dW_t^S \\ dv_t &=& \kappa(\theta - v_t) \, dt + \sqrt{v_t}(\xi_1\, dW_t^S+\xi_2 \, dW_t^v), \end{eqnarray*} gets discretized as follows: \begin{eqnarray*} dS_{i+1} &=& S_i\left(1+rh+\sqrt{v_i}\Delta W^S\right)+\frac{1}{2}r^2S_ih^2 \\ &+&\left(\left[r+\frac{\xi_1-\kappa}{4}\right]S_i\sqrt{v_i}+\left[\frac{\kappa\theta}{4}-\frac{\xi^2}{16}\right]\frac{S_i}{\sqrt{v_i}}\right)\Delta W^Sh \\ &+&\frac{1}{2}S_i\left(v_i+\frac{\xi_1}{2}\right)((\Delta W^S)^2-h)+\frac{1}{4}\xi_2S_i(\Delta W^v\Delta W^S+\varepsilon) \\ v_{i+1} &=& \kappa\theta h+(1-\kappa h)v_i + \sqrt{v_i}(\xi_1\Delta W^S+\xi_2\Delta W^v)-\frac{1}{2}\kappa^2(\theta-v_i)h^2 \\ &+&\left(\left[\frac{\kappa\theta}{4}-\frac{\xi^2}{16}\right]\frac{1}{\sqrt{v_i}} - \frac{3\kappa}{2}\sqrt{v_i}\right)(\xi_1\Delta W^S+\xi_2\Delta W^v)h \\ &+&\frac{1}{2}\xi_1^2((\Delta W^S)^2-h)+\frac{1}{4}\xi_2^2((\Delta W^v)^2-h)+\frac{1}{2}\xi_1\xi_2\Delta W^S\Delta W^v \end{eqnarray*} where $\xi^2 = \xi_1^2+\xi_2^2$ and $$\varepsilon = \begin{cases} h, & \mbox{with prob. } \frac{1}{2} \\ -h, & \mbox{with prob. } \frac{1}{2} \end{cases}$$ $\varepsilon$ and $\Delta W^S$ being independent. One can consider taking the absolute value of $v_i$.

This scheme can be used for instance by fixing $h=\frac{T}{n}$ with various sample sizes $n$ and 1e+06 repetition for each $n$. This method is known to produce smaller estimation bias but it has a slightly worse convergence

  • $\begingroup$ Thanks, I will implement that one as a reference. The most important thing for me is that it can be trusted to be bias free, or close enough. I am aware of the Broadie-Kaya exact scheme, which is theorethically bias free. But it seemed to be very calculation intensive and probably quite slow. I have implemented the Andersen QE-scheme with martingale corrections which seems pretty good in most cases. It seems a bit worse when the vol-vol is large and the correlation is positive though. $\endgroup$ Oct 22, 2020 at 16:04
  • $\begingroup$ Also, one could try to fine tune the Crack-Nicholson scheme, but I think the 2nd order method is very good in terms of bias error. In the Glaserman's book you can find even more details. Another approach would be to simulate the cosine method for Heston. Please let me know how it goes with it and if you are happy, feel free to accept this answer. $\endgroup$
    – FunnyBuzer
    Oct 22, 2020 at 16:17

Here are some references to the state-of-art simulation schemes of the Heston model. Roughly speaking there are two lines of algorithms: time-discretization v.s. exact schemes.

Time-discretization schemes:

  • Andersen (2008)'s QE method
  • Tse & Wan (2013)'s Inverse-Gaussian approximation

The computation cost for jumping a step is cheap in this line of methods, but you have to make multiple jumps (with a small-time step) in order to control the error. Therefore, these are more suitable if you need the price time series (i.e., path) as in pricing path-dependent (e.g., Asian or barrier) options.

See Van Haastrecht & Pelsser (2010) for performance comparison. They propose their own methods, but the conclusion is that the QE method is practically the best. Tse & Wan (2013) is relatively new. The cost for one jump is higher than the QE, but it's more accurate, so you don't need as many jumps as in QE. So, the overall cost can be lower.

Exact simulation schemes:

  • Broadie & Kaya (2006)'s original exact simulation scheme.
  • Glasserman & Kim (2011)'s Gamma series expansion.

These are more suitable when you only care about the terminal price (not the path). In these methods, you can jump any time step albeit with high computation cost. Regarding this line, read another post of mine



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.