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

我不是女神ヾ 2023-05-21 15:18 6阅读 0赞

系列前面的文章：

1.  [多元线性回归和显著性检验（参数估计、T检验、F检验、拟合优度）][T_F]
2.  [多重共线性（导致结果、检验——方差膨胀因子、补救措施——岭回归）][Link 1]
3.  [异方差（导致结果、检验——White、补救措施——广义线性回归）][White]
4.  [自相关（导致结果、检验——D-W、补救措施——广义线性回归）][D-W]

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

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

1.  模型符合线性模式
2.   X X X满秩（无多重共线）
3.  零均值价值： E ( ε i ∣ X i ) = 0 E(ε\_i|X\_i)=0 E(εi∣Xi)=0（自变量外生）
4.  同方差： V a r ( ε i ∣ X i ) = σ Var(ε\_i|X\_i)=σ Var(εi∣Xi)=σ
5.  无自相关： c o v ( ε i , ε j ) = 0 cov(ε\_i, ε\_j)=0 cov(εi,εj)=0
6.  球形扰动：ε\_i是正态分布

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

#### 一、异方差性 ####

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

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

**1.2.自相关的原因**

1.  变量本身带有趋势  
    我们通常研究的时间序列自身就是带有一定趋势的，比如 y t y\_t yt和 y t − 1 y\_\{t-1\} yt−1本身就存在一定关系（此时可试试AR模型），此时就会出现异方差问题。
2.  变量的之后效应  
    通常，一个变量对另一个变量的影响并不会在当期表现出来，往往需要之后一期或几期才会表现出来，如果忽视了这些之后效应，模型也可能出现自相关问题。
3.  模型设定偏差  
    如果我们的模型设定存在偏差，比如真实模型为：  
     Y i = β 0 + β 1 ∗ x 1 i + β 2 ∗ x 2 i + β 3 ∗ x 3 i 2 + u i Y\_i=β\_0+β\_1\*x\_\{1i\}+β\_2\*x\_\{2i\}+β\_3\*x\_\{3i\}^2+u\_i Yi=β0\+β1∗x1i\+β2∗x2i\+β3∗x3i2\+ui  
    我们建模时，建立的模型为  
     Y i = β 0 + β 1 ∗ x 1 i + β 2 ∗ x 2 i + + v i Y\_i=β\_0+β\_1\*x\_\{1i\}+β\_2\*x\_\{2i\}++v\_i Yi=β0\+β1∗x1i\+β2∗x2i\+\+vi  
    此时 v i = β 3 ∗ x 3 i 2 + u i v\_i=β\_3\*x\_\{3i\}^2+u\_i vi=β3∗x3i2\+ui， v i v\_i vi包含了 x 3 i 2 x\_\{3i\}^2 x3i2的趋势。此时可能出现自相关

**1.2.自相关的后果**

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

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

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

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

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

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

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

满足上面条件后，就可以开产D-W检验，因为我们假定 u t = ρ u t − 1 + ε t u\_t=ρu\_\{t-1\}+ε\_t ut=ρut−1\+εt，本质上我们检验的是 ρ ρ ρ是否等于0！因此我们的原假设为 H 0 : ρ = 0 H\_0: ρ=0 H0:ρ=0备择假设为 H 1 : ρ ≠ 0 H\_1: ρ≠0 H1:ρ=0。 ρ ρ ρ就是 u t u\_t ut与 u t − 1 u\_\{t-1\} ut−1的相关系数。所以，我们先估计出原始多元线性回归 Y i = β 0 + β 1 ∗ x 1 i + β 2 ∗ x 2 i + . . . + β k ∗ x k i + u i Y\_i=β\_0+β\_1\*x\_\{1i\}+β\_2\*x\_\{2i\}+...+β\_k\*x\_\{ki\}+u\_i Yi=β0\+β1∗x1i\+β2∗x2i\+...\+βk∗xki\+ui的残差 e t e\_t et，并研究 e t e\_t et与 e t − 1 e\_t-1 et−1的关系，具体的我们构造D-W检验统计量：  
![在这里插入图片描述][20200417160852458.png]  
因为 ∑ t = 2 n e t 2 ≈ ∑ t = 2 n e t − 1 2 \\sum\_\{t=2\}^ne\_t^2\\approx\\sum\_\{t=2\}^ne\_\{t-1\}^2 ∑t=2net2≈∑t=2net−12,所以有：  
 d = 2 ( ∑ t = 2 n e t 2 − ∑ t = 2 n e t e t − 1 ) ∑ t = 2 n e t 2 = 2 ( 1 − ∑ t = 2 n e t e t − 1 ∑ t = 2 n e t 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\}) d=∑t=2net22t=2(∑net2−∑t=2netet−1)=2(1−∑t=2net2∑t=2netet−1)  
因为 ρ ρ ρ的估计量为： ρ ^ = ∑ t = 2 n e t e t − 1 ∑ t = 2 n e t 2 \\widehat\\rho=\\frac\{\\sum\_\{t=2\}^ne\_te\_\{t-1\}\}\{\\sum\_\{t=2\}^ne\_t^2\} ρ=∑t=2net2∑t=2netet−1，所以实际上 d ≈ 2 ( 1 − ρ ^ ) d≈2(1-\\widehat\\rho) d≈2(1−ρ)。

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

1.   0 < d < D L 0<d<D\_L 0<d<DL，表明存在一阶正相关（ ρ ^ \\widehat\\rho ρ接近于1时）
2.   D L < d < D U D\_L<d<D\_U DL<d<DU，表明不确定是否存在自相关（ ρ ^ \\widehat\\rho ρ接近于0时）
3.   4 − D U < d < 4 − D L 4-D\_U<d<4-D\_L 4−DU<d<4−DL，表明不确定是否存在自相关（ ρ ^ \\widehat\\rho ρ接近于0时）
4.   4 − D L < d < 4 4-D\_L<d<4 4−DL<d<4，表明表明存在一阶负相关（ ρ ^ \\widehat\\rho ρ接近于-1时）

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

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

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

[T_F]: https://blog.csdn.net/u011517132/article/details/105286113
[Link 1]: https://blog.csdn.net/u011517132/article/details/105323135
[White]: https://blog.csdn.net/u011517132/article/details/105392242
[D-W]: https://blog.csdn.net/u011517132/article/details/105581371
[20200417160852458.png]: /images/20230521/367109beba7c4a4ebf3ec2d8fbde0ce3.png
[20200417185619691.png]: /images/20230521/56d7b57dbb024f0286e00abc580a75e2.png