rm(list = ls())

library(tidyverse)
library(ggpubr)

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

load("censos_proy.RData")

wppM <- readxl::read_xlsx("WPP2022_POP_F01_2_POPULATION_SINGLE_AGE_MALE.xlsx",
                  sheet = "Estimates", range = "A17:DH20541")
colnames(wppM)[3] <- "country"

wppF <- readxl::read_xlsx("WPP2022_POP_F01_3_POPULATION_SINGLE_AGE_FEMALE.xlsx",
                          sheet = "Estimates", range = "A17:DH20541")
colnames(wppF)[3] <- "country"

pobM <- wppM %>% 
  filter(country == "Paraguay") %>% 
  select(-c(1:10)) %>% 
  gather(key = "age", value = "pob", -Year) %>% 
  type_convert() %>% 
  mutate(pob = 1000*pob)
pobF <- wppF %>% 
  filter(country == "Paraguay") %>% 
  select(-c(1:10)) %>% 
  gather(key = "age", value = "pob", -Year) %>% 
  type_convert() %>% 
  mutate(pob = 1000*pob)

pobM$age[pobM$age == "100+"] <- 100
pobF$age[pobF$age == "100+"] <- 100

pobM <- pobM %>% 
  type_convert() %>% 
  mutate(source = "wpp",
    sex = "Hombres", .before = Year) %>% 
  select(source, sex, year = Year, age, pob)
pobF <- pobF %>% 
  type_convert() %>% 
  mutate(source = "wpp",
         sex = "Mujeres", .before = Year) %>% 
  select(source, sex, year = Year, age, pob)

pobine <- rbind(base0 %>% 
  filter(type == "proy") %>% 
  type_convert() %>% 
  group_by(year, sex, age) %>% 
  summarise(pob = sum(pob)) %>% 
  mutate(source = "ine_proy") %>% 
  select(source, sex, year, age, pob),
  base0 %>% 
    filter(type == "cens", age != "NE") %>% 
    type_convert() %>% 
    group_by(year, sex, age) %>% 
    summarise(pob = sum(pob)) %>% 
    mutate(source = "ine_censo") %>% 
    select(source, sex, year, age, pob),
  base0 %>% 
  filter(type == "cens", age == "NE") %>% 
  mutate(age = 999) %>% 
  type_convert() %>% 
  group_by(year, sex, age) %>% 
  summarise(pob = sum(pob)) %>% 
  mutate(source = "ine_censo") %>% 
  select(source, sex, year, age, pob))


base1 <- rbind(pobF, pobM, pobine)

wpp.trunc <- base1 %>% 
  filter(source == "wpp") %>% 
  group_by(year) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(year >= 2000, year <= 2008)

# library(forecast)
# 
wpp.trunc.fit <- auto.arima(wpp.trunc$pob, trace = T, d = 1)
wpp.trunc.for <- forecast(wpp.trunc.fit, h = 14, level = c(95))

ggplot(data = base1 %>% 
         filter(year >= 2000, year <= 2022) %>% 
         group_by(source, year) %>% 
         summarise(pob = sum(pob)) %>% 
         ungroup() %>% 
         filter(source != "ine_censo")) +
  geom_line(aes(x = year, y = pob/1000000, color = source)) +
  geom_line(data = tibble(year = c(2000:2022), 
                          pob = c(as.numeric(wpp.trunc.for$x)/1000000, 
                                  as.numeric(wpp.trunc.for$mean)/1000000)),
    aes(x = year, y = pob), 
            color = "black", linetype = "dashed") +
  geom_point(data = base1 %>% 
               filter(year >= 2000) %>% 
               group_by(source, year) %>% 
               summarise(pob = sum(pob)) %>% 
               ungroup() %>% 
               filter(source == "ine_censo"),
             aes(x = year, y = pob/1000000, color = "Censos")) +
  labs(x = "Año", y = "Población (miles)", color = "", fill = NULL)  +
  theme_light() +
  scale_color_discrete(labels = c("Censos","Proyección DGEEC", "WPP-22"),
                       type = c("black", "red", "blue")) +
  scale_x_continuous(n.breaks = 15) +
  theme(legend.position = "bottom") +
  annotate(geom = "text", x = 2020, y = 6.1, label = "Estimación propia",
           size = 3)
ggsave(filename = "compar_proys.png", dpi = 320)

base1 %>% 
  group_by(year) %>% 
  mutate(prop = ifelse(sex == "Hombres", -100*pob/sum(pob), 100*pob/sum(pob))) %>% 
  ungroup() %>% 
  filter(year == 2002, source == "ine_proy") %>% 
  ggplot(aes(x = age, y = prop)) +
  geom_bar(aes(fill = age), stat = "identity", show.legend = F) +
  annotate(geom ="text", x = 105, y = -0.4, label="Hombres",
           color ="black") + 
  annotate(geom ="text", x = 105, y = 0.4, label="Mujeres",
           color ="black") + 
  geom_hline(yintercept = 0) + 
  coord_flip() + 
  labs(y = "Proporción de población", x = "Edad", title = "Proyección INE 2002") +
  scale_y_continuous(limits = c(-0.6, 0.6), breaks = seq(-0.6, 0.6, 0.2), 
                     labels = as.character(c(seq(0.6, 0, -0.2), 
                                             seq(0.2, 0.6, 0.2)))) +
  scale_x_continuous(limits = c(-1, 110), breaks = seq(0, 110, 5),
                     labels = seq(0, 110, 5)) +
  scale_fill_viridis_c(option = "D")
  
ggsave(filename = "proy_ine_2002.png", dpi = 320, width = 16.5, height = 11.5)
