Mortality Prediction using Survival Energy Models with Functional Data Analysis

The Survival Energy Model (SEM), as originally introduced by Shimizu et al. (2020), is designed to characterize human bioenergetics by employing diffusion processes or inverse Gaussian processes. While parametric models have been employed to articulate the SEM, they exhibit inherent sensitivity in their parameters and hyperparameters, which in turn introduces issues of instability in estimation and prediction. In this paper, we demonstrate that the utilization of functional data analysis techniques for nonparametric estimation and prediction of critical functions within the SEM leads to a substantial enhancement in prediction performance.


Introduction
Survival Energy Model (SEM) is a stochastic model proposed by Shimizu et al. [11] for cohortwise prediction of mortality.It assumes that human death is caused by the disappearance of hypothetical 'survival energy' and expresses the dynamics as a stochastic process.
On a stochastic basis (Ω, F , P; F) with the right-continuous filtration F := (F t ) t ≥0 , a survival energy process of a certain cohort, say c, is defined as a càdlàg process X c = (X c t ) t ≥0 .Then the time of death, say τ, is defined as the first-hitting time of X to zero: which is a F t -stopping time.Thus, the probability of death up to time t > 0 of the corresponding cohort is given by q c (t ) := P(τ c ≤ t ), (1.1) which is called the mortality function at cohort c.The first SEM by Shimizu et al. [11] is an inhomogeneous diffusion process, called ID-SEM: where x c > 0 is the 'initial energy', U c , V c : [0, ∞) → R are deterministic functions, and W is a Wiener process.The model features intuitive modeling by interpreting U c (t )dt as representing the average drift at time t , and V 2 c (t )dt as representing its dispersion, which can cause sudden death if it values larger.Note that this SEM has a continuous path, and X c τ c = 0 almost surely.On the other hand, the second SEM by Shimizu et al. [12] is based on an inverse-Gaussian process, say Y c = (Y c t ) t ≥0 : where the probability density of Y c t is given by where σ c > 0 is a constant and Λ c is a nonnegative increasing function with Λ c (0) = 0.This SEM is called IG-SEM, and the path may have jumps such that X c τ c < 0. In particular, if Λ c (t ) = t , then X c is a spectrally negative Lévy process with infinite activity jumps.Although the interpretation of Λ c is not as intuitively easy as in the case of ID-SEM, Λ c (t ) approximately represents the degree of energy change at time t , and the larger Λ c , the higher the mortality rate in the model.
The utilization of the SEM offers numerous computational advantages in actuarial science, particularly in the context of cohort-specific premium calculations, sensitivity analyses for premiums with respect to the change of mortality, and longevity assessments.This is welldocumented in the works of Shimizu et al. [13] and Shirai and Shimizu [14].
In employing these models, a suitable choice of U c ,V c , and Λ c is important.Occasionally, the predictive performance can be influenced by the hyperparameters embedded within these functions.While a well-constructed parametric model has demonstrated the potential to outperform not only the classical Lee-Carter model [6] but also other established cohort models such as the CBD model by Cairns et al. [1] and the RH model by Renshaw and Haberman [8,9] -as detailed in Shimizu et al. [11,12] numerical comparisons -the converse outcome can easily transpire if executed improperly.The selection of such parametric models requires a universal measure, like an information criterion.However, constructing such a criterion poses challenges since X c represents a fictitious quantity that was originally unobservable.To overcome this difficulty, embracing the contemporary paradigm of Functional Data Analysis (FDA) holds great promise.
In this manuscript, we aim to estimate these functions and generate forecasts for prospective cohorts through the application of the FDA.First, we employ nonparametric estimation of the functions U c ,V c , and Λ c , using the least-squares method at specific time points t = t 1 , . . ., t d .Plotting these estimates with respect to t , achieved through appropriate smoothing techniques, furnishes us with the functional data for U c ,V c and Λ c .This process is undertaken for each cohort, facilitating our comprehension of how these functional data evolve with cohort progression.Subsequently, by applying the FDA's guidelines to these developments, we can predict the functions U c ,V c , and Λ c for forthcoming cohorts.These predictions are then complemented by our original refinements, culminating in the computation of a predictive mortality function within each SEM.
Our research shows that using SEM predictions with the FDA is better than using parametric frameworks.Importantly, the FDA technique we use is a widely accepted and easy-to-use method available in the R programming package.
In Section 2, we explain the specific steps for estimating and predicting using the FDA for some critical functions, say 'key' functions, that represent each SEM.It is important to note that the predictions from the 'key' function determine the projected mortality function.We also provide an overview of Functional Principal Component Analysis (FPCA) in Subsection 2.2, as it constitutes a primary tool in our forecasting procedure.Section 3 is dedicated to the analysis of real data, sourced from the Human Mortality Database (HMD) [5], which is a prominent open-source repository.Our analysis conclusively demonstrates the superior performance of our proposed methodology in contrast to parametric alternatives.The paper is closed with concluding remarks, presented in Section 4, in which we synthesize our findings and offer insights into the broader implications of our research.
2 Prediction procedure with FDA

Functional datatization of 'key' functions in SEM
Our initial undertaking involves the process of "datatization" of the 'key' functions, which constitute unobservable entities within the SEM.The fundamental concept is to estimate these functions at discrete time points for each cohort and subsequently apply a basis expansion to transform them into (continuous) functional data.
Recall the explicit expressions for mortality functions given in (1.1) in ID-and IG-SEM.

ID-SEM:
The mortality function of (1.2) for cohort c is given as follows (Shimizu et al. [11]).
where Φ is the cumulative distribution function of the standard normal distribution, M c (t ) Due to the last assumption, it suffices to know the negative constant κ c and, e.g., the function M c (or S c ) to have q I D c .Therefore, we consider M c a 'key' function in ID-SEM.IG-SEM: The mortality function of (1.3) for cohort c is given as follows (Shimizu et al. [12]).
In this model, we need to know the constant σ c and the function Λ c , which is the 'key' function in IG-SEM.
Remark 2.1.In this paper, we adopt a fixed approach with respect to the parameters κ c in the ID-SEM and σ c in the IG-SEM, setting them to specific values that are determined through empirical selection.These values are referenced in (3.1) within the data analysis in Section 3. The empirical selection of these values is informed by the outcomes of parametric inference, with the primary objective being the attainment of stability in the estimation of M c or Λ c .While the determination of κ c and σ c remains an important consideration, we do not discuss it further in this paper, but it is worth noting that even if these values are chosen in an ad hoc manner, the predictive efficacy of the mortality function remains good, provided that the 'key' functions can be accurately estimated.
Our primary objective is to acquire the "functional data" pertaining to the 'key' functions.We denote an empirical mortality function by q D AT A c (t ) (t = S, S + 1, . .., w ; c = c 1 , c 2 , . . ., c m ).These empirical functions are derived from data in the Human Mortality Database, as detailed in Shimizu et al. [11].Here, w denotes the final age recorded in the Human Mortality Database, and the value of S is determined on a case-by-case basis by the observer.
In the sequel, we will explain the procedure to obtain functional data for M c in the case of ID-SEM.The one for Λ c in the case of IG-SEM is the same.Moreover, in practice, we use the conditional version of the mortality function following the manner of Shimizu et al. [11]: e.g., in the case of ID-SEM, and its empirical version Remark 2.2.We assume we have mortality data for individuals born in cohort c m until they reach the age of w years.In other words, we have access to mortality data extending up to the calendar year c m+w .Consequently, for cohort c m+L , we have mortality data available up to the age of w − L years, and our objective for this specific cohort is the prediction of the mortality function q c (t |S) for t > (w − L) ∨ S

Smoothing for discretely estimated 'key' functions:
• Estimate M c (t ) for each t = S, S + 1, . .., w ; c = c 1 , c 2 , . . ., c m by the least squares method: Note that the error on the right-hand side can be zero.In this stage, our focus is on determining the value of M at the specific time point t where q I D c (t ; M |S) ≈ q D AT A c (t |S).In essence, we are seeking discrete observations of the latent function M at the moment when the identified quantile function closely approximates the data quantile function.
• For each cohort c = c 1 , c 2 , . . ., c m , we treat the values {M c (t ) | t = S, S + 1, . . ., w } as discrete data sampled from the trajectory of (M c (t )) S≤t ≤w , and proceed to interpolate these values using continuous functions h l : R+ → R, which form an orthogonal basis in the space L 2 ([S, w ]).In practical terms, we truncate the infinite series, limiting it to a suitably "large" value of L: Remark 2.3.Regarding the orthonormal basis {h l } l =1,2,... , the choice of the Fourier series expansion is typically suitable for handling periodic data, whereas B -splines are preferred for non-periodic data.In our research, we use the B -spline expansion method.
The coefficient α c,l in each summand is estimated by the least squares method: In this way, we have m-discrete obervations as functional data for the 'key' function M c in ID-SEM, and we obtain a basis expnsion of centerd data where with

Functional principle component analysis
Functional Principal Component Analysis (FPCA) is an extension of classical Principal Component Analysis (PCA) to handle random elements that take on functional values.Consider a stochastic process Y c = (Y c (t )) t ∈[0,T ] with continuous paths, where T > 0. These processes are indexed by c = 1, 2, . . ., m, and their values reside in a Hilbert space L 2 ([0, T ]).We assume that the random sequence Y 1 , Y 2 , . . ., Y m comes from weakly stationary functionvalued time series, with a mean function µ(t ) = E[Y 1 (t )].In this context, the covariance kernel, denoted as K (t , s) := Cov(Y c (t ), Y c (s)), is invariant with respect to c.It has a spectral-type decomposition, which can be characterized using Mercer's Theorem (see, e.g., Riesz and Sz.-Nagy [10]): where λ j is the j th eigenvalue and e j is the corresponding eigenfunction in L 2 ([0, T ]) space.Then the Karhunen-Lòeve expansion allows us to decompose the random function Y c (•) as where c and e j : which corresponds to the principal component score in the classical (finite-dimensional) PCA.
According to this manner, we can assume that where L is an integer suitably chosen and ǫ c is a noise process depending on the index c.Our aim is to find an approximation of the expansion (2.6).First, we shall estimate the covariance kernel K by a sample version: and we would like to find the eigenvalues λ j 's and eigenfunctions e j 's of K : This practical procedure is seen in Ramsay and Silverman [7], Section 8.4.2, or a more general argument is found in Hyndman and Ullah [4].With the expansion of our data using an orthonormal system, as detailed in (2.4) and (2.5), we can now consider our data in a centered form, expressed as follows: where Following the manner of Ramsay and Silverman [7], the eigenfunctions e j 's of the kernel function K has the following basis expansion: Then, the estimated principal component score is given by As a consequence, we have an approximation of (2.6) as follows: for L ∈ N large enough, Z c j e j (t ).
For the theoretical details for FPCA and the inference, see also Horváth and Kokoszka [2], or Hsing and Eubank [3].

Prediction of the key function for future cohorts
We implement the procedure outlined in the preceding section, applying it to the functional data Y c = M c and Y c = Λ c .These datasets take the form: as in (2.4) and (2.5), respectively.We explain mainly the procedures of estimation and prediction for M c .It is important to note that the procedures for Λ c follows precisely the same way.
Let us start by recalling that we have the smoothly estimated functional data: Step 1: Approximating the Karhunen-Lòeve expansion.
Following the steps in the previous subsection, we have the eigenfunctions e 1 j 's and the corresponding eigenvalues λ 1 j 's of Without loss of generality, we may assume that Z (1)  c j e 1 j (t ), where Z (1)  c j = T 0 M * c (t ) e 1 j (t ) dt is the j th FPC-score.The number L is numerically chosen as L = rank K 1 , where K 1 = ( K 1 (t , s)) t ,s=S,S+1,...,w , and we can reduce it so that the contribution rate of the score exceeds a certain threshold ϑ ∈ (0, 1) (usually close to 1): Then the estimated (approximated) M c is given by Z (1)  c j e 1 j , c = c 1 , . . ., c m .
Following the identical procedure used for M c (t ; α c ), we arrive at an approximation for Λ c as follows: Z (2)  c j e 2 j (t ), c = c 1 , . . ., c m . where The selection of K 2ϑ is determined in a manner similar to that outlined in (2.7).
Step 2: Prediction of mortality function.
We now possess a set of time series data for the principal component scores, represented as { Z (1)  c 1 j , Z (1) c 2 j , . . ., Z (1) c m j }.Considering this dataset as historical data, our objective is to make predictions for { Z (1)  c j | c > c m } for each j = 1, . . ., K 1ϑ .To accomplish this, we will employ time series models and forecast future values based on the distribution of these scores.This approach aligns with the methodology described in Hyndman and Ullah [4].
Remark 2.4.Subsequently, we will employ the ARIMA model to fit the data and utilize the 95% prediction intervals to forecast the SEM mortality function for future cohorts.The process of implementing ARIMA predictions is straightforward, particularly when using the R package forecast.
Once we obtain the predictor { Z (1)  c j | c > c m } for each j = 1, . . ., K 1ϑ , we have the (conditional) mortality functions of future cohorts c > c m : Regarding IG-SEM, we take the same procedures as above.We employ ARIMA models to forecast the future values of { Z (2)  c j | c > c m } for each j = 1, . . ., K 2ϑ and obtain

Step 3: Modification of predicted mortality function.
This modification step is a specialized technique unique to the SEM approach and deviates from conventional forecasting procedures.In this regard, it represents a novel contribution to the field.
When performing our estimations, we utilize (empirical) mortality data up to the age of w for the c m -cohort.Therefore, when we seek to estimate mortality for a cohort c (> c m ), individuals born in that cohort have already reached the age of c m + w − c, and we have access to mortality data up to this age.As a result, we can draw upon information from actual mortality data up to the age of c m +w −c for the cohort c under consideration to enhance our predictions.
In the previous step, where we predicted { Z (k)  c j | c > c m } using the ARIMA model, we also obtained the corresponding prediction intervals.These intervals, denoted as I c,k, j δ , represent the δ-prediction intervals: For instance, with δ = 0.95, the true score Z c k j falls into the I c,k, j δ with a probability of δ.This assumption allows us to fine-tune the fitting of the mortality function, as exemplified below for the case where k = 1: ( Z (1)  c1 , . . ., Z (1)   cK 1ϑ ) = arg min where arg min Z c is taken over all Z (1)  c satisfying that .
Thus, we have the modified version of the conditional mortality function: with the modified 'key' function Z (1)  c j e 1 j (t ).
In the same way, we also have (2.9) with the modified 'key' function Z (2)  c j e 2 j (t ).
Remark 2.5.The above modification step is applicable to cohorts c = c m+1 , c m+2 , . . ., c m+w−S because no mortality data are available beyond the age of S years for cohorts beyond c m+w−S .

Data analysis
In this section, we leverage the modified mortality functions described in (2.8) and (2.9) to make predictions for selected future cohorts.In our comparative analysis, we evaluate not only these two modified approaches but also consider parametric models proposed in previous works by Shimizu et al. [11,12].

Swedish mortality
The Swedish mortality data available in the Human Mortality Database are among the oldest records, rendering them well-suited for long-term projections.We make use of empirical mortality functions q D AT A It is important to note that, as discussed in Remark 2.2, we already have mortality data available up to the year c m + w = 1940 (for the purpose of this analysis, we assume the current year to be 1940).Consequently, individuals born in cohort c are already (1940−c) years old, and our forecasting takes into account the age (1940 − c) + 1, assuming the individual has lived up to age S. For example, when predicting the cohort c = 1910, our analysis pertains to individuals over the age of 31, and we compare the curve q c (t |30) (for 31 ≤ t ≤ 110).
The parametric models for key functions used for comparison below are the following setup, similar to Shimizu et al. [12]: and T c is a change point where the mortality changes drastically, fixed at T c = 50.Note that it is too hard to estimate and predict T c in practice; see Shimizu et al. [11].
All parameter estimation and forecasting methods were performed according to the methodologies by Shimizu et al. [12].Subsequently, we maintain the values of x c representing the initial energy, κ c in (2.2), and σ c in (2.3) at a fixed set of values as follows: These specific values have been empirically selected based on insights obtained from data analysis, as detailed in Shimizu et al. [12,13].They are chosen to ensure that the estimation and optimization processes remain numerically stable across all cohorts.

Forecasting of M c and Λ c
We employed the R package fda to estimate the eigenfunctions e k j and eigenvalues λ k j (k = 1, 2).In determining the number of terms, K 1ϑ in equation (2.7), we set ϑ = 0.995, which yielded K 1ϑ = 3.The output includes the mean function µ(t ) = M (t ) and the eigenfunctions e 1 j for j = 1, 2, 3, as displayed in Figure 1 (a).The same procedure for Λ c obtains Figure 1 (b).
The predictions for the PCA scores Z (1)   c j for cohorts c > c m are conducted using forecast with ARIMA models.The results of these predictions are given in Figure 2, (a)-(c), and the forecasted trajectories of M δ c (t ) modified based on the prediction intervals are illustrated in (d).
To obtain a modified key function, M δ c (t ), we employ the δ = 95%-prediction interval (represented by the grey range in the Figure 2 (a)-(c)).The graphs reveal that the function for M c exhibits a tendency to increase year by year.Given that M c originally represents a drift of survival energy, this suggests that the trajectory of survival energy is on an upward trajectory year by year.Consequently, future cohorts are less likely to experience mortality.
Similarly, we employ analogous procedures to forecast Λ c , and the corresponding results are presented in Figure 4.

Parametric vs. Nonparametric
Using the key function M c derived from the preceding step, we compute a predicted mortality function and compare it with the results of the parametric model proposed by Shimizu et al. [12] (using the same model setup as in their paper).Our predictions are cohorts born in c = 1870, 1890, and 1910.In our data configuration, we used data up to the cohort born in 1830 (covering ages 20-110), placing us in the year 1940.Consequently, our predictions target individuals from the c = 1870, 1890, and 1910 cohorts at ages 70+, 50+, and 30+, respectively.Figures 3 and 5 display q I D c (t |30) and q IG c (t |30) for the male cohort born in c = 1910, which represents the longest projection.Additionally, Tables 1 and 2 The outcomes clearly indicate that our FDA approach outperforms the parametric models (except for a few rare cases) and offers superior accuracy in mortality prediction.As illustrated in Figures 3 and 5, our methodology demonstrates a propensity for heightened accuracy in mortality prediction, particularly when applied to the elderly demographic.Furthermore, this nonparametric approach eliminates the need to select change points T c required in ID-SEM.This constitutes a significant practical contribution, as accurate forecasting for the elderly demographic carries substantial overarching significance.Occasionally, ID-SEM outperforms the FDA approach in terms of predictions (e.g., 1890cohort for females).Nevertheless, the selection of the hyperparameter T c in ID-SEM bears significant importance.An erroneous choice in this regard can lead to a sudden deterioration in predictions, rendering the parametric approach notably unstable.Remark 3.1.As discussed by Shimizu et al. [12], the introduced "Modification" in Section 2.3 does not invariably lead to improved predictions.This is attributed to the risk of entering an "over-learning" state, induced by excessive LSE-fitting to existing data, which subsequently diminishes the model's predictive capability for future data.One potential remedy involves the reduction of the prediction threshold, e.g., δ = 0.8, etc.Nevertheless, based on our extensive experience in numerous data analyses, it is evident that, in most cases, this "Modification" indeed results in a reduction of MSE in the predictions.

Concluding remarks
We used a nonparametric approach with functional data analysis for cohort-specific mortality rate prediction based on Survival Energy Modeling (SEM with parametric models, which lack objective model selection criteria and make it difficult to identify data-driven models, our FDA approach allows for easy estimation and prediction of key functions.In fact, this approach was able to significantly outperform predictions made by existing parametric models in both the ID-SEM and IG-SEM.Additionally, the convenient environment of R package for the FDA makes these approaches readily applicable in practice without the need for extensive theoretical expertise, making it a valuable tool.While finding more powerful candidates than ID-SEM or IG-SEM for future predictions is a potential challenge, our data analysis suggests that ID-SEM already possesses sufficient predictive power at this stage.However, amending forecasts through the utilization of the prediction interval associated with the PCA score can occasionally yield counterproductive outcomes.This issue bears a resemblance to the concept of "overlearning", where excessive adaptation to historical data results in inadequate performance when applied to future data.One potential remedy involves the reduction of the prediction threshold denoted as δ.Nonetheless, based on our empirical findings, even when errors worsen, they tend to be of minor magnitude, and, in the majority of instances, such adjustments lead to an enhancement in forecasting accuracy. As discussed by Shimizu et al. [13], cohort-specific mortality rate prediction with SEM can be a powerful tool for calculating insurance premiums for each cohort, predicting changes in pension amounts, and financial forecasting.Particularly in countries with increasing life expectancy, such as Japan, the analysis of longevity risk is a crucial task, and SEM is likely to become a valuable tool for addressing this challenge.▲▲▲▲ ▲▲▲ ▲▲▲▲ ▲▲▲ ▲▲▲ ▲▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▲ ▲▲ ▲▲▲▲▲▲▲▲▲▲▲▲▲     ▲▲▲▲ ▲▲▲ ▲▲▲▲ ▲▲▲ ▲▲▲ ▲▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▲ ▲▲ ▲▲▲▲▲▲▲▲▲▲▲▲▲  1910 (t |30) predictions using parametric (Shimizu et al. [12]) and nonparametric (FDA) methods.In this case, the FDA method outperforms the parametric one overall.
provide a comparison of the Mean Squared Error (MSE) of males and females predicted mortality functions for these cohorts: for c = 1870, 1890, 1910 with c = 1940 − c, and for w = 110,

Figure 1 :
Figure 1: Figure (a) is the estimated mean functions M (t ) and eigenfunctions for M * c (t ; α c ). Figure (b) is the one for Λ c (t ) and Λ * c (t ; β c ).

Figure 2 :
Figure 2: Figures (a)-(c) illustrate the predictions for the 1st to 3rd PCA scores.In each figure, the blue line represents the mean, while the purple and grey ranges denote the 80% and 95% prediction intervals, respectively.Figure (d) provides the forecast for M δ c for cohorts 1840 ≤ c ≤ 1910 with δ = 0.95.The grey curves correspond to the estimated past curves M c for cohorts 1781 ≤ c ≤ 1830.

Figure 3 :
Figure 3: Comparison between the Swedish Male ID-SEM q I D1910 (t |30) predictions using parametric (Shimizu et al.[12]) and nonparametric (FDA) methods.It's evident that the FDA approach outperforms the parametric one, particularly in predicting the elderly population.

Figure 4 :
Figure 4: Figures (a)-(c) illustrate the predictions for the 1st to 3rd PCA scores.In each figure, the blue line represents the mean, while the purple and grey ranges denote the 80% and 95% prediction intervals, respectively.Figure (d) provides the forecast for Λ δ c for cohorts 1840 ≤ c ≤ 1910 with δ = 0.95.The grey curves correspond to the estimated past curves Λ c for cohorts 1781 ≤ c ≤ 1830.

Figure 5 :
Figure 5: Comparison between the Swedish Male IG-SEM q IG1910 (t |30) predictions using parametric (Shimizu et al.[12]) and nonparametric (FDA) methods.In this case, the FDA method outperforms the parametric one overall.

Table 1 :
).In contrast to existing research Parametric vs. FDA (Males): MSEs for predictions of Swedish mortality functions of cohort c = 1870, 1890 and 1910, in ID-SEM and IG-SEM.

Table 2 :
Parametric vs. FDA (Females): MSEs for predictions of Swedish mortality functions of cohort c = 1870, 1890 and 1910, in ID-SEM and IG-SEM.