Build the Heston Model from scratch in Python — Part II: Calibration
In the previous section, we went over the intuition behind the Heston model.
If you already have options quotes and risk free rates from Bloomberg/Reuters/Excel, you can skip directly to the calibration section.
To get the .ipynb and .py files, you can find these on our Github!
Table of Contents
- Recap
- Why Calibrate?
- Calibration process: Intuition
- Web scrapping for quotes
- Web scrapping for risk free rates
- Define the objective function
- Calibrate
Recap
As a recap, we showed how to implement this formula:
where
Our resulting code was:
import numpy as np# Parallel computation using numba
from numba import jit, njit, prange,
from numba import cudai = complex(0,1)# To be used in the Heston pricer
@jit
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 (allow for parallel processing with numba)
@jit(forceobj=True)
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 prange(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 = P + ds *(numerator1 - numerator2)/denominator
element2 = P/np.pi
return np.real((element1 + element2))
I’ve added some numba elements because in calibration, we will need to speed up the computations. We use numba to optimize the evaluation of the functions.
Why Calibrate?
In practice, we actually don’t use the pricer to calculate prices. What we do is try to obtain the parameters for the Heston model based on the prices we observe in the markets.
These parameters are used in pricing complex/exotic options. For very complex options, we can’t really find a closed form formula for pricing exotic options.
What we do in this case is use the plain vanilla options to find parameters of the Heston model and use the parameters we obtain from the market to calculate the price of the exotic options using Monte Carlo.
The objective of this document is to show exactly how this is done.
Calibration process: Intuition
When calibrating you have to consider that you have options with different maturities and expiries.
You want to have one set of parameters for the entire set of options or at the very least, have one for the most important regions of the set of options.
If you are calculating sensitivities e.g. delta, vega, gamma. You want these to mean the same thing for all the options concerned. If you have a model for each point, you can’t aggregate sensitivities because it’s basically an apples to oranges comparison.
We start with a set of parameters, e.g.{𝜎=0.1 ,𝜅=0.1 ,𝜃=0.1 ,ρ=0.1 ,ν=0.1}. We calculate the Heston price for all our options and find the squared error between the market prices and the prces we obtained from our Heston.
We then use a solver to reduce this error to close 0. So this can actually be done on Excel. It’ll take a bit of time but it can be done.
Then using out optimal parameters we can price our options using the Monte Carlo Method.
Enough Talking, Let’s calibrate
Workflow
- Obtain options quotes and construct the surface (webscrapping or data API eg Bloomberg/Reuters)
- Obtain risk free rates (USSW for US or EONIA/EIOPA Curves for Europe)
- Create a vector function to calculate the prices of all the options on our surface
- Minimize the square error between all the prices calculated using our model & the market prices
First thing is to get option quotes. You can get these on Bloomberg if you have the terminal or scrape them from the web. We’ve done a section on this on this webscrapping repository on Github.
You can test it out to get option quotes from the internet.
I’ll do a quick sample here.
Web scrapping for option quotes
Options quotes typically don’t go out very long into the future. This is essentially because it would be very difficult to say where markets will be in the next 10 years for example.
On the other hand, with bonds, options on fixed income products do tend to go out very long into the future because of the maturities of bonds which can go as long as infinity (perpetuity bonds.)
NOTE: We’re only web scrapping because this is a private project. For work purposes, and with access to Bloomberg/Reuters, this should be a very quick set-up.
We’re using S&P 500 index options to calibrate our Heston. However, it’s possible to do this with any specific stock that has enough quoted options e.g. AAPL,MSFT.
Before importing packages in the python script, you will need to install selenium and beautifulsoup4 because they don’t come as default installed. In the command line, simply run:
pip install selenium
pip install beautifulsoup4Next we import our packages. BeautifulSoup allows us to make sense of the HTML pages we obtain from the internet.
Time and datetime are used to manage time and dates in python. We use time to pause Python for a few seconds as internet pages load.
Re is the default regular expression package in python. It allows us to search for, filter out and manipulate large bodies of strings.
Selenium is the key package used in web scrapping. We use Chromium to open a driver that will run all our automated commands.
from bs4 import BeautifulSoup
import re
import time
import pandas as pd
import os
from datetime import datetime as dt
import time
import itertools# Selenium Related packages
###############################
from selenium import webdriver
from selenium.webdriver.support.ui import Select
from selenium.webdriver.common.by import By
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
We have all the packages we need for getting our quotes.
The first thing we do is connect to the site we want to get quotes from. We’ll be using Barchart which provides a fairly large surface.
# We're using the barchart page to extract options quotes
url = 'https://www.barchart.com/stocks/quotes/$SPX/options'Next we open our Chrome driver. You can definitely use Firefox if you’re more comfortable with it. Next we connect to the url we defined above. We stop execution for 2 seconds to allow for the site to load. (You can definitely change or remove this depending on your internet speed)
# On ubuntu: sudo apt-get install chromium-chromedriver
# Connect to driver using Python
# Setting up options for WebDriver(Firefox, Safari, Edge)
driver = webdriver.Chrome()## Direct to Home Page
driver.get(url)time.sleep(2)
Small HTML detour
If you right-click on the Call quote table and select Inspect element you will obtain a window similar to this.
To get a specific element from all the mumbled up code, you can Inspect the element once again. This time, it should highlight the specific item you want to inspect.
For example, we are interested in the call option table. So we right click on the Strike text on top of the table, we will be directed to the following section:
At the very top, we’ll see a section called table. This is the part that interests us. We’re trying to extract details from this table.
We can also see a thead element which refers to the column heads.
Right after we see the tr element which represents a row. Similarly, we can seeth element which represents a column on the respective tr. For example, we can see that the th in view is given the class attribute as strikePrice.
So when we look for an element by xpath essentially what we’re doing is trying to find elements with the tags or names or classes defined in a certain way on our HTML text.
For example, adiv element in HTML is just a container with other text or buttons inside it. So we are essentially looking for a div element with the class attribute equal to bc-table-scrollable-inner. You should get both the call and put tables because they are the only scrollable tables on the page.
Now back to selenium
We accept terms and conditions:
# Accept terms and conditions
driver.find_elements_by_xpath("//*[contains(text(),
'Accept all')]")[0].click()We can start obtaining quotes. We filter out quotes that are not liquid by selecting only options that are close to the money. To do this, we search for a select element that contains a name tag called moneyness.
# We want to obtain quotes close to the money
rows = Select(driver.find_element_by_xpath(
"//select[@name = 'moneyness']"))
rows.select_by_value('inTheMoney')For different maturities, we will have different quotes for different strikes. On the website we obtain all the possible maturities provided in the drop-down with dates. Similarly here, we look for a select element containing a tag data-event-name called onExpirationDateChanged.
# Select a date on the options market
mButtons =Select(driver.find_element_by_xpath(
"//select[@data-event-name = 'onExpirationDateChanged']"))surface = []
dateIndex = []
We have an empty list surface that will contain all our quotes as a list of lists. Each component list will essentially represent a maturity date. The lists are the same length because we’re focusing on the same strikes.
The dateIndex will contain the maturity dates. We use these to calculate the maturity dates with respect to today in years.
So for each buttonin mButtons, we want to obtain the strike and the corresponding option price. We loop through each button element to obtain our quotes.
We add a try-except expression to ignore quotes we can’t find possibly due to different formats. This is just a quick and easy way to get quotes because either way, we will have more than we need.
for option in quotes.options:
# Ignore quotes we can't find
try:
option.click() # Wait for page to load
time.sleep(5) # Find the call and put tables
#Select first table (Call)
tableDiv = option.find_elements_by_xpath(
"//div[@class = 'bc-table-scrollable-inner']")[0] # We want the first table(Call quotes)
callQuotes = tableDiv.find_elements_by_xpath("//tbody")[0]
callQuotes = callQuotes.find_elements_by_xpath(
"//tr[@sly-repeat = 'row in rows.data']") # We want to obtain all the quotes from each row
regex = '[0-9,.]{2,}'
quotes = [re.findall(regex, elem.text)
for elem in callQuotes] #Strikes are element 0, mid quotes element 4 on table
strikes = [float(re.sub(',','', quote[0]))
for quote in quotes]
midQuotes = [float(re.sub(',','', quote[3]))
for quote in quotes] # Construct the dataframe
date = re.search('[0-9]{4}-[0-9]{2}-[0-9]{2}',
option.text).group()
date = dt.strptime(date, '%Y-%m-%d') # Create the date
surface.append(midQuotes)
dateIndex.append(date) except Exception:
pass
We break down each element of this loop. First, by clicking the option we are essentially selecting a date on the drop-down. We then wait a few seconds for the page to load.
option.click()# Wait for page to load
time.sleep(5)
We want to find the tables of quotes. In our case, we will only have the call and put tables. We select the first table as it appears on the page.
# Find the call and put tables
#Select first table (Call)
tableDiv = option.find_elements_by_xpath(
"//div[@class = 'bc-table-scrollable-inner']")[0]Similarly, we find items with the tag tbody which should essentially mean table body. Next, we find the tr which means a table row with the specific tag or name we used.
# We want the first table(Call quotes)
callQuotes = tableDiv.find_elements_by_xpath("//tbody")[0]
callQuotes = callQuotes.find_elements_by_xpath(
"//tr[@sly-repeat = 'row in rows.data']")Next, we want to obtain the quotes themselves. We find all quotes that are numeric containing a full stop and longer than 2 digits.
The numbers have commas so we remove them with the re.sub() function.
We add the quotes to the surface list. Similarly, we add the maturity to thedateIndex list.
# We want to obtain all the quotes from each row
regex = '[0-9,.]{2,}'
quotes = [re.findall(regex, elem.text) for elem in callQuotes]# Strikes are the first elements, mid quotes are the 5th element
strikes = [float(re.sub(',','', quote[0])) for quote in quotes]
midQuotes = [float(re.sub(',','', quote[3])) for quote in quotes]# Construct the dataframe
date = re.search('[0-9]{4}-[0-9]{2}-[0-9]{2}', option.text).group()
date = dt.strptime(date, '%Y-%m-%d')# Create the date
surface.append(midQuotes)
dateIndex.append(date)
Finally, we construct a volatility surface:
# Construct a full volatility surface
maturities = [(date - dt.today()).days/365.25
for date in dateIndex]volSurface = pd.DataFrame(surface,
index = maturities,
columns = strikes)volSurface = volSurface[volSurface.index > 0.1]
volSurface
The result is a fairly large surface:
Web scrapping for risk free rates
We use the same structure to find the risk free curve for the previous day’s date.
To do this, we use the treasury rates from this website. Again, this would be easier with some adjustments to the swap rates(EUSA or USSW in Bloomberg.)
We then apply an interpolation method to obtain our interest rate curve. This will be important to find the rate corresponding to each maturity.
To do this interpolation, we’ll use the Nelson Siegel Svensson method. This is one of the most common methods of interpolating a curve. For actuarial purposes in Europe, however, the Smith Wilson interpolation is used.
I’ll implement the Smith Wilson as it’s used in insurance at the end of the document.
ratesUrl = 'https://www.treasury.gov/resource-center/data-chart-center/interest-rates/Pages/TextView.aspx?data=yield'## Direct to Home Page
driver.get(ratesUrl)time.sleep(2)# Find the table of rates
rows = driver.find_elements_by_xpath("//tr")# Extract the last line (the latest possible rate)
rates = [elem.text for elem in rows[len(rows)-1].find_elements_by_xpath("//td[@class = 'text_view_data']")]# End the chromium driver
driver.close()# Convert the rates we obtained to float.
# Convert to a numpy array - this is the input format for the NSS
rates = np.array(rates[-12:]).astype(float)/100
maturities = np.array([1/12, 2/12, 3/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30])# Apply the nelson siegel svennson model
curve, status = calibrate_ns_ols(maturities, rates)
We can now view the extrapolated curve:
Define the Objective Function
We begin by initializing our parameters. We use these initial parameters to test our whether our heston pricer is coherent.
We convert our surface to a list of quotes containing the maturity, strike and quote for each row.
#Initialize parameters
sigma, kappa, theta, volvol, rho = 0.1, 0.1, 0.1, 0.1, 0.1# Convert our vol surface to long
volSurfaceLong = volSurface.melt(ignore_index =
False).reset_index()
volSurfaceLong.columns = ['maturity', 'strike', 'price']# Calculate the risk free rate for each maturity
volSurfaceLong['rate'] = volSurfaceLong['maturity'].apply(curve)# Define global variables to be used in optimization
maturities = volSurfaceLong['maturity'].to_numpy('float')
strikes = volSurfaceLong['strike'].to_numpy('float')
marketPrices = volSurfaceLong['price'].to_numpy('float')
rates = volSurfaceLong['rate'].to_numpy('float')
Calibrate
Another simple function OFHest should give us a vector of differences between the prices we obtained in our model and those quoted in the market.
When we fit our optimization algorithm, these errors will be squared and summed on each iteration. This is how the algorithm is able to continue minimizing the error until an optimal solution is found.
# This is the calibration function
def calibratorHeston(St, initialValues = [0.5,0.5,0.5,0.5,-0.5],
lowerBounds = [1e-2,1e-2,1e-2,1e-2,-1],
upperBounds = [10,10,10,10,0]): '''1) Define parameters
=========================='''
params = Parameters()
params.add('sigma',value = initialValues[0],
min = lowerBounds[0], max = upperBounds[0])
params.add('kappa',value = initialValues[1],
min = lowerBounds[1], max = upperBounds[1])
params.add('theta',value = initialValues[2],
min = lowerBounds[2], max = upperBounds[2])
params.add('volvol', value = initialValues[3],
min = lowerBounds[3], max = upperBounds[3])
params.add('rho', value = initialValues[4],
min = lowerBounds[4], max = upperBounds[4])
'''2) Define objective function
================================'''
OFHest = lambda param:(marketPrices-
priceHestonMid(St,
strikes,
rates,
maturities,
paramVect['sigma'].value,
paramVect['kappa'].value,
paramVect['theta'].value,
paramVect['volvol'].value,
paramVect['rho'].value))/
marketPrices'''3) Optimize parameters
=============================='''
result = minimize(OFHest,
params,
method = 'leastsq',
iter_cb = iter_cb,
ftol = 1e-6)
return(result)
The results should be:
Our optimal parameters should be:
Finally, we get correlations:
And the error color graph:
Annex: Smith Wilson
#Simple Smith Wilson Implementation'
'----------------------------------------'
class SmithWilson:
'''Implementation based on: https://www.scor.com/sites/default/files/2017_it_tesidi_laurea_magurean.pdf
Inputs
---------
ZCPrices = Zero coupon curve ie 0.99
maturities = maturities for each of the ZC prices
alpha = default parameter in Smith Wilson.
Determines the quickness of convergence
UFR = the theoretical ultimate forward rate
(The forward rate we suppose rates converge to)
'''#These are the default input maturities
maturities = np.arange(1, 41)
def __init__(self, ZCPrices,
maturities = maturities,
alpha = 0.134, UFR = 0.039):
self.UFR = np.log(1+UFR)
self.alpha = alpha
self.maturities = maturities
self.ZCPrices = ZCPrices
'Define the Wilson Function'
def wilsonFunction(self, t, uj):
'''
Exact same function on wikipedia: https://en.wikipedia.org/wiki/Smith%E2%80%93Wilson_method
Inputs
-------------
uj = time to maturity
t = maturity. This is a specific maturity for a ZC Bond
'''
mu = np.exp(- self.UFR*(t+uj))
maxT = self.alpha * np.maximum(t, uj)
minT = self.alpha * np.minimum(t, uj)
# We define the Wilson kernel function
wj = mu*(minT-0.5*np.exp(-maxT)*
(np.exp(minT)-np.exp(-minT)))
return (wj)
def wilsonVectorFunction(self, inputMaturities):
'''
Calculates a square matrix,
calculating the Wilson Function
for each maturity on each maturity
Inputs
---------------
inputMaturities = vector of maturities we input
'''
wilsonMatrix = np.zeros(len(inputMaturities),
len(self.maturities)))
for t in inputMaturities:
for j in self.maturities:
wilsonMatrix[t-1, j-1] = self.wilsonFunction(t, j)
return(wilsonMatrix)
'Obtain the parameter vector zeta'
def calibrate(self):
#Create the matrix W of Kernels
self.WMatrix = self.wilsonVectorFunction(self.maturities)
#Invert the matrix
#Recall that our parameters are W^-1 * (ZCPrices - muVector)
# Invert the kernel
WMatrixInv = np.linalg.inv(self.WMatrix)
muVector = np.exp(-self.UFR * self.maturities)
SWParams = WMatrixInv.dot(self.ZCPrices - muVector)
return (SWParams)
'Fit Curve'
def curve(self):
params = self.calibrate() #Obtain parameters
parametrizedWilson = self.wilsonVectorFunction(np.arange(1,
151)).dot(params)
result = np.exp(-self.UFR * np.arange(1, 151))+
parametrizedWilson
return(result)
