#### Observatorio Demográfico 2025                                        ####
#### Anexo: Cálculo del efecto tempo en las tasas globales de fecundidad  
#### CELADE - División de Población de la CEPAL                           
#### Fecha: 16/10/2025                                                      

#### 1. Instalación de paquetes             ####  

rm(list = ls( ))
graphics.off( )

# Inserte la ruta de la carpeta donde se encuentra este "Anexo" en su computador:
#setwd(".../Anexo/Datos")


# Lista de paquetes requeridos
.packages <-  
  c( "readr", "readxl", "dplyr", 
     "tidyr", "ggplot2", "scales", "zoo")

# Instalación y carga de paquetes
.inst <- 
  .packages %in% installed.packages()

if( length( .packages[ !.inst ] ) > 0 ){
  install.packages( 
    .packages[ !.inst ], 
    dependencies = TRUE 
  )
}

lapply( .packages, require, character.only = TRUE )

# Definir tema base para ggplot
wpp_theme <- 
  function( ){
    theme_bw( ) +
      theme(
        legend.position = 'bottom',     
        legend.text      = element_text( color = 'black', size = 11 ), 
        panel.grid.major.y = element_line( color = 'gray75', linewidth = 0.30,
                                           linetype = 'longdash' ), 
        panel.grid.major.x = element_blank(), 
        panel.grid.minor = element_blank(),
        plot.caption  = element_text( color = 'black', size = 10, hjust = 0 ),
        plot.title    = element_text( color = 'black', size = 17 ),
        plot.subtitle = element_text( color = 'black', size = 13 ),
        strip.text    = element_text( color = 'black', size = 15 ),
        strip.background = element_rect( color = 'white', fill = 'gray93' ),
        axis.text.x   = element_text( color = 'black', size = 11 ),
        axis.text.y   = element_text( color = 'black', size = 11 ),
        axis.title    = element_text( color = 'black', size = 13 )
      )
  }


#### 2. Estimación de la Tasa Global de Fecundidad ajustada por efecto tempo ####
#### Insumo: WPP 2024. Ajuste basado en la Edad Media de la Fecundidad

# Procesamiento principal de los datos
#
# Cálculo de indicadores clave:
#   * TFR: Tasa Global de Fecundidad (suma de las tasas específicas de fecundidad —ASFR— / 1000, multiplicada implícitamente por 1 año por grupo etario).
#   * MAC: Edad Media de la Fecundidad (Media ponderada de las edades centrales de los grupos etarios).
#
# Luego, por región o país, se estiman los efectos tempo y los ajustes:
#   * r: Variación interanual aproximada de la MAC (usando lead y lag para diferencias centradas).
#   * tfradj: TFR ajustada por efecto tempo (según Bongaarts: TFR / (1 - r)).
#   * Suavizados móviles (ma3 y ma5) para estabilizar las estimaciones.
#   * tfradj_hibrido: Combinación híbrida que permite transiciones suaves en proyecciones
#                     (histórico suavizado, reciente menos suavizado y año más reciente sin suavizar).

# Cargar datos base: tasas específicas de fecundidad (ASFR) por año, edad y región/país, provenientes de WPP 2024
asfr <- read_csv2('../Datos/asfr.csv')

ind3_dt <- asfr %>% 
  mutate(age = as.numeric(age)) %>% 
  filter(year >= 1950 & year < 2028) %>% 
  group_by(LocationID, LocAreaEs, LocationNameCepal, year) %>% 
  summarize(
    tfr = sum(asfr / 1000, na.rm = TRUE), 
    mac = sum((age + 0.5) * (asfr / 1000), na.rm = TRUE) / sum(asfr / 1000, na.rm = TRUE),
    .groups = "drop"
  ) %>% 
  arrange(LocationID, year) %>% 
  group_by(LocationID) %>% 
  mutate(
    r = (lead(mac) - lag(mac)) / 2, 
    tfradj = tfr / (1 - r), 
    tfradj_ma3 = zoo::rollmean(tfradj, k = 3, fill = NA, align = "center"), 
    tfradj_ma5 = zoo::rollmean(tfradj, k = 5, fill = NA, align = "center"), 
    tfradj_hibrido = case_when(
      year < 2023 ~ tfradj_ma5,
      year == 2023.5 ~ tfradj_ma3, 
      year >= 2024 ~ tfradj
    )
  ) %>% 
  ungroup()


#### 3. Elaboración del gráfico: América Latina y el Caribe: edad media de la fecundidad, 2000-2023 (Gráfico III.5 del Observatorio Demográfico 2025) ####

# Verificar si el directorio de salida existe; crearlo si no
output_dir <- "../Gráficos"
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

ind3_dta <- ind3_dt %>%
  filter(LocationID == 904 & year %in% c(2000.5:2024.5)) 
  
 g1 <-  ggplot(ind3_dta, aes(x = year, y = mac)) +
  geom_line(linewidth = 2, alpha = 0.9) +
  scale_y_continuous(
    name = 'Edad media de la fecundidad',
    labels = function(x) format(x, decimal.mark = ","),
    limits = c(26.5, 28),
    breaks = seq(26.5, 28, 0.5)) +
  scale_x_continuous(
    name = NULL,
    breaks = seq(1950, 2025, 5)) +
    labs(
    title = "América Latina y el Caribe: edad media de la fecundidad, 2000–2024  (En años)"
     ) +
  wpp_theme() +
  theme(
    plot.title = element_text(
      size = 16,          
      hjust = 0.5        
      )
  )

ggsave(filename = "../Gráficos/III.5_EMM_ALC.png",
       plot = g1, width = 16, height = 10, units = "in", dpi = 300, device = "png")

#### 4. Elaboración del gráfico: América Latina y el Caribe: tasa global de fecundidad (TGF) de período y ajustada por efecto tempo, 2000-2023 (Gráfico III.6 del Observatorio Demográfico 2025) ####

g2<- ggplot(ind3_dta %>% filter(year > 1999, year <= 2023.5)) +
  geom_line(aes(x = year, y = tfr, color = "de periodo"),
            linewidth = 2, alpha = 0.9) +
  geom_line(aes(x = year, y = tfradj_hibrido, color = "ajustada por tempo"),
            linewidth = 2, alpha = 0.9) +
  scale_y_continuous(name = 'Tasa global de fecundidad', 
                     labels = function(x) format(x, decimal.mark = ","),
                     breaks = pretty_breaks(n = 6)) +
  scale_x_continuous(
    name = NULL,
    breaks = seq(1950, 2025, 5)) +
  labs(
    title = "América Latina y el Caribe: tasa global de fecundidad (TGF) de período y ajustada por efecto tempo, 2000-2023"
  ) +
  scale_color_manual(name = "TGF",
                     values = c("ajustada por tempo" = "red",
                                "de periodo" = "black")) +
  guides(color = guide_legend(nrow = 1)) +
  wpp_theme() +
  theme(
    plot.title = element_text(
      size = 16,          
      hjust = 0.5        
    )
  )

ggsave(filename = "../Gráficos/III.6_TGFajustada_ALC.png",
       plot = g2, width = 12, height = 8, units = "in", dpi = 300, device = "png")

### 5. Elaboración del gráfico: América Latina (20 países): tasa global de fecundidad (TGF) de período y ajustada por efecto tempo, 2000-2023 (gráfico III.7 del Observatorio Demográfico 2025) ####
plot_data <- ind3_dt  %>%
filter(LocAreaEs %in% "América Latina" , 
      LocationID <= 899, 
      year >= 2000.5 & year <= 2023.5)

# Crear el gráfico con ggplot2
  g3 <- ggplot(plot_data) +
    # Línea para TGF de período (sin ajuste)
    geom_line(aes(x = year, y = tfr, color = "de periodo"),
              linewidth = 1.2, alpha = 0.9) +  
    
    # Línea para TGF ajustada por tempo (híbrida)
    geom_line(aes(x = year, y = tfradj_hibrido, color = "ajustada por efecto tempo"),
              linewidth = 1.2, alpha = 0.9) +
       
     scale_y_continuous(name = 'Tasa global de fecundidad', 
                       labels = function(x) format(x, decimal.mark = ","),
                       breaks = scales::pretty_breaks(n = 6)) +
        scale_x_continuous(name = NULL, 
                       breaks = seq(2000, 2020, 5)) +  
        scale_color_manual(name = "TGF",
                       values = c("de periodo" = "black",
                                  "ajustada por efecto tempo" = "red")) +
        guides(color = guide_legend(nrow = 1)) +
        facet_wrap(~ LocationNameCepal, scales = "free_x", ncol = 5) +
        wpp_theme() +
        theme(legend.position = "bottom", 
          strip.text = element_text(size = 10),
          axis.text.x = element_text(angle = 45, hjust = 1))
  
  # Guardar el gráfico en formato PNG
    ggsave(filename = "../Gráficos/III.7_TGFajustada_AL.png",
         plot = g3, width = 16, height = 14, units = "in", dpi = 300, device = "png")
  

### 6. Elaboración del gráfico: El Caribe (27 países y territorios): tasa global de fecundidad (TGF) de período y ajustada por efecto tempo, 2000-2023 (Gráfico III.8 del Observatorio Demográfico 2025) ####
  
# Gráfico para El Caribe
plot_data <- ind3_dt  %>%
  filter(LocAreaEs %in% "El Caribe" , 
         LocationID <= 899, 
         year >= 2000.5 & year <= 2023.5)

  # Crear el gráfico
 g4 <- ggplot(plot_data) +
    geom_line(aes(x = year, y = tfr, color = "de periodo"),
              linewidth = 1.2, alpha = 0.9) +  
    
    geom_line(aes(x = year, y = tfradj_hibrido, color = "ajustada por efecto tempo"),
              linewidth = 1.2, alpha = 0.9) +
    
    scale_y_continuous(name = 'Tasa global de fecundidad', 
                       labels = function(x) format(x, decimal.mark = ","),
                       breaks = scales::pretty_breaks(n = 6)) +
    
    scale_x_continuous(name = NULL, 
                       breaks = seq(2000, 2020, 5)) + 
    
    scale_color_manual(name = "TGF",
                       values = c("de periodo" = "black",
                                  "ajustada por efecto tempo" = "red")) +
    
    guides(color = guide_legend(nrow = 1)) +
    
    facet_wrap(~ LocationNameCepal, scales = "free_x", ncol = 5) +
    
    wpp_theme() +
    
    theme(legend.position = "bottom",
          strip.text = element_text(size = 10),
          axis.text.x = element_text(angle = 45, hjust = 1))
  
  # Guardar en PNG
  ggsave(filename = "../Gráficos/III.8_TGFajustada_CA.png",
         plot = g4, width = 16, height = 14, units = "in", dpi = 300, device = "png")
  
### 7. Estimación del efecto tempo en la TGF de Chile mediante el método basado en el orden de nacimiento y distintos métodos de suavizamiento ####
##     Datos utilizados: denominador: WPP 2024, numerador: nacidos vivos por orden de nacimiento de las estadísticas vitales

Bx_total <- read_csv2("../Datos/Bx_EV_Chile.csv")
Exp1 <- read_csv("../Datos/WPP_POPFEM_CHILE_1990_2024.csv")

Exp1 <- Exp1 %>%
  mutate(Age = ifelse(Age == "100+", 100, Age)) %>%
  mutate(Age = as.integer(Age)) %>%
  filter(Age > 12, Age < 55) %>%
  rename(Exposure = Exp)  

# 7.1 Cálculo de la TFR ajustada detallada por orden de nacimiento ####
# Procedimiento: 
# - Combinar los datos de nacimientos y exposiciones.
# - Calcular las tasas específicas por orden de nacimiento (Fxi).
# - Para cada año y orden: calcular TFRi (suma de Fxi) y MABi (edad media al nacimiento por orden).
# - Para cada orden de nacimiento: estimar ri (variación temporal de la MAB) y TFRadji (ajuste según Bongaarts-Feeney).
# - Redondear los resultados para su presentación.

TFRadj_detallado <- Bx_total %>%
  left_join(Exp1, by = join_by(Year, Age)) %>%  # Join eficiente con join_by (dplyr 1.1+)
  filter(!is.na(Age), !is.na(Exposure), !is.na(Bxi)) %>%  # Filtrar NA tempranamente para evitar propagación
  mutate(Fxi = Bxi / Exposure) %>%  # Tasa específica por orden
  group_by(Year, order) %>%
  summarise(
    TFRi = sum(Fxi, na.rm = TRUE),  # TFR por orden (suma tasas)
    MABi = if_else(TFRi > 0, sum(Fxi * Age, na.rm = TRUE) / TFRi, NA_real_),  # Edad media, condicional para evitar div/0
    .groups = "drop"
  ) %>%
  arrange(order, Year) %>%
  group_by(order) %>%
  mutate(
    ri = (lead(MABi) - lag(MABi)) / 2,  # Efecto tempo: cambio anual centrado en MAB
    TFRadji = if_else(abs(1 - ri) > 0.0001, TFRi / (1 - ri), NA_real_)  # Ajuste tempo, evitar div cerca de 0
  ) %>%
  ungroup() %>%
  mutate(
    TFRi = round(TFRi, 3),
    TFRadji = round(TFRadji, 3)
  )

# 7.2. Reorganización de la TFR ajustada por orden para columnas separadas (adjTFR1, adjTFR2, etc.) ####
TFRadj_por_orden <- TFRadj_detallado %>%
  select(Year, order, TFRadji) %>%  # dplyr::select para claridad
  mutate(nombre = paste0("adjTFR", order)) %>%
  pivot_wider(names_from = nombre, values_from = TFRadji) %>%
  arrange(Year)  # Asegurar orden cronológico

# 7.3. Calcular TFR total, observada y ajustada + suavizados promedios móviles ####
# - Sumar por año: TGF (observada total), TFRadj (ajustada total)
# - Híbrido: Sin ajuste en el último año, MA3 en penúltimo año, MA5 en históricos
TFRadj_total <- TFRadj_detallado %>%
  group_by(Year) %>%
  summarise(
    TFRadj = if_else(any(!is.na(TFRadji)), round(sum(TFRadji, na.rm = TRUE), 3), NA_real_),
    TGF = if_else(any(!is.na(TFRi)), round(sum(TFRi, na.rm = TRUE), 3), NA_real_),
    .groups = "drop"
  ) %>%
  arrange(Year) %>%
  mutate(
    TFRadj_3ma = round(zoo::rollmean(TFRadj, k = 3, fill = NA, align = "center"), 3),
    TFRadj_5ma = round(zoo::rollmean(TFRadj, k = 5, fill = NA, align = "center"), 3)
  ) %>%
  mutate(
    TFRadj_hibrido = case_when(
      Year == max(Year, na.rm = TRUE) ~ TFRadj,
      Year == max(Year, na.rm = TRUE) - 1 ~ TFRadj_3ma,
      TRUE ~ TFRadj_5ma
    )
  )

# 7.4. Unir TFR total (con suavizados) y por orden (columnas adjTFRx) ####
TFRadj_total_y_ordenes <- TFRadj_total %>%
  left_join(TFRadj_por_orden, by = "Year")


# 7.5. Calcular MAB completa: Total y por orden de nacimiento (agrupando órdenes >=5 como "5p") ####

MAB_por_orden <- Bx_total %>%
  left_join(Exp1, by = join_by(Year, Age)) %>%
  filter(!is.na(Exposure), !is.na(Bxi), !is.na(Age)) %>%  # Filtrar NA tempranamente
  mutate(
    Fxi = Bxi / Exposure,
    edad_media = Age + 0.5
  ) %>%
  group_by(Year, order) %>%
  summarise(
    TFRi = sum(Fxi, na.rm = TRUE),
    MABi = if_else(TFRi > 0, sum(Fxi * edad_media, na.rm = TRUE) / TFRi, NA_real_),  # Evitar div/0
    .groups = "drop"
  ) %>%
  mutate(
    order_group = if_else(order >= 5, "5p", as.character(order)),
    nombre = paste0("MAB", order_group),
    MABi = round(MABi, 2)
  ) %>%
  select(Year, nombre, MABi) %>%
  distinct(Year, nombre, .keep_all = TRUE) %>%  # Asegurar unicidad
  pivot_wider(names_from = nombre, values_from = MABi) %>%
  arrange(Year)

MAB_total <- Bx_total %>%
  left_join(Exp1, by = join_by(Year, Age)) %>%
  filter(!is.na(Exposure), !is.na(Bxi), !is.na(Age)) %>%
  mutate(
    Fxi = Bxi / Exposure,
    edad_media = Age + 0.5
  ) %>%
  group_by(Year) %>%
  summarise(
    TFR_total = sum(Fxi, na.rm = TRUE),
    MAB = if_else(TFR_total > 0, round(sum(Fxi * edad_media, na.rm = TRUE) / TFR_total, 2), NA_real_),
    .groups = "drop"
  ) %>%
  select(Year, MAB)  # Solo MAB total para join

MAB_completo <- MAB_total %>%
  left_join(MAB_por_orden, by = "Year") %>%
  arrange(Year)

# Asignar país (para merges posteriores)
MAB_completo$Pais <- "Chile"

# 7.6. Calcular TFR ajustada por método de tasas específicas totales (sin distinción por orden) ####
# - r: efecto tempo, TFRadj: ajuste Bongaarts-Feeney
# - Suavizado simple MA3 manual (lag/lead)
TFRadj_especificas <- Bx_total %>%
  left_join(Exp1, by = join_by(Year, Age)) %>%
  filter(!is.na(Age), !is.na(Exposure), !is.na(Bxi)) %>%
  mutate(Fx = Bxi / Exposure) %>%
  group_by(Year) %>%
  summarise(
    TFR = sum(Fx, na.rm = TRUE),
    MAC = if_else(TFR > 0, sum((Age + 0.5) * Fx, na.rm = TRUE) / TFR, NA_real_),
    .groups = "drop"
  ) %>%
  arrange(Year) %>%
  mutate(
    r = (lead(MAC) - lag(MAC)) / 2,
    TFRadj = if_else(abs(1 - r) > 0.0001, TFR / (1 - r), NA_real_),  # Evitar inestabilidad numérica
    TFRadj_ma = (lag(TFRadj) + TFRadj + lead(TFRadj)) / 3  # MA3 centrado manual
  ) %>%
  mutate(
    TFR = round(TFR, 3),
    TFRadj = round(TFRadj, 3),
    TFRadj_ma = round(TFRadj_ma, 3)
  )

# Usar TFRadj_total_y_ordenes del script previo (asumido calculado con ajuste por orden)
# Asignar país
TFRadj_total_y_ordenes$Pais <- "Chile"

# 7.7. Preparar la base de datos final para comparar métodos ####
# - metodo1: Basado en ajuste por orden (del script anterior: TFRadj_orden es suma ajustada por orden)
# - metodo2: Basado en tasas específicas totales
# - Unir por Year para permitir comparaciones y visualizaciones gráficas lado a lado

metodo1 <- TFRadj_total_y_ordenes %>%
  select(
    Year,
    TFRadj_orden = TFRadj,      # Ajustada total por órdenes
    TFRadj_3ma, 
    TFRadj_5ma, 
    TFRadj_hibrido,
    TGF_orden = TGF             # TGF observada por órdenes
  )

metodo2 <- TFRadj_especificas %>%
  select(
    Year,
    TFRadj_total = TFRadj,      # Ajustada por MAC total
    MAC,
    TFRadj_ma,
    TGF_total = TFR                 # TGF observada total
  )

comparacion <- left_join(metodo1, metodo2, by = "Year") %>%
  arrange(Year)

### 8. Elaboración del gráfico: Chile: comparación entre la tasa global de fecundidad (TGF) observada, ajustada por efecto tempo y suavizada, 1992-2022 (gráfico III.9 del Observatorio Demográfico 2025) ####

g5<- ggplot(comparacion, aes(x = Year)) +
  geom_line(aes(y = TGF_total, color = "TGF observada"), linewidth = 2) +
  geom_line(aes(y = TFRadj_5ma, color = "TGF ajustada por orden suavizada 5 años"), linewidth = 2) +
  geom_line(aes(y = TFRadj_3ma, color = "TGF ajustada por orden suavizada 3 años"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_total, color = "TGF ajustada por EMM"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_ma, color = "TGF ajustada por EMM suavizada 3 años"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_orden, color = "TGF ajustada por orden"), linewidth = 1.5) +
  geom_hline(yintercept = 2.1, linetype = "dashed", color = "gray50") +
  scale_y_continuous(
    name = "Tasa global de fecundidad",
    breaks = sort(unique(c(scales::pretty_breaks(n = 7)(range(comparacion[, grep("TGF|TFR", names(comparacion))], na.rm = TRUE)), 2.1))),
    labels = function(x) format(x, decimal.mark = ",")
  ) +
  scale_x_continuous(name = "Año", breaks = scales::pretty_breaks(n = 10)) +
  scale_color_manual(name = "",
                     values = c(
                       "TGF observada" = "black",
                       "TGF ajustada por orden suavizada 5 años" = "#9932CC",
                       "TGF ajustada por orden suavizada 3 años" = "#FF1493",
                       "TGF ajustada por EMM" = "lightblue",
                       "TGF ajustada por EMM suavizada 3 años" = "#FC8D62",
                       "TGF ajustada por orden" = "#1B9E77"
                     )
  ) +
  labs(
    title = "Chile: comparación entre la tasa global de fecundidad (TGF) observada, ajustada por efecto tempo y suavizada, 1992-2022",
    y = "Tasa global de fecundidad",
    x = "Año",
    color = ""
  ) +
  wpp_theme() +
  theme(
    plot.title = element_text(size=17, hjust=0.5),
    legend.position = "bottom",
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  guides(color = guide_legend(nrow = 2, byrow = TRUE))

ggsave(filename = "../Gráficos/III.9_TGFajustadasuavizada_CHL.png",plot = g5, width = 14, height = 10, units = "in", dpi = 300, device = "png")


### 9. Estimación del efecto tempo en la TGF de Argentina utilizando el método basado en el orden de nacimiento y distintos suavizamientos ####
## Datos utilizados: denominador: WPP 2024; numerador: nacidos vivos por orden de nacimiento según las estadísticas vitales
#  Fuente: Dirección de Estadísticas e Información en Salud del Ministerio de Salud de la Nación
#  Se revisaron y corrigieron inconsistencias evidentes entre la edad de la madre y la cantidad de nacidos vivos, especialmente en los años más antiguos.
#  Por esta razón, los totales de las variables podrían no coincidir exactamente con los publicados originalmente.

Bx <- read_excel("../Datos/Orden_nac_arg.xlsx")
Exp <- read.csv("../Datos/WPP_POPFEM_ARG_1990_2024.csv")

Bx <- Bx %>% 
  rename(Year = Año,
         Age = Edad_madre,
         Bxi = Cantidad,
         order = Orden_NacViv) %>% 
  filter(Age > 0) %>%
  filter(order != 0) %>%
  mutate(order = case_when(
    as.numeric(order) <= 4 ~ as.character(as.numeric(order)),
    as.numeric(order) >= 5 ~ "5",
    TRUE ~ NA_character_)) %>%
  group_by(Year, Age, order) %>%
  summarise(Bxi = sum(Bxi, na.rm = TRUE), .groups = "drop")

Exp <- Exp %>%
  mutate(Age = ifelse(Age == "100+", 100, Age))

Bx <- Bx %>% mutate(Year = as.integer(Year))
Exp <- Exp %>% mutate(Year = as.integer(Year))
Bx <- Bx %>% mutate(Age = as.numeric(Age))
Exp <- Exp %>% mutate(Age = as.numeric(Age))
Exp <- Exp %>% filter(Age>= 13 & Age <=55)

TFRadj_detallado <- Bx %>%
  left_join(Exp, by = join_by(Year, Age)) %>%
  filter(!is.na(Age)) %>%
  mutate(Fxi = Bxi / Exp) %>%
  group_by(Year, order) %>%
  summarise(
    TFRi = sum(Fxi),
    MABi = sum(Fxi * Age) / TFRi,
    .groups = "drop") %>%
  arrange(order, Year) %>%
  group_by(order) %>%
  mutate(
    ri = (lead(MABi) - lag(MABi)) / 2,
    TFRadji = TFRi / (1 - ri)
  ) %>%
  ungroup() %>%
  mutate(
    TFRi = round(TFRi, 3),
    TFRadji = round(TFRadji, 3)
  )

# Obtener tabla con TFR ajustada total + por orden
TFRadj_total_y_ordenes <- TFRadj_detallado %>%
  group_by(Year) %>%
  summarise(
    TFRadj = if (any(is.na(TFRadji))) NA else round(sum(TFRadji), 3),
    TGF = if (any(is.na(TFRi))) NA else round(sum(TFRi), 3),
    .groups = "drop"
  ) %>%
  left_join(
    TFRadj_detallado %>%
      dplyr::select(Year, order, TFRadji) %>%
      mutate(nombre = paste0("adjTFR", order)) %>%
      pivot_wider(names_from = nombre, values_from = TFRadji),
    by = "Year")


# 9.1. Calcular TFR ajustada por orden ####
TFRadj_detallado <- Bx %>%
  left_join(Exp, by = join_by(Year, Age)) %>%
  filter(!is.na(Age)) %>%
  mutate(Fxi = Bxi / Exp) %>%
  group_by(Year, order) %>%
  summarise(
    TFRi = sum(Fxi),
    MABi = sum(Fxi * Age) / TFRi,
    .groups = "drop"
  ) %>%
  arrange(order, Year) %>%
  group_by(order) %>%
  mutate(
    ri = (lead(MABi) - lag(MABi)) / 2,
    TFRadji = TFRi / (1 - ri)
  ) %>%
  ungroup() %>%
  mutate(
    TFRi = round(TFRi, 3),
    TFRadji = round(TFRadji, 3)
  )

# 9.2. Pivotear TFR ajustada por orden (adjTFR1, adjTFR2, etc.) ####
TFRadj_por_orden <- TFRadj_detallado %>%
  dplyr::select(Year, order, TFRadji) %>%
  mutate(nombre = paste0("adjTFR", order)) %>%
  pivot_wider(names_from = nombre, values_from = TFRadji)

# 9.3. Calcular TFRadj total + suavizados ####
TFRadj_total <- TFRadj_detallado %>%
  group_by(Year) %>%
  summarise(
    TFRadj = if (any(is.na(TFRadji))) NA else round(sum(TFRadji), 3),
    TGF = if (any(is.na(TFRi))) NA else round(sum(TFRi), 3),
    .groups = "drop") %>%
  arrange(Year) %>%
  mutate(
    TFRadj_3ma = round(
      rollapply(TFRadj, width = 3, align = "center", fill = NA,
                FUN = function(x) if (any(is.na(x))) NA else mean(x)),
      3),
    TFRadj_5ma = round(
      rollapply(TFRadj, width = 5, align = "center", fill = NA,
                FUN = function(x) if (any(is.na(x))) NA else mean(x)),
      3)) %>%
  mutate(
    TFRadj_hibrido = case_when(
      Year == max(Year) ~ TFRadj,
      Year == max(Year) - 1 ~ TFRadj_3ma,
      TRUE ~ TFRadj_5ma
    )
  )

# 9.4. Unir TFR total y por orden ####
TFRadj_total_y_ordenes <- TFRadj_total %>%
  left_join(TFRadj_por_orden, by = "Year") 

# 9.5 Estimación de la Edad media a la fecundidad total y por orden  ####
MAB_completo <- Bx %>% 
  left_join(Exp, by = join_by(Year, Age)) %>%
  filter(!is.na(Exp), !is.na(Bxi)) %>%
  mutate(
    Fxi = Bxi / Exp,
    edad_media = Age + 0.5) %>%
  group_by(Year) %>%
  mutate(TFR_total = sum(Fxi, na.rm = TRUE)) %>%
  ungroup() %>%
  # 9.5.1. MAB total
  group_by(Year) %>%
  summarise(
    MAB = round(sum(Fxi * edad_media, na.rm = TRUE) / sum(Fxi, na.rm = TRUE), 2),
    .groups = "drop") %>%
  # 9.5.2. Añadir MAB por orden
  left_join({
    Bx %>%
      left_join(Exp, by = join_by(Year, Age)) %>%
      filter(!is.na(Exp), !is.na(Bxi)) %>%
      mutate(Fxi = Bxi / Exp) %>%
      group_by(Year, order) %>%
      summarise(
        TFRi = sum(Fxi, na.rm = TRUE),
        MABi = sum(Fxi * (Age + 0.5), na.rm = TRUE) / TFRi,
        .groups = "drop") %>%
      mutate(order_group = ifelse(order >= 5, "5p", as.character(order)),
             nombre = paste0("MAB", order_group),
             MABi = round(MABi, 2)) %>%
      dplyr::select(Year, nombre, MABi) %>%
      distinct() %>%
      pivot_wider(names_from = nombre, values_from = MABi)
  }, by = "Year")

TFRadj_total_y_ordenes$Pais = "Argentina"
MAB_completo$Pais = "Argentina"

#### 9.6 Método con Tasas específicas ####
TFRadj <- Bx %>%
  left_join(Exp, by = join_by(Year, Age))

TFRadj2 <- TFRadj %>%
  filter(!is.na(Age)) %>%
  mutate(Fx = Bxi / Exp) %>%
  group_by(Year) %>%
  summarise(
    TFR = sum(Fx, na.rm = TRUE),
    MAC = sum((Age + 0.5) * Fx, na.rm = TRUE) / sum(Fx, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(Year) %>%
  mutate(
    r = (lead(MAC) - lag(MAC)) / 2,
    TFRadj = TFR / (1 - r),
    TFRadj_ma = (lag(TFRadj) + TFRadj + lead(TFRadj)) / 3
  )

#### 9.7 Unir bases de datos ####
metodo1 <- TFRadj_total_y_ordenes %>%
  dplyr::select(Year, TFRadj_orden = TFRadj, TFRadj_3ma, TFRadj_5ma, TFRadj_hibrido, TGF_orden = TGF)

metodo2 <- TFRadj2 %>%
  dplyr::select(Year, TFRadj_total = TFRadj, MAC, TFRadj_ma, TGF_total = TFR)

comparacion <- left_join(metodo1, metodo2, by = "Year")


#### 10. Elaboración del gráfico: Argentina: Comparación entre la TGF observada, ajustada por efecto tempo y suavizada, 1996 – 2022 (Gráfico III.10 del Observatorio Demográfico 2025) ####

g6<- ggplot(comparacion, aes(x = Year)) +
  geom_line(aes(y = TFRadj_orden, color = "TGF ajustada por orden"), linewidth = 1.5) +
  geom_line(aes(y = TGF_total, color = "TGF observada"), linewidth = 2) +
  geom_line(aes(y = TFRadj_total, color = "TGF ajustada por EMM"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_ma, color = "TGF ajustada por EMM suavizada 3 años"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_3ma, color = "TGF ajustada por orden suavizada 3 años"), linewidth = 1.5) +
  geom_line(aes(y = TFRadj_5ma, color = "TGF ajustada por orden suavizada 5 años"), linewidth = 2) +
  geom_hline(yintercept = 2.1, linetype = "dashed", color = "gray5") +
  scale_y_continuous(
    breaks = sort(unique(c(scales::pretty_breaks(n = 7)(range(comparacion$TGF_total, na.rm = TRUE)), 2.1)))) +
  scale_color_manual(values = c(
    "TGF ajustada por orden" = "#1B9E77",
    "TGF observada" = "black",
    "TGF ajustada por EMM" = "lightblue",
    "TGF ajustada por EMM suavizada 3 años" = "#FC8D62",
    "TGF ajustada por orden suavizada 3 años" = "#FF1493",
    "TGF ajustada por orden suavizada 5 años" = "#9932CC")) +
  labs(
    title = "Argentina: comparación entre la tasa global de fecundidad (TGF) observada, ajustada por efecto tempo y suavizada, 1996-2022",
    y = "Tasa Global de Fecundidad",
    x = "Año",
    color = "") +
  wpp_theme() +
  theme( plot.title = element_text(size=17, hjust=0.5),
         legend.position = "bottom",
         legend.title = element_blank(),
         axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  guides(color = guide_legend(nrow = 2, byrow = TRUE))


ggsave(filename = "../Gráficos/III.10_TGFajustadasuavizada_ARG.png",plot = g6, width = 14, height = 10, units = "in", dpi = 300, device = "png")


