线性回归—投资额（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_{}$