********************************************************************************
*
* Figure 5: Difference-in-differences evidence: Falsification in no-speech states
*
********************************************************************************
 
global program "Figure5"

cap log close
cap log using "$log_files/log_$program$today", replace

clear all
set mem 2g
set more off
set matsize 11000
set scheme s1color

***************************
***** 1. PREPARE DATA *****
***************************

// Get data

use "$data/Data_County", clear

// Controls definition

global geog 		"railroad_dens ldist_capital ldist_largest" // 
global demo 		"urban_dum perc_pop_dstr lpop_total perc_white perc_male21 perc_below21 perc_native perc_literate newsp_circ_elig" 
global demo_cont	"perc_pop_dstr lpop_total perc_white perc_male21 perc_below21 perc_native perc_literate newsp_circ_elig" 
global econ 		"perc_farm perc_manuf dlog_avbasket ldist_silver ldist_gold" 
global econ_cont	"perc_farm perc_manuf dlog_avbasket ldist_silver ldist_gold" 
global econ_p 		"perc_farm perc_manuf log_avbasket ldist_silver ldist_gold" 
global vote 		"lag_margin_pres_ttl" 			
global demo_add		"urban_dum lpop_total larea perc_white perc_below21 perc_native perc_literate newsp_circ_elig pop_male2144 pop_male4564 pop_male6599 perc_catholic perc_black"

// New variables

sort county year
g id=_n
g nospeech=1*(speech==0 & state_speech==1)

tab state_speech, miss
recode state_speech (0=1) (1=0), g(nostate_speech)
tab nostate_speech, miss

*********************************
***** 2. IDENTIFY NEIGHBORS *****
*********************************

// Fake treatment

qui tab state if state_speech==0
return list
qui tab state if state_speech==1
return list

sort id

qui teffects nnmatch (dempop_pres_ttl lag_dem_pres_ttl lag_pop_pres_ttl $demo $geog $econ) (nostate_speech) if year==1896 & railroad_dum==1, ///
	metric(euclidean) biasadj(lag_dem_pres_ttl lag_pop_pres_ttl $demo_cont $geog $econ_cont) vce(robust) nneighbor(1) generate(nn)

// Pair id
	
egen pair_id=group(id nn1)

tempfile  tmp
save `tmp.dta', replace

// Get pairs

use `tmp.dta', clear

keep if nostate_speech==1
keep id nn pair_id

drop if id==. | nn1==.

merge m:1 id using `tmp.dta', keepusing(county)
drop if _merge!=3
drop _merge id

rename county county_false

rename nn1 id
merge m:1 id using `tmp.dta', keepusing(county speech)
drop if _merge!=3
drop _merge id

rename county county_true
rename speech speech_false

// Reshape to keep pairs

reshape long county_, i(pair_id) j(sample) string
rename county_ county

tempfile  tmp
save `tmp.dta', replace

// Expand to get county characteristics

use "$data/Data_County", clear

keep year
duplicates drop

cross using `tmp.dta'

// Get county characteristics

merge m:1 county year using "$data/Data_County"
drop if _merge!=3
drop _merge

// State-year fixed effects

egen year_state=group(year state)

// Time varying controls

global X "c.perc_farm#i.year c.perc_manuf#i.year" 

// Treatment-year interaction

cap g x=1
forvalues y=1880(4)1900{
	g speech_`y'=speech_false*(year==`y')
	la var speech_`y' "`y'" 
}

// Set baseline to 1892

la var x "1892"

*******************************
***** 3. DIF-IN-DIF GRAPH *****
*******************************

// False regression

reghdfe dempop_pres_ttl speech_1884 speech_1888 x speech_1896 speech_1900 $X	railroad_dum $geog $demo $econ_p if year>1880 & year<=1900 & sample=="false", cluster(county) a(county year_state)  
estimates store DD_binary
test speech_1884=speech_1888=speech_1896=speech_1900=0

coefplot DD_binary, keep(speech_1884 speech_1888 x speech_1896 speech_1900) levels(95) ciopts(recast(rcap) lc(gs13))  vertical msymbol(S) msize(large)  mfcolor(white) pstyle(p3) base omitted /// 
	mlabel mlabsize(small) format(%9.3f) mlabposition(6) mlabcolor(navy) mlabgap(*2) mcolor(navy) /// 
	ytitle("Coefficient (95% C.I.)") /// 
	ylabel(-4(2)4, angle(0) format(%9.0f)) yscale(r(-4 4)) ///  
	xlabel(1 `" "1884" "(None)" "' 2 `" "1888" "(None)" "' 3 `" "1892" "(None)" "' 4 `" "1896" "(Bryan)" "' 5 `" "1900" "(Both D & R)" "') ///
	xline(3.5, lcolor(gs12) lpattern(dash)) xline(4.5, lcolor(gs12) lpattern(dash)) /// 
	yline(0, lpattern(dash)) yline(2, lpattern(dash) lcolor(gs15)) yline(4, lpattern(dash) lcolor(gs15)) yline(-2, lpattern(dash) lcolor(gs15)) yline(-4, lpattern(dash) lcolor(gs15))  

************************************************
***** 4. DID GRAPH WITH JOINT SIGNIFICANCE *****
************************************************
 
// Year specific coefficient

reghdfe dempop_pres_ttl speech_1884 speech_1888 x speech_1896 speech_1900 $X	railroad_dum $geog $demo $econ_p if year>1880 & year<=1900 & sample=="false", cluster(county) a(county year_state)
matrix b=e(b)
g b1884=b[1,1]
g b1888=b[1,2]
g b1892=b[1,3]
g b1896=b[1,4]
g b1900=b[1,5]
matrix V=e(V)
g se1884=sqrt(V[1,1])
g se1888=sqrt(V[2,2])
g se1892=sqrt(V[3,3])
g se1896=sqrt(V[4,4])
g se1900=sqrt(V[5,5])
g df_full=e(df_r)
forvalues y=1884(4)1900{
	g t`y'=b`y'/se`y'
	g p`y'=2*ttail(df_full, abs(t`y'))
}

// Grouped 1884-1888 and 1896-1900 coefficients

g speech_pre=speech_1884+speech_1888
g speech_post=speech_1896+speech_1900
reghdfe dempop_pres_ttl speech_pre x speech_post $X	railroad_dum $geog $demo $econ_p if year>1880 & year<=1900 & sample=="false", cluster(county) a(county year_state)
matrix b=e(b)
g bpre=b[1,1]
g bpost=b[1,3]
matrix V=e(V)
g sepre=sqrt(V[1,1])
g sepost=sqrt(V[3,3])
g df_group=e(df_r)
foreach y in pre post{
	g t`y'=b`y'/se`y'
	g p`y'=2*ttail(df_group, abs(t`y'))
}

// Keep results

keep b1* bp* se1* sep* df_* t1* tp* p1* pp*
duplicates drop

rename b1* b_full_1*
rename se1* se_full_1*
rename t1* t_full_1*
rename p1* p_full_1*

rename *pre *_pre
rename *post *_post

g id=1
reshape long b_full_ se_full_ t_full_ p_full_, i(id) j(year)

foreach var in b se t p{
	replace `var'_post=. if year<=1892
	replace `var'_pre=. if year>=1892
}

expand 2 if year!=1892
bysort year: egen rank=rank(year), unique

foreach var in b se t p{
	replace `var'_full=. if rank==2
}

replace year=1883 if year==1884 & rank==2
replace year=1889 if year==1888 & rank==2
replace year=1895 if year==1896 & rank==2
replace year=1901 if year==1900 & rank==2

drop rank
sort year

expand 2 if year==1883 | year==1901
bysort year: egen rank=rank(year), unique

replace year=1882 if year==1883 & rank==2
replace year=1902 if year==1901 & rank==2

replace b_post=. if year==1902
replace b_pre=. if year==1882

// Get bounds

foreach case in full pre post{
	gen ub_`case'= b_`case'+invttail(df_full,0.025)*se_`case'
	gen lb_`case'= b_`case'+invttail(df_full,0.975)*se_`case'
}

****************************
***** 3. EXPORT FIGURE *****
****************************

g yline=0

format b_full %9.3f

twoway (rarea ub_pre lb_pre year, fcolor(navy) fintensity(2) lcolor(navy*.5) lwidth(vthin)) ///
	(line b_pre year, lcolor(navy*.5)) ///  
	(line yline year, lcolor(cranberry) lpattern(dash)) ///
	(rcap ub_full lb_full year, color(gs13)) ///
	(scatter b_full year, mcolor(navy) msymbol(S) msize(large) mfcolor(white) mlabel(b_full) mlabsize(small) mlabposition(6) mlabcolor(navy) mlabgap(*2) mcolor(navy)), ///
	ytitle("Coefficient (95% C.I.)") xtitle("") ///
	ylabel(-4(2)4, angle(0) format(%9.0f)) yscale(r(-4.4 4.4)) ///
	xlabel(1884 `" "1884" "(None)" "' 1888 `" "1888" "(None)" "' 1892 `" "1892" "(None)" "' 1896 `" "1896" "(Bryan)" "' 1900 `" "1900" "(Both D & R)" "') ///
	yline(2, lpattern(dash) lcolor(gs15)) yline(4, lpattern(dash) lcolor(gs15)) yline(-2, lpattern(dash) lcolor(gs15)) yline(-4, lpattern(dash) lcolor(gs15)) /// 
	xline(1894, lcolor(gs12) lpattern(dash)) xline(1898, lcolor(gs12) lpattern(dash)) ///
	legend(off) graphregion(color(white))
graph export "$figures/$program.pdf", replace

cap log close