Proyecciones utilizando el modelo Lee-Carter
Nader Wizani
8 de enero de 2019
Introducción
ilc es un paquete abierto al público para el análisis de datos de mortalidad en un período de edades, que implementa regresiones especializadas y métodos descriptivos para ajustar modelos de Lee-Carter. El objetivo de la modelización de la mortalidad en este paquete es implementar una mejor estructura de modelización, que en sí, mejora el modelo estándar de Lee-Carter basado en el error de la Normal. Consecuentemente, el paquete se aleja de la forma tradicional de descomposición en valores singulares para el ajuste, que asume los residuales Gaussianos, en cambio, implementa una herramienta de regresión basada en procesos de máxima verosimilitud de Poisson. En particular, el paquete utiliza los métodos ilustrados en Renshaw y Haberman (2006), que generaliza la modelización básica de Lee-Carter y extiende el trabajo de Brouhns et al. (2002), para desarrollar un proceso iterativo adaptado para actualizar las estimaciones de parámetros.
Esta metodología de modelización generalizada, está implentada en el software estadístico R en la forma de una serie de funciones comando construidas expresamente para aplicar el ajuste del modelo iterativo mencionado anteriormente. El paquete contiene métodos para analizar 6 distintos tipos de modelos log-lineales en el marco de los GLM, con errores de Poisson que incluye el modelo básico de Lee-Carter. El paquete además, ofrece opciones de gráficos, correción preliminar de datos (incluyendo potenciales datos dispersos).
Datos de Mortalidad
Consideremos una experiencia de mortalidad observada en edades individuales (x) y años calendario (t), dado el aumento a un total de (k*n) celdas de data disponibles, y así podemos estimar la tasa de mortalidad central \((m_{xt})\) y la correspondiente fuerza de mortalidad \((\mu_{xt})\) con la siguiente expresión: \(\bar\mu_{xt} = \bar m_{xt} = \dfrac{y_{xt}}{e_{xt}}\), donde \(y_{xt}\) y \(e_{xt}\) representan el número de muertes y la correspondiente exposición central de cada sub-grupo, respectivamente. Además, para cada combinación de \(x\) y período \(t\), definimos el año de la cohorte \(z = t - x\), representando el año de nacimiento de cada sub-grupo en la data.
Período de años básiico de Lee-Carter
El modelo básico AP (período de años por sus siglas en inglés), fue propuesto en primer lugar por Lee y Carter (1992), y fue introducido como un tipo de modelos de componentes principales de la tasa de mortalidad \((m_{xt})\), dependiendo sólo de factores relacionados con la edad y el tiempo. El modelo es expresado como: \(log(m_{xt}) = \alpha_{x} + \beta_{x}k_{t} + \epsilon_{xt}\), donde los parámetros son interpretados de la siguiente manera:
- \(\alpha_{x}\) representa un patrón constante de una edad específica de mortalidad;
- \(k_{t}\) mide la tendencia en la mortalidad en el tiempo;
- \(\beta_{x}\) mide las desviaciones del cambio en la mortalidad en el tiempo de una edad específica;
- \(\epsilon_{xt}\) son efectos aleatorios Gaussianos distribuidos \(N(0, \sigma^{2})\) por la edad y el tiempo.
Debido a la construcción bilineal multiplicativa (\(\beta_{x}k_{t}\)) en la ecuación, hay un claro problema de identificabilidad que es resuelto tradicionalmente asegurándose que estos parámetros satisfagan dos restricciones:
- \(\sum_{x} \beta_{x} = 1\)
- \(\sum_{t = t_{1}}^{t_{n}} k{t} = 0\)
Entonces, el modelo estándar de Lee-Carter puede ser estimado usando la decomposición en valores singulares (SVD por sus siglas en inglés), que lleva al siguiente estimador de los efectos de edades específicas:
- \(\bar a_{x} = \dfrac{1}{n} \sum_{t = t_{1}}^{t_{n}} log(\bar m_{xt})\),
que minimiza el término de la suma de errores al cuadraro \((S = \sum_{xt} \epsilon^{2}_{xt})\). Lee y Carter también realizan algunos ajustes a las estimaciones \(\bar k_{t}\) para asegurarse de que en cada año, el total de muertes predichas por el modelo sea igual al total de muertes observadas \(\sum_{x} y_{xt}\).
Familia Generalizada de Modelos Lee-Carter
En un estudio más reciente, el modelo básico ha sido extendido para incluir un término bilineal adicional, conteniendo un segundo efecto de período (como en Renshaw y Haberman, 2003) o un efecto de una cohorte (como en Renshaw y Haberman, 2006). En particular, este acercamiento tuvo luz debido a los patrones de mortalidad a principios del siglo 20 en Gales e Inglaterra. Además, el modelo básico de Lee-Carter puede ser transformado en un esquema más general para poder analizar la relación entre el tiempo y la edad, y el impacto conjunto que tienen sobre las tasas de mortalidad. En este paquete, se sigue el esquema de la modelización APC (age-period-cohort) y el ajuste del modelo con una metodología propuesta por Renshaw y Haberman (2006), que especifica la fuerza de la mortalidad con una estructura generalizada escrita de la forma:
- \(\mu_{xt} = e^{\alpha_{x} + \beta^{0}_{xt}l_{t}-x + \beta^{1}_{x}k_{t}}\), donde \(\alpha_{x}\) mapea el perfil principal de mortalidad, \(l_{t}\) y \(k_{t}\) representan efectos de cohorte y período, respectivamente, donde \(\beta^{0}\) y \(\beta^{1}\) miden las correspondientes interacciones con la edad.
Notamos que el modeo representa una familia de seis modelos generalizados no-lineales de la estructura Lee-Carter con la función log-link. Las sub-categorías del modelo general pueden ser definidas ajustando independientemente los parámetros \(\beta^{0}\) y \(\beta^{1}\) a lo siguiente:
- desconocido (debe ser estimado);
- =1, (ajustado);
- =0, (vacío).
La función principal del paquete ilc las 6 sub-estructuras del modelo anterior, haciendo uso bien sea de los errores Gaussianos o de Poisson. La estimación general de esta clase de estructuras de modelo procede con la misma técnica de minimización iterativa descrita anteriormente.
Proyecciones
Las proyecciones de tasas de mortalidad en el caso de la familia de modelos Lee-Carter está basada en predicciones de series de tiempo de parámetros dependientes de años calendario (\(l_{t}-x, k_{t}\)). Lo podemos escribir de la siguiente forma:
- \(\mu_{x,t_{n}+s} = e^{\bar\alpha_{x} + \bar\beta^{0}_{xt}l_{t_{n}+s-x} + \bar\beta^{1}_{x}k_{t_{n+s}}}\), s > 0,
donde \(l_{t_{n}+s-x}\) y \(k_{t_{n+s}}\), representan las proyeccciones de los efectos de la cohorte y períodos, respectivamente. Observamos que en caso de las proyecciones de efectos de cohorte, los valores proyectados se revierten a los parámetros ajustados cuando el horizonte de la proyección es menor al rango de datos disponibles. Este método de proyección nos permite generar futuros valores promedio y evaluar la variabilidad futura de la tasa central de mortalidad. Además, la variabilidad de las predicciones pueden ser aplicadas para medir la incertidumbre en el riesgo de la longevidad.
El tipo más común de metódos de extrapolación de series de tiempo aplicados en el esquema de Lee-Carter es el ARIMA (Modelo autorregresivo integrado de media móvil) univariado, que son caracterizados por tres parámetros (p;d;q). El tipo de modelo ARIMA que se utiliza depende del perfil de los parámetros ajustados en el rango de datos disponibles. En la mayoría de los modelos Lee-Carter el camino aleatorio con deriva (0,1,0) es la opción más común para efectos de período (\(k_{t}\)), esto puede ser expresado como:
- \(k_{t} = k_{t-1} + d + e_{t}\),
donde \(d\) mide la deriva en forma de las desviaciones promedio anuales y \(e_{t}\) representa el ruido blanco en el proceso estocástico.
Metodología de Ajuste del Modelo
Como se mencionó anteriormente, la metodología de ajuste que implementa el paquete está basada en un algoritmo iterativo que minimiza la funcion de desviación. El paquete hace usa de un proceso cíclico de actualización de los parámetros estimados hasta que la diferencia mínima entre la verosimilitud del modelo ajustado y la verosimilitud del modelo saturado es alcanzada. Además, el mecanismo de actualización para un parámetro \(\theta\) es dado por el método de minimización de Newton-Raphson aplicado a la función de desviación, que puede ser expresada como:
- \(u(\bar\theta) = \bar\theta - \dfrac{\dfrac{\delta D}{\delta \theta}}{\dfrac{\delta^{2}D}{\delta \theta^{2}}}\)
Si observamos la estructura de la función de Poisson con errores, obtenemos que:
- \(\dfrac{\delta D}{\delta \theta} = \sum 2wa(\bar y - y)\), donde
- \(\bar y' = \dfrac{\delta \bar y}{\delta \theta} = a \bar y\)
Haciendo uso de las notaciones simplificadas anteriores, podemos expresar la segunda derivada parcial de la funcion de desviación como:
- \(\dfrac{\delta^{2}D}{\delta \theta^{2}} = \sum 2wa\bar y' = \sum 2wa^{2}\bar y'\)
Sustituyendo en la ecuación principal las dos igualdades obtenidas:
- \(u(\bar\theta) = \bar\theta - \dfrac{\sum 2wa(\bar y - y)}{\sum 2wa^{2}\bar y'} = \bar\theta + \dfrac{\sum 2wa(\bar y - y)}{\sum 2wa^{2}\bar y'}\)
Observamos que el paquete ilc implementa los algoritmos actualizadores correspondientes a los modelos con estructuras de errores de Poisson y Gaussianos.
Aplicación de los Modelos Generalizados de Lee-Carter en R con ilc
A continuación se presentarán las funciones más importantes del paquete ilc para ajustar y analizar modelos de mortalidad dependientes de la edad y el tiempo. Se extrajo la data de la experiencia de la mortalidad de hombres en Italia hasta el año 2014 de la HMD (Human Mortality Database). La mayoría de las funciones del paquete son complementadas con los paquetes demography y forecast.
Preparando los datos
Para ajustar un modelo generalizado de Lee-Carter, los datos de mortalidad deben ser de clase demogdata, formato del paquete demography. Afortunadamente, se extrajo los datos mediante una función que tuvo que ser modificada debido a que la función original del paquete no ha sido actualizada correctamente, y esta extrae los datos con formato demogdata. Al final encontrarán como anexo el código para poder extraer data de la HMD (Human Mortality Database).
italyDemo<-hmd.mx2(country="ITA", username="naderwizani@gmail.com", password="21759759", label="Italy")
par(mfrow=c(1,3))
plot(italyDemo,series="male",datatype="rate", main="Tasas de mortalidad hombres")
plot(italyDemo,series="female",datatype="rate", main="Tasas de mortalidad mujeres")
plot(italyDemo,"total",datatype="rate", main="Tasas de mortalidad totales")
Acá observamos la tasa de mortalidad por edades de hombres, mujeres, y el total. Lógicamente, a medida que la edad aumenta, la tasa de mortalidad es mayor tanto para hombres como para mujeres.
par(mfrow=c(1,3))
plot(italyDemo,series="male",datatype="rate", plot.type="time", main="Tasas de mortalidad hombres",xlab="Años")
plot(italyDemo,series="female",datatype="rate", plot.type="time", main="Tasas de mortalidad mujeres",xlab="Años")
plot(italyDemo,series="total",datatype="rate", plot.type="time", main="Tasas de mortalidad totales",xlab="Años")
Estas tres gráficas nos muestran cómo ha evolucionado la tasa de la mortalidad en el tiempo. Se observa que esta ha tenido un decrecimiento a medida que avanzan los años, esto debido a las mejores condiciones sanitarias, calidad de vida y avances en la tecnología en general.
El paquete nos ofrece hacer gráficas con leyendas y ajustar los años que queremos graficar:
par(mfrow=c(1,1))
plot_dd(italyDemo, xlim=c(40, 110), lpar=list(x.int=-0.2, y.int=0.9, cex=0.85))
par(mfrow=c(1,2))
plot_dd(italyDemo, year=1985:1995, transf=F)
plot_dd(italyDemo, year=1995:1997, transf=F, lty=1:3, col=1:3)
También podemos extraer el número de muertes de nuestros datos y graficarlos de la siguiente forma:
tmp_d <- extract.deaths(italyDemo, ages = 55:100, fill = "perks")
tmp_d$type <- "mortality"
par(mfrow=c(1,1))
plot_dd(tmp_d, year=1995:2003, transf=F, lty=1:8)
La función dd.rfp puede tomar un objeto de clase demogdata y de tipo ‘mortality’ y ajustar el logaritmo de las tasas de mortalidad observadas mediante un vector de efectos aditivos distribuidos Poisson con medias predeterminadas. Por ejemplo, con los datos que tenemos, podemos producir aleatoriamentes unos datos de mortalidad estratificados de la siguiente forma:
rfpitaly <- dd.rfp(italyDemo, rfp=c(0.5, 1.2, -0.7, 2.5))
rfpitaly ## RESUMEN## Multidimensional Mortality data for: Italy [female] [male] [total]
## Across covariates:
## years: 1872 - 2014
## ages: 0 - 110
## X: base, a, b, c, dPodemos graficar los perfiles de las exposiciones centrales y tasas de mortalidad logaritmicas de la siguiente forma:
## Exposiciones centrales
par(mfrow=c(1,2))
matplot(rfpitaly$age, rfpitaly$pop[,,1], type='l', xlab='Age',ylab='Ec', main='Nivel Base')
matplot(rfpitaly$age, rfpitaly$pop[,,2], type='l', xlab='Age',ylab='Ec', main='Nivel 1') 
## Tasas de mortalidad logaritmicas
par(mfrow=c(1,2))
matplot(rfpitaly$age, log(rfpitaly$mu[,,1]), type='l', xlab='Age',ylab='log(mu)', main='Nivel Base')
matplot(rfpitaly$age, log(rfpitaly$mu[,,2]), type='l', xlab='Age',ylab='log(mu)', main='Nivel 1') 
Ajustando el Modelo Generalizado de Lee-Carter y realizando Proyecciones
La función lca.rh es la función principal del paquete ilc, desarrollada para ajustar cualquiera de las seis variantes de las estructuras de modelos Lee-Carter (incluyendo el modelo base de Lee-Carter) utilizando el método iterativo de calibración. Veamos los argumentos:
args(lca.rh)## function (dat, year = dat$year, age = dat$age, series = 1, max.age = 100,
## dec.conv = 6, clip = 3, error = c("poisson", "gaussian"),
## model = c("m", "h0", "h1", "h2", "ac", "lc"), restype = c("logrates",
## "rates", "deaths", "deviance"), scale = F, interpolate = F,
## verbose = T, spar = NULL)
## NULLEstimación del modelo base de Lee-Carter (con errores de Poisson)
Aplicaremos el modelo base de Lee-Carter a la data que tenemos de Italia para los hombres.
modelo1 <- lca.rh(italyDemo, mod= 'lc', interpolate = T, verbose = F, series = "male")## Original sample: Mortality data for Italy
## Series: female male total
## Years: 1872 - 2014
## Ages: 0 - 110
## Applied sample: Mortality data for Italy
## Series: female male total
## Years: 1872 - 2014
## Ages: 0 - 100
##
## Fitting model: [ LC = a(x)+b1(x)*k(t) ]
## - with poisson error structure and with deaths as weights -
##
## Iterations finished in: 158 stepsmodelo1##
## ------------------------------------------------------------
## Iterative Lee-Carter Family Regression:
## Fitted Model: LC = a(x)+b1(x)*k(t)
## ------------------------------------------------------------
## Call:
## lca.rh(dat = italyDemo, series = "male", model = "lc", interpolate = T,
## verbose = F)
## Error Structure: poisson
## Data Source: Italy [male] over
## calendar years: (1872 - 2014) and ages: (0 - 100)
## Deviance convergence in: 158 iterations
## dev dev.c df df.c
## 1 Mean deviance base 157.163 df base 14100
## 2 Mean deviance total 815.785 df tot 14241par(mfrow=c(1,1))
plot(modelo1)
coef(modelo1) ##coeficientes del modelo## $ax
## 0 1 2 3 4 5
## -2.7521111 -4.4034962 -5.3982359 -5.9378847 -6.2658520 -6.5208848
## 6 7 8 9 10 11
## -6.6710019 -6.8039555 -6.9093709 -6.9788447 -7.0120240 -7.0033802
## 12 13 14 15 16 17
## -6.9694619 -6.9018931 -6.7252862 -6.5570825 -6.3628111 -6.2330180
## 18 19 20 21 22 23
## -6.0281115 -5.9365961 -5.8593570 -5.8159399 -5.7988894 -5.7918974
## 24 25 26 27 28 29
## -5.8014264 -5.8132436 -5.8222641 -5.8271073 -5.8157945 -5.8116391
## 30 31 32 33 34 35
## -5.7667582 -5.7692898 -5.7421232 -5.7138662 -5.6969797 -5.6657631
## 36 37 38 39 40 41
## -5.6313526 -5.5940729 -5.5410830 -5.4896487 -5.4140972 -5.3683037
## 42 43 44 45 46 47
## -5.3006613 -5.2386117 -5.1712865 -5.1034623 -5.0372589 -4.9699686
## 48 49 50 51 52 53
## -4.8863758 -4.8070321 -4.7079055 -4.6465355 -4.5629455 -4.4901351
## 54 55 56 57 58 59
## -4.4135185 -4.3457796 -4.2702319 -4.1910679 -4.1046704 -4.0129442
## 60 61 62 63 64 65
## -3.9012403 -3.8257775 -3.7286752 -3.6409477 -3.5550620 -3.4674873
## 66 67 68 69 70 71
## -3.3856245 -3.2966263 -3.2033959 -3.1109177 -3.0006088 -2.9223717
## 72 73 74 75 76 77
## -2.8197484 -2.7290455 -2.6348641 -2.5258407 -2.4308968 -2.3372639
## 78 79 80 81 82 83
## -2.2364200 -2.1366541 -2.0543923 -1.9722414 -1.8690416 -1.7721116
## 84 85 86 87 88 89
## -1.6789512 -1.5993543 -1.5171527 -1.4383923 -1.3576851 -1.2766662
## 90 91 92 93 94 95
## -1.2069870 -1.1413220 -1.0639729 -1.0073775 -0.9459818 -0.8959533
## 96 97 98 99 100+
## -0.8204520 -0.7655373 -0.6965941 -0.6269499 -0.5203481
##
## $bx1
## 0 1 2 3 4
## 0.0155369837 0.0236255664 0.0295523795 0.0306618907 0.0296276098
## 5 6 7 8 9
## 0.0283244684 0.0259804915 0.0230070265 0.0201990186 0.0184512164
## 10 11 12 13 14
## 0.0186890012 0.0186511785 0.0186838410 0.0178138412 0.0161747596
## 15 16 17 18 19
## 0.0147921927 0.0132673882 0.0129365337 0.0163903634 0.0183305412
## 20 21 22 23 24
## 0.0193893990 0.0196869398 0.0196461613 0.0192201402 0.0185544179
## 25 26 27 28 29
## 0.0176227026 0.0167620352 0.0160490016 0.0154037068 0.0149276735
## 30 31 32 33 34
## 0.0142994722 0.0139318126 0.0135501406 0.0130868694 0.0127283106
## 35 36 37 38 39
## 0.0122990012 0.0118471418 0.0114620812 0.0111037787 0.0107918302
## 40 41 42 43 44
## 0.0105568774 0.0099906817 0.0097574557 0.0093074601 0.0089188436
## 45 46 47 48 49
## 0.0084633501 0.0079631752 0.0076529911 0.0073559569 0.0070655348
## 50 51 52 53 54
## 0.0070241952 0.0065970072 0.0063545128 0.0060579132 0.0057162025
## 55 56 57 58 59
## 0.0053319059 0.0050806187 0.0048417120 0.0047877503 0.0047708709
## 60 61 62 63 64
## 0.0049202062 0.0046890529 0.0047613226 0.0047576680 0.0047131867
## 65 66 67 68 69
## 0.0046306944 0.0045252907 0.0044608738 0.0045100007 0.0044891816
## 70 71 72 73 74
## 0.0046681876 0.0044811306 0.0045104675 0.0044595896 0.0043939983
## 75 76 77 78 79
## 0.0045000395 0.0044177006 0.0043806250 0.0043548604 0.0042853926
## 80 81 82 83 84
## 0.0039998347 0.0038542099 0.0038277629 0.0037314309 0.0036439662
## 85 86 87 88 89
## 0.0034525886 0.0032620800 0.0030797534 0.0028765429 0.0027339332
## 90 91 92 93 94
## 0.0025451803 0.0023031984 0.0021612520 0.0019329941 0.0017354899
## 95 96 97 98 99
## 0.0015162594 0.0014445880 0.0011998533 0.0011602024 0.0011655827
## 100
## 0.0008099039
##
## $kt
## Time Series:
## Start = 1872
## End = 2014
## Frequency = 1
## 1872 1873 1874 1875 1876
## 97.4999655 93.0971546 92.6700393 93.6931957 89.4694046
## 1877 1878 1879 1880 1881
## 86.7265316 87.7593546 88.2087119 90.8018058 88.1250712
## 1882 1883 1884 1885 1886
## 88.5540972 86.2337416 83.1216216 82.9516542 86.2783649
## 1887 1888 1889 1890 1891
## 84.6902696 82.8330893 78.5613625 79.3639657 79.9863533
## 1892 1893 1894 1895 1896
## 78.5018276 76.7846465 76.4598854 77.4010875 75.0614512
## 1897 1898 1899 1900 1901
## 69.0962835 71.2798503 67.9020730 72.1192956 67.8627431
## 1902 1903 1904 1905 1906
## 69.7394510 69.0000042 66.0469113 67.1126862 64.3964981
## 1907 1908 1909 1910 1911
## 62.8758719 68.5021282 65.1443162 60.5638528 64.2629835
## 1912 1913 1914 1915 1916
## 55.2846462 56.2284299 52.0181569 76.5798940 87.2213392
## 1917 1918 1919 1920 1921
## 94.4478049 109.9645334 76.3103966 61.6020873 54.7132741
## 1922 1923 1924 1925 1926
## 51.6066074 48.1138105 47.0801430 48.5432915 49.2700753
## 1927 1928 1929 1930 1931
## 44.2545672 44.1703195 45.0790133 36.8693125 37.3063113
## 1932 1933 1934 1935 1936
## 37.7082091 32.8099385 31.1836363 33.2961984 31.6988883
## 1937 1938 1939 1940 1941
## 35.9000318 33.7342992 28.8045248 32.8504788 44.0283266
## 1942 1943 1944 1945 1946
## 54.5715656 66.0025170 54.7792900 45.4472664 24.4116508
## 1947 1948 1949 1950 1951
## 15.8920731 4.6324796 0.9131809 -7.0073415 -5.5385396
## 1952 1953 1954 1955 1956
## -8.7334754 -12.3594016 -18.7565550 -20.4295189 -16.7058917
## 1957 1958 1959 1960 1961
## -17.1880472 -23.3537182 -26.1022175 -24.8194324 -27.9203891
## 1962 1963 1964 1965 1966
## -23.9170322 -24.1317019 -31.0717797 -30.3766235 -35.5621027
## 1967 1968 1969 1970 1971
## -35.8730961 -33.8869163 -33.2481536 -39.0268559 -40.2852951
## 1972 1973 1974 1975 1976
## -41.8743628 -41.3559091 -46.5195388 -44.4482784 -46.7308076
## 1977 1978 1979 1980 1981
## -48.6964189 -50.7478916 -52.9224921 -51.4014758 -55.2839167
## 1982 1983 1984 1985 1986
## -60.0355064 -57.5752938 -65.4308017 -66.2584253 -69.7517306
## 1987 1988 1989 1990 1991
## -74.1838380 -76.2236723 -80.2720629 -80.4319743 -79.5214491
## 1992 1993 1994 1995 1996
## -84.7470575 -89.1968812 -91.7454799 -94.4642112 -100.1829216
## 1997 1998 1999 2000 2001
## -105.4322911 -106.5604482 -113.3569117 -120.6646340 -125.6952442
## 2002 2003 2004 2005 2006
## -130.4744906 -133.6769587 -148.1932467 -149.9401084 -158.3093287
## 2007 2008 2009 2010 2011
## -161.9590665 -166.1986594 -170.3776232 -178.3300009 -182.1992061
## 2012 2013 2014
## -183.8458642 -193.9178693 -200.6397358
##
## attr(,"class")
## [1] "coef"El modelo realizó 158 iteraciones hasta encontrar el valor mínimo de la suma de errores al cuadrado. Finalmente, podemos hacer proyecciones de la mortalidad y la correspondiente esperanza de vida futura, basada en el modelo ajustado. El paquete ilc se complementa con el paquete forecast a la hora de hacer las proyecciones para predecir los valores del parámetro de tendencia \(k_{t}\) usando un ARIMA(0,1,0). Si queremos hacer una proyección a 8 años lo hacemos de la siguiente forma:
forecast <- forecast(modelo1, h = 8, jump = 'fit', level = 90, shift = F)
plot_dd(forecast, xlim = c(45,100), lpar = list(x.int=-0.2, y.int=0.9, cex=0.95))
Además, el objeto forecast también contiene la proyección media de la esperanza de vida y su intervalo de confianza del 90%. Podemos hacer varias cosas con este objeto, como por ejemplo:
forecast$e0## Time Series:
## Start = 2015
## End = 2022
## Frequency = 1
## [1] 80.72140 80.81675 80.91141 81.00541 81.09876 81.19146 81.28354 81.37500lifeexpectancy <- life.expectancy(forecast, age=60)
flc.plot(modelo1, at = 60, h = 30, level = 90)
La función flc.plot nos muestra la gráfica de la proyección de la esperanza de vida en un horizonte dado, con un intervalo de confianza y a una edad en específico. En este caso, nuestro horizonte es 30 años, con un 90% de confianza y para una persona de edad 60.
Estimación del modelo APC de Lee-Carter (Age-Period-Cohort)
En esta aplicación utilizaremos un rango de edad restringida y el tipo de respuesta del modelo como ‘desviación’ para los residuales. Según Renshaw y Haberman (2006), el tipo preferido de residuales para realizar diagnósticos en el modelo, son los residuales estandarizados de desviación.
modelo_apc <- lca.rh(italyDemo, age = 55:95, res = 'dev', dec = 3, verb = F, series = "male")## Original sample: Mortality data for Italy
## Series: female male total
## Years: 1872 - 2014
## Ages: 0 - 110
## Applied sample: Mortality data for Italy
## Series: female male total
## Years: 1872 - 2014
## Ages: 55 - 95
##
## Fitting model: [ M = a(x)+b0(x)*i(t-x)+b1(x)*k(t) ]
## - with poisson error structure and with deaths as weights -
##
## Iterations finished in: 10440 stepsmodelo_apc##
## ------------------------------------------------------------
## Iterative Lee-Carter Family Regression:
## Fitted Model: M = a(x)+b0(x)*i(t-x)+b1(x)*k(t)
## ------------------------------------------------------------
## Call:
## lca.rh(dat = italyDemo, age = 55:95, series = "male", dec.conv = 3,
## restype = "dev", verbose = F)
## Error Structure: poisson
## Data Source: Italy [male] over
## calendar years: (1872 - 2014) and ages: (55 - 95)
## Deviance convergence in: 10440 iterations
## dev dev.c df df.c
## 1 Mean deviance base 9.169 df base 5445
## 2 Mean deviance total 101.920 df tot 5781plot(modelo_apc)
Realizamos las proyecciones tal y como lo hicimos anteriormente:
forecast_apc <- forecast(modelo_apc, h = 8, jump = 'fit', level = 90, shift = F)
plot_dd(forecast_apc)
Comparemos ahora las dos proyecciones que realizamos:
Si notamos las adeades de 55 a 95 años de la primera gráfica (modelo de Lee-Carter básico), observamos que es muy parecida a la del modelo APC que utilizamos con el tipo de respuesta ‘desviación’. El modelo que se recomendaría utilizar es el modelo APC, debido a que cuando utilizamos la desviación de los residuales, esto hace que los valores ajustados tengan mejor precisión (Renshaw y Haberman, 2006).
En conclusión, el paquete ilc nos ofrece una gran variedad de alternativas e innovaciones al modelo básico de Lee-Carter, la iteración normalmente suele ser un mecanismo eficaz a la hora de calibrar un modelo, y el paquete innova el modelo Lee-Carter básico, además de ofrecer gráficas especiales para graficar datos provenientes de tablas de mortalidad.
ANEXO
Anexo el código para arreglar la función de extracción de data de la Human Mortality Database
hmd.mx2 <- function(country, username, password, label=country)
{
path <- paste("https://www.mortality.org/hmd/", country, "/STATS/", "Mx_1x1.txt", sep = "")
userpwd <- paste(username, ":", password, sep = "")
txt <- RCurl::getURL(path, userpwd = userpwd)
con <- textConnection(txt)
mx <- try(utils::read.table(con, skip = 2, header = TRUE, na.strings = "."),TRUE)
close(con)
if(class(mx)=="try-error")
stop("Connection error at www.mortality.org. Please check username, password and country label.")
path <- paste("https://www.mortality.org/hmd/", country, "/STATS/", "Exposures_1x1.txt", sep = "")
userpwd <- paste(username, ":", password, sep = "")
txt <- RCurl::getURL(path, userpwd = userpwd)
con <- textConnection(txt)
pop <- try(utils::read.table(con, skip = 2, header = TRUE, na.strings = "."),TRUE)
close(con)
if(class(pop)=="try-error")
stop("Exposures file not found at www.mortality.org")
obj <- list(type="mortality",label=label,lambda=0)
obj$year <- sort(unique(mx[, 1]))
#obj$year <- ts(obj$year, start=min(obj$year))
n <- length(obj$year)
m <- length(unique(mx[, 2]))
obj$age <- mx[1:m, 2]
mnames <- names(mx)[-c(1, 2)]
n.mort <- length(mnames)
obj$rate <- obj$pop <- list()
for (i in 1:n.mort)
{
obj$rate[[i]] <- matrix(mx[, i + 2], nrow = m, ncol = n)
obj$rate[[i]][obj$rate[[i]] < 0] <- NA
obj$pop[[i]] <- matrix(pop[, i + 2], nrow = m, ncol = n)
obj$pop[[i]][obj$pop[[i]] < 0] <- NA
dimnames(obj$rate[[i]]) <- dimnames(obj$pop[[i]]) <- list(obj$age, obj$year)
}
names(obj$pop) = names(obj$rate) <- tolower(mnames)
obj$age <- as.numeric(as.character(obj$age))
if (is.na(obj$age[m]))
obj$age[m] <- 2 * obj$age[m - 1] - obj$age[m - 2]
return(structure(obj, class = "demogdata"))
}