温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者:黄欣怡 (中山大学) 邮箱:eugenehuangcheeks@163.com
「Source:Weights to address non-parallel trends in panel difference-in-differences models 」
目录
双重差分法 ( Difference in differences,简称为 DID 或 DD ) 是用于估计处理效应的最为常见的计量方法,广泛应用于评估政策实施、项目执行及其他外生事件产生的效应。DID 法通过衡量处理组与控制组在被处理前后发生的变化,将处理组的变动减去控制组的变动,得到处理效应。由于在估计过程中采用了两次差分,故名为双重差分法。
采用 DID 法的重要前提假设为平行趋势假设,即假定处理组在未受到处理情况下的变化趋势与控制组相同。该假设在具体实证研究中转化为检验处理组与控制组在接受处理前 ( pre-treatment period ) 是否具有相同的变化趋势。关于平行趋势假设的检验在连享会之前的推文「多期 DID:平行趋势检验图示」和「Stata: 多期倍分法 (DID) 详解及其图示」中均有涉猎。
实证研究中,当个体接受处理的可能性或接受处理产生的效应与该个体随时间变动的因素相关时,可能不再满足平行趋势假设。在违背平行趋势假设的情况下,不能直接采用 DID 估计处理效应。部分研究采用了增加控制变量的方法,在回归中进一步控制了交互固定效应 ( interactive fixed effects ) 或随时间变化的混合因素 ( time-varing confounding factors ) 。
与此同时,部分学者通过改变样本的权重,使处理组与新构造的控制组在接受处理前具有相似的变化趋势。这类方法包括合成控制法 ( synthetic control method ) 、倾向得分匹配法 ( propensity score matching ) 、逆概率加权法 ( inverse probability weighting ) 等。
Ahlfeldt (2018) 指出,这类计量方法只适用于估计单一处理效应的情形;在需要区分多类相互关联的处理效应的研究情境中,或面对连续型 DID 的情形时,这些方法不再适用。Ahlfeldt (2018) 由此提出 Weighted Parallel Trends DD ( 以下简称 WPT DD ) 的估计方法。下文主要对 Ahlfeldt (2018) 提出的网格搜索算法 ( Grid Search ) 和迭代法 ( Iterative Approach ) 进行介绍。
我们参照如下模型生成模拟数据:
其中,
在模型 (1) 中,
其中,
其中,
遵循这一思路,Ahlfeldt (2018) 首先生成总体
下文将利用样本接受处理前 ( pre-treatment period ) 的数据对权重
我们采用 Ahlfeldt (2018) 提供的 [数据和程序] 进行模拟。模拟数据由 1000 个体 ×10 年构成 (
在模拟数据生成过程中,我们设定
* 导入模拟数据
cd "E:\Data\GA__PTW_REP"
u "DATA/INPUT/MCDATA.dta", replace
* 生成不随时间变化的个体特征变量,即 b1, b2, b3
qui gen h1 = runiform(0,1)
qui gen h2 = runiform(0,1)
qui gen h3 = runiform(0,1)
* 生成模型(2)中对应的变量 Hi
qui gen H = (-1.5+h1+h2+h3)
* 生成两组 treatment variable (连续型 DID)
qui gen T1 = runiform(0,1)
qui gen T2 = runiform(0,1)
* 将 TREATMENT x H 加入至 TREND 变量中,从而生成模型 (2) 对应的个体变化趋势 OMEGA
qui replace TREND = TREND + (0.5 * T1 + 0.5 * T2) * H
* 生成抽样比率 Fi 与抽样权重 Si
* 生成 r1, r2, r3
gen RW = runiform(0,1)
* 生成抽样权重 Si
qui gen MODELW = (RW[1]*h1+RW[2]*h2+(RW[3])*h3)
* 对抽样权重作标准化处理
qui sum MODELW
qui replace MODELW = MODELW /r(mean)
* 生成抽样权重的倒数:抽样比率 Fi
qui gen IVW = 1/MODELW
* 依据抽样比率 Fi 从总体 j = 1,...,J 中抽取样本 i = 1,...,I
qui gsample `obs' [w=IVW]
* 将数据由宽型数据转换为长型数据
qui foreach num of numlist 1/10 {
preserve
keep ID LE YE Y`num' h* TREND T1 T2 MODELW IVW
gen PERIOD = `num'
ren Y`num' Y
gen panelID = _n
save "TEMP/re_temp`num'.dta", replace
restore
}
qui u "TEMP/re_temp1.dta", clear
qui erase "TEMP/re_temp1.dta"
qui foreach num of numlist 2/10 {
append using "TEMP/re_temp`num'.dta"
erase "TEMP/re_temp`num'.dta"
}
* 定义面板数据
qui xtset panelID PERIOD
* 将f(T) x OMEGA 加入至结果变量中
qui replace Y = Y + TREND*PERIOD
* 将处理效应 TREATMENT 加入至结果变量中 ,以t=5为接受处理的时点
qui replace Y = Y + T1 if PERIOD >=6
qui replace Y = Y + T2 if PERIOD >=6
* 生成 TREATMENT x POST DID 变量,其系数即我们考察的处理效应
qui gen T1_POST = T1 * (PERIOD >=6)
qui gen T2_POST = T2 * (PERIOD >=6)
* 在结果变量中加入随机扰动项
qui replace Y = Y + rnormal(0,0.1)
* 存储数据以进行下一步运算
save "TEMP/temp", replace
我们首先采用 OLS DD 估计模拟数据的处理效应。
** 估计 OLS DD
reg Y T1_POST T2_POST i.PERIOD, abs(panelID)
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
T1_POST | .6451472 .0786551 8.20 0.000 .4909652 .7993292
T2_POST | .3899474 .0813273 4.79 0.000 .2305274 .5493675
|
PERIOD |
2 | .2779982 .0528507 5.26 0.000 .1743988 .3815976
3 | .8783501 .0528507 16.62 0.000 .7747507 .9819495
4 | 1.648978 .0528507 31.20 0.000 1.545378 1.752577
5 | 1.991826 .0528507 37.69 0.000 1.888226 2.095425
6 | 2.832213 .0769469 36.81 0.000 2.68138 2.983046
7 | 2.538575 .0769469 32.99 0.000 2.387742 2.689409
8 | 3.607691 .0769469 46.89 0.000 3.456857 3.758524
9 | 3.472193 .0769469 45.12 0.000 3.32136 3.623027
10 | 4.313866 .0769469 56.06 0.000 4.163033 4.4647
|
_cons | 1.315558 .0373711 35.20 0.000 1.242303 1.388814
------------------------------------------------------------------------------
其次,我们已知生成样本数据所采用的抽样比率
** 采用抽样权重估计 WPT DD
reg Y T1_POST T2_POST i.PERIOD [w=MODELW], abs(panelID)
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
T1_POST | .8356594 .0766697 10.90 0.000 .6853693 .9859495
T2_POST | 1.0463 .0769499 13.60 0.000 .8954609 1.197139
|
PERIOD |
2 | .3825006 .0507956 7.53 0.000 .2829295 .4820716
3 | 1.083468 .0507956 21.33 0.000 .9838966 1.183039
4 | 1.959798 .0507956 38.58 0.000 1.860227 2.059369
5 | 2.404149 .0507956 47.33 0.000 2.304578 2.50372
6 | 2.931054 .0742346 39.48 0.000 2.785538 3.076571
7 | 2.739205 .0742346 36.90 0.000 2.593688 2.884722
8 | 3.910568 .0742346 52.68 0.000 3.765051 4.056084
9 | 3.879056 .0742346 52.25 0.000 3.733539 4.024573
10 | 4.824459 .0742346 64.99 0.000 4.678942 4.969976
|
_cons | 1.420531 .0359179 39.55 0.000 1.350124 1.490938
------------------------------------------------------------------------------
下文我们将分别通过网格搜索法 ( Grid Search ) 与迭代法 ( Iterative Approach ) 搜索
我们假定抽样权重
可采用的目标函数共分为 3 类,分别为:
在网格搜索中,我们选择叠加法进行演算,即最小化目标方程
* 导入数据并调整数据格式
qui u "TEMP/temp", clear
keep if PERIOD ==2
keep ID Y
ren Y Y2
save "TEMP/ptwtemp", replace
qui u "TEMP/temp", clear
keep if PERIOD ==1
merge m:m ID using "TEMP/ptwtemp.dta"
drop _m
ren Y Y1
erase "TEMP/ptwtemp.dta"
* 生成结果变量自第1期至第2期的变化值
gen DELTA = Y2- Y1
* 将个体特征变量标准化
foreach num of numlist 1/3 {
gen VAR`num'= h`num'
sum VAR`num',d
gen wa_VAR`num' = VAR`num'- r(min)
replace wa_VAR`num' = wa_VAR`num'/r(sd)
}
* 将TREATMENT VARIABLE标准化
sum T1
replace T1 = T1/ r(sd)
sum T2
replace T2 = T2 / r(sd)
* 生成网格搜索法中需要用到的暂时性变量
gen wa = .
gen OBJ = .
gen OBJmin = .
gen VAR1min = .
gen VAR2min = .
gen VAR3min = .
gen wa_min = .
* 确定搜索范围与搜索步长
local min = 0
local max = 1
local step = 0.1
* 运行网格搜索
* 开启对q1的循环
local om_VAR1 = `min'
while `om_VAR1' <= `max' {
* 开启对q2的循环
local om_VAR2 = `min'
while `om_VAR2' <= `max' {
* 开启对q3的循环
local om_VAR3 = `min'
while `om_VAR3' <= `max' {
* 生成权重Si
qui replace wa = `om_VAR1'*wa_VAR1+`om_VAR2'*wa_VAR2+`om_VAR3'*wa_VAR3
* 确保抽样权重Si>0
qui sum wa
if r(mean) > 0 {
* 对Si作标准化处理
qui replace wa = wa / r(mean)
* 依照模型(4)进行加权回归估计处理效应
qui reg DELTA T1 T2 [iw=wa]
* 计算目标方程
qui replace OBJ = _b[T1]^2+_b[T2]^2
* 若 qm 使目标方程较上一轮循环更小,则存储参数
qui replace VAR1min = `om_VAR1' if OBJ < OBJmin
qui replace VAR2min = `om_VAR2' if OBJ < OBJmin
qui replace VAR3min = `om_VAR3' if OBJ < OBJmin
* 若目标方程较上一轮循环更小,存储计算所得的 WPT 权重 Si
qui replace wa_min = wa if OBJ < OBJmin
qui replace OBJmin = OBJ if OBJ < OBJmin
display "VAR1 "`om_VAR1' " VAR2 " `om_VAR2' " 'VAR3 " `om_VAR3'
}
else {
}
* 结束对q3的循环
local om_VAR3 = `om_VAR3'+`step'
}
* 结束对q2的循环
local om_VAR2 = `om_VAR2'+`step'
}
* 结束对q1的循环
local om_VAR1 = `om_VAR1'+`step'
}
* 结束循环并存储样本权重
ren wa_min WA_MC
keep ID WA_MC
save "TEMP/WA_MC.dta", replace
* 将计算所得的样本权重与原模拟数据进行合并
qui u "TEMP/temp", clear
qui merge m:m ID using "TEMP/WA_MC.dta"
erase "TEMP/WA_MC.dta"
qui tab _m
qui drop _m
* 估计 WPT DD
reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)
* 去除权重变量
qui drop WA_MC
我们得到 WPT DD 的估计结果如下:
reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
T1_POST | .8385112 .077283 10.85 0.000 .6870188 .9900035
T2_POST | 1.090782 .0778193 14.02 0.000 .9382385 1.243326
|
PERIOD |
2 | .3796894 .0510905 7.43 0.000 .2795404 .4798385
3 | 1.079326 .0510905 21.13 0.000 .9791772 1.179475
4 | 1.952109 .0510905 38.21 0.000 1.85196 2.052258
5 | 2.392124 .0510905 46.82 0.000 2.291975 2.492273
6 | 2.893577 .0748875 38.64 0.000 2.74678 3.040373
7 | 2.69863 .0748875 36.04 0.000 2.551833 2.845426
8 | 3.865206 .0748875 51.61 0.000 3.718409 4.012002
9 | 3.832881 .0748875 51.18 0.000 3.686084 3.979677
10 | 4.77509 .0748875 63.76 0.000 4.628293 4.921886
|
_cons | 1.417661 .0361264 39.24 0.000 1.346845 1.488477
------------------------------------------------------------------------------
迭代法中,我们以
其中
在模型 (6) 中,估计得到的边际处理效应为
迭代的具体思路为:每一轮循环对权重
当模型 (5) 估计所得
迭代法中我们选取目标方程是
* 存储数据
save "TEMP/temp", replace
* 清除内存
set more off
estimates drop _all
program drop _all
* 设定迭代次数为100
local MAX = 100
* 设定目标方程的目标值为0.005
local TARGET = 0.005
* 导入数据并整理数据格式
u "TEMP/temp.dta", clear
keep if PERIOD ==2
keep ID Y
ren Y Y2
save "TEMP/ptwtemp.dta", replace
u "TEMP/temp.dta", clear
keep if PERIOD ==1
merge m:m ID using "TEMP/ptwtemp.dta"
drop _m
ren Y Y1
erase "TEMP/ptwtemp.dta"
* 生成结果变量自第1期至第2期的变化值
gen DELTA = Y2- Y1
* 对结果变量和处理变量作标准化处理
foreach var of varlist DELTA T1 T2 {
qui sum `var'
qui gen S_`var' = `var'/r(sd)
}
* 生成处理变量T1 T2与个体特征变量bm的交乘项
qui foreach var1 of varlist S_T1 S_T2 {
foreach var2 of varlist h1 h2 h3 {
gen `var1'_`var2' = `var1'*`var2'
}
}
* 生成迭代所需要的暂时性变量
gen W = 1
gen W1 = 1
gen W2 = 1
gen WOLD = 1
qui gen ME_T1 = .
qui gen ME_T2 = .
gen WBEST = 1
gen bBEST = 100
scalar TH = 0.00
gen R = .
* 定义迭代法所使用的程序 ************************************************
program ALGO
* 估计模型(5)
qui replace WOLD = W
qui reg S_DELTA S_T1 S_T2 [w=W]
* 存储估计所得的处理效应 (模型(5)中的c_s^n)
qui scalar bT1 = _b[S_T1]
qui scalar bT2 = _b[S_T2]
* 存储目标方程的函数值
qui scalar bsum = abs(bT1)+abs(bT2)
* 估计模型(6)
qui reg S_DELTA S_T1 S_T2 S_T1_h* S_T2_h* [w=W]
* 存储估计所得的处理效应 (模型(6)中的b_s^n,m)
qui replace ME_T1 = _b[S_T1]
qui replace ME_T2 = _b[S_T2]
qui foreach num of numlist 1/3 {
qui replace ME_T1 = ME_T1 + _b[S_T1_h`num']*h`num'
qui replace ME_T2 = ME_T2 + _b[S_T2_h`num']*h`num'
}
* 生成调整后的权重Si
* TREATMENT 1 (第一类处理效应)
* 若模型(5)中估计所得的处理效应为负
qui if bT1 < -1*TH {
display "neg"
* 增加处理效应为正的个体在样本中的权重
replace W1 = W + W*ME_T1 if ME_T1 > 0
}
else {
* 若模型(5)中估计所得的处理效应为正
if bT1 > TH {
display "pos"
* 增加处理效应为负的个体在样本中的权重
replace W1 = W - W*ME_T1 if ME_T1 < 0
}
else {
}
}
else {
}
* TREATMENT 2 (第二类处理效应)
* 若模型(5)中估计所得的处理效应为负
qui if bT2 < -1*TH {
display "neg"
* 增加处理效应为正的个体在样本中的权重
replace W2 = W + W*ME_T2 if ME_T2 > 0
}
else {
* 若模型(5)中估计所得的处理效应为正
if bT2 > TH {
display "pos"
* 增加处理效应为负的个体在样本中的权重
replace W2 = W - W*ME_T2 if ME_T2 < 0
}
else {
}
}
else {
}
* 根据两类处理效应的估计结果对权重进行调整
* 对于绝对值较大的处理效应给予更大幅度的调整
qui replace W =[exp(abs(bT1))* W1+exp(abs(bT2))*W2] / [exp(abs(bT1))+exp(abs(bT2))]
* 基于新生成的权重估计模型(5)
qui reg S_DELTA S_T1 S_T2 [w=W]
* 将处理效应存储为单值
qui scalar bT1new = _b[S_T1]
qui scalar bT2new = _b[S_T2]
* 计算新的目标方程函数值
qui scalar bsumnew = abs(_b[S_T1])+abs(_b[S_T2])
* 基于目标方程的函数值判定是否对权重进行更新
qui replace bBEST = bsumnew if bsumnew<bBEST
qui replace WBEST = W if bsumnew<bBEST
* 若上一轮循环得到的权重更优,则结合本次循环和上一轮循环得到的权重进行调整
if bsumnew > bBEST + 0.025 {
replace W = 0.8*WBEST+0.2*W
}
else {
}
end
* 程序ALGO定义结束 ****************************************************
* 迭代循环过程 *************************************************************
local it = 1
while `it' < `MAX' {
* 运行以上定义的程序 ALGO
qui ALGO
* 基于本轮循环得到的权重估计处理效应
qui reg S_DELTA S_T1 S_T2 [w=W]
* 存储估计所得处理效应
qui local T1CORR = _b[S_T1]
qui local T2CORR = _b[S_T2]
* 展示迭代结果
display "Iteration "`it' " T1 CORR " `T1CORR' " T2 CORR " `T2CORR'
* 检验是否达到目标方程函数值小于0.005的目标
if abs(`T1CORR') < `TARGET' & abs(`T2CORR') < `TARGET' {
local it = 1000
}
else {
local it = `it'+1
}
}
* 存储最终得到的pre-treatment period的处理效应
reg S_DELTA S_T1 S_T2 [w=W]
scalar T1CORR = _b[S_T1]
scalar T2CORR = _b[S_T2]
* 存储最终得到的最优权重
qui sum W
gen WA_MC = W / r(mean)
keep ID WA_MC
save "TEMP/WA_MC.dta", replace
* 将原模拟数据与权重数据合并
qui u "TEMP/temp", clear
qui merge m:m ID using "TEMP/WA_MC.dta"
erase "TEMP/WA_MC.dta"
qui tab _m
qui drop _m
* 估计 WPT DD
reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)
* 去除权重变量
qui drop WA_MC
erase "TEMP/temp.dta"
我们得到 WPT DD 的估计结果如下:
reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
T1_POST | .8385112 .077283 10.85 0.000 .6870188 .9900035
T2_POST | 1.090782 .0778193 14.02 0.000 .9382385 1.243326
|
PERIOD |
2 | .3796894 .0510905 7.43 0.000 .2795404 .4798385
3 | 1.079326 .0510905 21.13 0.000 .9791772 1.179475
4 | 1.952109 .0510905 38.21 0.000 1.85196 2.052258
5 | 2.392124 .0510905 46.82 0.000 2.291975 2.492273
6 | 2.893577 .0748875 38.64 0.000 2.74678 3.040373
7 | 2.69863 .0748875 36.04 0.000 2.551833 2.845426
8 | 3.865206 .0748875 51.61 0.000 3.718409 4.012002
9 | 3.832881 .0748875 51.18 0.000 3.686084 3.979677
10 | 4.77509 .0748875 63.76 0.000 4.628293 4.921886
|
_cons | 1.417661 .0361264 39.24 0.000 1.346845 1.488477
------------------------------------------------------------------------------
Ahlfeldt (2018) 进行了总计 1000 次蒙特卡洛模拟实验,图 1 呈现了蒙特卡洛模拟计算所得的处理效应 T1 和 T2 的分布,表 1 展示了处理效应分布的均值、中位数和标准差。
由图表可以看出,OLS DD 估计所得处理效应明显下偏,低于我们实际设定的处理效应
由此可以看出,WPT DD 有利于改善因违背平行趋势假设导致 DID 估计结果有偏的问题。
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看,所有课程可以随时购买观看。
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
⭕ Stata数据清洗 | 游万海 | 直播, 2 小时,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD