Heston Monte Carlo Variance Reversion and Path Interpretation

# How to fix my Monte Carlo simulation?


# How to fix my Monte Carlo simulation?












I hope that you are all having a blessed day,

I am working on calibrating the Heston model to observed data, and one of the steps (as proposed by Mikhalov, 2003) is to compare the performance of my Heston pricing formula with the results of a Monte Carlo simulation.

I decided to use the following discretization to do so:

$ S_{t+dt} = S_t * exp((r-\frac{1}{2} * v_t)*dt + \sqrt{v_t * dt}*Z_s) $ (Milstein)

$ v_{t+dt} = (\theta + (v_t - \theta)*e^{-\kappa * dt}) * exp(-\frac{1}{2}*\Gamma_t^{2} + \Gamma_t*Z_v)$

$ \Gamma_t = ln(1 + \frac{\sigma^2v_t(1-e^{-2\kappa dt})}{2\kappa (\theta + (v_t - \theta)*e^{-\kappa * dt})^2})$

With the equations for $ v_{t+dt}$ and $\Gamma_t$ coming from Andersen Moment matching scheme (Coming from F. Rouah - The Heston model and and its extensions in Matlab and C#).

I tried to implement this in python with the following code:

```
# Parameters used only for implementation, will generate new parameters when checking the accuracy of my pricer
S0 =100
K = 100
v0 = 0.1
r = 0.02
kappa = 1.5768
theta = 0.0398
sigma = 0.3
rho = -0.5711
tau = 1
N = 252
M = 1000

def heston_model_sim(S0, v0, rho, kappa, theta, sigma, t, N, M):
# Initializing other parameters
dt = t/N
mu = np.array([0,0])
cov = np.array([[1,rho],[rho,1]])
g0 = np.log(1 + (sigma ** 2 * v0 * (1 - np.exp(-2 * kappa * dt)))/(2 * kappa * (theta + (v0 - theta) * np.exp(- kappa * dt)) ** 2))

# Arrays for storing prices and variances
S = np.full(shape = (N+1, M), fill_value = S0)
v = np.full(shape = (N+1, M), fill_value = v0)
gamma = np.full(shape = (N+1, M), fill_value = g0)

# Sampling correlated brownian motions under risk-neutral measure
Z = np.random.multivariate_normal(mu, cov, (N, M))

for i in range(1, N+1):
    S[i] = S[i-1] * np.exp((r - 0.5 * v[i-1]) * dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,1])
    gamma[i-1] = np.log(1 + (sigma ** 2 * v[i-1] * (1 - np.exp(-2 * kappa * dt)))/(2 * kappa * (theta + (v[i-1] - theta) * np.exp(- kappa * dt)) ** 2)) 
    v[i] = (theta + (v[i-1] - theta) * np.exp(- kappa * dt)) * np.exp(-0.5 * gamma[i-1]**2 + gamma[i-1] * Z[i-1,:,0])
Return S,v
```

(Remark - I am not the best with pasting code on stackexchange but the indexation of the return statement is right in my code.)

And this is the output that I get:

It seems like there is a (big) problem with my volatility process (maybe as well with stock process but it is less eye catching) but I can not seem to find where it comes from and how can I fix it.

Do you see any mistakes in my code or in the discretization I am using ?

Thank you so much in advance

Update:

this is what happens when I lower the initial volatility:

It seems to be more correct, but why is there a problem when starting with a higher volatility ?

## Answer by KaiSqDist (score 1, accepted)

https://quant.stackexchange.com/a/79626

The problem is that the Heston process describes variance as a mean-reverting process (towards theta).

Therefore, if you specify variance as higher than 0.0398 (the theta level it reverts to), the variance will revert to that level. That is why you see in your first variance plot it is trending down.

In your second plot you started at roughly 0.04, which is why the variance stays stochastic around that level.