####################
### INTRODUCTION ###
####################
# This script calculates growth-emission elasticities
# It requires the following packages: dplyr, joyn, lme4, DescTools, tidyr
library(dplyr)
rm(list=ls())

####################
### PREPARE DATA ###
####################
load("02-intermediatedata/ghg.Rda")
load("02-intermediatedata/gdp.Rda")
load("02-intermediatedata/pop.Rda")

elasticities <- joyn::joyn(gdp,ghg,by=c("code","year"),match_type="1:1",keep="inner",reportvar=FALSE) |>
                filter(!is.na(gdppc) & year>=2014) |>
                mutate(codenr = as.numeric(as.factor((code)))) |>
                joyn::joyn(pop,by=c("code","year"),match_type = "1:1", keep="inner",reportvar=FALSE) |>
                mutate(co2pc = co2/pop,
                       ch4pc = ch4/pop,
                       n2opc = n2o/pop) |>
                group_by(code) |>
                arrange(year) |>
                mutate(ggdppc = log(gdppc)-lag(log(gdppc)),
                       gco2pc = log(co2pc)-lag(log(co2pc)),
                       gch4pc = log(ch4pc)-lag(log(ch4pc)),
                       gn2opc = log(n2opc)-lag(log(n2opc))) |>
                ungroup()

# Estimating and printing models
# This output lies behind Table S1
(model_co2 <- lme4::lmer(gco2pc ~ ggdppc + (1 +ggdppc || codenr), data=elasticities))
(model_ch4 <- lme4::lmer(gch4pc ~ ggdppc + (1 +ggdppc || codenr), data=elasticities))
(model_n2o <- lme4::lmer(gn2opc ~ ggdppc + (1 +ggdppc || codenr), data=elasticities))


codes <- elasticities |>
         select(code,codenr) |>
         distinct()

elasticities <- as.data.frame(lme4::ranef(model_co2)) |>
                 select(grp,term,condval) |>
                 rename("elasticity_co2" = "condval") |>
                 joyn::joyn(as.data.frame(lme4::ranef(model_ch4)),by=c("grp","term"),reportvar=FALSE) |>
                 rename("elasticity_ch4" = "condval") |>
                 joyn::joyn(as.data.frame(lme4::ranef(model_n2o)),by=c("grp","term"),reportvar=FALSE) |>
                 rename("elasticity_n2o" = "condval") |>
                 filter(term=="ggdppc") |>
                 select(grp,starts_with("elasticity")) |>
                 joyn::joyn(codes,by=c("grp=codenr"),match_type="1:1",reportvar=FALSE) |>
                 mutate(elasticity_co2_uw = elasticity_co2 + summary(model_co2)$coefficients[2,1],
                        elasticity_ch4_uw = elasticity_ch4 + summary(model_ch4)$coefficients[2,1],
                        elasticity_n2o_uw = elasticity_n2o + summary(model_n2o)$coefficients[2,1]) 
                 
elasticities$elasticity_co2 = DescTools::Winsorize(elasticities$elasticity_co2_uw, val = quantile(elasticities$elasticity_co2_uw, probs = c(0.1, 0.9), na.rm = TRUE))
elasticities$elasticity_ch4 = DescTools::Winsorize(elasticities$elasticity_ch4_uw, val = quantile(elasticities$elasticity_ch4_uw, probs = c(0.1, 0.9), na.rm = TRUE))
elasticities$elasticity_n2o = DescTools::Winsorize(elasticities$elasticity_n2o_uw, val = quantile(elasticities$elasticity_n2o_uw, probs = c(0.1, 0.9), na.rm = TRUE))

elasticities <- elasticities |>
                 select(-grp) |>
                 tidyr::pivot_longer(-code,values_to="elasticity",names_to="gas") |>
                 mutate(type = if_else(grepl("uw", gas), "elasticity_unwinsorized", "elasticity"),
                        gas  = gsub("elasticity_|_uw","",gas),
                        gas_label = factor(toupper(gas), levels = c("CO2", "CH4", "N2O"), # The internal factor levels 
                                          labels = c(expression(CO[2]),expression(CH[4]), expression(N[2]*O)))) |>
                 tidyr::pivot_wider(names_from=type,values_from=elasticity) 
  
save(elasticities, file="02-intermediatedata/elasticities.Rda")


