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

• 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. – JohnDoe Apr 6 at 9:04
• @johndoe, looking at the code, you're correct. – will Apr 7 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. – Daniele Rocchi Apr 7 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).