[源码分享] 【量化实验室q.datayes.com】量化分析师的Python日记【第4天：scipy篇】 [推广有奖]

上一篇介绍了numpy,本篇中着重介绍一下另一个量化金融中常用的库 scipy

一、SciPy概述

前篇已经大致介绍了NumPy，接下来让我们看看SciPy能做些什么。NumPy替我们搞定了向量和矩阵的相关操作，基本上算是一个高级的科学计算器。SciPy基于NumPy提供了更为丰富和高级的功能扩展，在统计、优化、插值、数值积分、时频转换等方面提供了大量的可用函数，基本覆盖了基础科学计算相关的问题。

在量化分析中，运用最广泛的是统计和优化的相关技术，本篇重点介绍SciPy中的统计和优化模块，其他模块在随后系列文章中用到时再做详述。

本篇会涉及到一些矩阵代数，如若感觉不适，可考虑跳过第三部分或者在理解时简单采用一维的标量代替高维的向量。

首先还是导入相关的模块，我们使用的是SciPy里面的统计和优化部分：

1

import numpy as np

2

import scipy.stats as stats

3

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的贝塔分布随机数：

1

rv_unif = stats.uniform.rvs(size=10)

2

print rv_unif

3

rv_beta = stats.beta.rvs(size=10, a=4, b=2)

4

print rv_beta

[ 0.6419336   0.48403001 0.89548809  0.73837498  0.65744886 0.41845577

[ 0.82164685  0.69563836 0.74207073  0.94348192  0.82979411 0.87013796

在每个随机分布的生成函数里，都内置了默认的参数，如均匀分布的上下界默认是0和1。可是一旦需要修改这些参数，每次生成随机都要敲这么老长一串有点麻烦，能不能简单点？SciPy里头有一个Freezing的功能，可以提供简便版本的命令。SciPy.stats支持定义出某个具体的分布的对象，我们可以做如下的定义，让beta直接指代具体参数a=4和b=2的贝塔分布。为让结果具有可比性，这里指定了随机数的生成种子，由NumPy提供。

1

np.random.seed(seed=2015)

2

rv_beta = stats.beta.rvs(size=10, a=4, b=2)

3

print "method 1:"

4

print rv_beta

5

6

np.random.seed(seed=2015)

7

beta = stats.beta(a=4, b=2)

8

print "method 2:"

9

print beta.rvs(size=10)

method 1:

[ 0.43857338  0.9411551  0.75116671  0.92002864  0.62030521 0.56585548

method 2:

[ 0.43857338  0.9411551  0.75116671  0.92002864  0.62030521 0.56585548

invalid syntax(line 2)

2.2 假设检验

好了，现在我们生成一组数据，并查看相关的统计量（相关分布的参数可以在这里查到：http://docs.scipy.org/doc/scipy/reference/stats.html）：

1

norm_dist = stats.norm(loc=0.5, scale=2)

2

n = 200

3

dat = norm_dist.rvs(size=n)

4

print "mean of data is: " + str(np.mean(dat))

5

print "median of data is: " + str(np.median(dat))

6

print "standard deviation of data is: " + str(np.std(dat))

mean of data is:0.383309149888

median of data is:0.394980561217

standard deviationof data is: 2.00589851641

假设这个数据是我们获取到的实际的某些数据，如股票日涨跌幅，我们对数据进行简单的分析。最简单的是检验这一组数据是否服从假设的分布，如正态分布。这个问题是典型的单样本假设检验问题，最为常见的解决方案是采用K-S检验（ Kolmogorov-Smirnov test）。单样本K-S检验的原假设是给定的数据来自和原假设分布相同的分布，在SciPy中提供了kstest函数，参数分别是数据、拟检验的分布名称和对应的参数：

1

mu = np.mean(dat)

2

sigma = np.std(dat)

3

stat_val, p_val = stats.kstest(dat, 'norm', (mu, sigma))

4

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：

1

stat_val, p_val = stats.ttest_1samp(dat, 0)

2

print 'One-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val)

One-samplet-statistic D =  2.696, p-value = 0.0076

我们看到p-value<0.05，即给定显著性水平0.05的前提下，我们应拒绝原假设：数据的均值为0。我们再生成一组数据，尝试一下双样本的t检验（ttest_ind）：

1

norm_dist2 = stats.norm(loc=-0.2, scale=1.2)

2

dat2 = norm_dist2.rvs(size=n/2)

3

stat_val, p_val = stats.ttest_ind(dat, dat2, equal_var=False)

4

print 'Two-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val)

Two-samplet-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_ksamp用于进行Anderson-Darling的K-样本检验等。

原文太长无法全部发表，请看以下链接