Python學習教程（Python學習路線）：Python——SciPy精講

import scipy

這樣你就可以用 scipy 里面所有的內置方法 (build-in methods) 了，比如插值、積分和優化。

numpy.interpolatenumpy.integratenumpy.optimize

但是每次寫 scipy 字數有點多，通常我們給 scipy 起個別名 sp，用以下語法，這樣所有出現 scipy 的地方都可以用 sp 替代。

import scipy as sp

SciPy 是建立 NumPy 基礎上的，很多關于線性代數的矩陣運算在 NumPy 都能做，因此就不重復在這里講了。此外在〖 數組計算之 NumPy (下) 〗也說過，數組計算比矩陣計算更通用，

本章換一種寫法，我們專門針對科學計算中三個具體問題來介紹 SciPy，它們就是

1. 插值 (interpolation)

2. 積分 (integration)

3. 優化 (optimization)

對于以上每一個知識點我都會介紹一個

• 簡單例子來明晰 SciPy 里各種函數的用法

• 和金融相關的實際例子

• 計算遠期利率：在零息曲線中插值折現因子

• 計算期權價格：將期望寫成積分并數值求解

• 配置資產權重：優化「風險平 價」模型權重

1.插值

給定一組 (xi, yi)，其中 i = 1, 2, ..., n，而且 xi 是有序的，稱為「標準點」。插值就是對于任何新點 xnew，計算出對應的 ynew。換句話來說，插值就是在標準點之間建立分段函數 (piecewise function)，把它們連起來。這樣給定任意連續 x 值，帶入函數就能計算出任意連續 y 值。

在 SciPy 中有個專門的函數 scipy.interpolate 是用來插值的，首先引進它并記為 spi。 

import scipy.interpolate as spi

首先定義 x 和函數 f(x)：

x = np.linspace(-2 * np.pi, 2 * np.pi, 11)f = lambda x: np.sin(x) + 0.5 * xf(x)

array

([

-3.14159265

, 

-1.56221761

, 

-1.29717034

, 

-1.84442231

, 



       -1.57937505

, 

0.

         , 

1.57937505

, 

1.84442231

,



        1.29717034

, 

1.56221761

, 

3.14159265

])

接下來介紹 scipy.interpolate 里面兩大殺器：splrep 和 splev。兩個函數名稱都是以 spl 開頭，全稱 spline (樣條)，可以理解這兩個函數都和樣條有關。不同的是，兩個函數名稱以 rep 和 ev 結尾，它們的意思分別是：

• rep：representation 的縮寫，那么 splrep 其實生成的是一個「表示樣條的對象」

• ev：evaluation 的縮寫，那么 splev 其實用于「在樣條上估值」

splrep 和 splev 像是組合拳 (one two punch)

• 前者將 x, y 和插值方式轉換成「樣條對象」tck

• 后者利用它在 xnew 上生成 ynew

一圖勝千言：

接下來仔細分析一下 tck。

tck = spi.splrep( x, f(x), k=1 )tck

tck 就是樣條對象，以元組形式返回，tck 這名字看起來很奇怪，實際指的是元組 (t, c, k) 里的三元素：

• t - vector of knots (節點)

• c - spline cofficients (系數)

• k - degree of spline (階數)

對照上圖，tck 確實一個元組，包含兩個 array 和一個標量 1，其中

1. 第一個 array 是節點，即標準點，注意到一開始 x 只有 11 個，但現在是 13 個，首尾都往外補了一個和首尾一樣的 x

2. 第二個 array 是系數，注意它就是 y 在尾部補了兩個 0

3. 標量 1 是階數，因為在調用 splrep 時就把 k 設成 1

注：前兩個 array 我只是發現這個規律，但解釋不清楚為什么這樣。它和 matlab 里面的 spline() 的產出不太一樣，希望懂的讀者可以留言區解釋一下。

雖然解釋不清楚前兩個 array，那就把 tck 當成是個黑匣子 (black-box) 直接用了。比如可用 PPoly.from_spline 來查看每個分段函數的系數。

pp = spi.PPoly.from_spline(tck)pp.c.T

array(

[[ 1.25682673, -3.14159265]

,



       [ 1.25682673, -3.14159265]

,



       [ 0.21091791, -1.56221761]

,



       [-0.43548928, -1.29717034]

,



       [ 0.21091791, -1.84442231]

,



       [ 1.25682673, -1.57937505]

,



       [ 1.25682673, 0. ]

,



       [ 0.21091791, 1.57937505]

,



       [-0.43548928, 1.84442231]

,



       [ 0.21091791, 1.29717034]

,



       [ 1.25682673, 1.56221761]

,



       [ 1.25682673, 3.14159265]

])

t ck 的系數數組里有 13 個元素，而上面數組大小是 (12, 2)，12 表示 12 段，2 表示每段線性函數由 2 個系數確定。

把 x 和 tck 丟進 splev 函數，我們可以插出在 x 點對應的值 iy。

iy = spi.splev( x, tck )iy

array

([

-3.14159265

, 

-1.56221761

, 

-1.29717034

, 

-1.84442231

, 



       -1.57937505

, 

0.

         , 

1.57937505

, 

1.84442231

,



        1.29717034

, 

1.56221761

, 

3.14159265

])

用 Matplotlib 來可視化插值的 iy 和原函數 f(x) 發現 iy 都在 f(x) 上。Matplotlib 是之后的課題，現在讀者可忽略其細節。

除了可視化，我們還可以用具體的數值結果來評估插值的效果：

np.allclose(f(x), iy)np.sum((f(x) - iy) ** 2) / len(x)

True




0.0

第一行 allclose 的結果都是 True 證明插值和原函數值完全吻合，第二行就是均方誤差 (mean square error, MSE)，0.0 也說明同樣結果。

上面其實做的是在「標準點 x」上插值，那得到的結果當然等于「標準點 y」了。這種插值確實意義不大，但舉這個例子只想讓大家

1. 明晰 splrep 和 splev 是怎么運作的

2. 如何可視化插出來的值和原函數的值

3. 如何用 allclose 來衡量插值和原函數值之間的差異

一旦弄明白了這些基礎，接下來就可以秒懂更實際的例子了。

正規例子

上面在「標準點 x」上插值有點作弊，現在我們在 50 個「非標準點 xd」上線性插值得到 iyd。

xd = np.linspace( 1.0, 3.0, 50 )iyd = spi.splev( xd, tck )print( xd, '\n\n', iyd )

密密麻麻的數字啥都看不出來，可視化一下把。

這插得是個什么鬼？紅色插值點在第二段和深青色原函數差別也太遠了吧 (MSE 也顯示有差異)。

np.sum((f(xd) - iyd) ** 2) / len(xd)

0.011206417290260647

問題出在哪兒呢？當「標準點 x」不密集時 (只有 11 個)，分段線性函數來擬合 sin(x) + 0.5x 函數當然不會太好啦。那我們試試分段三次樣條函數 (k = 3)。

tck = spi.splrep( x, f(x), k=3 )iyd = spi.splev( xd, tck )

可視化一下并計算 MSE 看看

np.sum((f(xd) - iyd) ** 2) / len(xd)

1.6073247177220542e-05

視覺效果好多了！誤差小多了！

金融例子

用 scipy.interpolate 來插值折現因子來計算遠期利率。

在金融市場中，每個貨幣都有自己相對應的折現曲線，簡單來說，就是在「標準日期」上一組折現因子 (discount factor) 或零息利率 (zero rate)。

那么在「非標準日期」上折現因子或零息利率怎么求呢？插值！

本小節的知識點內容來之〖 弄懂金融十大話題 (上) 〗。

對上面曲線插值有一個概念后，首先用 pandas 讀取數據。Pandas 是下帖內容，你就先把它當成一個可以用字符串來索引或切片的二維數據結構。

import pandas as pdcurve = pd.read_excel('CNY zero curve.xlsx')curve

該曲線用于估值日 2019-04-01，上圖第一個點的日期是 2019-04-03，通常稱為即期日，往后的日期分別是從即期日開始往后推 1W, 1M, 2M, 3M, 6M, 9M, 1Y 和 2Y 得到的。

用 Matplotlib 來可視化折現因子和零息利率。

這里用了雙 y 軸來區分折現因子和零息利率，左邊是折現因子，右邊是零息利率，其實通過觀察 y 軸的數值也可以區分出來兩者。

現在實際問題是要計算從起始日 2019-08-05 到終止日 2019-11-05 的 3M 遠期利率，根據其公式 (不推導)：

要計算遠期利率，核心問題就是計算 2019-08-05 和 2019-11-05 兩天的折現因子。為了簡化，我們把這兩天之間的年限差近似定為 0.25 ≈ 3個月/12個月。具體步驟如下：

1. 計算曲線上「標準日期」到「估值日」之間的天數差

2. 計算「起始日」和「終止日」到「估值日」之間的天數差
   

3. 插出「起始日」和「終止日」上的折現因子 (四種方法)

4. 將折現因子帶入公式計算遠期利率

第一步：計算曲線上「標準日期」到「估值日」之間的天數差

today = pd.Timestamp('2019-04-01')daydiff = curve['Date'] - todaydaydiff

 

2

 days



1

 

9

 days



2

 

32

 days



3

 

63

 days



4

 

93

 days



5

 

185

 days



6

 

277

 days



7

 

368

 days



8

 

733

 days


Name: 

Date

, dtype: timedelta64[ns]

上面結果不是數值型變量 (還帶個 days)，用 .dt.days.values 得到相應的 numpy 數組。

d = daydiff.dt.days.valuesd

array([ 

2

, 

9

, 

32

, 

63

, 

93

, 



      185

, 

277

, 

368

, 

733

], dtype=

int64

)

第二步：計算「起始日」和「終止日」到「估值日」之間的天數差

import datetimestart = datetime.datetime.strptime('2019-08-05','%Y-%m-%d')end = datetime.datetime.strptime('2019-11-05','%Y-%m-%d')d_s = (start - today).daysd_e = (end- today).daysprint( d_s, d_e )

126 218

需要引進 datetime 這個庫將字符型日期轉成真正的 date 格式。

第三步：插出「起始日」和「終止日」上的折現因子，有多種方法，不同數據商對不同曲線也有不同的設置，常見的四種有：

1. 在折現因子上線性插值

2. 在折現因子上三次樣條插值

3. 在 ln(折現因子) 上線性插值

4. 在零息曲線上線性插值，再計算折現因子

DF 上線性插值

tck = spi.splrep( d, curve['Discount Factor'], k=1 )DF_s = spi.splev( d_s, tck )DF_e = spi.splev( d_e, tck )print( DF_s, DF_e )

0.9909485166188177 0.9828538249018102

splrep 里面 k 設為 1 表示線性插值。

DF 上三次樣條插值

tck = spi.splrep( d, curve['Discount Factor'], k=3 )DF_s = spi.splev( d_s, tck )DF_e = spi.splev( d_e, tck )print( DF_s, DF_e )

0.9909572012597055 0.9827493083500931

splrep 里面 k 設為 3 表示三次樣條插值。

ln(DF) 上線性插值

tck = spi.splrep( d, np.log(curve['Discount Factor']), k=1 )DF_s = np.exp(spi.splev( d_s, tck ))DF_e = np.exp(spi.splev( d_e, tck ))print( DF_s, DF_e )

0.9909402218834239 0.9828472942639631

把 ln(DF) 放入 splrep 中，插出來也是 ln 形式，要最終得到折現因子，還需要用 exp 函數還原。

Rate 上線性插值

tck = spi.splrep( d, curve['Zero Rate (%)'], k=1 )r_s = spi.splev( d_s, tck )r_e = spi.splev( d_e, tck )DF_s = np.exp(-d_s/365 * r_s/100)DF_e = np.exp(-d_e/365 * r_e/100)print( DF_s, DF_e )

0.9921606726777862 0.9843810241053533

插出來的零息利率，需要用以下公式計算出折現因子

    DF = exp( -d/365 × r/100)

d 除以 365 轉換成年限，r 除以 100 是因為 r 單位是 %。

第四步：將折現因子帶入公式計算遠期利率

F = 0.25*(DF_s/DF_e - 1) * 100

第三步中四種方法計算出來的遠期利率 (%) 為

1. DF 上線性插值 - 2.059%

2. 折DF 上三次樣條插值 - 2.088%

3. ln(DF) 上線性插值 - 2.059%

4. Rate 上線性插值 - 1.976%

四個遠期利率差別都不大，業界使用較多的是第 3 和 4 種。

2.積分

在 SciPy 中有個專門的函數 scipy.integrate 是用來做數值積分的，首先引進它并記為 sci。 

import scipy.integrate as sci

簡單例子

用 scipy.integrate 來對函數 sin(x) + 0.5x 求積分。

首先定義被積函數 f(x)：

f = lambda x: np.sin(x) + 0.5 * x

假設我們想從 x= 0.5 到 9.5 對 f(x) 求積分，可以手推出

在 scipy.integrate 里還有些數值積分的函數：

• fixed_quad：fixed Gaussian quadrature (定點高斯積分)

• quad：adaptive quadrature (自適應積分)

• romberg：Romberg integration (龍貝格積分)

• trapz：用 trapezoidal 法則

• simps：用 Simpson’s 法則

前三個函數 fixed_quad, quad, romberg 的參數是被積函數、下界和上界。代碼如下：

sci.fixed_quad(f, a, b)[0]sci.quad(f, a, b)[0]sci.romberg(f, a, b)

24

.366995967084602



24

.374754718086752



24

.374754718086713

對后兩個函數 trapz 和 simps，首先在上下界之間取 n 個點 xi，再求出對應的函數值 f(xi)，再把當參數 f(xi) 和 xi 傳到函數中。代碼如下：

xi = np.linspace(a, b, 100)sci.trapz( f(xi), xi )sci.simps( f(xi), xi )

24

.373463386819818



24

.37474011548407

和解析解 24.37475471808675 比較，quad 的結果最精確。一般當被積函數不規則時 (某段函數值激增)，quad (自適應積分) 的結果也是最好。

金融例子

用 scipy.integrate 來以數值積分的形式給歐式期權定價。

注：本節主要將數值積分的用途，因此金融模型上的很多設置我們都用最簡單的，比如常數型的模型參數等等。

股票類的Black-Scholes (BS) 模型下的 SDE 是描述股票價格 (stock price) 的走勢：

其中

    S(t) = 股票價格

    r = 瞬時無風險利率

    σ = S(t)的瞬時波動率

    B(t) = 布朗運動

歐式看漲期權 (call option) 在 BS 模型下的解析解 (closed-form solution) 如下：

編寫一個程序計算 call 的解析解很容易：

這里需要引入 scipy.stats 下的 norm 庫，使用里面 cdf 函數來計算正態分布的累積分布概率。

假設股價 S0 = 100，行權價格 K = 95，利率為 5%，期限為 1 年，波動率為 10%，帶入寫好的 bscall 函數來計算期權的價值。

(S0, K, r, T, sigma) = (100, 95, 0.05, 1, 0.1)bscall( S0, K, r, T, sigma )

10.405284289598569

大概記注上面的期權值 10.4053。假設我們推導能力不強或者對于更復雜的期權沒有解析解，只要知道 ST 的分布，我們可以試著把「期望值」寫成「積分」形式，再用 x = lnST 做個轉換，最終可推出下式：

為了求數值積分，我們需要知道 x 是如何分布，也就是推出 x 的密度分布函數 fX(x)，推導如下 (不是本帖的重點，如無興趣可跳過下框內容)：

顯然 x 是個正態分布，均值為 lnS0 +(rT - 0.5σ2T)，方差是 σ2T。用 NPDF 代表正態分布 (Normal) 的密度分布函數 (PDF)，可把 call 的價值寫成積分形式，其中

• 被積函數是「支付函數」和「正態分布密度分布函數 」的乘積

• 下界和上界分別是 lnK 和 +∞

最終表達式如下：

跟著「被積函數」的表達式敲代碼

mu = np.log(S0) +(r*T-0.5*sigma**2*T)v = sigma*np.sqrt(T)f = lambda x: np.exp(-r*T) * (np.exp(x)-K) * norm.pdf(x,mu,v)

定義上界和下界

(lb, ub) = (np.log(K), 7)

注意上界不要定義成 +∞。稍微分析下 x = lnST，當 ST= e7 ≈ 1097 對于 S0 = 100 已經很大了，因此上界設為 7 比較合理。

看看三個數值積分的結果如何。

sci.quad(f, lb, ub)[0]xi = np.linspace(lb, ub, 1000)sci.trapz( f(xi), xi )sci.simps( f(xi), xi )

10

.405284289598615



10

.405170993379011



10

.405287100064612

結果和之前的解析解 10.4053 都相當接近。

用數值積分來求解歐式期權的確有點多此一舉 (ovekill)，但很多復雜的產品是沒有解析解的，除了用數值解的「偏微分方法有限差分法」和「蒙特卡洛法」，數值積分也是一種選擇。

3.優化

在 SciPy 中有個專門的函數 scipy.optimize 是用來優化的，首先引進它并記為 spo。  

import scipy.optimize as spo

優化問題可分為無約束優化 (unconstrained optimization) 和有約束優化 (constrained optimization)，我們用簡單例子來介紹前者，用金融例子來介紹后者。

首先定義函數

f = lambda x,y: np.sin(x) + 0.05 * x**2 + np.sin(y) + 0.05 * y**2

接著可視化函數

不難看出該函數有多個局部最小值 (local minimum) 和一個全局最小值 (global minimum)。我們目標是求后者，主要步驟如下：

1. 在 (x-y) 定義域上選點，求出函數值 f(x, y)，找出最小值對應的 x* 和 y*

2. 用 x* 和 y* 當初始值，求出函數全局最小值

第一步：用蠻力找函數最小值以及對應的參數

之所以用「蠻力」一詞，是因為接下來要用到 brute 函數，而 bruteforce 就是蠻力的意思。首先定義函數 fo (其實就是上面的 f)，只不過 brute 函數要求用一個元組把若干參數集合起來。此外我們添加一個 print() 語句，為了檢查中間產出。

將 x 和 y 在 -10 到 10 以步長為 5 來切片 (回顧切片是前閉后開的，因此切片完得到的是 -10, -5, 0, 5，而不包括 10 這個點)

output = Truerranges = ((-10, 10, 5), (-10, 10, 5))spo.brute(fo, rranges, finish=None)

從上面結果可看出，函數在 (0, 0) 是取最小值 0。真是最小值嗎？我也不知道，但是以 5 為步長是不是太粗糙了些，接下來用 0.1 為步長。這時把 output 設為 False 是因為不想看到打印的內容。

output = Falserranges = ((-10, 10, 0.1), (-10, 10, 0.1))opt1 = spo.brute(fo, rranges, finish=None)opt1

array([-1.4, -1.4])

fo(opt1)

-1.7748994599769203

當步長變小，我們能在更細的網格上計算函數值，這是函數在 (-1.4, -1.4) 取最小值 -1.7749，明顯比函數在 (0, 0) 上的值 0 要小。

第二步：把參數當初始值，求函數全局最小值

如果網格足夠密，上面 -1.7749 大概率是全局最小值而  (-1.4, -1.4) 是對應的最優解；如果網格不是足夠密，那么以 (-1.4, -1.4) 當初始值，也能很大概率找到全局最小值。

用 fmin 函數，將剛才 opt1 傳進去，并設定 x 和 f 的終止條件 xtol 和 ftol，和最多迭代次數 maxiter 和最多運行函數次數 maxfun。

output = Trueopt2 = spo.fmin( fo, opt1, xtol=0.001, ftol=0.001,                  maxiter=15, maxfun=20 )opt2

此時最優解為 (-1.42702972, -1.42876755)，而對應的函數值為

output = Falsefo(opt2)

-1.7757246992239009

比剛才函數在 (-1.4, -1.4) 取的最小值 -1.7749 又小了一些。

好的初始值對求函數的最優解影響比較大。假設我們無腦的用 (2, 2) 當初始值，看看會發生什么。

output = Falseopt3 = spo.fmin(fo, (2, 2), maxiter=250)opt3

Optimization terminated successfully.


Current 


function

 

value

: 0.015826



Iteration

s:

 

46



Function evaluation

s:

 

86





array([

4.2710728

 , 

4.27106945

])

求得函數在 (4.2710728, 4.2710728) 取的最小值 0.015826，是不是錯的太離譜了。

金融例子

本小節的知識點內容來自〖 資產配置 〗。

投資組合的資產配置是個很重要的課題，投資者為了最大化回報或最小化風險，可以給各種資產配置不同的權重。本節我們看一個很流行的資產配置方法 - 風險平價 (Risk Parity, RP)。但首先我們先來看看它的通用版本，風險預算 (Risk Budgeting, RB)。

風險平價 (RP) 就是等量的風險預算，即為投資組合中的所有資產分配相等的風險。

這是一個有約束 (constrainted) 的優化問題，我們可用 scipy.optimize 里的 minimize 函數來求解 RP 的權重。首先來定義 risk_parity 函數：

該函數的兩個參數 sigma 和 rho 是 n 個資產的波動率向量 (一維數組) 和相關系數矩陣 (二維數組)，其中

• obj 就是用 numpy 把上面目標函數用「匿名函數」的形式表示出來

• 限制條件有兩種形式，等式 (eq) 和不等式 (ineq)，分別用 dict 的形式表達，而限制條件的表達式也用「匿名函數」來表示

最后在 minimize 函數設定好目標函數、初始值、算法、限制條件和終止條件，得到一個 dict 類的結果 w。

兩個資產

先分析簡單的股票和債券兩個資產組合：

• 股票的預期超額回報為 10%，波動率為 20%

• 債券的預期超額回報為 5%，波動率為 10%

• 它們相關系數為 -10%

mu = np.array([0.1, 0.05])sigma = np.array([0.2, 0.1])rho = np.array([[1, -0.1], [-0.1, 1]])

運行 risk_partiy 函數

result = risk_parity( sigma, rho )result

    fun

: 

3.26901989274624

e-

15




    jac

: array([-

1.86742928

e-

07

, 

1.55459627

e-

07

])



message

: 

'Optimization terminated successfully.'




   nfev

: 

22




    nit

: 

5




   njev

: 

5




 status

: 





success

: True



      x

: array([

0.33333332

, 

0.66666668

])

result 是一個字典：

• 'fun' 對應的是目標函數在最優解下的值，非常小接近于零證明找到了最優值。

• 'nfev' 對應的 22 表示目標函數被運行了 22 次

如果只關注最優權重，只需看 ‘x’

result.x

array([0.33333332, 0.66666668])

股票和債券的最優權重為 w* =[33.33%, 66.66%]

三個資產

接著分析股票、債券和信貸三個資產組合：

• 股票的預期超額回報為 10%，波動率為 20%

• 債券的預期超額回報為 5%，波動率為 10%

• 信貸的預期超額回報為 10%，波動率為 15%

• 股票與債券、股票與信貸、債券與信貸的相關系數為 -10%, 30%, -30%

mu = np.array([0.1, 0.05, 0.1])sigma = np.array([0.2, 0.1, 0.15])rho = np.array([[1, -0.1, 0.3], [-0.1, 1, -0.3], [0.3, -0.3, 1]])

運行 risk_partiy 函數

result = risk_parity( sigma, rho )result.x

array([0.19117648, 0.5147059 , 0.29411762])

股票、債券和信貸的最優權重為 w* = [19.12%, 51.47%, 29.41%]

4總結

• 用 scipy.interpolate 來插值 (interpolation)

• 用 scipy.integrate 來積分 (integration)

• 用 scipy.optimize 來優化 (optimization)

學完此貼后，至少你可以解決以下具體金融問題 (training set  )：

1. 在折現曲線上插出折現因子和零息利率

2. 數值積分求解期權價值

3. 優化出風險平價模型的權重

舉一反三一下，你還可以解決新的金融問題 (test set  )：

1. 在波動平面上插出波動率

2. 數值積分求解而二維金融衍生品價值

3. 優化出各種資產配置模型的權重 (加各種約束)