用模拟退火算法估价heston期权定价模型的五个参数

特此声明：笔者发现这个帖子需求量较大，于是重新整理的更加完整和可操作，有需要请查看
https://blog.csdn.net/weixin_45590329/article/details/119610079

1973年BS期权定价模型的诞生标志着期权定价进入精确的数量化测度阶段。但是BS模型假设标的资产波动率为常数，这与现实市场观测到的“波动率微笑”曲线严重不符。
heston假设标的资产的价格服从如下过程，其中波动率为时变函数[1]：

并且求出了欧式看涨期权定价公式[2]：


本文使用python实现了上述定价公式。该公式需要输入一共九个参数，其中[v0,kappa,theta,sigma,rho]需要提前自行设置并填入。另外四个 [K,t,s0,r]：[执行价格，剩余时间，标的资产价格，无风险利率] 则是期权价格数据，BS模型定价同样需要这几个参数。
sv_model.py

"""heston期权定价模型计算期权价格"""
import numpy as np
import scipy.integrate as spi

class SV():
    """计算Heston随机波动率模型下的期权价格"""

    def __init__(self,K:float,t:float,s0:float,r:float,v0:float=0.01,kappa:float=2,theta:float=0.1,sigma:float=0.1,rho:float=-0.5):
        """输入期权数据和基础参数
        
        #期权数据
        k   -期权执行价格
        t  -期权剩余到期时间（年）
        s0 -标的资产初始价格
        r  -无风险收益率
        
        #初始参数
        v0
        kappa
        theta
        sigma
        rho
        """
        # 期权数据
        self.K = K 
        self.t = t 
        self.s0 = s0 
        self.r = r

        # 初始参数
        self.v0 = v0
        self.kappa = kappa
        self.theta = theta
        self.sigma = sigma
        self.rho = rho

    # 特征函数
    def characteristic_function(self,phi, type):

        if type == 1:
            u = 0.5
            b = self.kappa - self.rho * self.sigma
        else:
            u = -0.5
            b = self.kappa

        a = self.kappa * self.theta
        x = np.log(self.s0)
        d = np.sqrt((self.rho * self.sigma * phi * 1j - b) ** 2 - self.sigma ** 2 * (2 * u * phi * 1j - phi ** 2))
        g = (b - self.rho * self.sigma * phi * 1j + d) / (b - self.rho * self.sigma * phi * 1j - d)
        D = self.r * phi * 1j * self.t + (a / self.sigma ** 2) * (
                (b - self.rho * self.sigma * phi * 1j + d) * self.t - 2 * np.log((1 - g * np.exp(d * self.t)) / (1 - g)))
        E = ((b - self.rho * self.sigma * phi * 1j + d) / self.sigma ** 2) * (1 - np.exp(d * self.t)) / (1 - g * np.exp(d * self.t))

        return np.exp(D + E * self.v0 + 1j * phi * x)

    # 积分部分
    def integral_function(self,phi,type):

        integral = (np.exp(-1 * 1j * phi * np.log(self.K)) * self.characteristic_function(phi, type=type))

        return integral

    # p值函数
    def sv_P_Value(self,type):

        """

        """
        ifun = lambda phi: self.integral_function(phi,type=type) / (1j * phi)
        return 0.5 + (1 / np.pi) * spi.quad(ifun, 0, 1000)[0]

    def sv_Call_Value(self):

        """

        """
        P1 = self.s0 * self.sv_P_Value(type=1)
        P2 = self.K * np.exp(-self.r * self.t) * self.sv_P_Value(type=2)

        if np.isnan(P1-P2):return 1000000#如果初始参数使定价结果为nan，就返还巨大的值，视为报错！巨大的期权价格会拉大定价误差
        else:return P1 - P2

    #方便检测定价结果
    @classmethod
    def test_multiple(cls,option_params:list,init_params:list):
        """
        输入：
        option_params   -list,期权价格参数
                        [K,t,s0,r]
        init_params     -list,sv模型参数
                        [v0,kappa,theta,sigma,rho]
        输出：
            sv模型期权价格

        示例：
        sv=SV(K =2.15,t = 0.05924,s0 = 2.143896,r = 0.03043)
        sv.test_multiple([2.15,0.059,2.143896,0.03043],[0.01,2,0.1,0.1,-0.5])#通过直接输入两个列表，很方便的计算期权价格
        """
        sv=cls(option_params[0],option_params[1],option_params[2],option_params[3],
               init_params[0],init_params[1],init_params[2],init_params[3],init_params[4])
        c_sv=sv.sv_Call_Value()
        return c_sv


if __name__=='__main__':
    sv=SV(K =2.15,t = 0.05924,s0 = 2.143896,r = 0.03043)
    sv.sv_Call_Value()

由于heston模型中有五个参数需要估计，传统的线性最小二乘法失效。本文采用退火模拟算法估计这五个参数。说实话，这五个参数非常难估计。
首先写出模拟退火算法的类：
myNG.py

"""模拟退火算法"""
import numpy as np
import matplotlib.pyplot as plt
import copy
from randomrange import random_range

"""
部分名词：
新解：在某轮循环中通过自变量随机变动产生的新解
新接受解：由于新解导致函数值变小，或者由于概率而接受的新解
历史最优解：在整个循环中产生的最优解
"""


class NG():
    def __init__(self,func,x0):
        """
        x0  -列表：函数的初始值
                [v0,kappa,theta,sigma,rho]
        func    -待求解函数
        """
        self.x=x0
        self.dim_x=len(x0)#解的维度
        self.func=func
        self.f = func(self.x)#计算y值

        self.x_best=self.x#记录下来历史最优解，即所有循环中的最优解
        self.f_best=self.f

        self.times_stay=0#连续未接受新解的次数
        self.times_stay_max=400#当连续未接受新解超过这个次数，终止循环

        self.T=100#初始温度：初始温度越高，前期接受新解的概率越大
        self.speed=0.7#退火速度：温度下降速度越快，后期温度越低，接受新解的概率越小，从而最终稳定在某个解附件
        self.T_min=1e-6#最低温度：当低于该温度时，终止全部循环
        self.xf_best_T={}#记录下接受的所有新解

        # 最初若函数值变动为delta，则认为函数值变动很大，可以产生p_expec概率接受新解
        #若在初期便产生巨大的概率接受新解，则前期寻找新解的过程将变成盲目的随机漫步毫无意义，因此利用alpha调节概率
        self.p_expec = 0.9
        self.delta_standard = 0.7
        self.alpha =self.find_alpha()#调节概率因子

        self.times_delta_samller = 0#统计新旧最优值之差绝对值连续小于某值的次数
        self.delta_min=0.001#当新旧最优值之差绝对值连续小于此值达到某一次数时，终止该温度循环
        self.times_delta_min=100#当新旧最优值之差绝对值连续小于此值达到这个次数时，终止该温度循环

        self.times_max=500#当每个温度下循环超过这个次数，终止该温度循环
        self.times_cycle=0#记录算法循环总次数

        self.times_p=0#统计因为p值较大而接受新解的次数

        self.xf_all={self.times_cycle:[self.x,self.f]}#记录下来每一次循环产生的新解和函数值
        self.xf_best_all={self.times_cycle:[self.x,self.f]}#记录下来每一次循环接受的新解和函数值

    #温度下降，产生新温度
    def T_change(self):
        self.T=self.T*self.speed
        print('当前温度为{},大于最小温度{}'.format(self.T,self.T_min))#展示当前温度和最小温度

    # 将所有的x和f、循环次数存储下来
    def save_xy(self):
        self.xf_all[self.times_cycle]=[self.x,self.f]

    #将所有的最优x,y、循环次数存储下来
    def save_best_xy(self):
        self.xf_best_all[self.times_cycle]=[self.x,self.f]

    # 当调节因子为alpha时，函数值变动值为delta产生的接受新解概率
    def __p_delta(self,alpha):
        return np.exp(-self.delta_standard / (self.T * alpha))

    # 用二分法寻找方程解
    def __find_solver(self,func, f0):
        """
        输入：
        func    -待求解方程的函数
        f0  -float,预期函数值
        输出：
        mid -float,函数=预期函数值 的解
        """
        up = 100
        down = 0.00001
        mid = (up + down) / 2
        while abs(func(mid) - f0) > 0.0001:
            if func(down) < f0 < func(mid):
                up = mid
                mid = (mid + down) / 2
            elif func(up) > f0 > func(mid):
                down = mid
                mid = (up + down) / 2
            else:
                print('error!')
                break
        return mid

    # 最初若函数值变动为delta，则认为函数值变动很大，可以产生p_expec概率接受新解
    def find_alpha(self):
        return  self.__find_solver(self.__p_delta, self.p_expec)

    #获得新的x
    def get_x_new(self):
        random=np.random.normal(0,1,self.dim_x)#新的随机增量
        return   self.x+random

    #判断是否可以接受新的解
    def judge(self):
        if self.delta<0:#如果函数值变动幅度小于0,则接受新解
            self.x=self.x_new
            self.f_last=self.f#在最优解函数值更新之前将其记录下来
            self.f=self.f_new
            self.save_best_xy()#记录每次循环接受的新解
            self.get_history_best_xy()#更新历史最优解
            self.times_stay = 0  # 由于未接受新解，将连续未接受新解的次数归零
            print('由于函数值变小新接受解{}:{}'.format(self.f,self.x))#展示当前接受的新解
        else:
            p=np.exp(-self.delta/(self.T*self.alpha))#接受新解的概率
            p_=np.random.random()#判断标准概率
            if p>p_:#如果概率足够大，接受新解
                self.x = self.x_new
                self.f_last = self.f  # 在接受的新解更新之前将其记录下来
                self.f = self.f_new
                self.save_best_xy()#记录每次循环接受的新解
                self.get_history_best_xy()  # 更新历史最优解
                print('由于概率{}大于{}，新接受解{}:{}'.format(p,p_,self.f, self.x))#展示当前接受的新解
                self.times_p+=1#统计因为概率而接受新解的次数
                self.times_stay=0#由于未接受新解，将连续未接受新解的次数归零
            else:
                self.times_stay+=1#连续接受新解次数加1
                print('连续未接受新解{}次'.format(self.times_stay))

    #获得历史最优解
    def get_history_best_xy(self):
        x_array = list(np.array(list(self.xf_best_all.values()))[:, 0])  # 从历史所有的最优x和f中获得所有的x
        f_array=list(np.array(list(self.xf_best_all.values()))[:, 1])#从历史所有的最优x和f中获得所有的f
        self.f_best=min(f_array)#从每阶段最优的f中获得最优的f
        self.x_best=x_array[f_array.index(self.f_best)]#利用最优f反推最优x
        return self.x_best,self.f_best

    #绘制函数值降低记录
    def plot(self,x,y):
        plt.figure(figsize=(13,8))
        plt.plot(x,y)
        plt.xlabel('循环次数',fontsize=15)
        plt.ylabel('函数值', fontsize=15)
        plt.title('函数值变化过程')

    #绘制最优值变化过程
    def plot_best(self):
        times_cycle_array=list(self.xf_best_all.keys())#获得接受新解时的循环次数
        f_array = np.array(list(self.xf_best_all.values()))[:, 1]  # 从所有接受的新解中获得所有的函数值
        self.plot(times_cycle_array,f_array)#绘制循环次数与函数值的折线图

    #统计新旧函数值之差的绝对值连续小于此值的次数
    def count_times_delta_smaller(self):
        if self.delta_best<self.delta_min:
            self.times_delta_samller+=1#如果新旧函数值之差绝对值小于某值，则次数加1，否则归零
        else:
            self.times_delta_samller=0
        print('差值{}连续小于{}达到{}次'.format(self.delta_best,self.delta_min,self.times_delta_samller))

    #终止循环条件
    def condition_end(self):
        if self.times_delta_samller>self.times_delta_min:#如果新旧函数值之差绝对值连续小于某值次数超过某值，终止该温度循环
            return True
        elif self.times_stay>self.times_stay_max:#当连续未接受新解超过这个次数，终止循环
            return True

    #在某一特定温度下进行循环
    def run_T(self):

        for time_ in range(self.times_max):
            self.x_new=self.get_x_new()#获得新解
            self.f_new=self.func(self.x_new)#获得新的函数值
            self.save_xy()#将新解和函数值记录下来
            self.delta=self.f_new-self.f#计算函数值的变化值
            self.judge()#判断是否接受新解
            self.times_cycle+=1#统计循环次数
            self.delta_best=np.abs(self.f-self.f_last)#上次函数值与这次函数值的差值绝对值
            self.count_times_delta_smaller()#统计新旧函数值之差的绝对值连续小于此值的次数
            if self.condition_end()==True:#如果满足终止条件，终止该温度循环
                print('满足终止条件：接受新解后的函数值变化连续小于{}达到次数'.format(self.delta_min))
                break
            print('当前历史最优解{}：{}'.format(self.f_best,self.x_best))#展示当前最优值
            print('当前接受的新解{}：{}'.format(self.f, self.x))  # 展示当前接受的新解
            print('当前新解{}：{}'.format(self.f_new, self.x_new))  # 展示当前新产生的解
            print('当前温度为{}'.format(self.T))#展示当前温度

    #当每个温度下的循环结束时，有一定概率将当前接受的新解替换为历史最优解
    def accept_best_xf(self):
        if np.random.random()>0.75:
            self.x=self.x_best
            self.f=self.f_best

    def run(self):
        while self.T>self.T_min:
            self.run_T()#循环在该温度下的求解
            self.xf_best_T[self.T] = [self.get_history_best_xy()]#记录在每一个温度下的最优解
            self.T_change()#温度继续下降
            self.accept_best_xf()#当每个温度下的循环结束时，有一定概率将当前接受的新解替换为历史最优解
            if self.condition_end()==True:#如果满足终止条件，终止该温度循环
                break

if __name__=='__main__':
    # 待检测函数1
    def func_target(list_x):
        x1 = list_x[0]
        x2 = list_x[1]
        return 3 * (x1 - 5) ** 2 + 6 * (x2 - 6) ** 2 - 7

    ng=NG(func=func_target,x0=[8,9])
    ng.run()
    x,f=ng.get_history_best_xy()#查看历史最优解
    ng.plot_best()#绘制最优解变化过程

这是个一般性的模拟退火算法类，解决如上所示的二次凸函数最优化非常简单。但是在解决高度复杂的heston模型时就非常困难了。为了解决heston模型的参数估计，基于以上模型专门写一个模拟退火算法子类。该子类的作用主要是将五个参数限制在一定范围内，因为在参数估计过程中发现部分参数取值过大时会出现 nan 类型的python无法表示的浮点数

randomrange.py

"""某一个数在指定范围内随机变动"""
import matplotlib.pyplot as plt
import numpy as np


#将x增加一个随机变动量，但会把x限制在a,b之内
def random_range(x,a,b):
    random=np.random.normal(0,0.01*(b-a),1)[0]
    if x+random>b :random=-np.abs(random)#如果新值超越最大值，就将增量变为负数
    elif x+random<a:random=np.abs(random)#如果新值超越最小值，就将增量变为正数
    return x+random


if __name__=='__main__':

    x=0.3;a=0;b=4
    sum_False=0
    x_all=[]
    for i in range(1000):
        x=random_range(x,a,b)
        x_all.append(x)
        print(x)
        if x>b or x<a:
            sum_False+=1

    plt.plot(x_all)

myNG.py

#专门用于求解sv模型参数的退火算法
class NGSV(NG):
    def __init__(self,func,x0):
        super().__init__(func,x0)
        self.T=90
        self.T_min=1e-7#由于算法耗时太长，故小做一段模拟试试看
        self.times_max=500

    #sv模型的各个参数由于存在取值范围，因此在获得新的参数估计值时需要对其取值范围加以限制
    def get_x_new(self):
        """
        [v0,kappa,theta,sigma,rho]
        其中：
        v0,kappa,theta,sigma>0
        -1<rho<1
        2kappa*theta>sigma**2
        """
        x=copy.deepcopy(self.x)#使用深copy，否则self.x会随着x一起变动
        x[0]=random_range(x[0],0,5)
        x[1] =random_range(x[1],0,1)
        x[2] =random_range(x[2],0,1)
        x[3] =random_range(x[3],0,3)
        x[4]=random_range(x[4],-1,1)
        return x

利用上面的heston模型定价类和模拟退火算法类，写一个利用模拟退火算法解决heston模型参数估计的类
sv_sa.py

"""使用模拟退火算法求解期权价格"""
from sv_model import SV
#from sko.SA import SA
import pandas as pd
import numpy as np
from scipy.optimize import minimize
from myNG import NG,NGSV
import copy


#读取初始数据
# data=pd.read_excel('期权数据.xlsx')
# data_option=data.ix[(data['call_put']=='C')&(data['type_in_out']==4)&(data['type_time']==3)][['exercise_price','time','ETF50','shibor','close']]
#
# data_option.columns=['K','t','s0','r','c']
# data_option.index=range(len(data_option))
# data_option.to_excel('待测试数据.xlsx')
data_option=pd.read_excel('待测试数据.xlsx',index_col=0)
data_option=data_option.ix[:100,:]


class SV_SA():

    def __init__(self,data,v0:float=0.01,kappa:float=2,theta:float=0.1,sigma:float=0.1,rho:float=-0.5):
        """输入数据

        data    -pandas.core.frame.DataFrame格式数据，具体样式如下：

                           K         t        s0         r       c
                30     2.150  0.194444  2.111919  0.031060  0.0546
                31     2.150  0.198413  2.115158  0.031120  0.0666
                32     2.150  0.202381  2.107673  0.031210  0.0627
                33     2.150  0.214286  2.122269  0.031250  0.0531
                90     3.240  0.202381  3.181339  0.047446  0.0724

        """

        self.data=data
        self.init_params=[v0,kappa,theta,sigma,rho]# 初始参数列表
        self.cycle=0#计算模拟退火算法轮数
        self.error=0.000000

    def error_mean_percent(self,init_params:list):
        """计算heston模型期权定价的百分比误差均值

        百分比误差均值=绝对值（（理论值-实际值）/实际值）/样本个数

        输入：
        init_params -初始参数,列表格式
                     [v0,kappa,theta,sigma,rho]

        返回： -误差百分点数   例如：返回5，表示5%
        """
        v0,kappa,theta,sigma,rho=init_params
        list_p_sv=[]
        for i in self.data.index:
            K,t,s0,r,p_real=self.data.ix[i,:].tolist()
            sv = SV(K=K, t=t, s0=s0, r=r,
                    v0=v0, kappa=kappa, theta=theta, sigma=sigma, rho=rho)
            p_sv = sv.sv_Call_Value()  # sv模型期权价格
            list_p_sv.append(p_sv)

        self.error = np.average(np.abs((np.array(list_p_sv) - self.data['c']) / self.data['c']))#sv模型的期权价格和实际价格的百分比误差均值
        print('\n')
        print('第{}轮,误差：{}'.format(self.cycle, self.error))#展示本轮的误差
        self.cycle += 1

        return self.error


    def error_mean(self,init_params:list):
        """计算heston模型期权定价的均方误差
        init_params -初始参数,列表格式
                     [v0,kappa,theta,sigma,rho]
        """
        v0,kappa,theta,sigma,rho=init_params
        list_p_sv=[]
        for i in self.data.index:
            K,t,s0,r,p_real=self.data.ix[i,:].tolist()
            sv = SV(K=K, t=t, s0=s0, r=r,
                    v0=v0, kappa=kappa, theta=theta, sigma=sigma, rho=rho)
            p_sv = sv.sv_Call_Value()  # sv模型期权价格
            list_p_sv.append(p_sv)

        self.error=np.sqrt(np.sum((np.array(list_p_sv)-self.data['c'])**2)/len(self.data))#sv模型的期权价格和实际价格的均方误差
        print('\n')
        print('第{}轮,误差：{}'.format(self.cycle, self.error))#展示本轮的误差
        self.cycle += 1

        return self.error

    def test_error_mean(self,multiple_parmas:dict):
        """将多组初始参数输入，计算各组参数的均方误差
        multiple_parmas -dict,多组初始参数
                        {
                        1:[0.01,2,0.1,0.1,-0.5],
                        2:[0.01,2,0.1,0.1,-0.5],
                        3:[0.01,2,0.1,0.1,-0.5]
                        }
        返回： -dict,记录各组初始参数的均方误差
        """
        dict_={}#用于记录各组初始参数的均方误差
        for i in multiple_parmas.keys():
            dict_[i]=self.error_mean(multiple_parmas[i])
        return dict_

    def test_option_price(self,multiple_parmas:dict):
        """将多组期权数据和初始参数输入，将期权价格合并在表格旁边
                multiple_parmas -dict,多组初始参数
                                multiple_parmas={
                                1:[1.5932492661058346, 3.3803420203705365, 0.3333248435472669, 5.622092726036617, 0.044881506437356666],
                                2:[1.1070063457234607, 3.501301312245266, 0.6276009140316863, 9.383112611111134, -0.6092511548040354],
                                3:[0.5675877305927083, 3.736229838972323, 0.21803303626214537, 8.74231319248172, 0.09393882921335006]
                                }
                返回： -dict,记录各组初始参数的均方误差
                """
        data_option_=copy.deepcopy(data_option)
        for i in multiple_parmas.keys():
            #i=3
            data_option_['第{}组参数'.format(i)]=0.000000
            init_params = multiple_parmas[i]
            for j in data_option_.index:
                #j=8
                option_params=data_option_.ix[j,:4].tolist()
                sv=SV(option_params[0],option_params[1],option_params[2],option_params[3],
                 init_params[0],init_params[1],init_params[2],init_params[3],init_params[4])
                c_sv=sv.sv_Call_Value()
                data_option_.ix[j,'第{}组参数'.format(i)]=c_sv
            print('已经完成第{}组'.format(i))
        return  data_option_

    def sa(self):
        """对均方误差函数用模拟退火算法计算最优值
        """
        # 实例化算法，并加入初始解
        # opt = minimize(self.error_mean, self.init_params, method='Nelder-Mead', tol=1e-6)#单纯形法求最优解
        # self.best_params=opt.x

        #self.x_star, self.y_star,self.list_,self.info = MySA(self.error_mean, self.init_params)
        self.ng=NGSV(func=self.error_mean_percent,x0=self.init_params)
        self.ng.run()
        self.x_star,self.y_star=self.ng.get_history_best_xy()

        print(self.x_star, self.y_star)# 生成最优解x和最优值y


if __name__=='__main__':
    model=SV_SA(data=data_option)
    model.sa()
    model.ng.get_history_best_xy()

所有代码及数据：
https://download.csdn.net/download/weixin_45590329/13082047

笔者发现这个帖子需求量较大，于是重新整理的更加完整和可操作，有需要请查看
https://blog.csdn.net/weixin_45590329/article/details/119610079