用Python底层编写进行计量经济分析（四）：自相关（原因、结果、检验：DW检验、补救：广义线性回归）

最新推荐文章于 2024-03-14 14:49:16 发布

nbszg

最新推荐文章于 2024-03-14 14:49:16 发布

阅读量6.2k

点赞数 7

分类专栏：计量数据分析 python 文章标签：数据分析统计学统计模型

本文链接：https://blog.csdn.net/u011517132/article/details/105581371

版权

数据分析同时被 3 个专栏收录

13 篇文章 26 订阅

订阅专栏

python

8 篇文章 0 订阅

订阅专栏

计量

7 篇文章 2 订阅

订阅专栏

系列前面的文章：

写的比较仓促，代码中如有错误欢迎指正！

多元线性回归的基本假定：

模型符合线性模式
$X$ 满秩（无多重共线）
零均值价值： $E(ε_i|X_i)=0$ （自变量外生）
同方差： $Var(ε_i|X_i)=σ$
无自相关： $cov(ε_i, ε_j)=0$
球形扰动：ε_i是正态分布

若果模型违反了相应的假设就会犯对应的错误，我们在计量经济分析中的检验就是检验出是否可能犯了某一类错误，若果极有可能犯了一种错误时，我们应该怎么修正它，才能保证分析的结果是有效的。

一、异方差性

和前面一样，举个例子。如果，我们一个理论上的多元线性回归模型为： $Y_i=β_0+β_1*x_{1i}+β_2*x_{2i}+β_3*x_{3i}^2+u_i$
但我们再建模时只考虑 $x_{3}$ 的一阶情况，此时我们实际建立的模型变为 $Y_i=β_0+β_1*x_{1i}+β_2*x_{2i}+β_3*x_{3i}+u_i$
如果 $x_{3}$ 本身具有一定自相关，这样遗漏一个变量会对我们的估计产生后果？

1.1.自相关
多元线性回归要求无自相关，即 $cov(ε_i, ε_j)≠0， i不等于j$ ，如果我们估计的模型 $cov(ε_i, ε_j)≠0， i≠j$ 时，就出现了自相关问题。

1.2.自相关的原因

变量本身带有趋势
我们通常研究的时间序列自身就是带有一定趋势的，比如 $y_t$ 和 $y_{t-1}$ 本身就存在一定关系（此时可试试AR模型），此时就会出现异方差问题。
变量的之后效应
通常，一个变量对另一个变量的影响并不会在当期表现出来，往往需要之后一期或几期才会表现出来，如果忽视了这些之后效应，模型也可能出现自相关问题。
模型设定偏差
如果我们的模型设定存在偏差，比如真实模型为：
$Y_i=β_0+β_1*x_{1i}+β_2*x_{2i}+β_3*x_{3i}^2+u_i$
我们建模时，建立的模型为
$Y_i=β_0+β_1*x_{1i}+β_2*x_{2i}++v_i$
此时 $v_i=β_3*x_{3i}^2+u_i$ ， $v_i$ 包含了 $x_{3i}^2$ 的趋势。此时可能出现自相关

1.2.自相关的后果

自相关和异方差都是在扰动项的出现了问题（深究的话都是 $Ω$ 和原始形态不一样），因此自相关产生的后果与异方差产生的后果基本一样。

参数估计仍然是无偏的
因为 $u$ 虽然自身有一定趋势，但是和 $x$ 不相关，此时参数的估计仍然是无偏的（具体推导请查书，大体思想是是 $β$ 估计值等于 $β + f (x u)$ , $x u = 0 时 f (x u) = 0$ ，因此估计仍是无偏的）。具体几何解释可以理解为， $u$ 和 $x$ 仍然是正交的， $y$ 在 $x$ 空间的投影不会发生变化。
参数估计方差不再是最小的
之前最小二乘的估计具有无偏和方差最小的性质。 $β$ 估计值等于 $β + f (x u)$ ，当出现异方差时， $β$ 的方差估计会出现 $E(u_i^2)$ 项，由于 $u_i$ 的方差不再， $β$ 的方差估计不再是最小的（非有效的）。具体的非有效程度取决于异方差的程度。此时我们的t检验值会被高估，显著性被放大，此时的t检验是无效的！！！

1.4.自相关的检验（D-W检验）

对于自相关的检验，可画出残差 $e_t$ 和 $e_{t-1}$ 关系图，简单观察是否存在自相关，也可以使用自相关最经常使用的D-W检验。

D-W检验的方法如下（公式不是非常多，就手打吧。。。）：

首先D-W基于一系列的假设条件下，包括：

解释变量X是非随机的
随机误差 $u_t$ 满足一阶自回归模型 $u_t=ρu_{t-1}+ε_t$ ，且满足线性回归的古典假定
原始线性回归中没有滞后的内生变量作为解释变量（ $Y_t$ 的回归式中不包含 $Y_{t-1}$ 等滞后项）
截距项不为0（截距项在回归分析中非常的重要，拟合优度，带约束的检验等都要求带有截距项，如果不是必要的情况，都要加上截距项）
数据无缺失

满足上面条件后，就可以开产D-W检验，因为我们假定 $u_t=ρu_{t-1}+ε_t$ ，本质上我们检验的是 $ρ$ 是否等于0！因此我们的原假设为 $H_0: ρ=0$ 备择假设为 $H_1: ρ≠0$ 。 $ρ$ 就是 $u_t$ 与 $u_{t-1}$ 的相关系数。所以，我们先估计出原始多元线性回归 $Y_i=β_0+β_1*x_{1i}+β_2*x_{2i}+...+β_k*x_{ki}+u_i$ 的残差 $e_t$ ，并研究 $e_t$ 与 $e_t-1$ 的关系，具体的我们构造D-W检验统计量：
在这里插入图片描述
因为 $\sum_{t=2}^ne_t^2\approx\sum_{t=2}^ne_{t-1}^2$ ,所以有：
$d=\frac{2\overset n{\underset{t=2}{(\sum}}e_t^2-\sum_{t=2}^ne_te_{t-1})}{\sum_{t=2}^ne_t^2}=2(1-\frac{\sum_{t=2}^ne_te_{t-1}}{\sum_{t=2}^ne_t^2})$
因为 $ρ$ 的估计量为： $\widehat\rho=\frac{\sum_{t=2}^ne_te_{t-1}}{\sum_{t=2}^ne_t^2}$ ，所以实际上 $d≈2(1-\widehat\rho)$ 。

我们可以计算出D-W统计量 $d$ 的值，但d不服从正态分布、t分布、F分布、卡方分布中的任何一种。所以发明者建立一个 $d$ 统计量检验的上限临界值 $D_U$ 和下限临界值 $D_L$ 。根据我们计算的 $d$ 值，自变量个数 $k$ ，和如下对应关系判断是否存在自相关:

$0<d<D_L$ ，表明存在一阶正相关（ $\widehat\rho$ 接近于1时）
$D_L<d<D_U$ ，表明不确定是否存在自相关（ $\widehat\rho$ 接近于0时）
$4-D_U<d<4-D_L$ ，表明不确定是否存在自相关（ $\widehat\rho$ 接近于0时）
$4-D_L<d<4$ ，表明表明存在一阶负相关（ $\widehat\rho$ 接近于-1时）

1.5.自相关的补救——广义最小二乘
和自相关一样，都是因为 $Ω$ 和原始形态不一样，到时最小二乘估计不是最有效的，因此使用广义最小二乘估计。关于广义最小二乘，此处就不再重复了，具体请见用Python底层编写进行计量经济分析（三）：异方差（原因、结果、检验：White检验、补救：广义线性回归）

二、自回归检验和广义线性回归的python实现

我们沿用之前的框架，因为没有找到Python中查找D-W对应的临界值表的方法，此处值返回求出D-W统计量，可以自己根据求出的D-W统计量手动查表判断存在自相关。此处就直连带之前的代码全都写上了。

import numpy as np
import pandas as pd
import scipy.stats as st

class mul_linear_model(object):
    #intercept定义是否有截距项，在计量经济分析中，基本都要有截距项，否则计算出的拟合优度是没有意义的
    def __init__(self, y, X, intercept=True):
        data_x = X.copy()
        data_y = y.copy()
        self.data_x = data_x
        self.y = np.array(data_y)
        self.intercept = intercept
        
        if intercept == False:
            self.columns = data_x.columns
            self.X = np.mat(data_x)
        else:
            #插入截距项
            data_x.insert(0,'intercept',np.ones(len(data_x)))   
            #每个变量名称
            self.columns = data_x.columns
            #转换为矩阵，方便后面运算
            self.X = np.mat(data_x)
        
        #X的转置，方便后面运算
        self.XT=np.mat(self.X).T 
        self.N = self.X.shape[0]   #样本量
        self.K = self.X.shape[1]  #变量个数
    
    #拟合模型
    def fit(self, output=True):
        #使用最小二乘法(X'X)的逆成X'y
        self.b = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)),self.XT),np.mat(self.y).T))
        if output:
            for i in range(self.K):
                print("variable {0}'s cofe estimate is : {1}".format(self.columns[i], self.b[i][0]))
        
        self.S = np.array(np.linalg.inv(np.dot(self.XT,self.X)))
        return self.b
    
    #预测函数
    def predict(self, X_p):
        #直接使用条件期望作y=Xb为预测值
        y_hat =np.array(np.dot(X_p,self.b))
        return y_hat
    
    #参数T检验
    def T_test(self, alpha):
        #计算预测值
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算e'e,即残差平方和
        error_square = np.sum(np.square(y_hat-self.y))
        #计算扰动项的方差估计,即残差平方和
        sigma_hat = error_square/(self.N-self.K)
        #计算(X'X)的逆，用于估计参数b的方差
        #S = self.S
        for i in range(self.K):
            sk = self.S[i][i]
            #b的方差估计为sigma_hat * sk，构造t统计量，服从自由度为n-k的T分布
            Tk = self.b[i][0]/np.sqrt(sigma_hat * sk)
            #查表进行检验
            if abs(Tk) > abs(st.t.ppf(alpha/2, df = self.N-self.K)):
                print("T of variable {0} is: {1}, refuse H0".format(self.columns[i], Tk))
            else:
                print("T of variable {0} is: {1}, can't refuse H0".format(self.columns[i], Tk))
                
    #拟合优度计算
    def R_square(self):
        y_ba = np.mean(self.y)
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算总平方和
        self.SST = np.sum(np.square(self.y-y_ba))
        #计算误差平方和
        self.SSE = np.sum(np.square(self.y-y_hat))
        #计算离差平方和
        self.SSR = self.SST - self.SSE
        #利用拟合优度公式计算拟合优度
        self.R_2 = self.SSR/self.SST
        
        return self.R_2
    
    #整体F检验
    def F_test(self, alpha):
        #先求出拟合优度
        self.R_square()
        #利用公式求F统计量，附送自由优度为(K-1, N-K)的F分布
        F = (self.R_2/(self.K-1))/((1-self.R_2)/(self.N-self.K))
        #查表进行检验
        if abs(F) > abs(st.f.ppf(1-alpha, dfn = self.K-1, dfd = self.N-self.K)):
            print("F of function is: {0}, refuse H0".format(F))
        else:
            print("F of function is: {0}, can't refuse H0".format(F))
    
    def data_relation(self):
        if self.intercept:
            return self.data_x.loc[:,self.data_x.columns!='intercept'].corr()
        else:
            return self.data_x.corr()
        
        
    
    def VIF_cul(self, ind):
        #使用最小二乘法(X'X)的逆成X'y
        yk = np.array(self.X[:, ind].T)
        #有截距项时：
        if self.intercept:
            Xk = self.X[:, [i for i in range(1,ind)]+[j for j in range(ind+1, self.K)]]
        #无截距项时：
        else:
            Xk = self.X[:, [i for i in range(ind)]+[j for j in range(ind+1, self.K)]]
        
        vif_model = mul_linear_model(yk, pd.DataFrame(Xk), intercept=False)
        vif_model.fit(output=False)
        Rk_2 = vif_model.R_square()
        VIF_k = 1/(1-Rk_2)
        return VIF_k
        
    
    def VIF_test(self):
        for i in range(1, self.K):
            VIF_k = self.VIF_cul(i)
            print("variable {0}'s VIF estimate is : {1}".format(self.columns[i], VIF_k))
    
    def condition_num(self):
        if self.intercept:
            M = np.dot(self.XT[[i for i in range(1, self.K)],:], self.X[:, [j for j in range(1, self.K)]])
        else:
            M = np.dot(self.XT, self.X)
        eigenvalue, eigenvector = np.linalg.eig(M)
        cond_num = np.sqrt(eigenvalue[0]/eigenvalue[-1])
        print("data X's condition number is: {0}".format(cond_num))
    
    #岭回归
    def Ridge_fit(self, k = 0.1):
        #使用岭回归(X'X+kI)的逆成X'y
        self.b = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K))),self.XT),self.y).T)
        
        for i in range(self.K):
            print("variable {0}'s cofe estimate is: {1}".format(self.columns[i], self.b[i][0]))
        
        self.S = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K))), np.dot(self.XT,self.X)), np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K)))))
        
        return self.b
    
    #异方差怀特检验
    def White_test(self, alpha):
        #算出y_hat
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算偏差
        e = y_hat-self.y
        import matplotlib.pyplot as plt
        plt.plot(e)
        plt.title('ei')
        #计算white检验的自变量white_x,包括每一个x,和x的平方，以及两两之间的交叉项
        white_x = self.data_x.copy()
        x_columns = white_x.columns
        #计算平方项与交叉项，第一列是截距项，无需计算
        for i in range(1, len(x_columns)):
            for j in range(1, i+1):
                white_x[x_columns[i]+'*'+x_columns[j]] = white_x[x_columns[i]]*white_x[x_columns[j]]
        
        #调用自身类，进行回归求R2,white_x已经包含截距项，无需再添加截距项
        white_model = mul_linear_model(e, white_x, intercept=False)
        #求回归参数
        white_model.fit(output=False)
        #求拟合优度
        white_R_2 = white_model.R_square()
        #计算统计量自由度
        white_df = white_model.K-1
        #计算white统计量
        white = white_model.N*white_R_2
        #进行假设检验
        if white > st.chi2.ppf(1-alpha, df = white_df):
            print("white statistic is: {0}, refuse H0".format(white))
        else:
            print("white statistic is: {0}, can't refuse H0".format(white))
    
    #自相关DW检验
    def DW_test(self):
        #算出y_hat
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算偏差
        e = y_hat-self.y
        #得到当期残差e(t)
        e_t = e[1:len(e)]
        #得到前一期的残差e(t-1)
        e_tlag = e[0:len(e)-1]
        #进行回归求rou
        DW_model = mul_linear_model(e_t, pd.DataFrame(e_tlag), intercept=False)
        DW_model.fit(output=False)
        rou = DW_model.b[0][0]
        #利用dw统计量和rou的关系求出DW统计量
        d = 2*(1-rou)
        return d
    
    def gls_fit(self, output=True):
        #算出y_hat
        y_hat = np.array(np.dot(self.X,self.b).T)[0]
        #计算偏差
        e = y_hat-self.y
        #得到当期残差e(t)
        e_t = e[1:len(e)]
        #得到前一期的残差e(t-1)
        e_tlag = e[0:len(e)-1]
        #进行回归求rou
        DW_model = mul_linear_model(e_t, pd.DataFrame(e_tlag), intercept=False)
        DW_model.fit(output=False)
        rou = DW_model.b[0][0]
        
        #生成omiga矩阵过程
        relation_array=[]
        for i in range(len(self.y)):
            relation_array.append(pow(rou, i))
        
        omiga = [relation_array.copy()]
        
        for i in range(1, len(self.y)):
            relation_array = relation_array.copy()
            relation_array.insert(0, pow(rou, i))
            relation_array.pop()
            omiga.append(relation_array)
        #求出omiga矩阵
        omiga  = np.mat(omiga)
        #求出omiga矩阵的逆inv(omiga)
        omiga_inv = np.linalg.inv(omiga)
        
        #使用广义最小二乘法(X' inv(omiga) X)的逆成X'inv(omiga)y
        #np.dot(np.dot(np.dot(np.linalg.inv(np.dot(np.dot(self.XT,omiga_inv), self.X)), self.XT), omiga_inv), np.mat(self.y).T)
        self.b = np.array(np.dot(np.dot(np.dot(np.linalg.inv(np.dot(np.dot(self.XT,omiga_inv), self.X)), self.XT), omiga_inv), np.mat(self.y).T))
        
        if output:
            for i in range(self.K):
                print("variable {0}'s cofe estimate is : {1}".format(self.columns[i], self.b[i][0]))
        
        #需要注意的是，广义线性回归的方差会一起变
        self.S = np.array(np.linalg.inv(np.dot(np.dot(self.XT,omiga_inv), self.X)))
        return self.b

我们使用上一节的数据继续进行试验：

linear_model = mul_linear_model(house_data['value'],house_data[['HouseAge', 'AveRooms', 'AveBedrms', 'Population']])
linear_model.fit(output = False)
dw = linear_model.DW_test()
linear_model.gls_fit()

得到结果如下:
在这里插入图片描述
计算出的D-W统计量为0.912，参数估计如图。

到上面，多元线性回归的部分基本就写完了，如果代码中存在错误欢迎指正！

nbszg

关注

7
点赞
踩
66

收藏

觉得还不错? 一键收藏
1
评论
用Python底层编写进行计量经济分析（四）：自相关（原因、结果、检验：DW检验、补救：广义线性回归）

系列前面的文章：1.用Python底层编写进行计量经济分析（一）：多元线性回归（参数估计、T检验、拟合优度、F检验）2.用Python底层编写进行计量经济分析（二）：多重共线性（原因、结果、检验：条件数/方差膨胀因子、补救：岭回归）3用Python底层编写进行计量经济分析（三）：异方差（原因、结果、检验：White检验、补救：广义线性回归）多元线性回归的基本假定：模型符合线性模式XX...
复制链接

扫一扫