空间计量：地理加权回归(GWR)模型残差自相关检验

作者: stata连享会 | 来源:发表于2019-10-08 12:19 被阅读0次

作者：陈凤 (西安交通大学)

Stata连享会计量专题 || 精品课程 || 简书推文 || 公众号合集

点击查看完整推文列表

2020寒假Stata现场班(北京, 1月8-17日，连玉君-江艇主讲)

2020寒假Stata现场班

@[toc]

1. 简介

当一般线性模型的扰动项存在自相关时，即：
${\rm E} (\varepsilon_i,\varepsilon_j | \boldsymbol X) \neq 0$ 时, 虽然回归系数的OLS估计值依然是无偏且一致的，但是由于系数估计值的方差 ${\rm Var}(\boldsymbol \beta|\boldsymbol X)$ 不在是 $\sigma^2(\boldsymbol X^{\rm T}\boldsymbol X)^{-1}$ , 故常用的 $t$ 检验和 $F$ 检验也会相应的失效[陈强,2010]。此外，由于扰动项的自相关违背了高斯-马尔科夫定理，因此OLS也不再是BLUE。

图1：扰动项自相关对系数估计的影响
如上图1所示，实线表示真实的回归直线，扰动项之间存在正自相关。由该图可知由于扰动项自相关的存在，使得数据在回归直线上下波动的幅度增加，最终使得回归系数估计的精度降低。同理，在GWR模型中，一个重要的假设条件是误差项是独立同分布的。但是对于空间数据来说，空间自相关性是其重要的特征之一。因此，对空间数据建立GWR模型时，有可能误差项会存在一定的自相关性，从而如一般线性模型扰动项自相关一样对参数估计和检验等产生一定的影响。综上所述，在实际应用中，有必要对模型扰动项的自相关性进行检验。

2. GWR模型残差自相关性检验

2.1 GWR模型的残差

GWR模型如下：
$y_i=\beta_0(u_i,v_i)+\sum_{j=1}^p\beta_j(u_i,v_i)x_{ij}+\varepsilon_i, \tag{1}$ 其中， $\beta_j(u,v)\ (j=0,1,\cdots,p)$ 为空间地理位置函数。令 $\boldsymbol Y=(y_1,y_2,\cdots,y_n)^{\rm T}$ , 采用GWR估计方法对模型(1)进行估计，得到 $\boldsymbol Y$ 和 $\boldsymbol \varepsilon$ 的估计值 $\boldsymbol{ \hat{Y}}$ 和 $\boldsymbol{\hat{\varepsilon}}$ 如下：
$\boldsymbol{ \hat{Y}}=\boldsymbol L \boldsymbol Y, \tag{2}$
和 $\boldsymbol{\hat{\varepsilon}}=(\hat\varepsilon_1,\hat\varepsilon_2,\cdots,\hat\varepsilon_n)^{\rm T}=(\boldsymbol I-\boldsymbol L)\boldsymbol Y, \tag{3}$
其中， $\boldsymbol X=(\boldsymbol x_1,\boldsymbol x_2,\cdots,\boldsymbol x_p)$ , $\boldsymbol{x}_{i}=(x_{1j},x_{2j},\cdots,x_{nj}),$ $\boldsymbol{\tilde{x}}_{i}^{\rm T}=(x_{i1},x_{i2},\cdots,x_{ip}),$
$\boldsymbol L(h)= \left( \begin{array}{cccc} \boldsymbol{ \tilde{x}}_{1}^{\rm T}\left(\boldsymbol X^{\rm T}\boldsymbol W_{h}(u_1,v_1)\boldsymbol X\right)^{-1}\boldsymbol X^{\rm T} \boldsymbol W_{h}(u_1,v_1)\\ \boldsymbol{ \tilde{x}}_{2}^{\rm T}\left(\boldsymbol X^{\rm T}\boldsymbol W_{h}(u_2,v_2)\boldsymbol X\right)^{-1}\boldsymbol X^{\rm T} \boldsymbol W_{h}(u_2,v_2)\\ \vdots\\ \boldsymbol{ \tilde{x}}_{n}^{\rm T}\left(\boldsymbol X^{\rm T}\boldsymbol W_{h}(u_n,v_n)\boldsymbol X\right)^{-1}\boldsymbol X^{\rm T} \boldsymbol W_{h}(u_n,v_n) \end{array} \right),$ 和 $\boldsymbol W_h(u_i,v_i)={\rm Diag}\left(K(d_{i1}/h),K(d_{i2}/h),\cdots,K(d_{in}/h)\right),$ $K(\cdot)$ 为核函数, $d_{ij}$ 为点 $(u_i,v_i)$ 和 $(u_j,v_j)$ 之间的欧式距离, $h$ 为窗宽. 最优窗口可以采用 ${\rm AIC_c}$ 、 ${\rm CV}$ 和 ${\rm GCV}$ 等准则来选取。

2.2 Moran's I 检验空间自相关性

记 $\boldsymbol{\hat{\varepsilon}}=(\hat\varepsilon_1,\hat\varepsilon_2,\cdots,\hat\varepsilon_n)^{\rm T}$ 为GWR模型的残差， $\boldsymbol W=(w_{ij})$ 为已知的空间权重矩阵，则Moran's I 的表达式如下：
$I_0=\frac{n}{s}\frac{\sum_{i=1}^n\sum_{j=1}^nw_{ij}\hat{\varepsilon}_i\hat{\varepsilon}_j}{\sum_{i=1}^n\hat{\varepsilon}^2_i}=\frac{n}{s}\frac{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol W\boldsymbol{\hat{\varepsilon}}}{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol{\hat{\varepsilon}}}, \tag{4}$ 其中， $s=\sum_{i=1}^n\sum_{j=1}^nw_{ij},$ $\boldsymbol W$ 一般情况下为行标准化的空间全矩阵并且主对角线上的元素为0. 如果 $I_0$ 显著的大于0，则说明残差存在正的空间自相关性; $I_0$ 显著的小于0，则说明残差存在负的空间自相关性。当 $\boldsymbol W$ 为非对称矩阵时，可以构造新的对称矩阵如下：
$\boldsymbol W^{\ast}=\frac{1}{2}(\boldsymbol W^{\rm T}+\boldsymbol W),$
并且有
$\frac{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol W\boldsymbol{\hat{\varepsilon}}}{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol{\hat{\varepsilon}}}=\frac{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol W^{\ast}\boldsymbol{\hat{\varepsilon}}}{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol{\hat{\varepsilon}}}，$
因此，不失一般性，假设 $\boldsymbol W$ 为一个对称矩阵。另外，由于 $\frac{n}{s}$ 不会对Moran's I 的检验 $p-$ 值产生影响。故 Moran's I 的计算如下：
$I_0=\frac{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol W\boldsymbol{\hat{\varepsilon}}}{\boldsymbol{\hat{\varepsilon}}^{\rm T}\boldsymbol{\hat{\varepsilon}}}. \tag{5}$

2.3 $\chi^2$ 三阶矩近似 $p-$ 值

在一般的线性模型中，Moran's I 的零假设分布可以近似为正态分布，但是由于矩阵 $(\boldsymbol I-\boldsymbol L)$ 在GWR模型中一般不再是幂等矩阵，所以很难去研究 $I_0$ 的正态性。针对这一情况，Leung, Y等(2000)给出了近似计算Moran's I 检验 $p-$ 值的 $\chi^2$ 三阶矩。具体计算如下：
令 $Q={\boldsymbol \varepsilon}^{\rm T}(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)\boldsymbol \varepsilon,$ $r$ 为 Moran's I 的观测值。
(1)当 ${\rm E}[Q-{\rm E}(Q)]^3>0$ 时，
$p(I_0\leq r)=p(Q \leq 0)\approx p\left \{\chi_h^2 \leq h-\frac{\sqrt{2h}{\rm E}(Q)}{\sqrt{{\rm var(Q)}}}\right\} \tag{6}$
(2)当 ${\rm E}[Q-{\rm E}(Q)]^3<0$ 时，
$p(I_0\leq r)\approx 1-p\left \{\chi_h^2 \leq h+\frac{\sqrt{2h}{\rm E}(Q)}{\sqrt{{\rm var(Q)}}}\right\} \tag{7}$
其中， $h=\frac{8[{\rm var}(Q)]^3}{\{{\rm E}[Q-{\rm E}(Q)]^3\}^2}=\frac{\{ tr[(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)]^2\}^3}{\{ tr[(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)]^3\}^2},$
${\rm E}(Q)=tr[(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)],$ ${\rm var}(Q)=2tr[(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)]^2,$ ${{\rm E}[Q-{\rm E}(Q)]^3}=8tr[(\boldsymbol I-\boldsymbol L)^{\rm T}(\boldsymbol W-r \boldsymbol I)(\boldsymbol I-\boldsymbol L)]^3,$ 和 $\chi^2_h$ 为自由度是 $h$ 的一个 $\chi^2$ 分布。当得到的Moran's I 值小于等于0的时候，对应的检验 $p-$ 值为 $p(I_0 \leq r)$ ；当得到的Moran's I 值大于0的时候，对应的检验 $p-$ 值为 $p(I_0 \geq r)$ 。

连享会计量方法专题……

3. Matlab代码实现

1. 调用函数GWR_estimation进行GWR参数估计：
输入数据自变量Y，因变量X以及光滑参数的取值范围，得到回归系数估计值Parameter，残差Error，帽子矩阵L，以及最优窗宽bestbandwidth。该函数在估计过程中使用核函数为高斯函数，并采用$\rm AIC_c$准则选取最优窗宽。
function [Parameter,Error,L,bestbandwidth]=GWR_estimation(Y,X,Smooth)
        [p,n]=length(X);
        Parameter=zeros(n,p);
        W_kernel=zeros(n,n);%%计算权矩阵(核函数)
        L_h=zeros(n,n);
        E_error=zeros(1,length(Smooth));
        AIC=zeros(1,length(Smooth));
        for k=1:length(Smooth)
            for i=1:n
                for j=1:n
                    W_kernel(j,j)=1/sqrt(2*pi)*exp(-(Distance(i,j))^2/(2*Smooth(1,k)^2));
                end
                L_h(i,:)=X(i,:)*inv(X'*W_kernel*X)*X'*W_kernel;
            end
            E_error(1,k)=(Y(:,1)'*(eye(n,n)-L_h)'*(eye(n,n)-L_h)*Y(:,1))/n;
            AIC(1,k)=log(E_error(1,k))+(n+trace(L_h))/(n-2-trace(L_h));
        end
        index_Smooth=min(AIC);
        bestbandwidth=Smooth(1,index_Smooth);
        L=zeros(n,n);
        for i=1:n
            for j=1:n
                W_kernel(j,j)=1/sqrt(2*pi)*exp(-(Distance(i,j))^2/(2*(Smooth(1,index_Smooth))^2));
            end
            L(i,:)=X(i,:)*inv(X'*W_kernel*X)*X'*W_kernel;
            Parameter(i,:)=(inv(X'*W_kernel*X)*X'*W_kernel*Y(:,1))';
        end
        Error=(eye(n)-L)*Y(:,1);
end
2. 调用Moran_test函数计算GWR模型残差的Moran's I 值以及p-值：
输入空间权矩阵W，残差Error和帽子矩阵L，调用Moran_test，可以得到残差的Moran's I 值MoranI以及对应的检验p-值P。
function [MoranI,P]=Moran_test(W,Error,L)
W3=1/2*(W+W');
MoranI=Error'*W3*Error/(Error'*Error);%计算的Moran I值
n=length(Error);
N=eye(n)-L;
T=N'*(W3-MoranI.*eye(n))*N;
Q_expecation=trace(T);%E(Q)
Q_Var=2*trace(T^2);%Var(Q)
Q_three=8*trace(T^3);%E[Q-E(Q)]^3
h=8*(Q_Var)^3/((Q_three)^2);%卡方分布的自由度
v=(sqrt(2*h)*Q_expecation)/(sqrt(Q_Var));
if Q_three>0
    P=chi2cdf(h-v,h);
elseif Q_three<0
    P=1-chi2cdf(h+v,h);
end
if MoranI<=0
    P=P;
elseif MoranI>0
    P=1-P;
end
end

连享会计量方法专题……

4. 参考文献：

陈强. 高级计量经济学及Stata应用[M]. 高等教育出版社, 2010.
Leung, Y., Mei, C.-L., & Zhang, W.-X. (2000). Testing for Spatial Autocorrelation among the Residuals of the Geographically Weighted Regression. Environment and Planning A: Economy and Space, 32(5), 871–890.[[PDF]] (https://journals.sagepub.com/doi/10.1068/a32117#articleCitationDownloadContainer)

关于我们

「Stata 连享会」 由中山大学连玉君老师团队创办，定期分享实证分析经验，公众号：StataChina。
公众号推文同步发布于 CSDN 、简书和知乎Stata专栏。可在百度中搜索关键词「Stata连享会」查看往期推文。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
欢迎赐稿： 欢迎赐稿。录用稿件达三篇以上，即可免费获得一期 Stata 现场培训资格。
E-mail： StataChina@163.com
往期推文：计量专题 || 精品课程 || 简书推文 || 公众号合集

点击此处-查看完整推文列表

在这里插入图片描述

空间计量：地理加权回归(GWR)模型残差自相关检验

2020寒假Stata现场班(北京, 1月8-17日，连玉君-江艇主讲)

1. 简介