Stata:三维引力模型介绍与估计-ppmlhdfe-nbreg-reghdfe

发布时间:2022-01-10 阅读 212

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

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

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

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

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

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:左志勇 (中山大学)
邮箱zuozhy6@mail2.sysu.edu.cn


目录


1. 什么是引力模型?

引力模型的思想和概念源自经典物理学中的万有引力定律:

其中,F 为两物体之间的引力大小,其与两物体的质量乘积 m1m2 成正比,与两物体之间的距离平方 r2 成反比, G 为引力常数。

Tinbergen (1962) 和 Poyhonen (1963) 最早将引力模型用于国际贸易的研究,并得出相同的结果:两国双边贸易规模与他们的经济总量成正比,与两国之间的距离成反比。其中 Tinbergen (1962) 列出的引力模型标准式为:

其中,Fij 为 i 国与 j 国之间的贸易总额 (取绝对值),mi 与 mj 为两国经济总量,dij 代表两国之间的地理距离,G 为常数。为使 (2) 式能够进行线性估计,将其转化为对数形式:

相比于国际贸易理论中经典的 Heckscher-Ohlin 模型,引力模型在形式上可以被看作是一个因素模型,我们不妨写出引力模型的面板形式:

其中,ϵij,t 为模型扰动项,且 fij,t=lnFij。在这种情况下 lnmit 和 lnmjt 成为影响 fij,t 的共同因素 (Common Factors),而地理距离变量 dij,t 则控制单一维度的固定效应。

基于 (4) 式,我们可在确定经济规模 (如名义 GDP) 与国家间距离 (如两国首都直线距离) 的衡量指标后,利用简单的 OLS 回归得出系数估计结果,这也是最基本的三因素引力模型。该模型也被用于解释国际贸易以外的问题,如分析城市之间的贸易流量、区域之间的交通流量与人口迁移流量等。

在利用名义 GDP 与两国首都的直线距离作为解释变量之外,许多学者也对引力模型的因素选取进行了拓展,如考虑人口数量 (Lineman,1966)、人均收入 (Berstrand,1989)、制度因素 (Groot 等,2004) 与贸易协定关系 (Piaru 和 Kume,2000) 等作为模型的解释变量。

进一步地,如果我们认为 (4) 式是对现实经济结构的正确表述,那么 lnmit 和 lnmjt 并不需要被确定为某一数值变量 (如该国的经济规模等),而是可以被替换为国家-时间固定效应。同时 lndij 也可以与 (4) 式中的常数项相结合,成为国家对固定效应。通过上述过程,我们可以得到如下三维固定效应引力模型:

2. 三维引力模型的系数估计

2.1 高维固定效应面板 OLS 估计

直观上理解,我们以 (5) 式为基础进行高维固定效应面板的 OLS 回归似乎并无不妥,然而:

  • 考虑到 fij,t=lnFij,当某两国之间的贸易总额为 0 时,取自然对数后贸易数据会出现缺失,以致于估计结果出现样本选择偏差。尽管可以采取 fij,t=ln(Fij+1) 来解决此问题,但是这被证明估计量是不一致的 (Santos Silva 和 Tenreyro,2006)。并且如果零值很多,那么因变量微小的调整会导致模型估计系数和解释力变化很大;
  • 被解释变量存在着分布离散的特征,与一般 OLS 回归形式依赖样本服从正态分布的假设不同,这可能会造成较严重的异方差问题,从而导致参数估计失去有效性且影响显著性判断。

基于以上理由,我们考虑采用其他方法对引力模型进行系数估计。

2.2 偏移量估计

为估计 θit 和 θjt,对 (5) 式进行如下变换:

令上标 ~ 表示变量相对于其时间序列均值的偏离,例如 f~ij,t=fij,t1Tt=1Tfij,t,则由 (5) 式减去 (6) 式得到:

(7) 式再对 i 取截面平均:

将 (7) 式减去 (8) 式:

当 m 的值足够大时,令上式再对 j 取截面平均:

式中 θ^it=θit~1ni=1nθit~,表示 θit~ 在固定时间的条件下相对于其个体序列均值的偏移。作为 θit 的估计量,以 fij,t 为基础表示为:

完全对称地,我们也可以得到:

值得注意的是,(11) 和 (12) 式仅在 (5) 式条件满足,即贸易模型可被表征为一个三维固定效应模型时才成立。如果共同因素数量大于 2,上式所表示的估计就是错误的 (Anderson 和 Wincoop,2003)。当我们需要考察具体因素 (如贸易协定) 对贸易额的影响时,偏移量估计法并不能给我们提供有效的信息。

因此我们认为如何确定因素的数量是一个相当有价值的问题,然而遗憾的是引力模型本身并未对解决这个问题提供任何帮助。另外,偏移量估计法仅提供了模型参数的估计值,无法对变量进行显著性检验。

2.3 高维固定效应面板泊松伪最大似然估计

泊松回归通常意义上是估计「计数数据模型」的标准方法,Gourieroux 等 (1984) 通过放松因变量分布的假设,使泊松回归不再局限于计数数据。作为泊松回归的一种简单拓展,泊松伪极大似然 (Poisson Pseudo Maximum Likelihood,PPML) 估计法被广泛用于估算含有大量零值且存在异方差的贸易数据,即使因变量不服从泊松分布,PPML 回归也能得到一致无偏的估算结果。

Correia 等 (2020) 进一步修正了 PPML 估计法,在存在高维固定效应 (Multiple High-Dimensional Fixed Effects) 的前提下进行 PPML 回归,称为高维固定效应泊松伪极大似然估计法 (PPMLHDFE)。与 PPML 估计法相比,PPMLHDFE 估计法可以更为稳健地检验伪极大似然估计。

2.3.1 IRLS 算法

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

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

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

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

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

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

2.3.2 泊松伪极大似然估计

泊松回归定义如下:

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

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

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

2.3.3 处理 HDFE

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

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

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

2.3.4 加速 HDFE–GIRLS

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

2.3.5 最大似然估计的存在性

Santos Silva 和 Tenreyro (2010,2011) 指出,对于某些数据配置,可能不存在泊松回归的最大似然估计 (MLE),进而导致估计可能无法收敛或收敛到错误的值。在泊松回归的情况下,如果对数似然随着一个或多个系数趋于无穷大而单调增加,则会发生这种情况。

Santos Silva 和 Tenreyro (2010) 认为发生这种情况的主要原因是变量间的多重共线性。为此,他们建议排除有问题的变量。但是,排除哪个回归变量是一个模棱两可的决定,可能会影响其余参数的识别。此外,在具有多个 HDFE 的泊松模型中,该策略甚至不可行。

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

2.3.6 渐进偏误

Weidner 和 Zylkin (2021)指出,当时间维度固定时,以 PPMLHDFE 得出的点估计值和聚类稳健夹心估计量 (通常用于推断) 均具有阶数为 1N 的渐近伴随参数偏差 (Asymptotic Incidental Parameter Bias),且后者为下偏,其中 N 是数据中的贸易主体数。

为此 Weidner 和 Zylkin (2021) 在其论文中提出了一套偏差矫正算法 (具体数学运算过程较复杂,限于篇幅不详细介绍),在 Stata 中我们可以通过 ppml_fe_bias 命令实现对该渐进偏误的矫正。

2.4 高维固定效应负二项回归模型

注意到在泊松回归中,隐含着一个重要假设:

如果贸易额度 (被解释变量) 的期望和方差差距很大——方差明显大于期望,则数据集存在过度分散(Overdispersion)。此时,我们就有必要引入负二项分布回归法 (Negative Binominal Regression,NBREG)。相比于泊松回归,其主要区别在于对 (14) 式的修正:

即在泊松回归的基础上多了一个扰动项 ϵi 来捕捉个体异质性或不可观测部分。可以证明,负二项回归模型的条件期望仍为 E(yi|xi)=μi=exp(xiβ),而条件方差为:

这表明在负二项回归中,条件方差大于条件期望,且条件方差是 α 的增函数,故 α 称为过度分散参数 (Overdisperson Parameter)。当 α0 时,负二项回归就转化为了泊松回归。因此在进行负二项回归后,只要对原假设 H0:α=0 进行检验,即可确定是使用泊松回归还是负二项回归。

在某种意义上,泊松回归和负二项回归的关系有如线性模型中 OLS 与 WLS 的关系。即使数据中存在过度分散,“泊松回归+稳健标准误” 依然提供了对参数及标准误的一致估计,这类似于在异方差的情况下使用 “OLS+稳健标准误”。因此在样本量足够大的情况下,或许泊松模型就已经能够满足研究者的精度要求了。

此外,Kareem 等 (2016) 比较了存在大量零额贸易时不同估计方法的表现,认为基于PPML估计法的一系列估计方法表现得更胜一筹。基于此,实例中我们仍以泊松模型为主进行命令介绍与演示。

3. Stata 命令介绍

3.1 reghdfe 命令

* 命令安装
ssc install reghdfe, replace
ssc install ftools, replace // ftools 为使用 reghdfe 命令的必须工具包

其语法结构与下面的 ppmlhdfe 基本一致,具体可参考连享会推文「reghdfe:多维面板固定效应估计」

3.2 ppmlhdfe 命令

* 命令安装
ssc install ppmlhdfe, replace 
// reghdfe 和 ftools 是使用 ppmlhdfe 的必须工具包
ssc install reghdfe, replace 
ssc install ftools, replace
* 语法结构
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:控制回归表的选项,如置信水平、数字格式等。

3.3 ppml_fe_bias 命令

* 命令安装
ssc install ppml_fe_bias, replace
// 以下四个 packages 是使用 ppml_fe_bias 的必要工具箱
ssc install outreg, replace
ssc install hdfe, replace
ssc install gtools, replace
ssc install rowmat_utils, replace
* 语法结构
ppml_fe_bias depvar [indepvars] [if] [in],  ///
             lambda(varname) i(exp_id) j(imp_id) t(time_id) [options]

其中,

  • lambda(varname):输入代表固定效应之和的变量 varname;
  • i(exp_id) j(imp_id) t(time_id):分别输入代表出口国、进口国与时间的类别变量。

options 选项:

  • bias(name):将系数的误差矫正信息存储在矩阵 name 中;
  • v(name):将误差矫正后的方差矩阵存储在 name 中;
  • beta(name):输入由 ppmlhdfe 得到的系数估计矩阵;
  • approx:采用估计算法计算方差偏差;
  • exact:采用精确算法计算方差偏差。

3.4 nbreg 命令

* 语法结构
nbreg depvar [indepvars] [if] [in] [weight] [, nbreg_options]

该命令为 Stata 内置命令,语法结构与基础的 reg 命令相似,在处理高维固定效应模型时需采取手动设定虚拟变量的方式。

4. Stata 应用实例与方法比较

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

4.1 高维固定效应面板泊松伪极大似然估计

. lxhuse example_trade_fta.dta, clear
. keep if category=="TOTAL"
. egen imp = group(isoimp)
. egen exp = group(isoexp)
. ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) cluster(imp#exp) d nolog

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 | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         fta |      0.192      0.042     4.59   0.000        0.110       0.275
       _cons |     16.457      0.022   757.32   0.000       16.414      16.500
------------------------------------------------------------------------------
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

. // 提取固定效应之和与系数估计矩阵并进行偏差矫正
. predict lambda
. matrix beta = e(b)
. ppml_fe_bias trade fta, i(exp) j(imp) t(year) lambda(lambda) beta(beta)

  Adjusted SEs
  .0527904093
  bias corrections (to be subtracted from original coefficients)
  .0060619547
note: beta matrix will be shortened to the same length as the number of x-variables
-------------------------------------------------------------
        original      bias     adjusted SEs  bias-corrected 
-------------------------------------------------------------
fta    0.1924455   0.0060620   0.0527904      0.1863835    
      (0.0419527)                            (0.0527904)*** 
-------------------------------------------------------------
Standard errors clustered by pair, using a local de-biasing adjustment
to account for estimation noise in the exp-year and imp-year fixed effects.
* p<0.10; ** p<0.05; *** p<0.01

利用 ppmlhdfe 命令进行估计,并采用 ppml_fe_bias 进行偏差矫正,结果显示贸易协定变量 fta 的系数估计修正值为 0.186,修正稳健标准误为 0.053,对国家间贸易额在 1% 水平下有显著正向影响。

4.2 高维固定效应面板负二项回归估计与 OLS 估计

. * 采用负二项回归时,需采用创建交叉项的虚拟变量的方式控制相应的固定效应
. egen impexp=group(imp exp)
. nbreg trade fta i.(imp#year) i.(exp#year) i.impexp vce(cluster impexp)
. reghdfe trade fta, absorb(imp#year exp#year imp#exp) cluster(imp#exp)

在实际操作过程中,不难发现命令 nbreg 所需的运算时间远长于其他两种估计方法,在处理大样本数据时这将会成为 nbreg 一个相对而言较明显的缺点。将三种估计方式的结果汇总如下:

                        Regression Table
---------------------------------------------------------
       |       (1)             (2)             (3)
       |     PPMLHDFE         NBREG          REGHDFE
-------+-------------------------------------------------
   fta |      0.186***	     0.156***       2.34e+06***
       |      (3.53)         (2.80)           (2.96)
-------+-------------------------------------------------
     N |       5950           5950             5950
  r2_a |                                       0.84
  r2_p |       0.99           0.11	
---------------------------------------------------------

表格中 fta 一行内括号中的数字代表显著性检验值,结果显示无论选取何种估计方式,贸易协定变量均在 1% 水平上对两国间贸易额具有显著正向影响。但检验值可能存在一定差别,因此在变量作用相对不显著的情况下,选取不同的估计方式可能会对实证结论产生影响。

需注意的是,由于 ppmlhdfenbreg 采用的均是伪极大似然估计法,故仅存在伪 R2 可作为判断模型拟合效果的指标。在这一例子中,ppmlhdfe 的估计准确性明显高于 nbreg

5. 参考资料

  • Correia S, Guimaraes P, Zylkin T, et al. Fast Poisson estimation with high-dimensional fixed effects[J]. 2020. -PDF-
  • Sul D. Panel Data Econometrics: Common Factor Analysis for Empirical Researchers[M]. 2019.
  • Correia, S. Linear Models with High-Dimensional Fixed Effects: An Efficient and Feasible Estimator. Working Paper. 2016. -PDF-
  • Martin Weidner, Thomas Zylkin. Bias and consistency in three-way gravity models. Journal of International Economics[J]. 2021. -PDF-
  • Kareem F O, Martinez-Zarzoso I, Brümmer B. Fitting the gravity model when zero trade flows are frequent: a comparison of estimation techniques using Africa's trade data[R]. 2016. -PDF-
  • 陈强. 高级计量经济学及Stata应用(第二版)[M]. 高等教育出版社, 2014.

6. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 多维 高维, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

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

课程主页

课程主页

关于我们

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

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

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

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

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