[量化学堂-数学知识]多元线性回归

线性回归
逐步回归
标签: #<Tag:0x00007f5c00b8b8a8> #<Tag:0x00007f5c00b8b5b0>

(iQuant) #1

导语:多元线性回归的应用比较普遍,本文将对其做相关介绍。


金融理论从资本资产定价模型(CAPM)发展到套利定价理论(APT),在数理统计方面就是从应用一元线性回归发展到应用多元线性回归。在实际运用中,多元线性回归比较普遍

一元线性回归研究的是一个因变量和一个自变量的线性关系的模型,多元线性回归研究的是一个因变量和多个自变量的线性关系的模型。 多元线性回归模型表示为:

$$Y_i=\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+...++\beta_KX_{Ki}+u_i$$
其中,$i=1,2,...,n,n表示样本容量$,$K$表示自变量的个数。

与一元线性回归分析相同,其基本思想是根据 最小二乘($OLS$)原理,求解$\beta_0、\beta_1、\beta_2、...、\beta_K$使得全部观测值$Y_i$与回归值$\hat{Y_i}$的残差平方和达到最小值,该方法旨在最小化预测和观测之间的平方误差 ,$\sum_{i=0}^n\hat{u_i}^2$,平方项使得正的残差和负的残差同样被认为是糟糕,并且将其放大。残差平方和表示如下:

$$Q=\sum_{i=0}^n(Y_i-\hat{Y_i})^2=\sum_{i=0}^n(Y_i-(\hat{\beta_0}+\hat{\beta_1}X_{1i}+\hat{\beta_2}X_{2i}+...+\hat{\beta_K}X_{Ki}))^2$$

本文将介绍多元线性回归模型的相关理论和实际运用,感兴趣的小朋友可以 克隆到自己的策略研究平台进行研究。


附件:多元线性回归的使用

克隆策略
In [3]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels import regression
import matplotlib.pyplot as plt
In [1]:
Y = np.array([1, 3.5, 4, 8, 12])
Y_hat = np.array([1, 3, 5, 7, 9])

print('Error ' + str(Y_hat - Y))

# 计算残差平方
SE = (Y_hat - Y) ** 2
# 残差平方
print('Squared Error' + str(SE))
# 残差平方和
print('Sum Squared Error ' + str(np.sum(SE)))
Error [ 0.  -0.5  1.  -1.  -3. ]
Squared Error[ 0.    0.25  1.    1.    9.  ]
Sum Squared Error 11.25

一旦我们使用多元线性回归来确定回归的系数,我们将能够使用$X$的新观察值来预测$Y$的值。

每个系数$\beta_j$告诉我们,如果在保持所有其他因变量不变的情况下将$X_j$改变1个百分比,$Y_i$将会改变多少。这使我们可以分离不同自变量变化导致因变量变化所产生的边际贡献。

In [4]:
# 我们首先人为构建一个我们知道精确关系的Y,X1和X2
X1 = np.arange(100)

#  X2 = X1^2 + X1
X2 = np.array([i ** 2 for i in range(100)]) + X1

Y = X1 + X2

plt.plot(X1, label='X1')
plt.plot(X2, label='X2')
plt.plot(Y, label='Y')
plt.legend();
In [5]:
# 使用 column_stack连接X1和X2这两个变量, 然后将单位向量作为截距项
X = sm.add_constant( np.column_stack( (X1, X2) ) )

# 运行回归模型
results = regression.linear_model.OLS(Y, X).fit()

print('Beta_0:', results.params[0])
print('Beta_1:', results.params[1])
print('Beta_2:', results.params[2])
Beta_0: 1.36424205266e-12
Beta_1: 1.0
Beta_2: 1.0

可以看出,$X_1$的系数为1,这是因为如果在保持$X_2$不变的情况下将$X_1$增加1%,则$Y$也增加1%。 可以看出,多元线性回归能够分析不同变量的边际贡献。

如果我们使用一元线性回归来分析两个变量的关系,我们很可能会得到一个比较高的贝塔值,这样的分析是有失偏颇的,因此需要加入另外一个变量。请看下面这个例子

In [6]:
# 获得贵州茅台、中国平安、沪深300的数据
start_date = '2014-01-01'
end_date = '2015-01-01'

asset1 =  D.history_data('600519.SHA',start_date,end_date,fields=['close']).set_index('date')['close']
asset2 =  D.history_data('000001.SZA',start_date,end_date,fields=['close']).set_index('date')['close']
benchmark = D.history_data('000300.SHA',start_date,end_date,fields=['close']).set_index('date')['close']
# First, run a linear regression on the two assets
slr = regression.linear_model.OLS(asset1, sm.add_constant(asset2)).fit()
print('SLR beta of asset2:', slr.params[1])
SLR beta of asset2: 0.664051496894
In [8]:
# 将asset2和benchmark两个变量都看成自变量,然后进行多元回归
mlr = regression.linear_model.OLS(asset1, sm.add_constant(np.column_stack((asset2, benchmark)))).fit()

prediction = mlr.params[0] + mlr.params[1]*asset2 + mlr.params[2]*benchmark
prediction.name = 'Prediction'

print('MLR beta of asset2:', mlr.params[1], '\nMLR beta of 000300:', mlr.params[2])
MLR beta of asset2: -0.0133063437023 
MLR beta of 000300: 0.237182936686

可以看出,之前只是一元线性回归的时候获取了一个较高的$\beta$(贝塔值)

得到分析结果以后,下一步是看看我们是否可以相信结果。一个简单有效的方法就是将变量值和预测者绘制在图表进行观察

In [82]:
asset1.name = 'asset1'
asset2.name = 'asset2'
benchmark.name = 'benchmark'

asset1.plot()
asset2.plot()
benchmark.plot()
prediction.plot(color='y')
plt.xlabel('return')
plt.legend(bbox_to_anchor=(1,1), loc=2);
In [9]:
# 只看观测值和预测值,可以看出预测的走势还是比较接近的
asset1.plot()
prediction.plot(color='y')
plt.xlabel('Price')
plt.legend();

我们可以从回归返回的结果中得到详细统计信息

In [84]:
mlr.summary()
Out[84]:
OLS Regression Results
Dep. Variable: asset1 R-squared: 0.461
Model: OLS Adj. R-squared: 0.456
Method: Least Squares F-statistic: 103.4
Date: Fri, 14 Apr 2017 Prob (F-statistic): 3.53e-33
Time: 11:02:09 Log-Likelihood: -1421.9
No. Observations: 245 AIC: 2850.
Df Residuals: 242 BIC: 2860.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 355.0840 39.897 8.900 0.000 276.494 433.674
x1 -0.0133 0.189 -0.070 0.944 -0.385 0.359
x2 0.2372 0.064 3.706 0.000 0.111 0.363
Omnibus: 42.022 Durbin-Watson: 0.039
Prob(Omnibus): 0.000 Jarque-Bera (JB): 58.618
Skew: -1.087 Prob(JB): 1.87e-13
Kurtosis: 4.006 Cond. No. 1.94e+04

下面是一个关于如何 选择模型的一个例子

当你决定应该包括哪些自变量的最佳模型时,有几种不同的方法来确定。如果你使用太多的自变量,你可能会增加模型 过拟合的风险,但如果你使用的太少,你可能会 欠拟合。决定最佳模型的一个常用方法是 逐步回归。前向逐步回归从空模型开始,并测试每个变量,选择可以形成最佳模型的变量,通常用AIC或BIC判定(该值越低越好)。然后,依次添加一个其余变量,在回归中测试每个随后的变量组合,并在每个步骤中计算AIC或BIC值。在回归结束时,选择 具有AIC最低的模型,并将其确定为最终模型。但这确实也有局限性。它不能测试所有变量中的的每一个全部可能的组合,因此如果在逐步回归提前删除了某个特定的变量,那么最终模型可能不是理论最佳模型。因此,逐步回归应结合你对模型变量的初始判断。

In [10]:
X1 = np.arange(100)
X2 = [i**2 for i in range(100)] - X1
X3 = [np.log(i) for i in range(1, 101)] + X2
X4 = 5 * X1
Y = 2 * X1 + 0.5 * X2 + 10 * X3 + X4

plt.plot(X1, label='X1')
plt.plot(X2, label='X2')
plt.plot(X3, label='X3')
plt.plot(X4, label='X4')
plt.plot(Y, label='Y')
plt.legend();
In [11]:
results = regression.linear_model.OLS(Y, sm.add_constant(np.column_stack((X1,X2,X3,X4)))).fit()

print("Beta_0: ", results.params[0])
print("Beta_1: ", results.params[1])
print("Beta_2: ", results.params[2])
print("Beta_3: ", results.params[3])
print("Beta_4: ", results.params[4])
Beta_0:  -6.36646291241e-12
Beta_1:  0.269230769231
Beta_2:  0.499999999994
Beta_3:  10.0
Beta_4:  1.34615384615
In [12]:
data = pd.DataFrame(np.column_stack((X1,X2,X3,X4)), columns = ['X1','X2','X3','X4'])
response = pd.Series(Y, name='Y')
In [13]:
# 逐步回归
def forward_aic(response, data):
    # This function will work with pandas dataframes and series
    
    # Initialize some variables
    explanatory = list(data.columns)
    selected = pd.Series(np.ones(data.shape[0]), name="Intercept")
    current_score, best_new_score = np.inf, np.inf
    
    # Loop while we haven't found a better model
    while current_score == best_new_score and len(explanatory) != 0:
        
        scores_with_elements = []
        count = 0
        
        # For each explanatory variable
        for element in explanatory:
            # Make a set of explanatory variables including our current best and the new one
            tmp = pd.concat([selected, data[element]], axis=1)
            # Test the set
            result = regression.linear_model.OLS(Y, tmp).fit()
            score = result.aic
            scores_with_elements.append((score, element, count))
            count += 1
        
        # Sort the scoring list
        scores_with_elements.sort(reverse = True)
        # Get the best new variable
        best_new_score, best_element, index = scores_with_elements.pop()
        if current_score > best_new_score:
            # If it's better than the best add it to the set
            explanatory.pop(index)
            selected = pd.concat([selected, data[best_element]],axis=1)
            current_score = best_new_score
    # Return the final model
    model = regression.linear_model.OLS(Y, selected).fit()
    return model
In [14]:
result = forward_aic(Y, data)
result.summary()
Out[14]:
OLS Regression Results
Dep. Variable: y R-squared: 1.000
Model: OLS Adj. R-squared: 1.000
Method: Least Squares F-statistic: 3.092e+26
Date: Tue, 09 May 2017 Prob (F-statistic): 0.00
Time: 19:42:06 Log-Likelihood: 1700.7
No. Observations: 100 AIC: -3393.
Df Residuals: 96 BIC: -3383.
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -1.455e-11 7.01e-09 -0.002 0.998 -1.39e-08 1.39e-08
X3 10.0000 4.24e-09 2.36e+09 0.000 10.000 10.000
X1 0.2692 1.3e-11 2.08e+10 0.000 0.269 0.269
X2 0.5000 4.24e-09 1.18e+08 0.000 0.500 0.500
X4 1.3462 6.48e-11 2.08e+10 0.000 1.346 1.346
Omnibus: 14.070 Durbin-Watson: 0.001
Prob(Omnibus): 0.001 Jarque-Bera (JB): 9.981
Skew: -0.647 Prob(JB): 0.00680
Kurtosis: 2.152 Cond. No. 6.22e+17

在模型的构建中,可以很明显地看到变量$X_4$与变量$X_1$密切相关,只是将其乘以一个标量而已。 然而,逐步回归的方法并没有捕捉到这个信息,而是简单地调整了$X_1$项的系数。 但是我们自己的判断会说,我们建模时应该剔除变量$X_4$,因此这显示出逐步回归的局限性。

还有其他方法来诊断模型,比如给予更复杂模型的不同程度的 惩罚。 我们以后再做深入介绍。


本文由BigQuant宽客学院推出,版权归BigQuant所有,转载请注明出处。


社区干货与精选整理(持续更新中...)