量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】

来源:https://uqer.io/community/share/551cfa1ff9f06c8f339044ff

在本篇中,我们将介绍Q宽客常用工具之一:函数插值。接着将函数插值应用于一个实际的金融建模场景中:波动率曲面构造。

通过本篇的学习您将学习到:

  1. 如何在scipy中使用函数插值模块:interpolate
  2. 波动率曲面构造的原理;
  3. interpolate运用于波动率曲面构造。

1. 如何使用scipy做函数插值

函数插值,即在离散数据的基础上补插连续函数,估算出函数在其他点处的近似值的方法。在scipy中,所有的与函数插值相关的功能都在scipy.interpolate模块中

  1. from scipy import interpolate
  2. dir(interpolate)[:5]
  3. ['Akima1DInterpolator',
  4. 'BPoly',
  5. 'BarycentricInterpolator',
  6. 'BivariateSpline',
  7. 'CloughTocher2DInterpolator']

作为介绍性质的本篇,我们将只关注interpolate.spline的使用,即样条插值方法:

  • xk离散的自变量值,为序列
  • yk对应xk的函数值,为与xk长度相同的序列
  • xnew需要进行插值的自变量值序列
  • order样条插值使用的函数基德阶数,为1时使用线性函数
  1. print interpolate.spline.__doc__
  2. Interpolate a curve at new points using a spline fit
  3. Parameters
  4. ----------
  5. xk, yk : array_like
  6. The x and y values that define the curve.
  7. xnew : array_like
  8. The x values where spline should estimate the y values.
  9. order : int
  10. Default is 3.
  11. kind : string
  12. One of {'smoothest'}
  13. conds : Don't know
  14. Don't know
  15. Returns
  16. -------
  17. spline : ndarray
  18. An array of y values; the spline evaluated at the positions `xnew`.

1.1 三角函数(np.sin)插值

一例胜千言!让我们这里用实际的一个示例,来说明如何在scipy中使用函数插值。这里的目标函数是三角函数:

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图1

假设我们已经观测到的f(x)在离散点x=(1,3,5,7,9,11,13)的值:

  1. import numpy as np
  2. from matplotlib import pylab
  3. import seaborn as sns
  4. font.set_size(20)
  5. x = np.linspace(1.0, 13.0, 7)
  6. y = np.sin(x)
  7. pylab.figure(figsize = (12,6))
  8. pylab.scatter(x,y, s = 85, marker='x', color = 'r')
  9. pylab.title(u'$f(x)$离散点分布', fontproperties = font)
  10. <matplotlib.text.Text at 0x142cafd0>

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图2

首先我们使用最简单的线性插值算法,这里面只要将spline的参数order设置为1即可:

  1. xnew = np.linspace(1.0,13.0,500)
  2. ynewLinear = interpolate.spline(x,y,xnew,order = 1)
  3. ynewLinear[:5]
  4. array([ 0.84147098, 0.83304993, 0.82462888, 0.81620782, 0.80778677])

复杂一些的,也是spline函数默认的方法,即为样条插值,将order设置为3即可:

最后我们获得真实的sin(x)的值:

  1. ynewReal = np.sin(xnew)
  2. ynewReal[:5]
  3. array([ 0.84147098, 0.85421967, 0.86647437, 0.87822801, 0.88947378])

让我们把所有的函数画到一起,看一下插值的效果。对于我们这个例子中的目标函数而言,由于本身目标函数是光滑函数,则越高阶的样条插值的方法,插值效果越好。

  1. pylab.figure(figsize = (16,8))
  2. pylab.plot(xnew,ynewReal)
  3. pylab.plot(xnew,ynewLinear)
  4. pylab.plot(xnew,ynewCubicSpline)
  5. pylab.scatter(x,y, s = 160, marker='x', color = 'k')
  6. pylab.legend([u'真实曲线', u'线性插值', u'样条插值', u'$f(x)$离散点'], prop = font)
  7. pylab.title(u'$f(x)$不同插值方法拟合效果:线性插值 v.s 样条插值', fontproperties = font)
  8. <matplotlib.text.Text at 0x1424cd50>

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图3

2. 函数插值应用 —— 期权波动率曲面构造

市场上期权价格一般以隐含波动率的形式报出,一般来讲在市场交易时间,交易员可以看到类似的波动率矩阵(Volatilitie Matrix):

  1. import pandas as pd
  2. pd.options.display.float_format = '{:,>.2f}'.format
  3. dates = [Date(2015,3,25), Date(2015,4,25), Date(2015,6,25), Date(2015,9,25)]
  4. strikes = [2.2, 2.3, 2.4, 2.5, 2.6]
  5. blackVolMatrix = np.array([[ 0.32562851, 0.29746885, 0.29260648, 0.27679993],
  6. [ 0.28841840, 0.29196629, 0.27385023, 0.26511898],
  7. [ 0.27659511, 0.27350773, 0.25887604, 0.25283775],
  8. [ 0.26969754, 0.25565971, 0.25803327, 0.25407669],
  9. [ 0.27773032, 0.24823248, 0.27340796, 0.24814975]])
  10. table = pd.DataFrame(blackVolMatrix * 100, index = strikes, columns = dates, )
  11. table.index.name = u'行权价'
  12. table.columns.name = u'到期时间'
  13. print u'2015年3月3日10时波动率矩阵'
  14. table
  15. 20153310时波动率矩阵
到期时间March 25th, 2015April 25th, 2015June 25th, 2015September 25th, 2015
行权价
2.2032.5629.7529.2627.68
2.3028.8429.2027.3926.51
2.4027.6627.3525.8925.28
2.5026.9725.5725.8025.41
2.6027.7724.8227.3424.81

交易员可以看到市场上离散值的信息,但是如果可以获得一些隐含的信息更好:例如,在2015年6月25日以及2015年9月25日之间,波动率的形状会是怎么样的?

2.1 方差曲面插值

我们并不是直接在波动率上进行插值,而是在方差矩阵上面进行插值。方差和波动率的关系如下:

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图4

所以下面我们将通过处理,获取方差矩阵(Variance Matrix):

  1. evaluationDate = Date(2015,3,3)
  2. ttm = np.array([(d - evaluationDate) / 365.0 for d in dates])
  3. varianceMatrix = (blackVolMatrix**2) * ttm
  4. varianceMatrix
  5. array([[ 0.00639109, 0.0128489 , 0.02674114, 0.04324205],
  6. [ 0.0050139 , 0.01237794, 0.02342277, 0.03966943],
  7. [ 0.00461125, 0.01086231, 0.02093128, 0.03607931],
  8. [ 0.00438413, 0.0094909 , 0.02079521, 0.03643376],
  9. [ 0.00464918, 0.00894747, 0.02334717, 0.03475378]])

这里的值varianceMatrix就是变换而得的方差矩阵。

下面我们将在行权价方向以及时间方向同时进行线性插值,具体地,行权价方向:

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图5

时间方向:

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图6

这个过程在scipy中可以直接通过interpolate模块下interp2d来实现:

  • ttm 时间方向离散点
  • strikes 行权价方向离散点
  • varianceMatrix 方差矩阵,列对应时间维度;行对应行权价维度
  • kind = 'linear' 指示插值以线性方式进行
  1. interp = interpolate.interp2d(ttm, strikes, varianceMatrix, kind = 'linear')

返回的interp对象可以用于获取任意点上插值获取的方差值:

  1. interp(ttm[0], strikes[0])
  2. array([ 0.00639109])

最后我们获取整个平面上所有点的方差值,再转换为波动率曲面。

  1. sMeshes = np.linspace(strikes[0], strikes[-1], 400)
  2. tMeshes = np.linspace(ttm[0], ttm[-1], 200)
  3. interpolatedVarianceSurface = np.zeros((len(sMeshes), len(tMeshes)))
  4. for i, s in enumerate(sMeshes):
  5. for j, t in enumerate(tMeshes):
  6. interpolatedVarianceSurface[i][j] = interp(t,s)
  7. interpolatedVolatilitySurface = np.sqrt((interpolatedVarianceSurface / tMeshes))
  8. print u'行权价方向网格数:', np.size(interpolatedVolatilitySurface, 0)
  9. print u'到期时间方向网格数:', np.size(interpolatedVolatilitySurface, 1)
  10. 行权价方向网格数: 400
  11. 到期时间方向网格数: 200

选取某一个到期时间上的波动率点,看一下插值的效果。这里我们选择到期时间最近的点:2015年3月25日:

  1. pylab.figure(figsize = (16,8))
  2. pylab.plot(sMeshes, interpolatedVolatilitySurface[:, 0])
  3. pylab.scatter(x = strikes, y = blackVolMatrix[:,0], s = 160,marker = 'x', color = 'r')
  4. pylab.legend([u'波动率(线性插值)', u'波动率(离散)'], prop = font)
  5. pylab.title(u'到期时间为2015年3月25日期权波动率', fontproperties = font)
  6. <matplotlib.text.Text at 0xea27f90>

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图7

最终,我们把整个曲面的图像画出来看看:

  1. from mpl_toolkits.mplot3d import Axes3D
  2. from matplotlib import cm
  3. maturityMesher, strikeMesher = np.meshgrid(tMeshes, sMeshes)
  4. pylab.figure(figsize = (16,9))
  5. ax = pylab.gca(projection = '3d')
  6. surface = ax.plot_surface(strikeMesher, maturityMesher, interpolatedVolatilitySurface*100, cmap = cm.jet)
  7. pylab.colorbar(surface,shrink=0.75)
  8. pylab.title(u'2015年3月3日10时波动率曲面', fontproperties = font)
  9. pylab.xlabel("strike")
  10. pylab.ylabel("maturity")
  11. ax.set_zlabel(r"volatility(%)")
  12. <matplotlib.text.Text at 0x14e03050>

量化分析师的Python日记【第8天 Q Quant兵器谱之函数插值】 - 图8