一文读懂空间计量及stata应用（附lr检验、空间可视化、权重矩阵、检验、模型dofile等）

本文来源：Spatial Econometrics Methods using Stata，更多信息详见https://www.researchgate.net/publication/313736408

转载请注明来源

Tutorial to replicate results of Workshop:

Spatial Econometrics Methods using Stata"

LISER (Luxembourg Institute of Socio-Economic Research)

Luxembourg, 15th February 2017.

Author: Marcos Herrera (CONICET-IELDE, UNSa, Argentina)

e-mail: [e-mail]

主要包括如下内容：

1、空间数据的操作与可视化: 使用地图表示。

2、空间权重矩阵的建立和空间自相关测试:Moran's I test, Geary's c test and Getis-Ord's G test. - LISA: Local Moran's I test.

3、基本空间计量经济学(截面数据):SLM、SEM、SARAR、SDM。

MLE和GMM下的空间参数估计。 

测试以确定适当的计量经济设定。

异方差性和内生性校正。

4、高级空间计量经济学: 

面板数据中空间效应的介绍，

静态和动态空间面板数据的估计。

1

空间计量经济学命令

需要下载安装外部命令清单如下：

version 14.0*ssc install spmap*ssc install shp2dta*net install sg162, from(http://www.stata.com/stb/stb60)*net install st0292, from(http://www.stata-journal.com/software/sj13-2)*net install spwmatrix, from(http://fmwww.bc.edu/RePEc/bocode/s)*net install splagvar, from(http://fmwww.bc.edu/RePEc/bocode/s)*ssc install xsmle.pkg*ssc install xtcsd

2

空间数据处理

1、导入数据Read the information shape in Stata

shp2dta using nuts2_164, database(data_shp) coordinates(coord) ///genid(id) genc(c) replace

结果为：

2、查看数据

use data_shp, cleardescribe

3、Themeless map (without information)

spmap using coord, id(id) note("Europe, EU15")

4、If there are some polygons without information: for example, Sweden and Finland

drop if id==3|id==5|id==6|id==164|id==7|id==8|id==12|id==4|id==2|id==1|id==11|id==12spmap using coord, id(id) note("Europe without Finland and Sweden, EU15")

5、 Import from Excel and save in Stata format 

clear allimport excel migr_unemp07_12.xls, firstrowsave migr_unemp.dta, replaceuse migr_unempdescribe

结果为：

6、 合并数据

use data_shp, clearmerge 1:1 POLY_ID using migr_unemp, gen(union) forceassert union==3

7、 合并数据

drop unionsave migr_unemp_shp.dta, replace

8、Showing the information using maps

use migr_unemp_shp.dta

9、Quantile map:

format U2012 %12.1fspmap U2012 using coord, id(id) clmethod(q) title("Unemployment rate") ///legend(size(medium) position(5)) fcolor(Blues2) note("Europe, 2012" "Source: Eurostat")

format NM2012 %12.1fspmap NM2012 using coord, id(id) clmethod(q) title("Net migration rate") ///legend(size(medium) position(5)) fcolor(BuRd) note("Europe, 2012" "Source: Eurostat")

10、Equal interval maps

spmap U2012 using coord, id(id) clmethod(e) title("Unemployment rate") ///legend(size(medium) position(5)) fcolor(Blues2) note("Europe, 2012" "Source: Eurostat")

spmap NM2012 using coord, id(id) clmethod(e) title("Net migration rate") ///legend(size(medium) position(5)) fcolor(BuRd) note("Europe, 2012" "Source: Eurostat")

11、Box maps

spmap U2012 using coord, id(id) clmethod(boxplot) title("Unemployment rate") ///legend(size(medium) position(5)) fcolor(Heat) note("Europe, 2012" "Source: Eurostat")
spmap NM2012 using coord, id(id) clmethod(boxplot) title("Net migration rate") ///legend(size(medium) position(5)) fcolor(Rainbow) note("Europe, 2012" "Source: Eurostat")
graph hbox U2012, asyvars ytitle("")graph hbox NM2012, asyvars ytitle("")

12、 Deviation maps

spmap U2012 using coord, id(id) clmethod(s) title("Unemployment rate") ///legend(size(medium) position(5)) fcolor(Blues2) note("Europe, 2012" "Source: Eurostat")
spmap NM2012 using coord, id(id) clmethod(s) title("Net migration rate") ///legend(size(medium) position(5)) fcolor(BuRd) note("Europe, 2012" "Source: Eurostat")


13、  Combination of points and polygons using both variables: 

spmap U2012 using coord, id(id) fcolor(RdYlBu) cln(8) point(data(migr_unemp_shp) xcoord(x_c) ///ycoord(y_c) deviation(NM2012) sh(T) fcolor(dknavy) size(*0.3)) legend(size(medium) position(5)) legt(Unemployment) ///note("Solid triangles indicate values over the mean of net-migration." "Europa, 2012. Source: Eurostat")


spmap NM2012 using coord, id(id) fcolor(RdYlBu) cln(8) diagram(var(U2012) xcoord(x_c) ycoord(y_c) ///fcolor(gs2) size(1)) legend(size(medium) position(5)) legstyle(3) legt(Net migration) ///note(" " "Boxes indicate values of unemployment." "Europe, 2012. Source: Eurostat")

14、  空间权重矩阵以及空间相关检验

spmat contiguity Wcontig using "coord.dta", id(id)* Problem with conguity criterion: 5 islands.
* We choose k-nn: 5 nearest neighbours row-standardizedspwmatrix gecon y_c x_c, wn(W5st) knn(5) row con
* We need the spatial W as a SPMAT object:* First, we generate W 5nn binary and then we export as txt
spwmatrix gecon y_c x_c, wn(W5bin) knn(5) xport(W5bin,txt) replace* Read the txt file and to adapt format for SPMATinsheet using "W5bin.txt", delim(" ") cleardrop in 1rename v1 idsave "W5bin.dta", replace* Generate SPMAT object: W5 row-standardize spmat dta W5_st v*, id(id) norm(row)spmat summarize W5_st, linksspmat graph W5_st

15、TIPS FOR MATRICES

1. From .GAL file to Stata format  //gal格式导入stata中* spwmatrix import using matrix.gal, wname(W_geoda)

* 2. From .GWT file to Stata format //gwt格式导入stata中* spmat import W_knn using knn.gwt, geoda
* 3. From SPMAT object to SPATWMAT object * spmat export Wknn using "Wknn_noid.txt", noid replace* insheet using "Wknn_noid.txt", delim(" ") clear* drop in 1* save "Wknn_noid.dta", replace*    spatwmat using "Wcont_noid.dta", name(Wks) standardize

16、空间相关检验

Moran I test, Geary's c test and Getis-Ord G test.
use migr_unemp_shp.dta, clear
spatgsa U2012, w(W5st) moran geary twospatgsa NM2012, w(W5st) moran geary two

17、空间相关检验

* For Getis-Ord test we need a binary matrix:* spatwmat using "W5bin_noid.dta", name(W5b)* This binary matrix has been created previously by spwmatrixspatgsa U2012, w(W5bin) go twospatgsa NM2012, w(W5bin) go two

18、空间相关检验

* Moran I scatterplotsplagvar U2012, wname(W5st) wfrom(Stata) ind(U2012) order(1) plot(U2012) moran(U2012)splagvar NM2012, wname(W5st) wfrom(Stata) ind(NM2012) order(1) plot(NM2012) moran(NM2012)
* Local Moran I (LISA)genmsp_v0 U2012, w(W5st)graph twoway (scatter Wstd_U2012 std_U2012 if pval_U2012>=0.05, msymbol(i) mlabel ///(id) mlabsize(*0.6) mlabpos(c)) (scatter Wstd_U2012 std_U2012 if pval_U2012<0.05, ///msymbol(i) mlabel (id) mlabsize(*0.6) mlabpos(c) mlabcol(red)) (lfit Wstd_U2012 ///std_U2012), yline(0, lpattern(--)) xline(0, lpattern(--)) xlabel(-1.5(1)4.5, ///labsize(*0.8)) xtitle("{it:z}") ylabel(-1.5(1)3.5, angle(0) labsize(*0.8)) ///ytitle("{it:Wz}") legend(off) scheme(s1color) title("Local Moran I of Unemployment rate")
spmap msp_U2012 using coord, id(id) clmethod(unique) title("Unemployment rate") ///legend(size(medium) position(4)) ndl("No signif.") fcolor(blue red) ///note("Europe, 2012" "Source: Eurostat")


19、空间计量分析

OLS estimation************************************************************************************use migr_unemp_shp, clearreg U2012 NM2012

20、空间计量分析

* Spatial tests
spwmatrix gecon y_c x_c, wn(W5st) knn(5) rowspatdiag, weights(W5st)

21、空间计量分析

************************************************************************************* Spatial models using Maximum Likelihood (ML)************************************************************************************* Spatial Lag Model (SLM) with W5_st spmat objectspreg ml U2012 NM2012, id(id) dlmat(W5_st)estimates store SLM_ml
* Spatial Error Model (SEM)spreg ml U2012 NM2012, id(id) elmat(W5_st)estimates store SEM_ml
* Spatial autoregressive SARAR model: combine SLM-SEMspreg ml U2012 NM2012, id(id) dlmat(W5_st) elmat(W5_st)estimates store SARAR_ml
* Spatial Durbin model (SDM)spmat lag wx_NM2012 W5_st NM2012
spreg ml U2012 NM2012 wx_NM2012, id(id) dlmat(W5_st)estimates store SDM_ml

结果为：

22、LR检验，SEM模型与SDM模型

* Selecting between SDM and SEM: LR_comfaclrtest SDM_ml SEM_ml
estimates table SLM_ml SEM_ml SARAR_ml SDM_ml CLIFFORD_ml, b(%7.2f) star(0.1 0.05 0.01)* Others alternative commands: "spmlreg" de Jeanty o "spatreg" de Pisati

23、空间计量分析


************************************************************************************* Spatial model using Instrumental Variables / Generalized method of moments(IV-GMM)************************************************************************************* Spatial Lag Model (SLM)spivreg U2012 NM2012, dl(W5_st) id(id)
* SLM could be estimated using habitual commands in Stataspmat lag wx_U2012 W5_st U2012spmat lag wx2_NM2012 W5_st wx_NM2012ivregress 2sls U2012 NM2012 (wx_U2012 = wx_NM2012 wx2_NM2012)
* Spatial Error Model (SEM)spivreg U2012 NM2012, el(W5_st) id(id)
* SARAR Modelspivreg U2012 NM2012, dl(W5_st) el(W5_st) id(id)
* Spatial Durbin Model (SDM)spivreg U2012 NM2012 wx_NM2012, dl(W5_st) id(id)ereturn list* Same result of SDM using ivregress:spmat lag wx3_NM2012 W5_st wx2_NM2012ivregress 2sls U2012 NM2012 wx_NM2012 (wx_U2012 = wx2_NM2012 wx3_NM2012)
* Cliff-Ord Modelspivreg U2012 NM2012 wx_NM2012, dl(W5_st) el(W5_st) id(id)
************************************************************************************* Corrections for heteroskedasticity and endogeneity
* SLM assuming endogeneity in net migration and heteroskedasticityspivreg U2012 NM2012 (NM2012 = NM2009), dl(W5_st) id(id) het
* SEM assuming endogeneity in net migration and heteroskedasticityspivreg U2012 NM2012 (NM2012 = NM2009), el(W5_st) id(id) het
* SARAR with heteroskedasticity correctionspivreg U2012 NM2012, el(W5_st) dl(W5_st) id(id) het
* Alternative command: "spreg gs2sls"
************************************************************************************* Interpretation of spatial estimation************************************************************************************
* Remembering the SLM estimated under MLuse migr_unemp_shp, clearspreg ml U2012 NM2012, dl(W5_st) id(id)
* Read W and betas in MATA language spmat getmatrix W5_st Wmata:b = st_matrix("e(b)")blambda = b[1,3]lambdaS = luinv(I(rows(W))-lambda*W)end
* Using formulas of effects
* Total effectsmata: (b[1,1]/rows(W))*sum(S)
* Direct effectsmata: (b[1,1]/rows(W))*trace(S)
* Indirect effects (spatial spillovers)mata: (b[1,1]/rows(W))*sum(S) - (b[1,1]/rows(W))*trace(S)



24、高级空间计量分析

* Reshaping format: from wide to long formatreshape long NM U, i(id) j(year)tsset id year
xtset id yearxtdes
************************************************************************************* STATIC MODELS************************************************************************************
************************************************************************************* Spatial deteccionxtreg U NM, fextcsd, pes abs
xtreg U NM, rextcsd, pes abs
* Fixed Effects
* SLM modelxsmle U NM, fe wmat(W5_st) mod(sar) hausmanxsmle U NM, fe type(ind, leeyu) wmat(W5_st) mod(sar) effectsestimates store SLM_fe
* SEM modelxsmle U NM, fe emat(W5_st) mod(sem) hausmanestimates store SEM_fe
* SARAR modelxsmle U NM, fe wmat(W5_st) emat(W5_st) mod(sac) effectsestimates store SARAR_fe
* SDM modelxsmle U NM, fe type(ind) wmat(W5_st) mod(sdm) hausman effectsestimates store SDM_fe * Common factor Testtestnl ([Wx]NM = -[Spatial]rho*[Main]NM)
estimates table SLM_fe SEM_fe SARAR_fe SDM_fe, b(%7.2f) star(0.1 0.05 0.01) statistics(aic)
* Random Effects
* SLM modelxsmle U NM, wmat(W5_st) mod(sar)* SEM modelxsmle U NM, emat(W5_st) mod(sem)* SDM modelxsmle U NM, wmat(W5_st) mod(sdm)
************************************************************************************* DYNAMIC MODELS************************************************************************************
************************************************************************************* Detecting serial dependencextserial U NM
* dynSLM model (named as SAR)xsmle U NM, dlag(1) fe wmat(W5_st) type(ind) mod(sar) effects nsim(499)estimates store dynSLM_1 xsmle U NM, dlag(2) fe wmat(W5_st) type(ind) mod(sar) effects nsim(499)estimates store dynSLM_2xsmle U NM, dlag(3) fe wmat(W5_st) type(ind) mod(sar) effects nsim(499)estimates store dynSLM_3estimates table dynSLM_1 dynSLM_2 dynSLM_3, b(%7.2f) star(0.1 0.05 0.01) statistics(aic)
* dynSDM model*xsmle U NM, dlag(1) fe wmat(W5_st) type(ind) mod(sdm) effects nsim(499)*estimates store dynSDM_1*xsmle U NM, dlag(2) fe wmat(W5_st) type(ind) mod(sdm) effects*estimates store dynSDM_2*xsmle U NM, dlag(3) fe wmat(W5_st) type(ind) mod(sdm) effects*estimates store dynSDM_3