Simulating Heston Variance with Correct Time Scaling and QE Methods

# Numerical simulation of Heston model


# 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[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 by Cantaro (score 5)

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

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 by jaehyukchoi49 (score 1)

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

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