Numerical simulation of Heston model

Question

I am trying to simulate on Python random paths for a general asset price as described by the Heston model:

\begin{equation} \begin{aligned} dS_t &= \mu S_t dt + \sqrt{\nu_t} S_t dW^S_t \\ d\nu_t &= \kappa(\theta - \nu_t) dt + \xi \sqrt{\nu_t} dW^{\nu}_t \\ \textrm{Corr}[W^S_t, W^{\nu}_t] &= \rho \end{aligned} \end{equation}

where:

• $W_t^{S}$ and $W_t^{\nu}$ are two standard Brownian motions with correlation $\rho$.
• $\nu _{t}$ is the instantaneous variance.
• $\mu$ is the rate of return of the asset.
• $\theta$ is the long variance.
• $\kappa$ is the rate at which $\nu_t$ reverts to $\theta$.
• $\xi$ is the volatility of the instantaneous volatility.

Hence I implemented the following function:

def HeMC (S0, mu, v0, rho, kappa, theta, xi, T, dt):

    # Generate a Monte Carlo simulation for the Heston model

    # Generate random Brownian Motion
    MU  = np.array([0, 0])
    COV = np.matrix([[1, rho], [rho, 1]])
    W   = np.random.multivariate_normal(MU, COV, T)
    W_S = W[:,0]
    W_v = W[:,1]

    # Generate paths
    vt    = np.zeros(T)
    vt[0] = v0
    St    = np.zeros(T)
    St[0] = S0
    for t in range(1,T):
        vt[t] = np.abs(vt[t-1] + kappa*(theta-np.abs(vt[t-1]))*dt + xi*np.sqrt(np.abs(vt[t-1]))*W_v[t])
        St[t] = St[t-1]*np.exp((mu - 0.5*vt[t])*dt + np.sqrt(vt[t]*dt)*W_S[t])

    return St, vt

The problem is that when I run this function with the following parameters I often obtain paths which seems not to make sense. Especially the instantaneous volatility does not seem mean-reverting but it often follows wild paths.

T     = 252
dt    = 1/252
S0    = 100 # Initial price
mu    = 0.1 # Expected return
sigma = 0.2 # Volatility
rho   = -0.2 # Correlation
kappa = 0.3 # Revert rate
theta = 0.2 # Long-term volatility
xi    = 0.2 # Volatility of instantaneous volatility
v0    = 0.2 # Initial instantaneous volatility

I believe the problem is in the way I discretised the process or maybe there are some bugs in my code but I could not find anything.

Thanks for your help.

Answer 1

There are two mistakes in the code:

1) In the line

vt[t] = np.abs(vt[t-1] + kappa*(theta-np.abs(vt[t-1]))*dt + xi*np.sqrt(np.abs(vt[t-1]))*W_v[t])

you forgot to multiply W_v[t] by np.sqrt(dt).
This is the reason the volatility increases so much.

2) The line

St[t] = St[t-1]*np.exp((mu - 0.5*vt[t])*dt + np.sqrt(vt[t]*dt)*W_S[t])

should be

St[t] = St[t-1]*np.exp((mu - 0.5*vt[t-1])*dt + np.sqrt(vt[t-1]*dt)*W_S[t])

Also there is no need to use three times the np.abs function. One is enough. (The more external).

Answer 2

The variance process in the Heston model (i.e., CIR process) is notorious for the simulation. A typical Euler/Milstein scheme ends up with negative variance in a high portion of paths. Making it positive by np.abs() is distorting the true distribution.

Consider Andersen (2008)'s QE scheme to resolve this issue:

• Andersen L (2008) Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189

Alternatively, consider the noncentral chi-squared distribution to simulate vt, which is the analytic property of the CIR process. The code below is for idea only. (Please modify properly if there's syntax error.)

chi_df = 4 * theta * kappa / xi**2
exph = np.exp(-kappa*dt/2)
phi = 4*kappa / xi**2 / (1/exph - exph)
chi_nonc = vt[t-1] * exph * phi
vt[t] = (exph / phi) * np.random.noncentral_chisquare(df=chi_df, nonc=chi_nonc)