I am trying to fit a time dependent Heston model using Quantlib Python. I'm getting the following runtime error: Boost assertion failed : px !=0.
Can somebody help in this or is there an example of time dependent heston ?
import numpy as np
#% matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.optimize import root
from scipy.optimize import root, least_squares, differential_evolution, basinhopping
from pandas import DataFrame
day_count = ql.Actual365Fixed()
calendar = ql.UnitedStates()
calculation_date = ql.Date(7, 6, 2023)
spot = 4271.21
ql.Settings.instance().evaluationDate = calculation_date
dividend_yield = ql.QuoteHandle(ql.SimpleQuote(0.0))
risk_free_rate = 0.0553
dividend_rate = 0.01283
flat_ts = ql.YieldTermStructureHandle(
ql.FlatForward(calculation_date, risk_free_rate, day_count))
yield_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(7,7,2023), ql.Date(7,9,2023),
ql.Date(7,12,2023), ql.Date(7,6,2024)]
ttm=[ql.Period(m-calculation_date, ql.Days) for m in expiration_dates]
strikes = [3844.09, 4057.65, 4164.43, 4271.21, 4377.99, 4484.77, 4698.33]
data = [
[ 0.22004, 0.16338, 0.13716, 0.11289 ,0.10017, 0.10047, 0.13063 ],
#[, 0.20691, 0.16693, 0.14635, 0.12616 ,0.11087, 0.10298, 0.10272 ,0.17100],
[ 0.20358, 0.16993, 0.15206, 0.13394 ,0.11860, 0.10861, 0.10283 ],
[ 0.20554, 0.18071, 0.16752, 0.15401 ,0.14106, 0.12976, 0.11482 ],
#[0.24909, 0.20919, 0.18866, 0.17759, 0.16614 ,0.15474, 0.14404, 0.12691 ,0.11426],
[ 0.21129, 0.19361, 0.18420, 0.17443 ,0.16453, 0.15487, 0.13788 ]]
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)
strike = 4300
expiry = 0.6 # years
black_var_surface.blackVol(expiry, strike)
v0 = 0.02; kappa = 3; theta = 0.02; rho = -0.8; sigma = 2;
times = [0.0833,0.25,0.5,1]
for i, time in enumerate(times):
kappaTS.setParam(i, kappa[i])
#thetaTS.setParam(i, theta[i])
#rhoTS.setParam(i, rho[i])
#sigmaTS.setParam(i, sigma[i])
grid = ql.TimeGrid(times)
# process = ql.HestonProcess(flat_ts, dividend_ts,
# ql.QuoteHandle(ql.SimpleQuote(spot)),
# v0, kappa, theta, sigma, rho,grid)
model = ql.PiecewiseTimeDependentHestonModel(flat_ts, dividend_ts,
ql.QuoteHandle(ql.SimpleQuote(spot)),
v0, ql.Parameter(),ql.Parameter(),
ql.Parameter(), ql.Parameter(),grid)
engine = ql.AnalyticPTDHestonEngine(model)
# engine = ql.FdHestonVanillaEngine(model)
heston_helpers = []
black_var_surface.setInterpolation("bicubic")
for i, date in enumerate(expiration_dates):
for j, s in enumerate(strikes):
t = (date - calculation_date )
p = ql.Period(t, ql.Days)
sigma = data[i][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()
model.calibrate(heston_helpers, lm, ql.EndCriteria(500, 300, 1.0e-8, 1.0e-8, 1.0e-8))
