Introduction to demofit

library(demofit)

1. Introduction

The demofit package provides a unified framework for fitting parametric mortality curves and forecasting mortality rates with stochastic models.

The package implements multiple mortality curves and stochastic mortality models with a consistent interface and common S3 methods.

2. Preparing mortality data

For mortality curves, demofit accepts mortality rates in vector form.

For stochastic mortality models, demofit accepts mortality data in matrix form. Rows represent calendar years and columns represent ages.

The function mortalitymatrix() converts mortality data in data frame or demogdata format into the required numeric matrix format.

# generate synthetic mortality rates in data frame
a <- c(-4.8499,-4.7676,-4.6719,-4.5722,-4.4847,-4.3841,-4.2813,-4.1863,-4.0861,-3.9962,
-3.8885,-3.7896,-3.6853,-3.5737,-3.4728,-3.3718,-3.2586,-3.1474,-3.0371,-2.9206,
-2.7998,-2.6845,-2.5653,-2.4581,-2.3367,-2.2159,-2.1017,-1.9941,-1.8821, -1.7697)
b <- c(0.0283,0.0321,0.0335,0.0336,0.0341,0.0358,0.0368,0.0403,0.0392,0.0395,
0.0396,0.0399,0.0397,0.0386,0.039,0.0375,0.0367,0.0368,0.035,0.0354,
0.0336,0.0323,0.0313,0.0295,0.0282,0.0265,0.024,0.0226,0.0219,0.0183)
k <- c(12.11,10.69,11.18,9.64,9.35,8.21,6.89,5.74,4.56,3.6,
3.27,2.04,1.11,-0.44,-1.05,-1.03,-1.84,-2.9,-4.03,-4.12,
-5.18,-5.64,-6,-6.51,-6.91,-6.9,-8.32,-8.53,-9.69,-9.31)
set.seed(123)
Raw <- exp(outer(k,b)+matrix(a,nrow=30,ncol=30,byrow=TRUE)+rnorm(900,0,0.035))
Data <- as.data.frame(Raw)
# convert data frame into numeric matrix
M <- mortalitymatrix(Data)

The resulting matrix can be passed directly to the stochastic mortality model functions.

3. Mortality curve fitting

The function MC() provides a unified interface for fitting parametric mortality curves.

# generate synthetic mortality rates in vector form
x <- 60:89
set.seed(123); m <- 0.0000082*exp(0.10771*c(60:89)+rnorm(30,0,0.1))
# fit Gompertz curve to mortality rates
fit <- MC(x=x,m=m,curve="gompertz")
# print estimated parameters
coef(fit)
##            B            C 
## 9.803215e-06 1.052498e-01
# plot fitted Gompertz curve
plot(fit)

The fitted Gompertz curve closely follows the underlying mortality pattern.

# generate synthetic mortality rates in vector form
x <- 60:89
set.seed(123); m <- 0.005+0.0000082*exp(0.10771*c(60:89)+rnorm(30,0,0.1))
# fit Makeham curve to mortality rates
fit <- MC(x=x,m=m,curve="makeham")
# print estimated parameters
coef(fit)
##            A            B            C 
## 5.339115e-03 8.531247e-06 1.068591e-01
# plot fitted Makeham curve
plot(fit)

The fitted Makeham curve closely follows the underlying mortality pattern.

Different mortality curves can be fitted using the same interface by specifying the curve argument.

4. Mortality forecasting

The function FCS() provides a unified interface for forecasting mortality rates with stochastic models.

# generate synthetic mortality rates in matrix form
x <- 60:89
a <- c(-4.8499,-4.7676,-4.6719,-4.5722,-4.4847,-4.3841,-4.2813,-4.1863,-4.0861,-3.9962,
-3.8885,-3.7896,-3.6853,-3.5737,-3.4728,-3.3718,-3.2586,-3.1474,-3.0371,-2.9206,
-2.7998,-2.6845,-2.5653,-2.4581,-2.3367,-2.2159,-2.1017,-1.9941,-1.8821, -1.7697)
b <- c(0.0283,0.0321,0.0335,0.0336,0.0341,0.0358,0.0368,0.0403,0.0392,0.0395,
0.0396,0.0399,0.0397,0.0386,0.039,0.0375,0.0367,0.0368,0.035,0.0354,
0.0336,0.0323,0.0313,0.0295,0.0282,0.0265,0.024,0.0226,0.0219,0.0183)
k <- c(12.11,10.69,11.18,9.64,9.35,8.21,6.89,5.74,4.56,3.6,
3.27,2.04,1.11,-0.44,-1.05,-1.03,-1.84,-2.9,-4.03,-4.12,
-5.18,-5.64,-6,-6.51,-6.91,-6.9,-8.32,-8.53,-9.69,-9.31)
set.seed(123)
M <- exp(outer(k,b)+matrix(a,nrow=30,ncol=30,byrow=TRUE)+rnorm(900,0,0.035))
# fit Lee-Carter model to mortality rates and smooth forecasted rates by Makeham curve
fit <- FCS(x=x,M=M,model="LC",curve="makeham",h=30,jumpoff=2)
# print estimated parameters
coef(fit)
## $alpha
##  [1] -4.851558 -4.761369 -4.671056 -4.575497 -4.491137 -4.378732 -4.280798
##  [8] -4.189451 -4.084336 -3.984429 -3.886095 -3.793503 -3.684144 -3.582582
## [15] -3.471694 -3.361932 -3.257354 -3.142917 -3.034809 -2.927308 -2.806363
## [22] -2.688123 -2.573719 -2.455479 -2.325837 -2.212565 -2.100037 -2.004359
## [29] -1.873074 -1.762622
## 
## $beta
##  [1] 0.02950412 0.03241615 0.03260187 0.03536899 0.03473039 0.03589613
##  [7] 0.03637652 0.04058335 0.03869840 0.03862533 0.03786024 0.04002463
## [13] 0.04050987 0.03798978 0.04087331 0.03772439 0.03861324 0.03617972
## [19] 0.03529944 0.03536673 0.03326351 0.03089149 0.03328026 0.02890509
## [25] 0.02585061 0.02693123 0.02238475 0.02194416 0.02222336 0.01908294
## 
## $kappa
##  [1] 12.2305818 10.6042761 11.5356417  9.3733652  9.3130432  8.3209665
##  [7]  7.0833041  5.6544742  4.6653718  3.3634376  3.4474689  1.8589739
## [13]  1.1647736 -0.3772328 -1.1142009 -1.1582084 -1.8767680 -3.2091821
## [19] -4.0391981 -4.1294601 -5.1228184 -5.5461830 -6.1201753 -6.4149479
## [25] -6.9697116 -6.8832021 -8.1439493 -8.6889046 -9.6572739 -9.1642622
# plot fitting and forecasting results
plot(fit)

The function CBDS() fits and forecasts mortality rates using the CBD model specifically.

# generate synthetic mortality rates in matrix form
x <- 60:89
k1 <- -2.97-0.0245*(0:29)
k2 <- 0.101+0.000345*(0:29)
set.seed(123)
M <- exp(matrix(k1,nrow=30,ncol=30,byrow=FALSE)+outer(k2,(x-mean(x)))+rnorm(900,0,0.035))
# fit CBD model to mortality rates and smooth forecasted rates by Makeham curve
fit <- CBDS(x=x,M=M,curve="makeham",h=30,jumpoff=2)
# print estimated parameters
coef(fit)
## $kappa1
##  [1] -2.966968 -2.995883 -3.006593 -3.051805 -3.068170 -3.086121 -3.109082
##  [8] -3.142045 -3.162883 -3.197979 -3.208270 -3.245587 -3.261703 -3.285007
## [15] -3.313835 -3.339804 -3.363685 -3.398987 -3.411299 -3.433770 -3.458551
## [22] -3.481077 -3.512675 -3.528915 -3.559368 -3.581160 -3.601869 -3.635392
## [29] -3.653522 -3.674440
## 
## $kappa2
##  [1] 0.1002979 0.1016178 0.1012783 0.1023293 0.1016872 0.1028631 0.1026888
##  [8] 0.1042037 0.1036287 0.1036629 0.1043472 0.1048233 0.1054915 0.1050585
## [15] 0.1069324 0.1052939 0.1063750 0.1068979 0.1060689 0.1089243 0.1071847
## [22] 0.1080932 0.1096365 0.1089425 0.1098004 0.1103041 0.1100144 0.1104469
## [29] 0.1110232 0.1113381
# plot fitting and forecasting results
plot(fit)

The two fitted models forecast mortality rates for the next 30 years, anchored to observed mortality rates in the last data year.

The stochastic mortality models can be fitted using either the unified interface or individual model-specific functions.

5. Advanced features

The function ENI() generates ensemble interval forecast of mortality rates.

# generate synthetic mortality rates in matrix form
x <- 60:69
a <- c(-4.8499,-4.7676,-4.6719,-4.5722,-4.4847,-4.3841,-4.2813,-4.1863,-4.0861,-3.9962)
b <- c(0.0801,0.0909,0.0948,0.0951,0.0965,0.1014,0.1042,0.1141,0.1110,0.1118)
k <- c(12.11,10.69,11.18,9.64,9.35,8.21,6.89,5.74,4.56,3.60,
3.27,2.04,1.11,-0.44,-1.05,-1.03,-1.84,-2.90,-4.03,-4.12,
-5.18,-5.64,-6.00,-6.51,-6.91,-6.90,-8.32,-8.53,-9.69,-9.31)
set.seed(123)
M <- exp(outer(k,b)+matrix(a,nrow=30,ncol=10,byrow=TRUE)+rnorm(300,0,0.035))
# fit Lee-Carter, Renshaw-Haberman, age-period-cohort, and CBD models to mortality rates
fit1 <- LCS(x=x,M=M,curve="makeham",h=30,jumpoff=2)
fit2 <- RHS(x=x,M=M,curve="makeham",h=30,jumpoff=2)
fit3 <- APCS(x=x,M=M,curve="makeham",h=30,jumpoff=2)
fit4 <- CBDS(x=x,M=M,curve="makeham",h=30,jumpoff=2)
# generate ensemble interval forecast of mortality rates
fit <- ENI(fit1,fit2,fit3,fit4)
# plot ensemble results
plot(fit)

The function CFMS() fits and forecasts mortality rates of two populations using the common factor model.

# generate synthetic mortality rates of two populations in matrix form
x <- 60:89
a1 <- c(-5.18,-5.12,-4.98,-4.92,-4.82,-4.73,-4.66,-4.53,-4.45,-4.35,
-4.26,-4.17,-4.05,-3.95,-3.84,-3.73,-3.65,-3.52,-3.40,-3.29,
-3.14,-3.02,-2.88,-2.76,-2.64,-2.49,-2.37,-2.25,-2.12,-2.00)
a2 <- c(-4.78,-4.68,-4.57,-4.49,-4.39,-4.29,-4.19,-4.10,-4.00,-3.89,
-3.80,-3.69,-3.60,-3.49,-3.39,-3.29,-3.17,-3.07,-2.96,-2.85,
-2.71,-2.62,-2.49,-2.37,-2.26,-2.14,-2.04,-1.91,-1.82,-1.72)
B <- c(0.0381,0.0340,0.0420,0.0389,0.0423,0.0414,0.0406,0.0393,0.0415,0.0400,
0.0411,0.0362,0.0387,0.0381,0.0384,0.0385,0.0356,0.0314,0.0317,0.0337,
0.0316,0.0298,0.0284,0.0270,0.0248,0.0262,0.0205,0.0215,0.0142,0.0145)
K <- c(9.66,9.89,10.66,9.83,9.52,7.39,7.64,6.36,2.32,4.18,
2.91,-0.61,0.28,-0.38,-1.79,-3.34,-1.74,-3.50,-4.28,-4.77,
-4.98,-7.13,-5.09,-6.41,-5.56,-5.65,-6.12,-5.64,-7.35,-6.28)
b1 <- c(0.0012,-0.0033,0.0523,0.0161,0.0529,0.0220,0.0312,0.0437,0.0709,0.0444,
0.0398,0.0361,0.0403,0.0396,0.0506,0.0315,0.0428,0.0261,0.0384,0.0388,
0.0300,0.0269,0.0275,0.0256,0.0239,0.0421,0.0314,0.0284,0.0174,0.0314)
k1 <- c(-1.24,-1.38,-3.48,-2.51,-1.32,-1.90,-3.42,-0.94,0.24,-0.48,
-0.26,2.70,1.39,-0.46,1.74,2.53,0.90,1.43,0.76,2.48,
0.74,2.32,0.42,1.69,-0.64,1.30,0.19,-0.69,-1.11,-1.01)
b2 <- c(-0.0014,0.0272,0.0083,0.0273,0.0209,0.0253,0.0144,0.0333,0.0460,0.0439,
0.0439,0.0674,0.0331,0.0443,0.0312,0.0240,0.0570,0.0312,0.0403,0.0376,
0.0500,0.0289,0.0466,0.0418,0.0349,0.0149,0.0366,0.0178,0.0361,0.0372)
k2 <- c(2.35,0.62,-0.38,0.12,0.00,0.80,-1.39,0.38,2.47,0.40,
0.76,3.06,1.42,-0.73,0.79,1.94,0.12,0.60,-0.43,0.29,
0.17,0.98,-1.01,-0.13,-2.46,-1.24,-1.65,-2.48,-2.32,-3.06)
set.seed(123)
M1 <- exp(outer(k1,b1)+outer(K,B)+matrix(a1,nrow=30,ncol=30,byrow=TRUE)+rnorm(900,0,0.07))
M2 <- exp(outer(k2,b2)+outer(K,B)+matrix(a2,nrow=30,ncol=30,byrow=TRUE)+rnorm(900,0,0.07))
# fit common factor model to mortality rates and smooth forecasted rates by Makeham curve
fit <- CFMS(x=x,M1=M1,M2=M2,curve="makeham",h=30,jumpoff=2)
# print estimated parameters
coef(fit)
## $alpha1
##  [1] -5.183272 -5.107493 -4.978280 -4.926552 -4.832840 -4.719220 -4.658954
##  [8] -4.536264 -4.446442 -4.326420 -4.255149 -4.177768 -4.047650 -3.967727
## [15] -3.837754 -3.710224 -3.647474 -3.510997 -3.395386 -3.303383 -3.153093
## [22] -3.027214 -2.896807 -2.754728 -2.618246 -2.483309 -2.366654 -2.270498
## [29] -2.101929 -1.985832
## 
## $alpha2
##  [1] -4.757664 -4.695171 -4.579278 -4.480465 -4.381814 -4.283737 -4.194741
##  [8] -4.093199 -4.007107 -3.887608 -3.817219 -3.682068 -3.601522 -3.486368
## [15] -3.402580 -3.288480 -3.166292 -3.067173 -2.944959 -2.846177 -2.702362
## [22] -2.608148 -2.494114 -2.364741 -2.247449 -2.135749 -2.042906 -1.919588
## [29] -1.813195 -1.715405
## 
## $B
##  [1] 0.03725862 0.03205288 0.04116460 0.03793507 0.04302767 0.04118460
##  [7] 0.04012376 0.03984389 0.04183725 0.03888903 0.04261718 0.03436459
## [13] 0.04077983 0.03847557 0.04434605 0.03942397 0.03862935 0.03095216
## [19] 0.03053224 0.02904238 0.03301609 0.02974130 0.03059406 0.02318300
## [25] 0.02429096 0.02503781 0.01955552 0.02117862 0.01503119 0.01589075
## 
## $K
##  [1] 10.3992626  9.5956005 11.5641104 10.4974471  9.9041023  7.7638157
##  [7]  7.5377179  6.0258711  2.7177135  3.1073968  4.1318204 -0.9711885
## [13]  0.2525531 -0.4277492 -1.7497507 -2.8831810 -1.6697076 -4.2475739
## [19] -3.9644383 -4.2705815 -4.3867120 -6.3033589 -6.3745216 -6.7367001
## [25] -6.4704871 -5.9948342 -6.7553510 -6.2579858 -7.9720398 -6.0612500
## 
## $beta1
##  [1] -0.009900327 -0.014478270  0.053632523  0.007436720  0.048632697
##  [6]  0.021060265  0.033469774  0.041973534  0.075229853  0.052415982
## [11]  0.060662150  0.033327084  0.034589692  0.039336346  0.064633964
## [16]  0.035861990  0.048438923  0.027051593  0.038055718  0.016843407
## [21]  0.040934898  0.033940238  0.023453597  0.009931392  0.032128780
## [26]  0.039750679  0.034584042  0.027785852  0.022005787  0.027211119
## 
## $kappa1
##  [1] -1.7113062 -1.3923368 -3.6400712 -3.5650618 -1.7437159 -2.0315534
##  [7] -2.9719948 -0.8627244  0.1496509  0.2106459 -0.9777645  2.5295305
## [13]  1.5573598 -0.6265752  1.5669343  1.8173200  0.7888943  1.7692606
## [19]  0.3546920  1.9883226  0.3425878  1.7261205  1.4410241  2.0591940
## [25]  0.1367043  1.6573201  1.0685545 -0.3835887 -0.3271824 -0.9302410
## 
## $beta2
##  [1] -0.024684636  0.032907912 -0.003942276  0.021151359  0.013811921
##  [6]  0.028599604  0.015788987  0.035239642  0.046158141  0.046592993
## [11]  0.038736752  0.084228320  0.054691271  0.048800715  0.007631359
## [16]  0.027181650  0.048246364  0.059863260  0.055446827  0.032687367
## [21]  0.046859865  0.045157175  0.037310969  0.043323854  0.042090376
## [26]  0.001484669  0.016063012  0.020225882  0.039406992  0.038939677
## 
## $kappa2
##  [1]  1.402438e+00  1.521767e+00 -1.042946e+00 -2.771619e-01 -1.824959e-01
##  [6]  3.279962e-01 -1.414447e+00  4.489068e-01  1.853083e+00  1.462887e+00
## [11] -1.044947e-01  2.690480e+00  1.528225e+00 -3.663089e-01  2.080287e-01
## [16]  2.077982e+00  4.572661e-05  1.893909e+00 -5.876895e-01 -4.226409e-01
## [21] -5.747285e-02 -5.512255e-02  1.263464e-01 -4.307235e-02 -1.812958e+00
## [26] -1.429187e+00 -1.562694e+00 -1.619645e+00 -1.831580e+00 -2.732180e+00
# plot fitting (residuals) and forecasting results
plot(fit,which=2)

The package also supports other advanced features such as mortality ensemble and modelling sparse data.