rm(list = ls())

library(tidyverse)
library(ggpubr)

setwd("~/Dropbox/C\ O\ N\ S\ U\ T\ O\ R\ I\ A\ S/UNFPA-INE_Paraguay/parity")

base02 <- readxl::read_xlsx("03 Datos_Nacidos vivos x edad de la madre 2002-2012-2022_.xlsx", sheet = "2002_CLEAN2") %>% 
  replace_na(list(nac_viv = 0, muertos = 0, hnv12 = 0)) %>% 
  mutate(year = 2002, .before = depto)

# base12 <- readxl::read_xlsx("03 Datos_Nacidos vivos x edad de la madre 2002-2012-2022_.xlsx", sheet = "2012_CLEAN") %>% 
#   replace_na(list(hnv_1 = 0)) %>% 
#   mutate(year = 2012, .before = depto)

base22 <- readxl::read_xlsx("03 Datos_Nacidos vivos x edad de la madre 2002-2012-2022_.xlsx", sheet = "2022_CLEAN2") %>% 
  replace_na(list(nac_viv = 0, muertos = 0, hnv12 = 0)) %>% 
  mutate(year = 2022, .before = depto)

base02$hnv12[base02$hnv12 == "."] <- 0
base22$hnv12[base22$hnv12 == "."] <- 0

nacim <- rbind(base02, base22) %>% 
  type_convert()

nacim$depto[nacim$depto == "ASUNCION"] <- "ASUNCIÓN"
nacim$depto[nacim$depto == "ALTO PARANA"] <- "ALTO PARANÁ"
nacim$depto[nacim$depto == "BOQUERON"] <- "BOQUERÓN"
nacim$depto[nacim$depto == "CAAGUAZU"] <- "CAAGUAZÚ"
nacim$depto[nacim$depto == "CAAZAPA"] <- "CAAZAPÁ"
nacim$depto[nacim$depto == "CANINDEYU"] <- "CANINDEYÚ"
nacim$depto[nacim$depto == "CONCEPCION"] <- "CONCEPCIÓN"
nacim$depto[nacim$depto == "GUAIRA"] <- "GUAIRÁ"
nacim$depto[nacim$depto == "ITAPUA"] <- "ITAPÚA"
nacim$depto[nacim$depto == "ÑEEMBUCU"] <- "ÑEEMBUCÚ"
nacim$depto[nacim$depto == "PARAGUARI"] <- "PARAGUARÍ"
nacim$depto[nacim$depto == "PDTE. HAYES"] <- "PRESIDENTE HAYES"

nacim_q <- nacim %>% 
  mutate(
    age5 = case_when(
      edad < 15 ~ "10-14",
      edad >= 15 & edad < 20 ~ "15-19",
      edad >= 20 & edad < 25 ~ "20-24",
      edad >= 25 & edad < 30 ~ "25-29",
      edad >= 30 & edad < 35 ~ "30-34",
      edad >= 35 & edad < 40 ~ "35-39",
      edad >= 40 & edad < 45 ~ "40-44",
      edad >= 45 & edad < 50 ~ "45-49",
      edad >= 50 ~ "50+",
    )
  ) %>% 
  group_by(year, depto, age5) %>% 
  summarise(mujeres = sum(mujeres), muj_fec = sum(muj_fec), muj_nfec = sum(muj_nfec),
            nac_viv = sum(nac_viv), muertos = sum(muertos), hnv12 = sum(hnv12)) %>% 
  ungroup()

correct <- readxl::read_xlsx("P_F_correct.xlsx")

base_nal <- nacim_q %>% 
  group_by(year, age5) %>% 
  summarise(mujeres = sum(mujeres), muj_fec = sum(muj_fec), mujnfec = sum(muj_nfec),
            nac_viv = sum(nac_viv), muertos = sum(muertos), hnv12 = sum(hnv12)) %>% 
  ungroup() %>% 
  filter(!age5 %in% c("10-14", "50+")) %>% 
  group_by(year) %>% 
  mutate(agemid = seq(17.5, 47.5, 5),
    Ei = hnv12/mujeres,
         Pi = nac_viv/muj_fec,
         cum_asfr = 5*cumsum(Ei)) %>% 
  select(year, age5, agemid, Ei, Pi, cum_asfr) %>% 
  ungroup()

m <- base_nal %>% 
  group_by(year) %>% 
  summarise(m = sum(agemid*Ei)/sum(Ei)) %>% 
  ungroup()

f1_f2_O <- base_nal %>% 
  group_by(year) %>% 
  summarise(f1f2 = Ei[1]/Ei[2]) %>% 
  ungroup()

yr <- as.numeric(names(table(base_nal$year)))

nacim_est <- list()
for(i in 1:2){
nacim_est[[i]] <- cbind(base_nal %>% 
                        filter(year == yr[i]),
          correct %>% 
            filter(abs(f1_f2 - f1_f2_O$f1f2[i]) == min(abs(f1_f2 - f1_f2_O$f1f2[i])))) %>% 
  select(year, age5, agemid, Ei, Pi, cum_asfr, wi) %>% 
  mutate(Fi = cum_asfr + wi*Ei,
         PF = Pi/Fi, 
         fi_fit = PF[2]*Ei)
}

base_nal_est <- do.call(rbind, nacim_est)

tgf.nal <- base_nal_est %>% 
  group_by(year) %>% 
  summarise(tgf = 5*sum(fi_fit)) %>% 
  ungroup() %>% 
  mutate(source = "censo", .before = "year")

tgfp <- readxl::read_xlsx("03 Datos_Nacidos vivos x edad de la madre 2002-2012-2022_.xlsx", sheet = "TGF") %>% 
  gather(key = "depto", value = "tgf", -Period) %>% 
  rename(year = Period) %>% 
  mutate(source = "proy", .before = "year")

tgf.fin.nal <- rbind(tgf.nal, tgfp)

# departamentos ----

base_dep <- nacim_q %>% 
  filter(!age5 %in% c("10-14", "50+")) %>% 
  group_by(year, depto) %>% 
  mutate(agemid = seq(17.5, 47.5, 5),
         Ei = hnv12/mujeres,
         Pi = nac_viv/muj_fec,
         cum_asfr = 5*cumsum(Ei)) %>% 
  select(year, depto, age5, agemid, Ei, Pi, cum_asfr) %>% 
  ungroup()

m <- base_dep %>% 
  group_by(year, depto) %>% 
  summarise(m = sum(agemid*Ei)/sum(Ei)) %>% 
  ungroup()

f1_f2_O <- base_dep %>% 
  group_by(year, depto) %>% 
  summarise(f1f2 = Ei[1]/Ei[2]) %>% 
  ungroup()

yr <- as.numeric(names(table(base_dep$year)))
dep <- names(table(base_dep$depto))

fec_est_dep1 <- list()
fec_est_dep2 <- list()
for (j in 1:18) {
  for(i in 1:2){
    fec_est_dep1[[i]] <- cbind(base_dep %>% 
                              filter(year == yr[i], depto == dep[j]),
                            correct %>% 
                              filter(abs(f1_f2 - f1_f2_O %>% 
                                           filter(year == yr[i], depto == dep[j]) %>% 
                                           select(f1f2) %>% 
                                           as.numeric()) == min(abs(f1_f2 - 
                                                                      f1_f2_O %>%
                                                                      filter(year == yr[i], 
                                                                             depto == dep[j]) %>% 
                                                                      select(f1f2) %>% 
                                                                      as.numeric())))) %>% 
      select(year, depto, age5, agemid, Ei, Pi, cum_asfr, wi) %>% 
      mutate(Fi = cum_asfr + wi*Ei,
             PF = Pi/Fi, 
             fi_fit = PF[2]*Ei)
  }
  fec_est_dep2[[j]] <- do.call(rbind, fec_est_dep1)
}

base_dep_est <- do.call(rbind, fec_est_dep2)

tgf <- base_dep_est %>% 
  group_by(year, depto) %>% 
  summarise(tgf = 5*sum(fi_fit)) %>% 
  ungroup() %>% 
  mutate(source = "censo", .before = "year")

tgfp <- readxl::read_xlsx("03 Datos_Nacidos vivos x edad de la madre 2002-2012-2022_.xlsx", sheet = "TGF") %>% 
  gather(key = "depto", value = "tgf", -Period) %>% 
  rename(year = Period) %>% 
  mutate(source = "proy", .before = "year")

tgf.fin <- rbind(tgf, tgfp)

tgf.fin$depto[tgf.fin$depto == "ALTO PARANA"] <- "ALTO PARANÁ"
tgf.fin$depto[tgf.fin$depto == "BOQUERON"] <- "BOQUERÓN"
tgf.fin$depto[tgf.fin$depto == "CAAGUAZU"] <- "CAAGUAZÚ"
tgf.fin$depto[tgf.fin$depto == "CAAZAPA"] <- "CAAZAPÁ"
tgf.fin$depto[tgf.fin$depto == "CANINDEYU"] <- "CANINDEYÚ"
tgf.fin$depto[tgf.fin$depto == "CONCEPCION"] <- "CONCEPCIÓN"
tgf.fin$depto[tgf.fin$depto == "GUAIRA"] <- "GUAIRÁ"
tgf.fin$depto[tgf.fin$depto == "ITAPUA"] <- "ITAPÚA"
tgf.fin$depto[tgf.fin$depto == "ÑEEMBUCU"] <- "ÑEEMBUCÚ"
tgf.fin$depto[tgf.fin$depto == "PARAGUARI"] <- "PARAGUARÍ"
tgf.fin$depto[tgf.fin$depto == "PDTE. HAYES"] <- "PRESIDENTE HAYES"

dep <- names(table(tgf.fin$depto))

for (i in 1:18) {
tgf.fin %>% 
  filter(depto == dep[i]) %>% 
  ggplot(aes(x = year, y = tgf, color = source)) +
  geom_point() +
  labs(x = "Año", y = "Tasa Global de Fecundidad", title = dep[i]) +
    scale_color_discrete(name = "Fuente", labels = c("Censos", "Proyecciones"),
                         type = c("blue", "red")) +
    scale_x_continuous(n.breaks = 15) +
    scale_y_continuous(n.breaks = 12) +
  theme_light()
  
ggsave(filename = paste("tgf_", dep[i], ".png", sep  = ""), dpi = 320)
}
