0

我正在尝试使用库 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

有人可以帮助我以使这个棘手的优化工作吗?

谢谢,

4

0 回答 0