Implementations - 2023.1 English

The discrete form (Euler-Maruyama Form) of Heston Model is:

\[\ln{S(j\Delta)} = \ln{S((j-1)\Delta)} + \mu\Delta + \sqrt[2]{\nu((j-1)\Delta)} \Delta W_t^S\]

\[\nu(j\Delta) = \nu((j-1)\Delta) + \kappa(\theta - \nu((j-1)\Delta))\Delta + \sigma\sqrt[2]{\nu((j-1)\Delta)}\Delta W_t^\nu\]

Where \(\Delta\) stands for unit timestep length. \(S(j\Delta)\) stands for S in j th timesteps, AKA \(S(j * \Delta t)\). \(\nu(j\Delta)\) has a similar meaning. To simplify the process, we use \(\ln{S(j\Delta)}\) instead of \(S\) since multiplication becomes addition after logarithm.

\(\Delta W_t^S\) and \(\Delta W_t^\nu\) could be calculated by:

\[\Delta W_t^S = \Delta Z_S\]

\[\Delta W_t^\nu = \Delta Z_\nu\]

Where \(Z_S\) and \(Z_\nu\) are two uniform distributed random numbers that have correlation \(\rho\)

Although \(\nu\) is non-negative in continuous form, it would become negative if we use the Euler-Maruyama form above directly. There are several variations to solve this issue, here we provide 5 of most commonly used.

kDTReflection: use absolute value of the volatility of the last iteration.

\[\ln{S(j\Delta)} = \ln{S((j-1)\Delta)} + \mu\Delta + \sqrt[2]{|\nu((j-1)\Delta)|} \Delta W_t^S\]

\[\nu(j\Delta) = |\nu((j-1)\Delta)| + \kappa(\theta - |\nu((j-1)\Delta)|)\Delta + \sigma\sqrt[2]{|\nu((j-1)\Delta)|}\Delta W_t^\nu\]

kDTPartialTruncation: use only absolute value of volatility of the last iteration in sqrt root part.

\[\ln{S(j\Delta)} = \ln{S((j-1)\Delta)} + \mu\Delta + \sqrt[2]{|\nu((j-1)\Delta)|} \Delta W_t^S\]

\[\nu(j\Delta) = \nu((j-1)\Delta) + \kappa(\theta - \nu((j-1)\Delta))\Delta + \sigma\sqrt[2]{|\nu((j-1)\Delta)|}\Delta W_t^\nu\]

kDTFullTruncation: use only positive part or zero of volatility of the last iteration.

\(\nu((j-1)\Delta)^+ = \nu((j-1)\Delta)\) if \(\nu((j-1)\Delta) > 0\) \(\nu((j-1)\Delta)^+ = 0\) if \(\nu((j-1)\Delta) \leq 0\).

\[\ln{S(j\Delta)} = \ln{S((j-1)\Delta)} + \mu\Delta + \sqrt[2]{\nu((j-1)\Delta)^+} \Delta W_t^S\]

\[\nu(j\Delta) = \nu((j-1)\Delta) + \kappa(\theta - \nu((j-1)\Delta)^+)\Delta + \sigma\sqrt[2]{\nu((j-1)\Delta)^+}\Delta W_t^\nu\]

kDTQuadraticExponential and kDTQuadraticExponentialMartingale are more accurate variation, details could be found in reference papers [ANDERSON2005]. They use a different approximation method to calculate \(\nu\), here’s brief on its algorithm

Step 1: we calculate first order moment and second order moment of \(\nu\).

\[m = \theta + (\nu((j-1)\Delta) - \theta)e^{-\kappa \theta}\]

\[s^2 = \frac{\nu((j-1)\Delta)\sigma^2 e^{-\kappa \Delta} }{\kappa} (1 - e^{-\kappa \Delta}) + \frac{\theta \sigma^2}{2 \kappa}(1 - e^{-\kappa \Delta})^2\]

Step 2: Calculate \(\Psi = s^2 / m^2\)

Step 3: If \(\Psi \leq \Psi_{sw}, \Psi_{sw} = 1.5\), Then

Step 3.1: Calculate \(a\) and \(b^2\)

\[b^2 = 2\Psi^{-1} - 1 + \sqrt[2]{2\Psi^{-1}} \sqrt[2]{2\Psi^{-1}-1}\]

\[a = \frac{m}{1+b^2}\]

Step 3.2: Calculate \(\nu(j\Delta)\)

\[\nu(j\Delta) = a(b+Z_\nu)^2\]

Step 4: If Step 3 does not hold, Then

Step 4.1: Calculate \(\beta\) and \(p\)

\[p = \frac{\Psi - 1}{\Psi + 1}\]

\[\beta = \frac{2}{m(\Psi + 1)}\]

Step 4.2: Calculate \(U_\nu = \Psi(Z_\nu)\)

Step 4.3: Calculate \(\nu(j\Delta)\)

\[\nu(j\Delta) = 0 \:\: if \:\: U_\nu \leq p\]

\[\nu(j\Delta) = \frac{1}{\beta}(\frac{\log{(1-p)}}{1-U_\nu}) \:\: if \:\: p < U_\nu\]

It should be noticed that they both have two branches for value in different range. These two branches have a similar calculation process. Furthermore, only one branch is active at the same time. By merging these two branches into one branch and manually binding calculations to DSPs, it will cut off DSP cost. This won’t change its performance and accuracy.

In Monte Carlo Simulation, we need to compute stock prices of multiple paths at multiple time steps. Therefore we need two loops to calculate prices and volatilities, the inner loop is either timestep loop or path loop. Price at each time step is calculated using last time step’s price and volatility as input. And we use 1-D array to store price and volatility of each path’s history (last timestep).

If the inner loop is timestep loop, as red arrows demonstrate in the diagram below, it will keep update the same array element until reaches max timesteps. Such operation can not achieve initiation interval (II)=1 and will greatly slow down the calculation process. If the inner loop is path loop, as green arrows demonstrate in the diagram below, it will keep updating different array element each time. Such operation will avoid dependency issue and reach II=1, which is used in this implementation.