clear
set more off

********************************************************************************
* March 26 estimations
********************************************************************************

local end_time = 26

use "$temp/final", clear

* Create groups of departements as a function of epidemicson March 15 (+ 11 days lag in measurement)
*******************************************************************************

* Bottom 25% and 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(25,  `third') 
	gen temp_high25 = (nbre_hospit_corona_cum_pc > r(r1)) if date == mdy(03,`end_time',2020)
	gen temp_high33 = (nbre_hospit_corona_cum_pc > r(r2)) if date == mdy(03,`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(03,`end_time',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 - prediction_`end_time'_log3)/population *100000

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

* 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(03,`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

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

local treat = "high33"

* Create interaction terms
forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(02,29,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'

local date_text : di %td_Month_dd `i'
label var coeff_turnout_`j'  "`date_text' $\times$ Turnout $\times$ Low circulation"
label var coeff_high_turnout_`j' "`date_text' $\times$ Turnout $\times$  High circulation"

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

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

* Full sample
********************************************************************************

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

keep if e(sample)

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

gen b_turnout_high = .
gen b_turnout_all = .

levelsof id_dep, local(levels) 
 foreach id of local levels {	
 
reghdfe diff_`end_time'_log3_pc  coeff_turnout_* coeff_high* coeff_date_*  coeff_density_* coeff_population_* coeff_meteo_* coeff_mobility_* [aw = weight] if id_dep != `id', ///
									absorb(id_dep) vce(cluster date id_de)

forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(02,29,2020)
if `j' != `exclude'{
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)

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) ///
		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
		