绝对实用! 空间计量方法使用指南(附代码和实例)
稿件：econometrics666@126.com
所有计量经济圈方法论丛的code程序, 宏微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.
“空间计量研究小组”将为引荐空间计量全套方法和相关说明。之前，已经推荐过空间计量方面introductory性质的文章“空间计量模型选择、估计、权重、检验(Spatial effect)”和“空间计量经济学最新进展和理论框架，comprehensive track”。
空间计量什么鬼？
空间计量也叫做spatial econometrics，实际上就是把咱们平时用到的那些方法加入一些空间效应后做的系列回归。空间效应，实质上就是一种网络效应，证明了万事万物之间都有关联性，当然，越是靠得近的两个objects也更可能有强关联性。要不然，你真的会相信，一只拉丁美洲的蝴蝶振动几下翅膀，然后中国就必定发生地震或者洪灾了。再考虑下，当年美国准备attack North Korea，因为我们知道“唇亡齿寒”的道理，所以中国派遣了那么多的志愿军跨过鸭绿江去援助老金家。现在的物理距离可能还不如经济距离实在，比如，我们与韩国是那么近距离的邻居，可是在经济军事上韩国与美国可能反而走得更近(如果不是如此，请忽略该陈述)。
举二个例子来看看空间计量在经济事务中的运用。北京市的发展会通过neighbourhood效应影响河北和天津的发展，比如，导致河北和天津发展相对滞后一些(大众直觉)。如果想要研究政府出台的“限购政策”对房价的影响，那么，我们可以搜集全国出台限购政策的城市的几年面板数据，然后做一个普通的xtreg回归就行。然而，我们要晓得，各个城市出台的限购政策不仅会影响当地城市的房价，还会通过诸如“人口流动”等影响另一个城市的房价，而且，这个城市的房价还可能直接影响另一个城市的房价。所以，我们就需要考虑这样一种剪不断、理还乱的网状关系，让限购政策对房价的影响分成两部分：直接影响和间接影响。
空间计量经济学可以处理截面数据，也可以处理面板数数据，而且还可以处理变量内生性问题(目前，处理类型是受限的)。空间计量属于比较新发展起来的学科，很多计量理论和方法都还在探索中，因此，下面的部分主要集中于几个用得比较多的空间计量模型。
空间计量回归数据
对于空间计量回归，我们首先需要找到的是shapefile，一种可以装着各种经纬度等地理信息的file包。当我们拿到一个shapefile后，就需要通过spshape2dta转换成我们当前软件可以识别的数据格式。我们要晓得，这里面只有一些地理信息相关的数据，因此我们首先需要生成一个地方唯一的标识——它的ID，这样才能与之后咱们要关联的其他经济社会数据进行合并使用。比如，我们想通过广东的县域数据了解平均受教育水平与经济发展之间的关系。那么，我们需要把广东各个县域的shapefile找到，然后让每个县都有一个唯一的ID去标识(一般根据城市+县编制)。然后，我们再把广东县域的经济社会统计年鉴数据库与上面的shapefile根据唯一的ID进行一一合并。最终得到的数据库就可以用来作为空间计量分析的样本了。
以下就是一个shapefile数据，里面只有这种各个地方的经纬度和名称的地理信息，因此，我们首先需要把这些信息根据需要转换成可以使用的数据(前面有提到)。
下面，我们直接通过一些例子来理解今天将要引荐的空间计量方法。
截面空间计量回归
紧接着的这个是截面数据，_ID是识别这个观测值的唯一identity，_CX, _CY是shapefile里该点的经纬度，cname, sname分别为县名和州名。当_ID=8时，这个点就代表了Texas州Menard县的观测值，所以你可以把把这些想象成习以为常的”某个人“。所以，我们也可以用空间计量来研究一下邻里效应或者网络效应或同伴效应，毕竟这些模型里都彼此存在着溢出效应。
前面，我们已经说过，需要通过唯一的ID(下面我们用的是fips)去标识一个县(州的code+县的code)。然后，再通过这个ID把我们感兴趣的数据库与shapefile进行一一合并。下面这个数据库就是合并之后的截面数据库(去掉了与咱们主题不相关的其他变量)。在紧接着的例子里，咱们只保留了三个变量：college(县里的大学毕业率), income(县的平均收入), unemployment(县的失业率)。接下来，我们想要研究的是，县的大学毕业率与县的失业率之间的关系，即是否大学毕业率降低了失业率水平——县里的大学生越多会降低失业率水平。
既然我们想要研究的是失业率水平，那么我们索性把数据库里所有县的失业率水平通过地图展示出来(颜色越深代表失业率越高)。
下面这张图是各个县的大学毕业率水平(颜色越深代表大学毕业率越高)，你可以看看与上面那张图片在颜色层次上重复度高不高？
我们首先做一个普通的OLS回归，把空间相关性忽视掉。从这个结果中，我们可以看出，当控制了收入后，大学毕业率确实降低了一个地方的失业率(college系数为负)。
我们现在想看看，这个回归中到底存不存在空间相关性呢？因此，我们考虑用Moran test去检验一下。检验结果表明，在我们的这个回归中存在很强的空间效应，因为它拒绝了H0: error is iid。
在做Moran test之前，实际上我们需要先创建一个空间权重矩阵spatial weighting matrix。SWM的创建方法很多种，但是我们常用的是二种方式：近邻空间权重矩阵(contiguity)和距离倒数空间权重矩阵(inverse distance)(翻译的种类不同，自己理解就行)。下面这个是咱们用近邻方式创建的SWM，可以看到的是，这个SWM里面大部分都是为0，表明二个县之间不存在邻里关系，因此彼此的变量不影响另一方。这些都可以自动让程序去执行，下方这个是输出的SWM结果。
空间计量回归有二种algorithms，他们分别是广义空间二阶段最小二乘法(gs2sls)和极大似然估计法(ml)。两者的主要区别：后者在error服从正态分布的时候更加有效，否则就不如gs2sls那么robust(比如，在有异方差情形下)。gs2sls实质上就是咱们习以为常的GMM估计，找到各种矩条件然后进行优化操作寻找到arg minimization时的估计值。
下面是我们用gs2sls做的空间自回归模型，即我们把其他县的因变量unemployment作为这个县的自变量放进模型，这有点像咱们平时遇到的Autoregressive regression (AR模型)。画黄线的wald检验表明，咱们的空间效应确实存在的，空间项是显著的，因此应该使用spatial回归。
下面是我们使用极大似然估计得到的空间自回归模型。如果直接比较这两个方法所得到的各个系数，那么我们还是能够看出一些差别的。
下面这个是我们把空间计量中常用的三个模型类别全部揉在一起的结果。上面，我们只用到了空间自回归模型，即把其他县的unemployment放在咱们县作为自变量进行的回归。下面这个允许三种形式的空间依存(实质上分别对应着三种不同类型的空间计量模型)：1. 其他县的大学毕业率college(spatial lag of independent variables)，2. 其他县的失业率unemployment(spatial lag of dependent variable)，即空间自相关模型， 3. error的相关性，即空间误差模型(spatial lag of error)，这有点像咱们习以为常的moving average模型(MA model)。
注意，我们一般不是通过上面的系数去汇报结果，毕竟我们压根解读不出来直接效应和间接效应，原因在于这会存在recursive问题。想一想，A影响B，B影响C，A也影响C，C影响B，C也影响A，我们就不能直接分清这里面的各种效应大小。因此，我们通过下面这张图获得各个变量的平均效应大小。比如，大学毕业率每上升１个百分点，那么失业率将会下降0.16%，而大学毕业率的间接效应竟然占据了大半(自己算一下，0.1 divides 0.16)，所以我们不得不考虑空间溢出效应问题。
在上面的回归模型基础上，下面我们稍微改变一下回归形式，即增加另一个空间权重矩阵SWM(两个黄颜色的点代表了二个不同的空间权重矩阵)。W为近邻空间权重矩阵(contiguity matrix)，M为距离倒数空间权重矩阵(inverse distance matrix)。这就是说，在空间回归模型中，我们还可以同时添加不同的空间权重矩阵，再用这些空间权重矩阵分别去做空间自回归或空间误差回归等模型。
通过获得平均效应，我们可以发现，现在，咱们的大学毕业率每上升1%，那么失业率会下降0.20%，要比前面的那个结果的效应大一些。具体要怎么选择，我们认为在研究中可以把这两个都放上，反正最后得出的效果是负数，符合你研究的hypothesis。
空间计量内生性问题
我们紧接着讨论一下内生性问题，也将通过一个例子来看看空间计量中如何去处理内生性问题。我们想要研究一下，与饮酒相关的被捕率是否会受到一个地方警察人数的影响——一个地方警察人数多可能会导致与饮酒相关的被捕率上升或下降。下方是美国南部这些县的与饮酒相关的被捕率分布。
与饮酒相关的被捕率(dui)相关的变量有警察人数(police)，其他种类的被捕率(nondui)，交通工具数量(vehicles)，是否禁酒(dry)。但是，我们要晓得，与饮酒相关的被捕率是与警察人数存在双向因果关系的，之前的被捕率高的话当然应当增加警察人数，即police from (to) dui。所以，我们通过工具变量“是否有县选举”作为police的工具变量(符合相关性和外生性条件)。
通过平均效应图，我们可以发现，通过工具变量处理了内生性问题后，警察人数的增加降低一个地方与饮酒相关的被捕率。这有力的证明了，一个地方要想降低逮捕率需要增加警察人数。这个结论好像与咱们的直觉有点相左，但这与中国现实是相符合的。在中国，越是警察人数多的地方，相对而言，与饮酒相关的被捕率要低很多(你可以不认同)。
面板空间计量回归
下面这个是一个面板数据，panel为_ID编号，year为1960-1990四个时期，因此(_ID=876, year=1960)表示West Virginia州Hancock县在1960年的观测值。面板数据的好处，请看看这篇文章“面板数据是怎样处理内生性的,一篇让人豁然明朗的文章”，主要是可以帮助我们识别出一个观测值的趋势，更好地处理那些不可观测因素的混淆效应。
我们下面想要研究县里的gini系数(gini)是否直接影响着县里的自杀率(hrate)，毕竟一个地方贫富差距过大在直觉上就会导致人们心理出现问题(比如，抑郁)。当然，我们把县的人口密度(ln_density)和人口规模(ln_population)作为控制变量添加进去，因为他们会影响一个地方的自杀率。
我们先做一个普通的随机效应回归，从下面可以看出，gini系数的增加会导致一个地方的自杀率上升(gini前面的系数为正)。把年份效应加入进去之后，我们发现相对于1960年，1970、1980、1990年的自杀率都要高一些。sigma_u: panel-level standard deviation; sigma_e: standard deviation of epsilon_it。
下面是咱们通过近邻权重矩阵创建的空间权重矩阵SWM(WX)，注意，这里我们使用的year=1990年里所有县的数据创建的SWM，因为我们认为在1990年是邻居的县在1960、1970和1980年里应当还是邻居(当然有例外，需要把例外找到)。
下面是我们使用空间回归模型来做一下，里面包括了spatial lag of dependent variable (hrate)和spatial lag of error (e.hrate)(如果你不清楚这两处的英文意思，请看文章的前面部分的解释)。我们仍然得到了相似的结论，即gini系数的增大会导致自杀率上升。
我们之前说过，由于上面的那个系数没有告诉我们gini的直接效应和间接效应，因此需要通过平均效应这个表格来读出相关信息。从下面的图片，我们可以看出，gini系数每上升1个百分点，那么县里每10万人中的自杀人数会增加46.85人。此时的直接效应与间接效应差不多，表明其他县的自杀率实际上会对我们这个县的自杀率产生影响的(这个需要引起重视)。至于具体的传递机制，那是“公说公有理，婆说婆有理”，“可到底谁是正确的”好像不是我们关心的。
在前面，我们发现1970、1980和1990三个时期的自杀率都比1960年的要高，从而让我们对于“是否随着时间推移gini系数对自杀率的影响增强”话题产生兴趣了。下面这张图片显示，我们用year对gini系数做了一个交叉项，以1960年作为基期后，发现似乎gini系数对hrate的影响与日剧增。
我们通过下图看出，year与gini的交叉项对自杀率的影响是显著的。
下面，我们通过四张按照年份展示的平均效应图，来看看gini系数对自杀率的影响是否在逐渐增强。
紧接着的是1960年的gini对自杀率的平均效应，它是54.9。
下面是1970年gini系数对自杀率的平均效应，它是75.4。
下面是1980年gini系数对自杀率的平均效应，它是131.5。
下面是1990年gini系数对自杀率的平均效应，它是231.8。
通过上面四张平均效应图，我们可以得到：gini系数对自杀率的影响是与日剧增的。
当errorlag出现时,添加了sarpanel允许个体效应(在这里是县域效应)出现像误差项那样的空间回归形式，即我们常说的ai也出现像error那用的空间回归形式。在这个数据集里，整体而言，当我们把个体效应也放入空间回归中，得到的结果与前面没有考虑个体效应的空间回归也差不多。
下面展示的是空间面板数据的固定效应回归结果，fixed effect estimator与random effect estimator估计方法是有差别的(这里不详述)。我们前面运行的re估计采用的是极大似然估计方法，因此运行所花费的时间就多，而当我们采用fe时运行所花费的时间就短很多(这与使用的算法有关)。从中我们可以得到的结论是，gini系数的上升会导致自杀率的升高，这个与前面的回归结果是一致的。
通过看下面平均效应图结果，我们发现gini系数对自杀率的总效应是37.7。这表明gini系数每上升1个百分点，那么县里每10万人中的自杀人数会增加37.7人。这比通过随机空间面板效应方法得到的效应要小一些。
**written by @空间计量研究小组(spatial econometrics research team)，作者: 甲壳虫
**-------------------------------------------------
**估计大学生毕业率对失业率的影响
//------------------------------------------------
**下载并处理成空间计量形式的数据----------
use texas_ue.dta, clear //载入这个数据
describe //描述一下这个数据
merge 1:1 fips using tl_2017_us_county //把前面的数据与shapefile合并
keep if _merge==3 //保留匹配好了的样本
drop _merge
rename NAME countyname
drop STATEFP COUNTYFP COUNTYNS GEOID
drop NAMELSAD LSAD CLASSFP MTFCC CSAFP
drop CBSAFP METDIVFP FUNCSTAT
drop ALAND AWATER INTPTLAT INTPTLON
save texas_ue.dta, replace //生成新的包含shapefile的样本
**数据整理好之后的空间计量估计-------------
summarize unemployment
grmap, activate
grmap unemployment //画出unemployment的地图
regress unemployment college income //做一个ols回归
spmatrix create contiguity W //生成相邻权重矩阵W
estat moran, errorlag(W) //检测一下是否有spatial correlation
spmatrix export W using contig.xls //把空间权重矩阵W输出
spmatrix create idistance M //生成距离倒数的权重矩阵M
**---------------------------------------------
**spregress方法考虑空间效应
//-----------------------------------------------
**如果没有加入iv,dv,error spatial lag做回归------------
regress unemployment college income
spregress unemployment college income, gs2sls
spregress unemployment college income, ml
**1.spregress有两种估计方法, spatial lag of the DV-------
spregress unemployment college income, gs2sls dvarlag(W) //用广义空间2sls估计
spregress unemployment college income, ml dvarlag(W) //用ml估计spatial lag of the dependent variable
**2.spatial lag of the IV-------------
spregress unemployment college income, gs2sls ivarlag(W:college)
estat impact //平均的效应
**3.spatial lag of error------------
spregress unemployment college income, gs2sls errorlag(W)
estat impact
**4.spatial lag of DV, IV and error--------
spregress unemployment college income, gs2sls dvarlag(W) heteroskedastic ///
errorlag(W) ivarlag(W: college) //同时把dv,iv和error的spatial都放进去
estat impact
spregress unemployment college income, gs2sls dvarlag(W) errorlag(W) ///
ivarlag(M: college) ///用了两个空间权重矩阵
estat impact
spregress unemployment college income, ml vce(robust) dvarlag(W) errorlag(W) ///
ivarlag(W: college) ivarlag(M: college) //ml方法可以有两个spatial lag of iv
estat impact
**----------------------------------------------
** spivregress方法处理内生性问题
//----------------------------------------------
use dui_southern, clear //用数据库来运行spivregress
spset //设置成空间回归形式的数据
**检测是否有sparital correlation-------
regress dui nondui vehicles i.dry
spmatrix create contiguity WI //生成相邻权重矩阵WI
spmatrix create idistance MI //生成距离倒数权重矩阵MI
estat moran, errorlag(WI) //检测是否有spatial correlation
estat moran, errorlag(WI) errorlag(MI)
grmap, activate
grmap dui //酗酒导致的被捕率(10wan)地图
spivregress dui nondui vehicles i.dry (police = elect), dvarlag(WI) ///
errorlag(WI) //用elect变量作为police的工具变量
estat impact
spivregress dui nondui vehicles i.dry (police = elect), dvarlag(WI) ///
errorlag(WI) ivarlag(WI: i.dry) //
**------------------------------------------
** spxtregress方法处理面板数据
//-------------------------------------------
use homicide_1960_1990, clear //用数据库来运行spxtregress
list _ID year in 1/8, sepby(_ID)
xtset _ID year //需要定义面板数据形式
spset //设置成空间回归形式的数据
bysort _ID: gen npanel = _N //检测一下是否_ID有重复
tab npanel //如果有重复，那需要使用spbalance
**空间随机效应-----------------------
xtreg hrate ln_population ln_pdensity gini i.year, re //普通的随机效应回归
spmatrix create contiguity WX if year == 1990 //生成紧邻权重矩阵WX(注意year==1990)
spmatrix export WX using contigx.xls //把空间权重矩阵W输出
spxtregress hrate ln_population ln_pdensity gini i.year, re dvarlag(WX) ///
errorlag(WX) //空间随机效应spatial lag of dv and error
estat impact gini //获得gini的平均效应
spmatrix create idistance MX if year == 1990 //生成距离倒数权重矩阵MX(注意year==1990)
spmatrix dir //显示一下所有的权重矩阵
spxtregress hrate ln_population ln_pdensity c.gini##i.year, re dvarlag(MX) ///
errorlag(MX) //空间随机效应,更换了一个权重矩阵且加入了交叉项
contrasts c.gini#year //检验gini与year交叉项的显著性
estat impact gini if year == 1960 //交叉项显著因此可以得到每一个year的gini效应
estat impact gini if year == 1970
estat impact gini if year == 1980
estat impact gini if year == 1990 //发现gini影响效应是在随着时间增强的
spxtregress hrate ln_population ln_pdensity c.gini##i.year, re sarpanel ///
dvarlag(MX) errorlag(MX) //当errorlag出现时,添加了sarpanel允许个体效应也像误差项那样的空间回归形式
**空间固定效应---------
xtreg hrate ln_population ln_pdensity gini, fe //普通的固定效应回归
spxtregress hrate ln_population ln_pdensity gini, fe dvarlag(WX) //空间固定效应spatial lag of dv
estat impact gini
spxtregress hrate ln_population ln_pdensity gini c.gini#i.year , ///
fe //如果不在每个panel里变动的话,就不能包括在fe模型里, 因为会被省略去
** 我们可以把内存里的权重矩阵删除掉---
local m "W M WX MX MI WI"
foreach x of local m {
spmatrix drop `x'
} //删掉之前的空间权重矩阵