##### R CODE
##### The Link between Regional Temperature and Regional Incomes: Econometric Evidence with Sub-National Data
##### December 2020
##### Christina Greßer, Daniel Meierrieks, David Stadelmann
##### main replication files

library(lfe)
library(stargazer)
library(miceadds)
options(scipen=50)
library(ggplot2)
library(mgcv)
library(coefplot)
library(pastecs)
library(lme4)
library(plotrix)

load("data_climate_Gennaioli.Rda")
load("data_climate_Gennaioli_growth.Rda")
load("data_climate_DHS2005.Rda")
load("data_climate_DHS2015.Rda")





#####FIGURES#####
#####Manuscript - Figure 1#####
plotincome1<-ggplot(data_climate_Gennaioli,aes(x=Temperature,y=Ln_GDP_region))+geom_point()+geom_smooth()
plotincome1 + theme(panel.grid.major = element_blank(), axis.text = element_text(size = 15), axis.title=element_blank(),
                    panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
cor.test(data_climate_Gennaioli$Ln_GDP_region, data_climate_Gennaioli$Temperature)

plotnightlights1<-ggplot(data_climate_DHS2015,aes(x=TEMP15,y=ln_nightlights))+geom_point()+geom_smooth()
plotnightlights1 + theme(panel.grid.major = element_blank(), axis.text = element_text(size = 15), axis.title=element_blank(),
                         panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
cor.test(data_climate_DHS2015$ln_nightlights, data_climate_DHS2015$TEMP15)

plotGCP<-ggplot(data_climate_DHS2005,aes(x=TEMP05,y=ln_gross_cell_production))+geom_point()+geom_smooth()
plotGCP + theme(panel.grid.major = element_blank(), axis.text = element_text(size = 15), axis.title=element_blank(),
                panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
cor.test(data_climate_DHS2005$ln_gross_cell_production, data_climate_DHS2005$TEMP05)




#####Manuscript - Figure 2#####
scatterplotincome2<-subset(data_climate_Gennaioli,Country%in%c("United States","Colombia","Russian Federation","China"))
plotincome2<-ggplot(scatterplotincome2,aes(x=Temperature,y=Ln_GDP_region,shape=Country,color=Country))+geom_point()+geom_smooth()
plotincome2 + theme(panel.grid.major = element_blank(), axis.text = element_text(size = 15), axis.title=element_blank(),
                    panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

scatterplotnightlights2<-subset(data_climate_DHS2015,country%in%c("Zimbabwe","Malawi","Namibia","Zambia"))
plotnightlights2<-ggplot(scatterplotnightlights2,aes(x=TEMP15,y=ln_nightlights,shape=country,color=country))+geom_point()+geom_smooth()
plotnightlights2 + theme(panel.grid.major = element_blank(), axis.text = element_text(size = 15), axis.title=element_blank(),
                         panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))





#####TABLES#####
#####Manuscript - Table 1: Baseline Income#####
attach(data_climate_Gennaioli)
income1<-felm(Ln_GDP_region~Temperature|0|0|Country)
income2<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
income3<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)+factor(year)|0|Country)
income4<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)+factor(year)|0|Country)
stargazer(income1,income2,income3,income4,type="text",out="table1.txt")
detach(data_climate_Gennaioli)




#####Table 2: Baseline Nightlights#####
attach(data_climate_DHS2015)
nightlights1<-felm(ln_nightlights~TEMP15|0|0|country)
nightlights2<-felm(ln_nightlights~TEMP15|factor(country)|0|country)
nightlights3<-felm(ln_nightlights~TEMP15+Poorlights+(Poorlights*TEMP15)|factor(country)|0|country)
nightlights4<-felm(ln_nightlights~TEMP15+TEMP152|factor(country)|0|country)
stargazer(nightlights1,nightlights2,nightlights3,nightlights4,type="text",out="table2.txt")
detach(data_climate_DHS2015)




#####Table 3: Baseline GCP#####
attach(data_climate_DHS2005)
GCP1<-felm(ln_gross_cell_production~TEMP05|0|0|country)
GCP2<-felm(ln_gross_cell_production~TEMP05|factor(country)|0|country)
GCP3<-felm(ln_gross_cell_production~TEMP05+PoorGCP+(PoorGCP*TEMP05)|factor(country)|0|country)
GCP4<-felm(ln_gross_cell_production~TEMP05+TEMP052|factor(country)|0|country)
stargazer(GCP1,GCP2,GCP3,GCP4,type="text",out="table3.txt")
detach(data_climate_DHS2005)












#####SUPPLEMENTARY INFORMATION#####

#####Figure 3: Within country analysis#####
fits <- lmList(Ln_GDP_region ~ Temperature | country, data=data_climate_Gennaioli)
df <- summary(fits)$coef[,,"Temperature"]
View(df)
hist(df[,"Estimate"], breaks=10, main="Histogram of coefficient estimates", xlab = "coefficients estimates")




#####Figure 4#####
attach(data_climate_Gennaioli)
income5<-gam(Ln_GDP_region~s(Temperature,bs="cr")+factor(country)+factor(year))
summary(income5)
plot(income5)
detach(data_climate_Gennaioli)




######Table 4: Descriptives#####
attach(data_climate_Gennaioli)
descriptives_G<-stat.desc(data_climate_Gennaioli)
View(descriptives_G)
detach(data_climate_Gennaioli)

attach(data_climate_Gennaioli_growth)
descriptives_Gro<-stat.desc(data_climate_Gennaioli_growth)
View(descriptives_Gro)
detach(data_climate_Gennaioli_growth)

attach(data_climate_DHS2005)
descriptives_DHS05<-stat.desc(data_climate_DHS2005)
View(descriptives_DHS05)
detach(data_climate_DHS2005)

attach(data_climate_DHS2015)
descriptives_DHS15<-stat.desc(data_climate_DHS2015)
View(descriptives_DHS15)
detach(data_climate_DHS2015)




#####Table 5: see Figure 3#####
View(df)




#####Table 6: Income with controls#####
attach(data_climate_Gennaioli)
incomecontrols1<-felm(Ln_GDP_region~Temperature+Landlockedregion+nbr+nbr_nr+Latitude+ln_Area_sqkm+Malaria_ecology+Ln_Cum_Oil_Gas_Prod+Ln_Pop_density+CapitalisinRegion+YearsofEducation|factor(Country)+factor(year)|0|Country)
incomecontrols2<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)+Landlockedregion+nbr+nbr_nr+Latitude+ln_Area_sqkm+Malaria_ecology+Ln_Cum_Oil_Gas_Prod+Ln_Pop_density+CapitalisinRegion+YearsofEducation|factor(Country)+factor(year)|0|Country)
incomecontrols3<-felm(Ln_GDP_region~Temperature+Temperature2+Landlockedregion+nbr+nbr_nr+Latitude+ln_Area_sqkm+Malaria_ecology+Ln_Cum_Oil_Gas_Prod+Ln_Pop_density+CapitalisinRegion+YearsofEducation|factor(Country)+factor(year)|0|Country)
stargazer(incomecontrols1,incomecontrols2,incomecontrols3,type="text",out="table6.txt")
detach(data_climate_Gennaioli)




#####Table 7: nightlights with controls#####
attach(data_climate_DHS2015)
rega<-felm(ln_nightlights~TEMP15+Latitude+ln_pop_2015+Aridity_2015+drought_episodes+Enhanced_Vegetation_Index_2015+Frost_Days_2015+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2015+Malaria_Incidence_2015+Malaria_Prevalence_2015+PET_2015+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2015+Slope+Wet_Days_2015|factor(country)|0|country)
regb<-felm(ln_nightlights~TEMP15+Poorlights+(Poorlights*TEMP15)+Latitude+ln_pop_2015+Aridity_2015+drought_episodes+Enhanced_Vegetation_Index_2015+Frost_Days_2015+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2015+Malaria_Incidence_2015+Malaria_Prevalence_2015+PET_2015+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2015+Slope+Wet_Days_2015|factor(country)|0|country)
regc<-felm(ln_nightlights~TEMP15+TEMP152+Latitude+ln_pop_2015+Aridity_2015+drought_episodes+Enhanced_Vegetation_Index_2015+Frost_Days_2015+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2015+Malaria_Incidence_2015+Malaria_Prevalence_2015+PET_2015+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2015+Slope+Wet_Days_2015|factor(country)|0|country)
stargazer(rega,regb,regc,type="text",out="table7.txt")
detach(data_climate_DHS2015)



#####Table 8: GCP with controls#####
attach(data_climate_DHS2005)
regA<-felm(ln_gross_cell_production~TEMP05+Latitude+ln_pop_2005+Aridity_2005+drought_episodes+Enhanced_Vegetation_Index_2005+Frost_Days_2005+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2005+Malaria_Incidence_2005+Malaria_Prevalence_2005+PET_2005+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2005+Slope+Wet_Days_2005|factor(country)|0|country)
regB<-felm(ln_gross_cell_production~TEMP05+PoorGCP+(PoorGCP*TEMP05)+Latitude+ln_pop_2005+Aridity_2005+drought_episodes+Enhanced_Vegetation_Index_2005+Frost_Days_2005+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2005+Malaria_Incidence_2005+Malaria_Prevalence_2005+PET_2005+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2005+Slope+Wet_Days_2005|factor(country)|0|country)
regC<-felm(ln_gross_cell_production~TEMP05+TEMP052+Latitude+ln_pop_2005+Aridity_2005+drought_episodes+Enhanced_Vegetation_Index_2005+Frost_Days_2005+global_human_footprint+growing_season_length+Irrigation+ITN_Coverage_2005+Malaria_Incidence_2005+Malaria_Prevalence_2005+PET_2005+ln_proximity_to_national_borders+ln_Proximity_to_Protected_Areas+ln_proximity_to_water+Rainfall_2005+Slope+Wet_Days_2005|factor(country)|0|country)
stargazer(regA,regB,regC,type="text",out="table8.txt")
detach(data_climate_DHS2005)




#####Table 9: Regional Income in 7 year subsamples#####

year1<-subset(data_climate_Gennaioli,year%in%c("1950"))
attach(year1)
income1950a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income1950b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income1950c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year1)

year2<-subset(data_climate_Gennaioli,year%in%c("1960"))
attach(year2)
income1960a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income1960b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income1960c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year2)

year3<-subset(data_climate_Gennaioli,year%in%c("1970"))
attach(year3)
income1970a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income1970b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income1970c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year3)

year4<-subset(data_climate_Gennaioli,year%in%c("1980"))
attach(year4)
income1980a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income1980b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income1980c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year4)

year5<-subset(data_climate_Gennaioli,year%in%c("1990"))
attach(year5)
income1990a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income1990b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income1990c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year5)

year6<-subset(data_climate_Gennaioli,year%in%c("2000"))
attach(year6)
income2000a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income2000b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income2000c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year6)

year7<-subset(data_climate_Gennaioli,year%in%c("2010"))
attach(year7)
income2010a<-felm(Ln_GDP_region~Temperature|factor(Country)|0|Country)
income2010b<-felm(Ln_GDP_region~Temperature+Poor+(Poor*Temperature)|factor(Country)|0|Country)
income2010c<-felm(Ln_GDP_region~Temperature+Temperature2|factor(Country)|0|Country)
detach(year7)

stargazer(income1950a,income1950b,income1950c,income1960a,income1960b,income1960c,income1970a,income1970b,income1970c,income1980a,income1980b,income1980c,income1990a,income1990b,income1990c,income2000a,income2000b,income2000c,income2010a,income2010b,income2010c,type="text",out="table9.txt")




##### Table 10: Diff_max_min temperatures and december/july temperatures for nightlights#####

attach(data_climate_DHS2015)
regj<-felm(ln_nightlights~Diff_Max_Min15|factor(country)|0|country)
regk<-felm(ln_nightlights~Diff_Max_Min15+Poorlights+(Poorlights*Diff_Max_Min15)|factor(country)|0|country)
regd<-felm(ln_nightlights~Temperature15_December|factor(country)|0|country)
rege<-felm(ln_nightlights~Temperature15_December+Poorlights+(Poorlights*Temperature15_December)|factor(country)|0|country)
regf<-felm(ln_nightlights~Temperature15_December+Temperature152_December|factor(country)|0|country)
regg<-felm(ln_nightlights~Temperature15_July|factor(country)|0|country)
regh<-felm(ln_nightlights~Temperature15_July+Poorlights+(Poorlights*Temperature15_July)|factor(country)|0|country)
regi<-felm(ln_nightlights~Temperature15_July+Temperature152_July|factor(country)|0|country)
detach(data_climate_DHS2015)
stargazer(regj,regk,regd,rege,regf,regg,regh,regi,type="text",out="table10.txt")




#####Table 11: Diff_max_min temperatures and december/july temperatures for GCP#####

attach(data_climate_DHS2005)
regD<-felm(ln_gross_cell_production~Temperature05_December|factor(country)|0|country)
regE<-felm(ln_gross_cell_production~Temperature05_December+PoorGCP+(PoorGCP*Temperature05_December)|factor(country)|0|country)
regF<-felm(ln_gross_cell_production~Temperature05_December+Temperature052_December|factor(country)|0|country)
regG<-felm(ln_gross_cell_production~Temperature05_July|factor(country)|0|country)
regH<-felm(ln_gross_cell_production~Temperature05_July+PoorGCP+(PoorGCP*Temperature05_July)|factor(country)|0|country)
regI<-felm(ln_gross_cell_production~Temperature05_July+Temperature052_July|factor(country)|0|country)
regJ<-felm(ln_gross_cell_production~Diff_Max_Min05|factor(country)|0|country)
regK<-felm(ln_gross_cell_production~Diff_Max_Min05+PoorGCP+(PoorGCP*Diff_Max_Min05)|factor(country)|0|country)
detach(data_climate_DHS2005)
stargazer(regJ,regK,regD,regE,regF,regG,regH,regI,type="text",out="table11.txt")




#####Table 12: Various robustness tests (outliers, educ, weighted, continents, income level)#####

outlier_gdpreg2stdev<-subset(data_climate_Gennaioli,outlier_gdpreg2stdev%in%c("0"))
attach(outlier_gdpreg2stdev)
outlier1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(outlier_gdpreg2stdev)
outlier_conflict<-subset(data_climate_Gennaioli,outlier_conflict%in%c("0"))
attach(outlier_conflict)
outlier2<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(outlier_conflict)
outlier_hottestcoldestreg<-subset(data_climate_Gennaioli,outlier_hottestcoldestreg%in%c("0"))
attach(outlier_hottestcoldestreg)
outlier3<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(outlier_hottestcoldestreg)
outlier_temp2stdev<-subset(data_climate_Gennaioli,outlier_temp2stdev%in%c("0"))
attach(outlier_temp2stdev)
outlier4<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(outlier_temp2stdev)
outlier_nightlights2stdev<-subset(data_climate_DHS2015,outlier_nightlights2stdev%in%c("0"))
attach(outlier_nightlights2stdev)
outliera<-felm(ln_nightlights~TEMP15|factor(country)|0|country)
detach(outlier_nightlights2stdev)
outlier_hottestcoldestreg15<-subset(data_climate_DHS2015,outlier_hottestcoldestreg%in%c("0"))
attach(outlier_hottestcoldestreg15)
outlierb<-felm(ln_nightlights~TEMP15|factor(country)|0|country)
detach(outlier_hottestcoldestreg15)
outlier_temp2stdev15<-subset(data_climate_DHS2015,outlier_temp2stdev15%in%c("0"))
attach(outlier_temp2stdev15)
outlierc<-felm(ln_nightlights~TEMP15|factor(country)|0|country)
detach(outlier_temp2stdev15)
outlier_GCP2stdev<-subset(data_climate_DHS2005,outlier_gcp2stdev%in%c("0"))
attach(outlier_GCP2stdev)
outlierA<-felm(ln_gross_cell_production~TEMP05|factor(country)|0|country)
detach(outlier_GCP2stdev)
outlier_conflict05<-subset(data_climate_DHS2005,outlier_conflict05%in%c("0"))
attach(outlier_conflict05)
outlierB<-felm(ln_gross_cell_production~TEMP05|factor(country)|0|country)
detach(outlier_conflict05)
outlier_hottestcoldestreg05<-subset(data_climate_DHS2005,outlier_hottestcoldestreg05%in%c("0"))
attach(outlier_hottestcoldestreg05)
outlierC<-felm(ln_gross_cell_production~TEMP05|factor(country)|0|country)
detach(outlier_hottestcoldestreg05)
outlier_temp2stdev05<-subset(data_climate_DHS2005,outlier_temp2stdev%in%c("0"))
attach(outlier_temp2stdev05)
outlierD<-felm(ln_gross_cell_production~TEMP05|factor(country)|0|country)
detach(outlier_temp2stdev05)
stargazer(outlier1,outlier2,outlier3,outlier4,outliera,outlierb,outlierc,outlierA,outlierB,outlierC,outlierD,type="text",out="table12outlier.txt")

attach(data_climate_Gennaioli)
incomeweighted<-felm(Ln_GDP_region~Temperature_weighted|factor(Country)+factor(year)|0|Country)
detach(data_climate_Gennaioli)
attach(data_climate_DHS2015)
nightlightsweighted<-felm(ln_nightlights~TEMP15_weighted|factor(country)|0|country)
detach(data_climate_DHS2015)
attach(data_climate_DHS2005)
GCPweighted<-felm(ln_gross_cell_production~TEMP05_weighted|factor(country)|0|country)
detach(data_climate_DHS2005)
stargazer(incomeweighted,nightlightsweighted,GCPweighted,type="text",out="table12weighted.txt")

Africaincome<-subset(data_climate_Gennaioli,continent%in%c("Africa"))
attach(Africaincome)
Africa1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(Africaincome)
Asiaincome<-subset(data_climate_Gennaioli,continent%in%c("Asia"))
attach(Asiaincome)
Asia1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(Asiaincome)
Europeincome<-subset(data_climate_Gennaioli,continent%in%c("Europe"))
attach(Europeincome)
Europe1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(Europeincome)
Africanightlights<-subset(data_climate_DHS2015,continent%in%c("Africa"))
attach(Africanightlights)
Africaa<-felm(ln_nightlights~TEMP15|factor(country)+factor(year)|0|country)
detach(Africanightlights)
Asianightlights<-subset(data_climate_DHS2015,continent%in%c("Asia"))
attach(Asianightlights)
Asiaa<-felm(ln_nightlights~TEMP15|factor(country)+factor(year)|0|country)
detach(Asianightlights)
Europenightlights<-subset(data_climate_DHS2015,continent%in%c("Europe"))
attach(Europenightlights)
Europea<-felm(ln_nightlights~TEMP15|factor(country)+factor(year)|0|country)
detach(Europenightlights)
AfricaGCP<-subset(data_climate_DHS2005,continent%in%c("Africa"))
attach(AfricaGCP)
AfricaA<-felm(ln_gross_cell_production~TEMP05|factor(country)+factor(year)|0|country)
detach(AfricaGCP)
AsiaGCP<-subset(data_climate_DHS2005,continent%in%c("Asia"))
attach(AsiaGCP)
AsiaA<-felm(ln_gross_cell_production~TEMP05|factor(country)+factor(year)|0|country)
detach(AsiaGCP)
EuropeGCP<-subset(data_climate_DHS2005,continent%in%c("Europe"))
attach(EuropeGCP)
EuropeA<-felm(ln_gross_cell_production~TEMP05|factor(country)+factor(year)|0|country)
detach(EuropeGCP)
stargazer(Africa1,Asia1,Europe1,Africaa,Asiaa,Europea,AfricaA,AsiaA,EuropeA,type="text",out="table12continents.txt")

Lowerincome<-subset(data_climate_Gennaioli,WB_classification%in%c("Lower-income"))
attach(Lowerincome)
Lowerincome1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(Lowerincome)
Middleincome<-subset(data_climate_Gennaioli,WB_classification%in%c("Lower-middle-income","Upper-middle-income"))
attach(Middleincome)
Middleincome1<-felm(Ln_GDP_region~Temperature|factor(Country)+factor(year)|0|Country)
detach(Middleincome)
Lowernightlights<-subset(data_climate_DHS2015,WB_classification%in%c("Lower-income"))
attach(Lowernightlights)
Lowernightlightsa<-felm(ln_nightlights~TEMP15|factor(country)+factor(year)|0|country)
detach(Lowernightlights)
Middlenightlights<-subset(data_climate_DHS2015,WB_classification%in%c("Lower-middle-income","Upper-middle-income"))
attach(Middlenightlights)
Middlenightlightsa<-felm(ln_nightlights~TEMP15|factor(country)+factor(year)|0|country)
detach(Middlenightlights)
LowerGCP<-subset(data_climate_DHS2005,WB_classification%in%c("Lower-income"))
attach(LowerGCP)
LowerGCPA<-felm(ln_gross_cell_production~TEMP05|factor(country)+factor(year)|0|country)
detach(LowerGCP)
MiddleGCP<-subset(data_climate_DHS2005,WB_classification%in%c("Lower-middle-income","Upper-middle-income"))
attach(MiddleGCP)
MiddleGCPA<-felm(ln_gross_cell_production~TEMP05|factor(country)+factor(year)|0|country)
detach(MiddleGCP)
stargazer(Lowerincome1,Middleincome1,Lowernightlightsa,Middlenightlightsa,LowerGCPA,MiddleGCPA,type="text",out="table12incomelevels.txt")

attach(data_climate_Gennaioli)
incomeeduc<-felm(Ln_GDP_region~Temperature+YearsofEducation+(Temperature*YearsofEducation)|factor(Country)+factor(year)|0|Country)
detach(data_climate_Gennaioli)
stargazer(incomeeduc,type="text",out="table12educ.txt")




#####Table 13: Growth#####

attach(data_climate_Gennaioli_growth)
growth1<-lm(ggrowth2~gTemperature)
growth2<-felm(ggrowth2~gTemperature|factor(gCountry)|0|gCountry)
growth3<-felm(ggrowth2~gTemperature+gPoor+(gPoor*gTemperature)|factor(gCountry)|0|gCountry)
growth4<-felm(ggrowth2~gTemperature+gTemperature2|factor(gCountry)|0|gCountry)
growth5<-felm(ggrowth2~gTemperature+gLandlockedcountry+gLandlockedregion+ln_gdist_gg+gnbr+gnbr_nr+gLatitude+ln_gArea_sqkm+gMalaria_ecology+gLn_Cum_Oil_Gas_Prod+gLn_Pop_density+gCapitalisinRegion+gLn_GDP_country+gYearsofEducation|factor(gCountry)|0|gCountry)
detach(data_climate_Gennaioli_growth)
goutlier_growth2stdev<-subset(data_climate_Gennaioli_growth,goutlier_growth2stdev%in%c("0"))
attach(goutlier_growth2stdev)
growth6<-felm(ggrowth2~gTemperature|factor(gCountry)|0|gCountry)
detach(goutlier_growth2stdev)
goutlier_hottestcoldestreg<-subset(data_climate_Gennaioli_growth,goutlier_hottestcoldestreg%in%c("0"))
attach(goutlier_hottestcoldestreg)
growth7<-felm(ggrowth2~gTemperature|factor(gCountry)|0|gCountry)
detach(goutlier_hottestcoldestreg)
goutlier_temp2stdev<-subset(data_climate_Gennaioli_growth,goutlier_temp2stdev%in%c("0"))
attach(goutlier_temp2stdev)
growth8<-felm(ggrowth2~gTemperature|factor(gCountry)|0|gCountry)
detach(goutlier_temp2stdev)
stargazer(growth1,growth2,growth3,growth4,growth5,growth6,growth7,growth8,type="text",out="table13.txt")





#####Table 14: Precipitation#####

attach(data_climate_DHS2015)
prec0<-lm(ln_nightlights~PREC15)
prec1<-felm(ln_nightlights~PREC15|factor(country)|0|country)
detach(data_climate_DHS2015)
attach(data_climate_DHS2005)
precA<-lm(ln_gross_cell_production~PREC05)
precB<-felm(ln_gross_cell_production~PREC05|factor(country)|0|country)
detach(data_climate_DHS2005)
stargazer(prec0,prec1,precA,precB,type="text",out="table14.txt")





#####Table 15: ""Replication" of Dell (2009) with regional data from Gennaioli#####

dellreplication1<-subset(data_climate_Gennaioli,Dell_2000>0)
attach(dellreplication1)
dell1<-lm.cluster(data=dellreplication1,Ln_GDP_region~Temperature,cluster="country")
summary(dell1)
dell2<-lm.cluster(data=dellreplication1,Ln_GDP_region~Temperature+ln_dist_gg2+ln_Area_sqkm+Ln_Pop_density,cluster="country")
summary(dell2)
dell3<-felm(Ln_GDP_region~Temperature+ln_dist_gg2+ln_Area_sqkm+Ln_Pop_density|factor(country)|0|country)
detach(dellreplication1)
dellreplication2<-subset(data_climate_Gennaioli,Dell_countries>0)
attach(dellreplication2)
dell4<-felm(Ln_GDP_region~Temperature+ln_dist_gg2+ln_Area_sqkm+Ln_Pop_density|factor(year)+factor(country)|0|country)
detach(dellreplication2)
stargazer(dell3,dell4,type="text",out="table15_dell.txt")




#####Table 16: Replication of Li (2019) with regional data from Gennaioli#####

lireplication1<-subset(data_climate_Gennaioli,country%in%c("China"))
lireplication2<-subset(lireplication1,year>1965)
attach(lireplication2)
li1<-felm(Ln_GDP_region~Temperature|factor(year)|0|region)
li2<-felm(Ln_GDP_region~Temperature+Temperature2|factor(year)|0|region)
detach(lireplication2)
attach(lireplication1)
li3<-felm(Ln_GDP_region~Temperature|factor(year)|0|region)
li4<-felm(Ln_GDP_region~Temperature+Temperature2|factor(year)|0|region)
detach(lireplication1)

stargazer(li1,li2,li3,li4,type="text",out="table16.txt")