clear
set more off

local end_time = 9

use "$temp/final_round_2", clear

* Add second round turnout data
merge n:1 departement_code using "$temp/departement_turnout_2020_round_2"
drop if _merge == 2 // Drop departements out of metropolitan France
drop _merge

* Get weather conditions on the day of the second round
merge 1:1 departement_code date using "$temp/departement_meteo"
drop if _merge == 2
drop _merge
pca precip temp_mean // if date == mdy(03,15,2020) // wind_mean humid_mean
predict meteo if date == mdy(06,28,2020), score
drop meteo_elec
egen meteo_elec = max(meteo), by(departement_code)
drop meteo

* Document departements for which predictions are not available
********************************************************************************
preserve

distinct id_dep if prediction_`end_time'_log3 == .

gen insample = (prediction_`end_time'_log3 != .)

summ population if date == mdy(07,`end_time',2020) & insample == 0
local a = r(sum)
summ population if date == mdy(07,`end_time',2020)
di `a'/r(sum)

restore
********************************************************************************

distinct id_dep if prediction_`end_time'_log3 != .

* Create groups of departement as a function of epidemics on June 28 (+ 11 days lag in measurement)
*******************************************************************************

gen nbre_hospit_corona_cum_9_pc = (nbre_hospit_corona_cum_9/population)*100000
	
* Bottom 25% and bottom third with 32-day cumulated hospitalizations
	local third = 1/3*100
	_pctile nbre_hospit_corona_cum_9_pc if date == mdy(07,09,2020) & prediction_`end_time'_log3!=.  , p(25,  `third') 
	gen temp_high25 = (nbre_hospit_corona_cum_9_pc > r(r1)) if date == mdy(07,`end_time',2020)
	gen temp_high33 = (nbre_hospit_corona_cum_9_pc > r(r2)) if date == mdy(07,`end_time',2020)
	bys id_dep: egen high25 = mean(temp_high25)
	bys id_dep: egen high33 = mean(temp_high33)
	drop temp*

* At least 7 days of increase before end of fit
	tsset id_dep date
	local depth = 7
	gen temp0 = (nbre_hospit_corona != 0) if nbre_hospit_corona != .
	gen temp1 = 0
	replace temp1 = L.temp1 + temp0 if L.temp1 != .
	gen temp2 = (temp1>`depth') if date == mdy(07,09,2020)
	bys id_dep: egen act`depth' = mean(temp2)
	drop temp*

* Define sample period and construct dependent variable
********************************************************************************	
	
gen diff_`end_time'_log3_pc = (nbre_hospit_corona_cum_9 - prediction_`end_time'_log3)/population *100000

local date_start = mdy(06,14,2020)
local date_end = mdy(07,`end_time',2020) + 7

keep if inrange(date,`date_start',`date_end')

replace turnout_2020 = votes_round2/voters_round2

* Construct weights
********************************************************************************

* Per capita standard error of prediction
gen se_pc = (prediction_`end_time'_log3_se/population)*100000

gen weight = 1
replace weight = 1/(se_pc) if date > mdy(07,`end_time',2020)

* Normalize sum of weights to 1 within each day (useless given the estimation framework)
egen h = sum(weight), by(date)
replace weight = weight/h
drop  h
	
local treat = "high33"

* Create interaction terms
forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(06,13,2020)

gen coeff_date_`j' = (date == `i') 
gen coeff_turnout_`j' = turnout_2020 * (date == `i')   * (1 - `treat')
gen coeff_high_`j' = (date == `i') * `treat'
gen coeff_high_turnout_`j' = turnout_2020 * (date == `i') * `treat'

* Controls
gen coeff_density_`j' = density * (date == `i') 
gen coeff_population_`j' = prop_60_more * (date == `i') 
gen coeff_meteo_`j' = meteo_elec * (date == `i')
	gen coeff_high_meteo_`j' = meteo_elec * (date == `i') * `treat'
}

* Exclude reference period
local exclude = 1
foreach var of varlist coeff_*_`exclude'{
drop `var'
}

* Estimate model
********************************************************************************
									
reghdfe diff_`end_time'_log3_pc  coeff_turnout_* coeff_high* coeff_date_* coeff_density_* coeff_population_* coeff_meteo_* [aw = weight], ///
									absorb(id_dep) vce(cluster date id_de)

* Save and plot estimated coefficients
********************************************************************************

gen b_turnout_high = 0
gen b_turnout_high_ci_up = 0
gen b_turnout_high_ci_low = 0
gen b_turnout_all = 0
gen b_turnout_all_ci_up = 0
gen b_turnout_all_ci_low = 0

forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(06,13,2020)
if `j' != `exclude'{
replace b_turnout_all = _b[coeff_turnout_`j'] if date == `i'
replace b_turnout_all_ci_up = _b[coeff_turnout_`j'] + invttail(e(df_r),0.025)*_se[coeff_turnout_`j'] if date == `i'
replace b_turnout_all_ci_lo = _b[coeff_turnout_`j'] - invttail(e(df_r),0.025)*_se[coeff_turnout_`j'] if date == `i'

replace b_turnout_high = _b[coeff_high_turnout_`j'] if date == `i'
replace b_turnout_high_ci_up = _b[coeff_high_turnout_`j'] + invttail(e(df_r),0.025)*_se[coeff_high_turnout_`j'] if date == `i'
replace b_turnout_high_ci_lo = _b[coeff_high_turnout_`j'] - invttail(e(df_r),0.025)*_se[coeff_high_turnout_`j'] if date == `i'
}
}

gen date1 = date + 0.2
gen date2 = date - 0.2

local calib_pos = mdy(07,`end_time',2020)

local x1 = mdy(06,14,2020)
local x2 = mdy(06,21,2020)
local elec = mdy(06,28,2020)
local x3 = mdy(07,06,2020)
local x4 = mdy(07,13,2020)

local y_calib = 100
local y_elec = 100

twoway rcap b_turnout_all_ci_up b_turnout_all_ci_low date2, color(gs9) lpattern(dash) msize(0) || ///
		scatter b_turnout_all date2, mcolor(gs9) msymbol(Oh) || ///
		rcap b_turnout_high_ci_up b_turnout_high_ci_low date1, color(gs7) lpattern(solid)  msize(0) || ///
		scatter b_turnout_high date1, mcolor(gs7) msymbol(O)  || ///
		, ///
		yline(0, lcolor(black)) xline(`calib_pos', lcolor(black) lpattern(dash)) xline(`elec', lcolor(black) lpattern(dash)) /// 
		xlabel(`x1' `x2' `elec' `x3' `x4', format(%td_Month_dd) labsize(medsmall)) ///
		text(`y_calib' `calib_pos' "End of calibration" "period (July `end_time')", placement(w) justification(right) margin(vsmall) size(small)) ///
		text(`y_elec' `elec' "Second round of" "municipal elections", placement(w) justification(right) margin(vsmall) size(small)) ///
		ylabel(-50(50)150, angle(0) nogrid) ///
		xtitle("") ///
		bgcolor(white) graphregion(color(white)) ///
		legend(col(1) order (2 "Estimate for départements with {bf:low} COVID-19 circulation by June 28" 4 "Estimate for départements with {bf:high} COVID-19 circulation by June 28") region(lstyle(none) lcolor(white))) ///
		ytitle("Effect of turnout on excess cumulated" "hospitalizations per 100,000 inhabitants" " ")
		
		graph export "$graphs/effect_het_`treat'_july`end_time'.eps", replace
		graph export "$graphs/effect_het_`treat'_july`end_time'.png", replace
