【空间计量】（非原创）

https://www.sohu.com/a/386218186_698752

目前，空间计量经济学研究包括以下四个感兴趣的领域：

计量经济模型中空间效应的确定； 合并了空间影响的模型的估计；空间效应存在的说明、检验和诊断；空间预测。

空间计量经济学模型有多种类型（Anselin，et al. 2004）。 首先介绍纳入了空间效应（空间相关和空间差异）、适用于截面数据的空间常系数回归模型，包括空间滞后模型（Spatial Lag Model，SLM）与空间误差模型（Spatial Error Model，SEM）两种，以及空间变系数回归模型——地理加权回归模型（Geographical Weighted Regression，GWR）。适用于时间序列和截面数据合成的空间面板数据计量经济学模型将在以后予以介绍。

空间滞后模型（Spatial Lag Model，SLM）主要是探讨各变量在一地区是否有扩散现象（溢出效应）。其模型表达式为：参数 反映了自变量对因变量的影响，空间滞后因变量 是一内生变量，反映了空间距离对区域行为的作用。区域行为受到文化环境及与空间距离有关的迁移成本的影响，具有很强的地域性（Anselin et al.，1996）。由于SLM模型与时间序列中自回归模型相类似，因此SLM也被称作空间自回归模型（Spatial Autoregressive Model，SAR）。

空间误差模型（Spatial Error Model，SEM）存在于扰动误差项之中的空间依赖作用，度量了邻近地区关于因变量的误差冲击对本地区观察值的影响程度。由于SEM模型与时间序列中的序列相关问题类似，也被称为空间自相关模型（Spatial Autocorrelation Model，SAC）。

估计技术：鉴于空间回归模型由于自变量的内生性，对于上述两种模型的估计如果仍采用OLS，系数估计值会有偏或者无效，需要通过IV、ML或GLS、GMM等其他方法来进行估计。Anselin（1988）建议采用极大似然法估计空间滞后模型（SLM）和空间误差模型（SEM）的参数。

空间自相关检验与SLM、SEM的选择：判断地区间创新产出行为的空间相关性是否存在，以及SLM和SEM那个模型更恰当，一般可通过包括Moran’s I检验、两个拉格朗日乘数（Lagrange Multiplier）形式LMERR、LMLAG及其稳健（Robust）的R-LMERR、R-LMLAG）等形式来实现。由于事先无法根据先验经验推断在SLM和SEM模型中是否存在空间依赖性，有必要构建一种判别准则，以决定哪种空间模型更加符合客观实际。Anselin和Florax（1995）提出了如下判别准则：如果在空间依赖性的检验中发现LMLAG较之LMERR在统计上更加显著，且R-LMLAG显著而R-LMERR不显著，则可以断定适合的模型是空间滞后模型；相反，如果LMERR比LMLAG在统计上更加显著，且R-LMERR显著而R-LMLAG不显著，则可以断定空间误差模型是恰当的模型。

除了拟合优度R2检验以外，常用的检验准则还有：自然对数似然函数值（Log likelihood，LogL）、似然比率（Likelihood Ratio，LR）、赤池信息准则（Akaike information criterion，AIC）、施瓦茨准则（Schwartz criterion，SC）。对数似然值越大，AIC和SC值越小，模型拟合效果越好。这几个指标也用来比较OLS估计的经典线性回归模型和SLM、SEM，似然值的自然对数最大的模型最好。

空间变系数回归模型及估计：就目前国内外的研究来看，大多直接假定横截面单元是同质的，即地区或企业之间没有差异。传统的OLS只是对参数进行“平均”或“全域”估计，不能反映参数在不同空间的空间非稳定性（吴玉鸣，李建霞，2006；苏方林，2007）。 当用横截面数据建立计量经济学模型时，由于这种数据在空间上表现出的复杂性、自相关性和变异性，使得解释变量对被解释变量的影响在不同区域之间可能是不同的，假定区域之间的经济行为在空间上具有异质性的差异可能更加符合现实。空间变系数回归模型（Spatial Varying-Coefficient Regression Model）中的地理加权回归模型（Geographical Weighted Regression，GWR）是一种解决这种问题的有效方法。 、空间计量主要命令

spmat 生成空间权重矩阵

spatwmat 用于定义空间权重矩阵

spatgsa 用于全局空间自相关检验

gsa表示global spatial autocorrelation

spatlsa 进行局部空间自相关检验

lsa表示local spatial autocorrelation

spatcorr 考察空间自相关指标对距离临界值d的依赖性

spatdiag 针对ols回归结果，考察是否存在空间效应

spatreg 估计空间滞后与空间误差模型

空间面板主要命令为：help xsmle

Spatial Autoregressive (SAR) model

xsmle depvar [indepvars] [if] [in] [weight] , wmat(name) model(sar) [SAR_options]

Spatial Durbin (SDM) model

xsmle depvar [indepvars] [if] [in] [weight] , wmat(name) model(sdm) [SDM_options]

Spatial Autocorrelation (SAC) model

xsmle depvar [indepvars] [if] [in] [weight] , wmat(name) emat(name) model(sac) [SAC_options]

Spatial Error (SEM) model

xsmle depvar [indepvars] [if] [in] [weight] , emat(name) model(sem) [SEM_options]

Generalized Spatial Panel Random Effects (GSPRE) model

xsmle depvar [indepvars] [if] [in] [weight] , wmat(name) model(gspre) [emat(name) GSPRE_options]

潜变量模型的类型

近年来，潜变量分析模型领域最重要的革新是将连续型潜变量和类别型潜变量整合到一个模型中即混合模型(Hybrids Model, Muthén, 2008)。表1根据连续-类别潜变量和横断面-纵向研究设计将潜变量模型划分为不同的类别。

表1潜变量类型—扩展

连续潜变量

类别潜变量

混合

横断面模型

Cross-section Models

因子分析模型，SEM

潜类别分析(Latent Class Analysis, LCA)，

回归混合模型(Regression Mixture Modeling, RMM)

因子混合模型(Factor Mixture Modeling, FMM)

纵向模型

Longitudinal Models

潜变量增长模型

(Latent Growth Curve Model, LGCM)

潜在转换分析(Latent Transition Analysis, LTA)；

潜类别增长模型(latent class growth modeling, LCGM)

增长混合模型(Growth Mixture Model, GMM)

stata命令大全

********* 面板数据计量分析与软件实现 *********

说明：以下do文件相当一部分内容来自于中山大学连玉君STATA教程，感谢他的贡献。本人做了一定的修改与筛选。

*----------面板数据模型

* 1.静态面板模型：FE 和RE

* 2.模型选择：FE vs POLS, RE vs POLS, FE vs RE （pols混合最小二乘估计）

* 3.异方差、序列相关和截面相关检验

* 4.动态面板模型（DID-GMM,SYS-GMM）

* 5.面板随机前沿模型

* 6.面板协整分析（FMOLS,DOLS）

*** 说明：1-5均用STATA软件实现， 6用GAUSS软件实现。

* 生产效率分析（尤其指TFP）：数据包络分析（DEA）与随机前沿分析（SFA）

***

说明：DEA由DEAP2.1软件实现，SFA由Frontier4.1实现，尤其后者，侧重于比较C-D与Translog生产函数，一步法与两步法的区别。常应用于地区经济差异、FDI溢出效应（Spillovers

Effect）、工业行业效率状况等。

* 空间计量分析：SLM模型与SEM模型

*说明：STATA与Matlab结合使用。常应用于空间溢出效应（R&D）、财政分权、地方政府公共行为等。

* ---------------------------------

* -------- 一、常用的数据处理与作图 -----------

* ---------------------------------

* 指定面板格式

xtset id year （id为截面名称，year为时间名称）

xtdes /*数据特征*/

xtsum logy h /*数据统计特征*/

sum logy h /*数据统计特征*/

*添加标签或更改变量名

label var h "人力资本"

rename h hum

*排序

sort id year /*是以STATA面板数据格式出现*/

sort year id /*是以DEA格式出现*/

*删除个别年份或省份

drop if year<1992

drop if id==2 /*注意用==*/

*如何得到连续year或id编号（当完成上述操作时，year或id就不连续，为形成panel格式，需要用egen命令）

egen year_new=group(year)

xtset id year_new

**保留变量或保留观测值

keep inv /*删除变量*/

**或

keep if year==2000

**排序

sort id year /*是以STATA面板数据格式出现

sort year id /*是以DEA格式出现

**长数据和宽数据的转换

*长>>>宽数据

reshape wide logy,i(id) j(year)

*宽>>>长数据

reshape logy,i(id) j(year)

**追加数据（用于面板数据和时间序列）

xtset id year

*或者

xtdes

tsappend,add(5) /表示在每个省份再追加5年，用于面板数据/

tsset

*或者

tsdes

.tsappend,add(8) /表示追加8年，用于时间序列/

*方差分解，比如三个变量Y,X,Z都是面板格式的数据，且满足Y=X+Z，求方差var(Y),协方差Cov(X,Y)和Cov（Z,Y）

bysort year:corr Y X Z,cov

**生产虚拟变量

*生成年份虚拟变量

tab year,gen(yr)

*生成省份虚拟变量

tab id,gen(dum)

**生成滞后项和差分项

xtset id year

gen ylag=l.y /*产生一阶滞后项)，同样可产生二阶滞后项*/

gen ylag2=L2.y

gen dy=D.y /*产生差分项*/

*求出各省2000年以前的open inv的平均增长率

collapse (mean) open inv if year<2000,by(id)

变量排序，当变量太多，按规律排列。可用命令

aorder

或者

order fdi open insti

*-----------------

* 二、静态面板模型

*-----------------

*--------- 简介 -----------

* 面板数据的结构(兼具截面资料和时间序列资料的特征)

use product.dta, clear

browse

xtset id year

xtdes

* ---------------------------------

* -------- 固定效应模型 -----------

* ---------------------------------

* 实质上就是在传统的线性回归模型中加入 N-1 个虚拟变量，

* 使得每个截面都有自己的截距项，

* 截距项的不同反映了个体的某些不随时间改变的特征

*

* 例如： lny = a_i + b1*lnK + b2*lnL + e_it

* 考虑中国29个省份的C-D生产函数

*******-------画图------*

*散点图+线性拟合直线

twoway (scatter logy h) (lfit logy h)

*散点图+二次拟合曲线

twoway (scatter logy h) (qfit logy h)

*散点图+线性拟合直线+置信区间

twoway (scatter logy h) (lfit logy h) (lfitci logy h)

*按不同个体画出散点图和拟合线，可以以做出fe vs re的初判断*

twoway (scatter logy h if id<4) (lfit logy h if id<4) (lfit logy h if

id==1) (lfit logy h if id==2) (lfit logy h if id==3)

*按不同个体画散点图,so beautiful!!!*

graph twoway scatter logy h if id==1 || scatter logy h if id==2,msymbol(Sh)

|| scatter logy h if id==3,msymbol(T) || scatter logy h if id==4,msymbol(d) || ,

legend(position(11) ring(0) label(1 "北京") label(2 "天津") label(3 "河北") label(4

"山西"))

**每个省份logy与h的散点图，并将各个图形合并

twoway scatter logy h,by(id) ylabel(,format(%3.0f))

xlabel(,format(%3.0f))

*每个个体的时间趋势图*

xtline h if id<11,overlay legend(on)

* 一个例子：中国29个省份的C-D生产函数的估计

tab id, gen(dum)

list

* 回归分析

reg logy logk logl dum*,

est store m_ols

xtreg logy logk logl, fe

est store m_fe

est table m_ols m_fe, b(%6.3f) star(0.1 0.05 0.01)

* Wald 检验

test logk=logl=0

test logk=logl

* stata的估计方法解析

* 目的：如果截面的个数非常多，那么采用虚拟变量的方式运算量过大

* 因此，要寻求合理的方式去除掉个体效应

* 因为，我们关注的是 x 的系数，而非每个截面的截距项

* 处理方法：

*

* y_it = u_i + x_it*b + e_it (1)

* ym_i = u_i + xm_i*b + em_i (2) 组内平均

* ym = um + xm*b + em (3) 样本平均

* (1) - (2), 可得：

* (y_it - ym_i) = (x_it - xm_i)*b + (e_it - em_i) （4） /*within estimator*/ *

(4)+(3), 可得：

* (y_it-ym_i+ym) = um + (x_it-xm_i+xm)*b + (e_it-em_i+em)

* 可重新表示为：

* Y_it = a_0 + X_it*b + E_it

* 对该模型执行 OLS 估计，即可得到 b 的无偏估计量

**stata后台操作，揭开fe估计的神秘面纱！！！

egen y_meanw = mean(logy), by(id) /*个体内部平均*/

egen y_mean = mean(logy) /*样本平均*/

egen k_meanw = mean(logk), by(id)

egen k_mean = mean(logk)

egen l_meanw = mean(logl), by(id)

egen l_mean = mean(logl)

gen dyw = logy - y_meanw

gen dkw = logk - k_meanw

gen dlw=logl-l_meanw

reg dyw dkw dlw,nocons

est store m_stata

gen dy = logy - y_meanw + y_mean

gen dk = logk - k_meanw +k_mean

gen dl=logl-l_meanw+l_mean

reg dy dk dl

est store m_stata

est table m_*, b(%6.3f) star(0.1 0.05 0.01)

* 解读 xtreg,fe 的估计结果

xtreg logy h inv gov open,fe

*-- R^2

* y_it = a_0 + x_it*b_o + e_it (1) pooled OLS

* y_it = u_i + x_it*b_w + e_it (2) within estimator

* ym_i = a_0 + xm_i*b_b + em_i (3) between estimator

*

* -->R-sq: within 模型(2)对应的R2，是一个真正意义上的R2

* -->R-sq: between corr{xm_i*b_w,ym_i}^2

* -->R-sq: overall corr{x_it*b_w,y_it}^2

*

*-- F(4,373) = 855.93检验除常数项外其他解释变量的联合显著性

*

*

*-- corr(u_i, Xb) = -0.2347

*

*-- sigma_u, sigma_e, rho

* rho = sigma_u^2 / (sigma_u^2 + sigma_e^2)

dis e(sigma_u)^2 / (e(sigma_u)^2 + e(sigma_e)^2)

*

* 个体效应是否显著？

* F(28, 373) = 338.86 H0: a1 = a2 = a3 = a4 = a29

* Prob >F = 0.0000 表明，固定效应高度显著

*---如何得到调整后的 R2,即 adj-R2 ？

ereturn list

reg logy h inv gov open dum*

*---拟合值和残差

* y_it = u_i + x_it*b + e_it

* predict newvar, [option]

/*

xb xb, fitted valuesthe default

stdp calculate standard error of the fitted values

ue u_i + e_it, the combined residual

xbu xb + u_i, prediction including effect

u u_i, the fixed- or random-error component

e e_it, the overall error component */

xtreg logy logk logl, fe

predict y_hat

predict a , u

predict res,e

predict cres, ue

gen ares = a + res

list ares cres in 1/10

* ---------------------------------

* ---------- 随机效应模型 ---------

* ---------------------------------

* y_it = x_it*b + (a_i + u_it)

* = x_it*b + v_it

* 基本思想：将随机干扰项分成两种

* 一种是不随时间改变的，即个体效应 a_i

* 另一种是随时间改变的，即通常意义上的干扰项 u_it

* 估计方法：FGLS

* Var(v_it) = sigma_a^2 + sigma_u^2

* Cov(v_it,v_is) = sigma_a^2

* Cov(v_it,v_js) = 0

* 利用Pooled OLS，Within Estimator, Between Estimator

* 可以估计出sigma_a^2和sigma_u^2,进而采用GLS或FGLS

* Re估计量是Fe估计量和Be估计量的加权平均

* yr_it = y_it - theta*ym_i

* xr_it = x_it - theta*xm_i

* theta = 1 - sigma_u / sqrt[(T*sigma_a^2 + sigma_u^2)]

* 解读 xtreg,re 的估计结果

use product.dta, clear

xtreg logy logk logl, re

*-- R2

* -->R-sq: within corr{(x_it-xm_i)*b_r, y_it-ym_i}^2

* -->R-sq: between corr{xm_i*b_r,ym_i}^2

* -->R-sq: overall corr{x_it*b_r,y_it}^2

* 上述R2都不是真正意义上的R2，因为Re模型采用的是GLS估计。

*

* rho = sigma_u^2 / (sigma_u^2 + sigma_e^2)

dis e(sigma_u)^2 / (e(sigma_u)^2 + e(sigma_e)^2)

*

* corr(u_i, X) = 0 (assumed)

* 这是随机效应模型的一个最重要，也限制该模型应用的一个重要假设

* 然而，采用固定效应模型，我们可以粗略估计出corr(u_i, X)

xtreg market invest stock, fe

*

* Wald chi2(2) = 10962.50 Prob>chi2 = 0.0000

*-------- 时间效应、模型的筛选和常见问题

*---------目录--------

* 7.2.1 时间效应（双向固定(随机)效应模型）

* 7.2.2 模型的筛选

* 7.2.3 面板数据常见问题

* 7.2.4 面板数据的转换

* ----------------------------------

* ------------时间效应--------------

* ----------------------------------

* 单向固定效应模型

* y_it = u_i + x_it*b + e_it

* 双向固定效应模型

* y_it = u_i + f_t + x_it*b + e_it

qui tab year, gen(yr)

drop yr1

xtreg logy logk logl yr*, fe

* 随机效应模型中的时间效应

xtreg logy logk logl yr*, fe

* ---------------------------------

* ----------- 模型的筛选 ----------

* ---------------------------------

* 固定效应模型还是Pooled OLS？

xtreg logy logk logl yr*, fe /*Wald 检验*/

qui tab id, gen(dum) /*LR检验*/

reg logy logk logl /*POLS*/

est store m_ols

reg logy logk logl dum*,nocons

est store m_fe

lrtest m_ols m_fe

est table m_*, b(%6.3f) star(0.1 0.05 0.01)

* RE vs Pooled OLS？

* H0: Var(u) = 0

* 方法一：B-P 检验

xtreg logy logk logl, re

xttest0

* FE vs RE?

* y_it = u_i + x_it*b + e_it

*--- Hausman 检验 ---

* 基本思想：如果 Corr(u_i,x_it) = 0, Fe 和 Re 都是一致的，但Re更有效

* 如果 Corr(u_i,x_it)!= 0, Fe 仍然有效，但Re是有偏的

* 基本步骤

***情形1：huasman为正数

xtreg logy logk logl, fe

est store m_fe

xtreg logy logk logl, re

est store m_re

hausman m_fe m_re

*** 情形2：

qui xtreg logy h inv gov open,fe

est store fe

qui xtreg logy h inv gov open,re

est store re

hausman fe re

* Hausman 检验值为负怎么办？

* 通常是因为RE模型的基本假设 Corr(x,u_i)=0 无法得到满足

* 检验过程中两个模型的方差-协方差矩阵都采用Fe模型的

hausman fe re, sigmaless

* 两个模型的方差-协方差矩阵都采用Re模型的

hausman fe re, sigmamore

*== 为何有些变量会被drop掉？

use nlswork.dta, clear

tsset idcode year

xtreg ln_wage hours tenure ttl_exp, fe /*正常执行*/

* 产生种族虚拟变量

tab race, gen(dum_race)

xtreg ln_wage hours tenure ttl_exp dum_race2 dum_race3, fe

* 为何 dum_race2 和 dum_race3 会被 dropped ?

* 固定效应模型的设定：y_it = u_i + x_it*b + e_it (1)

* 由于个体效应 u_i 不随时间改变，

* 因此若 x_it 包含了任何不随时间改变的变量，

* 都会与 u_i 构成多重共线性，Stata会自动删除之。

*******异方差、序列相关和截面相关问题

* ---------------- 简 介 -------------

* y_it = x_it*b + u_i + e_it

*

* 由于面板数据同时兼顾了截面数据和时间序列的特征，

* 所以异方差和序列相关必然会存在于面板数据中；

* 同时，由于面板数据中每个截面（公司、个人、国家、地区）之间还可能存在内在的联系， * 所以，截面相关性也是一个需要考虑的问题。

*

* 此前的分析依赖三个假设条件：

* （1） Var[e_it] = sigma^2 同方差假设

* (2) Corr[e_it, e_it-s] = 0 序列无关假设

* (3) Corr[e_it, e_jt] = 0 截面不相关假设

*

* 当这三个假设无法得到满足时，便分别出现 异方差、序列相关和截面相关问题； * 我们一方面要采用各种方法来检验这些假设是否得到了满足；

* 另一方面，也要在这些假设无法满足时寻求合理的估计方法。

* ---------------- 假设检验 -------------

*== 组间异方差检验（截面数据的特征）

* Var(e_i) = sigma_i^2

* Fe 模型

xtreg logy logk logl, fe

xttest3

* Re 模型

* Re本身已经较大程度的考虑了异方差问题，主要体现在sigma_u^2上

*== 序列相关检验

* Fe 模型

* xtserial Wooldridge(2002)，若无序列相关，则一阶差分后残差相关系数应为-0.5

xtserial logy logk logl

xtserial logy logk logl, output

* Re 模型

xtreg logy logk logl, re

xttest1 /*提供多个统计检验量*/

*== 截面相关检验

* xttest2命令 H0: 所有截面残差的相关系数都相等

xtreg logy logk logl, fe

xttest2

* 由于检验过程中执行了SUE估计，所以要求T>N

xtreg logy logk logl if id<6, fe

xttest2

* xtcsd 命令（提供了三种检验方法）

xtreg logy logk logl, fe

xtcsd , pesaran /*Pesaran(2004)*/

xtcsd , friedman /*Friedman(1937)*/

xtreg logy logk logl, re

xtcsd , pesaran

* ----------------- 估计方法 ---------------------

*== 异方差稳健型估计

xtreg logy h inv gov open, fe robust

est store fe_rb

xtreg logy h inv gov open, fe robust

est store fe

* 结果对比

esttab fe_rb fe, b(%6.3f) se(%6.3f) mtitle(fe_rb fe)

*== 序列相关估计

* 一阶自相关 xtregar, fe/re

* 模型： y_it = u_i + x_it*b + v_it (1)

* v_it = rho*v_it-1 + z_it (2)

xtregar logy h inv gov open, fe

est store fe_ar1

xtregar logy h inv gov open,fe lbi /*Baltagi-Wu LBI test*/

---