SVI模型拟合

Posted on September 9, 2018 by laofish

波动率模型是学界搞期权的重点，BS模型中，使用的是固定的波动率。现在提的较多的是局部波动率模型和随机波动率模型。局部波动率模型是随机波动率模型的一种简化，将波动率定义为标的价格S和时间t的函数sigma(S,t)，随机波动率模型中，对波动率描述存在随机项（独立于价格中的随机项）。

局部波动率模型中，常用的是SVI模型(stochastic volatility inspired)，SVI模型本质上是描述了波动率微笑曲线，并且增加了一些期权性质上的约束。

我前面研究过一次SVI模型，但由于对python里面最优化函数理解的不透，对外层使用了是穷举法，这个方法显然是不好的，只是让我初次了解这个模型的性质。最近部门一个实习生用matlab实现了一下这个模型，虽然拟合的结果并不是很好，但让我发现原来matlab里面的几个函数，如lsqlin就是二次非线性最优化，而fminsearch的函数，使用的就是Nelder-Mead Simplex方法。

实际上我以前虽然也用matlab做最优化，但是对里面的求解算法实在是没研究，反正直接用就好了。但前面自己用python摸索求解函数的时候，突然觉得自己根本不知道在python里面用什么函数去求解自己要解的最优化方程。

最初以为scipy里面应该有对应的最优化函数，后来发现matlab里面基本的lsqlin在scipy里面就没有实现，scipy最接近的是lsq_liner，但这个函数只能接受上下限，无法接受不等式方程约束，后来发现用凸优化报cvxopt才有希望。但是这个包里面也没有lsqlin函数，不过知识点到这里，搜索一下lsqlin+python关键词就能找到一个俄罗斯人贡献的用python实现的lsqlin代码了。

在这些知识点的基础上，我重新将自己的svi模型实现了一下。

首先是数据源格式

  1. 合约编码    行权日 行权日 剩余到期时间(日历日) 期权类型    标的收盘价   期权行权价   收盘中间价 
  2. 10000887    20171227    2017-12-27  0.052082192 C   2.865   2.209   0.6624 
  3. 10000888    20171227    2017-12-27  0.052082192 C   2.865   2.258   0.6135 
  4. 10000889    20171227    2017-12-27  0.052082192 C   2.865   2.307   0.56465 
  5. 10000890    20171227    2017-12-27  0.052082192 C   2.865   2.356   0.5158 
  6. 10000891    20171227    2017-12-27  0.052082192 C   2.865   2.405   0.46705 
  7. 10000892    20171227    2017-12-27  0.052082192 P   2.865   2.209   0.00015 
  8. 10000893    20171227    2017-12-27  0.052082192 P   2.865   2.258   0.00015 
  9. 10000894    20171227    2017-12-27  0.052082192 P   2.865   2.307   0.00015 
  10. 10000895    20171227    2017-12-27  0.052082192 P   2.865   2.356   0.0002 
  11. 10000896    20171227    2017-12-27  0.052082192 P   2.865   2.405   0.0003 
  12. 10000897    20171227    2017-12-27  0.052082192 C   2.865   2.16    0.7117 

其次是拟合代码

  1. #!/usr/bin/env python 
  2. # encoding: utf-8 
  3. """ 
  4. @version: ?? 
  5. @author: laofish 
  6. @contact: laofish@outlook.com 
  7. @site: http://www.laofish.com 
  8. @file: fitSVI.py 
  9. @time: 2018-09-07 23:48 
  10.  
  11. 拟合SVI曲线 
  12.  
  13. """ 
  14.  
  15. import traceback 
  16.  
  17. from scipy.optimize import lsq_linear 
  18. from scipy.optimize import fmin 
  19.  
  20. from lsqlin import * 
  21.  
  22. from derivaties_fun2 import * 
  23.  
  24. from scipy.optimize import least_squares 
  25.  
  26. # 中文和负号的正常显示 
  27. import socket 
  28.  
  29. hostName = socket.gethostname() 
  30.  
  31. if hostName == 'laofish-home-PC': 
  32.     plt.rcParams['font.sans-serif'] = ['Yahei Mono']  # mat 2.1 
  33. else: 
  34.     plt.rcParams['font.sans-serif'] = ['Microsoft Yahei Mono'] 
  35.  
  36. plt.rcParams['axes.unicode_minus'] = False 
  37.  
  38. # 设置图形的显示风格 
  39. plt.style.use('ggplot') 
  40.  
  41.  
  42. def neicengOptimization(guess, myargs): 
  43.     '''内侧优化使用''' 
  44.     # SVI模型内层优化 
  45.     # m,sigma 为给定值，直接从参数里面得到 
  46.     # x 为自变量,即最优化问题里面的x 
  47.     # a d c 为系数，即最优化问题里面的c 
  48.  
  49.     m, sigma = guess 
  50.  
  51.     x = myargs[0] 
  52.  
  53.     omg_i = myargs[1] 
  54.  
  55.     # x = lsqlin(C,d,A,b) solves the linear system C*x = d 
  56.     # in the least-squares sense subject to A*x ≤ b, where C is m-by-n. 
  57.  
  58.     yx = (x - m) / sigma 
  59.     zx = np.sqrt(yx ** 2 + 1) 
  60.  
  61.     omega = max(omg_i) 
  62.  
  63.     # xz = np.array([np.ones([len(x), 1]).T, yx, zx]) 
  64.     xz = np.array([np.array(np.ones([len(x), 1]).T.tolist()[0]), yx, zx]) 
  65.  
  66.     # 实际上可能可以用cvxpy这个凸优化包的使用更为简单 
  67.     A = [[0, 0, -1], [0, 0, 1], [0, -1, -1], [0, 1, -1], [0, 1, 1], [0, -1, 1], [-1, 0, 0], [1, 0, 0]] 
  68.     # 约束边界b 
  69.  
  70.     b = [0, 4 * sigma, 0, 0, 4 * sigma, 4 * sigma, 0, omega] 
  71.     # lsqlin(C, d, reg=0, A=None） 
  72.     xmatric = lsqlin(xz.T, omg_i, 0, np.array(A), np.array(b)) 
  73.     acap = np.array(xmatric['x'])[0][0] 
  74.     d = np.array(xmatric['x'])[1][0] 
  75.     c = np.array(xmatric['x'])[2][0] 
  76.     # [acap, d, c]=lsqlin(xz.T, omg_i, 0, np.array(A), np.array(b)) 
  77.  
  78.     # [acap,d,c] = lsqlin(xz, omg_i,A,b) 
  79.  
  80.     a = acap 
  81.     b = c / sigma 
  82.     rho = d / c 
  83.  
  84.     # 目标函数 
  85.     sigma2 = np.sum(np.array(acap + d * yx + c * zx - omg_i) ** 2) 
  86.  
  87.     # return [a,b,rho,sigma2] 
  88.  
  89.     return sigma2 
  90.  
  91.  
  92. def neicengOptimization2(guess, myargs): 
  93.     '''返回a d c ''' 
  94.     m, sigma = guess 
  95.     x = myargs[0] 
  96.     omg_i = myargs[1] 
  97.  
  98.     yx = (x - m) / sigma 
  99.     zx = np.sqrt(yx ** 2 + 1) 
  100.  
  101.     omega = max(omg_i) 
  102.  
  103.     # xz = np.array([np.ones([len(x), 1]).T, yx, zx]) 
  104.     xz = np.array([np.array(np.ones([len(x), 1]).T.tolist()[0]), yx, zx]) 
  105.  
  106.     A = [[0, 0, -1], [0, 0, 1], [0, -1, -1], [0, 1, -1], [0, 1, 1], [0, -1, 1], [-1, 0, 0], [1, 0, 0]] 
  107.     # 约束边界b 
  108.  
  109.     b = [0, 4 * sigma, 0, 0, 4 * sigma, 4 * sigma, 0, omega] 
  110.     # lsqlin(C, d, reg=0, A=None） 
  111.     xmatric = lsqlin(xz.T, omg_i, 0, np.array(A), np.array(b)) 
  112.     acap = np.array(xmatric['x'])[0][0] 
  113.     d = np.array(xmatric['x'])[1][0] 
  114.     c = np.array(xmatric['x'])[2][0] 
  115.     # [acap, d, c]=lsqlin(xz.T, omg_i, 0, np.array(A), np.array(b)) 
  116.  
  117.     # [acap,d,c] = lsqlin(xz, omg_i,A,b) 
  118.  
  119.     a = acap 
  120.     b = c / sigma 
  121.     rho = d / c 
  122.  
  123.     # 目标函数 
  124.     # sigma2 = np.sum(np.array(acap + d * yx + c * zx - omg_i) ** 2) 
  125.     return [a, b, rho] 
  126.  
  127.  
  128. def getfitsiv(x, t2M, a_star, b_star, rho_star, m_star, sigma_star): 
  129.     '''生成拟合后的SVI波动率''' 
  130.  
  131.     # 计算拟合值 
  132.     fit_omg = a_star + b_star * (rho_star * (x - m_star) + np.sqrt((x - m_star) ** 2 + sigma_star ** 2)) 
  133.  
  134.     fsigma = np.sqrt(fit_omg / t2M) 
  135.  
  136.     return fsigma 
  137.  
  138.  
  139. if __name__ == '__main__': 
  140.  
  141.     optionData = pd.read_excel('rawOptionData.xlsx') 
  142.  
  143.     strikeArr = pd.unique(optionData['期权行权价']) 
  144.     time2Mature = pd.unique(optionData['剩余到期时间(日历日)']) 
  145.     ulprice = pd.unique(optionData['标的收盘价']) 
  146.  
  147.     r = 0.04 
  148.     S = ulprice[0] 
  149.     q = 0 
  150.  
  151.     # 鉴于国内特色，只计算call 
  152.     oD = optionData 
  153.     vArr = [] 
  154.     strikeArrE = [] 
  155.     time2MatureE = [] 
  156.     cData = pd.DataFrame() 
  157.     try: 
  158.         # 处理错误 
  159.  
  160.         # step 1 
  161.         # 计算期权市场的隐含波动率 
  162.         for strike in strikeArr: 
  163.             for t2M in time2Mature: 
  164.                 oprice = oD[(oD['期权行权价'] == strike) & (oD['期权类型'] == 'C') 
  165.                             & (oD['剩余到期时间(日历日)'] == t2M)]['收盘中间价'] 
  166.                 if oprice.empty == False: 
  167.                     cimpv = bsmImpVol(S, strike, t2M, r, q, oprice.tolist()[0], 'C') 
  168.                 else: 
  169.                     cimpv = 0 
  170.  
  171.                 # 这个就是隐含方差 
  172.                 # total implied variance 
  173.                 v = cimpv ** 2 * t2M 
  174.  
  175.                 vArr.append(v) 
  176.                 strikeArrE.append(strike) 
  177.                 time2MatureE.append(t2M) 
  178.  
  179.         cData['K'] = strikeArrE 
  180.         cData['t2M'] = time2MatureE 
  181.         cData['v'] = vArr 
  182.  
  183.         # 去掉v等于0的元素 
  184.         cData2 = cData[cData['v'] != 0] 
  185.  
  186.         # 论文中没提到的是，这里应该一个到期时间一个到期时间的循环 
  187.         for t2M in cData2['t2M'].unique(): 
  188.             # 扣除所需数据 
  189.             thisData = cData2[cData2['t2M'] == t2M] 
  190.             x = np.log(thisData['K'] / (S * np.exp(r * t2M))) 
  191.  
  192.             omg_i = thisData['v'].values 
  193.  
  194.             # %para为m和sigma的初值 
  195.             m = 0.1 
  196.             sigma = 0.2 
  197.  
  198.             myargs = np.array([x, omg_i])  # 需要传进去的参数 
  199.  
  200.             # arg_out= fmin(neicengOptimization, para) 
  201.             # neicengOptimization(m, sigma, myargs) 
  202.             # arg_out = fmin(neicengOptimization, x0=[m, sigma],args=(myargs,)) 
  203.             guess = [m, sigma] 
  204.             arg_out = fmin(neicengOptimization, guess, args=(myargs,)) 
  205.  
  206.             # m sigma 
  207.             m_star, sigma_star = arg_out 
  208.             # 获取a,b,rho 
  209.             a_star, b_star, rho_star = neicengOptimization2(arg_out.tolist(), myargs) 
  210.  
  211.             # =================================== 
  212.             # 绘制拟合和真是曲线 
  213.             k1 = thisData['K'].tolist() 
  214.  
  215.             # 换成隐含波动率比较 
  216.             tsigma = np.sqrt(omg_i / t2M) 
  217.  
  218.             fig = plt.figure() 
  219.             ax1 = fig.add_subplot(111)  # 在图表1中创建子图1 
  220.  
  221.             ax1.scatter(k1, tsigma, color='.25', label='Value') 
  222.  
  223.             # 生成拟合后的曲线 
  224.             # xaxis = np.linspace(min(k1), max(k1), 100) 
  225.             xaxis = np.linspace(min(k1), max(k1), 100) 
  226.             # np.log(thisData['K'] / (S * np.exp(r * t2M))) 
  227.             # 行权价转化为实际计算用x 
  228.             xbx=[np.log(x / (S * np.exp(r * t2M))) for x in xaxis] 
  229.  
  230.             yaxis = [getfitsiv(xk, t2M, a_star, b_star, rho_star, m_star, sigma_star) for xk in xbx] 
  231.  
  232.             # ax1.plot(xaxis.tolist(), yaxis, color='r--',label='fitValue') 
  233.             ax1.plot(xaxis.tolist(), yaxis, 'r--', label='fitValue') 
  234.  
  235.             # 去除图形顶部边界和右边界的刻度 
  236.             # ax1.tick_params(top='off', right='off') mat 2.1 
  237.             ax1.tick_params(top=False, right=False) 
  238.  
  239.             ax1.set_title('svi模型拟合对比') 
  240.  
  241.             ax1.legend() 
  242.  
  243.             plt.savefig('fitsvi_t2m=' + str(round(t2M, 3)) + '.png') 
  244.  
  245.             plt.close() 
  246.             # plt.show() 
  247.  
  248.         # 检测是否在要求阈值内 
  249.  
  250.         print(cData) 
  251.     except:  # 处理异常 
  252.  
  253.         traceback.print_exc() 
  254.         pass 
  255.  
  256.     print(optionData) 

最后展示一下拟合的结果

This entry was posted in 学术研究 and tagged , , . Bookmark the permalink.

5 Responses to SVI模型拟合

  1. Jojo says:
    February 10, 2020 at 9:49 pm

    请问python对应的lsqlin代码方便分享吗？或者相应的网页链接。搜索lsqlin+python没有对应结果，谢谢

    Reply
  2. L says:
    July 27, 2023 at 11:14 pm

    你好
    文中提到鉴于国内情况，只计算call ,国内股指期权出现 深度实值put 隐含波动率偏低有什么比较深的原因么 ？

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *