clear
set more off

********************************************************************************
* March 25 and March 27 estimations
********************************************************************************

foreach end_time in 25 27{

use "$temp/final", clear

* Add March 25 predictions
********************************************************************************
merge 1:1 date departement_code using "$temp/predictions_`end_time'", update
drop _merge

* 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(03,`end_time',2020) & insample == 0
local a = r(sum)
summ population if date == mdy(03,`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 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"

	preserve
	
* 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 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_* coeff_mobility_* [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(02,29,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_low = _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_low = _b[coeff_high_turnout_`j'] - invttail(e(df_r),0.025)*_se[coeff_high_turnout_`j'] if date == `i'
}
}

foreach var of varlist *_up{
replace `var' =150 if `var' > 150 & `var' != .
}

foreach var of varlist *_low{
replace `var' =-50 if `var' < -50 & `var' != .
}

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

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 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 (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(col(1) order (2 "Estimate for départements with {bf:low} COVID-19 circulation by March 15" 4 "Estimate for départements with {bf:high} COVID-19 circulation 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'_march`end_time'.eps", replace
		graph export "$graphs/effect_het_`treat'_march`end_time'.png", replace

* Calculate total elections-related excess hospitalizations	
********************************************************************************

keep if e(sample)
			
gen excess = .
gen excess_up = .
gen excess_lo = .

* Get average turnout_high coefficient over calibration period
summ b_turnout_high if date <= mdy(03,`end_time',2020)
local bench = r(mean)

forvalues i =`date_start'(1)`date_end'{
local j = `i' - mdy(02,29,2020)
if `j' != `exclude'{
replace excess = turnout_2020 * (_b[coeff_high_turnout_`j'] - (`bench') )  if date == `i' & `treat' == 1
replace excess_up = turnout_2020 * (_b[coeff_high_turnout_`j']  + invttail(e(df_r),0.025)*_se[coeff_high_turnout_`j'] - (`bench') )  if date == `i' & `treat' == 1
replace excess_lo = turnout_2020 * (_b[coeff_high_turnout_`j']  - invttail(e(df_r),0.025)*_se[coeff_high_turnout_`j'] - (`bench') )  if date == `i' & `treat' == 1
}
}

foreach var of varlist excess*{
replace `var' = 0 if `var' == .
gen `var'_tot = `var' * population / 100000
}

collapse (sum) excess*_tot nbre_hospit_corona_cum, by(date)

* Truncate for graphical representation
replace excess_lo_tot = 0 if excess_lo_tot < 0

local calib_pos = mdy(03,`end_time',2020)
local y_calib = 35000

local y_elec = 35000

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 pos1 = `date_end'
summ nbre_hospit_corona_cum if date == `pos1'
local y_pos1 =r(mean)

local pos2 = `date_end'
summ excess_tot if date > mdy(03,`end_time',2020)
local y_pos2 =r(max)

twoway scatter nbre_hospit_corona_cum date, msymbol(none) connect(line) lcolor(gs9) lwidth(thick) || ///
		scatter excess_tot date if date >= mdy(03,`end_time',2020), msymbol(none) connect(line) lcolor(gs6) lwidth(thick) lpattern(dash) || ///
		scatter excess_up_tot date if date >= mdy(03,`end_time',2020), msymbol(none) connect(line) lcolor(gs6) lwidth(medthin ) lpattern(shortdash) || ///
		scatter excess_lo_tot date if date >= mdy(03,`end_time',2020), msymbol(none) connect(line) lcolor(gs6) lwidth(medthin ) lpattern(shortdash) || ///
		, ///
		xline(`calib_pos', lcolor(black) lpattern(dash)) xline(`elec', lcolor(black) lpattern(dash)) /// 
		ylabel(0(10000)40000 0 "0", angle(0) nogrid format(%6.0fc)) ///
		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)) ///
		text(`y_pos1' `pos1' "Total cumulated" "hospitalizations", placement(nw) justification(left) margin(vsmall) size(small)) ///
		text(`y_pos2' `pos2' "Elections-related" "excess" "hospitalizations", box fcolor(white) lcolor(white) placement(nw) justification(right) margin(vsmall) size(small)) ///
		xtitle("") ///
		bgcolor(white) graphregion(color(white)) plotregion(margin(1 0 0 1)) ///
		legend(off)
		
		graph export "$graphs/effect_het_`treat'_march`end_time'_excess.eps", replace
		graph export "$graphs/effect_het_`treat'_march`end_time'_excess.png", replace
		
	restore
}
