线性回归—投资额（python、OLS最小二乘、残差图、DW检验）

年份序号	投资额	国民生产总值	物价指数
1	90.9	596.7	0.7167
2	97.4	637.7	0.7277
3	113.5	691.1	0.7436
4	125.7	756	0.7676
5	122.8	799	0.7906
6	133.3	873.4	0.8254
7	149.3	944	0.8679
8	144.2	992.7	0.9145
9	166.4	1077.6	0.9601
10	195	1185.9	1.0000
11	229.8	1326.4	1.0575
12	228.7	1434.2	1.1508
13	206.1	1594.2	1.2579
14	257.9	1718	1.3234
15	324.1	1918.3	1.4005
16	386.6	2163.9	1.5042
17	423	2417.8	1.6342
18	401.9	2631.7	1.7842
19	474.9	2954.7	1.9514
20	424.5	3073	2.0688

二、问题分析：

    以时间为序的数据称为时间序列。时间序列中同一变量的顺序观测值之间存在自相关，若采用普通回归模型直接处理，将会出现不良后果。因此，需要诊断并消除数据的自相关性，建立新的模型。另外许多经济数据在时间上有一定的滞后性，也会影响模型效果。
    本文按照正常建模流程来处理数据，分析每个模型的优缺点并进行比较。
三、基本模型：
变量描述：
    $t$~年份序号
    $y_t$~投资额
    $x_{1t}$~GNP
    $x_{2t}$~物价指数
画散点图：
    拿到数据的第一步，画散点图，观察投资额与GNP和物价指数的基本关系。

图1 投资额与GNP散点图

图2 投资额与物价指数散点图

可以看到，投资额与GNP及物价指数间均有很强的线性关系。所以，构建基本回归模型：

$$y_t={\beta }_0+{\beta }_1{x}_{1t}+{\beta }_2{x}_{2t}+{ϵ}_t$$

四、模型求解和代码展示：
    第一步，导入相关库，并解决画图中文问题。

#导入库
import cmath
import numpy as np 
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
#matplotlib画图中中文显示会有问题，需要这两行设置默认字体

    第二步，导入数据，变量描述。

#导入数据
touzi = pd.read_csv("touzi.csv")
touzi = pd.DataFrame(touzi)

#变量描述
t = touzi['年份序号']
yt = touzi['投资额']
x1t = touzi['国民生产总值']
x2t = touzi['物价指数']

    第三步，分析 $y_t$(投资额) 和 $x_{1t}$(GNP)及$x_{2t}$(物价指数) 之间的关系，通过散点图得以直观表达。分别是图1、图2。

    第四步(最重要一步),求解模型回归系数：${\beta }_0、{\beta }_1、{\beta }_2$.
    本文用OLS普通最小二乘法求解，在下列代码 (‘yt ~ x1t + x2t’, data = touzi)中，yt表示因变量。自变量则在yt~后边表示成变量相加的格式，模型会自动求解出变量前的回归系数。

#yt = b0 + b1x1t + b2x2t + e
model_1 = sm.formula.ols('yt ~ x1t + x2t', data = touzi).fit()  #OLS普通最小二乘法
model_1.summary()

图3 拟合参数

    下面是模型结果中部分参数的解释。

左边:
Dep.Variable:    输出变量的名称
Model:         模型名称
Method:        方法，其中 Least Squares 表示最小二乘法
Date:         日期
Time:         时间
No.Observations: 样本数目
Df Residuals:    残差自由度
Df Model:        模型参数个数，相当于输入的X的元素个数
右边:
R- squared:        可决系数,用来判断估计的准确性,范围在[0,1]越接近1,说明对y的解释能力越强,拟合越好
Adj-R- squared:    通过样本数量与模型数量对R-squared进行修正,奥卡姆剃刀原理,避免描述冗杂
F-statistic:       衡量拟合的显著性,重要程度
Prob(F-statistic): 当prob（F-statistic）<α时,表示拒绝原假设,即认为模型是显著的
Log likelihood:    对数似然比LLR
AIC:               衡量拟合优良性
BIC:               贝叶斯信息准则

    主要看此处的结果。Intercept表示 ${\beta }_0$ 的数值，$x_{1t}、x_{2t}$分别表示求得自身系数 ${\beta }_1、{\beta }_2$ 的数值。

coef:     系数
std err:  系数估计的基本标准误差
t:        t统计值,衡量系数统计显著程度的指标
P>|t|:    P值
[0.025,0.975]:   95％置信区间的下限和上限值

五、基本回归模型结果分析：

equation_yt = 'yt方程为:' + '\tyt' + '  =  ' + '363.7029' + '  +  ' + '0.6792' + '*x1t  ' + '-  ' + '972.7857' + '*x2t'
print(equation_yt)
pre_yt = 363.7029 + 0.6792*x1t - 972.7857*x2t
pre_yt = list(pre_yt)
MSE = cmath.sqrt(sum((yt - pre_yt)**2/(20)))
print('剩余标准差:\ts   = ',MSE)
print('R2 = ',0.990,',拟合度高')

    模型优点：$R^2=0.990$，拟合度高
    模型缺点：没有考虑时间序列数据的滞后性影响，可能忽视了随机误差存在自相关，如果存在自相关，用此模型会有不良后果。

六、自相关的定性诊断：
1、残差诊断法
    模型残差$e_t=y_t-{\hat{y}}_t$
    $e_t$为随机误差${ϵ}_t$的估计值。
    作残差$e_t$~$e_{t-1}$散点图，如果大部分点落在第1，3象限，则${ϵ}_t$存在正的自相关；如果大部分点落在第2，4象限，则${ϵ}_t$存在负的自相关。

#残差诊断法
et = yt - pre_yt    #残差
et = list(et)

#残差et 与 et_1 散点图
plt.ylim(-30, 20)
plt.xlim(-30, 20)
plt.plot([-30,20],[0,0],'b',linewidth=1)
plt.plot([0,0],[-30,20],'b',linewidth=1)
et_1 = et[1:20]
et_change = et[0:19]
plt.scatter(et_1,et_change,marker = '+',color = 'g')

图4 基本回归模型残差图

    自相关直观判断，基本回归模型的随机误差${ϵ}_t$存在正的自相关。注：此图不能直观看出大部分点落在第1，3象限，与案例中的图不太相同，不知原因，如有解决的朋友欢迎在评论区留言。下面展示案例图：

图4 案例中基本回归模型残差图

2、D-W统计量与D-W检验

$$DW=\frac{\sum _{t=2}^n{(e_t-e_{t-1})}^2}{\sum _{t=2}^n{e}_t^2}\approx 2[1-\frac{\sum _{t=2}^n{e}_t{e}_{t-1}}{\sum _{t=2}^n{e}_t^2}]=2(1-\hat{\rho })$$
    其中，$\hat{\rho }=\frac{\sum _{t=2}^n{e}_t{e}_{t-1}}{\sum _{t=2}^n{e}_t^2}$,此处$\approx$的前提是$n$较大。

    其中，$d_L$和$d_U$需要查询D-W分布表。$-1$ $\le$ $\hat{\rho }$$\le$ $1$ 即 $0$ $\le$ $DW$ $\le$ $4$ .
    注：以下DW计算数据均用课件给出数据，不再计算，有兴趣的朋友可以计算，因为用最小二乘估计参数与课件数据有偏差，计算出来的DW值不符合正自相关性，如果我解决了会及时修改。不过程序计算步骤都是没有问题的。

    $DW$计算程序：

#D-W统计量
DW = sum([(m-n)**2 for m,n in zip(et_change,et_1)])/sum([m**2 for m in et_change])
print('DW =',DW,' 解释:因误差，假设原模型有正自相关')
p = 1-DW/2
print('p =',p)

    原模型$D{W}_{old}=0.8754$，根据样本容量$n=20$,$\alpha =0.05$,通过查表得到$d_L=1.10$,$d_U=1.54$，因为$D{W}_{old}<d_L$，所以原模型有正自相关。
七、广义差分变换：
    原模型：$y_t={\beta }_0+{\beta }_1{x}_{1t}+{\beta }_2{x}_{2t}+{\epsilon }_t,{\epsilon }_t=\rho {\epsilon }_{t-1}+u_t$
    变换：$y_t^{\ast }=y_t-\rho y_{t-1},x_{it}^{\ast }=x_{it}-\rho x_{i,t-1},i=1,2$
    新模型：$y_t^{\ast }={\beta }_0^{\ast }+{\beta }_1{x}_{1t}^{\ast }+{\beta }_2{x}_{2t}^{\ast }+u_t,{\beta }_0^{\ast }={\beta }_0(1-\rho )$
    以 ${\beta }_0^{\ast },{\beta }_1,{\beta }_2$ 为回归系数的普通回归模型。

    DW检验后的处理步骤如下图：

    $\hat{\rho }=1-DW/2=0.5623$
    作变换：
    $y_i^{\ast }=y_t-0.5623y_{t-1}$
    $x_{it}^{\ast }=x_{it}-0.5623x_{i,t-1},i=1,2$

    给出程序：

#去掉数据后
yt_1 = [97.4,113.5,125.7,122.8,133.3,149.3,144.2,166.4,195.0,229.8,228.7,206.1,257.9,324.1,386.6,423.0,401.9,474.9,424.5]
x1t_1 = [637.7,691.1,756.0,799.0,873.4,944.0,992.7,1077.6,1185.9,1326.4,1434.2,1594.2,1718.0,1918.3,2163.9,2417.8,2631.7,2954.7,3073.0]
x2t_1 = [0.7277,0.7436,0.7676,0.7906,0.8254,0.8679,0.9145,0.9601,1.0,1.0575,1.1508,1.2579,1.3234,1.4005,1.5042,1.6342,1.7842,1.9514,2.0688]

#广义差分变换
#变换
yt_change = [a - p*b for a,b in zip(yt[1:20],yt_1)]
x1t_change = [c - p*d for c,d in zip(x1t[1:20],x1t_1)]
x2t_change = [e - p*f for e,f in zip(x2t[1:20],x2t_1)]

yt_change = pd.DataFrame(yt_change)
x1t_change = pd.DataFrame(x1t_change)
x2t_change = pd.DataFrame(x2t_change)
touzi['yt_change'] = yt_change
touzi['x1t_change'] = x1t_change
touzi['x2t_change'] = x2t_change

八、新模型求解及验证：
1、求解模型
$$y_t^{\ast }={\beta }_0^{\ast }+{\beta }_1{x}_{1t}^{\ast }+{\beta }_2{x}_{2t}^{\ast }+u_t$$.

#新模型
#yt_change = b0_change + b1_1x1t_change + b2_1x2t_change + u
model_2 = sm.formula.ols('yt_change ~ x1t_change + x2t_change', data = touzi).fit()  #OLS普通最小二乘法
model_2.summary()

#新模型
equation_yt_change = 'yt_change方程为:' + '  yt_change' + '  =  ' + '252.0659' + '  +  ' + '0.7176' + '*x1t_change  ' + '-  ' + '1040.8011' + '*x2t_change'
print(equation_yt_change)
pre_yt_change = 252.0659 + 0.7176*x1t_change - 1040.8011*x2t_change
MSE_change_1 = (yt_change - pre_yt_change)**2/(19)
MSE_change_2 = MSE_change_1.apply(lambda x: x.sum())
MSE_change_3 = cmath.sqrt(MSE_change_2)
print('剩余标准差:\t  s_change   = ',MSE_change_3,'< s = 12.2853')
print('R2_change  = ',0.991)

    总体效果良好，$s_{new}=7.6100<s_{old}=12.2853$

2、新模型的自相关检验
    新模型$D{W}_{new}=1.5751$,样本容量$n=19$,回归变量数目$k=3$,$\alpha =0.05$,通过查表得，临界值$d_L=1.08,d_U=1.53$.
    $d_U<D{W}_{new}<4-d_U$,新模型无自相关性。
    新模型：${\hat{y}}_t^{\ast }=252.0659+0.7176x_{1t}^{\ast }-1040.8011x_{2t}^{\ast }$.
    还原为原始变量：${\hat{y}}_t=252.0659+0.5623y_{t-1}+0.7176x_{1,t}-1040.8011x_{2,t}-0.4.35x_{1,t-1}+585.2425x_{2,t-1}$,称为一阶自回归模型。
    其中，$0.4.35x_{1,t-1}$为 $0.5623y_{t-1}\times \rho$;$585.2425x_{2,t-1}=1040.8011x_{2t}\times \rho$.

九、基本回归模型与一阶自回归模型比较：
1、残差图比较

#一阶自回归模型
onelevel_pre_yt =np.array([p*q for q in yt_1])  + np.array([0.7176*r for r in x1t[0:19]]) + np.array([p*(-0.7176)*s for s in x1t_1]) + np.array([-1040.8011*t for t in x2t[0:19]]) + np.array([p*1040.8011*v for v in x2t_1])

#新残差诊断法
et_one = list(yt[0:19]) - onelevel_pre_yt   #残差
et_one = list(et_one)

    将原模型与新模型（一阶自回归模型）残差归一化，把数据规模转换到[-20,20]，进行比较。

#对 et 数据归一化到[-20,20]区间范围
a_min = -20
b_max = 20
et_min = min(et[0:19])
et_max = max(et[0:19])
k=(b_max-a_min)/(et_max - et_min)
et_sta=[(a_min + k*(l - et_min)) for l in et[0:19]]

#对 et_one 数据归一化到[-20,20]区间范围
et_min_one = min(et_one)
et_max_one = max(et_one)
k2=(b_max-a_min)/(et_max_one - et_min_one)
et_sta_one=[(a_min + k2*(ll - et_min_one)) for ll in et_one]

    GNP（国民生产总值）的新旧模型残差比较。

#残差与国民生产总值的关系
plt.scatter(x1t[0:19],et_sta,marker = '+')
plt.scatter(x1t[0:19],et_sta_one,marker = '+',color='r')
plt.legend(["原模型","新模型"])

图5 GNP新旧模型残差比较图

    物价指数新旧模型残差比较。

#残差与物价指数的关系
plt.scatter(x2t[0:19],et_sta,marker = '+')
plt.scatter(x2t[0:19],et_sta_one,marker = '+',color='r')
plt.legend(["原模型","新模型"])

图6 物价指数新旧模型残差比较图

2、拟合图比较

#对 yt 数据归一化到[0,500]区间范围
c_min = 0
d_max = 500
yt_min = min(yt[0:19])
yt_max = max(yt[0:19])
k3=(d_max-c_min)/(yt_max - yt_min)
yt_sta=[(c_min + k3*(lll - yt_min)) for lll in yt[0:19]]

#对 onelevel_pre_yt 数据归一化到[-20,20]区间范围
onelevel_pre_yt_min = min(onelevel_pre_yt)
onelevel_pre_yt_max = max(onelevel_pre_yt)
k4=(d_max-c_min)/(onelevel_pre_yt_max - onelevel_pre_yt_min)
onelevel_pre_yt_sta=[(c_min + k4*(llll - onelevel_pre_yt_min)) for llll in onelevel_pre_yt]

#对 x1t 数据归一化到[0,20]区间范围
e_min = 0
f_max = 20
x1t_min = min(x1t[0:19])
x1t_max = max(x1t[0:19])
k5=(f_max-e_min)/(x1t_max - x1t_min)
x1t_sta=[(e_min + k5*(lllll - x1t_min)) for lllll in x1t[0:19]]

#对 et_one 数据归一化到[-20,20]区间范围
x2t_min = min(x2t[0:19])
x2t_max = max(x2t[0:19])
k6=(f_max-e_min)/(x2t_max - x2t_min)
x2t_sta=[(e_min + k6*(llllll - x2t_min)) for llllll in x2t[0:19]]

#拟合图比较
plt.scatter(x1t_sta[0:19],yt_sta[0:19],marker='o',c='',edgecolors='b')
plt.scatter(x2t_sta[0:19],yt_sta[0:19],marker='+',color='r',edgecolors='b')
plt.legend(["原模型","新模型"])

    新旧模型拟合图比较。

图7 新旧模型拟合图比较

    按道理，这里是要看出：一阶自回归模型残差$e$比基本回归模型要小，这里我也没有明确看出来。
    接下来就可以根据方程式进行投资额预测，本文不再展开描述。