量化分析師的Python日記【第4天：一大波金融Library來襲之scipy篇】 · Python 量化交易教程

# 量化分析師的Python日記【第4天：一大波金融Library來襲之scipy篇】 > 來源：https://uqer.io/community/share/54d83bb3f9f06c276f651a6e 上一篇介紹了`numpy`,本篇中著重介紹一下另一個量化金融中常用的庫 `scipy` ## 一、SciPy概述前篇已經大致介紹了NumPy，接下來讓我們看看SciPy能做些什么。NumPy替我們搞定了向量和矩陣的相關操作，基本上算是一個高級的科學計算器。SciPy基于NumPy提供了更為豐富和高級的功能擴展，在統計、優化、插值、數值積分、時頻轉換等方面提供了大量的可用函數，基本覆蓋了基礎科學計算相關的問題。在量化分析中，運用最廣泛的是統計和優化的相關技術，本篇重點介紹SciPy中的統計和優化模塊，其他模塊在隨后系列文章中用到時再做詳述。本篇會涉及到一些矩陣代數，如若感覺不適，可考慮跳過第三部分或者在理解時簡單采用一維的標量代替高維的向量。首先還是導入相關的模塊，我們使用的是SciPy里面的統計和優化部分： ```py import numpy as np import scipy.stats as stats import scipy.optimize as opt ``` ## 二、統計部分 ### 2.1 生成隨機數我們從生成隨機數開始，這樣方便后面的介紹。生成n個隨機數可用`rv_continuous.rvs(size=n)`或`rv_discrete.rvs(size=n)`，其中`rv_continuous`表示連續型的隨機分布，如均勻分布（uniform）、正態分布（norm）、貝塔分布（beta）等；`rv_discrete`表示離散型的隨機分布，如伯努利分布（bernoulli）、幾何分布（geom）、泊松分布（poisson）等。我們生成10個`[0, 1]`區間上的隨機數和10個服從參數`a=4`，`b=2`的貝塔分布隨機數： ```py rv_unif = stats.uniform.rvs(size=10) print rv_unif rv_beta = stats.beta.rvs(size=10, a=4, b=2) print rv_beta [ 0.6419336 0.48403001 0.89548809 0.73837498 0.65744886 0.41845577 0.3823512 0.0985301 0.66785949 0.73163835] [ 0.82164685 0.69563836 0.74207073 0.94348192 0.82979411 0.87013796 0.78412952 0.47508183 0.29296073 0.52551156] ``` 在每個隨機分布的生成函數里，都內置了默認的參數，如均勻分布的上下界默認是0和1。可是一旦需要修改這些參數，每次生成隨機都要敲這么老長一串有點麻煩，能不能簡單點？SciPy里頭有一個Freezing的功能，可以提供簡便版本的命令。`SciPy.stats`支持定義出某個具體的分布的對象，我們可以做如下的定義，讓`beta`直接指代具體參數`a=4`和`b=2`的貝塔分布。為讓結果具有可比性，這里指定了隨機數的生成種子，由NumPy提供。 ```py np.random.seed(seed=2015) rv_beta = stats.beta.rvs(size=10, a=4, b=2) print "method 1:" print rv_beta np.random.seed(seed=2015) beta = stats.beta(a=4, b=2) print "method 2:" print beta.rvs(size=10) method 1: [ 0.43857338 0.9411551 0.75116671 0.92002864 0.62030521 0.56585548 0.41843548 0.5953096 0.88983036 0.94675351] method 2: [ 0.43857338 0.9411551 0.75116671 0.92002864 0.62030521 0.56585548 0.41843548 0.5953096 0.88983036 0.94675351] ``` ### 2.2 假設檢驗好了，現在我們生成一組數據，并查看相關的統計量（相關分布的參數可以在[這里](http://docs.scipy.org/doc/scipy/reference/stats.html)查到）： ```py norm_dist = stats.norm(loc=0.5, scale=2) n = 200 dat = norm_dist.rvs(size=n) print "mean of data is: " + str(np.mean(dat)) print "median of data is: " + str(np.median(dat)) print "standard deviation of data is: " + str(np.std(dat)) mean of data is: 0.383309149888 median of data is: 0.394980561217 standard deviation of data is: 2.00589851641 ``` 假設這個數據是我們獲取到的實際的某些數據，如股票日漲跌幅，我們對數據進行簡單的分析。最簡單的是檢驗這一組數據是否服從假設的分布，如正態分布。這個問題是典型的單樣本假設檢驗問題，最為常見的解決方案是采用K-S檢驗（ Kolmogorov-Smirnov test）。單樣本K-S檢驗的原假設是給定的數據來自和原假設分布相同的分布，在SciPy中提供了`kstest`函數，參數分別是數據、擬檢驗的分布名稱和對應的參數： ```py mu = np.mean(dat) sigma = np.std(dat) stat_val, p_val = stats.kstest(dat, 'norm', (mu, sigma)) print 'KS-statistic D = %6.3f p-value = %6.4f' % (stat_val, p_val) KS-statistic D = 0.037 p-value = 0.9428 ``` 假設檢驗的p-value值很大（在原假設下，p-value是服從`[0, 1]`區間上的均勻分布的隨機變量，可參考 http://en.wikipedia.org/wiki/P-value ），因此我們接受原假設，即該數據通過了正態性的檢驗。在正態性的前提下，我們可進一步檢驗這組數據的均值是不是0。典型的方法是t檢驗（t-test），其中單樣本的t檢驗函數為`ttest_1samp`： ```py stat_val, p_val = stats.ttest_1samp(dat, 0) print 'One-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val) One-sample t-statistic D = 2.696, p-value = 0.0076 ``` 我們看到`p-value<0.05`，即給定顯著性水平0.05的前提下，我們應拒絕原假設：數據的均值為0。我們再生成一組數據，嘗試一下雙樣本的t檢驗（`ttest_ind`）： ```py norm_dist2 = stats.norm(loc=-0.2, scale=1.2) dat2 = norm_dist2.rvs(size=n/2) stat_val, p_val = stats.ttest_ind(dat, dat2, equal_var=False) print 'Two-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val) Two-sample t-statistic D = 3.572, p-value = 0.0004 ``` 注意，這里我們生成的第二組數據樣本大小、方差和第一組均不相等，在運用t檢驗時需要使用Welch's t-test，即指定`ttest_ind`中的`equal_var=False`。我們同樣得到了比較小的`p-value$`，在顯著性水平0.05的前提下拒絕原假設，即認為兩組數據均值不等。 `stats`還提供其他大量的假設檢驗函數，如`bartlett`和`levene`用于檢驗方差是否相等；`anderson_ksam`p用于進行Anderson-Darling的K-樣本檢驗等。 ### 2.3 其他函數有時需要知道某數值在一個分布中的分位，或者給定了一個分布，求某分位上的數值。這可以通過`cdf`和`ppf`函數完成： ```py g_dist = stats.gamma(a=2) print "quantiles of 2, 4 and 5:" print g_dist.cdf([2, 4, 5]) print "Values of 25%, 50% and 90%:" print g_dist.pdf([0.25, 0.5, 0.95]) quantiles of 2, 4 and 5: [ 0.59399415 0.90842181 0.95957232] Values of 25%, 50% and 90%: [ 0.1947002 0.30326533 0.36740397] ``` 對于一個給定的分布，可以用`moment`很方便的查看分布的矩信息，例如我們查看`N(0,1)`的六階原點矩： ```py stats.norm.moment(6, loc=0, scale=1) 15.000000000895332 ``` `describe`函數提供對數據集的統計描述分析，包括數據樣本大小，極值，均值，方差，偏度和峰度： ```py norm_dist = stats.norm(loc=0, scale=1.8) dat = norm_dist.rvs(size=100) info = stats.describe(dat) print "Data size is: " + str(info[0]) print "Minimum value is: " + str(info[1][0]) print "Maximum value is: " + str(info[1][1]) print "Arithmetic mean is: " + str(info[2]) print "Unbiased variance is: " + str(info[3]) print "Biased skewness is: " + str(info[4]) print "Biased kurtosis is: " + str(info[5]) Data size is: 100 Minimum value is: -5.73556523159 Maximum value is: 3.77439818033 Arithmetic mean is: -0.00559348382755 Unbiased variance is: 3.64113204268 Biased skewness is: -0.600615731841 Biased kurtosis is: 0.432147856587 ``` 當我們知道一組數據服從某些分布的時候，可以調用`fit`函數來得到對應分布參數的極大似然估計（MLE, maximum-likelihood estimation）。以下代碼示例了假設數據服從正態分布，用極大似然估計分布參數： ```py norm_dist = stats.norm(loc=0, scale=1.8) dat = norm_dist.rvs(size=100) mu, sigma = stats.norm.fit(dat) print "MLE of data mean:" + str(mu) print "MLE of data standard deviation:" + str(sigma) MLE of data mean:0.00712958665203 MLE of data standard deviation:1.71228079199 ``` `pearsonr`和`spearmanr`可以計算`Pearson`和`Spearman`相關系數，這兩個相關系數度量了兩組數據的相互線性關聯程度： ```py norm_dist = stats.norm() dat1 = norm_dist.rvs(size=100) exp_dist = stats.expon() dat2 = exp_dist.rvs(size=100) cor, pval = stats.pearsonr(dat1, dat2) print "Pearson correlation coefficient: " + str(cor) cor, pval = stats.pearsonr(dat1, dat2) print "Spearman's rank correlation coefficient: " + str(cor) Pearson correlation coefficient: -0.0345336831321 Spearman's rank correlation coefficient: -0.0345336831321 ``` 其中的p-value表示原假設（兩組數據不相關）下，相關系數的顯著性。最后，在分析金融數據中使用頻繁的線性回歸在SciPy中也有提供，我們來看一個例子： ```py x = stats.chi2.rvs(3, size=50) y = 2.5 + 1.2 * x + stats.norm.rvs(size=50, loc=0, scale=1.5) slope, intercept, r_value, p_value, std_err = stats.linregress(x, y) print "Slope of fitted model is:" , slope print "Intercept of fitted model is:", intercept print "R-squared:", r_value**2 Slope of fitted model is: 1.20010505908 Intercept of fitted model is: 2.04778311819 R-squared: 0.781316678233 ``` 在前面的鏈接中，可以查到大部分`stat`中的函數，本節權作簡單介紹，挖掘更多功能的最好方法還是直接讀原始的文檔。另外，StatsModels（ http://statsmodels.sourceforge.net ）模塊提供了更為專業，更多的統計相關函數。若在SciPy沒有滿足需求，可以采用StatsModels。 ## 三、優化部分優化問題在投資中可謂是根本問題，如果手上有眾多可選的策略，應如何從中選擇一個“最好”的策略進行投資呢？這時就需要用到一些優化技術針對給定的指標進行尋優。隨著越來越多金融數據的出現，機器學習逐漸應用在投資領域，在機器學習中，優化也是十分重要的一個部分。以下介紹一些常見的優化方法，雖然例子是人工生成的，不直接應用于實際金融數據，我們希望讀者在后面遇到優化問題時，能夠從這些簡單例子迅速上手解決。 ### 3.1 無約束優化問題所謂的無約束優化問題指的是一個優化問題的尋優可行集合是目標函數自變量的定義域，即沒有外部的限制條件。例如，求解優化問題 ![](https://box.kancloud.cn/2016-07-30_579cb72dc0dbd.jpg) 就是一個無約束優化問題，而求解 ![](https://box.kancloud.cn/2016-07-30_579cb72dd7f00.jpg) 則是一個帶約束的優化問題。更進一步，我們假設考慮的問題全部是凸優化問題，即目標函數是凸函數，其自變量的可行集是凸集。（詳細定義可參考斯坦福大學Stephen Boyd教授的教材convex optimization，下載鏈接：http://stanford.edu/~boyd/cvxbook ）我們以Rosenbrock函數 ![](https://box.kancloud.cn/2016-07-30_579cb72dea8d4.jpg) 作為尋優的目標函數來簡要介紹在SciPy中使用優化模塊`scipy.optimize`。首先需要定義一下這個Rosenbrock函數： ```py def rosen(x): """The Rosenbrock function""" return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) ``` 3.1.1 Nelder-Mead單純形法單純形法是運籌學中介紹的求解線性規劃問題的通用方法，這里的Nelder-Mead單純形法與其并不相同，只是用到單純形的概念。設定起始點`x0=(1.3,0.7,0.8,1.9,1.2)`，并進行最小化的尋優。這里`xtol`表示迭代收斂的容忍誤差上界： ```py x_0 = np.array([0.5, 1.6, 1.1, 0.8, 1.2]) res = opt.minimize(rosen, x_0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True}) print "Result of minimizing Rosenbrock function via Nelder-Mead Simplex algorithm:" print res Optimization terminated successfully. Current function value: 0.000000 Iterations: 436 Function evaluations: 706 Result of minimizing Rosenbrock function via Nelder-Mead Simplex algorithm: status: 0 nfev: 706 success: True fun: 1.6614969876635003e-17 x: array([ 1., 1., 1., 1., 1.]) message: 'Optimization terminated successfully.' nit: 436 ``` Rosenbrock函數的性質比較好，簡單的優化方法就可以處理了，還可以在`minimize`中使用`method='powell'`來指定使用Powell's method。這兩種簡單的方法并不使用函數的梯度，在略微復雜的情形下收斂速度比較慢，下面讓我們來看一下用到函數梯度進行尋優的方法。 3.1.2 Broyden-Fletcher-Goldfarb-Shanno法 Broyden-Fletcher-Goldfarb-Shanno（BFGS）法用到了梯度信息，首先求一下Rosenbrock函數的梯度： ![](https://box.kancloud.cn/2016-07-30_579cb72e09bac.jpg) 其中當`i=j`時，`δi,j=1`，否則`δi,j=0`。邊界的梯度是特例，有如下形式： ![](https://box.kancloud.cn/2016-07-30_579cb72e1c464.jpg) 我們可以如下定義梯度向量的計算函數了： ```py def rosen_der(x): xm = x[1:-1] xm_m1 = x[:-2] xm_p1 = x[2:] der = np.zeros_like(x) der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm) der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0]) der[-1] = 200*(x[-1]-x[-2]**2) return der ``` 梯度信息的引入在`minimize`函數中通過參數`jac`指定： ```py res = opt.minimize(rosen, x_0, method='BFGS', jac=rosen_der, options={'disp': True}) print "Result of minimizing Rosenbrock function via Broyden-Fletcher-Goldfarb-Shanno algorithm:" print res Optimization terminated successfully. Current function value: 0.000000 Iterations: 52 Function evaluations: 63 Gradient evaluations: 63 Result of minimizing Rosenbrock function via Broyden-Fletcher-Goldfarb-Shanno algorithm: status: 0 success: True njev: 63 nfev: 63 hess_inv: array([[ 0.00726515, 0.01195827, 0.0225785 , 0.04460906, 0.08923649], [ 0.01195827, 0.02417936, 0.04591135, 0.09086889, 0.18165604], [ 0.0225785 , 0.04591135, 0.09208689, 0.18237695, 0.36445491], [ 0.04460906, 0.09086889, 0.18237695, 0.36609277, 0.73152922], [ 0.08923649, 0.18165604, 0.36445491, 0.73152922, 1.46680958]]) fun: 3.179561068096293e-14 x: array([ 1. , 0.99999998, 0.99999996, 0.99999992, 0.99999983]) message: 'Optimization terminated successfully.' jac: array([ 4.47207141e-06, 1.30357917e-06, -1.86454207e-07, -2.00564982e-06, 4.98799446e-07]) ``` 3.1.3 牛頓共軛梯度法（Newton-Conjugate-Gradient algorithm）用到梯度的方法還有牛頓法，牛頓法是收斂速度最快的方法，其缺點在于要求Hessian矩陣（二階導數矩陣）。牛頓法大致的思路是采用泰勒展開的二階近似： ![](https://box.kancloud.cn/2016-07-30_579cb72e2daf5.jpg) 其中`H(x0)`表示二階導數矩陣。若Hessian矩陣是正定的，函數的局部最小值可以通過使上面的二次型的一階導數等于0來獲取，我們有： ![](https://box.kancloud.cn/2016-07-30_579cb72e413d5.jpg) 這里可使用共軛梯度近似Hessian矩陣的逆矩陣。下面給出Rosenbrock函數的Hessian矩陣元素通式： ![](https://box.kancloud.cn/2016-07-30_579cb72e56c07.jpg) 其中`i,j∈[1,N?2]`。其他邊界上的元素通式為： ![](https://box.kancloud.cn/2016-07-30_579cb72e68e54.jpg) 例如，當`N=5`時的Hessian矩陣為： ![](https://box.kancloud.cn/2016-07-30_579cb72e7bb11.jpg) 為使用牛頓共軛梯度法，我們需要提供一個計算Hessian矩陣的函數： ```py def rosen_hess(x): x = np.asarray(x) H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1) diagonal = np.zeros_like(x) diagonal[0] = 1200*x[0]**2-400*x[1]+2 diagonal[-1] = 200 diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:] H = H + np.diag(diagonal) return H ``` ```py res = opt.minimize(rosen, x_0, method='Newton-CG', jac=rosen_der, hess=rosen_hess, options={'xtol': 1e-8, 'disp': True}) print "Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian):" print res Optimization terminated successfully. Current function value: 0.000000 Iterations: 20 Function evaluations: 22 Gradient evaluations: 41 Hessian evaluations: 20 Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm: status: 0 success: True njev: 41 nfev: 22 fun: 1.47606641102778e-19 x: array([ 1., 1., 1., 1., 1.]) message: 'Optimization terminated successfully.' nhev: 20 jac: array([ -3.62847530e-11, 2.68148992e-09, 1.16637362e-08, 4.81693414e-08, -2.76999090e-08]) ``` 對于一些大型的優化問題，Hessian矩陣將異常大，牛頓共軛梯度法用到的僅是Hessian矩陣和一個任意向量的乘積，為此，用戶可以提供兩個向量，一個是Hessian矩陣和一個任意向量`p`的乘積，另一個是向量`p`，這就減少了存儲的開銷。記向量`p=(p1,…,pN?1)`，可有 ![](https://box.kancloud.cn/2016-07-30_579cb72e8f479.jpg) 我們定義如下函數并使用牛頓共軛梯度方法尋優： ```py def rosen_hess_p(x, p): x = np.asarray(x) Hp = np.zeros_like(x) Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1] Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \ -400*x[1:-1]*p[2:] Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1] return Hp res = opt.minimize(rosen, x_0, method='Newton-CG', jac=rosen_der, hessp=rosen_hess_p, options={'xtol': 1e-8, 'disp': True}) print "Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian times arbitrary vector):" print res Optimization terminated successfully. Current function value: 0.000000 Iterations: 20 Function evaluations: 22 Gradient evaluations: 41 Hessian evaluations: 58 Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian times arbitrary vector): status: 0 success: True njev: 41 nfev: 22 fun: 1.47606641102778e-19 x: array([ 1., 1., 1., 1., 1.]) message: 'Optimization terminated successfully.' nhev: 58 jac: array([ -3.62847530e-11, 2.68148992e-09, 1.16637362e-08, 4.81693414e-08, -2.76999090e-08]) ``` ### 3.2. 約束優化問題無約束優化問題的一種標準形式為： ![](https://box.kancloud.cn/2016-07-30_579cb72ea3c20.jpg) 其中![](https://box.kancloud.cn/2016-07-30_579cb72eb6da3.jpg)為![](https://box.kancloud.cn/2016-07-30_579cb72ec97c8.jpg)空間上的二次可微的凸函數；`A`為`p×n`矩陣且秩`rankA=p<n`。我們考察如下一個例子： ![](https://box.kancloud.cn/2016-07-30_579cb72edd929.jpg) 定義目標函數及其導數為： ```py def func(x, sign=1.0): """ Objective function """ return sign*(2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2) def func_deriv(x, sign=1.0): """ Derivative of objective function """ dfdx0 = sign*(-2*x[0] + 2*x[1] + 2) dfdx1 = sign*(2*x[0] - 4*x[1]) return np.array([ dfdx0, dfdx1 ]) ``` 其中`sign`表示求解最小或者最大值，我們進一步定義約束條件： ```py cons = ({'type': 'eq', 'fun': lambda x: np.array([x[0]**3 - x[1]]), 'jac': lambda x: np.array([3.0*(x[0]**2.0), -1.0])}, {'type': 'ineq', 'fun': lambda x: np.array([x[1] - 1]), 'jac': lambda x: np.array([0.0, 1.0])}) ``` 最后我們使用SLSQP（Sequential Least SQuares Programming optimization algorithm）方法進行約束問題的求解（作為比較，同時列出了無約束優化的求解）： ```py res = opt.minimize(func, [-1.0, 1.0], args=(-1.0,), jac=func_deriv, method='SLSQP', options={'disp': True}) print "Result of unconstrained optimization:" print res res = opt.minimize(func, [-1.0, 1.0], args=(-1.0,), jac=func_deriv, constraints=cons, method='SLSQP', options={'disp': True}) print "Result of constrained optimization:" print res Optimization terminated successfully. (Exit mode 0) Current function value: -2.0 Iterations: 4 Function evaluations: 5 Gradient evaluations: 4 Result of unconstrained optimization: status: 0 success: True njev: 4 nfev: 5 fun: -1.9999999999999996 x: array([ 2., 1.]) message: 'Optimization terminated successfully.' jac: array([ -2.22044605e-16, -0.00000000e+00, 0.00000000e+00]) nit: 4 Optimization terminated successfully. (Exit mode 0) Current function value: -1.00000018311 Iterations: 9 Function evaluations: 14 Gradient evaluations: 9 Result of constrained optimization: status: 0 success: True njev: 9 nfev: 14 fun: -1.0000001831052137 x: array([ 1.00000009, 1. ]) message: 'Optimization terminated successfully.' jac: array([-1.99999982, 1.99999982, 0. ]) nit: 9 ``` 和統計部分一樣，Python也有專門的優化擴展模塊，CVXOPT（ http://cvxopt.org ）專門用于處理凸優化問題，在約束優化問題上提供了更多的備選方法。CVXOPT是著名的凸優化教材convex optimization的作者之一，加州大學洛杉磯分校Lieven Vandenberghe教授的大作，是處理優化問題的利器。 SciPy中的優化模塊還有一些特殊定制的函數，專門處理能夠轉化為優化求解的一些問題，如方程求根、最小方差擬合等，可到SciPy優化部分的指引頁面查看。