# 量化分析師的Python日記【第8天 Q Quant兵器譜之函數插值】 > 來源：https://uqer.io/community/share/551cfa1ff9f06c8f339044ff > 在本篇中，我們將介紹Q寬客常用工具之一：函數插值。接著將函數插值應用于一個實際的金融建模場景中：波動率曲面構造。 > 通過本篇的學習您將學習到： > 1. 如何在`scipy`中使用函數插值模塊：`interpolate`； > 2. 波動率曲面構造的原理； > 3. 將`interpolate`運用于波動率曲面構造。 ## 1. 如何使用`scipy`做函數插值 函數插值，即在離散數據的基礎上補插連續函數，估算出函數在其他點處的近似值的方法。在`scipy`中，所有的與函數插值相關的功能都在`scipy.interpolate`模塊中 ```py from scipy import interpolate dir(interpolate)[:5] ['Akima1DInterpolator', 'BPoly', 'BarycentricInterpolator', 'BivariateSpline', 'CloughTocher2DInterpolator'] ``` 作為介紹性質的本篇，我們將只關注`interpolate.spline`的使用，即樣條插值方法： + `xk`離散的自變量值，為序列 + `yk`對應`xk`的函數值，為與`xk`長度相同的序列 + `xnew`需要進行插值的自變量值序列 + `order`樣條插值使用的函數基德階數，為1時使用線性函數 ```py 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`）插值 一例勝千言！讓我們這里用實際的一個示例，來說明如何在`scipy`中使用函數插值。這里的目標函數是三角函數：  假設我們已經觀測到的`f(x)`在離散點`x=(1,3,5,7,9,11,13)`的值： ```py 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> ```  首先我們使用最簡單的線性插值算法，這里面只要將`spline`的參數`order`設置為1即可： ```py 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]) ``` 復雜一些的，也是`spline`函數默認的方法，即為樣條插值，將`order`設置為3即可： 最后我們獲得真實的`sin(x)`的值： ```py ynewReal = np.sin(xnew) ynewReal[:5] array([ 0.84147098, 0.85421967, 0.86647437, 0.87822801, 0.88947378]) ``` 讓我們把所有的函數畫到一起，看一下插值的效果。對于我們這個例子中的目標函數而言，由于本身目標函數是光滑函數，則越高階的樣條插值的方法，插值效果越好。 ```py 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. 函數插值應用 —— 期權波動率曲面構造 市場上期權價格一般以隱含波動率的形式報出，一般來講在市場交易時間，交易員可以看到類似的波動率矩陣（Volatilitie Matrix): ```py 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時波動率矩陣 ``` | 到期時間 | March 25th, 2015 | April 25th, 2015 | June 25th, 2015 | September 25th, 2015 | | --- | --- | | 行權價 | | | | | | 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 | 交易員可以看到市場上離散值的信息，但是如果可以獲得一些隱含的信息更好：例如，在2015年6月25日以及2015年9月25日之間，波動率的形狀會是怎么樣的？ ### 2.1 方差曲面插值 我們并不是直接在波動率上進行插值，而是在方差矩陣上面進行插值。方差和波動率的關系如下：  所以下面我們將通過處理，獲取方差矩陣（Variance Matrix): ```py 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]]) ``` 這里的值`varianceMatrix`就是變換而得的方差矩陣。 下面我們將在行權價方向以及時間方向同時進行線性插值，具體地，行權價方向：  時間方向：  這個過程在`scipy`中可以直接通過`interpolate`模塊下`interp2d`來實現： + `ttm` 時間方向離散點 + `strikes` 行權價方向離散點 + `varianceMatrix` 方差矩陣，列對應時間維度；行對應行權價維度 + `kind = 'linear'` 指示插值以線性方式進行 ```py interp = interpolate.interp2d(ttm, strikes, varianceMatrix, kind = 'linear') ``` 返回的`interp`對象可以用于獲取任意點上插值獲取的方差值： ```py interp(ttm[0], strikes[0]) array([ 0.00639109]) ``` 最后我們獲取整個平面上所有點的方差值，再轉換為波動率曲面。 ```py 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) 行權價方向網格數： 400 到期時間方向網格數： 200 ``` 選取某一個到期時間上的波動率點，看一下插值的效果。這里我們選擇到期時間最近的點：2015年3月25日： ```py 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> ```  最終，我們把整個曲面的圖像畫出來看看： ```py 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> ```