# Simulation of Heston model, best reference?

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.

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