事件研究法用于DID的经典文献"环境规制"论文数据和程序

凡是搞计量经济的，都关注这个号了

邮箱：[e-mail]

所有计量经济圈方法论丛的code程序, 宏微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

这篇文章主要研究环保法令（政策）有没有改善空气、水以及婴儿死亡率，作者采用的是与事件研究法类似的DID估计策略。

作者采用的是两步法：1. 用事件研究法得到在事件窗口期空气、水等平均污染程度，2. 然后就用第一步得到的平均污染度对政策变量+时间趋势+政策变量*时间趋势做回归。

部分程序，值得深究，比如这里使用的reg做DID

完整程序、数据放在社群，有需要可以参考


// An highlighted block

Objective: The objective of this program is to generate Table 3 of the paper. This program
	performs mean shift and trend break analysis of air pollution policies using both a
	single-stage standard DiD approach as well as a two-stage event study and trend break
	approach. 


Structure:
0. Stata Setup
1. Data prep: Generate policy event-years on air dataset
2. Run single-stage and two-stage regressions
	a. CAT policy
	b. SCAP policy
3. Clean-up

Specifications:
1. Single-stage approach (standard difference-in-differences regressions):
	a. Mean shift (includes post*treatment only)
	b. Mean shift (includes post*treatment and city-year time trends)
	c. Trendbreak (included post*treatment, city-year time trends, and post*treatment*time-trend)
2. Two-stage approach (event study followed by mean-shift or trend-break regressions)
	a. Mean shift: reg taub post*treatment [aw=1/tause]
	b. Mean shift w/ time trend: reg taub post*treatment timetrend [aw=1/tause]
	c. Trend break: reg taub post*treatment timetrend post*treatment*timetrend [1/tause]
*/

***********************************************************************************
***********************************************************************************
* 0. Setup
***********************************************************************************

* Reset memory options
clear all
set more off
set mem 500m
set matsize 10000

* Establish log
cap log close
log using "${mypath}Log/Table3_01102014.txt", replace t

***********************************************************************************
* 1. Data prep: Generate policy event-years on air dataset
***********************************************************************************

* Supreme Court Action Plan and Catalytic Converter policies *
foreach pollutant in spm so2 no2 { 
	use "${mypath}Data/Air Data/Final Data/combined.dta", clear
	egen city_id = group(city)
	bys state city: egen count = count(e_`pollutant'_mean)

	* Generate Supreme Court Action Plan taus
	cap drop temp*
	g temp = year if actionplan_sc == 1
	egen temp2 = min(temp), by(city_id)
	egen temp3 = min(year), by(city_id)
	g tauSC = year - temp2 if temp2 > temp3 & e_`pollutant'_mean<.
	g neveradoptSC = temp2 >= .
	replace tauSC = 0 if neveradoptSC == 1 & e_`pollutant'_mean < . // We code the non-adopting cities so that they are always in event year 0
	g tempx = 1 if actionplan_sc==1 & e_`pollutant'_mean<.		// the following 3 lines of code determine whether there are any cities which enacted the policy but are missing pollution data for ALL post-policy years
	egen tempy = min(tempx) if temp2<., by(city_id)
	ta tempy if temp2<., missing		// no cities fit this description, so no more changes need be made

	g temp4 = tauSC >= 3 & tauSC < .
	g temp5 = tauSC <= -3 & tauSC < .
	egen Mtemp4 = max(temp4), by(city_id) 
	egen Mtemp5 = max(temp5), by(city_id) 
	g useSC = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptSC == 1 & count > 1)

	* Generate Catalytic Converter Policy taus
	cap drop Mtemp* temp*
	g temp = year if catconverter == 1
	egen temp2 = min(temp), by(city_id) 
	egen temp3 = min(year), by(city_id) 
	g tauCAT = year - temp2 if temp2 > temp3 & e_`pollutant'_mean<.
	g neveradoptCAT = temp2 >= .
	replace tauCAT = 0 if neveradoptCAT == 1 & e_`pollutant'_mean < . // We code the non-adopting cities so that they are always in event year 0
	g tempx = 1 if catconverter==1 & e_`pollutant'_mean<.
	egen tempy = min(tempx) if temp2<., by(city_id)
	ta tempy if temp2<., missing
	replace neveradoptCAT = 1 if tempy>=. & temp2<.
	replace tauCAT = . if tempy>=. & temp2<.

	* The following lines restrict the sample to only those cities with at an observation at least three years prior to and after
	*	policy adoption. We also want to estimate Eqns (1) and (2) with an unrestricted sample. For this reason, we'll save restricted
	*	and unrestricted versions of the dataset.
	g temp4 = tauCAT >= 3 & tauCAT < .
	g temp5 = tauCAT <= -3 & tauCAT < .
	egen Mtemp4 = max(temp4), by(city_id)
	egen Mtemp5 = max(temp5), by(city_id)
	g useCAT = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptCAT == 1 & count > 1)

	drop if useSC==0 | useCAT==0
	drop if e_`pollutant'_mean >= .
	ta tauSC
	ta tauCAT
	g one = 1

	* Generate individual event-year dummies for Supreme Court Action Plans
	forv tau = 7(-1)1 {
		g tauSCm`tau' = tauSC == -`tau'
		la var tauSCm`tau' "This obs is `tau' years before action plan began"
	}
	forv tau = 0/3 {
		g tauSC`tau' = tauSC == `tau'
		la var tauSC`tau' "This obs is `tau' years after action plan began"
	}
	g tauSCL = tauSC<-7
	g tauSCR = tauSC>3 & tauSC<.

	* Generate individual event-year dummies for Catalytic Converter Policy
	forv tau = 7(-1)1 {
		g tauCATm`tau' = tauCAT == -`tau'
		la var tauCATm`tau' "This obs is `tau' years before catalytic converters mandated"
	}
	forv tau = 0/9 {
		g tauCAT`tau' = tauCAT == `tau'
		la var tauCAT`tau' "This obs is `tau' years after catalytic converters mandated"
	}
	g tauCATL = tauCAT<-7
	g tauCATR = tauCAT>9 & tauCAT<.
	
	* Save pollutant specific dataset
	save "${mypath}Data/Air Data/Final Data/combined_taus_`pollutant'_stddid.dta", replace	
}

***********************************************************************************
* 2. Run single-stage and two-stage regressions
***********************************************************************************

* Catalytic Converter policy *
local placeholder "replace" // we use this to control the outreg command
foreach pollutant in spm so2 no2 {
	use "${mypath}Data/Air Data/Final Data/combined_taus_`pollutant'_stddid.dta", clear
* Single-Stage Regressions:
	gen scaprange = (tauSC >= -7 & tauSC <= 3) // the two-stage regressions only include SC event years between -7 & 3
	gen scap = (tauSC >= 0) & (neveradoptSC == 0)
	gen scappolicy = scap*scaprange
	gen scaptau = tauSC*scaprange
	gen scap_trend = scap*scaprange*tauSC
	replace scap_trend = scap_trend + 1 if (tauSC >= 0 & tauSC <= 3)
	gen catrange = (tauCAT >= -7 & tauCAT <= 9) // the two-stage regressions only include CAT event years between -7 & 9
	gen catconv = (tauCAT >= 0) & (neveradoptCAT == 0)
	gen catconvpolicy = catconv*catrange
	gen catconvtau = tauCAT*catrange
	gen catconv_trend = catconv*catrange*tauCAT
	replace catconv_trend = catconv_trend + 1 if (tauCAT >= 0 & tauCAT <= 9)

*	Mean shift (includes post*policy plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR catrange catconvpolicy tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban], cluster(city)
		outreg2 catconvpolicy using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 3a") aster(coef) replace
*	Mean shift (includes post*policy and city-year time trend, plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR scaptau catrange catconvpolicy tauCATL tauCATR catconvtau lit_urban mean i.city i.year [aw = pop_urban], cluster(city)
		outreg2 catconvpolicy catconvtau using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 3b") aster(coef) append
*	Trendbreak (includes post*treated, city-year time trend, and post*treated*time-trend, plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR scaptau scap_trend catrange catconvpolicy tauCATL tauCATR catconvtau catconv_trend lit_urban mean i.city i.year [aw = pop_urban], cluster(city)
		lincom catconvpolicy + 5*catconv_trend
		local coef = r(estimate)
		local coef = round(`coef', .01)
		local coef = "`coef'"
		if length("`coef'") < 3 {
			local coef = "`coef'" + "0"
		}
		local coef = substr("`coef'", 1, 6)
		local se = r(se)
		local tstat = r(estimate)/r(se)
		local pval = tprob(r(df), abs(`tstat'))
		local pval = round(`pval', .01)
		if `pval' <= .01 {
			local coef = "`coef'" + "***"
		}
		if `pval' > .01 & `pval' <= .05 {
			local coef = "`coef'" + "**"
		}
		if `pval' > .05 & `pval' <= .1 {
			local coef = "`coef'" + "*"
		}
		if length("`pval'") < 3 {
			local pval = "`pval'" + "0"
		}
		outreg2 catconvpolicy catconvtau catconv_trend using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 3c") aster(coef) addtext("5-Year Effect", "`coef'", "p-value", "[`pval']") append
		
* Two-Stage Regressions:
	* Recall that, for the purposes of the single-stage regressions, control cities were coded as event year 0. Revert them to
	*	event year missing.
	replace tauCAT = . if neveradoptCAT == 1 & e_`pollutant'_mean < .
	replace tauSC = . if neveradoptSC == 1 & e_`pollutant'_mean < .
	replace tauCAT0 = 0 if neveradoptCAT == 1 & e_`pollutant'_mean < .
	replace tauSC0 = 0 if neveradoptSC == 1 & e_`pollutant'_mean < .	
*	First stage (Event study)
	qui xi: reg e_`pollutant'_mean tauSCm7-tauSCm1 tauSC0-tauSC3 tauSCL tauSCR tauCATm7-tauCATm1 tauCAT0-tauCAT9 tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban]

	* Prepare sigmat-hats
	gen taub = .
	gen tause = .
	collapse (mean) taub tause tauCATm7-tauCAT9, by(tauCAT)
	drop if tauCAT > 9 | tauCAT < -7
	foreach k of varlist tauCATm7-tauCATm1 tauCAT0-tauCAT9 {
		replace taub = _b[`k'] if `k' == 1
		replace tause = _se[`k'] if `k' == 1
	}
	drop tauCATm7-tauCAT9
	gen catconv = (tauCAT >= 0)
	gen tau_trend = tauCAT + 8
	gen catconv_trend = catconv*tauCAT
	replace catconv_trend = catconv_trend + 1 if tauCAT >= 0
	
* Save the dataset of sigma-hats (for use with breakpoint tests)
	save "${mypath}Data/Air Data/Final Data/sigmas_`pollutant'.dta", replace		
	
*	Second stage (mean shift and trendbreak regressions)	
	* Trendbreak regressions - Eqn (2a)
	display "Outcome: `pollutant'; Policy: Catalytic Converters; Sample: all city*years"
	regress taub catconv [aw = 1/tause]
		outreg2 catconv using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 2a") aster(coef) append	

	* Trendbreak regressions - Eqn (2b)
	reg taub tauCAT catconv [aw=1/tause]
		outreg2 catconv tauCAT using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 2b") aster(coef) append
	
	* Trendbreak regressions - Eqn (2c)
	reg taub tauCAT catconv catconv_trend [aw=1/tause]
		lincom catconv + 5*catconv_trend
			local coef = r(estimate)
			local coef = round(`coef', .01)
			local coef = "`coef'"
			if length("`coef'") < 3 {
				local coef = "`coef'" + "0"
			}
			local coef = substr("`coef'", 1, 6)
			local tstat = r(estimate)/r(se)
			local pval = tprob(r(df), abs(`tstat'))
			local pval = round(`pval', .01)
			if `pval' <= .01 {
				local coef = "`coef'" + "***"
			}
			if `pval' > .01 & `pval' <= .05 {
				local coef = "`coef'" + "**"
			}
			if `pval' > .05 & `pval' <= .1 {
				local coef = "`coef'" + "*"
			}
			if length("`pval'") < 3 {
				local pval = "`pval'" + "0"
			}
			outreg2 catconv tauCAT catconv_trend using "${mypath}Output/T3/Table3-CAT-`pollutant'.xls", dec(2) se par ctitle("Eqn 2c") aster(coef) addtext("5-Year Effect", "`coef'", "p-value", "[`pval']") append
}
*/
*SCAP policy *
local placeholder "replace" // we use this to control the outreg command
foreach pollutant in spm so2 no2 {
	use "${mypath}Data/Air Data/Final Data/combined_taus_`pollutant'_stddid.dta", clear
* Single-Stage Regressions:
	gen scaprange = (tauSC >= -7 & tauSC <= 3) // the two-stage regressions only include SC event years between -7 & 3
	gen scap = (tauSC >= 0) & (neveradoptSC == 0)
	gen scappolicy = scap*scaprange
	gen scaptau = tauSC*scaprange
	gen scap_trend = scap*scaprange*tauSC
	replace scap_trend = scap_trend + 1 if (tauSC >= 0 & tauSC <= 3)
	gen catrange = (tauCAT >= -7 & tauCAT <= 9) // the two-stage regressions only include CAT event years between -7 & 9
	gen catconv = (tauCAT >= 0) & (neveradoptCAT == 0)
	gen catconvpolicy = catconv*catrange
	gen catconvtau = tauCAT*catrange
	gen catconv_trend = catconv*catrange*tauCAT
	replace catconv_trend = catconv_trend + 1 if (tauCAT >= 0 & tauCAT <= 9)

*	Mean shift (includes post*policy plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR catrange catconvpolicy tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban], cluster(city)
		outreg2 scappolicy using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 3a") aster(coef) replace
*	Mean shift (includes post*policy and city-year time trend, plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR scaptau catrange catconvpolicy tauCATL tauCATR catconvtau lit_urban mean i.city i.year [aw = pop_urban], cluster(city)
		outreg2 scappolicy scaptau using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 3b") aster(coef) append
*	Trendbreak (includes post*treated, city-year time trend, and post*treated*time-trend, plus the usual controls)
	xi: reg e_`pollutant'_mean scaprange scappolicy tauSCL tauSCR scaptau scap_trend catrange catconvpolicy tauCATL tauCATR catconvtau catconv_trend lit_urban mean i.city i.year [aw = pop_urban], cluster(city)	
		lincom scappolicy + 5*scap_trend
		local coef = r(estimate)
		local coef = round(`coef', .01)
		local coef = "`coef'"
		if length("`coef'") < 3 {
			local coef = "`coef'" + "0"
		}
		local coef = substr("`coef'", 1, 6)
		local se = r(se)
		local tstat = r(estimate)/r(se)
		local pval = tprob(r(df), abs(`tstat'))
		local pval = round(`pval', .01)
		if `pval' <= .01 {
			local coef = "`coef'" + "***"
		}
		if `pval' > .01 & `pval' <= .05 {
			local coef = "`coef'" + "**"
		}
		if `pval' > .05 & `pval' <= .1 {
			local coef = "`coef'" + "*"
		}
		if length("`pval'") < 3 {
			local pval = "`pval'" + "0"
		}
		outreg2 scappolicy scaptau scap_trend using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 3c") aster(coef) addtext("5-Year Effect", "`coef'", "p-value", "[`pval']") append
		
* Two-Stage Regressions:
	* Recall that, for the purposes of the single-stage regressions, control cities were coded as event year 0. Revert them to
	*	event year missing.
	replace tauCAT = . if neveradoptCAT == 1 & e_`pollutant'_mean < .
	replace tauSC = . if neveradoptSC == 1 & e_`pollutant'_mean < .
	replace tauCAT0 = 0 if neveradoptCAT == 1 & e_`pollutant'_mean < .
	replace tauSC0 = 0 if neveradoptSC == 1 & e_`pollutant'_mean < .	
*	First stage (Event study)
	qui xi: reg e_`pollutant'_mean tauSCm7-tauSCm1 tauSC0-tauSC3 tauSCL tauSCR tauCATm7-tauCATm1 tauCAT0-tauCAT9 tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban]

	* Prepare sigmat-hats
	gen taub = .
	gen tause = .
	collapse (mean) taub tause tauSCm7-tauSC3, by(tauSC)
	drop if tauSC > 3 | tauSC < -7
	foreach k of varlist tauSCm7-tauSCm1 tauSC0-tauSC3 {
		replace taub = _b[`k'] if `k' == 1
		replace tause = _se[`k'] if `k' == 1
	}
	drop tauSCm7-tauSC3
	gen scap = (tauSC >= 0)
	gen tau_trend = tauSC + 8
	gen scap_trend = scap*tauSC
	replace scap_trend = scap_trend + 1 if tauSC >= 0
	
*	Second stage (mean shift and trendbreak regressions)	
	* Trendbreak regressions - Eqn (2a)
	display "Outcome: `pollutant'; Policy: SCAP; Sample: all city*years"
	regress taub scap [aw = 1/tause] 
		outreg2 scap using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 2a") aster(coef) append	

	* Trendbreak regressions - Eqn (2b)
	reg taub tauSC scap [aw=1/tause]
		outreg2 scap tauSC using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 2b") aster(coef) append
	
	* Trendbreak regressions - Eqn (2c)
	reg taub tauSC scap scap_trend [aw=1/tause]
		lincom scap + 5*scap_trend
			local coef = r(estimate)
			local coef = round(`coef', .01)
			local coef = "`coef'"
			if length("`coef'") < 3 {
				local coef = "`coef'" + "0"
			}
			local coef = substr("`coef'", 1, 6)
			local tstat = r(estimate)/r(se)
			local pval = tprob(r(df), abs(`tstat'))
			local pval = round(`pval', .01)
			if `pval' <= .01 {
				local coef = "`coef'" + "***"
			}
			if `pval' > .01 & `pval' <= .05 {
				local coef = "`coef'" + "**"
			}
			if `pval' > .05 & `pval' <= .1 {
				local coef = "`coef'" + "*"
			}
			if length("`pval'") < 3 {
				local pval = "`pval'" + "0"
			}
			outreg2 scap tauSC scap_trend using "${mypath}Output/T3/Table3-SCAP-`pollutant'.xls", dec(2) se par ctitle("Eqn 2c") aster(coef) addtext("5-Year Effect", "`coef'", "p-value", "[`pval']") append
}

***********************************************************************************
* 3. Clean-up
***********************************************************************************

* Erase intermediate datasets
foreach pollutant in spm so2 no2 {
	erase "${mypath}Data/Air Data/Final Data/combined_taus_`pollutant'_stddid.dta"
}

log close

事件图Event plot

部分程序，值得深究

// An highlighted block
/*
Objective: The objective of this program is to generate Figure 5 in the paper. Figure 5
	plots the air policy event-study estimates developed in the first stage of the two-stage
	regression analysis. 
	
Structure:
1. Figure 5. Air Pollution Event Study Figure

For each of the Air Pollution measures, the graphs are generated by:

1. Cleaning and generating necessary variables for analysis
2. Generating individual graphs
3. Combining graphs into one

Notes:
1. The policies considered are SCAP and Catalytic Converters
*/

***********************************************************************************

clear all
set more off
set mem 500m

* Establish log
cap log close
log using "${mypath}Log/Figure5_01102014.txt", replace t

*ssc install estout, replace

***********************************************
* Event Studies

*** Air Pollution ***

/* PM */
u "${mypath}Data/Air Data/Final Data/combined.dta", clear
cd "${mypath}Output/F5"

egen city_id = group(city)
bys state city: egen count = count(e_spm_mean)

* Preparing variables for analysis

cap drop temp*
g temp = year if actionplan_sc == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauSC = year - temp2 if temp2 > temp3 & e_spm_mean<.
g neveradoptSC = temp2 >= .
g tempx = 1 if actionplan_sc==1 & e_spm_mean<.
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing

g temp4 = tauSC >= 3 & tauSC < .
g temp5 = tauSC <= -3 & tauSC < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id) 
g useSC = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptSC == 1 & count > 1)

drop Mtemp* temp*
g temp = year if catconverter == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauCAT = year - temp2 if temp2 > temp3 & e_spm_mean<.
g neveradoptCAT = temp2 >= .
g tempx = 1 if catconverter==1 & e_spm_mean<.
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing
replace neveradoptCAT = 1 if tempy>=. & temp2<.
replace tauCAT = . if tempy>=. & temp2<.

g temp4 = tauCAT >= 3 & tauCAT < .
g temp5 = tauCAT <= -3 & tauCAT < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id)
g useCAT = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptCAT == 1 & count > 1)

drop if useSC==0 | useCAT==0
drop if e_spm_mean >= .
ta tauSC
ta tauCAT
g one = 1

forv tau = 7(-1)1 {
	g tauSCm`tau' = tauSC == -`tau'
	la var tauSCm`tau' "This obs is `tau' years before action plan began"
}
forv tau = 0/3 {
	g tauSC`tau' = tauSC == `tau'
	la var tauSC`tau' "This obs is `tau' years after action plan began"
}
g tauSCL = tauSC<-7
g tauSCR = tauSC>3 & tauSC<.
forv tau = 7(-1)1 {
	g tauCATm`tau' = tauCAT == -`tau'
	la var tauCATm`tau' "This obs is `tau' years before catalytic converters mandated"
}
forv tau = 0/9 {
	g tauCAT`tau' = tauCAT == `tau'
	la var tauCAT`tau' "This obs is `tau' years after catalytic converters mandated"
}
g tauCATL = tauCAT<-7
g tauCATR = tauCAT>9 & tauCAT<.

xi: reg e_spm_mean tauSCm7-tauSCm2 tauSC0-tauSC3 tauSCL tauSCR tauCATm7-tauCATm2 tauCAT0-tauCAT9 tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban]
mat M = e(b)'
svmat M
replace M1 = . if _n >= 11
g n = _n - 8 in 1/11
tsset one n 
g temp6 = L.M1
replace M1 = temp6 if _n >= 8
replace M1 = 0 if _n == 7
tw connected M1 n if _n <= 11, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on PM", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Action Plan Mandate", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Supreme Court Action Plan") ///
	xlab(-7 -6 -3 0 3) ///
	saving(taus_apPM, replace)
drop n temp6
mat N = e(b)'
svmat N
replace N1 = . if _n >= 29 | _n <= 12
g n = _n - 20 in 13/29
tsset one n 
g temp6 = L.N1
replace N1 = temp6 if _n >= 8
replace N1 = 0 if _n == 7
tw connected N1 n if _n <= 17, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on PM", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Catalytic Converters Mandated", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Catalytic Converters") ///
	xlab(-7 -6 -3 0 3 6 9) ///
	saving(taus_catconPM, replace)			
			
* Combine PM graphs
graph combine taus_apPM.gph taus_catconPM.gph, cols(1) scheme(s1color) ysize(5) xsize(3) iscale(*1.2)
graph export F5a.eps, replace	

/* SO2 */
u "${mypath}Data/Air Data/Final Data/combined.dta", clear

egen city_id = group(city)
bys state city: egen count = count(e_so2_mean)

cap drop temp*
g temp = year if actionplan_sc == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauSC = year - temp2 if temp2 > temp3 & e_so2_mean<.
g neveradoptSC = temp2 >= .
g tempx = 1 if actionplan_sc==1 & e_so2_mean<.		// the following 3 lines of code determine whether there are any cities which enacted the policy but are missing pollution data for ALL post-policy years
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing		// no cities fit this description, so no more changes need be made

g temp4 = tauSC >= 3 & tauSC < .
g temp5 = tauSC <= -3 & tauSC < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id) 
g useSC = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptSC == 1 & count > 1)

drop Mtemp* temp*
g temp = year if catconverter == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauCAT = year - temp2 if temp2 > temp3 & e_so2_mean<.
g neveradoptCAT = temp2 >= .
g tempx = 1 if catconverter==1 & e_so2_mean<.		// the following 3 lines of code determine whether there are any cities which enacted the policy but are missing pollution data for ALL post-policy years
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing		// Silvassa again. Include it as a 'neveradopt'.
replace neveradoptCAT = 1 if tempy>=. & temp2<.
replace tauCAT = . if tempy>=. & temp2<.

g temp4 = tauCAT >= 3 & tauCAT < .
g temp5 = tauCAT <= -3 & tauCAT < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id)
g useCAT = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptCAT == 1 & count > 1)

drop if useSC==0 | useCAT==0
drop if e_so2_mean >= .
ta tauSC
ta tauCAT
g one = 1

forv tau = 7(-1)1 {
	g tauSCm`tau' = tauSC == -`tau'
	la var tauSCm`tau' "This obs is `tau' years before action plan began"
}
forv tau = 0/3 {
	g tauSC`tau' = tauSC == `tau'
	la var tauSC`tau' "This obs is `tau' years after action plan began"
}
g tauSCL = tauSC<-7
g tauSCR = tauSC>3 & tauSC<.
forv tau = 7(-1)1 {
	g tauCATm`tau' = tauCAT == -`tau'
	la var tauCATm`tau' "This obs is `tau' years before catalytic converters made mandatory"
}
forv tau = 0/9 {
	g tauCAT`tau' = tauCAT == `tau'
	la var tauCAT`tau' "This obs is `tau' years after catalytic converters made mandatory"
}
g tauCATL = tauCAT<-7
g tauCATR = tauCAT>9 & tauCAT<.

xi: reg e_so2_mean tauSCm7-tauSCm2 tauSC0-tauSC3 tauSCL tauSCR tauCATm7-tauCATm2 tauCAT0-tauCAT9 tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban]
mat M = e(b)'
svmat M
replace M1 = . if _n >= 11
g n = _n - 8 in 1/11
tsset one n 
g temp6 = L.M1
replace M1 = temp6 if _n >= 8
replace M1 = 0 if _n == 7
tw connected M1 n if _n <= 11, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on SO2", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Action Plan Mandate", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Supreme Court Action Plan") ///
	xlab(-7 -6 -3 0 3) ///
	saving(taus_apSO, replace)
drop n temp6
mat N = e(b)'
svmat N
replace N1 = . if _n >= 29 | _n <= 12
g n = _n - 20 in 13/29
tsset one n 
g temp6 = L.N1
replace N1 = temp6 if _n >= 8
replace N1 = 0 if _n == 7
tw connected N1 n if _n <= 17, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on SO2", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Catalytic Converters Mandated", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Catalytic Converters") ///
	xlab(-7 -6 -3 0 3 6 9) ///
	saving(taus_catconSO, replace)			

* Combine so2 graphs
graph combine taus_apSO.gph taus_catconSO.gph, cols(1) scheme(s1color) ysize(5) xsize(3) iscale(*1.2)
graph export F5b.eps, replace	

/* NO2 */
u "${mypath}Data/Air Data/Final Data/combined.dta", clear

egen city_id = group(city)
bys state city: egen count = count(e_no2_mean)

cap drop temp*
g temp = year if actionplan_sc == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauSC = year - temp2 if temp2 > temp3 & e_no2_mean<.
g neveradoptSC = temp2 >= .
g tempx = 1 if actionplan_sc==1 & e_no2_mean<.		// the following 3 lines of code determine whether there are any cities which enacted the policy but are missing pollution data for ALL post-policy years
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing		// no cities fit this description, so no more changes need be made

g temp4 = tauSC >= 3 & tauSC < .
g temp5 = tauSC <= -3 & tauSC < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id) 
g useSC = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptSC == 1 & count > 1)
	
drop Mtemp* temp*
g temp = year if catconverter == 1
egen temp2 = min(temp), by(city_id) 
egen temp3 = min(year), by(city_id) 
g tauCAT = year - temp2 if temp2 > temp3 & e_no2_mean<.
g neveradoptCAT = temp2 >= .
g tempx = 1 if catconverter==1 & e_no2_mean<.		// the following 3 lines of code determine whether there are any cities which enacted the policy but are missing pollution data for ALL post-policy years
egen tempy = min(tempx) if temp2<., by(city_id)
ta tempy if temp2<., missing		// Silvassa.
replace neveradoptCAT = 1 if tempy>=. & temp2<.
replace tauCAT = . if tempy>=. & temp2<.

g temp4 = tauCAT >= 3 & tauCAT < .
g temp5 = tauCAT <= -3 & tauCAT < .
egen Mtemp4 = max(temp4), by(city_id) 
egen Mtemp5 = max(temp5), by(city_id)
g useCAT = (Mtemp4 == 1 & Mtemp5 == 1 & count > 1) | (neveradoptCAT == 1 & count > 1)

drop if useSC==0 | useCAT==0
drop if e_no2_mean >= .
ta tauSC
ta tauCAT
g one = 1

forv tau = 7(-1)1 {
	g tauSCm`tau' = tauSC == -`tau'
	la var tauSCm`tau' "This obs is `tau' years before action plan began"
}
forv tau = 0/3 {
	g tauSC`tau' = tauSC == `tau'
	la var tauSC`tau' "This obs is `tau' years after action plan began"
}
g tauSCL = tauSC<-7
g tauSCR = tauSC>3 & tauSC<.
forv tau = 7(-1)1 {
	g tauCATm`tau' = tauCAT == -`tau'
	la var tauCATm`tau' "This obs is `tau' years before catalytic converters made mandatory"
}
forv tau = 0/9 {
	g tauCAT`tau' = tauCAT == `tau'
	la var tauCAT`tau' "This obs is `tau' years after catalytic converters made mandatory"
}
g tauCATL = tauCAT<-7
g tauCATR = tauCAT>9 & tauCAT<.

xi: reg e_no2_mean tauSCm7-tauSCm2 tauSC0-tauSC3 tauSCL tauSCR tauCATm7-tauCATm2 tauCAT0-tauCAT9 tauCATL tauCATR lit_urban mean i.city i.year [aw = pop_urban]
mat M = e(b)'
svmat M
replace M1 = . if _n >= 11
g n = _n - 8 in 1/11
tsset one n 
g temp6 = L.M1
replace M1 = temp6 if _n >= 8
replace M1 = 0 if _n == 7
tw connected M1 n if _n <= 11, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on NO2", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Action Plan Mandate", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Supreme Court Action Plan") ///
	xlab(-7 -6 -3 0 3) ///
	saving(taus_apNO, replace)
drop n temp6
mat N = e(b)'
svmat N
replace N1 = . if _n >= 29 | _n <= 12
g n = _n - 20 in 13/29
tsset one n 
g temp6 = L.N1
replace N1 = temp6 if _n >= 8
replace N1 = 0 if _n == 7
tw connected N1 n if _n <= 17, lwidth(medthick) ///
	sort yline(0, lpattern(dash)) ///
	ytit("Effect on NO2", margin(medsmall)) ///
	xline(0, lpattern(dash)) ///
	xtit("Years since Catalytic Converters Mandated", margin(medsmall) size(small)) ///
	scheme(s1color) ///
	subtit("Policy:  Catalytic Converters") ///
	xlab(-7 -6 -3 0 3 6 9) ///
	saving(taus_catconNO, replace)			

* Combine no2 graphs
graph combine taus_apNO.gph taus_catconNO.gph, cols(1) scheme(s1color) ysize(5) xsize(3) iscale(*1.2)
graph export F5c.eps, replace	

log close

太长了，数据和程序都放到了咱们社群里，有需要可以下载参看。

下面这些短链接文章属于合集，可以收藏起来阅读，不然以后都找不到了。

2年，计量经济圈公众号近1000篇文章，

Econometrics Circle

数据系列：空间矩阵 | 工企数据 | PM2.5 | 市场化指数 | CO2数据 |  夜间灯光 | 官员方言  | 微观数据

计量系列：匹配方法 | 内生性 | 工具变量 | DID | 面板数据 | 常用TOOL | 中介调节  | 时间序列 | RDD断点 | 合成控制 | 

数据处理：Stata | R | Python | 缺失值 | 

干货系列：能源环境 | 效率研究 | 空间计量 | 国际经贸 | 计量软件 | 商科研究 | 机器学习 | SSCI | CSSCI

计量经济圈组织了一个计量社群，有如下特征：热情互助最多、前沿趋势最多、社科资料最多、社科数据最多、科研牛人最多、海外名校最多。因此，建议积极进取和有强烈研习激情的中青年学者到社群交流探讨，始终坚信优秀是通过感染优秀而互相成就彼此的。