我正在尝试使用库 quantlib 和全局优化器作为 DE 来校准 heston 模型的参数。
首先，我们有以下隐含波动率、执行价格和到期日列表：
data = [[0.22, 0.22, 0.19, 0.14, 0.14],[0.16, 0.19, 0.2, 0.15, 0.14, 0.14, 0.14, 0.11, 0.16, 0.1, 0.23, 0.1, 0.17, 0.08], [0.16, 0.22, 0.21, 0.2, 0.19, 0.18, 0.14, 0.11, 0.17, 0.08, 0.11, 0.11, 0.13, 0.16, 0.11, 0.12, 0.18, 0.12, 0.12, 0.14, 0.16, 0.15, 0.15, 0.16, 0.17, 0.14, 0.17, 0.15, 0.16, 0.16, 0.16, 0.16, 0.18, 0.17, 0.14, 0.16, 0.17, 0.16, 0.16, 0.18, 0.16], [0.2, 0.14], [0.16, 0.15, 0.18, 0.19, 0.16], [0.25, 0.2], [0.14], [0.21, 0.23, 0.17, 0.18, 0.11, 0.11, 0.11], [0.14], [0.16, 0.14, 0.14], [0.23, 0.25], [0.2, 0.19, 0.14, 0.12, 0.18, 0.17, 0.1, 0.12, 0.16, 0.16, 0.17, 0.16, 0.13, 0.16, 0.19], [0.12, 0.19, 0.13, 0.17], [0.1, 0.14, 0.13, 0.18], [0.2], [0.18, 0.11, 0.11, 0.1, 0.12, 0.12, 0.13, 0.17, 0.14], [0.16], [0.16, 0.11, 0.17, 0.12, 0.16, 0.16, 0.17], [0.18, 0.19, 0.15, 0.14, 0.14, 0.14, 0.18, 0.13], [0.23, 0.18, 0.16, 0.18, 0.14, 0.14], [0.21, 0.17, 0.14, 0.16, 0.1, 0.08, 0.11, 0.12, 0.17], [0.23, 0.26]]
strikes = [[13210.0, 13220.0, 13775.0, 13825.0, 13875.0], [12250.0, 13425.0, 13460.0, 13790.0, 13810.0, 13820.0, 13830.0, 13840.0, 13860.0, 13870.0, 13880.0, 13930.0, 14010.0, 14050.0], [8100.0, 10750.0, 11250.0, 11600.0, 12075.0, 12175.0, 12525.0, 13075.0, 13230.0, 13270.0, 13290.0, 13310.0, 13320.0, 13340.0, 13360.0, 13375.0, 13380.0, 13410.0, 13440.0, 13590.0, 13630.0, 13640.0, 13660.0, 13710.0, 13720.0, 13730.0, 13740.0, 13760.0, 13770.0, 13780.0, 13890.0, 13920.0, 13940.0, 13975.0, 13980.0, 13990.0, 14020.0, 14040.0, 14060.0, 14070.0, 14075.0], [13700.0, 14300.0], [11050.0, 13925.0, 14225.0, 14275.0, 14400.0], [13125.0, 13675.0], [15250.0], [12700.0, 12825.0, 13100.0, 14100.0, 14425.0, 14575.0, 15350.0], [13325.0], [10975.0, 14475.0, 15700.0], [12975.0, 13725.0], [12200.0, 12300.0, 12950.0, 13050.0, 13225.0, 13850.0, 13950.0, 14950.0, 15375.0, 15450.0, 15475.0, 15550.0, 15575.0, 15975.0, 16500.0], [9350.0, 10225.0, 12500.0, 13575.0], [14200.0, 15400.0, 17000.0, 17700.0], [16900.0], [13750.0, 14700.0, 15025.0, 15050.0, 16200.0, 17200.0, 17900.0, 18400.0, 19000.0], [12400.0], [8600.0, 11800.0, 13900.0, 14000.0, 14625.0, 15300.0, 15775.0], [15200.0, 15600.0, 16000.0, 17100.0, 18200.0, 20000.0, 20200.0, 20400.0], [13600.0, 16300.0, 16400.0, 17600.0, 19400.0, 19800.0], [7000.0, 8450.0, 10300.0, 11700.0, 12800.0, 13500.0, 14800.0, 16700.0, 19600.0], [7100.0, 13200.0]]
expiration_dates = [Date(16,12,2022), Date(15,12,2023), Date(16,4,2021), Date(17,12,2021), Date(18,6,2021), Date(7,5,2021), Date(23,4,2021), Date(21,5,2021), Date(14,4,2021), Date(15,10,2021), Date(30,4,2021), Date(14,5,2021), Date(26,4,2021), Date(12,4,2021), Date(3,5,2021), Date(17,6,2022), Date(19,4,2021), Date(17,9,2021), Date(16,7,2021), Date(18,3,2022), Date(28,4,2021), Date(20,8,2021)]
让我们从实现基础开始：
import QuantLib as ql
import numpy as np
import pandas as pd
from datetime import timedelta, date
day_count = ql.Actual365Fixed()
calendar = ql.UnitedStates()
tdy = date.today()
calculation_date = ql.Date(tdy.day,tdy.month,tdy.year)
spot = 13780.79296875
ql.Settings.instance().evaluationDate = calculation_date
risk_free_rate = 0.0
dividend_rate = 0.0
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))
成本函数和辅助函数由 表示：
def setup_helpers(engine, expiration_dates, strikes,
data, ref_date, spot, yield_ts,
dividend_ts,i):
heston_helpers = []
grid_data = []
date = expiration_dates[i]
l = 0
for s in strikes[i]:
t = (date - ref_date )
p = ql.Period(t, ql.Days)
vols = data[i][l]
helper = ql.HestonModelHelper(
p, calendar, spot, s,
ql.QuoteHandle(ql.SimpleQuote(vols)),
yield_ts, dividend_ts)
helper.setPricingEngine(engine)
heston_helpers.append(helper)
grid_data.append((date, s))
l =+ 1
return heston_helpers, grid_data
def cost_function_generator(model, helpers,norm=False):
def cost_function(params):
params_ = ql.Array(list(params))
model.setParams(params_)
avg = 0
for i, opt in enumerate(helpers):
err = (opt.modelValue()-opt.marketValue())**2
avg += err
if norm:
return avg/len(helpers)#MSE
else:
return error
return cost_function
def calibration_report(helpers, grid_data, detailed=False):
avg = 0.0
for i, opt in enumerate(helpers):
err = (opt.modelValue()-opt.marketValue())**2
avg += err
avg = avg/len(helpers) * 100 #in %
if detailed: print("-"*100)
summary = "Average MSE Error (%%) : %5.9f" % (avg)
print(summary)
return avg
def setup_model(_yield_ts, _dividend_ts, _spot,
init_condition=(0.02,0.2,0.5,0.1,0.01)):
theta, kappa, sigma, rho, v0 = init_condition
process = ql.HestonProcess(_yield_ts, _dividend_ts,
ql.QuoteHandle(ql.SimpleQuote(_spot)),
v0, kappa, theta, sigma, rho)
model = ql.HestonModel(process)
engine = ql.AnalyticHestonEngine(model)
return model, engine
此外，由于 Quantlib 不包含 Feller 条件，我将其合并为：
from scipy.optimize import NonlinearConstraint
def Cst(x):
theta, kappa, sigma, rho, v0 = x
return 2*kappa*theta - sigma**2
const = NonlinearConstraint(Cst,0.01,np.inf)
以下是对每对数据的优化：
from scipy.optimize import differential_evolution
summary = []
bounds = [(0.01,1),(0.01,20), (0.01,5.), (-1,1), (0.01,1.)]
for i in range(0,len(expiration_dates)):
model4, engine4 = setup_model(yield_ts, dividend_ts, spot)
heston_helpers4, grid_data4 = setup_helpers(
engine4, expiration_dates, strikes, data,
calculation_date, spot, yield_ts, dividend_ts,i)
initial_condition = list(model4.params())
cost_function = cost_function_generator(
model4, heston_helpers4, norm=True)
sol = differential_evolution(cost_function, bounds,constraints = const)
print(sol)
params_ = ql.Array(sol.x.tolist())
model4.setParams(params_)
theta, kappa, sigma, rho, v0 = model4.params()
print("theta = %f, kappa = %f, sigma = %f, rho = %f, v0 = %f" % \
(theta, kappa, sigma, rho, v0))
error = calibration_report(heston_helpers4, grid_data4)
theta, kappa, sigma, rho, v0 = list(model4.params())
summary.append([L['Time'].unique()[i], error] + list(model4.params()))
这里的问题是参数不稳定（我知道优化器是随机的。但是即使容差为 1e-6 ，参数每次也有显着差异）并且需要大量时间才能收敛（我什至没有确定确实如此）。此外，有时我会收到以下带有错误的消息（错误并不总是出现并且与约束有关，第二个奇怪的事情）：
/usr/local/lib/python3.7/dist-packages/scipy/optimize/_hessian_update_strategy.py:187: UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
'approximations.', UserWarning)
错误 ：
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-247-c3efacba5a02> in <module>()
----> 1 get_ipython().run_cell_magic('time', '', 'summary = []\ncost_function = cost_function_generator(\n model4, heston_helpers4, norm=True)\nsol = differential_evolution(cost_function, bounds, constraints = const,maxiter = 2000)\nparams_ = ql.Array(sol.x.tolist())\nprint(sol.x.tolist())\nmodel4.setParams(params_)\ntheta, kappa, sigma, rho, v0 = model4.params()\nprint("theta = %f, kappa = %f, sigma = %f, rho = %f, v0 = %f" % \\\n (theta, kappa, sigma, rho, v0))\nerror = calibration_report(heston_helpers4, grid_data4)\nsummary.append(["Scipy DE1", error] + list(model4.params()))')
16 frames
<decorator-gen-53> in time(self, line, cell, local_ns)
<timed exec> in <module>()
/usr/local/lib/python3.7/dist-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
484 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
485 raise ValueError(
--> 486 "array must not contain infs or NaNs")
487 return a
488
ValueError: array must not contain infs or NaNs
有人可以帮助我以使这个棘手的优化工作吗？
谢谢，