# install.packages("MortalityLaws")
library(MortalityLaws)
# =============================================================================
# Extrapolation of Argentine Life Tables using Parametric Mortality Laws
# =============================================================================
# Description:
#   This script reads the official Argentine life tables from INDEC for
#   two reference periods (2000-01 and 2008-10), and extrapolates mortality
#   rates beyond the observed ages using two parametric laws:
#     - Kannisto (1992)
#     - Gompertz (1825)
#
#   The extrapolated life tables are then used as inputs for the r-variable
#   method to estimate centenarian population.
#
# Repository structure expected:
#   inputs/   <- raw INDEC life table Excel files (this script reads from here)
#   outputs/  <- extrapolated life tables in CSV (this script writes here)
#   codes/    <- this script lives here
#
# Inputs (inputs/):
#   - tm_0810_totalpais.xlsx : INDEC life tables 2008-10 (sheet 2=men, sheet 3=women)
#   - tm_0001_totalpais.xlsx : INDEC life tables 2000-01 (sheet 2=men, sheet 3=women)
#
# Outputs (outputs/life_tables/):
#   Kannisto extrapolation:
#     tm_0810_varones_kannisto.csv, tm_0810_mujeres_kannisto.csv
#     tm_0001_varones_kannisto.csv, tm_0001_mujeres_kannisto.csv
#   Gompertz extrapolation (robustness check):
#     tm_0810_varones_gompertz.csv, tm_0810_mujeres_gompertz.csv
#     tm_0001_varones_gompertz.csv, tm_0001_mujeres_gompertz.csv
#
# Execution order:
#   1. codes/extrapolar_tablas_v2.R   <- this script
#   2. codes/r_variable_estimation.py <- reads the CSVs produced here
#
# Dependencies: MortalityLaws, tidyverse, readxl, dplyr, magrittr
# =============================================================================
# install.packages("MortalityLaws")
library(MortalityLaws)
library(tidyverse)
library(readxl)
library(dplyr)
library(magrittr)
# -----------------------------------------------------------------------------
# 0. Paths
# -----------------------------------------------------------------------------
# All paths are relative to the repository root.
# Run this script from the root of the repository (not from codes/).
input_dir  <- file.path("inputs")
output_dir <- file.path("outputs", "life_tables")
# Create output subdirectory if it does not exist
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# -----------------------------------------------------------------------------
# 1. Read official INDEC life tables
# -----------------------------------------------------------------------------
# Life tables 2008-10
tm_0810_varones <- read_excel(file.path(input_dir, "tm_0810_totalpais.xlsx"), sheet = 2)
# =============================================================================
# Extrapolation of Argentine Life Tables using Parametric Mortality Laws
# =============================================================================
# Description:
#   This script reads the official Argentine life tables from INDEC for
#   two reference periods (2000-01 and 2008-10), and extrapolates mortality
#   rates beyond the observed ages using two parametric laws:
#     - Kannisto (1992)
#     - Gompertz (1825)
#
#   The extrapolated life tables are then used as inputs for the r-variable
#   method to estimate centenarian population.
#
# Repository structure expected:
#   inputs/   <- raw INDEC life table Excel files (this script reads from here)
#   outputs/  <- extrapolated life tables in CSV (this script writes here)
#   codes/    <- this script lives here
#
# Inputs (inputs/):
#   - tm_0810_totalpais.xlsx : INDEC life tables 2008-10 (sheet 2=men, sheet 3=women)
#   - tm_0001_totalpais.xlsx : INDEC life tables 2000-01 (sheet 2=men, sheet 3=women)
#
# Outputs (outputs/life_tables/):
#   Kannisto extrapolation:
#     tm_0810_varones_kannisto.csv, tm_0810_mujeres_kannisto.csv
#     tm_0001_varones_kannisto.csv, tm_0001_mujeres_kannisto.csv
#   Gompertz extrapolation (robustness check):
#     tm_0810_varones_gompertz.csv, tm_0810_mujeres_gompertz.csv
#     tm_0001_varones_gompertz.csv, tm_0001_mujeres_gompertz.csv
#
# Execution order:
#   1. codes/extrapolar_tablas_v2.R   <- this script
#   2. codes/r_variable_estimation.py <- reads the CSVs produced here
#
# Dependencies: MortalityLaws, tidyverse, readxl, dplyr, magrittr
# =============================================================================
# install.packages("MortalityLaws")
library(MortalityLaws)
library(tidyverse)
library(readxl)
library(dplyr)
library(magrittr)
# -----------------------------------------------------------------------------
# 0. Paths
# -----------------------------------------------------------------------------
# All paths are relative to the repository root.
# Run this script from codes/
input_dir  <- file.path("..", "inputs")
output_dir <- file.path("..", "outputs", "life_tables")
# Create output subdirectory if it does not exist
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# -----------------------------------------------------------------------------
# 1. Read official INDEC life tables
# -----------------------------------------------------------------------------
# Life tables 2008-10
tm_0810_varones <- read_excel(file.path(input_dir, "tm_0810_totalpais.xlsx"), sheet = 2)
tm_0810_mujeres <- read_excel(file.path(input_dir, "tm_0810_totalpais.xlsx"), sheet = 3)
# Life tables 2000-01
tm_0001_varones <- read_excel(file.path(input_dir, "tm_0001_totalpais.xlsx"), sheet = 2)
tm_0001_mujeres <- read_excel(file.path(input_dir, "tm_0001_totalpais.xlsx"), sheet = 3)
# -----------------------------------------------------------------------------
# 2. Function: extrapolate a life table using a parametric mortality law
# -----------------------------------------------------------------------------
#
# Arguments:
#   tm          : data frame with the observed life table. Expected columns:
#                   Edad (age), nmx (central death rate), ndx (deaths),
#                   nLx (person-years lived), lx (survivors), n (interval width),
#                   Ax (average fraction of interval lived at death).
#   edad_ajuste : numeric vector of length 2 — age range used to fit the law,
#                 e.g. c(70, 95). The model is fit on observed data within
#                 this interval.
#   edad_predic : scalar — age from which the extrapolation takes over.
#                 Ages below this value are kept as observed; from this age
#                 onwards, rates are replaced by model predictions.
#   ley         : character — mortality law to apply. Options: "kannisto"
#                 (logistic, recommended for 80+) or "gompertz" (exponential).
#                 Passed directly to MortalityLaw().
#
# Returns: a list with two elements:
#   [[1]] tm_predic — full life table (observed + extrapolated), single-year ages
#   [[2]] L_predic  — abridged version with nLx summed into 5-year age groups
#                     (0, 5, ..., 120), used as input for the r-variable method.
#
extrapolar_tabla <- function(tm, edad_ajuste = c(70, 90), edad_predic = 100,
ley = "kannisto") {
# --- 2.1 Fit the mortality law on the specified age range -----------------
# MortalityLaw() estimates the model parameters via maximum likelihood.
# Inputs are:
#   x  = midpoint of each 5-year age group (Edad + 2.5)
#   Dx = observed deaths in that age group (ndx)
#   Ex = person-years of exposure, approximated by nLx
ages_fit <- tm$Edad >= edad_ajuste[1] & tm$Edad <= edad_ajuste[2]
ajuste <- MortalityLaw(
x   = tm$Edad[ages_fit] + 2.5,   # midpoint of each 5-year age group
Dx  = tm$ndx[ages_fit],           # deaths
Ex  = tm$nLx[ages_fit],           # exposure (person-years)
law = ley
)
summary(ajuste)
# --- 2.2 Predict mortality rates for ages from edad_predic to 150 ---------
# Single-year age groups from edad_predic onwards.
# nmx is predicted at each age midpoint (age + 0.5).
ages_pred <- edad_predic:150
tm_predic <- data.frame(
Edad = ages_pred,
nmx  = predict(ajuste, ages_pred + 0.5)
)
# nqx: probability of dying in [x, x+1), derived from the central death rate
# under the constant-force-of-mortality assumption: nqx = 1 - exp(-nmx)
tm_predic$nqx <- 1 - exp(-tm_predic$nmx)
# Ax: average fraction of the year lived by those who die in [x, x+1).
# Under constant force: Ax = 1/nmx - 1/(exp(nmx) - 1)
# which simplifies to: Ax = 1 - (nmx - nqx) / (nmx * nqx)
tm_predic$Ax <- 1 - (tm_predic$nmx - tm_predic$nqx) / (tm_predic$nmx * tm_predic$nqx)
# Single-year interval width
tm_predic$n <- rep(1, nrow(tm_predic))
# --- 2.3 Merge observed (below edad_predic) with extrapolated rows --------
tm$n <- as.numeric(tm$n)
tm_predic <- bind_rows(tm[tm$Edad < edad_predic, ], tm_predic)
# --- 2.4 Rebuild all life table columns from nqx --------------------------
# lx: survivors at exact age x, with radix l(0) = 100,000.
# l(x+1) = l(x) * (1 - nqx); equivalently, l(x) = 100,000 * prod_{i<x}(1 - nq_i).
# We keep l(0) from the observed table and recompute the rest.
tm_predic$lx[2:nrow(tm_predic)] <- head(cumprod(1 - tm_predic$nqx) * 100000, -1)
# ndx: deaths between ages x and x+n
tm_predic$ndx <- tm_predic$lx * tm_predic$nqx
# nLx: person-years lived in [x, x+n).
# nLx = l(x+n) * n  +  ndx * Ax
tm_predic$nLx <- lead(tm_predic$lx, n = 1L) * tm_predic$n + tm_predic$ndx * tm_predic$Ax
tm_predic$nLx[is.na(tm_predic$nLx)] <- 0   # terminal open-ended group
# Tx: total person-years lived above exact age x (tail sum of nLx)
tm_predic$Tx <- rev(cumsum(rev(tm_predic$nLx)))
# exo: remaining life expectancy at age x
tm_predic$exo <- tm_predic$Tx / tm_predic$lx
# --- 2.5 Build the 5-year abridged nLx table (for r-variable inputs) ------
# The r-variable method requires nLx aggregated into standard 5-year groups.
# Age groups: [0,5), [5,10), ..., [115,120), [120+].
x    <- seq(from = 0, to = 120, by = 5)
L_ab <- vapply(
x,
function(i) sum(tm_predic$nLx[tm_predic$Edad >= i & tm_predic$Edad < (i + 5)]),
numeric(1)
)
# The last group (120+) is open-ended: sum all remaining person-years
L_ab[length(L_ab)] <- sum(tm_predic$nLx[tm_predic$Edad >= 120])
L_predic <- data.frame(Edad = x, nLx = L_ab)
return(list(tm_predic, L_predic))
}
# -----------------------------------------------------------------------------
# 3. Apply Kannisto extrapolation
# -----------------------------------------------------------------------------
# 2008-10: fit on ages 70–95, predict from 100 onwards.
# 2000-01: fit on ages 70–90, predict from 95 onwards (the INDEC table for
#   this period has reliable data only up to ~90).
#
# Both [[1]] (full single-year table) and [[2]] (5-year abridged nLx) are
# extracted. The abridged tables are saved separately as they are the direct
# inputs to the r-variable estimation script.
# 2008-10 — Kannisto
results_0810_varones_k <- extrapolar_tabla(tm_0810_varones, edad_ajuste = c(70, 95), edad_predic = 100, ley = "kannisto")
results_0810_mujeres_k <- extrapolar_tabla(tm_0810_mujeres, edad_ajuste = c(70, 95), edad_predic = 100, ley = "kannisto")
tm_0810_varones_kannisto <- results_0810_varones_k[[1]]
tm_0810_mujeres_kannisto <- results_0810_mujeres_k[[1]]
nLx_0810_varones_kannisto <- results_0810_varones_k[[2]]
nLx_0810_mujeres_kannisto <- results_0810_mujeres_k[[2]]
# 2000-01 — Kannisto
results_0001_varones_k <- extrapolar_tabla(tm_0001_varones, edad_ajuste = c(70, 90), edad_predic = 95, ley = "kannisto")
results_0001_mujeres_k <- extrapolar_tabla(tm_0001_mujeres, edad_ajuste = c(70, 90), edad_predic = 95, ley = "kannisto")
tm_0001_varones_kannisto <- results_0001_varones_k[[1]]
tm_0001_mujeres_kannisto <- results_0001_mujeres_k[[1]]
nLx_0001_varones_kannisto <- results_0001_varones_k[[2]]
nLx_0001_mujeres_kannisto <- results_0001_mujeres_k[[2]]
# -----------------------------------------------------------------------------
# 4. Diagnostic plot: observed vs. Kannisto-extrapolated log(nmx)
# -----------------------------------------------------------------------------
# Visually verify that the extrapolation follows a plausible trajectory.
# Kannisto should produce a decelerating rise in log mortality at the oldest
# ages — a flatter slope than Gompertz beyond the fit range.
ggplot() +
geom_line(data = tm_0810_varones,
aes(x = Edad, y = log(nmx), color = "INDEC (observed)")) +
geom_line(data = tm_0810_varones_kannisto,
aes(x = Edad, y = log(nmx), color = "INDEC + Kannisto")) +
scale_color_manual(
name   = "",
values = c("INDEC (observed)" = "black", "INDEC + Kannisto" = "#2ca25f")
) +
coord_cartesian(xlim = c(60, 120)) +
labs(
title    = "Log mortality rates: observed vs. Kannisto extrapolation",
subtitle = "Men, Argentina 2008-10",
x        = "Age",
y        = "log(nmx)"
) +
theme_minimal() +
theme(legend.position = "bottom")
# -----------------------------------------------------------------------------
# 5. Export Kannisto extrapolated life tables
# -----------------------------------------------------------------------------
# Full single-year tables (used for diagnostics and documentation)
write.csv(tm_0810_varones_kannisto,  file.path(output_dir, "tm_0810_varones_kannisto.csv"),  row.names = FALSE)
write.csv(tm_0810_mujeres_kannisto,  file.path(output_dir, "tm_0810_mujeres_kannisto.csv"),  row.names = FALSE)
write.csv(tm_0001_varones_kannisto,  file.path(output_dir, "tm_0001_varones_kannisto.csv"),  row.names = FALSE)
write.csv(tm_0001_mujeres_kannisto,  file.path(output_dir, "tm_0001_mujeres_kannisto.csv"),  row.names = FALSE)
# 5-year abridged nLx tables (direct inputs for r_variable_estimation.py)
write.csv(nLx_0810_varones_kannisto, file.path(output_dir, "nLx_0810_varones_kannisto.csv"), row.names = FALSE)
write.csv(nLx_0810_mujeres_kannisto, file.path(output_dir, "nLx_0810_mujeres_kannisto.csv"), row.names = FALSE)
write.csv(nLx_0001_varones_kannisto, file.path(output_dir, "nLx_0001_varones_kannisto.csv"), row.names = FALSE)
write.csv(nLx_0001_mujeres_kannisto, file.path(output_dir, "nLx_0001_mujeres_kannisto.csv"), row.names = FALSE)
# -----------------------------------------------------------------------------
# 6. Apply Gompertz extrapolation (robustness check)
# -----------------------------------------------------------------------------
# Gompertz assumes exponential increase in mortality at old ages, which tends
# to overestimate death rates relative to Kannisto. Used for sensitivity analysis.
# 2008-10 — Gompertz
results_0810_varones_g <- extrapolar_tabla(tm_0810_varones, edad_ajuste = c(70, 95), edad_predic = 100, ley = "gompertz")
results_0810_mujeres_g <- extrapolar_tabla(tm_0810_mujeres, edad_ajuste = c(70, 95), edad_predic = 100, ley = "gompertz")
tm_0810_varones_gompertz  <- results_0810_varones_g[[1]]
tm_0810_mujeres_gompertz  <- results_0810_mujeres_g[[1]]
nLx_0810_varones_gompertz <- results_0810_varones_g[[2]]
nLx_0810_mujeres_gompertz <- results_0810_mujeres_g[[2]]
# 2000-01 — Gompertz
results_0001_varones_g <- extrapolar_tabla(tm_0001_varones, edad_ajuste = c(70, 90), edad_predic = 95, ley = "gompertz")
results_0001_mujeres_g <- extrapolar_tabla(tm_0001_mujeres, edad_ajuste = c(70, 90), edad_predic = 95, ley = "gompertz")
tm_0001_varones_gompertz  <- results_0001_varones_g[[1]]
tm_0001_mujeres_gompertz  <- results_0001_mujeres_g[[1]]
nLx_0001_varones_gompertz <- results_0001_varones_g[[2]]
nLx_0001_mujeres_gompertz <- results_0001_mujeres_g[[2]]
# -----------------------------------------------------------------------------
# 7. Export Gompertz extrapolated life tables
# -----------------------------------------------------------------------------
write.csv(tm_0810_varones_gompertz,  file.path(output_dir, "tm_0810_varones_gompertz.csv"),  row.names = FALSE)
write.csv(tm_0810_mujeres_gompertz,  file.path(output_dir, "tm_0810_mujeres_gompertz.csv"),  row.names = FALSE)
write.csv(tm_0001_varones_gompertz,  file.path(output_dir, "tm_0001_varones_gompertz.csv"),  row.names = FALSE)
write.csv(tm_0001_mujeres_gompertz,  file.path(output_dir, "tm_0001_mujeres_gompertz.csv"),  row.names = FALSE)
write.csv(nLx_0810_varones_gompertz, file.path(output_dir, "nLx_0810_varones_gompertz.csv"), row.names = FALSE)
write.csv(nLx_0810_mujeres_gompertz, file.path(output_dir, "nLx_0810_mujeres_gompertz.csv"), row.names = FALSE)
write.csv(nLx_0001_varones_gompertz, file.path(output_dir, "nLx_0001_varones_gompertz.csv"), row.names = FALSE)
write.csv(nLx_0001_mujeres_gompertz, file.path(output_dir, "nLx_0001_mujeres_gompertz.csv"), row.names = FALSE)
setwd("~/centenarios-argentina")
source("codes/logquad_estimation.R")
source("codes/logquad_estimation.R")
