# ============================================================
#  CÁLCULO DE e_xy A PARTIR DE l_x  (Moivre, x=40, y=50, ω=100)
# ============================================================

x     <- 40
y     <- 50
omega <- 100
l0    <- 100

# Función de cohorte de Moivre: l(z) = l0 * (1 - z/omega)
l <- function(z) l0 * (1 - z/omega)

# Horizonte del grupo
t_max <- omega - max(x, y)   # = 50

# ============================================================
#  TABLA CON l_x, l_y, l_{x+t}, l_{y+t}
# ============================================================
t          <- 0:t_max
edad_x     <- x + t
edad_y     <- y + t

l_x        <- l(x)            # l_40 = 60 (constante en la tabla)
l_y        <- l(y)            # l_50 = 50 (constante en la tabla)
l_xt       <- l(edad_x)       # l_{40+t}
l_yt       <- l(edad_y)       # l_{50+t}

# Supervivencias residuales a partir de los cocientes de l
p_x        <- l_xt / l_x
p_y        <- l_yt / l_y
p_xy       <- p_x * p_y

tabla <- data.frame(
  t       = t,
  edad_x  = edad_x,
  edad_y  = edad_y,
  l_x     = l_x,
  l_y     = l_y,
  l_xt    = l_xt,
  l_yt    = l_yt,
  p_x     = round(p_x, 6),
  p_y     = round(p_y, 6),
  p_xy    = round(p_xy, 6)
)

print(tabla)

# ============================================================
#  e_xy POR REGLA DEL TRAPECIO
# ============================================================
e_xy_trap <- sum( (p_xy[-length(p_xy)] + p_xy[-1]) / 2 )
cat("\ne_xy (trapecio) =", round(e_xy_trap, 6), "\n")

# Suma directa (lo que hacías antes): NO es e_xy
suma_directa <- sum(p_xy)
cat("Suma directa de p_xy(t) =", round(suma_directa, 6), "\n")

# Suma corregida: quitar medio valor en los extremos
e_xy_corr <- suma_directa - 0.5 * p_xy[1] - 0.5 * p_xy[length(p_xy)]
cat("e_xy (suma corregida) =", round(e_xy_corr, 6), "\n")

# ============================================================
#  FÓRMULA CERRADA PARA COMPROBAR
# ============================================================
dx <- omega - x
dy <- omega - y
e_xy_exacto <- t_max -
  (1/dx + 1/dy) * t_max^2 / 2 +
  (1/(dx*dy)) * t_max^3 / 3
cat("e_xy (fórmula cerrada) =", round(e_xy_exacto, 6), "\n")

# ============================================================
#  ESPERANZAS INDIVIDUALES Y EXTINCIÓN
# ============================================================
e_x <- (omega - x) / 2
e_y <- (omega - y) / 2
cat("e_x =", e_x, "  e_y =", e_y, "\n")
cat("e_xy_barra = e_x + e_y - e_xy =", round(e_x + e_y - e_xy_exacto, 6), "\n")
# install.packages("openxlsx")
library(openxlsx)
write.xlsx(tabla, "tabla_exy.xlsx")