一、SciPy概述
前篇已經(jīng)大致介紹了NumPy，接下來讓我們看看SciPy能做些什么。NumPy替我們搞定了向量和矩陣的相關(guān)操作，基本上算是一個高級的科學計算器。SciPy基于NumPy提供了更為豐富和高級的功能擴展，在統(tǒng)計、優(yōu)化、插值、數(shù)值積分、時頻轉(zhuǎn)換等方面提供了大量的可用函數(shù)，基本覆蓋了基礎科學計算相關(guān)的問題。

在量化分析中，運用最廣泛的是統(tǒng)計和優(yōu)化的相關(guān)技術(shù)，本篇重點介紹SciPy中的統(tǒng)計和優(yōu)化模塊，其他模塊在隨后系列文章中用到時再做詳述。
本篇會涉及到一些矩陣代數(shù)，如若感覺不適，可考慮跳過第三部分或者在理解時簡單采用一維的標量代替高維的向量。
首先還是導入相關(guān)的模塊，我們使用的是SciPy里面的統(tǒng)計和優(yōu)化部分：
In [1]:

import numpy as np
import scipy.stats as stats
import scipy.optimize as opt
二、統(tǒng)計部分?
2.1 生成隨機數(shù)?
我們從生成隨機數(shù)開始，這樣方便后面的介紹。生成n個隨機數(shù)可用rv_continuous.rvs(size=n)或rv_discrete.rvs(size=n)，其中rv_continuous表示連續(xù)型的隨機分布，如均勻分布（uniform）、正態(tài)分布（norm）、貝塔分布（beta）等；rv_discrete表示離散型的隨機分布，如伯努利分布（bernoulli）、幾何分布（geom）、泊松分布（poisson）等。我們生成10個[0, 1]區(qū)間上的隨機數(shù)和10個服從參數(shù)$a = 4$，$b = 2$的貝塔分布隨機數(shù)：
In [2]:

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.20630272  0.25929204  0.16859206  0.92573462  0.16383319  0.3475617
  0.83792048  0.79574153  0.37945051  0.23439682]
[ 0.71216492  0.85688464  0.70310131  0.3783662   0.69507561  0.78626586
  0.54529967  0.4261079   0.26646767  0.8519046 ]
在每個隨機分布的生成函數(shù)里，都內(nèi)置了默認的參數(shù)，如均勻分布的上下界默認是0和1。可是一旦需要修改這些參數(shù)，每次生成隨機都要敲這么老長一串有點麻煩，能不能簡單點？SciPy里頭有一個Freezing的功能，可以提供簡便版本的命令。SciPy.stats支持定義出某個具體的分布的對象，我們可以做如下的定義，讓beta直接指代具體參數(shù)$a = 4$和$b = 2$的貝塔分布。為讓結(jié)果具有可比性，這里指定了隨機數(shù)的生成種子，由NumPy提供。
In [3]:

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 假設檢驗?
好了，現(xiàn)在我們生成一組數(shù)據(jù)，并查看相關(guān)的統(tǒng)計量（相關(guān)分布的參數(shù)可以在這里查到：http://docs.scipy.org/doc/scipy/reference/stats.html）：
In [4]:

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.705195138069
median of data is: 0.658167882933
standard deviation of data is: 2.08967006905
假設這個數(shù)據(jù)是我們獲取到的實際的某些數(shù)據(jù)，如股票日漲跌幅，我們對數(shù)據(jù)進行簡單的分析。最簡單的是檢驗這一組數(shù)據(jù)是否服從假設的分布，如正態(tài)分布。這個問題是典型的單樣本假設檢驗問題，最為常見的解決方案是采用K-S檢驗（ Kolmogorov-Smirnov test）。單樣本K-S檢驗的原假設是給定的數(shù)據(jù)來自和原假設分布相同的分布，在SciPy中提供了kstest函數(shù)，參數(shù)分別是數(shù)據(jù)、擬檢驗的分布名稱和對應的參數(shù)：
In [5]:

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.045 p-value = 0.8195

假設檢驗的$p$-value值很大（在原假設下，$p$-value是服從[0, 1]區(qū)間上的均勻分布的隨機變量，可參考http://en.wikipedia.org/wiki/P-value ），因此我們接受原假設，即該數(shù)據(jù)通過了正態(tài)性的檢驗。在正態(tài)性的前提下，我們可進一步檢驗這組數(shù)據(jù)的均值是不是0。典型的方法是$t$檢驗（$t$-test），其中單樣本的$t$檢驗函數(shù)為ttest_1samp：
In [6]:

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 =  4.761, p-value = 0.0000

我們看到$p$-value$ < 0.05$，即給定顯著性水平0.05的前提下，我們應拒絕原假設：數(shù)據(jù)的均值為0。我們再生成一組數(shù)據(jù)，嘗試一下雙樣本的$t$檢驗（ttest_ind）：
In [7]:

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 =  5.565, p-value = 0.0000
注意，這里我們生成的第二組數(shù)據(jù)樣本大小、方差和第一組均不相等，在運用$t$檢驗時需要使用Welch's $t$-test，即指定ttest_ind中的equal_var=False。我們同樣得到了比較小的$p$-value$，在顯著性水平0.05的前提下拒絕原假設，即認為兩組數(shù)據(jù)均值不等。

stats還提供其他大量的假設檢驗函數(shù)，如bartlett和levene用于檢驗方差是否相等；anderson_ksamp用于進行Anderson-Darling的K-樣本檢驗等。

2.3 其他函數(shù)?

有時需要知道某數(shù)值在一個分布中的分位，或者給定了一個分布，求某分位上的數(shù)值。這可以通過cdf和ppf函數(shù)完成：
In [8]:

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)$的六階原點矩：
In [9]:

stats.norm.moment(6, loc=0, scale=1)

Out[9]:

15.000000000895332

describe函數(shù)提供對數(shù)據(jù)集的統(tǒng)計描述分析，包括數(shù)據(jù)樣本大小，極值，均值，方差，偏度和峰度：
In [10]:

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: -4.12414564687
Maximum value is: 4.82577602489
Arithmetic mean is: 0.0962913592209
Unbiased variance is: 2.88719292463
Biased skewness is: -0.00256548794681
Biased kurtosis is: -0.317463421177

當我們知道一組數(shù)據(jù)服從某些分布的時候，可以調(diào)用fit函數(shù)來得到對應分布參數(shù)的極大似然估計（MLE, maximum-likelihood estimation）。以下代碼示例了假設數(shù)據(jù)服從正態(tài)分布，用極大似然估計分布參數(shù)：
In [11]:

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.249880829912
MLE of data standard deviation:1.89195303507

pearsonr和spearmanr可以計算Pearson和Spearman相關(guān)系數(shù)，這兩個相關(guān)系數(shù)度量了兩組數(shù)據(jù)的相互線性關(guān)聯(lián)程度：
In [12]:

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.0262911931014
Spearman's rank correlation coefficient: -0.0262911931014

其中的$p$-value表示原假設（兩組數(shù)據(jù)不相關(guān)）下，相關(guān)系數(shù)的顯著性。

最后，在分析金融數(shù)據(jù)中使用頻繁的線性回歸在SciPy中也有提供，我們來看一個例子：
In [13]:

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.44515601191
Intercept of fitted model is: 1.91080684516
R-squared: 0.798786910173
在前面的鏈接中，可以查到大部分stat中的函數(shù)，本節(jié)權(quán)作簡單介紹，挖掘更多功能的最好方法還是直接讀原始的文檔。另外，StatsModels（http://statsmodels.sourceforge.net ）模塊提供了更為專業(yè)，更多的統(tǒng)計相關(guān)函數(shù)。若在SciPy沒有滿足需求，可以采用StatsModels。

三、優(yōu)化部分?

優(yōu)化問題在投資中可謂是根本問題，如果手上有眾多可選的策略，應如何從中選擇一個“最好”的策略進行投資呢？這時就需要用到一些優(yōu)化技術(shù)針對給定的指標進行尋優(yōu)。隨著越來越多金融數(shù)據(jù)的出現(xiàn)，機器學習逐漸應用在投資領域，在機器學習中，優(yōu)化也是十分重要的一個部分。以下介紹一些常見的優(yōu)化方法，雖然例子是人工生成的，不直接應用于實際金融數(shù)據(jù)，我們希望讀者在后面遇到優(yōu)化問題時，能夠從這些簡單例子迅速上手解決。
 
3.1 無約束優(yōu)化問題?

所謂的無約束優(yōu)化問題指的是一個優(yōu)化問題的尋優(yōu)可行集合是目標函數(shù)自變量的定義域，即沒有外部的限制條件。例如，求解優(yōu)化問題 [
minimizef(x)=x2?4.8x+1.2
] 就是一個無約束優(yōu)化問題，而求解 [
minimizef(x)=x2?4.8x+1.2 subject tox≥0

]

則是一個帶約束的優(yōu)化問題。更進一步，我們假設考慮的問題全部是凸優(yōu)化問題，即目標函數(shù)是凸函數(shù)，其自變量的可行集是凸集。（詳細定義可參考斯坦福大學Stephen Boyd教授的教材convex optimization，下載鏈接：http://stanford.edu/~boyd/cvxbook ）

我們以Rosenbrock函數(shù) [ f(\mathbf{x}) = \sum{i=1}^{N-1} 100 (x_i - x{i-1}^2)^2 + (1 - x_{i-1})^2 ] 作為尋優(yōu)的目標函數(shù)來簡要介紹在SciPy中使用優(yōu)化模塊scipy.optimize。

首先需要定義一下這個Rosenbrock函數(shù)：
In [14]:

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單純形法?

單純形法是運籌學中介紹的求解線性規(guī)劃問題的通用方法，這里的Nelder-Mead單純形法與其并不相同，只是用到單純形的概念。設定起始點$\mathbf{x}_0 = (1.3, 0.7, 0.8, 1.9, 1.2)$，并進行最小化的尋優(yōu)。這里‘xtol’表示迭代收斂的容忍誤差上界：
In [15]:

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函數(shù)的性質(zhì)比較好，簡單的優(yōu)化方法就可以處理了，還可以在minimize中使用method='powell'來指定使用Powell's method。這兩種簡單的方法并不使用函數(shù)的梯度，在略微復雜的情形下收斂速度比較慢，下面讓我們來看一下用到函數(shù)梯度進行尋優(yōu)的方法。
 
3.1.2 Broyden-Fletcher-Goldfarb-Shanno法?

Broyden-Fletcher-Goldfarb-Shanno（BFGS）法用到了梯度信息，首先求一下Rosenbrock函數(shù)的梯度：

[ \begin{split} \frac{\partial f}{\partial xj} &= \sum{i=1}^N 200(xi - x{i-1}^2)(\delta{i,j} - 2x{i-1}\delta{i-1,j}) -2(1 - x{i-1})\delta_{i-1,j} \ &= 200(xj - x{j-1}^2) - 400xj(x{j+1} - x_j^2) - 2(1 - x_j) \end{split}] 其中當$i=j$時，$\delta_{i,j} = 1$，否則$\delta_{i,j} = 0$。

邊界的梯度是特例，有如下形式： [ \begin{split} \frac{\partial f}{\partial x_0} &= -400x_0(x_1 - x_0^2) - 2(1 - x_0), \ \frac{\partial f}{\partial x{N-1}} &= 200(x{N-1} - x_{N-2}^2) \end{split}]

我們可以如下定義梯度向量的計算函數(shù)了：
In [16]:

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函數(shù)中通過參數(shù)jac指定：
In [17]:

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矩陣（二階導數(shù)矩陣）。牛頓法大致的思路是采用泰勒展開的二階近似： [ f(\mathbf{x}) \approx f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0) + \frac{1}{2}(\mathbf{x} - \mathbf{x}_0)^T\mathbf{H}(\mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0) ] 其中$\mathbf{H}(\mathbf{x}_0)$表示二階導數(shù)矩陣。若Hessian矩陣是正定的，函數(shù)的局部最小值可以通過使上面的二次型的一階導數(shù)等于0來獲取，我們有： [ \mathbf{x}_{\mathrm{opt}} = \mathbf{x}_0 - \mathbf{H}^{-1}\nabla f ]

這里可使用共軛梯度近似Hessian矩陣的逆矩陣。下面給出Rosenbrock函數(shù)的Hessian矩陣元素通式：

[ \begin{split} H{i,j} = \frac{\partial^2 f}{\partial x_i \partial x_j} &= 200(\delta{i,j} - 2x{i-1}\delta{i-1,j}) - 400xi(\delta{i+1,j} - 2xi\delta{i,j}) - 400\delta{i,j}(x{i+1} - xi^2) + 2\delta{i,j}, \ &= (202 + 1200xi^2 - 400x{i+1}) \delta{i,j} - 400x_i\delta{i+1,j} - 400x{i-1}\delta{i-1,j} \end{split}] 其中$i,j \in [1, N-2]$。其他邊界上的元素通式為： [ \begin{split} \frac{\partial^2 f}{\partial x_0^2} &= 1200x_0^2 - 400x_1 + 2, \ \frac{\partial^2 f}{\partial x_0 \partial x_1} = \frac{\partial^2 f}{\partial x_1 \partial x_0} &= -400x_0, \ \frac{\partial^2 f}{\partial x{N-1} \partial x{N-2}} = \frac{\partial^2 f}{\partial x{N-2} \partial x{N-1}} &= -400x_{N-2}, \ \frac{\partial^2 f}{\partial x_{N-1}^2} &= 200. \end{split}]

例如，當$N=5$時的Hessian矩陣為：

[ \mathbf{H} =
[1200x20?400x1+2?400x0000 ?400x0202+1200x21?400x2?400x100 0?400x1202+1200x22?400x3?400x20 00?400x2202+1200x23?400x4?400x3 000?400x3200]

]

為使用牛頓共軛梯度法，我們需要提供一個計算Hessian矩陣的函數(shù)：
In [18]:

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

In [19]:

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 (Hessian):
  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])

對于一些大型的優(yōu)化問題，Hessian矩陣將異常大，牛頓共軛梯度法用到的僅是Hessian矩陣和一個任意向量的乘積，為此，用戶可以提供兩個向量，一個是Hessian矩陣和一個任意向量$\mathbf{p}$的乘積，另一個是向量$\mathbf{p}$，這就減少了存儲的開銷。記向量$\mathbf{p} = (p_1, \ldots, p_{N-1})$，可有

[ \mathbf{H(x)p} = \begin{bmatrix} (1200x0^2 - 400x_1 + 2)p_0 -400x_0p_1 \ \vdots \ -400x{i-1}p{i-1} + (202 + 1200x_i^2 - 400x{i+1})pi - 400x_ip{i+1} \ \vdots \ -400x{N-2}p{N-2} + 200p_{N-1} \end{bmatrix} ]

我們定義如下函數(shù)并使用牛頓共軛梯度方法尋優(yōu)：
In [20]:

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