PythonStatsmodels统计包之OLS回归

PythonStatsmodels统计包之OLS回归
Statsmodels 是 Python 中⼀个强⼤的统计分析包，包含了回归分析、时间序列分析、假设检
验等等的功能。

Statsmodels 在计量的简便性上是远远不及 Stata 等软件的，但它的优点在于可以与 Python 的其他的任务（如 NumPy、Pandas）有效结合，提⾼⼯作效率。

在本⽂中，我们重点介绍最回归分析中最常⽤的 OLS（ordinary least square）功能。

当你需要在 Python 中进⾏回归分析时……
import statsmodels.api as sm！！！
在⼀切开始之前
上帝导⼊了 NumPy（⼤家都叫它囊派？我叫它囊辟），
import numpy as np
便有了时间。

上帝导⼊了 matplotlib，
import matplotlib.pyplot as plt
便有了空间。

上帝导⼊了 Statsmodels，
import statsmodels.api as sm
世界开始了。

简单 OLS 回归
假设我们有回归模型
Y=β0+β1X1+⋯+βnXn+ε,
并且有 k 组数据。

OLS 回归⽤于计算回归系数βi 的估值 b0,b1,…,bn，使误差平⽅
最⼩化。

statsmodels.OLS 的输⼊有 (endog, exog, missing, hasconst) 四个，我们现在只考虑前两个。

第⼀个输⼊ endog 是回归中的反应变量（也称因变量），是上⾯模型中的 y(t), 输⼊是⼀个长度为 k 的 array。

第⼆个输⼊ exog 则是回归变量（也称⾃变量）的值，即模型中的x1(t),…,xn(t)。

但是要注意，statsmodels.OLS 不会假设回归模型有常数项，所以我们应该假设模型是
并且在数据中，对于所有 t=1,…,k，设置 x0(t)=1。

因此，exog的输⼊是⼀个 k×(n+1) 的array，其中最左⼀列的数值全为 1。

往往输⼊的数据中，没有专门的数值全为1的⼀列，Statmodels 有直接解决这个问题的函数：sm.add_constant()。

它会在⼀个 array 左侧加上⼀列 1。

（本⽂中所有输⼊ array 的情况也可以使⽤同等的 list、pd.Series 或 pd.DataFrame。

）
确切地说，statsmodels.OLS 是 statsmodels.regression.linear_model ⾥的⼀个函数（从这个命名也能看出，statsmodel 有很多很多功能，其中的⼀项叫回归）。

它的输出结果是⼀个 statsmodels.regression.linear_model.OLS，只是⼀个类，并没有进⾏任何运算。

在 OLS 的模型之上调⽤拟合函数 fit()，才进⾏回归运算，并且得到 statsmodels.regression.linear_model.RegressionResultsWrapper，它包含了这组数据进⾏回归拟合的结果摘要。

调⽤ params 可以查看计算出的回归系数 b0,b1,…,bn。

简单的线性回归
上⾯的介绍绕了⼀个⼤圈圈，现在我们来看⼀个例⼦，假设回归公式是:
我们从最简单的⼀元模型开始，虚构⼀组数据。

⾸先设定数据量，也就是上⾯的 k 值。

nsample = 100
然后创建⼀个 array，是上⾯的 x1 的数据。

这⾥，我们想要 x1 的值从 0 到 10 等差排列。

x = np.linspace(0, 10, nsample)
使⽤ sm.add_constant() 在 array 上加⼊⼀列常项1。

X = sm.add_constant(x)
然后设置模型⾥的β0,β1，这⾥要设置成 1,10。

beta = np.array([1, 10])
然后还要在数据中加上误差项，所以⽣成⼀个长度为k的正态分布样本。

e = np.random.normal(size=nsample)
由此，我们⽣成反应项 y(t)。

y = np.dot(X, beta) + e
好嘞，在反应变量和回归变量上使⽤ OLS() 函数。

model = sm.OLS(y,X)
然后获取拟合结果。

results = model.fit()
再调取计算出的回归系数。

print(results.params)
得到
[ 1.04510666, 9.97239799]
和实际的回归系数⾮常接近。

当然，也可以将回归拟合的摘要全部打印出来。

print(results.summary())
得到
中间偏下的 coef 列就是计算出的回归系数。

我们还可以将拟合结果画出来。

先调⽤拟合结果的 fittedvalues 得到拟合的 y 值。

y_fitted = results.fittedvalues
然后使⽤ matplotlib.pyploft 画图。

⾸先设定图轴，图⽚⼤⼩为 8×6。

fig, ax = plt.subplots(figsize=(8,6))
画出原数据，图像为圆点，默认颜⾊为蓝。

ax.plot(x, y, 'o', label='data')
画出拟合数据，图像为红⾊带点间断线。

ax.plot(x, y_fitted, 'r--.',label='OLS')
放置注解。

ax.legend(loc='best')
得到
在⼤图中看不清细节，我们在 0 到 2 的区间放⼤⼀下，可以见数据和拟合的关系。

加⼊改变坐标轴区间的指令
ax.axis((-0.05, 2, -1, 25))
⾼次模型的回归
假设反应变量 Y 和回归变量 X 的关系是⾼次的多项式，即
我们依然可以使⽤ OLS 进⾏线性回归。

但前提条件是，我们必须知道 X 在这个关系中的所有次⽅数；⽐如，如果这个公式⾥有⼀个.5项，但我们对此并不知道，那么⽤线性回归的⽅法就不能得到准确的拟合。

虽然 X 和 Y 的关系不是线性的，但是 Y 和的关系是⾼元线性的。

也就是说，只要我们把⾼次项当做其他的⾃变量，即设。

那么，对于线性公式
可以进⾏常规的 OLS 回归，估测出每⼀个回归系数βi。

可以理解为把⼀元⾮线性的问题映射到⾼元，从⽽变成⼀个线性关系。

下⾯以
为例做⼀次演⽰。

⾸先设定数据量，也就是上⾯的 k 值。

nsample = 100
然后创建⼀个 array，是上⾯的 x1 的数据。

这⾥，我们想要 x1 的值从 0 到 10 等差排列。

x = np.linspace(0, 10, nsample)
再创建⼀个 k×2 的 array，两列分别为 x1 和 x2。

我们需要 x2 为 x1 的平⽅。

X = np.column_stack((x, x**2))
使⽤ sm.add_constant() 在 array 上加⼊⼀列常项 1。

X = sm.add_constant(X)
然后设置模型⾥的β0,β1,β2，我们想设置成 1,0.1,10。

beta = np.array([1, 0.1, 10])
然后还要在数据中加上误差项，所以⽣成⼀个长度为k的正态分布样本。

e = np.random.normal(size=nsample)
由此，我们⽣成反应项 y(t)，它与 x1(t) 是⼆次多项式关系。

y = np.dot(X, beta) + e
在反应变量和回归变量上使⽤ OLS() 函数。

model = sm.OLS(y,X)
然后获取拟合结果。

results = model.fit()
再调取计算出的回归系数。

print(results.params)
得到
[ 0.95119465, 0.10235581, 9.9998477]
获取全部摘要
print(results.summary())
得到
拟合结果图如下在横轴的 [0,2] 区间放⼤，可以看到
哑变量
⼀般⽽⾔，有连续取值的变量叫做连续变量，它们的取值可以是任何的实数，或者是某⼀区间⾥的任何实数，⽐如股价、时间、⾝⾼。

但有些性质不是连续的，只有有限个取值的可能性，⼀般是⽤于分辨类别，⽐如性别、婚姻情况、股票所属⾏业，表达这些变量叫做分类变量。

在回归分析中，我们需要将分类变量转化为哑变量(dummy variable)。

如果我们想表达⼀个有 d 种取值的分类变量，那么它所对应的哑变量的取值是⼀个 d 元组（可以看成⼀个长度为 d 的向量），其中有⼀个元素为 1，其他都是0。

元素呈现出 1 的位置就是变量所取的类别。

⽐如说，某个分类变量的取值是 {a,b,c,d}，那么类别 a 对应的哑变量是(1,0,0,0)，b 对应 (0,1,0,0)，c 对应(0,0,1,0)，d 对应 (0,0,0,1)。

这么做的⽤处是，假如 a、b、c、d 四种情况分别对应四个系数β0,β1,β2,β3，设 (x0,x1,x2,x3) 是⼀个取值所对应的哑变量，那么
可以直接得出相应的系数。

可以理解为，分类变量的取值本⾝只是分类，⽆法构成任何线性关系，但是若映射到⾼元的 0,1 点上，便可以⽤线性关系表达，从⽽进⾏回归。

Statsmodels ⾥有⼀个函数 categorical() 可以直接把类别 {0,1,…,d-1} 转换成所对应的元组。

确切地说，sm.categorical() 的输⼊有 (data, col, dictnames, drop) 四个。

其中，data 是⼀个 k×1 或 k×2 的 array，其中记录每⼀个样本的分类变量取值。

drop 是⼀个 Bool值，意义为是否在输出中丢掉样本变量的值。

中间两个输⼊可以不⽤在意。

这个函数的输出是⼀个k×d 的 array（如果 drop=False，则是k×(d+1)），其中每⼀⾏是所对应的样本的哑变量；这⾥ d 是 data 中分类变量的类别总数。

我们来举⼀个例⼦。

这⾥假设⼀个反应变量 Y 对应连续⾃变量 X 和⼀个分类变量 Z。

常项系数为 10，XX 的系数为 1；Z 有 {a,b,c}三个种类，其中 a 类有系数1，b 类有系数 3，c 类有系数 8。

也就是说，将 Z 转换为哑变量 (Z1,Z2,Z3)，其中 Zi 取值于 0,1，有线性公式
可以⽤常规的⽅法进⾏ OLS 回归。

我们按照这个关系⽣成⼀组数据来做⼀次演⽰。

先定义样本数量为 50。

nsample = 50
设定分类变量的 array。

前 20 个样本分类为 a。

groups = np.zeros(nsample, int)
之后的 20 个样本分类为 b。

groups[20:40] = 1
最后 10 个是 c 类。

groups[40:] = 2
转变成哑变量。

dummy = sm.categorical(groups, drop=True)
创建⼀组连续变量，是 50 个从 0 到 20 递增的值。

x = np.linspace(0, 20, nsample)
将连续变量和哑变量的 array 合并，并加上⼀列常项。

X = np.column_stack((x, dummy))
X = sm.add_constant(X)
定义回归系数。

我们想设定常项系数为 10，唯⼀的连续变量的系数为 1，并且分类变量的三种分类 a、b、c 的系数分别为 1,3,8。

beta = [10, 1, 1, 3, 8]
再⽣成⼀个正态分布的噪⾳样本。

e = np.random.normal(size=nsample)
最后，⽣成反映变量。

y = np.dot(X, beta) + e
得到了虚构数据后，放⼊ OLS 模型并进⾏拟合运算。

result = sm.OLS(y,X).fit()
print(result.summary())
再画图出来
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, result.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best')
得出
这⾥要指出，哑变量是和其他⾃变量并⾏的影响因素，也就是说，哑变量和原先的 x 同时影响了回归的结果。

初学者往往会误解这⼀点，认为哑变量是⼀个选择变量：也就是说，上图中给出的回归结果，是在只做了⼀次回归的情况下完成的，⽽不是分成３段进⾏３次回归。

哑变量的取值藏在其他的三个维度中。

可以理解成：上图其实是将⾼元的回归结果映射到平⾯上之后得到的图。

简单应⽤
我们来做⼀个⾮常简单的实际应⽤。

设 x 为上证指数的⽇收益率，y 为深证成指的⽇收益率。

通过对股票市场的认知，我们认为 x 和 y 有很强的线性关系。

因此
可以假设模型
并使⽤ Statsmodels 包进⾏ OLS 回归分析。

我们取上证指数和深证成指⼀年中的收盘价。

data = get_price(['000001.XSHG', '399001.XSHE'], start_date='2015-01-01', end_date='2016-01-01', frequency='daily', fields=['close'])['close'] x_price = data['000001.XSHG'].values
y_price = data['399001.XSHE'].values
计算两个指数⼀年内的⽇收益率，记载于 x_pct 和 y_pct 两个 list 中。

x_pct, y_pct = [], []
for i in range(1, len(x_price)):
x_pct.append(x_price[i]/x_price[i-1]-1)
for i in range(1, len(y_price)):
y_pct.append(y_price[i]/y_price[i-1]-1)
将数据转化为 array 的形式；不要忘记添加常数项。

x = np.array(x_pct)
X = sm.add_constant(x)
y = np.array(y_pct)
上吧，λu.λv.(sm.OLS(u,v).fit())！全靠你了！
results = sm.OLS(y, X).fit()
print(results.summary())
得到
恩，y=0.002+0.9991x，合情合理，或者⼲脆直接四舍五⼊到 y=x。

最后，画出数据和拟合线。

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, results.fittedvalues, 'r--', label="OLS")
ax.legend(loc='best')
结语
本篇⽂章中，我们介绍了 Statsmodels 中很常⽤ OLS 回归功能，并展⽰了⼀些使⽤⽅法。

线性回归的应⽤场景⾮常⼴泛。

在我们量化课堂应⽤类的内容中，也有相当多的策略内容采⽤线性回归的内容。

我们会将应⽤类⽂章中涉及线性回归的部分加上链接，链接到本篇⽂章中来，形成体系。

量化课堂在未来还会介绍Statsmodel 包其他的⼀些功能，敬请期待。

到JoinQuant查看代码并与作者交流讨论：。