优化，使用 Quantlib 和 Scipy 差分进化，不收敛

如何解决优化，使用 Quantlib 和 Scipy 差分进化，不收敛

我正在尝试使用库 quantlib 和作为 DE 的全局优化器校准 heston 模型的参数。

首先，我们有以下隐含波动率、执行价格和到期日期列表：

data = [[0.22,0.22,0.19,0.14,0.14],[0.16,0.2,0.15,0.11,0.16,0.1,0.23,0.17,0.08],0.21,0.18,0.08,0.13,0.12,0.16],[0.2,[0.25,0.2],[0.14],[0.21,0.11],[0.23,0.25],0.19],[0.12,0.17],[0.1,0.18],[0.2],[0.18,[0.16],0.13],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,2023),Date(16,4,2021),Date(17,Date(18,6,Date(7,5,Date(23,Date(21,Date(14,10,Date(30,Date(26,Date(12,Date(3,Date(19,9,7,3,Date(28,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,s,ql.QuoteHandle(ql.SimpleQuote(vols)),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.5,0.01)):
    theta,kappa,sigma,rho,v0 = init_condition
    process = ql.HestonProcess(_yield_ts,ql.QuoteHandle(ql.SimpleQuote(_spot)),v0,theta,rho)
    model = ql.HestonModel(process)
    engine = ql.AnalyticHestonEngine(model) 
    return model,engine

此外，由于 Quantlib 没有包含 Feller 条件，因此我将其合并，例如：

from scipy.optimize import NonlinearConstraint

def Cst(x):
    theta,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),5.),(-1,1.)]

for i in range(0,len(expiration_dates)):
    model4,engine4 = setup_model(yield_ts,spot)
    heston_helpers4,grid_data4 = setup_helpers(
    engine4,calculation_date,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,v0 = model4.params()
    print("theta = %f,kappa = %f,sigma = %f,rho = %f,v0 = %f" % \
    (theta,v0))
    error = calibration_report(heston_helpers4,grid_data4)
    theta,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,norm=True)\nsol = differential_evolution(cost_function,constraints = const,maxiter = 2000)\nparams_ = ql.Array(sol.x.tolist())\nprint(sol.x.tolist())\nmodel4.setParams(params_)\ntheta,v0 = model4.params()\nprint("theta = %f,v0 = %f" % \\\n    (theta,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

有人可以帮我实现这个棘手的优化吗？

谢谢，