# Numerical simulation of Heston model

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 = v0
St    = np.zeros(T)
St = 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.

• Without looking at your code, the initial volatility and the long run mean seems to be pretty high. Remember, that you model the variance $v_t$ and not the standard deviation, so try to use 0.04 as the respective values and see what happens. Apr 6, 2019 at 9:04
• @johndoe, looking at the code, you're correct.
– will
Apr 7, 2019 at 9:30
• @JohnDoe yes you are right, I badly calibrated the model. Recalling that its actually variance and not volatility, I adjusted the values for v0, theta and kappa and now works perfectly fine. Apr 7, 2019 at 22:40

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).

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:

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)