We want to try and get the intuition behind the model so that we can implement and use it. The derivation is less important to us in this document.
We also provide the base pricing formula and the different problems that could occur while pricing with the Heston model.
You only need to know about probability distributions and integrals. The rest is just adjustments to what you already know.
This document is also on GitHub. So if you want to get the .ipynb or .py files you can feel free to do so here!
Table of Contents
Black Scholes Model
We begin by introducing the Black and Scholes pricer and make parallels to the base options pricing model in comparison to the Heston.
We know that the Black and Scholes pricing model gives the following final formula for a call option:
𝑁 is the cumulative distribution function for a normal random variable. 𝑁() is just the NORM.DIST function in Excel. This function takes the following inputs:NORM.DIST(x, mean, standard_dev)and calculates, for any input x, the area under the bell curve starting from the left.
Say we have 0 as the mean and 1 as the standard deviation.
This means that at 0 the area of the graph below starting from the furthest point on the left will be exactly half.
This area represents the probability of a random value selected from a distribution of mean 0 and standard deviation 1 being less than or equal to 𝑥.
Therefore, the larger the number 𝑥x the higher the probability that any random number selected will be less than 𝑥.
From Black Scholes to the Heston Model
First have to get a bird’s eye view before we venture into the details. The formula seems daunting but in a few steps, we’ll show you that it’s not really that far from the Black & Scholes model in practice.
Our example is just a simple call option. Lets say you have a European option expiring in 2 years with a strike at 100 and the price of the stock (say Apple) is at 110 today. You’d expect the value of the option to be about 10. We also just assume that the risk free interest rate is 5%.
We can use some general mathematical notation and define:
- 𝐶 as the value of the call option
- 𝑆𝑡 the stock price today
- 𝐾 the strike price
- 𝑇 the maturity (therefore 𝑇−𝑡 = τ is the time to maturity)
- 𝑟 the risk free rate
Our Heston pricing formula for a call option is:
In the Black Scholes Model P1 and P2 functions are actually the N(d1) and N(d2).
In the Heston, we replace N(d1) with:
and N(d2) with:
Don’t be intimidated by the formula. Remember that the use is just the same as N(d1) and N(d2).
So our full formula is:
Heston Model Breakdown: component by component
We will break this down piece by piece so we have everything we need. After that, we will start with the simplest elements and build up to the full integral.
We have two integrals. There is an Re encapsulating the entire integral. This just means that we have some complex numbers and when we calculate the price, we’ll only consider the real number part. So don’t worry if you haven’t dealt with complex numbers because you won’t have to!
Also, we use i to define the complex number i = √-1
So we go back to our P1 or P2. We can shorten these to Pj
Each integral has a function: fj. These functions are defined as:
and the functions C and D are defined as:
We also define the final elements x, d, g and u that are used in fj:
We also have a BRS element that is defined as:
Simple Heston Version
Instead of doing two integrals, you can put everything under one integral and one function f. (Just to be sure, I’ve shown how to do this in the Annex of this jupyter notebook document)
and so we define one function f:
Code: Step by step build
- Defining f (fHeston)
So we start by breaking down this monster of a function. We’ll start by defining 𝑓 as fHeston with inputs s, St, K, r, T, sigma, kappa, theta, volvol, rho representing our inputs and the model parameters {𝜎,𝜅,𝜃,𝜈,ρ} respectively.
Because our function will be using complex numbers, we can just define i as a global variable. Every complex number will essentially just be an operation on i.
You’ll also note that in the implementation, you can just treat i like any other variable and you should be fine. So no worries if you're not sure how to work with complex numbers!
The expression ρ.𝜎.𝑖.𝑠 is used quite a lot. We can define this beforehand and just call it prod
import numpy as npi = complex(0,1)# To be used in the Heston pricer
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
# To be used a lot
prod = rho * sigma *i *s
We extend our function. We split and calculate d and g
# To be used in the Heston pricer
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
# To be used a lot
prod = rho * sigma *i *s
# Calculate d
d1 = (prod - kappa)**2
d2 = (sigma**2) * (i*s + s**2)
d = np.sqrt(d1 + d2)
# Calculate g
g1 = kappa - prod - d
g2 = kappa - prod + d
g = g1/g2Next break down the large exponential. We deal with the first part of the exponential:
# To be used in the Heston pricer
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
# To be used a lot
prod = rho * sigma *i *s
# Calculate d
d1 = (prod - kappa)**2
d2 = (sigma**2) * (i*s + s**2)
d = np.sqrt(d1 + d2)
# Calculate g
g1 = kappa - prod - d
g2 = kappa - prod + d
g = g1/g2
# Calculate first exponential
exp1 = np.exp(np.log(St) * i *s) * np.exp(i * s* r* T)
exp2 = 1 - g * np.exp(-d *T)
exp3 = 1- g
mainExp1 = exp1*np.power(exp2/exp3, -2*theta*kappa/(sigma **2))We then deal with the second integral to obtain the full fHeston function
import numpy as npi = complex(0,1)
u = 1# To be used in the Heston pricer
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
# To be used a lot
prod = rho * sigma *i *s
# Calculate d
d1 = (prod - kappa)**2
d2 = (sigma**2) * (i*s + s**2)
d = np.sqrt(d1 + d2)
# Calculate g
g1 = kappa - prod - d
g2 = kappa - prod + d
g = g1/g2
# Calculate first exponential
exp1 = np.exp(np.log(St) * i *s) * np.exp(i * s* r* T)
exp2 = 1 - g * np.exp(-d *T)
exp3 = 1- g
mainExp1 = exp1*np.power(exp2/exp3, -2*theta*kappa/(sigma**2))
# Calculate second exponential
exp4 = theta * kappa * T/(sigma **2)
exp5 = volvol/(sigma **2)
exp6 = (1 - np.exp(-d * T))/(1 - g * np.exp(-d * T))
mainExp2 = np.exp((exp4 * g1) + (exp5 *g1 * exp6))
return (mainExp1 * mainExp2)
2. Calculate the integral
Setting the limit
Because our integral goes up to infinity, we cannot have something running continuously.
So what we do, is select a large limit that allows our integral to converge. By convergence, we simply mean that at some point, the additional areas under the curve we’re calculating are so small that they’re insignificant.
In our case, we selected 100. You can try this with different values to see how it works.
Selecting an integration scheme
To compute our integral, we’ll use a very simple integration scheme. This is the scheme we all used in elementary school. We simply split our area into rectangles and find the mid point. We use this mid point to calculate the area of each rectangle and sum them together.
In our case, we initialize P, the final price, as 0, we split our area into 1000 rectangles and fix the limit to 100.
So in essence, the width of each rectangle will be du = 100/1000.
We also calculate the first part before the integral:
# Heston Pricer
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
P, iterations, maxNumber = 0,1000,100
ds = maxNumber/iterations
element1 = 0.5 * (St - K * np.exp(-r * T))We define s1 and s2 as defined in the formula.
You observe that we increment s1 at each run of the loop. For example, at j=1, we make s1 = du * (2*1 +1)/2 which i just s1 = du * 1.5.
When j=2, then s1 = du * (2*2 +1)/2 which is s1 = 2.5 * du. This verifies our mid point rule for the integration.
You’ll also note that we don’t start from 0. This is because going close to 0 may mean dividing a number by 0. Of course that'll be a problem.
# Heston Pricer
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
P, iterations, maxNumber = 0,1000,100
ds = maxNumber/iterations
element1 = 0.5 * (St - K * np.exp(-r * T))
# Calculate the complex integral
# Using j instead of i to avoid confusion
for j in range(1, iterations):
s1 = ds * (2*j + 1)/2
s2 = s1 - i
numerator1 = fHeston(s2, St, K, r, T,
sigma, kappa, theta, volvol, rho)
numerator2 = K * fHeston(s1, St, K, r, T,
sigma, kappa, theta, volvol, rho)
denominator = np.exp(np.log(K) * i * s1) *i *s1So at this point, we have the width of all our rectangles, and the height is simply the value from the fHeston as we move along our x axis.
The sum of the results from fHeston divided by the denominator are our height. So we can now just calculate our areas and sum them together.
# Heston Pricer
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
P, iterations, maxNumber = 0,1000,100
du = maxNumber/iterations
element1 = 0.5 * (St - K * np.exp(-r * T))
# Calculate the complex integral
# Using j instead of i to avoid confusion
for j in range(1, iterations):
s1 = du * (2*j + 1)/2
s2 = u1 - i
numerator1 = fHeston(s2, St, K, r, T,
sigma, kappa, theta, volvol, rho)
numerator2 = K * fHeston(s1, St, K, r, T,
sigma, kappa, theta, volvol, rho)
denominator = np.exp(np.log(K) * i * s1) *i *s1
P += ds *(numerator1 - numerator2)/denominator
element2 = P/np.pi
return np.real((element1 + element2))Full Code
So we’ve managed to implement the Heston model in less than 50 lines of code!
import numpy as npi = complex(0, 1)# To be used in the Heston pricer
def fHeston(s, St, K, r, T, sigma, kappa, theta, volvol, rho):
# To be used a lot
prod = rho * sigma * i * s # Calculate d
d1 = (prod - kappa)**2
d2 = (sigma**2) * (i * s + s**2)
d = np.sqrt(d1 + d2) # Calculate g
g1 = kappa - prod - d
g2 = kappa - prod + d
g = g1 / g2 # Calculate first exponential
exp1 = np.exp(np.log(St) * i * s) * np.exp(i * s * r * T)
exp2 = 1 - g * np.exp(-d * T)
exp3 = 1 - g
mainExp1 = exp1 * np.power(exp2 / exp3, -2 * theta * kappa/ (sigma**2)) # Calculate second exponential
exp4 = theta * kappa * T / (sigma**2)
exp5 = volvol / (sigma**2)
exp6 = (1 - np.exp(-d * T)) / (1 - g * np.exp(-d * T))
mainExp2 = np.exp((exp4 * g1) + (exp5 * g1 * exp6)) return (mainExp1 * mainExp2)# Heston Pricer
def priceHestonMid(St, K, r, T, sigma, kappa, theta, volvol, rho):
P, iterations, maxNumber = 0, 1000, 100
ds = maxNumber / iterations element1 = 0.5 * (St - K * np.exp(-r * T)) # Calculate the complex integral
# Using j instead of i to avoid confusion
for j in range(1, iterations):
s1 = ds * (2 * j + 1) / 2
s2 = s1 - i numerator1 = fHeston(s2, St, K, r, T,
sigma, kappa, theta, volvol, rho)
numerator2 = K * fHeston(s1, St, K, r, T,
sigma, kappa, theta, volvol, rho)
denominator = np.exp(np.log(K) * i * s1) * i * s1 P += ds * (numerator1 - numerator2) / denominator element2 = P / np.pi return np.real((element1 + element2))
