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

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

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

``````from scipy import interpolate
dir(interpolate)[:5]

['Akima1DInterpolator',
'BPoly',
'BarycentricInterpolator',
'BivariateSpline',
'CloughTocher2DInterpolator']
``````

• `xk`离散的自变量值，为序列
• `yk`对应`xk`的函数值，为与`xk`长度相同的序列
• `xnew`需要进行插值的自变量值序列
• `order`样条插值使用的函数基德阶数，为1时使用线性函数
``````print interpolate.spline.__doc__

Interpolate a curve at new points using a spline fit

Parameters
----------
xk, yk : array_like
The x and y values that define the curve.
xnew : array_like
The x values where spline should estimate the y values.
order : int
Default is 3.
kind : string
One of {'smoothest'}
conds : Don't know
Don't know

Returns
-------
spline : ndarray
An array of y values; the spline evaluated at the positions `xnew`.
``````

### 1.1 三角函数（`np.sin`）插值

``````import numpy as np
from matplotlib import pylab
import seaborn as sns
font.set_size(20)
x = np.linspace(1.0, 13.0, 7)
y = np.sin(x)
pylab.figure(figsize = (12,6))
pylab.scatter(x,y, s = 85, marker='x', color = 'r')
pylab.title(u'\$f(x)\$离散点分布', fontproperties = font)

<matplotlib.text.Text at 0x142cafd0>
``````

``````xnew = np.linspace(1.0,13.0,500)
ynewLinear = interpolate.spline(x,y,xnew,order = 1)
ynewLinear[:5]

array([ 0.84147098,  0.83304993,  0.82462888,  0.81620782,  0.80778677])
``````

``````ynewReal = np.sin(xnew)
ynewReal[:5]

array([ 0.84147098,  0.85421967,  0.86647437,  0.87822801,  0.88947378])
``````

``````pylab.figure(figsize = (16,8))
pylab.plot(xnew,ynewReal)
pylab.plot(xnew,ynewLinear)
pylab.plot(xnew,ynewCubicSpline)
pylab.scatter(x,y, s = 160, marker='x', color = 'k')
pylab.legend([u'真实曲线', u'线性插值', u'样条插值', u'\$f(x)\$离散点'], prop = font)
pylab.title(u'\$f(x)\$不同插值方法拟合效果：线性插值 v.s 样条插值', fontproperties = font)

<matplotlib.text.Text at 0x1424cd50>
``````

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

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

2015年3月3日10时波动率矩阵
``````

2.20 32.56 29.75 29.26 27.68
2.30 28.84 29.20 27.39 26.51
2.40 27.66 27.35 25.89 25.28
2.50 26.97 25.57 25.80 25.41
2.60 27.77 24.82 27.34 24.81

### 2.1 方差曲面插值

``````evaluationDate = Date(2015,3,3)
ttm = np.array([(d - evaluationDate) / 365.0 for d in dates])
varianceMatrix = (blackVolMatrix**2) * ttm
varianceMatrix

array([[ 0.00639109,  0.0128489 ,  0.02674114,  0.04324205],
[ 0.0050139 ,  0.01237794,  0.02342277,  0.03966943],
[ 0.00461125,  0.01086231,  0.02093128,  0.03607931],
[ 0.00438413,  0.0094909 ,  0.02079521,  0.03643376],
[ 0.00464918,  0.00894747,  0.02334717,  0.03475378]])
``````

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

``````interp(ttm[0], strikes[0])

array([ 0.00639109])
``````

``````sMeshes = np.linspace(strikes[0], strikes[-1], 400)
tMeshes = np.linspace(ttm[0], ttm[-1], 200)
interpolatedVarianceSurface = np.zeros((len(sMeshes), len(tMeshes)))
for i, s in enumerate(sMeshes):
for j, t in enumerate(tMeshes):
interpolatedVarianceSurface[i][j] = interp(t,s)

interpolatedVolatilitySurface = np.sqrt((interpolatedVarianceSurface / tMeshes))
print u'行权价方向网格数：', np.size(interpolatedVolatilitySurface, 0)
print u'到期时间方向网格数：', np.size(interpolatedVolatilitySurface, 1)

``````

``````pylab.figure(figsize = (16,8))
pylab.plot(sMeshes, interpolatedVolatilitySurface[:, 0])
pylab.scatter(x = strikes, y = blackVolMatrix[:,0], s = 160,marker = 'x', color = 'r')
pylab.legend([u'波动率（线性插值）', u'波动率（离散）'], prop = font)
pylab.title(u'到期时间为2015年3月25日期权波动率', fontproperties = font)

<matplotlib.text.Text at 0xea27f90>
``````

``````from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

maturityMesher, strikeMesher = np.meshgrid(tMeshes, sMeshes)
pylab.figure(figsize = (16,9))
ax = pylab.gca(projection = '3d')
surface = ax.plot_surface(strikeMesher, maturityMesher, interpolatedVolatilitySurface*100, cmap = cm.jet)
pylab.colorbar(surface,shrink=0.75)
pylab.title(u'2015年3月3日10时波动率曲面', fontproperties = font)
pylab.xlabel("strike")
pylab.ylabel("maturity")
ax.set_zlabel(r"volatility(%)")

<matplotlib.text.Text at 0x14e03050>
``````