OLS回归分析原理实战及结果解析-python3

OLS 回归

变量之间存在着相关关系，比如，人的身高和体重之间存在着关系，一般来说，人高一些，体重要重一些，身高和体重之间存在的是不确定性的相关关系。回归分析是研究相关关系的一种数学工具，它能帮助我们从一个变量的取值区估计另一个变量的取值。

OLS(最小二乘法)主要用于线性回归的参数估计，它的思路很简单，就是求一些使得实际值和模型估值之差的平方和达到最小的值，将其作为参数估计值。就是说，通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法可用于曲线拟合，其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

简单的线性回归

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.0, 9.]

和实际的回归系数非常接近。

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

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)) plt.show() 

高次模型的回归

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, x2)) 

使用 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., 0., 9.]

获取全部摘要

print(results.summary()) 

哑变量

一般而言，有连续取值的变量叫做连续变量，它们的取值可以是任何的实数，或者是某一区间里的任何实数，比如股价、时间、身高。但有些性质不是连续的，只有有限个取值的可能性，一般是用于分辨类别，比如性别、婚姻情况、股票所属行业，表达这些变量叫做分类变量。在回归分析中，我们需要将分类变量转化为哑变量(dummy variable)。

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 中分类变量的类别总数。

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') plt.show() 

这里要指出，哑变量是和其他自变量并行的影响因素，也就是说，哑变量和原先的 x 同时影响了回归的结果。初学者往往会误解这一点，认为哑变量是一个选择变量：也就是说，上图中给出的回归结果，是在只做了一次回归的情况下完成的，而不是分成３段进行３次回归。哑变量的取值藏在其他的三个维度中。可以理解成：上图其实是将高元的回归结果映射到平面上之后得到的图。

简单应用

data = get_price(['000001.XSHG', '.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['.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()) 

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') 

全文代码

# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm #OLS 回归 nsample = 100 #x1 的值从 0 到 10 等差排列 x = np.linspace(0, 10, nsample) #在 array 上加入一列常项1。 X = sm.add_constant(x) beta = np.array([1, 10]) e = np.random.normal(size=nsample) y = np.dot(X, beta) + e model = sm.OLS(y,X) results = model.fit() print(results.params) print(results.summary()) y_fitted = results.fittedvalues 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') ax.axis((-0.05, 2, -1, 25)) plt.show() #高次模型的回归 nsample = 100 X = np.column_stack((x, x2)) X = sm.add_constant(X) beta = np.array([1, 0.1, 10]) e = np.random.normal(size=nsample) y = np.dot(X, beta) + e model = sm.OLS(y,X) results = model.fit() print(results.params) print(results.summary()) y_fitted = results.fittedvalues 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') ax.axis((-0.05, 2, -1, 25)) plt.show() #哑变量 nsample = 50 groups = np.zeros(nsample, int) groups[20:40] = 1 groups[40:] = 2 dummy = sm.categorical(groups, drop=True) x = np.linspace(0, 20, nsample) X = np.column_stack((x, dummy)) X = sm.add_constant(X) beta = [10, 1, 1, 3, 8] e = np.random.normal(size=nsample) y = np.dot(X, beta) + e 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') plt.show() data = get_price(['000001.XSHG', '.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['.XSHE'].values 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) x = np.array(x_pct) X = sm.add_constant(x) y = np.array(y_pct) results = sm.OLS(y, X).fit() 

summary结果中提供了很多关于拟合的信息，下面是这些描述信息的含义：
第一个表左边部分是关于拟合的基本信息:

第一个表的右边部分显示的是拟合的好坏情况：

第二个表显示的是拟合系数信息：

最后一个表显示的是对残差分布的统计检验评估：

部分较常用的结果数值提取具体操作示例如下

import statsmodels.api as sm # 模型训练 model = sm.OLS(y, x).fit() # 查看模型结果 print(model.summary()) 

提取元素-回归系数类

# 提取回归系数 model.params # 提取回归系数标准差 model.bse # 提取回归系数p值 model.pvalues # 提取回归系数t值 model.tvalues # 提取回归系数置信区间 默认5%，括号中可填具体数字 比如0.05, 0.1 model.conf_int() # 提取模型预测值 model.fittedvalues # 提取残差 model.resid # 模型自由度（系数自由度） model.df_model # 残差自由度（样本自由度） model.df_resid # 模型样本数量 model.nobs 

模型评价类

# 提取R方 model.rsquared # 提取调整R方 model.rsquared_adj # 提取AIC model.aic # 提取BIC model.bic # 提取F-statistic model.fvalue # 提取F-statistic 的pvalue model.f_pvalue # 模型mse model.mse_model # 残差mse model.mse_resid # 总体mse model.mse_total 

下面是不太常用的计量经济学方面的系数

# 协方差矩阵比例因子 model.scale # White异方差稳健标准误 model.HC0_se # MacKinnon和White（1985）的异方差稳健标准误 model.HC1_se # White异方差矩阵 model.cov_HC0 # MacKinnon和White（1985）的异方差矩阵 model.cov_HC1