引力模型-高维固定效应面板泊松模型

发布时间:2021-03-20 阅读 668

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装命令如下:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

作者:李胜胜 (安徽大学)
邮箱lishengsheng2016@126.com

编者按:本文摘译自下文,特此致谢!
Source:Correia S, Guimarães P, Zylkin T. Fast Poisson estimation with high-dimensional fixed effects[J]. The Stata Journal, 2020, 20(1): 95-115. -PDF-


目录


1. 引言

泊松回归现已成为「计数数据模型」的标准方法。与普通最小二乘法一样,泊松回归的一致性估计也需要正确说明因变量的条件均值假设 (Gourieroux 等,1984)。在该条件下,泊松回归就能转化为泊松伪最大似然 (PPML) 回归。Gourieroux 等 (1984) 通过放松因变量分布的假设,使泊松回归不再局限于计数数据。也因此,泊松伪最大似然 (PPML) 回归可以应用于任何非负的因变量。

实际上,当非负数据中存在较多 0 的时候,PPML 似乎是最好的选择。并且,这种情况常发生于以下几个领域,例如公司研发支出、专利引用计数、日常产品销售、医生就诊次数、公司信贷量、拍卖竞拍者数量和跨地区通勤者数量。

但是,学者更多使用线性对数回归,即使 PPML 更合理。一种可能的解释是,学者可以轻松估计控制多个固定效应的线性回归。尤其是伴随着大型面板数据集的可用性不断提高,以及具有高维固定效应 (HDFE) 的线性回归模型估计技术的进步,使得学者能够控制更多异质性来源。

需要指出的是,PPML 也可以轻松的实现 HDFE。为此,本文将介绍新命令 ppmlhdfe。相对于用于现有 HDFE 非线性估计算法,ppmlhdfe 消除了一些不必要的步骤,提高了参数的计算速度。

2. 估计方法

2.1 IRLS 算法

GLM 是基于 Nelder 和 Wedderburn (1972) 引入的指数分布族的一类回归模型,主要包括流行的非线性回归模型,例如 Logit、Probit、Cloglog 和 Poisson。Hardin 和 Hilbe (2018) 将指数族定义如下:

其中,a()b() 和 c()是特定函数,ϕ 和 θ是参数。对于这些模型而言:

给定一组具有 n 个独立的观察值,每个观察值均由 i 索引,可以通过链接函数 g(.) 将期望值与一组协变量 (xi) 关联。更具体地说,假设下式:

则 GLM 的似然函数可以写为:

通过求解 (伪) 似然最大化的一阶条件,得到 β 的估计。应用高斯-牛顿算法和期望 Hessian 函数导出修正方程为:

其中 X 是解释变量矩阵,W(r1) 是权重矩阵,Z(r1) 是因变量的变换,r 是迭代索引。式 (1) 表明 β 的估计值是通过加权最小二乘法递归而获得,这种方法称之为 IRLS。

2.2 泊松回归

泊松回归定义如下:

实现 IRLS 的回归权重可以简化为:

而中间回归的因变量则为:

因此,在 (1) 中实现 IRLS 更新回归仅需要计算在先前迭代中获得拟合值xiβ(r1)的向量即可。

处理 HDFE

在 HDFE 下,IRLS 的困难是 X 可能包含很多固定效应,导致无法直接计算 (XW(r1)X)。解决方案是使用一个替代的公式,即只估计非固定效应协变量的系数 (例如,δ),从而降低维数。由于式 (1) 是加权线性回归,因此可以依靠 FWL 定理来探讨固定效应。这意味着代替式 (1),可以使用下式:

其中,X~ 和 Z~ 分别是主协变量矩阵 X 和因变量 z 的变换后的加权。此外,FWL 定理还暗示从式 (4)计算的残差与从式 (1) 计算的残差相同。这意味着可以使用以下方法对 W 和 z 执行所需的更新:

其中,e(r) 是一个向量,该向量收集使用式 (4) 计算的残差。然后,与原始 IRLS 循环一样,W(r) 和 z(r) 的新值直接从式 (2) 和 (3) 开始。其次,这也意味着,一旦 δ(r) 收敛到正确的估计值 δ^,式 (4) 中加权最小二乘回归的估计方差-协方差矩阵将成为 δ^ 正确的方差-协方差矩阵。

加速 HDFE–GIRLS

命令 poi2hdfe 可以实现式 (4),同时使用 reghdfe 作为运行 HDFE 加权最小二乘回归。这是一个密集型计算过程,需要在每次 IRLS 迭代中估算 HDFE 回归模型。但是,ppmlhdfe 中有多种变通办法,可以使其效率更高。例如,ppmlhdfe 直接嵌入 reghdfe 的 Mata 中。因此利用了这样一个事实,它们在每个 IRLS 迭代中保持不变,某些计算只需要执行一次。但最显著的速度改进来自于对标准 HDFE–IRLS 算法的改进,该算法旨在减少对 reghdfe 的调用次数。

最大似然估计的存在性

Santos Silva 和 Tenreyro (2010,2011) 指出,对于某些数据配置,可能不存在泊松回归的最大似然估计 (MLE),进而导致估计可能无法收敛或收敛到错误的值。在泊松回归的情况下,如果对数似然随着一个或多个系数趋于无穷大而单调增加,则会发生这种情况。Santos Silva 和 Tenreyro (2010) 认为发生这种情况的主要原因是变量间的多重共线性。为此,他们建议排除有问题的变量。但是,排除哪个回归变量是一个模棱两可的决定,可能会影响其余参数的识别。此外,在具有多个 HDFE 的泊松模型中,该策略甚至不可行。

Correia 等 2019 讨论了各种 GLM 模型估计中的必要条件和充分条件,并表明在泊松回归情况下,如果从样本中删除一些观察值,可能得到总 MLE 估计。这些单独的观测值不传递估计过程的相关信息,可以安全地丢弃。同时,删除这些观察值后,某些回归变量将产生共线性,因此也必须删除。此外,Correia 等 2019 提出了一种识别分离观察结果的方法,并且即使在 HDFE 存在情况下,也可以成功运行。

3. ppmlhdfe 命令

3.1 命令安装:

ssc install ppmlhdfe, replace

3.2 语法结构:

ppmlhdfe depvar [indepvars] [if] [in] [weight] ,  ///
         [absorb(absvars)] [options]
  • depvar:被解释变量;

  • indepvars:解释变量;

  • absorb(absvars):要吸收的分类变量 (固定效应),也允许单独的斜率;

  • absorb(..., savefe):使用 hdfe # 保存所有固定效应估计值。

  • [options] 选项:

    • exposure(varname):在系数约束为 1 的模型中包含 ln(varname);
    • offset(varname):在系数约束为 1 的模型中包含 varname;
    • d(newvar):将固定效应之和另存为newvar;
    • d:如上,但变量将另存为 _ppmlhdfe_d;
    • vce(vcetype):vcetype 可以是 robust 或聚类;
    • verbose(#):显示调试信息量;
    • nolog:隐藏迭代日志;
    • tolerance(#):收敛标准,默认为 tolerance(1e-8)
    • guess(string):设置用于设置初始值的规则;
    • eform:报告指数系数;
    • irr:eform 的同义词;
    • separation(string):用于删除分离的观测值及其相关回归变量的算法;
    • maxiteration(#):指定最大迭代次数;
    • keepsingletons:不要删除单例组;
    • version:报告 ppmlhdfe 的版本号和日期;
    • display_options:控制回归表的选项,如置信水平、数字格式等。

4. Stata Stata 范例

4.1 范例 1:模拟数据

从一个简单的例子开始,我们来说明 ppmlhdfe 方法处理 MLE 估计不存在时的优势。如前所述,ppmlhdfe 注意识别分离的观测值,然后以保证存在有意义的 MLE 方式限制样本。

以下数据包括六个观察值和三个解释变量:

input y x1 x2 x3
0       1       2       1
0       0       0       2
0       2       3       3
1       1       2       4
2       2       4       5
3       1       2       6
end

如果试图用 glm 命令估计泊松回归,Stata 就无法收敛。poisson 命令可以得到三个解释变量系数的估计值,但是会发现估计结果并不可靠。

. poisson y x1 x2 x3, nolog

Poisson regression                       Number of obs     =          6
                                         LR chi2(3)        =       8.89
                                         Prob > chi2       =     0.0308
Log likelihood = -4.0415302              Pseudo R2         =     0.5237

-----------------------------------------------------------------------
    y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------+----------------------------------------------------------------
   x1 |  -31.29521   8467.059    -0.00   0.997    -16626.43    16563.83
   x2 |   15.84339   4233.529     0.00   0.997    -8281.722    8313.408
   x3 |   .7970409   .4608654     1.73   0.084    -.1062388     1.70032
_cons |  -4.032453   2.868066    -1.41   0.160    -9.653759    1.588853
-----------------------------------------------------------------------

使用命令 ppml 识别存在问题的数据,并删除变量 x3。但是,拟合回归仍然有问题。

. ppml y x1 x2 x3

note: checking the existence of the estimates
note: starting ppml estimation

Iteration 1:   deviance =  5.984675
Iteration 2:   deviance =  5.625005
Iteration 3:   deviance =  5.538839
Iteration 4:   deviance =  5.507367
Iteration 5:   deviance =   5.49579
Iteration 6:   deviance =   5.49153
Iteration 7:   deviance =  5.489963
Iteration 8:   deviance =  5.489387
Iteration 9:   deviance =  5.489175
Iteration 10:  deviance =  5.489097
Iteration 11:  deviance =  5.489068
Iteration 12:  deviance =  5.489058
Iteration 13:  deviance =  5.489054
Iteration 14:  deviance =  5.489052
Iteration 15:  deviance =  5.489052
Warning:  variance matrix is nonsymmetric or highly singular

Number of parameters: 3
Number of observations: 6
Number of observations dropped: 0
Pseudo log-likelihood: -6.5473014
R-squared: .3399843
WARNING: The model appears to overfit some observations with y=0
-----------------------------------------------------------------------
      |               Robust
    y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------+----------------------------------------------------------------
   x1 |  -32.31708          .        .       .            .           .
   x2 |   16.58591          .        .       .            .           .
_cons |  -.8167293          .        .       .            .           .
-----------------------------------------------------------------------
Number of regressors dropped to ensure that the estimates exist: 1
Dropped variables:  x3
Option strict is off

如果使用 strict 选项运行,则 ppml 会删除所有回归变量。而 ppmlhdfe 方法不同,与其寻找有问题的回归变量,不如寻找有问题的观察值。在这种情况下,ppmlhdfe 删除了第三个观察值,然后丢弃 x2 以避免完全多重共线,此时估计结果会更合理。

. ppmlhdfe y x1 x2 x3, nolog

(simplex method dropped 1 separated observation)
note: 1 variable omitted because of collinearity: x2
Converged in 6 iterations and 6 HDFE sub-iterations (tol = 1.0e-08)

PPML regression                            No. of obs      =          5
                                           Residual df     =          2
                                           Wald chi2(2)    =      50.78
Deviance             =  .4775093816        Prob > chi2     =     0.0000
Log pseudolikelihood = -4.041530113        Pseudo R2       =     0.4532
-----------------------------------------------------------------------
      |               Robust
    y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------+----------------------------------------------------------------
   x1 |   .3914642   .1733025     2.26   0.024     .0517976    .7311309
   x2 |          0  (omitted)
   x3 |   .7969293   .1582404     5.04   0.000     .4867839    1.107075
_cons |  -4.031679   1.119577    -3.60   0.000    -6.226011   -1.837348
-----------------------------------------------------------------------

4.2 范例 2:Stata 手册例子

接下来,我们使用固定效应选项的命令 xtpoisson 来复制 Stata 16 手册中的示例 1。本示例使用的数据为 ships.dta,并以此来估计多个回归变量上船舶事故数量的泊松回归。它将可变船舶作为固定效应,即控制五种类型的船舶。通过控制 exposure(service) 估计回归,并将系数报告为发生率。

. webuse ships, clear
. xtpoisson acc op_75_79 co_65_69 co_70_74 co_75_79, ///
            exposure(service) irr fe nolog

Conditional fixed-effects Poisson regression    Number of obs     =     34
Group variable: ship                            Number of groups  =      5

                                                Obs per group:
                                                              min =      6
                                                              avg =    6.8
                                                              max =      7

                                                Wald chi2(4)      =  48.44
Log likelihood  = -54.641859                    Prob > chi2       = 0.0000

--------------------------------------------------------------------------
    accident |        IRR   Std. Err.    z    P>|z|   [95% Conf. Interval]
-------------+------------------------------------------------------------
    op_75_79 |   1.468831   .1737218   3.25   0.001   1.164926    1.852019
    co_65_69 |   2.008002   .3004803   4.66   0.000   1.497577    2.692398
    co_70_74 |    2.26693    .384865   4.82   0.000   1.625274    3.161912
    co_75_79 |   1.573695   .3669393   1.94   0.052   .9964273    2.485397
 ln(service) |          1  (exposure)
--------------------------------------------------------------------------

ppmlhdfe 可以获得等效的结果,语法为:

. ppmlhdfe acc op_75_79 co_65_69 co_70_74 co_75_79,  ///
           absorb(ship) exposure(service) irr nolog

Converged in 6 iterations and 6 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs   =     34
Absorbing 1 HDFE group                            Residual df  =     25
                                                  Wald chi2(4) = 111.06
Deviance             =  38.69505154               Prob > chi2  = 0.0000
Log pseudolikelihood = -68.28077143               Pseudo R2    = 0.8083
-----------------------------------------------------------------------
             |             Robust
    accident |   exp(b)   Std. Err.      z    P>|z|        [95% CI]
-------------+---------------------------------------------------------
    op_75_79 | 1.468831   .1484359     3.80   0.000  1.204902  1.790572
    co_65_69 | 2.008002   .2202475     6.36   0.000  1.619572  2.489592
    co_70_74 |  2.26693   .3256501     5.70   0.000  1.710649  3.004107
    co_75_79 | 1.573695   .3117262     2.29   0.022  1.067358  2.320232
       _cons | .0011254   .0001061   -72.03   0.000  .0009356  .0013538
 ln(service) |        1  (exposure)
-----------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
        ship |         5           0           5     |
-----------------------------------------------------+

在这里吸收船舶作为固定效应,有几点值得一提,如预期的那样,变量的估计系数与使用 xtpoisson 命令获得的系数相同。但是,标准误的估计结果有所不同,因为默认情况下,ppmlhdfe 报告的是可靠的标准误。这两个命令显示的对数似然值也不同,命令 xtpoisson 报告条件对数似然值,而 ppmlhdfe 报告实际的泊松对数似然值。

考虑到我们使用的是一个小数据集,可以使用 poisson 命令复制ppmlhdfe 获得的结果,如下所示:

. ppmlhdfe acc op_75_79 co_65_69,    ///
     absorb(ship co_70_74 co_75_79)  ///
     exposure(service) irr nolog

值得注意的是,ppmlhdfe 可以吸收任何类别变量作为固定效应。例如,如果仅对 op_75_79 和 co_65_69 的系数感兴趣,则可以将 ship、co_70_74 和 co_ 75_79 吸收。

. ppmlhdfe acc op_75_79 co_65_69,    ///
     absorb(ship co_70_74 co_75_79)  ///
     exposure(service) irr nolog

Converged in 7 iterations and 23 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs    =      34
Absorbing 3 HDFE groups                           Residual df   =      25
                                                  Wald chi2(2)  =   71.60
Deviance             =  38.69505154               Prob > chi2   =  0.0000
Log pseudolikelihood = -68.28077143               Pseudo R2     =  0.8083
-------------------------------------------------------------------------
             |               Robust
    accident |     exp(b)   Std. Err.      z    P>|z|       [95% CI]
-------------+-----------------------------------------------------------
    op_75_79 |   1.468831   .1484359     3.80   0.000  1.204902  1.790572
    co_65_69 |   2.008002   .2202475     6.36   0.000  1.619572  2.489592
       _cons |   .0015435   .0001204   -83.00   0.000  .0013247  .0017984
 ln(service) |          1  (exposure)
-------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
        ship |         5           0           5     |
    co_70_74 |         2           1           1     |
    co_75_79 |         2           1           1    ?|
-----------------------------------------------------+
? = number of redundant parameters may be higher

4.3 范例 3:贸易引力模型

使用 ppml_ panel_sg 命令提供的辅助数据和示例来拟合引力模型。该数据集包含 35 个国家 1986 年至 2004 年的年度双边贸易数据。目的是估计自由贸易协定变量 fta 对贸易的影响。在本例中,我们希望控制国家对固定效应 (country-pair fixed effects) 和国家时间固定效应 (对于进口国和出口国)。此外,希望标准误聚类在国家对的级别上。

* 数据下载地址:
* https://gitee.com/arlionn/data/blob/master/data01/EXAMPLE_TRADE_FTA.dta

use http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA if category=="TOTAL", clear
egen imp=group(isoimp)
egen exp=group(isoexp)
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) cluster(imp#exp) nolog
Converged in 11 iterations and 35 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                         No. of obs   =   5,950
Absorbing 3 HDFE groups                      Residual df  =   1,189
Statistics robust to heteroskedasticity      Wald chi2(1) =   21.04
Deviance             =  377332502.3          Prob > chi2  =  0.0000
Log pseudolikelihood = -188710931.7          Pseudo R2    =  0.9938

Number of clusters (imp#exp)=      1,190
                 (Std. Err. adjusted for 1,190 clusters in imp#exp)
-------------------------------------------------------------------
       |            Robust
 trade |    Coef.  Std. Err.      z    P>|z|   [95% Conf. Interval]
-------+-----------------------------------------------------------
   fta | .1924455  .0419527     4.59   0.000   .1102197    .2746713
 _cons | 16.45706  .0217308   757.32   0.000   16.41447    16.49965
-------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
    imp#year |       175           0         175     |
    exp#year |       175           5         170     |
     imp#exp |      1190        1190           0    *|
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

因为交易示例包括 HDFE,所以它也为我们提供了一个机会,可以从编程的角度解开 ppmlhdfe 的工作原理和背后运行机制。ppmlhdfe 本身是一个复杂的命令,主要在 Mata 中进行编码,并尽可能利用 reghdfe 现有内部工作原理。但是,用于实现 HDFE–IRLS 基本算法具有很高的可移植性,并且可以使用任何已经为内部转换和标准加权回归提供可靠算法的语言进行编程。

4.4 范例 4:结果输出及展示

本示例将展示如何将 ppmlhdfeesttab 命令结合来构建发布高质量的回归表格。为了构造固定效应的标签,使用了 reghdfe 附带的 estfe 命令。该示例由一个包含 1,000,000 多个观察值的数据组成,其中包含 1991 年至 2009 年间收集的所有 SSC 引文。该数据具有期刊 ID、文章标题、作者人数、两个主要的 JEL 分类代码、出版年份、文章类型、以及每年收集的引用次数。但是,所有观测值中有将近 58% 为零。因变量是引文的数量,我们主要的关注点在于一旦控制了各种可以控制的变量,一篇文章中作者数量是否真的提升了引文数量。

* 数据下载地址:
* https://gitee.com/arlionn/data/blob/master/data01/citations_example.rar

use citations_example, clear
estimates clear
ppmlhdfe cit nbaut, absorb(issn type jel2 pubyear)
eststo
ppmlhdfe cit nbaut, absorb(issn#c.year type jel2 pubyear)
eststo
ppmlhdfe cit nbaut, absorb(issn#year type jel2 pubyear)
eststo
estfe *, labels(type "Article type FE" jel2 "JEL code FE" pubyear "Publication year FE" ///
issn "ISSN FE" issn#c.year "Year trend by ISSN" issn#year "ISSN-Year FE")

esttab, indicate(`r(indicate_fe)',   ///
        labels("Yes" "")) b(3) se(3) ///
        varwidth(25) label compress  ///
        stat(N ll, fmt(%9.0fc %9.2fc)) ///
        se starlevels(* 0.1 ** 0.05 *** 0.01) 

输出结果如下:

----------------------------------------------------------------
                                (1)          (2)          (3)   
                                cit          cit          cit   
----------------------------------------------------------------
nbaut                         0.190***     0.190***     0.189***
                            (0.003)      (0.003)      (0.003)   
Constant                      0.054***     0.081***     0.110***
                            (0.006)      (0.006)      (0.006)   
Article type FE                 Yes          Yes          Yes   
JEL code FE                     Yes          Yes          Yes   
Publication year FE             Yes          Yes          Yes   
ISSN FE                         Yes                             
Year trend by ISSN                           Yes                
ISSN-Year FE                                              Yes   
----------------------------------------------------------------
N                         1,083,701    1,083,701    1,080,051   
ll                        -1.71e+06    -1.69e+06    -1.66e+06   
----------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

结果表明作者数量对被引次数具有显著正影响。

5. 参考文献

  • Gourieroux C, Monfort A, Trognon A. Pseudo maximum likelihood methods: Theory[J]. Econometrica: journal of the Econometric Society, 1984: 681-700. -PDF-
  • Correia S, Guimarães P, Zylkin T. Fast Poisson estimation with high-dimensional fixed effects[J]. The Stata Journal, 2020, 20(1): 95-115. -PDF-

6. 相关推文

Note:产生如下推文列表的命令为:
lianxh 固定效应 泊松
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,400+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh