温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者:游万海 (福州大学)
目录
图片来源:Golgher 和 Voss (2016)
由于面板数据模型所具有的众多优点 (刻画个体异质性,减弱模型共线性和增加自由度等),其被广泛应用于实证计量中。在 「 Stata: 面板数据模型-一文读懂」 文中,我们已对面板数据模型进行了介绍。
然而,当研究样本涉及到多个单元时 (如多个国家),研究单元间的空间相关性不可忽略。例如,在利用跨国样本研究环境污染 (如
正如地理学第一定律 ( Tobler's First Law ) 所说:任何事物都是与其他事物相关的,只不过相近的事物关联更紧密。空间计量计量学的发展为这一理论研究提供了有力的工具。 本推文将对面板数据框架下的空间计量模型 (空间面板数据模型) 进行讨论。本文主要分四部分内容进行介绍:
在介绍空间面板模型之前,我们先简单回顾一下普通的面板数据模型。这里以双向固定效应模型为例进行说明,详情参见「 Stata: 面板数据模型-一文读懂」。
面板数据模型同时包含了截面和时间两个维度,后文用
其中,
显然,
面板数据模型的主要用途之一就在于处理这些不可观测的个体效应或时间效应。当对所有的 reg y x
命令进行参数估计。
上述模型中解释变量仅仅纳入自身的因素,未考虑其他地区的一些因素。例如,一个国家的
空间面板数据模型进一步加入了空间滞后被解释变量
进一步,可将该模型写成向量形式:
其中,
与普通面板数据模型相比,空间面板数据模型通过引入空间权重矩阵来定义研究单元之间的关联方式和关联程度。从现有文献来看,常用的空间权重矩阵可以归纳为两类。一类为基于地理位置构造的,如 Rook 、 Queen、K-nearest 和 距离倒数次方权重矩阵;一类为基于社会经济因素的空间权重矩阵。前者的优点为直观且满足空间权重矩阵外生性假定,而后者的优点为有较强的经济含义且更符合实际应用背景。基于社会经济因素的空间权重矩阵,由于不同年份的社会经济指标值不同,其通常是随着时间变化而改变,但是此时模型通常较为复杂,实证中常取平均值构造空间权重矩阵 (此外,时变空间权重矩阵也被提出)。
为了对模型进行正确识别,必须保证空间权重矩阵的外生性。基于经济社会因素构造的空间权重矩阵有时经济含义较为明显,但是通常都不满足外生性假设。因此,在实际应用中最常用的空间权重矩阵是基于地理信息或者地理位置的,如 Rook,Queen 以及距离权重矩阵。
例如,在省域环境问题研究中,通常认为"邻近"省份的环境治理行为存在相似。这里的"邻近"可以认为从地理位置来看,两个省份存在地理接壤。若存在共同边界,元素设定为 1,否则设定为 0。假设我们有 5 个省份 (下文统称为 单元#),其地理分布如下图所示
则 Rook 空间权重矩阵是一个 5 行 5 列的方阵:
具体解释如下:
对于这种 孤零零的单元,另一种处理办法是放宽假设,可以认为只要存在顶点相接,就认为两地区为"邻居"关系,这就是 Queen 空间权重矩阵。按照这一设定,空间权重矩阵可定义为
前文已经提及,空间权重矩阵不仅可以定义 关联方式 (是否关联),而且可以定义 关联程度 (关联的紧密程度)。然而基于 Rook 和 Queen 的空间权重矩阵假定两个地区只要存在关联,其关联程度都是相同的 (等权重)。然而,这种假定往往有悖常理。例如,在环境污染研究中,一般认为距离近的地区,关联程度较高;又如,在国际贸易研究中,往来密切的地区之间关联程度较高。在这种情况下,基于空间距离的权重矩阵较为合适。基于距离的空间权重矩阵定义如下:
其中:
与基于地理信息构造的空间权重矩阵相比,经济权重矩阵的实际应用背景较强,其没有统一的数学表达式。例如,研究经济增长影响因素时,研究样本为世界各国,一般认为两国间的贸易往来越多,它们之间的经济增长水平联系越紧密,此时可根据两国的贸易往来数据构造空间权重矩阵,空间权重矩阵定义如下:
其中:
再如,基于 GDP 的空间权重可以定义为
其中,
当假设研究单元不相依时 (independence),在时刻
根据上述式子可知,模型系数和边际效应大小相同,解释变量
空间计量模型因为包含
也就是说,某一单元
为了解决空间计量模型系数难以解释问题,LeSage 和 Pace (2009) 提出直接效应 (direct effects, DE), 间接效应 (indirect effects, IE) 和 总效应 (total effects, TE)。对于空间 Durbin 模型,模型可以表示为
其中,在特定稳定条件下,
为了更为直观理解上述公式,以下用截面数据框架下 SAR 模型 (仅考虑
若研究样本为中国 34 个省级行政区域,经济发展水平为解释变量,记为
直接效应为上述矩阵对角元素的平均值,若
直接效应 表示的是某地区自变量对该因变量的影响大小,其包括反馈效应 (feedback effects, FE),即对其他地区的影响又会反过来影响该地区。从数值上来看,直接效应等于空间 Durbin 模型系数与反馈效应之和。具体地,反馈效应 指的是某个地区的解释变量 (如:技术创新) 会对其"邻近"地区的被解释变量 (如:可再生能源生产量),反过来又会影响该地区的被解释变量 (如:可再生能源生产量)。例如:某地区技术会被其"邻近"地区模仿与采用,用来生产可再生能源,这些地区又可将生产的再生能源出口到本地区,进而影响本地区的可再生能源量。
间接效应 为矩阵非对角元素的行平均,若
间接效应,又称为空间溢出效应 (spatial spillover effects),用于度量“邻近”地区的某个解释变量(如:技术创新)对本地区的被解释变量的影响 (如:可再生能源生产)。值得注意的是,不能仅仅依靠空间自回归系数
总效应为直接效应和间接效应之和,即
空间面板数据模型中包括了 splagvar
或者 spgenerate
命令生成变量 xtabond2
命令进行动态面板数据模型估计。在贝叶斯估计方法,更加详细可以参考 LeSage 和 Pace (2009)。本部分所介绍的 Stata 命令,主要是基于极大似然法进行参数估计。
Stata 15 版增加相应的空间计量模型估计 (help spxtregress
)和相关检验命令 (help spxtregress postestimation
),在估计模型后可以通过 estat impact
计算直接效应,间接效应和总效应。详细的手册可以参考「 Stata: Stata 15 PDF 电子手册」 。目前,spxtregress
还不支持动态空间面板数据模型估计,这里主要介绍 Federico Belotti 等人开发的 xsmle
命令,该命令同时支持空间动态面板数据模型的估计,该命令已发表在 Stata Journal,更为详细的可以参见
「 Spatial panel-data models using Stata 」。
以空间杜宾模型为例, xsmle
语法如下:
xsmle depvar [indepvars] [if] [in] [weight] , ///
wmat(name) model(sdm) [SDM_options]
为了方便说明,同时列出普通面板数据模型的命令语法,
xtreg depvar [indepvars] [if] [in] [weight] , fe [FE_options]
对比可以发现,xsmle
命令增加了 wmat() 和 model() 这两个选项。model() 用于定义空间面板数据模型形式,包括空间杜宾模型 (SDM),空间滞后模型 (Spatial Autoregressive model, SAR),空间误差模型 (Spatial Error model, SEM),空间自相关模型 (Spatial Autocorrelation model, SAC),以及广义空间面板随机效应模型(Generalized Spatial Panel Random Effects model, GSRE)。
前文已经介绍,利用空间计量模型的一个重要前提是空间权重矩阵的合理定义,wmat 选项就是用于空间权重矩阵,其可以为 Stata 中矩阵或者 spmat
目标(即利用 spmat
产生, help spmat
)。下文基于地理位置信息构造空间权重矩阵进行说明。
为了计算该类空间权重矩阵,首先需要获得研究单元的地理位置信息,如经度和维度。当研究样本为国家时,通常以国家的首都或者最大城市位置信息作为国家位置的代表。获取经度和维度有多种途径,比较方便的是下载到研究单元对象的 ESRI Shapefile (shapefile) 文件。中国的 shape 文件,包括省级,市级和县级的可以在国家基础信息中心申请下载 (各大论坛应该也有)。若研究的样本为世界各个国家,可以到 (http://www.diva-gis.org/gdata) 下载 。 例如,下载的 shape 文件包括:
首先,生成样本数据,这里假设研究样本为 10 个国家,时间为 1935 到 1954 年。
clear
input str16 countries
"China"
"United States"
"Germany"
"Japan"
"India"
"Brazil"
"Russia"
"South Africa"
"United Kingdom"
"Spain"
end
expand 20
sort countries
bys countries:gen year = _n + 1934
egen company = group(countries)
merge 1:1 company year using grunfeld
keep if _merge==3
drop _merge time
save G:\sample_sp,replace
若以 invest 为因变量,mvalue 和 kstock 为自变量。若采用普通面板固定效应模型,在声明面板数据结构后,可利用如下命令进行估计:
xtreg invest mvalue kstock,fe
若采用空间面板数据模型进行建模,则还需要找出这些国家的经度和维度等地理信息。前文已经说明了如何下载 shapefile 文件,该文件中即包含了各个国家的地理信息。那么可以将研究单元数据和 shapefile 进行合并。这里需要用到 shp2dta
命令读取 shapefile 文件,该命令为外部命令,安装方式参考「 Stata: 外部命令的搜索、安装与使用」 。
. cd "c:\"
. shp2dta using "TM_WORLD_BORDERS_SIMPL-0.2.shp", ///
database(data_db) coordinates(data_xy) ///
genid(id) gencentroids(c) replace
. use "data_db", clear
. rename NAME countries
. list countries LAT LON in 1/10
+-----------------------------------------+
| countries LAT LON |
|-----------------------------------------|
1. | Antigua and Barbuda 17.078 -61.783 |
2. | Algeria 28.163 2.632 |
3. | Azerbaijan 40.43 47.395 |
4. | Albania 41.143 20.068 |
5. | Armenia 40.534 44.563 |
|-----------------------------------------|
6. | Angola -12.296 17.544 |
7. | American Samoa -14.318 -170.73 |
8. | Argentina -35.377 -65.167 |
9. | Australia -24.973 136.189 |
10. | Bahrain 26.019 50.562 |
+-----------------------------------------+
读取 shapefile 后,可以发现有 LAT 和 LON 两个变量,分别代表了对应国家的纬度和经度。
接下来,利用该数据和样本数据 sample_sp.dta 根据 countries 这一变量进行合并,并将数据保存为 spatial_weight.dta。
merge 1:m countries using "G:\sample_sp.dta"
keep if _merge==3
drop _merge
keep countries year LAT LON invest mvalue kstock
sort countries year
by countries, sort: list, constant
bys countries: keep if _n==1
save G:\spatial_weight.dta,replace
为了构造空间权重矩阵,利用上述生成的 spatial_weight.dta 数据中的纬度和经度信息进行计算。这里可以利用 spwmatrix
生成空间权重矩阵。该命令可以构造多种形式的空间权重矩阵,包括基于地理位置信息的距离权重矩阵,基于经济社会变量的经济权重矩阵。此外,也可以导入 Geoda 生成的 gal 文件,或直接导入 txt 文件。
例如,可利用上面代码产生距离倒数空间权重矩阵
use "G:\spatial_weight.dta",clear
spwmatrix gecon LAT LON, wn(example_w2) ///
wtype(inv) r(3958.761) alpha(1) ///
xport(example_w2,txt) row replace
spmat import example_w2 using example_w2.txt,replace
spmat save example_w2 using example_w2m.spmat,replace
spmat use WW using example_w2m.spmat,replace
我们来看 spwmatrix
命令的语法:
spwmatrix gecon varlist [if] [in], ///
wname(wght_name) ///
[wtype(bin|inv|econ|invecon) ///
cart r(#) dband(numlist) ///
alpha(#) knn(#) ///
econvar(varname1) beta(#) Other_options]
wname() 填写生成的矩阵名称,wtype() 内填写生成矩阵类型,如 bin 表示 binary 二元权重矩阵,inv 表示生成逆距离矩阵,econ 为经济权重矩阵等。r() 内填写地球半径,单位可以是英里或者千米。alpha() 内填写数字,如 1 表示距离 1 次方倒数,2 表示距离平方的倒数。knn() 用于定义最邻近单元个数。row 表示对该矩阵进行行标准化。
在运行 spwmatrix
后,在本地会生成 example_w2.txt 文件,可以将其导入到 Stata 中查看,
insheet using "example_w2.txt", clear delimiter(" ") names
drop v1
rename (v#) (v#),renumber
export excel using "spatial_weight.xls",replace
list in 1/10
上述生成的 spatial_weight.xls文件就可以导入到 R 或 Matlab 中运行。
此外,也可以生成 k-nearest 空间权重矩阵,代码如下:
use G:\spatial_weight.dta,clear
spwmatrix gecon LAT LON, ///
wn(knear_w) knn(5) ///
xport(distant_nearest,txt) ///
row replace /*这里定义最近的5个地区为邻居*/
spmat import nearest_w using distant_nearest.txt,replace
spmat save nearest_w using nearest_w.spmat,replace
spmat use WW using nearest_w.spmat,replace
根据步骤 (2) ,我们已经生成了空间权重矩阵,接下来按照第 (1) 部分生成的 sample_sp.dta 数据进行实证建模
cd "H:\"
use spatial_weight.dta,clear
spwmatrix gecon LAT LON, ///
wn(example_w2) wtype(inv) ///
r(3958.761) alpha(1) ///
xport(example_w2,txt) row replace
spmat import example_w2 using "example_w2.txt", replace
spmat save example_w2 using "example_w2m.spmat", replace
spmat use WW using "example_w2m.spmat", replace
use sample_sp.dta,clear
rename company id
xtset id year,y
local xorg "invest mvalue kstock"
foreach f of local xorg {
gen ln`f' = log(`f')
}
est clear
eststo: xtreg lninvest lnmvalue lnkstock, fe /*普通面板数据模型*/
eststo:xsmle lninvest lnmvalue lnkstock, fe ///
model(sar) wmat(WW) type(both) ///
effects nolog /*空间面板数据模型:SAR model*/
估计结果显示如下:
convergence not achieved
Computing marginal effects standard errors using MC simulation...
SAR with spatial and time fixed-effects Number of obs = 200
Group variable: id Number of groups = 10
Time variable: year Panel length = 20
R-sq: within = 0.4258
between = 0.8555
overall = 0.8165
Mean of fixed-effects = 3.5886
Log-likelihood = 31.1740
------------------------------------------------------------------------------
lninvest | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Main |
lnmvalue | .3087483 .0771167 4.00 0.000 .1576024 .4598942
lnkstock | .0788984 .0306623 2.57 0.010 .0188015 .1389953
-------------+----------------------------------------------------------------
Spatial |
rho | -.4325611 .167113 -2.59 0.010 -.7600965 -.1050257
-------------+----------------------------------------------------------------
Variance |
sigma2_e | .0415918 .0042036 9.89 0.000 .0333528 .0498308
-------------+----------------------------------------------------------------
LR_Direct |
lnmvalue | .3189515 .0808705 3.94 0.000 .1604482 .4774549
lnkstock | .079349 .0306296 2.59 0.010 .019316 .1393819
-------------+----------------------------------------------------------------
LR_Indirect |
lnmvalue | -.0952446 .0387047 -2.46 0.014 -.1711044 -.0193847
lnkstock | -.0239403 .0128027 -1.87 0.061 -.0490331 .0011526
-------------+----------------------------------------------------------------
LR_Total |
lnmvalue | .223707 .064382 3.47 0.001 .0975205 .3498934
lnkstock | .0554087 .0221398 2.50 0.012 .0120155 .098802
------------------------------------------------------------------------------
从结果来看,出现了 convergence not achieved,表明在规定迭代次数内,模型未收敛。这时候可以通过 copy() 选项进行调整,大致的思路是先根据现有估计结果,将模型参数估计系数作为初始值进行迭代。代码如下:
eststo:xsmle lninvest lnmvalue lnkstock, ///
fe model(sar) wmat(WW) type(both) effects nolog ///
from(0.308 0.078 0.5 0.0418,copy)
结果显示如下:
Computing marginal effects standard errors using MC simulation...
SAR with spatial and time fixed-effects Number of obs = 200
Group variable: id Number of groups = 10
Time variable: year Panel length = 20
R-sq: within = 0.4213
between = 0.8549
overall = 0.8155
Mean of fixed-effects = 3.6617
Log-likelihood = 31.1653
------------------------------------------------------------------------------
lninvest | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Main |
lnmvalue | .3079884 .0772697 3.99 0.000 .1565426 .4594342
lnkstock | .0789791 .0307231 2.57 0.010 .0187628 .1391954
-------------+----------------------------------------------------------------
Spatial |
rho | -.4491497 .1674137 -2.68 0.007 -.7772745 -.1210249
-------------+----------------------------------------------------------------
Variance |
sigma2_e | .0417572 .0042414 9.85 0.000 .0334441 .0500703
-------------+----------------------------------------------------------------
LR_Direct |
lnmvalue | .3186822 .0811596 3.93 0.000 .1596123 .4777522
lnkstock | .079556 .0307443 2.59 0.010 .0192983 .1398136
-------------+----------------------------------------------------------------
LR_Indirect |
lnmvalue | -.0981599 .03883 -2.53 0.011 -.1742652 -.0220546
lnkstock | -.0247473 .0129782 -1.91 0.057 -.0501842 .0006895
-------------+----------------------------------------------------------------
LR_Total |
lnmvalue | .2205223 .0635574 3.47 0.001 .0959521 .3450926
lnkstock | .0548086 .0218984 2.50 0.012 .0118885 .0977288
------------------------------------------------------------------------------
(est20 stored)
可以发现,在调整初始值后,模型达到收敛 ( convergence )。在这里需要特别强调,有时候估计模型时也会出现 initial values not feasible ,这时候必须通过**copy()**选项进行设置调整。
现在我们来看下结果,Main 为模型估计系数
接下来,将 model(sar) 中修改为 sdm,即估计空间杜宾模型,结果如下:
eststo:xsmle lninvest lnmvalue lnkstock, fe ///
model(sdm) wmat(WW) type(both) effects nolog
Warning: All regressors will be spatially lagged
**convergence not achieved**
Computing marginal effects standard errors using MC simulation...
SDM with spatial and time fixed-effects Number of obs = 200
Group variable: id Number of groups = 10
Time variable: year Panel length = 20
R-sq: within = 0.6420
between = 0.8481
overall = 0.6868
Mean of fixed-effects = -3.0978
Log-likelihood = 34.2413
------------------------------------------------------------------------------
lninvest | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Main |
lnmvalue | .3992177 .0907614 4.40 0.000 .2213286 .5771068
lnkstock | .1580089 .0429793 3.68 0.000 .073771 .2422469
-------------+----------------------------------------------------------------
Wx |
lnmvalue | .60829 .3927921 1.55 0.121 -.1615683 1.378148
lnkstock | .5086279 .209802 2.42 0.015 .0974236 .9198323
-------------+----------------------------------------------------------------
Spatial |
rho | -.5926545 .1769285 -3.35 0.001 -.939428 -.245881
-------------+----------------------------------------------------------------
Variance |
sigma2_e | .0396852 .004072 9.75 0.000 .0317042 .0476662
-------------+----------------------------------------------------------------
LR_Direct |
lnmvalue | .3747412 .0854161 4.39 0.000 .2073287 .5421536
lnkstock | .1274166 .0350247 3.64 0.000 .0587695 .1960638
-------------+----------------------------------------------------------------
LR_Indirect |
lnmvalue | .2894559 .254896 1.14 0.256 -.210131 .7890428
lnkstock | .2878004 .1273755 2.26 0.024 .038149 .5374517
-------------+----------------------------------------------------------------
LR_Total |
lnmvalue | .6641971 .2879684 2.31 0.021 .0997893 1.228605
lnkstock | .415217 .1438093 2.89 0.004 .1333558 .6970781
------------------------------------------------------------------------------
(est31 stored)
模型不收敛,通过 copy() 调整初始值,代码如下
. eststo:xsmle lninvest lnmvalue lnkstock,fe ///
model(sdm) wmat(WW) type(both) effects nolog ///
from(0.3995296 0.1598841 0.6188835 0.5201294 0.689 0.0396793,copy)
. Warning: All regressors will be spatially lagged
Computing marginal effects standard errors using MC simulation...
SDM with spatial and time fixed-effects Number of obs = 200
Group variable: id Number of groups = 10
Time variable: year Panel length = 20
R-sq: within = 0.6420
between = 0.8481
overall = 0.6868
Mean of fixed-effects = -3.0978
Log-likelihood = 34.2412
------------------------------------------------------------------------------
lninvest | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Main |
lnmvalue | .3992206 .09085 4.39 0.000 .2211578 .5772834
lnkstock | .1580265 .0430199 3.67 0.000 .0737089 .242344
-------------+----------------------------------------------------------------
Wx |
lnmvalue | .6083891 .3931709 1.55 0.122 -.1622116 1.37899
lnkstock | .5087355 .2099962 2.42 0.015 .0971505 .9203204
-------------+----------------------------------------------------------------
Spatial |
rho | -.5929535 .1769972 -3.35 0.001 -.9398616 -.2460454
-------------+----------------------------------------------------------------
Variance |
sigma2_e | .0397628 .0040881 9.73 0.000 .0317502 .0477754
-------------+----------------------------------------------------------------
LR_Direct |
lnmvalue | .3747347 .0854988 4.38 0.000 .2071601 .5423092
lnkstock | .1274166 .0350578 3.63 0.000 .0587046 .1961286
-------------+----------------------------------------------------------------
LR_Indirect |
lnmvalue | .2894315 .2550993 1.13 0.257 -.2105539 .7894169
lnkstock | .2877979 .1274783 2.26 0.024 .0379451 .5376507
-------------+----------------------------------------------------------------
LR_Total |
lnmvalue | .6641662 .2881908 2.30 0.021 .0993227 1.22901
lnkstock | .4152145 .1439223 2.88 0.004 .133132 .697297
------------------------------------------------------------------------------
(est32 stored)
至此,本文讲解了空间权重矩阵构造和静态空间面板数据模型估计。此外,还有模型的选择(SDM vs SAR vs SEM),大家可以参考「 Spatial panel-data models using Stata 」。
此外,动态空间面板模型估计可以通过调整dlag选项进行定义:
dlag(1):
dlag(2):
dlag(3):