Fitting two parameters Feller process to mortality data.

Model we want to fit:

dm = (am)dt + ssqrt(m)dW

Mortality data

A two parameter Feller process was calibrated to the mortality rates provided by the Human Mortality Database www.mortality.org

The period considered in this study is from year 1985 to 2014. The tables used were those of an Italian male with age 50 in 1985 (79 in 2014).

The 30 year period considered was splitted into a 20 years training period (from 1985 to 2004) and a 10 year test period (from 2005 to 2014). The training data is used to fit the model while the test data is used to measure its predictive performance.

## $estimation
##       Year Age      mx      qx  ax    lx   dx    Lx      Tx    ex
## 12594 1985  50 0.00568 0.00566 0.5 92856  526 92593 2372731 25.55
## 12706 1986  51 0.00596 0.00594 0.5 92527  549 92253 2306244 24.92
## 12818 1987  52 0.00684 0.00682 0.5 92134  628 91820 2252285 24.45
## 12930 1988  53 0.00688 0.00686 0.5 91660  629 91346 2180202 23.79
## 13042 1989  54 0.00734 0.00732 0.5 91119  667 90785 2122297 23.29
## 13154 1990  55 0.00796 0.00793 0.5 90376  716 90018 2040315 22.58
## 13266 1991  56 0.00901 0.00897 0.5 89464  802 89063 1953427 21.83
## 13378 1992  57 0.00935 0.00930 0.5 88899  827 88486 1899469 21.37
## 13490 1993  58 0.01018 0.01013 0.5 88648  898 88199 1832408 20.67
## 13602 1994  59 0.01073 0.01067 0.5 87876  938 87407 1759825 20.03
## 13714 1995  60 0.01168 0.01161 0.5 87176 1012 86670 1693766 19.43
## 13826 1996  61 0.01252 0.01244 0.5 86681 1078 86142 1635711 18.87
## 13938 1997  62 0.01349 0.01340 0.5 86154 1154 85577 1569527 18.22
## 14050 1998  63 0.01449 0.01439 0.5 85343 1228 84729 1487494 17.43
## 14162 1999  64 0.01585 0.01573 0.5 84627 1331 83961 1434837 16.95
## 14274 2000  65 0.01614 0.01601 0.5 83972 1344 83300 1387851 16.53
## 14386 2001  66 0.01730 0.01715 0.5 83001 1423 82289 1328543 16.01
## 14498 2002  67 0.01831 0.01814 0.5 82151 1490 81406 1266799 15.42
## 14610 2003  68 0.01999 0.01979 0.5 81231 1608 80427 1188049 14.63
## 14722 2004  69 0.02153 0.02130 0.5 80655 1718 79796 1169861 14.50
## 
## $test
##       Year Age      mx      qx  ax    lx   dx    Lx      Tx    ex
## 14834 2005  70 0.02285 0.02259 0.5 79476 1796 78578 1092671 13.75
## 14946 2006  71 0.02352 0.02324 0.5 78366 1821 77456 1047623 13.37
## 15058 2007  72 0.02604 0.02570 0.5 77114 1982 76123  980748 12.72
## 15170 2008  73 0.02830 0.02790 0.5 75684 2112 74628  916786 12.11
## 15282 2009  74 0.03071 0.03025 0.5 74007 2238 72888  854051 11.54
## 15394 2010  75 0.03255 0.03203 0.5 72645 2327 71481  803337 11.06
## 15506 2011  76 0.03570 0.03507 0.5 70818 2484 69576  742319 10.48
## 15618 2012  77 0.03977 0.03900 0.5 68616 2676 67278  674416  9.83
## 15730 2013  78 0.04300 0.04210 0.5 67069 2823 65657  635021  9.47
## 15842 2014  79 0.04739 0.04630 0.5 64968 3008 63464  586162  9.02

Fitting procedure.

The survival function implied by the Feller process describing the force of mortality was fitted to the observed survival function calculated according to the mortality rates (column mx) reported in these tables.

The fitting was computed minimizing the root mean squared error by means of the NLopt library (https://nlopt.readthedocs.io/en/latest/), specifically the R interface to the library provided by the nloptr package (https://CRAN.R-project.org/package=nloptr).

The algorithm used was NLOPT_LN_COBYLA without gradient information and the stop threshold “xtol_rel” was chosen so that to avoid overfitting while keeping the error of prediction over the test period under control. In fact, if the threshold is too low, the algorithm will seek for an almost perfect match to the training data, setting the estimated volatility to almost zero thus reducing the Feller stochastic model to a deterministic Gompertz model. Instead if the threshold is set too high the predictive performance will suffer and we will have a poor prediction of the mortality rates over the test period.

## Post check, all these must be positive

##  [1] 0.008618330 0.008051815 0.007522539 0.007028055 0.006566075
##  [6] 0.006134462 0.005731221 0.005354487 0.005002517 0.004673683
## [11] 0.004366464 0.004079440 0.003811284 0.003560754 0.003326692
## [16] 0.003108016 0.002903715 0.002712843 0.002534518

## Stop threshold:  0.001

## The estimated parameters  were: a= 0.06767605  s= 0.004642473

Plot of simulated paths

Quality of prediction

The quality of the prediction is measured by means of the mean relative errors calculated over the test period.

Specifically year by year the mean relative errors are calculated by taking the mean of the relative errors of each simulated path, where the relative error is calculated with the following formula:

(estimated mx - actual mx) / actual mx

The yearly mean relative errors will be positive if the estimated mortality rate is greater than the actual.

## Mean relative errors by year

##  [1] -0.06976440 -0.03491483 -0.06936170 -0.08626497 -0.10038807
##  [6] -0.09289752 -0.11716088 -0.15421615 -0.16446189 -0.19070978

## Mean of 10 years mean relative errors

## [1] -0.108014

Expected residual lifetime

The expected residual lifetime is calculated using the doubly stochastic model and Monte Carlo simulation where the paths are generated from the calibrated Feller model.

## The expected residual life time at  50 :  32.4351

## Standard error:  0.4509209