Stata:局部回归分布估计量-lpdensity

发布时间:2022-01-12 阅读 123

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下载 - 推文合集

作者:陈卓然 (中山大学岭南学院)
邮箱chenzhr25@mail2.sysu.edu.cn

编者按:本文主要摘译自下文,特此致谢!
Source: Cattaneo M D, Jansson M, Ma X. Local regression distribution estimators[J]. Journal of Econometrics, 2021. -PDF-


目录


1. 问题背景

自从 Whitney Newey 开创性地将核平滑应用到非参数或半参数估计中以后,分布函数和密度函数的非参数核估计及其高阶导数在计量经济学中就一直扮演着重要的角色。本文在此基础上进一步研究局部回归分布估计量的大样本性质,以便更好地将其应用于非参/半参数的估计和推断当中。

局部回归分布估计量是通过对一个随机变量 xXR 经验分布函数的局部最小二乘近似得到的,在点 xX 的局部化是通过一个核函数和一个带宽参数来实现的。此外我们通过使用一个有限维基准函数来得到局部函数近似。如果基准函数包含 pN 阶多项式,那么相应的最小二乘系数将会给出分布函数、密度函数、以及更高阶导数 (p1 阶以下) 在 xX 处的估计量。

相比于其他类似的估计量,局部回归分布估计量不需要进行预先数据的平滑,因此避开了调整参数选择问题,并且很容易实施和解释。

2. 局部回归分布估计量

假定 x1,x2,,xn 是随机变量 x 的一族随机样本,x 是定义在实数集 (XR) 上,且具有绝对连续累积分布函数 F() 和相应的勒贝格密度 f() 的随机变量。本文提出一种新的 F()f() 及其导数的非参数估计量。

局部回归分布估计量应用于当 F() 在 x 附近光滑并且当 h 很小时,下式也很小:

其中,R() 是已知的局部基准函数。那么,θ(x)=(F(x),f(x),,f(p1)(x)) 的局部回归估计量为:

其中,Wi=K((xix)/h)/hRi=R(xix)F^i=1nj=1n1(xjxi) 是 xi 处的经验分布函数。

局部回归分布估计量是通过经验分布函数的一个最小二乘映射而得到的,这种映射赋予每一个观测值相等的权重。因此将前述的局部回归估计量应用 L2(F^)-映射,得到局部回归分布估计量为:

3. Stata 应用

本文的案例演示都以 Stata 为例,关于 Python 和 R 语言的应用,详见复制文件「CJM_2021_JOE」

3.1 命令介绍

*命令安装
lxhget lpdensity.pkg, install replace
// 或者 net install lpdensity, from(https://raw.githubusercontent.com/nppackages/lpdensity/master/stata) replace
// 关于 Python 和 R 的命令安装包请参考:https://nppackages.github.io/lpdensity
*命令语法
lpdensity Var [if] [in] [,
              grid(Var) bw(Var or #) p(#) q(#) v(#) kernel(KernelFn) 
              scale(#) nomasspoints
              bwselect(BwMethod) nlocalmin(#) nuniquemin(#) noregularize nostdvar
              cweights(Var) pweights(Var)
              genvars(NewVarName)
              rgrid(Var) rindex(Var) level(#) ciuniform cisimul(#) separator(#)
              plot
              estype(ESOpts) 
              esline_options(ESLineOpts) espoint_options(ESPointOpts)
              citype(CIOpts) 
              ciregion_options(CIRegionOpts) 
              ciline_options(CILineOpts)   ciebar_options(CIEbarOpts)
              histogram hiplot_options(HistOpts)
              graph_options(GraphOpts)
              ]

其中,

  • grid(var):指明密度被估计的网格,默认是选择数据集的 5% 到 95%,步长为 0.05;
  • bw(var or #):指明用作估计的带宽 (可以是包含带宽的变量,也可以是一个数字),当被省略时,带宽将被通过 bwselect(BwMethod) 来计算;
  • p(#):指明用作构建点估计的局部多项式的阶数,默认为 2 阶;
  • q(#):指明用来构建置信区间的局部多项式的阶数,默认为 p(#)+1,通常情况下它应该严格大于 p(#)
  • v(#):指明待估计的分布函数的导数,v(0) 表示分布函数,v(1) 表示密度函数 (默认);
  • kernel(KernelFn):指明用来构建局部多项式估计量的核函数:
    • triangularK(u)=(1|u|)(|u|<=1) 这是默认情形;
    • epanechnikovK(u)=0.75(1u2)(|u|<=1)
    • uniformK(u)=0.5(|u|<=1)
  • scale(#):控制估计量的成倍增加或减少;
  • nomassponts:将不会调整点估计或者标准误即使数据集中有质点;
  • bwselect(BwMethod):指明带宽选择的方法:
    • mse-dpi:每一个网格点的均方误最优带宽;
    • imse-dpi:加总均方误最优带宽;
    • mse-rot:基于高斯推断模型的经验带宽;
    • imse-rot:加总的基于高斯推断模型的经验带宽。
  • nlocalmin:指明在每一个局部邻域内不同观测值的最小个数;
  • noregularize:将每一个局部邻域内不同观测值的最小个数设置为 0;
  • nostdvar:将不会根据带宽选择来标准化数据;
  • cweights(var):指明用于反事实分布构建的权重;
  • pweights(var):指明用于抽样的权重,应该是非负的;
  • genvars(NewVarName):生成一个用来存储估计结果的变量;
  • rgrid(var):指明用来展示结果的网格点;
  • rindex(var):指明用来呈现结果的指数;
  • level(#):控制置信区间的水平;
  • ciuniform:计算统一的置信区间,而非逐点置信区间;
  • cisimul(#):指明用来构建临界值的模拟次数,默认是 2000 次;
  • separator(#):在每 # 个变量之后画一条分割线;
  • plot:如果被指明的话表示点估计和区间估计都将被绘制;
  • estype(ESOpts):指明点估计的绘制形状:
    • line:一条曲线 (默认);
    • points:单独的点;
    • both:点和曲线;
    • none:将不会绘制点估计。
  • esline_options(ESLineOpts):指明twoway line 的相关选项;
  • espoint_options(ESPointOpts):指明twoway scatter 的相关选项;
  • citype(CIOpts):指明置信区间的绘制形状:
    • region:绘制阴影区域;
    • line:绘制上下边界;
    • ebar:误差柱 (error bars);
    • all:以上全部;
    • none:将不会绘制置信区间。
  • ciregion_options(CILineOpts):指明 twoway rarea 的有关选项;
  • ciline_options(CILineOpts):指明 twoway rline 的有关选项;
  • ciebr_options(VIEbarOpts):指明 twoway rcap 的有关选项;
  • histgram:如被指明,将绘制一个直方图;
  • hiplot_options(HistOpts):指明额外的 twoway histogram 的选项;
  • graph_options(GraphOpts):指明额外绘制选项,例如图例和标签等。

3.2 反事实密度

3.2.1 理论介绍

反事实分析中的一个关键特征是使用被估计出来的权重体系去估计在全部支撑集 (包括边界和边界邻域) 上的密度函数。为了构建一个反事实密度估计量,我们需要合理设定权重体系 (w1,w2,,wn)。此外,我们还需要估计这些权重的预先一致估计量。假设被观察到的数据为 (xi,ti,zi),i=1,2,,n,这些数据分为两组 (控制组和处理组), 其中 xi 是主要的结果变量,zi 是其他的协变量,ti 是指明 i 属于哪一组的二元变量。

结果变量 xi 全样本下的边际分布在无权重的条件下 (wi=1) 可以很容易被估计出来。控制组和处理组中的条件密度 f^0(x) 和 f^1(x) 分别可以通过 wi0=(1ti)/P[ti=0] 和 wi1=ti/P[ti=1] 估计出来。在随机控制实验中,这两种密度函数分别可以用来刻画控制组和处理组中结果变量的分布特征。

那么一个更具挑战性的问题出现了:如果处理组中的个体具有和控制组中的个体相同的协变量分布时,结果变量分布将会是怎样的呢?这时的处理组中结果变量的密度函数被称为反事实密度,记为 f10(x)。反事实密度对于理解 f1(x) 和 f0(x) 至关重要。

此外应该注意到的是,反事实分布还有一个更为有用的解释:假定结果变量是从潜在产出 xi=tixi(1)+(1ti)xi(0) 中产生的,并且假定条件于 ziti 是独立于 (xi(0),xi(1)) 的,那么 f10(x) 就是控制组的反事实分布:在 ti=0 的条件下 xi(1) 分布的密度函数。

f10(x) 可以很容易地通过上述的局部回归分布估计量以及下面的权重估计出来:

在实践当中,这种权重体系是未知的,因为倾向得分 P[ti=1zi] 是不可观测到的。因此研究者们往往选择 Probit 或者 Logit 模型去估计,我们的局部回归分布估计量中允许使用这类估计出的权重去替换理论上的权重,只要这种估计的权重可以收敛地足够快。

3.2.2 具体案例

在这一节中我们主要展示边际密度,条件密度以及反事实密度是如何通过使用局部分布回归的方法去估计的。具体地,我们采用 Abadie 等 (2002) 研究教育对收入影响数据集中的一个子数据集。这份数据中包括没有登记在 Job Training Partnership Act (JTPA) 的全部个体。

主要的结果变量是在 30 个月的期限内收入的加总,这些个体根据其教育程度被分为两组:t1=1 代表那些有高中学历的个体;t1=0 表示其他。这份数据中还包含性别、种族、年龄、婚姻状况、收到的 AFDC (Aid to Families with Dependent Children)、以及是否该个体在一年内工作至少 12 周。

样本量共有 5447,其中有 3927 个高中毕业生。该子样本的描述性统计我们放在下一个应用案例 (用到的是全样本) 中一起展示。众所周知的是教育对于劳动收入具有显著的影响,我们首先分别画出具有高中学历和不具有高中学历两个子样本的收入分布情况,即画出 f^1(x) 和 f^0(x)

// 调用 dta 格式数据
lxhuse jtpa.dta, clear
// 下载 csv 格式数据
lxhget jtpa.csv, replace 

* Figure 2
capture drop grid // set up grid points
qui gen grid = 1.97 + 0.03 * _n if _n <= 101
capture drop eduFull_*
capture drop edu1_*
capture drop edu0_*

// full sample
qui lpdensity logincome if treatment == 0, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(eduFull)
// with HS or GED
qui lpdensity logincome if treatment == 0 & hsorged > 0, bwselect(imse-dpi) grid(grid) ///
    cweights(hsorged) ciuniform	genvars(edu1)
// no HS or GED
qui lpdensity logincome if treatment == 0 & hsorged == 0, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(edu0)
// Figure 2(a) * NewVarName_f_p 表示用 p(#) 多项式得到的点估计
twoway (line eduFull_f_p eduFull_grid, col(red)) ///
    (line edu1_f_p edu1_grid, col(blue))		 ///
    (line edu0_f_p edu0_grid, col(green))		 ///
    (rarea eduFull_CI_l eduFull_CI_r eduFull_grid, col(red%20)   lcolor(red%0))	 ///
    (rarea edu1_CI_l edu1_CI_r edu1_grid, col(blue%20)  lcolor(blue%0))	         ///
    (rarea edu0_CI_l edu0_CI_r edu0_grid, col(green%20) lcolor(green%0)),        ///
    legend(cols(3) order(1 "Full" 2 "HS or GED" 3 "No HS or GED"))               ///
    xtitle("log(income)") ytitle("")

从上图可以看出高中毕业生和没有取得高中学历的毕业生之间的收入分布是非常不同的。f^1(x) 的均值和中位数都要高于 f^0(x)。但是正如前文所述,直接比较 f^1(x) 和 f^0(x) 是不能揭示拥有高中学历对收入的影响的,因为除高中学历以外的个体特征并没有被我们控制,所以下面我们使用局部回归分布估计量来估计反事实分布 f10(x)

// counterfactual weights estimation
capture drop male2nonwhite
capture drop male2wkless13
capture drop nonwhite2wkless13
capture drop afdc2male
capture drop cweights
qui gen male2nonwhite = male * nonwhite
qui gen male2wkless13 = male * wkless13
qui gen nonwhite2wkless13 = nonwhite * wkless13
qui gen afdc2male = afdc * male
qui reg hsorged male nonwhite wkless13 married male2nonwhite male2wkless13 nonwhite2wkless13 afdc2male /// 
    age2225 age2629 age3035 age3644 age4554 if treatment == 0
qui predict cweights
qui replace cweights = (1 - cweights) / cweights
qui replace cweights = cweights * hsorged if hsorged > 0

// estimate counterfactual density
capture drop educ_*
qui lpdensity logincome if treatment == 0 & hsorged > 0, bwselect(imse-dpi) grid(grid) ///
    cweights(cweights) ciuniform genvars(educ)

// Figure 2(b)
twoway (line educ_f_p educ_grid, col(red))		///
    (line edu1_f_p edu1_grid, col(blue))		///
    (line edu0_f_p edu0_grid, col(green))		///
    (rarea educ_CI_l educ_CI_r educ_grid, col(red%20)   lcolor(red%0))	     ///
    (rarea edu1_CI_l edu1_CI_r edu1_grid, col(blue%20)  lcolor(blue%0))      ///
    (rarea edu0_CI_l edu0_CI_r edu0_grid, col(green%20) lcolor(green%0)),    ///
    legend(cols(3) order(1 "Counterfactual" 2 "HS or GED" 3 "No HS or GED")) ///
    xtitle("log(income)") ytitle("")

从上图我们可以看到 f^10(x) 和 f^1(x) 的差别并不是很显著,尽管看上去 f^10(x) 的中位数和均值要小一些。另一方面,可以看到的是 f^0(x) 和 f^10(x) 的差距仍然非常显著。基于此我们可以得出结论:教育确实可以导致显著的人力资本积累进而提高劳动收入,所以教育水平确实是收入差距的一个极为重要的解释变量。

3.3 工具变量的识别和异质性

自选择和处理效应异质性是因果推断和社会经济学中需要重点考虑的方面。众所周知的是,包括平均处理效应或处理组的处理效应在内的古典处理效应,由于存在不完全服从个体从而导致即便在充分随机分组的条件下也不能被识别。基于此,我们首先应用局部回归分布估计量的密度函数来进行工具变量的识别检验,然后估计一下服从者潜在产出的密度函数。

我们采用和 3.2 节中相同的数据集。JIPA 是一项大型的对经济困难或存在就业障碍的个体提供再就业培训的公共福利计划。个体被随机提供 JIPA 培训,但是在被提供 JIPA 培训的样本当中只有 67% 的个体接受。因此 JTPA 项目是一个有效的工具变量用以研究就业培训对收入的影响。

我们首先通过局部分布估计量 f^(x) 进行工具变量识别检验,下面我们画出两种估计出来的密度函数。

lxhuse jtpa.dta, clear
qui gen grid = 1.97 + 0.03 * _n if _n <= 101
// full sample
capture drop full_*
qui lpdensity logincome, bwselect(imse-dpi) grid(grid) genvars(full) ciuniform

// JTPA offer
capture drop offer1_* 
capture drop offer0_* 
// JTPA offer
qui lpdensity logincome if instrument == 1, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(offer1)
// no JTPA offer
qui lpdensity logincome if instrument == 0, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(offer0)

// JTPA Enroll
capture drop treat1_* 
capture drop treat0_* 
// Enrolled in JTPA
qui lpdensity logincome if treatment == 1, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(treat1)
// not Enrolled in JTPA
qui lpdensity logincome if treatment == 0, bwselect(imse-dpi) grid(grid) ///
    ciuniform genvars(treat0)

// weights for compliers, Abadie (2003)
capture drop cweights
capture drop cweights0
capture drop cweights1
capture drop hsorged2male
capture drop hsorged2nonwhite
capture drop hsorged2wkless13
capture drop male2nonwhite
capture drop male2wkless13
capture drop nonwhite2wkless13
qui gen hsorged2male = hsorged * male
qui gen hsorged2nonwhite = hsorged * nonwhite
qui gen hsorged2wkless13 = hsorged * wkless13
qui gen male2nonwhite = male * nonwhite
qui gen male2wkless13 = male * wkless13
qui gen nonwhite2wkless13 = nonwhite * wkless13
qui reg instrument male hsorged nonwhite wkless13 married hsorged2male  hsorged2nonwhite ///
    hsorged2wkless13 male2nonwhite male2wkless13 nonwhite2wkless13
qui predict cweights
qui gen cweights0 = cweights
qui gen cweights1 = cweights
qui replace cweights = 1 - instrument * (1 - treatment) / cweights - ///
    (1 - instrument) * treatment / (1 - cweights)
qui replace cweights1 = treatment * (instrument - cweights1) / (cweights1 * (1 - cweights1))
qui replace cweights0 = (1 - treatment) * (1 - instrument - 1 +     ///
    cweights0) / (cweights0 * (1 - cweights0))

// density for compliers
capture drop complier_* 
qui lpdensity logincome, bwselect(imse-dpi) grid(grid) cweights(cweights) ///
    ciuniform genvars(complier)
capture drop complier0_* 
qui lpdensity logincome, bwselect(imse-dpi) grid(grid) cweights(cweights0) ///
    ciuniform genvars(complier0)
capture drop complier1_* 
qui lpdensity logincome, bwselect(imse-dpi) grid(grid) cweights(cweights1) ///
    ciuniform genvars(complier1)

// Figure 4(a)
twoway (line full_f_p full_grid, col(red))        ///
    (line complier_f_p complier_grid, col(blue))  ///
    (rarea full_CI_l full_CI_r full_grid, col(red%20) lcolor(red%0))                ///
    (rarea complier_CI_l complier_CI_r complier_grid, col(blue%20) lcolor(blue%0)), ///
    legend(cols(3) order(1 "Full" 2 "Compliers")) ///
    xtitle("log(income)") ytitle("") title("a") saving(a)
// Figure 4(b)
twoway (line offer0_f_p offer0_grid, col(red))	///
    (line offer1_f_p offer1_grid, col(blue))	///
    (rarea offer0_CI_l offer0_CI_r offer0_grid, col(red%20) lcolor(red%0))	  ///
    (rarea offer1_CI_l offer1_CI_r offer1_grid, col(blue%20) lcolor(blue%0)), ///
    legend(cols(3) order(1 "JTPA Offer = 0" 2 "1")) ///
    xtitle("log(income)") ytitle("") title("b") saving(b)
// Figure 4(c)
twoway (line treat0_f_p treat0_grid, col(red))	///
    (line treat1_f_p treat1_grid, col(blue))	///
    (rarea treat0_CI_l treat0_CI_r treat0_grid, col(red%20)   lcolor(red%0))   ///
    (rarea treat1_CI_l treat1_CI_r treat1_grid, col(blue%20)  lcolor(blue%0)), ///
    legend(cols(3) order(1 "JTPA Enrollment = 0" 2 "1")) ///
    xtitle("log(income)") ytitle("") title("c") saving(c)
// Figure 4(d)
twoway (line complier0_f_p complier0_grid, col(red))    ///
    (line complier1_f_p complier1_grid, col(blue))      ///
    (rarea complier0_CI_l complier0_CI_r complier0_grid, col(red%20) lcolor(red%0))    ///
    (rarea complier1_CI_l complier1_CI_r complier1_grid, col(blue%20) lcolor(blue%0)), ///
    legend(cols(3) order(1 "x(0) for compliers" 2 "x(1) for compliers")) ///
    xtitle("log(income)") ytitle("") title("d") saving(d)
graph combine a.gph b.gph c.gph d.gph

在图 (a) 中,我们绘制了全样本和只针对服从者的收入分布,两者看起来很相似,尽管后者有着更高的均值并且左尾更瘦。从图 (b) 中可以看到被提供 JIPA 和没有被提供 JIPA 的收入分布是显著不同的。在图 (c) 中,我们也绘制出是否接受 JIPA 对个体收入分布的影响,不出意外,两者的收入分布差距更大。

在图 (d) 中,我们构建出对于服从者而言的潜在收入分布,并估计了在这样的工具变量设定背景下被识别的分布处理效应。可以看出,JIPA 确实对收入有着正的显著的影响,而且它证实了自选择的存在:那些参与 JIPA 的个体平均而言会受益最多,其次是那些被认为在 “在无差异边缘” 的服从者。

4. 相关推文

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