reghdfe：多维面板固定效应估计

Posted 2023-04-09 a useful man

reghdfe：多维面板固定效应估计| 连享会主页

实证分析中，我们经常需要控制各个维度的个体效应，以便尽可能减轻 遗漏变量 导致的偏误。在最常用的二维面板数据中，我们通常会采用 xtreg y x i.year, fe 的形式来控制 公司个体效应 和 年度效应。然而，在有些情况下，我们需要对三维甚至更高维度的数据进行分析 (例如，公司-年度-高管，省份-城市-行业-年度)，此时，一方面要考虑估计的可行性，另一方面还需兼顾计算速度问题。

本文介绍的 reghdfe 命令可以很好地达成上述目的。reghdfe 主要用于实现多维固定效应线性回归。该命令类似于 areg 及 xtreg,fe，但允许引入多维固定效应。此外，该命令在运行速度方面远远优于 areg 及 xtreg， 因此倍受研究者青睐。

本文对该命令的介绍基于 A Feasible Estimator for Linear Models with Multi-Way Fixed Effects (Correia, 2016）。

1. 命令的安装

我们可以使用 ssc install 命令安装最新版本的 reghdfe 相关程序文件。需要特别说明的是，该命令依赖于最新版的 gtools 命令，因此，需要同时安装后者，否则执行 reghdfe 时可能会提示 错误信息。

ssc install gtools, replace  
ssc install reghdfe, replace // 安装最新版命令

可以使用如下命令查看你使用的版本是否为最新版本 (如下是写作本文stata版本信息)：

. which ftools
D:\\stata15/ado\\plus\\f\\ftools.ado
*! version 2.37.0 16aug2019

. which reghdfe
D:\\stata15/ado\\plus\\r\\reghdfe.ado
*! version 5.7.3 13nov2019

2. 命令的语法

该命令的具体语法如下：

reghdfe depvar [indepvars] [if] [in] [weight], absorb(absvars) [options]

其中，

• depvar: 因变量

• indepvars: 解释变量

• absorb(absvars)：引入固定效应

  • 可以包含多维固定效应，即 absorb (var1,var2,var3,...)。若想保存对某变量的固定效应，则运行命令absorb (var1,var2,FE3=var3) ， 变量 FE3 将保存对 var3 的固定效应估计结果。

  • 可以包含不同效应间的交互影响，即absorb(var1#var2)。

值得注意的是，reghdfe 允许定类变量 (categorical variable) 与连续性变量 (continuous variable) 进行交互，即 absorb(i.var1#c.var2) 。实证中很少引入这样的交互项。但如果对该问题感兴趣，可参考 Duflo (2014) 。

3. 命令的操作

这一部分用两个实证的例子介绍如何运用 reghdfe 。

3.1 估计双重差分的固定效应模型（DID)

该命令可用于估计双重差分的固定效应模型（DID）。估计 DID 模型的三个命令：reg, areg, 以及 xtreg。reghdfe 也可实现同样的估计结果，而且运行速度优于其他命令。

使用的数据请参考之前推文Stata: 双重差分的固定效应模型。该数据模拟的情况为，政策冲击发生在t=14时，对照组为i=1，控制组为i=0 。模型为y=0.3+0.19*i+1.67*d+0.56*i*d+e 。

. set obs 400
. gen firm=_n  ///生成企业数量
. expand 24
. bysort firm: gen t=_n  ///时间跨度设定为24个季度（6年）
. gen d=(t>=14)
. label var d "=1 if post-treatment" ///设定事件冲击发生在第14期
. gen r=rnormal()
. qui sum r, d
. bysort firm: gen i=(r>=r(p50)) if _n==1
. bysort firm: replace i=i[_n-1] if i==. & _n!=1 ///设定处理组和对照组
. drop r
. label var i "=1 if treated group, =0 if untreated group" 
. gen e = rnormal() ///设定随机变量
. label var e "normal random variable"
. gen y = 0.3 + 0.19*i + 1.67*d + 0.56*i*d + e ///模型设置

首先，回顾双重差分模型的设定形式，

其中， 为分组虚拟变量（处理组=1，控制组=0）； 为分期虚拟变量（政策实施后=1，政策实施前=0）；交互项  表示处理组在政策实施后的效应。 与  分别为个体固定效应和时间固定效应。

具体用于估计政策冲击对公司的影响的命令如下。

gen did = i*d ///生成交互项
reghdfe y did, absorb(firm t) vce(cluster firm)

变量 did 即为交互项，其系数为双重差分模型重点考察的处理效应。命令 absorb(firm t) 同时引入了公司固定效应以及时间固定效应。结果如下。

. reghdfe y did, absorb(firm t) vce(cluster firm)
(MWFE estimator converged in 2 iterations)

HDFE Linear regression                            Number of obs   =      9,600
Absorbing 2 HDFE groups                           F(   1,    399) =     175.80
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
                                                  R-squared       =     0.5102
                                                  Adj R-squared   =     0.4875
                                                  Within R-sq.    =     0.0198
Number of clusters (firm)    =        400         Root MSE        =     1.0043

                                 (Std. Err. adjusted for 400 clusters in firm)
------------------------------------------------------------------------------
             |               Robust
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         did |   .5656247   .0426601    13.26   0.000     .4817581    .6494914
       _cons |   1.143579   .0084565   135.23   0.000     1.126954    1.160204
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+

 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
        firm |       400         400           0    *|
           t |        24           0          24     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

 

3.2 估计多维固定效应的线性模型（复制一篇 AER 论文）

这一小节将介绍如何运用 reghdfe 估计多维固定效应的线性模型。American Economic Review一篇文章，The Costs of Patronage: Evidence from the British Empire (Xu, 2018)， 提供的可供复制的代码中出现了大量 reghdfe 命令。本小节介绍该作者如何用 reghdfe 命令输出其文章Table 2第六列的结果。

Source: Xu, G. (2018). The Costs of Patronage: Evidence from the British Empire. American Economic Review, 108 (11): 3170-98.

作者在这篇文章中想要探究 任命制 (patronage) 对英国 政治体系 的影响。具体于 Table2，作者意图研究社会联系（social connections) 是否会影响政府官员的工资水平。Table 2中，第六列所估计的回归为：

其中，

• 为政府官员  于时间  在  州执政时的对数工资水平。Stata 命令中，该变量名为 log_salary_governor_gbp

•  为虚拟变量（Dummy Variable)，当政府官员与其上任官员存在社会联系时，该变量取1。如若不然，则取0。社会联系包括：共同祖先，贵族身份以及教育背景。Stata 命令中，该变量名为 connected

•  为政府官员固定效应。该部分的设置为了解决政府官员的异质性 (heterogeneity) 问题。例如，具有较强能力的政府官员更有可能建立更多的社会关系。Stata 命令中, aid 为不同官员的 unique ID 变量。

•  为政府官员执政时长固定效应。设置该部分是因为，执政时间的长短可能也会对社会关系产生影响。Stata 命令中, duration 为官员执政时长变量。

•  为控制变量。作者选用了执政者在历史上执政过的州的数目。Stata 命令中，该变量名为 no_colonies 。

•  为年份固定效应。该部分的设置是为了吸收执政者们在不同时期受到的共同时间冲击。Stata 命令中, year 为年份变量。

•  为残差。作者使用了聚类标准误的方法。

该回归的原假设为，: 社会联系 (connected) 与政府官员的工资水平 (log_salary_governor_gbp) 无关。若 connected 的系数  不显著，则不拒绝原假设。若  显著，则拒绝原假设，并可以判定社会联系对政府官员的工资水平显著相关。

用 Stata 实现该回归的命令如下。

reghdfe log_salary_governor_gbp no_colonies connected, ///
        absorb(aid year duration) vce(cluster bilateral)

其中，absorb(aid year duration) 同时引入了官员固定效应、时间固定效应以及执政时长固定效应。

命令运行后的结果如下所示。数据请于 AER 官网 「下载」。

. quietly use "analysis.dta", replace
. reghdfe log_salary_governor_gbp no_colonies connected, ///
          absorb(aid year duration) vce(cluster bilateral)
(MWFE estimator converged in 26 iterations)

HDFE Linear regression                            Number of obs   =      3,510
Absorbing 3 HDFE groups                           F(   2,   1517) =      25.45
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
                                                  R-squared       =     0.9255
                                                  Adj R-squared   =     0.9109
                                                  Within R-sq.    =     0.0978
Number of clusters (bilateral) =      1,518       Root MSE        =     0.2374

                          (Std. Err. adjusted for 1,518 clusters in bilateral)
------------------------------------------------------------------------------
             |               Robust
log_salary~p |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 no_colonies |   .2234767   .0347473     6.43   0.000     .1553189    .2916346
   connected |   .0972969   .0355508     2.74   0.006     .0275628    .1670309
       _cons |   7.485619    .065766   113.82   0.000     7.356617    7.614621
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
         aid |       456           0         456     |
        year |       110           1         109     |
    duration |         7           1           6    ?|
-----------------------------------------------------+
? = number of redundant parameters may be higher

上述结果表明，变量 connected 的系数为 0.097, 标准误为 0.036。这说明该变量在 1% 的水平上显著大于 0 。其经济学含义为，与上一任官员存在社会联系的官员，相较于无社会联系的官员，工资水平要高出 9.7%。也就是说，官员的工资水平和其社会关系显著相关。

4.结语

这篇推文主要介绍了如何在实证中运用 reghdfe。具体而言，本推文列举了两个例子。其一，为运用该命令对 DID 模型进行估计。其二，为运用该命令进行多维固定效应线性模型的估计。

文献来源

• Correia, S. (2016). Linear Models with High-Dimensional Fixed Effects: An Efficient and Feasible Estimator, Working Paper. [PDF]

• Duflo, E. (2004). The medium run effects of educational expansion: Evidence from a large school construction program in Indonesia. Journal of Development Economics, 74(1), 163-197. [PDF]

• Xu, G. (2018). The Costs of Patronage: Evidence from the British Empire. American Economic Review, 108 (11): 3170-98. [PDF]  

从固定效应空间面板模型估计结果中提取AIC

我正在使用Millo和Piras（2012）的“splm”软件包来估算具有国家（个人）固定效应的SPatial Durbin模型。

 ks = log(spa.sak$index_ret+2) ~ log(spa.sak$lagindex+2) + log(spa.sak$cred_def_rate+2) +   log(spa.sak$ch_in_ex_rate+2) +  
            log(spa.sak$un_inf+2) + log(spa.sak$gdp_growth+2) + log(spa.sak$cha_in_int_rate+2) + country

    ERVM<- spml(formula = ks, data = spa.sak, index = c('id','months'), listw = ERV, model = "within", lag = TRUE, effect = "individual", spatial.error = "none") 
> summary(ERVM)
Spatial panel fixed effects lag model


Call:
spml(formula = ks, data = spa.sak, index = c("id", "months"), 
    listw = ERV, model = "within", effect = "individual", lag = TRUE, 
    spatial.error = "none")

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-0.783343 -0.028334 -0.003056  0.020040  1.199906 

Spatial autoregressive coefficient:
         Estimate Std. Error t-value Pr(>|t|)
lambda -0.0040442  0.0061785 -0.6546   0.5127

Coefficients:
                                    Estimate  Std. Error  t-value  Pr(>|t|)    
log(spa.sak$lagindex + 2)        -2.5577e-01  1.6593e-02 -15.4142 < 2.2e-16 ***
log(spa.sak$cred_def_rate + 2)    1.3427e-02  2.3244e-02   0.5777    0.5635    
log(spa.sak$ch_in_ex_rate + 2)   -5.1946e-04  2.2338e-03  -0.2325    0.8161    
log(spa.sak$un_inf + 2)           1.1360e-01  2.2112e-01   0.5138    0.6074    
log(spa.sak$gdp_growth + 2)       3.1093e+00  2.1339e+00   1.4571    0.1451    
log(spa.sak$cha_in_int_rate + 2) -1.4148e+01  3.1507e+00  -4.4905 7.105e-06 ***
countryEGYPT                      1.0659e-02  1.7940e-02   0.5942    0.5524    
countryGHANA                      2.3214e-02  2.0567e-02   1.1287    0.2590    
countryKENYA                      1.6923e-02  1.9599e-02   0.8634    0.3879    
countryMALAWI                     1.0190e-02  3.2028e-02   0.3182    0.7504    
countryMAURITIUS                  2.5843e-03  1.3530e-02   0.1910    0.8485    
countryMOROCCO                    4.2243e-03  1.4848e-02   0.2845    0.7760    
countryNAMIBIA                   -1.0416e-02  1.4408e-02  -0.7229    0.4697    
countryNIGERIA                    1.2794e-02  1.8327e-02   0.6981    0.4851    
countrySOUTH AFRICA              -9.4671e-04  1.3665e-02  -0.0693    0.9448    
countryTANZANIA                   1.9698e-02  1.9822e-02   0.9937    0.3204    
countryTUNISIA                    1.4579e-02  1.4895e-02   0.9788    0.3277    
countryUGANDA                     3.1850e-02  2.1042e-02   1.5137    0.1301    
countryZAMBIA                     3.3671e-02  2.1202e-02   1.5881    0.1123    
countryZIMBABWE                   1.4299e-02  3.1798e-02   0.4497    0.6529    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 > str(ERVM)
    List of 15
     $ coefficients : Named num [1:21] -0.004044 -0.255767 0.013427 -0.000519 0.113601 ...
      ..- attr(*, "names")= chr [1:21] "lambda" "log(spa.sak$lagindex + 2)" "log(spa.sak$cred_def_rate + 2)" "log(spa.sak$ch_in_ex_rate + 2)" ...
     $ errcomp      : NULL
     $ vcov         : num [1:21, 1:21] 3.82e-05 0.00 0.00 0.00 0.00 ...
     $ spat.coef    : Named num -0.00404
      ..- attr(*, "names")= chr "lambda"
     $ vcov.errcomp : NULL
     $ residuals    : Named num [1:3390] -0.1195 0.0395 0.0151 -0.0186 0.0285 ...
      ..- attr(*, "names")= chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
     $ fitted.values: num [1:3390, 1] 0.681 0.712 0.704 0.725 0.719 ...
      ..- attr(*, "dimnames")=List of 2
      .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. ..$ : NULL
      ..- attr(*, "names")= chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
     $ sigma2       : num 0.0192
     $ type         : chr "fixed effects lag"
     $ model        :'data.frame':  3390 obs. of  21 variables:
      ..$ y                               : num [1:3390] 0.537 0.731 0.697 0.685 0.728 ...
      ..$ log.spa.sak.lagindex...2.       : num [1:3390] 0.851 0.71 0.752 0.674 0.731 ...
      ..$ log.spa.sak.cred_def_rate...2.  : num [1:3390] 2.83 2.83 2.83 2.83 2.83 ...
      ..$ log.spa.sak.ch_in_ex_rate...2.  : num [1:3390] 0.696 0.69 0.711 0.768 0.695 ...
      ..$ log.spa.sak.un_inf...2.         : num [1:3390] 0.697 0.694 0.693 0.692 0.694 ...
      ..$ log.spa.sak.gdp_growth...2.     : num [1:3390] 0.693 0.693 0.693 0.693 0.693 ...
      ..$ log.spa.sak.cha_in_int_rate...2.: num [1:3390] 0.693 0.693 0.693 0.693 0.693 ...
      ..$ countryEGYPT                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryGHANA                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryKENYA                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMALAWI                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMAURITIUS                : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMOROCCO                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryNAMIBIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryNIGERIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countrySOUTH.AFRICA             : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryTANZANIA                 : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryTUNISIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryUGANDA                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryZAMBIA                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryZIMBABWE                 : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
     $ call         : language spml(formula = ks, data = spa.sak, index = c("id", "months"), listw = ERV, model = "within", effect = "individual| __truncated__
     $ logLik       : num 1895
     $ method       : chr "eigen"
     $ effects      : chr "spfe"
     $ res.eff      :List of 2
      ..$         :List of 7
      .. ..$ res.sfe   : num [1:15, 1] 0.00264 -0.00152 0.0017 0.0023 0.0106 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:15] "1" "2" "3" "4" ...
      .. .. .. ..$ : NULL
      .. ..$ res.se.sfe: Named num [1:15] 2.55 2.55 2.55 2.55 2.55 ...
      .. .. ..- attr(*, "names")= chr [1:15] "1" "2" "3" "4" ...
      .. ..$ intercept : num [1, 1] 8.43
      .. ..$ res.se.con: num [1, 1] 2.55
      .. ..$ xhat      : num [1:3390, 1] 0.681 0.712 0.704 0.725 0.719 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. .. .. ..$ : NULL
      .. ..$ N.vars    : int 35
      .. ..$ res.e     : num [1:3390, 1] -0.1195 0.0395 0.0151 -0.0186 0.0285 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. .. .. ..$ : NULL
      ..$ res.corr: NULL
     - attr(*, "class")= chr "splm"

除了以下两个，我成功获得估计结果。

1）如上所示，我只能得到lambda和eta系数的值。我无法得到 heta估计结果或与空间滞后解释变量相关的值，这些变量显示了时空关系。如何提取它？

2）我无法提取AIC和BIC进行模型比较。在第一种情况下，我的试验没有成功，而对于第二种情况，我试过了

AICsplm(ERVMX, criterion="AIC")

但是我收到了一个错误。

我该如何解决这些问题？

更新：我已经解决了这个问题，请参阅下面的内容。

答案

提前构建响应和预测变量的空间滞后，并添加到公式并运行。

就AIC而言，需要从github download AICsplm.R下载AICsplm.R来运行代码。

source(AICsplm.R) 
AICsplm(ERVM)