cd "C:\Users\hpv239\Dropbox\Coping\Data2021\Replication_files"

*Figure 1, panel a

use weekly_countries, clear

sort date2

label var prayer "Internet search intensity for prayer"

line prayer date2 if country=="Worldwide", xla(20820[60]22276, format(%td) angle(45) labsize(vsmall)) xtitle("") xline(20967 21317 21674 22028, lcolor(erose))

*Figure 1, panel b

reg gprayer ib36.month2 ramadan lramadan easter leaster ash lash pent lpent if country=="Worldwide", r
*coefplot, drop(_cons ramadan lramadan easter leaster ash lash pent lpent) keep(*:) omitted baselevels vertical yline(0) coeflabels(13.month2="Jan 17-18" 14.month2="Feb 17-18" 15.month2="Mar 17-18" 16.month2="Apr 17-18" 17.month2="May 17-18" 18.month2="Jun 17-18" 19.month2="Jul 17-18" 20.month2="Aug 17-18" 21.month2="Sep 17-18" 22.month2="Oct 17-18" 23.month2="Nov 17-18" 24.month2="Dec 17-18" 25.month2="Jan 18-19" 26.month2="Feb 18-19" 27.month2="Mar 18-19" 28.month2="Apr 18-19" 29.month2="May 18-19" 30.month2="Jun 18-19" 31.month2="Jul 18-19" 32.month2="Aug 18-19" 33.month2="Sep 18-19" 34.month2="Oct 18-19" 35.month2="Nov 18-19" 36.month2="Dec 18-19" 37.month2="Jan 19-20" 38.month2="Feb 19-20" 39.month2="Mar 19-20" 40.month2="Apr 19-20" 41.month2="May 19-20" 42.month2="Jun 19-20" 43.month2="Jul 19-20" 44.month2="Aug 19-20" 45.month2="Sep 19-20" 46.month2="Oct 19-20" 47.month2="Nov 19-20" 48.month2="Dec 19-20", labsize(vsmall) angle(45)) ylabel(-0.2(0.1)0.6, labsize(vsmall)) ytitle("Avg prayer search growth relative to the year before",size(small))
*graph export "fig1.png", as(png) replace

*Table 1
use weekly_countries, clear

gen p5=0
replace p5=1 if year==2020 & ( month==1 | month==2)
replace p5=2 if year==2020 & ( month==3 | month==4)
replace p5=3 if year==2020 & month>=5 & month!=.

gen p6=0
replace p6=1 if year==2019 & ( month==1 | month==2)
replace p6=2 if year==2019 & ( month==3 | month==4)
replace p6=3 if year==2019 & month>=5 & month!=.
replace p6=. if year==2020

label define cat2 1 "Jan-Feb 2020" 2 "Mar-Apr 2020" 3 "May-Dec 2020" 
label define cat4 1 "Jan-Feb 2019" 2 "Mar-Apr 2019" 3 "May-Dec 2019" 

label values p5 cat2
label values p6 cat4

est clear

local pray gprayer

reghdfe `pray' i.p5 if country=="Worldwide" & year>2018, absorb(ramadan lramadan easter leaster ash lash pent lpent) 
estadd loc samp "2018-2020"
estadd loc sample "World"
estadd loc trend "N"
estadd loc ctrend "N"
estadd ysumm
eststo r1

local pray gprayer lprayer

reghdfe `pray' i.p5 if country=="Worldwide" & year>2018, absorb(ramadan lramadan easter leaster ash lash pent lpent)
estadd loc samp "2018-2020"
estadd loc sample "World"
estadd loc trend "N"
estadd loc ctrend "N"
estadd ysumm
eststo r2

reghdfe `pray' i.p5 date2 if country=="Worldwide" & year>2018, absorb(ramadan lramadan easter leaster ash lash pent lpent)
estadd loc samp "2018-2020"
estadd loc sample "World"
estadd loc trend "Y"
estadd loc ctrend "N"
estadd ysumm
eststo r3

reghdfe `pray' i.p6 date2 if country=="Worldwide", absorb(ramadan lramadan easter leaster ash lash pent lpent) 
estadd loc samp "2017-2019"
estadd loc sample "World"
estadd loc trend "Y"
estadd loc ctrend "N"
estadd ysumm
eststo r4

reghdfe `pray' i.p5 date2 if country!="Worldwide" & year>2018, absorb(ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc samp "2018-2020"
estadd loc sample "Country"
estadd loc trend "Y"
estadd loc ctrend "N"
estadd ysumm
eststo r5

reghdfe `pray' i.p5 date2 if country!="Worldwide" & year>2018, absorb(code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc samp "2018-2020"
estadd loc sample "Country"
estadd loc trend "Y"
estadd loc ctrend "Y"
estadd ysumm
eststo r6

reghdfe `pray' i.p6 date2 if country!="Worldwide", absorb(code#c.date2 ramadan lramadan easter leaster ash lash pent lpent)  cl(code)
estadd loc samp "2017-2019"
estadd loc sample "Country"
estadd loc trend "Y"
estadd loc ctrend "Y"
estadd ysumm
eststo r7

esttab r1 r2 r3 r4 r5 r6 r7 
*using "week_regs_n1.tex", replace fragment drop(_cons 0.p5 0.p6 date2) order(lprayer) starlevels(* .10 ** .05 *** .01) se(3) b(a2) label stats(r2 N samp sample trend ctrend ymean, labels(R-squared Observations "Period" "Panel" "Trend" "Country-trend" "MeanDepVar")) nogaps nonotes  nocons nonum mlabels("(1)" "(2)" "(3)" "(4)" "(5)" "(6)" "(7)")

*Figure 2
use map_data, clear

spmap gprayer using cntr_coord if neg==0, id(_ID) fcolor(Greens2) clnumber(8) ndfcolor(gs15) legend(size(medium)) polygon(data("cntr_coord2.dta") by(neg) fcolor(gold) legenda(on)) 
*graph export "prayer_map_growth.png", as(png) replace


*Table 2
use weekly_countries, clear

gen p5=0
replace p5=1 if year==2020 & ( month==1 | month==2)
replace p5=2 if year==2020 & ( month==3 | month==4)
replace p5=3 if year==2020 & ( month==5 | month==6)
replace p5=4 if year==2020 & ( month==7 | month==8)
replace p5=5 if year==2020 & ( month==9 | month==10)
replace p5=6 if year==2020 & ( month==11 | month==12)

label define cat3 1 "Jan-Feb 2020" 2 "Mar-Apr 2020" 3 "May-Jun 2020" 4 "Jul-Aug 2020" 5 "Sep-Oct 2020" 6 "Nov-Dec 2020" 

label values p5 cat3

est clear

local pray gprayer lprayer

local if & year>2018

local period p5

reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & year>2018, absorb(code#c.date2  ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "N"
estadd loc trend "Y"
estadd loc cont "All"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r1

reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & year>2018, absorb(code code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "Y"
estadd loc trend "Y"
estadd loc cont "All"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r3

reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & (con==1 | con==4) & year>2018, absorb(code code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "Y"
estadd loc trend "Y"
estadd loc cont "Asia"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r5

reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & con==3 & year>2018, absorb(code code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "Y"
estadd loc trend "Y"
estadd loc cont "Africa"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r6

quietly reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & con==5 & year>2018, absorb(code code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "Y"
estadd loc trend "Y"
estadd loc cont "Americas"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r7

quietly reghdfe gprayer lprayer i.p5 date2 if country!="Worldwide" & con==2 & year>2018, absorb(code code#c.date2 ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc cntr "Y"
estadd loc trend "Y"
estadd loc cont "Europe"
sum prayer2019 if e(sample)
estadd scalar mean=r(mean)
eststo r8

esttab r1 r3 r5 r6 r7 r8 
*using "week_regs_more.tex", replace fragment drop(_cons 0.p5 date2 lprayer) order(lprayer) starlevels(* .10 ** .05 *** .01) se(3) b(a2) label stats(r2 N N_clust cntr cont, labels(R-squared Observations "Countries" "Country FE" "Region")) nogaps mlabels("(1)" "(2)" "(3)" "(4)" "(5)" "(6)") nocons nonotes nonum

*Table 3
use weekly_countries, clear

quietly reghdfe gprayer lprayer time_pan  if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg0
quietly reghdfe gprayer lprayer casedum10 if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd ysumm
eststo reg1
quietly reghdfe gprayer lprayer deathdum if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg2
quietly reghdfe gprayer lprayer anosmiadum if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg3
quietly reghdfe gprayer lprayer casedum10 time_pan if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg4
quietly reghdfe gprayer lprayer deathdum time_pan if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg5
quietly reghdfe gprayer lprayer anosmiadum time_pan if country!="Worldwide" & year==2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
estadd loc lag "Y"
estadd ysumm
eststo reg6
quietly reghdfe gprayer lprayer casedum10 if country!="Worldwide"  & year==2020, absorb(date2 code#c.date2 code) cl(code)
estadd loc timefe "Y"
estadd loc lag "Y"
estadd ysumm
eststo reg7
quietly reghdfe gprayer lprayer deathdum if country!="Worldwide" & year==2020, absorb(date2 code#c.date2 code) cl(code)
estadd loc timefe "Y"
estadd loc lag "Y"
estadd ysumm
eststo reg8
quietly reghdfe gprayer lprayer anosmiadum if country!="Worldwide" & year==2020, absorb(date2 code#c.date2 code) cl(code)
estadd loc timefe "Y"
estadd loc lag "Y"
estadd ysumm
eststo reg9

esttab reg0 reg1 reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9 
*using "week_regs.tex", replace fragment drop( _cons lprayer) starlevels(* .10 ** .05 *** .01) se(3) b(a2) label stats(r2 N N_clust timefe, labels(R-squared Observations Countries "Time Fe")) nogaps  nonum mlabels("(1)" "(2)" "(3)" "(4)" "(5)" "(6)" "(7)" "(8)" "(9)" "(10)") nonotes  nocons

*Table 4
use weekly_countries_bans, clear

reghdfe gprayer lprayer mass_all_avg if country!="Worldwide", absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
eststo r1
reghdfe gprayer lprayer mass_all_avg time_pan if country!="Worldwide", absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
eststo r2
reghdfe gchurch lchurch mass_all_avg if country!="Worldwide", absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
eststo r3
reghdfe gmosque lmosque mass_all_avg if country!="Worldwide", absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd loc timefe "N"
eststo r4
esttab r1 r2 r3 r4 
*using "week_bans.tex", replace fragment drop( _cons lprayer lchurch lmosque) starlevels(* .10 ** .05 *** .01) se(3) b(a2) label stats(r2 N N_clust, labels(R-squared Observations Countries)) nogaps  nonum mlabels("(1)" "(2)" "(3)" "(4)" "(5)") nonotes  nocons

*Figure 3, panel a
use internetchurch, clear
tsset date
line internetchurch date, lw(medium medium medium thick) xscale(range(21915[10]22010)) xlabel(21915[10]22010, format(%tdD_m) angle(45) labsize(vsmall)) xtitle("") ytitle("") xtitle("") legend(on)
*graph export "internetchurch.png", as(png) replace

*Figure 3, panel b
use internetchurch_prayer, clear
tsset date
line internetchurch prayer date, lw(medium medium medium thick) xscale(range(21915[10]22010)) xlabel(21915[10]22010, format(%tdD_m) angle(45) labsize(vsmall)) xtitle("") legend(order(2 "Prayer" 1 "Internet church") size(vsmall))

* Table 5
use daily_countries, clear
quietly: reghdfe prayer time_pan if NAMES_STD!="worldwide", absorb(pent ash code#c.date2 code ramadan easter) cl(code)
estadd ysumm
estadd loc cntrfe "Yes"
eststo reg00

foreach num in 0 1 2 3 4 5 6 {
use daily_countries, clear
keep if dow==`num'
quietly: reghdfe prayer time_pan if NAMES_STD!="worldwide", absorb(pent ash code#c.date2 code ramadan easter) cl(code)
estadd ysumm
estadd loc cntrfe "Yes"
test time_pan=10.6793
estadd scalar rp=r(p)
eststo reg`num'
}
esttab reg00 reg0 reg1 reg2 reg3 reg4 reg5 reg6 
*using "main.tex", replace fragment drop(_cons) starlevels(* .10 ** .05 *** .01) se(3) b(a2) label stats(r2 N N_clust, labels(R-squared Observations Countries)) nogaps mlabels("All" "Sun" "Mon" "Tue" "Wed" "Thu" "Fri" "Sat", lhs("Sample:")) nonotes  nocons

*Figure 5, panel a
use weekly_countries, clear

xtile prayer2019_perc=prayer2019_low, nq(4)

replace prayer2019_low=prayer2019_perc

foreach num in 0 1 2 3 {
gen pray`num'=(prayer2019_low>`num' & prayer2019_low<=`num'+1) & !missing(prayer2019_low)
gen prayl`num'=(prayer2019_low<=`num') & !missing(prayer2019_low)
gen prayh`num'=(prayer2019_low>`num'+1) & !missing(prayer2019_low)
}

gen pray19=.
gen prayl19=.
gen prayh19=.
gen time_int=.
gen time_intl=. 
gen time_inth=. 
tempname myresults_ras3
postfile `myresults_ras3' str12 mean_num nr mean_df_short mean_coeff_short mean_se_short using myresults_ras3.dta, replace
foreach j in 0 1 2 3 {
replace pray19=pray`j'
replace time_int=pray`j'*time_pan
replace prayl19=pray`j'
replace time_intl=prayl`j'*time_pan
replace prayh19=prayh`j'
replace time_inth=prayh`j'*time_pan
reghdfe gprayer lprayer time_int time_intl time_inth if country!="Worldwide" & year>=2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
post `myresults_ras3' (`"pray`j'"') (`j') (e(df_r)) (`=_b[time_int]') (`=_se[time_int]')
}

postclose `myresults_ras3' 
use myresults_ras3 , clear
import excel labels1e.xlsx, sheet("Sheet1") firstrow clear
rename number nr
merge 1:1 nr using myresults_ras3
drop _merge
gen p_val  = (2 * ttail( mean_df_short , abs( mean_coeff_short / mean_se_short )))
replace p_val=round(p_val, 0.001)
gen stars = 0 if p_val >=0.1
replace stars = 1 if p_val <0.1
replace stars = 2 if p_val <0.05
replace stars = 3 if p_val <0.01
sort mean_coeff_short
gen laws_order = _n
label var laws_order "The rise in prayer for different initial prayer levels"
save myresults_ras3, replace

set scheme plotplain

gen lowerbound_95 =  mean_coeff_short - invttail(mean_df_short,0.025)*mean_se_short
gen upperbound_95 =  mean_coeff_short + invttail(mean_df_short,0.025)*mean_se_short
gsort -nr
sencode name, g(nam2) 

eclplot mean_coeff_short lowerbound_95 upperbound_95 nam2,  /// 
horizontal xline(0) aspectratio(0.5) estopts(msize(vsmall) msymbol(circle) ) ciopts(msize(vsmall) lwidth(thin))  /// 
ytitle("Prayer search shares in 2019",size(small)) xtitle("The rise in prayer search shares after March 11 2020",size(small)) ylabel(1(1)4, labsize(vsmall)) xscale(range(-0.2(0.1)0.4)) xlabel(-0.2(0.1)0.4, labsize(vsmall))
*graph export "impact_prayerlevels2.png", as(png) replace

*Figure 5, panel b
use weekly_countries, clear

drop if dist_eq_34>1500000 

xtile prayer2019_perc=risk_earthq, nq(4)

replace prayer2019_low=prayer2019_perc

foreach num in 0 1 2 3 {
gen pray`num'=(prayer2019_low>`num' & prayer2019_low<=`num'+1) & !missing(prayer2019_low)
gen prayl`num'=(prayer2019_low<=`num') & !missing(prayer2019_low)
gen prayh`num'=(prayer2019_low>`num'+1) & !missing(prayer2019_low)
}

gen pray19=.
gen prayl19=.
gen prayh19=.
gen time_int=.
gen time_intl=. 
gen time_inth=. 
tempname myresults_ras3
postfile `myresults_ras3' str12 mean_num nr mean_df_short mean_coeff_short mean_se_short using myresults_ras3.dta, replace
foreach j in 0 1 2 3 {
replace pray19=pray`j'
replace time_int=pray`j'*time_pan
replace prayl19=pray`j'
replace time_intl=prayl`j'*time_pan
replace prayh19=prayh`j'
replace time_inth=prayh`j'*time_pan
reghdfe gprayer lprayer time_int time_intl time_inth if country!="Worldwide" & year>=2020 & dist_eq_34<=1500000, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
post `myresults_ras3' (`"pray`j'"') (`j') (e(df_r)) (`=_b[time_int]') (`=_se[time_int]')
}

postclose `myresults_ras3' 
use myresults_ras3 , clear
import excel labels1e.xlsx, sheet("Sheet1") firstrow clear
rename number nr
merge 1:1 nr using myresults_ras3
drop _merge
gen p_val  = (2 * ttail( mean_df_short , abs( mean_coeff_short / mean_se_short )))
replace p_val=round(p_val, 0.001)
gen stars = 0 if p_val >=0.1
replace stars = 1 if p_val <0.1
replace stars = 2 if p_val <0.05
replace stars = 3 if p_val <0.01
sort mean_coeff_short
gen laws_order = _n
label var laws_order "The rise in prayer for different initial prayer levels"
save myresults_ras3, replace

set scheme plotplain

gen lowerbound_95 =  mean_coeff_short - invttail(mean_df_short,0.025)*mean_se_short
gen upperbound_95 =  mean_coeff_short + invttail(mean_df_short,0.025)*mean_se_short
gsort -nr
sencode name, g(nam2) 

eclplot mean_coeff_short lowerbound_95 upperbound_95 nam2,  /// 
horizontal xline(0) aspectratio(0.5) estopts(msize(vsmall) msymbol(circle) ) ciopts(msize(vsmall) lwidth(thin))  /// 
ytitle("Earthquake risk",size(small)) xtitle("The rise in prayer search shares after March 11 2020",size(small)) ylabel(1(1)4, labsize(vsmall)) xscale(range(-0.2(0.1)0.6)) xlabel(-0.2(0.1)0.6, labsize(vsmall))
*graph export "impact_prayerlevels3.png", as(png) replace

*Figure 5, panel c
use weekly_countries, clear

gen time_int=.
gen time_intl=. 
tempname myresults_ras3
postfile `myresults_ras3' str12 mean_num mean_df_short mean_coeff_short mean_se_short using myresults_ras3.dta, replace
foreach j in chstian cath00 prot00 muslims hindus buddhist{
replace time_int=denom`j'*time_pan
replace time_intl=notdenom`j'*time_pan
quietly reghdfe gprayer lprayer time_int time_intl if country!="Worldwide" & year>=2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
post `myresults_ras3' (`"denom`j'"') (e(df_r)) (`=_b[time_int]') (`=_se[time_int]')
}

postclose `myresults_ras3' 
use myresults_ras3 , clear
gen num = _n
save myresults_ras3, replace
import excel labels2.xlsx, sheet("Sheet1") firstrow clear
merge 1:1 mean_num using myresults_ras3
drop _merge
gen p_val  = (2 * ttail( mean_df_short , abs( mean_coeff_short / mean_se_short )))
replace p_val=round(p_val, 0.001)
gen stars = 0 if p_val >=0.1
replace stars = 1 if p_val <0.1
replace stars = 2 if p_val <0.05
replace stars = 3 if p_val <0.01
sort mean_coeff_short
gen laws_order = _n
label var laws_order "The rise in prayer for different initial prayer levels"
save myresults_ras3, replace

set scheme plotplain

gen lowerbound_95 =  mean_coeff_short - invttail(mean_df_short,0.025)*mean_se_short
gen upperbound_95 =  mean_coeff_short + invttail(mean_df_short,0.025)*mean_se_short
sort num
sencode name, g(nam2) 

eclplot mean_coeff_short lowerbound_95 upperbound_95 nam2,  /// 
horizontal xline(0) aspectratio(0.5) estopts(msize(vsmall) msymbol(circle) ) ciopts(msize(vsmall) lwidth(thin))  /// 
ytitle("Denominations",size(small)) xtitle("The rise in prayer search intensity after March 11 2020",size(small)) ylabel(1(1)6, labsize(vsmall)) xlabel(-0.4(0.1)0.6, labsize(vsmall))
*graph export "impact_denoms.png", as(png) replace

*Table 6
use weekly_countries, clear

gen time_int=.
gen time_int_rel=.
gen time_rel=.
label var time_int "Pandemic x Variable"
label var time_int_rel "Pandemic x Variable x Prayer 2019"
label var time_rel "Pandemic x Prayer 2019"

label var time_int "Case dummy x Variable"
label var time_int_rel "Case dummy x Variable x Prayer 2019"
label var time_rel "Case dummy x Prayer 2019"

eststo clear

foreach var in wdi_povgap190 bl_lhmf gini ffp_ued ffp_fsi who_mrt hf_prights{
replace `var'=`var'/10
}

/*Panel A.*/
foreach var in lnrgdp wdi_povgap190 bl_lhmf gini ffp_ued ffp_fsi who_mrt hf_prights{
gen `var'_pan=`var'*time_pan
replace time_int=`var'_pan
quietly reghdfe gprayer lprayer time_pan time_int if country!="Worldwide" & year>=2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd ysumm
quietly sum `var' if e(sample), det
estadd loc cntrfe "Yes"
eststo r`var'
}

esttab rlnrgdp rwdi_povgap190 rbl_lhmf rgini rffp_ued rffp_fsi rwho_mrt rhf_prights 
*using "who.tex", fragment replace starlevels(* .10 ** .05 *** .01) se(3) b(a2) label drop(lprayer _cons) stats(r2 N N_clust , labels(R-squared Observations Countries)) mlabels("GDP" "Poverty" "Educ" "Gini"  "Uneven" "Fragile" "AdultMort" "PropRights", lhs("\textbf{Panel A}"))  nogaps 

/*Panel B.*/
eststo clear
foreach var in  lnrgdp wdi_povgap190 bl_lhmf gini ffp_ued ffp_fsi who_mrt hf_prights{
replace `var'_pan=`var'*time_pan
replace time_int=`var'_pan
replace time_rel=time_pan*prayer2019_low
quietly reghdfe gprayer lprayer time_pan time_int time_rel if country!="Worldwide" & year>=2020, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd ysumm
estadd loc cntrfe "Yes"
eststo r`var'

}
esttab rlnrgdp rwdi_povgap190 rbl_lhmf rgini rffp_ued rffp_fsi rwho_mrt rhf_prights 
*using "who.tex", fragment append starlevels(* .10 ** .05 *** .01) se(3) b(a2) label drop(lprayer _cons) stats(r2 N N_clust , labels(R-squared Observations Countries)) mlabels("" "" "" ""  "" "" "" "", lhs("\textbf{Panel B}"))  nogaps nonum

/*Panel C.*/
foreach var in lnrgdp wdi_povgap190 bl_lhmf gini ffp_ued ffp_fsi who_mrt hf_prights{
replace time_int=`var'_pan
quietly reghdfe gprayer lprayer time_pan time_int if country!="Worldwide" & year>=2020 & dist_eq_34<=1500000, absorb(code#c.date2 code ramadan lramadan easter leaster ash lash pent lpent) cl(code)
estadd ysumm
quietly sum `var' if e(sample), det
estadd loc cntrfe "Yes"
eststo r`var'
}
esttab rlnrgdp rwdi_povgap190 rbl_lhmf rgini rffp_ued rffp_fsi rwho_mrt rhf_prights 
*using "who.tex", fragment append starlevels(* .10 ** .05 *** .01) se(3) b(a2) label drop(lprayer _cons) stats(r2 N N_clust , labels(R-squared Observations Countries)) mlabels("" "" "" ""  "" "" "" "", lhs("\textbf{Panel C}"))  nogaps nonum

/*Panel D. IV*/
gen time_eq=time_pan*risk_earthq

eststo clear
foreach var in  lnrgdp wdi_povgap190 bl_lhmf gini ffp_ued ffp_fsi who_mrt hf_prights{
replace `var'_pan=`var'*time_pan
replace time_int=`var'_pan
replace time_rel=time_pan*prayer2019_low
replace time_eq=time_pan*risk_earthq
quietly ivreg2 gprayer lprayer time_pan time_int ( time_rel = time_eq ) i.code#c.date2 i.code ramadan lramadan easter leaster ash lash pent lpent if country!="Worldwide"  & year>=2020 & dist_eq_34<=1500000, r cl(code) partial(i.code#c.date2 i.code ramadan lramadan easter leaster ash lash pent lpent)
estadd ysumm
estadd loc cntrfe "Yes"
eststo r`var'
}
esttab rlnrgdp rwdi_povgap190 rbl_lhmf rgini rffp_ued rffp_fsi rwho_mrt rhf_prights 
*using "who.tex", fragment append starlevels(* .10 ** .05 *** .01) se(3) b(a2) label keep(time_int time_rel time_pan) order(time_pan time_int time_rel) stats(r2 N N_clust widstat, labels(R-squared Observations Countries FirstStageF)) mlabels("" "" "" ""  "" "" "" "", lhs("\textbf{Panel D}"))  nogaps nonum





