# Option Prices And Calibrating The Heston Model Code Question

I'm trying to understand this Python code that uses Quantlib to calibrate the parameters of the Heston model. The data that is provided in the code is the spot price, the risk free interest rate, the dividends of the underlining, the day of calibration, expiration dates of options, their strike prices, and their implied volatilities.

What isn't given is the the market prices, yet the code generates a comparison between the market prices and the price computed by the Heston model. How are those prices generated? How can the parameters be fit if the market prices of the options aren't explicitly given? Isn't the point to choose the parameters so that difference between predicted market prices via the Heston model and the actual market prices are minimized?

This is the original source of the code.

Any help clarifying this will be much appreciated.

https://github.com/AIMLModeling/Heston-Model-Calibration/blob/main/HestonVolVisuals.py

import QuantLib as ql
import math
import numpy as np from mpl_toolkits.mplot3d
import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
day_count = ql.Actual365Fixed()
calendar = ql.UnitedStates()

calculation_date = ql.Date(9, 11, 2021)
spot = 659.37
ql.Settings.instance().evaluationDate = calculation_date

dividend_yield = ql.QuoteHandle(ql.SimpleQuote(0.0))
risk_free_rate = 0.01
dividend_rate = 0.00
flat_ts = ql.YieldTermStructureHandle(
ql.FlatForward(calculation_date, risk_free_rate, day_count))
dividend_ts = ql.YieldTermStructureHandle(
ql.FlatForward(calculation_date, dividend_rate, day_count))

expiration_dates = [ql.Date(9,12,2021), ql.Date(9,1,2022), ql.Date(9,2,2022),
ql.Date(9,3,2022), ql.Date(9,4,2022), ql.Date(9,5,2022),
ql.Date(9,6,2022), ql.Date(9,7,2022), ql.Date(9,8,2022),
ql.Date(9,9,2022), ql.Date(9,10,2022), ql.Date(9,11,2022),
ql.Date(9,12,2022), ql.Date(9,1,2023), ql.Date(9,2,2023),
ql.Date(9,3,2023), ql.Date(9,4,2023), ql.Date(9,5,2023),
ql.Date(9,6,2023), ql.Date(9,7,2023), ql.Date(9,8,2023),
ql.Date(9,9,2023), ql.Date(9,10,2023), ql.Date(9,11,2023)]
strikes = [527.50, 560.46, 593.43, 626.40, 659.37, 692.34, 725.31, 758.28]
data = [
[0.37819, 0.34177, 0.30394, 0.27832, 0.26453, 0.25916, 0.25941, 0.26127],
[0.3445, 0.31769, 0.2933, 0.27614, 0.26575, 0.25729, 0.25228, 0.25202],
[0.37419, 0.35372, 0.33729, 0.32492, 0.31601, 0.30883, 0.30036, 0.29568],
[0.37498, 0.35847, 0.34475, 0.33399, 0.32715, 0.31943, 0.31098, 0.30506],
[0.35941, 0.34516, 0.33296, 0.32275, 0.31867, 0.30969, 0.30239, 0.29631],
[0.35521, 0.34242, 0.33154, 0.3219, 0.31948, 0.31096, 0.30424, 0.2984],
[0.35442, 0.34267, 0.33288, 0.32374, 0.32245, 0.31474, 0.30838, 0.30283],
[0.35384, 0.34286, 0.33386, 0.32507, 0.3246, 0.31745, 0.31135, 0.306],
[0.35338, 0.343, 0.33464, 0.32614, 0.3263, 0.31961, 0.31371, 0.30852],
[0.35301, 0.34312, 0.33526, 0.32698, 0.32766, 0.32132, 0.31558, 0.31052],
[0.35272, 0.34322, 0.33574, 0.32765, 0.32873, 0.32267, 0.31705, 0.31209],
[0.35246, 0.3433, 0.33617, 0.32822, 0.32965, 0.32383, 0.31831, 0.31344],
[0.35226, 0.34336, 0.33651, 0.32869, 0.3304, 0.32477, 0.31934, 0.31453],
[0.35207, 0.34342, 0.33681, 0.32911, 0.33106, 0.32561, 0.32025, 0.3155],
[0.35171, 0.34327, 0.33679, 0.32931, 0.3319, 0.32665, 0.32139, 0.31675],
[0.35128, 0.343, 0.33658, 0.32937, 0.33276, 0.32769, 0.32255, 0.31802],
[0.35086, 0.34274, 0.33637, 0.32943, 0.3336, 0.32872, 0.32368, 0.31927],
[0.35049, 0.34252, 0.33618, 0.32948, 0.33432, 0.32959, 0.32465, 0.32034],
[0.35016, 0.34231, 0.33602, 0.32953, 0.33498, 0.3304, 0.32554, 0.32132],
[0.34986, 0.34213, 0.33587, 0.32957, 0.33556, 0.3311, 0.32631, 0.32217],
[0.34959, 0.34196, 0.33573, 0.32961, 0.3361, 0.33176, 0.32704, 0.32296],
[0.34934, 0.34181, 0.33561, 0.32964, 0.33658, 0.33235, 0.32769, 0.32368],
[0.34912, 0.34167, 0.3355, 0.32967, 0.33701, 0.33288, 0.32827, 0.32432],
[0.34891, 0.34154, 0.33539, 0.3297, 0.33742, 0.33337, 0.32881, 0.32492]]
implied_vols = ql.Matrix(len(strikes), len(expiration_dates))
for i in range(implied_vols.rows()):
for j in range(implied_vols.columns()):
implied_vols[i][j] = data[j][i]
black_var_surface = ql.BlackVarianceSurface(
calculation_date, calendar,
expiration_dates, strikes,
implied_vols, day_count)

strikes_grid = np.arange(strikes[0], strikes[-1],10)
expiry = 1.0 # years
implied_vols = [black_var_surface.blackVol(expiry, s)
for s in strikes_grid] # can interpolate here
actual_data = data[11] # cherry picked the data for given expiry (1 year)

fig, ax = plt.subplots()
ax.plot(strikes_grid, implied_vols, label="Black Surface")
ax.plot(strikes, actual_data, "o", label="Actual")
ax.set_xlabel("Strikes", size=12)
ax.set_ylabel("Vols", size=12)
legend = ax.legend(loc="upper right")
plot_years = np.arange(0, 2, 0.1)
plot_strikes = np.arange(535.0, 750.0, 1.0)
fig = plt.figure()
ax = fig.subplots(subplot_kw={'projection': '3d'})
X, Y = np.meshgrid(plot_strikes, plot_years)
X, Y = np.meshgrid(plot_strikes, plot_years)

Z = np.array([black_var_surface.blackVol(y, x)
for xr, yr in zip(X, Y)
for x, y in zip(xr,yr) ]
).reshape(len(X), len(X[0]))

surf = ax.plot_surface(X,Y,Z, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0.1)
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()


# dummy parameters

v0 = 0.01; kappa = 0.2; theta = 0.02; rho = -0.75; sigma = 0.5;

process = ql.HestonProcess(flat_ts, dividend_ts,
ql.QuoteHandle(ql.SimpleQuote(spot)),
v0, kappa, theta, sigma, rho)
model = ql.HestonModel(process)
engine = ql.AnalyticHestonEngine(model)
heston_helpers = []
black_var_surface.setInterpolation("bicubic")
one_year_idx = 11 # 12th row in data is for 1 year expiry
date = expiration_dates[one_year_idx]
for j, s in enumerate(strikes):
t = (date - calculation_date )
p = ql.Period(t, ql.Days)
sigma = data[one_year_idx][j]
#sigma = black_var_surface.blackVol(t/365.25, s)
helper = ql.HestonModelHelper(p, calendar, spot, s,
ql.QuoteHandle(ql.SimpleQuote(sigma)),
flat_ts,
dividend_ts)
helper.setPricingEngine(engine)
heston_helpers.append(helper)
lm = ql.LevenbergMarquardt(1e-8, 1e-8, 1e-8)
model.calibrate(heston_helpers, lm,
ql.EndCriteria(500, 50, 1.0e-8,1.0e-8, 1.0e-8))
theta, kappa, sigma, rho, v0 = model.params()

print ("\ntheta = %f, kappa = %f, sigma = %f, rho = %f, v0 = %f" % (theta, kappa, sigma, rho, v0))
avg = 0.0

print ("%15s %15s %15s %20s" % (
"Strikes", "Market Value",
"Model Value", "Relative Error (%)"))
print ("="*70)
for i, opt in enumerate(heston_helpers):
err = (opt.modelValue()/opt.marketValue() - 1.0)
print ("%15.2f %14.5f %15.5f %20.7f " % (
strikes[i], opt.marketValue(),
opt.modelValue(),
100.0*(opt.modelValue()/opt.marketValue() - 1.0)))
avg += abs(err)
avg = avg*100.0/len(heston_helpers)
print ("-"*70)
print ("Average Abs Error (%%) : %5.3f" % (avg))

theta = 0.131545, kappa = 10.957884, sigma = 3.991680, rho = -0.351518, v0 =
0.076321
Strikes    Market Value     Model Value   Relative Error (%)
======================================================================
527.50       44.53343        44.32072           -0.4776460
560.46       54.89815        55.07795            0.3275233
593.43       67.20964        67.50343            0.4371238
626.40       80.76865        81.66110            1.1049427
659.37       98.71891        97.58279           -1.1508670
692.34       93.09155        92.34694           -0.7998741
725.31       79.44659        79.24588           -0.2526358
758.28       67.42434        67.77431            0.5190524


Average Abs Error (%) : 0.634