-
Notifications
You must be signed in to change notification settings - Fork 34
/
Copy pathTWFE_vs_SDID.do
296 lines (233 loc) · 10.1 KB
/
TWFE_vs_SDID.do
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
* This do file finishes four tasks:
* (1) It shows some Stata commands to obtain difference-in-differences estimates by the two-way fixed effects (TWFE) regression.
* (2) Following Arkhangelsky et al. (2021) and Clarke et al. (2023), it runs synthetic DID. The data are from Orzechowski & Walker (2005) and data from Bhalotra et al. (2022).
* (3) It shows how to run dynamic DID, using data from Serrato & Zidar (2018).
* (4) It briefly instroduces "xthdidregress", a new command introduced in Stata 18.
* Author: Ian He
* Date: June 29, 2023
********************************************************************************
clear all
global localdir "D:\research\DID Example"
global dtadir "$localdir\Data"
global figdir "$localdir\Figure"
********************************************************************************
**# Block Treatment: Orzechowski & Walker (2005)
use "$dtadir\OW05_prop99.dta", clear
encode state, gen(state_code)
xtset state_code year
table year treated
* Balanced panel: 39 states, years from 1970 to 2000.
* Treated group: California (1 unit).
* Control group: 38 other states.
* TWFE DID
xtdidregress (packspercapita) (treated), group(state_code) time(year) vce(cluster state_code)
estat trendplots // visualization
estat ptrends // parallel-trends test
estat granger // anticipation test
estat grangerplot // time-specific treatment effects
eststo twfe1: qui xtreg packspercapita treated i.year, fe cluster(state_code)
eststo twfe2: qui areg packspercapita treated i.year, absorb(state_code) cluster(state_code)
eststo twfe3: qui reghdfe packspercapita treated, absorb(state_code year) cluster(state_code)
estout twfe*, keep(treated) ///
coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///
starlevels(* .1 ** .05 *** .01) legend ///
stats(N r2_a, nostar labels("Observations" "Adjusted R2") fmt("%9.0fc"3))
* Sythetic DID
eststo syn_did: sdid packspercapita state year treated, vce(placebo) seed(1) ///
graph g1on ///
g1_opt(xtitle("") ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
) ///
g2_opt(ylabel(0(25)150) ytitle("Packs per capita") ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
) ///
graph_export("$figdir\prop99_did_", .pdf)
ereturn list
matrix list e(omega) // unit-specific weights
matrix list e(lambda) // time-specific weights
estout twfe* syn_did, keep(treated) ///
coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///
starlevels(* .1 ** .05 *** .01) legend ///
stats(N r2_a, nostar labels("Observations" "Adjusted R2") fmt("%9.0fc"3))
* Compare SDID, DID, and SC
foreach m in sdid did sc {
if "`m'"=="did" {
local g_opt = "msize(small)"
}
else {
local g_opt = ""
}
sdid packspercapita state year treated, ///
method(`m') vce(noinference) graph g1on `g_opt' ///
g1_opt(xtitle("") ytitle("") ///
xlabel(, labsize(tiny)) ///
ylabel(-100(25)50, labsize(small)) ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
) ///
g2_opt(title("`m'") xtitle("") ytitle("") ///
xlabel(, labsize(small)) ///
ylabel(0(25)150, labsize(small)) ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
)
graph save g1_1989 "$figdir\\`m'_1.gph", replace
graph save g2_1989 "$figdir\\`m'_2.gph", replace
}
graph combine "$figdir\sdid_2.gph" "$figdir\did_2.gph" "$figdir\sc_2.gph" "$figdir\sdid_1.gph" "$figdir\did_1.gph" "$figdir\sc_1.gph", ///
cols(3) xsize(3.5) ysize(2) ///
graphregion(fcolor(white) lcolor(white))
graph export "$figdir\compare_sdid_did_sc.pdf", replace
* Visualize the year-by-year effects
* We need to run the sdid as many times as post treatment times, and each time we use all pre-treatment years (1970-1988) and only one post-treatment year.
matrix te = J(12, 3, .) // generate a matrix to collect the estimates
local j = 1
forval t = 1989(1)2000 {
sdid packspercapita state year treated if year<=1988 | year==`t', vce(placebo) seed(1)
* store point and interval estimates
local b = e(tau)[1,1]
local se = e(se)
local lb = `b' + invnormal(0.025)*`se'
local ub = `b' + invnormal(0.975)*`se'
matrix te[`j',1] = `b'
matrix te[`j',2] = `lb'
matrix te[`j',3] = `ub'
local ++j
}
matlist te
matrix coln te = estimate lower upper
clear
svmat te, n(col) // covert the matrix to variables
gen year = _n + 1988
twoway (line estimate year, lc(black)) ///
(rcap lower upper year, lc(black)) ///
(scatter estimate year, mc(black) msize(medium)), ///
yline(0, lc(gs12) lp(dash)) ///
xlabel(1989(1)2000) ///
title("Year-by-Year Treatment Effects") xtitle("Year") ytitle("") ///
xlab(, nogrid) ylabel(-50(10)20) legend(off)
graph export "$figdir\sdid_att_over_time.pdf", replace
********************************************************************************
**# Staggered Treatment: Bhalotra et al. (2022)
use "$dtadir\BCGV22_gender_quota.dta", clear
drop if lngdp==.
table year quota
* Balanced panel: 115 countries, years from 1990 to 2015.
* Treated group: 9 countries.
* Control group: 106 other states.
eststo stagg1: sdid womparl country year quota, vce(bootstrap) seed(3) ///
graph g1on ///
g1_opt(xtitle("") ytitle("") xlabel(, labsize(tiny)) ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
) ///
g2_opt(ytitle("") ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
)
ereturn list
matrix list e(tau) // adoption-period specific estimate
matrix list e(lambda)
matrix list e(omega)
matrix list e(adoption) // different adoption years
* Covariate adjustments
eststo stagg2: sdid womparl country year quota, covariates(lngdp, optimized) vce(bootstrap) seed(3)
eststo stagg3: sdid womparl country year quota, covariates(lngdp, projected) vce(bootstrap) seed(3)
matrix list e(beta) // coefficient estimates of covariates
* Obtain the std. err. and p-values of coefficients of covariates
bysort country: egen treat_num = total(quota)
reghdfe womparl lngdp if treat_num==0, abs(country year) cluster(country)
* Compare three SDID results
estout stagg*, ///
coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///
starlevels(* .1 ** .05 *** .01) legend ///
stats(N, nostar labels("Observations") fmt("%9.0fc"))
********************************************************************************
**# Staggered Treatment: Serrato & Zidar (2018)
use "$dtadir\SZ18_state_taxes.dta", replace
keep if year >= 1980 & year <= 2010
drop if fips_state == 11 | fips_state == 0 | fips_state > 56
xtset fips_state year
* Balanced panel: 50 states, years from 1980 to 2010.
* Treated group: 15 states.
* Control group: 35 other states.
* Construct dependent variables
gen log_rev_corptax = 100*ln(rev_corptax+1)
gen log_gdp = 100*ln(GDP)
g r_g = 100*rev_corptax/GDP
* Screen out the tax change with a pre-determined threshold
local threshold = 0.5
gen ch_corporate_rate = corporate_rate - L1.corporate_rate
replace ch_corporate_rate = 0 if abs(ch_corporate_rate) <= `threshold'
gen ch_corporate_rate_inc = (ch_corporate_rate > 0 & !missing(ch_corporate_rate))
gen ch_corporate_rate_dec = (ch_corporate_rate < 0 & !missing(ch_corporate_rate))
* Contruct treatment dummies
** Static
gen change_year = year if ch_corporate_rate_dec==1
bysort fips_state: egen tchange_year = min(change_year)
gen treated = (year >= tchange_year)
** Dynamic
gen period = year - tchange_year
gen Dn5 = (period < -4)
forvalues i = 4(-1)1 {
gen Dn`i' = (period == -`i')
}
forvalues i = 0(1)5 {
gen D`i' = (period == `i')
}
gen D6 = (period >= 6 & period != .)
* Static DID
local ylist = "log_rev_corptax log_gdp r_g"
local i = 1
foreach yvar in `ylist' {
quietly{
eststo areg`i': areg `yvar' treated i.fips_state, absorb(year) cluster(fips_state)
eststo hdreg`i': reghdfe `yvar' treated, absorb(fips_state year) cluster(fips_state)
eststo sdid`i': sdid `yvar' fips_state year treated, vce(bootstrap) seed(2018)
local ++i
}
}
estout areg1 hdreg1 sdid1 areg2 hdreg2 sdid2 areg3 hdreg3 sdid3, keep(treated) ///
coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///
starlevels(* .1 ** .05 *** .01) legend ///
stats(N r2_a, nostar labels("Observations" "Adj. R2") fmt("%9.0fc" 3))
* Dynamic DID (sdid cannot do this)
local ylist = "log_rev_corptax log_gdp r_g"
local i = 1
foreach yvar in `ylist' {
quietly{
eststo areg`i': areg `yvar' Dn5-Dn2 D0-D6 i.fips_state, absorb(year) cluster(fips_state)
eststo hdreg`i': reghdfe `yvar' Dn5-Dn2 D0-D6, absorb(fips_state year) cluster(fips_state)
local ++i
}
}
estout areg1 hdreg1 areg2 hdreg2 areg3 hdreg3, keep(D*) ///
coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///
starlevels(* .1 ** .05 *** .01) legend ///
stats(N r2_a, nostar labels("Observations" "Adj. R2") fmt("%9.0fc" 3))
**# Stata 18 new command: xthdidregress
local controls = "FederalIncomeasStateTaxBase sales_wgt throwback FedIncomeTaxDeductible Losscarryforward FranchiseTax"
local ps_var = "FederalIncomeasStateTaxBase sales_wgt throwback FedIncomeTaxDeductible FranchiseTax"
xthdidregress aipw (log_gdp `controls') (treated `ps_var'), group(fips_state) vce(cluster fips_state)
* Visualizing ATT for each cohort
estat atetplot
graph export "$figdir\SZ18_cohort_ATT.pdf", replace
* Visualizing ATT over cohort
estat aggregation, cohort ///
graph(xlab(, angle(45) labsize(small)) legend(rows(1) position(6)))
graph export "$figdir\SZ18_cohort_agg_ATT.pdf", replace
* Visualizing ATT over time
estat aggregation, time ///
graph(xlab(, angle(45) labsize(small)) legend(rows(1) position(6)))
graph export "$figdir\SZ18_time_agg_ATT.pdf", replace
* Visualizing dynamic effects
estat aggregation, dynamic(-5/6) ///
graph( ///
title("Dynamic Effects on Log State GDP", size(medlarge)) ///
xlab(, labsize(small) nogrid) ///
legend(rows(1) position(6)) ///
xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///
)
graph export "$figdir\SZ18_dynamic_ATT.pdf", replace