clear
set more off

use "$temp/predictions_final", clear

local end_time = 26

distinct id_dep if prediction_26_log3 != .

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

* Bottom third
	local third = 1/3*100
	_pctile nbre_hospit_corona_cum_pc if date == mdy(03,`end_time',2020) & prediction_`end_time'_log3!=.  , p(`third') 
	gen temp_high33 = (nbre_hospit_corona_cum_pc > r(r1)) if date == mdy(03,`end_time',2020)
	bys id_dep: egen high33 = mean(temp_high33)
	drop temp*

* Define sample period and dependent variable
********************************************************************************

gen diff_`end_time'_log3_pc = (nbre_hospit_corona_cum - prediction_`end_time'_log3)/population *100000

local date_start = mdy(03,01,2020)
local date_end = mdy(03,`end_time',2020) + 7
	
keep if inrange(date,`date_start',`date_end')


* Estimations
********************************************************************************

local treat = "high33"

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

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

if `i' == `date_start'{
replace coeff_turnout_`j' = 0
replace coeff_high_`j' = 0
replace coeff_high_turnout_`j' = 0
}
}

gen b_turnout_high = .
gen b_turnout_all = .

* Exclude departements one-by-one
********************************************************************************

levelsof id_dep, local(levels) 
 foreach id of local levels {	
 
reghdfe diff_`end_time'_log3_pc  coeff_turnout_* coeff_high* if id_dep != `id', ///
									a( id_dep date  date#c.density  date#c.more60_prop) vce(cluster id_dep date)

forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(02,29,2020)

replace b_turnout_all = _b[coeff_turnout_`j'] if date == `i' & id_dep == `id'

replace b_turnout_high = _b[coeff_high_turnout_`j'] if date == `i' & id_dep == `id'
}
}

* Graph
********************************************************************************
sort id_dep date

local calib_pos = mdy(03,`end_time',2020) + 0.5

local x1 = mdy(03,01,2020)
local x2 = mdy(03,08,2020)
local elec = mdy(03,15,2020)
local x3 = mdy(03,22,2020)
local x4 = mdy(03,29,2020)

local y_calib = 100
local y_elec = 100

twoway ///
scatter b_turnout_all date, connect(L) msymbol(none) lcolor(gs12) || ///
scatter b_turnout_high date, connect(L) msymbol(none) lcolor(gs4) || ///
, ///
		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 (March `end_time')", placement(w) justification(right) margin(vsmall) size(small)) ///
		text(`y_elec' `elec' "First 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(off) ///
		legend(col(1) order (1 "Effect for all départements" 2 "Additional effect for départements with advanced COVID-19" "epidemics by March 15") region(lstyle(none) lcolor(white))) ///
		ytitle("Effect of turnout on excess cumulated" "hospitalizations per 100,000 inhabitants" " ")

		graph export "$graphs/effect_het_`treat'_exclude_departements_one_by_one.eps", replace
		graph export "$graphs/effect_het_`treat'_exclude_departements_one_by_one.png", replace
		