Python學習教程（Python學習路線）：Python——SciPy精講
Python學習教程（Python學習路線）：Python——SciPy精講
SciPy 是 Python 裡處理科學計算 (scientific computing) 的包，使用它遇到問題可訪問它的官網 (https://www.scipy.org/). 去找答案。 在使用 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，它們就是
- 插值 (interpolation)
- 積分 (integration)
- 優化 (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
簡單例 子
用 scipy.interpolate 來插值函數 sin(x) + 0.5x。
基本概念
首先定義 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，其中
- 第一個 array 是節點，即標準點，注意到一開始 x 只有 11 個，但現在是 13 個，首尾都往外補了一個和首尾一樣的 x
- 第二個 array 是係數，注意它就是 y 在尾部補了兩個 0
- 標量 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」了。這種插值確實意義不大，但舉這個例子只想讓大家
- 明晰 splrep 和 splev 是怎麼運作的
- 如何可視化插出來的值和原函數的值
- 如何用 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)。
那麼在「非標準日期」上折現因子或零息利率怎麼求呢？插值！
本小節的知識點內容來之〖 弄懂金融十大話題 (上) 〗。
知識點
這裡面說的插值是分段 (piecewise) 插值！對於線性插值，不是說一條直線擬合上表的 9 個點，這樣也是不可能做到的。但是分段線性插值就可以完美解決這個問題，因為 9 個點，有 8 段，每一段首尾兩個點，可以連一條直線，全部點之間連起來不就是分段線性插值嗎？三種最常見的插值方法
- 分段常函數
- 分段線性函數
- 分段三次樣條函數
首先給出數學符號。給定 N 數據點 (xi, fi), i = 1, 2, …, N，其中 x1 < x2 < ... < xN 。我們希望找到一個函數 f(x) 來擬合這 N 個數據點，對於分段函數，因為有 N 個數據點，需要 N -1 段函數。
分段常 (piecewise constant) 函數
在這種情況，每一段函數都是一個常數，這種插值方法
- 優點是簡單
- 缺點是在數據點上不連續，更不可導
- 適用於在某些模型的參數 (比如 Heston 模型中的均值迴歸率和波動率的波動率) 上插值 (模型參數通常只用常數和分段常函數，但後者比前者能更好的擬合市場數據，因為它有更多自由度)。
- 不適用於曲線和波動率插值
分段常函數不連續，通常稱作 C-1 函數。
分段線性 (piecewise linear) 函數
在這種情況，每一段函數都是一個線性函數，這種插值方法
- 優點是簡單，在數據點上連續，而且形狀保持性很好 (插出的值只和它相鄰兩個數據點有關，別的數據怎麼動都不影響它的插值)
- 缺點是在數據點上不可導
- 適用於曲線和波動率插值
- 不適用於在 Hull-White 模型下的曲線插值 (Hull-White 模型需要對曲線求二階導)
分段線性函數連續但是不可導，通常稱作 C0 函數。
分段三次樣條 (piecewise cubic spline) 函數
在這種情況，每一段函數都是一個三次多項式函數，這種插值方法
- 優點是在數據點上可導甚至可導三次 (非常平滑)
- 缺點是有些複雜，而且形狀保持性不好 (插出的值和整個數據點有關，別的數據動以下都會影響它的插值)
- 適用於曲線的插值
分段三次樣條函數連續而且二階可導，通常稱作 C2 函數。
對上面曲線插值有一個概念後，首先用 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個月。具體步驟如下：
- 計算曲線上「標準日期」到「估值日」之間的天數差
- 計算「起始日」和「終止日」到「估值日」之間的天數差
- 插出「起始日」和「終止日」上的折現因子 (四種方法)
- 將折現因子帶入公式計算遠期利率
第一步：計算曲線上「標準日期」到「估值日」之間的天數差
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 格式。
第三步：插出「起始日」和「終止日」上的折現因子，有多種方法，不同數據商對不同曲線也有不同的設置，常見的四種有：
- 在折現因子上線性插值
- 在折現因子上三次樣條插值
- 在 ln(折現因子) 上線性插值
- 在零息曲線上線性插值，再計算折現因子
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
第三步中四種方法計算出來的遠期利率 (%) 為
- DF 上線性插值 - 2.059%
- 折DF 上三次樣條插值 - 2.088%
- ln(DF) 上線性插值 - 2.059%
- 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)，推導如下 (不是本帖的重點，如無興趣可跳過下框內容)：
給定 S 的隨機微分方程，首先用伊藤公式推出 lnS 的隨機微分方程
在 0 到 T 兩邊求積分，整理得到 S T 的解。
其中 z 是標準正態分佈變量 z ~ N(0, 1)。
用之前的變量轉化 x = lnS T 得到 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)，我們用簡單例子來介紹前者，用金融例子來介紹後者。
簡單例子
用 scipy.optimize 來求出函數 sin(x) + 0.05x2 + sin(y) + 0.05y2 的最小值。
首先定義函數
f = lambda x,y: np.sin(x) + 0.05 * x**2 + np.sin(y) + 0.05 * y**2
接著可視化函數
不難看出該函數有多個局部最小值 (local minimum) 和一個全局最小值 (global minimum)。我們目標是求後者，主要步驟如下：
- 在 (x-y) 定義域上選點，求出函數值 f(x, y)，找出最小值對應的 x* 和 y*
- 用 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，是不是錯的太離譜了。
金融例子
用 scipy.optimize 來用「風險平價」模型為資產配置最優權重。
本小節的知識點內容來自〖 資產配置 〗。
投資組合的資產配置是個很重要的課題，投資者為了最大化回報或最小化風險，可以給各種資產配置不同的權重。本節我們看一個很流行的資產配置方法 - 風險平價 (Risk Parity, RP)。但首先我們先來看看它的通用版本，風險預算 (Risk Budgeting, RB)。
知識點
風險預算 (RB) 可以基於投資者對資產未來表現 (主要是風險) 的具體看法，或一些通用原則來給資產來分配風險預算，而不是給資產分配權重。下圖畫出兩者的區別。
傳統的 FW 模型把 60% 分給股票而 40% 分給債券，但是這樣的一個投資組合 90% 的風險都來自股票只有 10% 的風險來自債券。那麼這個組合更容易出現股票尾部風險 (tail risk)。一個風險更均衡的投資組合應該選擇配置更多債券 (比如 75%) 和更少股票 (比如 25%)，如下圖所示。
RB 模型的思路就是通過分配風險 (上圖的風險比例) 來影響權重 (上圖的資產權重)，通常是給風險低的資產 (如債券) 高風險配額，而風險高的資產 (如股票) 低風險配額。
接下來我們看看 RB 模型的數學公式吧，首先回顧組合風險
對於第 i 個資產，其邊際風險貢獻 (Marginal Risk Contribution, MRC) 是該資產權重 wi 的微小變化對組合風險 σp 所帶來的影響。用數學公式表示就是對 σp 求 wi 的偏導數。
第 i 個資產的總體風險貢獻 (Total Risk Contribution, TRC) 是其 MRC 乘以其權重，顧名思義，這個總體貢獻一方面來自 MRC，一方面來自權重，數學表達式為：
根據 TRCi 的定義，即第 i 個資產對總體風險的貢獻，可推出它們總和應該等於組合風險 sp，從數學上也可證實此關係
上式兩邊同時除以 σp，並定義風險預算 si 為 TRCi 的佔比，可得 sT1 = 1
由上式看出 si 也類似於權重，只不過是風險上的權重，而 wi 是資產上的權重。下面給出 si 和 wi 之間的關係
在 RB 模型中，股票權重等於風險預算除以貝塔，因此，RB 模型依賴於貝塔的預測質量。歸一化之後的權重等於
事先將一組風險預算分配好，例如 s = [20%, 30%, 50%]，再數值解下面序列二次規劃 (Sequential Quadratic Programming, SQP) 問題可以得到 RB 模型下的最佳權重
風險平價 (RP) 就是等量的風險預算，即為投資組合中的所有資產分配相等的風險。
知識點
類比 RB，RP 給每個風險配額 si 的分配 1/n 的權重，這時組合權重為
同樣可得到 RP 模型下的優化問題 (用 1/n 替代 si )
這是一個有約束 (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 可以實現的三個應用
- 用 scipy.interpolate 來插值 (interpolation)
- 用 scipy.integrate 來積分 (integration)
- 用 scipy.optimize 來優化 (optimization)
學完此貼後，至少你可以解決以下具體金融問題 (training set
)：
- 在折現曲線上插出折現因子和零息利率
- 數值積分求解期權價值
- 優化出風險平價模型的權重
舉一反三一下，你還可以解決新的金融問題
- 在波動平面上插出波動率
- 數值積分求解而二維金融衍生品價值
- 優化出各種資產配置模型的權重 (加各種約束)