Heston Model Simulation with Python

by John | May 14, 2023



The Heston model is a useful model for simulating stochastic volatility and its effect on the potential paths an asset can take over the life of an option. The Heston model also allows modeling the statistical dependence between the asset returns and the volatility which have been empirically shown to have an inverse relationship. This allows for more effective modeling than the Black-Scholes formula allows due to its restrictive assumption of constant volatility. 


Heston Model SDE

The SDE shown below is similar to the standard geometric Brownian motion apart from the fact that the volatility is a stochastic, mean-reverting process. The notation we use here is from Wikipedia


\(dS_t = \mu S_tdt +\sqrt{v_t}S_tdW_{t}^S\)

\(dv_t = \kappa(\theta - v_t)dt +\xi\sqrt{v_t}dW_{t}^v\)


\(dW_{t}^S\) and \(dW_{t}^v\) are Weiner processes with instantaneous correlation \(\rho\) 


\(\mu = \text{drift term}\\ \theta = \text{long run average variance}\\ \kappa= \text{rate of mean reversion}\\ \xi = \text{vol of vol}\\ \)


Generating Correlated Random Normal Variables


To simulate Heston paths we first need to examine how to generate correlated random variables in Python. Let's take a quick example first. Let \(\Sigma\) be the covariance matrix between two random normal variables:


\(\Sigma = \begin{pmatrix}1 & \rho\\\ \rho & 1\end{pmatrix}\)



 The means for each random variable are shown below as a vector.


\(\mu = \begin{pmatrix}\mu_1 \\ \mu_2 \end{pmatrix}\)



We can then draw from the following distribution ;


 \(W_t\sim\mathcal{N}(\mu, \Sigma)\)


Let's try that in Python and using a rho of -0.70. 


import numpy as np
import matplotlib.pyplot as plt

rho = -0.7
Ndraws = 1000
mu = np.array([0,0])
cov = np.array([[1, rho] , [rho , 1]])

W = np.random.multivariate_normal(mu, cov, size=Ndraws)

plt.title('Correlated Random Variables')



correlated random variables python


It appears from the plot above that the two processes do in fact move inversely as desired. We can now check statistically whether the desired relationship was created among the variables. 



#array([[ 1.        , -0.70388879],
#       [-0.70388879,  1.        ]])


Ok now that we are satisfied that it works we can introduce the concept of stochastic volatility and its relationship to asset returns. In a previous article on the Black-Scholes formula we discussed the problem we assuming that volatility is constant. The Heston model treats the volatility as a random variable which follows an Ornstein–Uhlenbeck process which is nice because it is able to replicate empirically observed properties of volatility. 



Properties of Volatility


Mean Reversion

Notice in the plot below for Apple that the volatility never jumps without eventually returning to a certain level, which looks to be approximately 25%.

Recall the notation we introduced at the beginning of the document \( \ \theta \) in this case would be approximately \(0.25^2 = 0.0625\)  the rate at which it is pulled towards this long-term variance is \(\kappa\) this behavior makes sense intuitively also since volatility can't go to zero or infinity.

 Notice the large spike in volatility in the first quarter of 2020 , it is followed by a return towards some long-term mean that we would need to estimate. 


instability of volatility example for AAPL



Clustering relates to the fact that if volatility is high today it is likely it will be high tomorrow also and vice versa for low volatility periods. In technical terms this means the volatility of a stock exhibits autocorrelation. 



Simulating a Heston Process


Essentially what we are doing here is just simulating a standard geometric Brownian motion with non-constant volatility, where the change in S has relationship \(\rho\) with the change in volatility. Here we set the drift term \(\mu\) from above to equal the risk free rate \(r\)


\(S_{t+i} = S_t e^{(r - \frac{v_t}{2})dt + \sqrt{v_i} dW_{t}^S }\)

\(v_{t+1} = v_t + \kappa(\theta-v_t)dt + \xi \sqrt{v_t} dW_t^v\)


So clearly we need to decide on an initial value of \(v_0\) which possibly could be estimated by using a 30 period window. The initial value for \(S_0\) is obviously just the price we observe currently in the market. 


Note that in that code we have assured that \(v\) is strictly positive as negative variance values aren't logical. Therefore where you see \(v \ \text{above think} \ |v|\)


def generate_heston_paths(S, T, r, kappa, theta, v_0, rho, xi, 
                          steps, Npaths, return_vol=False):
    dt = T/steps
    size = (Npaths, steps)
    prices = np.zeros(size)
    sigs = np.zeros(size)
    S_t = S
    v_t = v_0
    for t in range(steps):
        WT = np.random.multivariate_normal(np.array([0,0]), 
                                           cov = np.array([[1,rho],
                                           size=Npaths) * np.sqrt(dt) 
        S_t = S_t*(np.exp( (r- 0.5*v_t)*dt+ np.sqrt(v_t) *WT[:,0] ) ) 
        v_t = np.abs(v_t + kappa*(theta-v_t)*dt + xi*np.sqrt(v_t)*WT[:,1])
        prices[:, t] = S_t
        sigs[:, t] = v_t
    if return_vol:
        return prices, sigs
    return prices



We can now generate paths for a Heston process, the paths below were generated with an annual volatility of 20%. The Python snippet to plot paths for a Heston process is given at the end of the document to avoid clutter, see 'plotting paths of a Heston process' for the parameters used for this demonstration. 


Heston Model Simulation for 5 paths


We used a \(\kappa\) of 3 here, chosen for convenience. In practice these parameters would mostly likely be estimated from market prices. Perhaps we will create a separate article on how that is done.  


Heston Stochastic Volatility simulation



Let's take a look at different distributions of \(S_T\) when varying the correlation measure \(\rho\) between the Weiner processes. 


kappa =4
theta = 0.02
v_0 =  0.02
xi = 0.9
r = 0.02
S = 100
paths =50000
steps = 2000
T = 1

prices_pos = generate_heston_paths(S, T, r, kappa, theta,
                                    v_0, rho=0.9, xi=xi, steps=steps, Npaths=paths,
prices_neg  = generate_heston_paths(S, T, r, kappa, theta,
                                    v_0, rho=-0.9, xi=xi, steps=steps, Npaths=paths,
gbm_bench = S*np.exp( np.random.normal((r - v_0/2)*T , 
                                np.sqrt(theta)*np.sqrt(T), size=paths))

import seaborn as sns

fig, ax = plt.subplots()

ax = sns.kdeplot(data=prices_pos, label=r"$\rho = 0.9$", ax=ax)
ax = sns.kdeplot(data=prices_neg, label=r"$\rho= -0.9$ ", ax=ax)
ax = sns.kdeplot(data=gbm_bench, label="GBM", ax=ax)

ax.set_title(r'Tail Density by Varying $\rho$')
plt.axis([40, 180, 0, 0.055])



Heston Model skew in distribution by varying rho



Pricing Options with Heston Model


Let's take the terminal prices we got from the simulation above when \(\rho = 0.9\) and price options for a range of strikes. We will price a chain of puts between 30 - 200$. And investigate whether we get a volatility smile. To run the script below you will need the BS_CALL, BS_PUT and implied_vol function from this article on the subject



strikes =np.arange(30, 200,1)

puts = [] 

for K in strikes:
    P = np.mean(np.maximum(K-prices_neg,0))*np.exp(-r*T)

ivs = [implied_vol(P, S, K, T, r, type_ = 'put' ) for P, K in zip(puts,strikes)]

plt.plot(strikes, ivs)
plt.ylabel('Implied Volatility')
plt.axvline(S, color='black',linestyle='--',
            label='Spot Price')
plt.title('Implied Volatility Smile from Heston Model')


It appears we do in fact get a volatility smile as expected, albeit a rather weird looking one. This is probably due to the unrealistically high correlation we have used in this example. We will leave it to the reader to play around with these parameters and inspect the resulting smiles. 


Heston Model Volatility Smile


As mentioned previously Heston derived a closed form solution, which is necessary in order to calibrate the model to market prices as we have shown for the Merton jump diffusion model



So to summarize, clearly the Heston model is an improvement on the Black-Scholes model, due to the ability to model volatility as a mean-reverting random variable. 


It may be interesting to try to incorporate jumps to both the stock return and the volatility, perhaps trying to model stochastic correlation may be an interesting topic also. All of which we will endeavor to cover in future articles. 


Plotting Paths of Heston Process

kappa =3
theta = 0.04
v_0 =  0.04
xi = 0.6
r = 0.05
S = 100
paths =3
steps = 10000
T = 1
rho = -0.8
prices,sigs = generate_heston_paths(S, T, r, kappa, theta,
                                    v_0, rho, xi, steps, paths,
plt.title('Heston Price Paths Simulation')
plt.xlabel('Time Steps')
plt.ylabel('Stock Price')

plt.axhline(np.sqrt(theta), color='black', label=r'$\sqrt{\theta}$')
plt.title('Heston Stochastic Vol Simulation')
plt.xlabel('Time Steps')




Join the discussion

Share this post with your friends!