【应用计量系列151】Sant' Anna的DID规定动作（4-4）：处理配置机制——基于可观测协变量的处理

当我们第一次听说“标准”双向固定效应TWF声明在交叠处理场景中存在偏误时，经常会听到或读到作者说，由于他们没有交叠处理，所以他们将使用TWFE，且包含各种协变量。今天，我将说明当你的潜在结果具有特定协变量平行趋势或其他形式的基于协变量的异质性趋势时，双向固定效应模型也存在偏误。前面的内容，我们呈现过另一种处理配置机制——基于潜在结果变量的处理：

第四步：处理配置机制---理论篇

第四步：处理配置机制---应用篇

第四步：处理配置机制---阿什费尔特沉降（Ashenfelter's Dip）

基于可观测变量的处理和条件平行趋势

我会用模拟数据来说明我想要展示的问题。我将在此数据生成过程中介绍两件事。首先，我让个体根据其可观察到的协变量来选择处理。在我的案例中，它将基于年龄和高中 GPA，这是一项关于读大学---收入回报的研究。从这个意义上讲，我想继续探索Ghanem 等人（2024 年）概述的各种选择处理机制，以描述什么机制不会影响平行趋势，什么机制会打破平行趋势。

但是，该数据生成过程的第二部分将是平行趋势仅基于协变量的条件下成立。这意味着什么？考虑平行趋势成立，但处理组和对照组的个体构成不同。例如，处理组可能由城市组成，但对照组由郊区或农村地区组成。也许你想知道在城市建造大型体育场对税收的影响。平行趋势意味着，建造体育场的城市与没有体育场的农村和郊区的发展类似。

如果你怀疑处理组和对照组是否会有相似的反事实趋势，因为处理组和对照组的协变量分布非常不同，那么你就是在怀疑你的平行趋势假设。但是，如果你相信一旦你以人口规模为条件，平行趋势对于处理组和对照组而言是相似的，那么你就在说平行趋势确实成立——只是它在某种程度上适用于人口密集的地区，而对于人口稀少的地区则以另一种方式适用。

这就是“条件平行趋势”的含义，我第一次看到它的方程是在Heckman、Ichimura 和 Todd (1997) 的论文中（第 613 页上方程 A-4' 和 A-5 之间的未编号方程），这也是Abadie (2005)以及最近的Sant'Anna 和 Zhou (2020)中的假设。

然而，当我们估计带有协变量的回归模型时，我们通常不会想到条件平行趋势。我们通常考虑的是遗漏变量偏差或无混杂性。但是，当你估计DID模型时，你明确地诉诸于假设没有混杂因素的平行趋势。换句话说，你说的是，平均潜在结果 Y(0) 对于两个组来说会类似演变，因此，如果它们对于两个组来说类似演变，那么我们对混杂因素该说些什么呢？它们要么不存在，要么不重要，要么对于处理和对照组来说都相同。

换句话说，条件平行趋势并不是意味着无混淆性。它只是说不同的个体根据其协变量具有不同的趋势，但具有相似协变量的个体具有相似的反事实趋势，仅此而已。那么，是什么产生了我的条件平行趋势呢？

首先，我考察基于可观察变量来选择进入处理。GPA 较高且年龄较大的个体更有可能属于处理组，而不是对照组。我使用以下方法选择进入处理的机制。

* Generate Covariates (Baseline values)
  gen age = rnormal(35, 10)
  gen gpa = rnormal(2.0, 0.5)

  * Center Covariates (Baseline)
  sum age, meanonly
  qui replace age = age - r(mean)
  sum gpa, meanonly
  qui replace gpa = gpa - r(mean)

  * Generate Polynomial and Interaction Terms (Baseline)
  gen age_sq = age^2
  gen gpa_sq = gpa^2
  gen interaction = age * gpa
  
  * Treatment probability increases with age and decrease with gpa
  gen propensity = 0.3 + 0.3 * (age > 0) + 0.2 * (gpa > 0)
  gen treat = runiform() < propensity

在上述模拟数据中，处理的倾向得分很简单——如果年龄或 GPA 高于平均水平，那么接受处理的概率更高，我们使用另一个随机变量（“runiform()）”来比较。这会将人们纳入和排除在处理组之外，如果我们愿意，我们可以估计倾向得分并查看两组的可比性。虽然处理组的平均年龄比对照组大 4 岁左右，并且 GPA 更高，但情况还不错，有重叠：

这很好，因为对于在下面使用的双重稳健方法，我需要在这些协变量中至少有一些处理和控制个体。这是因为条件平行趋势基本上是取两个看起来不同的组，并或多或少地“修复它们”，以便平行趋势可以通过协变量发挥作用，即使平行趋势平均不适用于两组比较。以下是条件平行趋势：

其中b是处理前或者基期，t是处理后时期，结果是潜在结果，上标0是未处理的潜在结果。这意味着，如果城市没有建造体育场，那么，对于人口规模相当的地区，城市收税变化与农村/郊区类似。

我们用下列代码，看看实践的含义：

* Generate fixed effect with control group making 10,000 more at baseline
  qui gen unit_fe = 40000 + 10000 * (treat == 0) 
  
  * Generate Potential Outcomes with Baseline and Year Difference
  gen e = rnormal(0, 1500)
  qui gen     y0 = unit_fe        + 100 * age + 1000 * gpa + e if year == 1987
  qui replace y0 = unit_fe + 1000 + 150 * age + 1500 * gpa + e if year == 1988
  qui replace y0 = unit_fe + 2000 + 200 * age + 2000 * gpa + e if year == 1989
  qui replace y0 = unit_fe + 3000 + 250 * age + 2500 * gpa + e if year == 1990
  qui replace y0 = unit_fe + 4000 + 300 * age + 3000 * gpa + e if year == 1991
  qui replace y0 = unit_fe + 5000 + 350 * age + 3500 * gpa + e if year == 1992

  * Covariate-based treatment effect heterogeneity
  gen y1 = y0
  replace y1 = y0 + 1000 + 250 * age + 1000 * gpa if year >= 1991

我做了两件事。

首先，请注意 Y(0)（DID中的关键潜在结果）每年增长 1000，但也会根据个人的基期年龄和 GPA 增长不同的规模。

然后我做的第二件事是引入基于协变量的异质处理效应。解释比看起来要简单：如果一个人从未获得大学学位，那么他们的收入会随着年龄的增长而增加。但如果他们确实获得了大学学位，那么他们的收入会随着年龄的增长而增加得更多。这是一个声明，即年龄对收入的影响取决于你是否接受过处理。

使用这种数据生成过程，1991 年和 1992 年的 ATT 只是接受处理的个人及其个人处理效应的加总。这取决于抽签，因为它会随着年龄和 GPA 而变化，而年龄和 GPA 有点随机，但最终计算结果为 1991 年约为 1620 美元，1992 年约为 1620 美元（即没有动态处理效应）。

使用OLS估计

首先展示来自 OLS 的估计值，其中我控制了年龄、gpa、年龄平方、gpa 平方和年龄-gpa 交互项。

代码如下：

reg earnings treat##ib1990.year age gpa age_sq gpa_sq interaction, robust

我运行了 1,000 次，每次都估算该方程并保存处理变量的系数。然后，绘制了这 1,000 次蒙特卡罗试验的平均值和 95% 的上限和下限。当我们绘制事件研究图并将其与“真实效应”叠加时，我们发现了以下结果：

上述结果中，有几个值得注意的地方：

第一，1991年开始才有处理效应；

第二，1992年和1991年的处理效应相同，即没有动态处理效应；

第三，TWFE估计量是有偏的。处理前趋势存在偏误，处理后处理效应也存在偏误。如果使我们在跑TWFE事件研究，看到处理前系数显著，可能我们要疯掉了。今天还有人问我，如果事件研究的处理前趋势检验通不过，是不是DID就不能做了。在这里，我们的模拟数据处理前效应为0，所以TWFE的处理前系数时有偏误的，这并不是因为我们遗漏了协变量，我们把所有的协变量都控制住了（正如绝大多数人在实践中作的那样），而是因为TWFE要求处理效应在协变量上是不变的，并且要求协变量对潜在结果Y0的影响也是恒定的。但是，我们上面的模拟数据恰恰是潜在结果Y0随着年龄的增长而提高，而获得大学学位的收入回报随着年龄、GPA的变化而变化。

因此，TWFE 的这种偏误并非来自异质性处理效应，也不是来自交叠处理。相反，它来自不同的协变量趋势，以及与协变量相关的异质性处理效应。

我在我的DID课程中，讲过很多关于协变量的问题。

还有一个解决方案就是双重稳健估计量（Sant'Anna 和 Zhao，2021）。我们使用STATA的csdid命令。

双重稳健估计量具有以下表达式：

因为只有一组，所以就变成了 1991 年和 1992 年的 ATT，以及三个处理前时期，即 1987 年、1988 年和 1989 年。括号内部项是权重，第二个外部项是从基期到 1991 年或 1992 年的收入变化 (Delta Y) 减去如果经过处理D=1 的预测反事实。双重稳健结合了 Abadie 的倾向得分加权法和 Heckman、Ichimura 和 Todd 早期的结果回归方法。只需要其中一个模型正确，该方法就是无偏的。并且它对处理效应异质性和影响 Y(0)（没有上过大学的收入）的协变量在不同年份（我们也有）上具有不同的稳健性。估计这个的命令使用 - csdid -。

csdid earnings age gpa age_sq gpa_sq interaction, ivar(id) time(year) gvar(treat_date)

我也运行了 1,000 次，保存了 1987、1988、1989、1991 和 1992 年的估计值，并将它们绘制为以下事件研究：

结果显示，CS 对 1991 年和 1992 年的处理前系数以及处理后 ATT 参数进行了无偏估计。现在让我们将它们放在同一张图上：

我们既然已经在TWFE回归中包含了所有协变量，那为什么OLS仍然有偏呢？

首先，正如我所说，协变量存在处理效应异质性。

然后，年龄和GPA影响收入，而且是存在不同年份有不同的影响。

这个时候，重新思考协变量的异质性趋势就很容易了。许多经济学家会绘制不同年龄段工人的平均收入。我们在这里要说的是，收入会随着年龄而变化。我所做的是，在没有人获得大学学位的世界里，年龄对收入有自己的影响。这被称为“特定协变量趋势”，但需要注意的是，这是 Y(0) 的趋势，如果你还记得的话，平行趋势就是 Y(0) 的变化对于两个组来说都是相同的。这是我制作的一张图，显示了不同年龄段的平均 Y(0)。

正如我所说，我运行的 OLS 模型要求处理效应恒定，并且要求年龄和收入不会对 Y(0) 产生趋势影响。如果这些数据中两者都不成立，那么 OLS 就是有偏的。

在这里我们再次看到“外生性”对 DGP 施加了强有力的假设。它对 Y(0) 和恒定处理效应施加了函数形式假设。我们本可以在 DID中使用回归调整来解决这个问题，但我只是运行了简单的 OLS，因为绝大部分应用研究者要控制协变量时都会这样做。但是上面的模拟数据告诉我们，这还不够。

我们学到了什么？

第一件事，基于可观测变量来选择进入处理并不会影响平行趋势。

第二件事，即使我们在TWFE中控制了所有协变量，OLS估计量仍然可能有偏。我们可以用两种方式来尝试稳健性检验：

① 控制处理前协变量与（不同）时间趋势交乘项；
② 使用双重稳健估计量

    
  
    
  
* name: step4d.do
* Selection on observables and conditional parallel trends

clear all
capture log close
set seed 12345

********************************************************************************
* Define dgp
********************************************************************************
cap program drop dgp
program define dgp

  * 1,000 workers (25 per state), 40 states, 4 groups (250 per group), 6 years
  * First create the states
  quietly set obs 40
  gen state = _n

  * Generate 1000 workers. These are in each state. So 25 per state.
  quietly expand 25
  bysort state: gen worker=runiform(0,5)
  label variable worker "Unique worker fixed effect per state"
  quietly egen id = group(state worker)

  * Generate Covariates (Baseline values)
  gen age = rnormal(35, 10)
  gen gpa = rnormal(2.0, 0.5)

  * Center Covariates (Baseline)
  sum age, meanonly
  qui replace age = age - r(mean)
  sum gpa, meanonly
  qui replace gpa = gpa - r(mean)

  * Generate Polynomial and Interaction Terms (Baseline)
  gen age_sq = age^2
  gen gpa_sq = gpa^2
  gen interaction = age * gpa
  
  * Treatment probability increases with age and decrease with gpa
  gen propensity = 0.3 + 0.3 * (age > 0) + 0.2 * (gpa > 0)
  gen treat = runiform() < propensity

  * Generate the years
  quietly expand 6
  sort state
  bysort state worker: gen year = _n
  * years 1987 -- 1992
  replace year = 1986 + year
  
  * Post-treatment
  gen post = 0  
  qui replace post = 1 if year >= 1991

  * Generate treatment dates
  gen treat_date = 0
  replace treat_date = 1991 if treat==1

  * Generate fixed effect with control group making 10,000 more at baseline
  qui gen unit_fe = 40000 + 10000 * (treat == 0) 
  
  * Generate Potential Outcomes with Baseline and Year Difference
  gen e = rnormal(0, 1500)
  qui gen     y0 = unit_fe        + 100 * age + 1000 * gpa + e if year == 1987
  qui replace y0 = unit_fe + 1000 + 150 * age + 1500 * gpa + e if year == 1988
  qui replace y0 = unit_fe + 2000 + 200 * age + 2000 * gpa + e if year == 1989
  qui replace y0 = unit_fe + 3000 + 250 * age + 2500 * gpa + e if year == 1990
  qui replace y0 = unit_fe + 4000 + 300 * age + 3000 * gpa + e if year == 1991
  qui replace y0 = unit_fe + 5000 + 350 * age + 3500 * gpa + e if year == 1992

  * Covariate-based treatment effect heterogeneity
  gen y1 = y0
  replace y1 = y0 + 1000 + 250 * age + 1000 * gpa if year >= 1991

  * Treatment effect
  gen delta = y1 - y0
  label var delta "Treatment effect for unit i (unobservable in the real world)"
  
  * Generate observed outcome based on treatment assignment
  gen earnings = y0
  qui replace earnings = y1 if treat == 1 & year >= 1991
end

********************************************************************************
* Generate a sample
********************************************************************************
clear
quietly dgp

* Regression breaks
reg earnings treat##ib1990.year age gpa age_sq gpa_sq interaction, robust

* CSDID works
csdid earnings age gpa age_sq gpa_sq interaction, ivar(id) time(year) gvar(treat_date)
csdid_plot, group(1991)


********************************************************************************
* Monte-carlo simulation
********************************************************************************
cap program drop simulation
program define simulation, rclass
    clear 
	  quietly dgp
    
    // True ATT
    gen true_att = y1 - y0
    qui sum true_att if treat == 1 & year == 1991
    return scalar att_1991 = r(mean)
    qui sum true_att if treat == 1 & year == 1992
    return scalar att_1992 = r(mean)

    // CSDID
    qui csdid earnings age gpa age_sq gpa_sq interaction, ivar(id) time(year) gvar(treat_date)
    matrix b = e(b)
    return scalar cs_pre1987 = b[1,1]
    return scalar cs_pre1988 = b[1,2]
    return scalar cs_pre1989 = b[1,3]
    return scalar cs_post1991 = b[1,4]
    return scalar cs_post1992 = b[1,5]
    return scalar cs_att = (b[1,4] + b[1,5]) / 2
    
    // OLS
    qui reg earnings treat##ib1990.year age gpa age_sq gpa_sq interaction, robust
    return scalar ols_pre1987 = _b[1.treat#1987.year]
    return scalar ols_pre1988 = _b[1.treat#1988.year]
    return scalar ols_pre1989 = _b[1.treat#1989.year]
    return scalar ols_post1991 = _b[1.treat#1991.year]
    return scalar ols_post1992 = _b[1.treat#1992.year]
    
    // OLS overall ATT
    qui reg earnings post##treat age gpa age_sq gpa_sq interaction, robust
    return scalar ols_att = _b[1.post#1.treat]  
end

simulate att_1991 = r(att_1991) ///
         att_1992 = r(att_1992) ///
         cs_pre1987 = r(cs_pre1987) ///
         cs_pre1988 = r(cs_pre1988) ///
         cs_pre1989 = r(cs_pre1989) ///
         cs_post1991 = r(cs_post1991) ///
         cs_post1992 = r(cs_post1992) ///
         cs_att = r(cs_att) ///
         ols_pre1987 = r(ols_pre1987) ///
         ols_pre1988 = r(ols_pre1988) ///
         ols_pre1989 = r(ols_pre1989) ///
         ols_post1991 = r(ols_post1991) ///
         ols_post1992 = r(ols_post1992) ///
         ols_att = r(ols_att), ///
         reps(100) seed(12345): simulation


// Summarize results
sum

// Store results
save ./step4d.dta, replace





********************************************************************************
* Plot results
********************************************************************************

use ./step4d.dta, clear

* Calculate means and standard deviations for OLS variables
summarize ols_pre1987
local ols_pre1987_mean = r(mean)
local ols_pre1987_sd = r(sd)

summarize ols_pre1988
local ols_pre1988_mean = r(mean)
local ols_pre1988_sd = r(sd)

summarize ols_pre1989
local ols_pre1989_mean = r(mean)
local ols_pre1989_sd = r(sd)

summarize ols_post1991
local ols_post1991_mean = r(mean)
local ols_post1991_sd = r(sd)

summarize ols_post1992
local ols_post1992_mean = r(mean)
local ols_post1992_sd = r(sd)

* Calculate means and standard deviations for CSDID variables
summarize cs_pre1987
local cs_pre1987_mean = r(mean)
local cs_pre1987_sd = r(sd)

summarize cs_pre1988
local cs_pre1988_mean = r(mean)
local cs_pre1988_sd = r(sd)

summarize cs_pre1989
local cs_pre1989_mean = r(mean)
local cs_pre1989_sd = r(sd)

summarize cs_post1991
local cs_post1991_mean = r(mean)
local cs_post1991_sd = r(sd)

summarize cs_post1992
local cs_post1992_mean = r(mean)
local cs_post1992_sd = r(sd)

summarize att_1992
local true_att_1991 = r(mean)
summarize att_1991
local true_att_1992 = r(mean)


* Create a new dataset for plotting
clear
set obs 5

* Define the years
gen year = 1987 + _n - 1
replace year = 1991 if _n == 4
replace year = 1992 if _n == 5

* True ATT values
gen truth = 0
replace truth = `true_att_1991' if year == 1991
replace truth = `true_att_1992' if year == 1992

* OLS means and confidence intervals
gen ols_mean = .
gen ols_ci_lower = .
gen ols_ci_upper = .
replace ols_mean = `ols_pre1987_mean' in 1
replace ols_mean = `ols_pre1988_mean' in 2
replace ols_mean = `ols_pre1989_mean' in 3
replace ols_mean = `ols_post1991_mean' in 4
replace ols_mean = `ols_post1992_mean' in 5

replace ols_ci_lower = ols_mean - 1.96 * `ols_pre1987_sd' in 1
replace ols_ci_lower = ols_mean - 1.96 * `ols_pre1988_sd' in 2
replace ols_ci_lower = ols_mean - 1.96 * `ols_pre1989_sd' in 3
replace ols_ci_lower = ols_mean - 1.96 * `ols_post1991_sd' in 4
replace ols_ci_lower = ols_mean - 1.96 * `ols_post1992_sd' in 5

replace ols_ci_upper = ols_mean + 1.96 * `ols_pre1987_sd' in 1
replace ols_ci_upper = ols_mean + 1.96 * `ols_pre1988_sd' in 2
replace ols_ci_upper = ols_mean + 1.96 * `ols_pre1989_sd' in 3
replace ols_ci_upper = ols_mean + 1.96 * `ols_post1991_sd' in 4
replace ols_ci_upper = ols_mean + 1.96 * `ols_post1992_sd' in 5

* CSDID means and confidence intervals
gen csdid_mean = .
gen csdid_ci_lower = .
gen csdid_ci_upper = .
replace csdid_mean = `cs_pre1987_mean' in 1
replace csdid_mean = `cs_pre1988_mean' in 2
replace csdid_mean = `cs_pre1989_mean' in 3
replace csdid_mean = `cs_post1991_mean' in 4
replace csdid_mean = `cs_post1992_mean' in 5

replace csdid_ci_lower = csdid_mean - 1.96 * `cs_pre1987_sd' in 1
replace csdid_ci_lower = csdid_mean - 1.96 * `cs_pre1988_sd' in 2
replace csdid_ci_lower = csdid_mean - 1.96 * `cs_pre1989_sd' in 3
replace csdid_ci_lower = csdid_mean - 1.96 * `cs_post1991_sd' in 4
replace csdid_ci_lower = csdid_mean - 1.96 * `cs_post1992_sd' in 5

replace csdid_ci_upper = csdid_mean + 1.96 * `cs_pre1987_sd' in 1
replace csdid_ci_upper = csdid_mean + 1.96 * `cs_pre1988_sd' in 2
replace csdid_ci_upper = csdid_mean + 1.96 * `cs_pre1989_sd' in 3
replace csdid_ci_upper = csdid_mean + 1.96 * `cs_post1991_sd' in 4
replace csdid_ci_upper = csdid_mean + 1.96 * `cs_post1992_sd' in 5

twoway (scatter truth year, mcolor(maroon) msize(6-pt) msymbol(lgx) mlabcolor() mfcolor(cranberry) mlwidth(medthick)) ///
       (scatter ols_mean year, mcolor(navy) msize(6-pt)) ///
       (line ols_mean year, lcolor(blue) lwidth(medthick)) ///
       (rcap ols_ci_lower ols_ci_upper year, lcolor(blue) lpattern(dash)), ///
       title("OLS vs Truth Event Study") ///
       subtitle("Additive controls specification") ///
       note("DGP uses conditional parallel trends and all controls included. 1000 Monte Carlo simulations.") ///
       legend(order(1 "Truth" 2 "OLS")) ///
	   xline(1990.5, lpattern(dash) lcolor(gray))


graph export "/Users/scunning/Library/CloudStorage/Dropbox-MixtapeConsulting/scott cunningham/0.1 Mixtape Consulting/Substack/Code/step4d_ols.png", as(png) name("Graph")	   replace
	   
twoway (scatter truth year, mcolor(maroon) msize(6-pt) msymbol(lgx) mlabcolor() mfcolor(cranberry) mlwidth(medthick)) ///
       (scatter csdid_mean year, mcolor(saddle) msize(6-pt)) ///
       (line csdid_mean year, lcolor(brown) lwidth(medthick)) ///
       (rcap csdid_ci_lower csdid_ci_upper year, lcolor(brown) lpattern(dash)), ///
       title("CS vs Truth Event Study") ///
       subtitle("Doubly Robust specification") ///
       note("DGP uses conditional parallel trends and all controls included. Estimated with csdid." "1000 Monte Carlo simulations.") ///
       legend(order(1 "Truth" 2 "CS")) ///
	   xline(1990.5, lpattern(dash) lcolor(gray))
	   
graph export "/Users/scunning/Library/CloudStorage/Dropbox-MixtapeConsulting/scott cunningham/0.1 Mixtape Consulting/Substack/Code/step4d_cs.png", as(png) name("Graph")	   replace

* Shift years slightly to avoid overlap
gen year_ols = year - 0.1
gen year_csdid = year

* Plotting
twoway (scatter truth year, mcolor(maroon) msize(6-pt) msymbol(lgx) mlabcolor() mfcolor(cranberry) mlwidth(medthick)) ///
       (scatter ols_mean year_ols, mcolor(navy) msize(6-pt)) ///
       (line ols_mean year_ols, lcolor(blue) lwidth(medthick)) ///
       (rcap ols_ci_lower ols_ci_upper year_ols, lcolor(blue)) ///
       (scatter csdid_mean year_csdid, mcolor(saddle) msize(6-pt)) ///
       (line csdid_mean year_csdid, lcolor(brown) lwidth(medthick) lpattern(dash)) ///
       (rcap csdid_ci_lower csdid_ci_upper year_csdid, lcolor(brown) lpattern(dash)), ///
       title("Event Study: OLS, CSDID, and Truth") ///
       note("DGP uses conditional parallel trends. OLS includes additive controls; CS uses double robust." "No differential timing. 1000 Monte Carlo simulations.") ///
       legend(order(1 "Truth" 3 "OLS" 6 "CS") ///
              label(1 "Truth" ) ///
              label(2 "OLS" ) ///
              label(3 "CS" )) ///
       xline(1990.5, lpattern(dash) lcolor(gray))

* Export the graph
graph export "/Users/scunning/Library/CloudStorage/Dropbox-MixtapeConsulting/scott cunningham/0.1 Mixtape Consulting/Substack/Code/step4d_combined.png", as(png) name("Graph") replace

宏观研学会

关注宏观经济与宏观经济学研究前沿问题，推广、普及DSGE及其他宏观经济研究方法。