Forecasting Age- and Sex-Specific Survival Functions: Application to Annuity Pricing

Abstract

We introduce the function principal component regression (FPCR) forecasting method to model and forecast age-specific survival functions observed over time. The age distribution of survival functions is an example of constrained data whose values lie within a unit interval. Because of the constraint, such data do not reside in a linear vector space. A natural way to deal with such a constraint is through an invertible logit transformation that maps constrained onto unconstrained data in a linear space. With a time series of unconstrained data, we apply a functional time-series forecasting method to produce point and interval forecasts. The forecasts are then converted back to the original scale via the inverse logit transformation. Using the age- and sex-specific survival functions for Australia, we investigate the point and interval forecast accuracies for various horizons. We conclude that the functional principal component regression (FPCR) provides better forecast accuracy than the Lee–Carter (LC) method. Therefore, we apply FPCR to calculate annuity pricing and compare it with the market annuity price.

1. Introduction

Accurate mortality forecasts are crucial for social, economic, financial, and medical decisions, but no consensus exists about what mortality variable to use and how to produce its forecasts (see Basellini et al. 2023; Booth and Tickle 2008, for detailed reviews). A standard life table contains age-specific death rates, probabilities of death, survival probabilities, life-table deaths, and life expectancy. The most commonly modeled mortality variable is the age-specific mortality rate, popularized by the Lee and Carter (1992) model. However, the LC model assumes a constant mortality improvement rate over time, which would under-estimate life expectancy if there is an accelerating trend in mortality improvements (Booth et al. 2002). This becomes a significant issue when health advancements and socioeconomic factors have accelerated mortality improvements in recent decades (Siu et al. 2011).

This deficiency prompts many researchers to consider other mortality variables, such as life-table deaths (see, e.g., Bergeron-Boucher et al. 2017, 2018; Oeppen 2008) or age-specific probabilities of death (see, e.g., Cairns et al. 2006, 2009). The life-table deaths resemble a probability density function with a non-negativity constraint and summability to the radix ${10}^5$. Similarly, the age-specific probabilities of death have two constraints: the non-negativity constraint and summability to one. When forecasting the probabilities of deaths, the logit transformation is often preferred due to the constraints (see, e.g., Li et al. 2015; Sweeting 2011). We study the age-specific survival function, similar to the cumulative probability of deaths, whose values are bounded in a unit interval. Compared to the life-table deaths or probabilities of death, this mortality variable has only one constraint, which is that its values lie within the unit interval.

With a logit transformation, a time series of survival functions is transformed to the unbounded data without constraints except for one in the last age group. For all age groups but the last one, we introduce a functional time-series forecasting method to model the historical patterns and extrapolate them into the future. The forecasting method performs a functional principal component analysis to extract the first few functional principal components and their associated principal component scores. The principal component scores can be modeled and forecasted by a univariate time-series forecasting method. By multiplying the forecasted scores with the estimated functional principal components, we obtain h-step-ahead forecasts within the transformed space. By taking the inverse logit transformation, the forecast survival functions are obtained by adding zeros for the last age group.

In the insurance industry, insurers rely on accurate mortality forecasts to price annuities or set reserves for mortality-related products. Annuity sales in voluntary markets have historically been low in Australia, primarily because of the availability of publicly funded age pensions and mandatory private retirement savings, such as superannuation. The demand for and profitability of annuity products is impacted by the misalignment that exists between retirees’ perceptions of their remaining lifespan and insurance companies’ assumptions and projections. In addition, heterogeneity in mortality risk across various factors further dampens demand for generic annuities (Brown and McDaid 2003).

However, the landscape has evolved in recent years. The desire for greater income, increased choice, and enhanced flexibility during retirement has surged, particularly fueled by rising inflation rates. Lifetime or fixed-term income stream products, particularly those offering investment-linked options, have witnessed growing popularity. With investment risk effectively transferred from insurers to policyholders, the level of confidence insurers place in their mortality forecasts significantly influences annuity pricing.

In this paper, we employ functional principal component analysis to forecast age-specific survival functions. We begin by introducing age- and sex-specific survival functions for Australia with logit transformation in Section 2. In Section 3, we apply Hyndman and Ullah (2007)’s method to model and forecast a time series of age- and sex-specific survival functions and then compare it with the widely used Lee–Carter method. In Section 4, we develop a single-premium temporary immediate annuity pricing model based on survival functions and subsequently compare the results with annuity prices on sales. Section 5 concludes with reflections on broader applications for the methods presented herein.

2. Age-Specific Survival Functions

We consider age- and sex-specific survival functions for Australia. These data can be publicly obtained from the Human Mortality Database (2024). For a given year t, the survival rates reduce as age increases. There are 111 ages, which are age $0,1,\dots ,109,110+$. At the last age group of 110+, the survival rates reach their minimum value of zero.

The survival function represents the probability that an individual survives to age x. According to this conditional feature, we construct survival function $S\left(x\right)$ as follows:

$$S(x+1)=S\left(x\right)\times (1-q_x)$$

where $q_x$ is the mortality rate at which people living to age x are expected to die before reaching age x + 1 and $S\left(0\right)=1$. In another way, we can obtain the survival function $S\left(x\right)$ from the following expression:

$$S\left(x\right)=\prod _{i=0}^{x-1}(1-q_i)$$

where i represents all time points from age 0 to $x-1$.

To understand the principal features of the data, Figure 1 presents rainbow plots of the age- and sex-specific survival functions in Australia from 1921 to 2020. The time ordering of the curves follows the color code of a rainbow, where data from the distant past years are shown in red, and data from more recent years are shown in purple (see, e.g., Hyndman and Shang 2010). For both males and females, the increasing ‘rectangularization’ (Wilmoth and Horiuchi 1999) of the survival curve indicates improving mortality across the lifespan. In contrast to the death rates, the survival functions are smooth.

The logit function can map data within bounded intervals onto unbounded data. The inverse logit function can ensure that the predicted values remain within the unit interval, which is important for probabilities. By taking the logit transformation, we obtain a time series of transformed data in Figure 2. When the survival functions reach zero in the last age group, the logit transformation produces undefined values. So, we removed the ones from the logit transformation and added them to the forecasts.

3. Functional Time-Series Forecasting Method

3.1. Functional Principal Component Analysis

The model proposed by Lee and Carter (1992) is a widely utilized statistical tool in actuarial science and demography for forecasting mortality rates and life expectancy. The two-factor (Age–Period) Lee–Carter (LC) model can be expressed as follows:

$$ln\left(m_{x,t}\right)=a_x+b_x{k}_t+{ϵ}_{x,t},$$

(1)

where $m_{x,t}$ represents the central death rate at age x in year t, $a_x$ is the average log mortality at age x over time, $b_x$ measures the relative change in the log death rate at each age, $k_t$ is the general level of the log death rate in year t, and ${\epsilon }_{x,t}$ is the residual error term at age x and year t.

The Lee–Carter model in Equation (1) is under-determined. Therefore, two constraints were introduced to ensure the Lee–Carter model yields consistent results before re-estimating $k_t$ based on observed total deaths in each year. The first constraint limits the sum of $b_x$ to 1, while the second constraint requires that the sum of $k_t$ equals zero.

As an extension to the LC method, Hyndman and Ullah (2007) (HU) introduced the HU method, which employs Functional Principal Component Analysis to decompose a set of curves into orthogonal functional principal components and their independent principal component scores. The HU method enhances the LC method through three key points:

(1)

Smoothing the log mortality rates is achieved using weighted penalized regression splines, with a specific emphasis on preserving monotonicity at older ages.

(2)

Utilizing multiple principal components.

(3)

A wider array of univariate time series models can be employed to forecast the principal component scores.

Let $S\left(x\right)=[S_1\left(x\right),S_2\left(x\right),\dots ,S_n\left(x\right)]$ be a time series of cumulative probabilities of deaths, and let $\mathcal{X}\left(u\right)=[{\mathcal{X}}_1\left(u\right),{\mathcal{X}}_2\left(u\right),\dots ,{\mathcal{X}}_n\left(u\right)]$ be a time series of logit-transformed data defined on a common function support $u\in \mathcal{I}\subset R$; $\mathcal{I}$ is a function support and R is the real line. $S_t\left(x\right)$ and ${\mathcal{X}}_t\left(u\right)$, where $t\in \{1,\dots ,n\}$ are elements of the Hilbert space $\mathcal{H}$ equipped with the inner product $\langle w,v\rangle ={\int }_{\mathcal{I}}w\left(u\right)v\left(u\right)du$, where u represents a continuum. Each function is a square-integrable function with a finite second moment. A non-negative definite covariance function is given by

$$\begin{array}{cc}\hfill c(u,v)& :=\mathrm{Cov}\left[\mathcal{X}\right(u),\mathcal{X}(v\left)\right]\hfill \\ \hfill & =\mathrm{E}\left\{\right[\mathcal{X}\left(u\right)-\mu \left(u\right)\left]\right[\mathcal{X}\left(v\right)-\mu \left(v\right)\left]\right\},\hfill \end{array}$$

(2)

for all $u,v\in \mathcal{I}$. The covariance function $c(u,v)$ in Equation (2) allows the covariance operator of $\mathcal{X}$, denoted by $\mathcal{K}$, to be defined as

$$\mathcal{K}\left(\varphi \right)\left(u\right)={\int }_{\mathcal{I}}c(u,v)\varphi \left(v\right)dv.$$

Via Mercer (1909)’s lemma, there exists an orthonormal sequence $\left({\varphi }_k\right)$ of continuous functions in the square-integrable space and a non-increasing sequence ${\lambda }_k$ of positive numbers, such that

$$c(u,v)=\sum _{k=1}^{\infty }{\lambda }_k{\varphi }_k\left(u\right){\varphi }_k\left(v\right).$$

Via Karhunen (1947)–Loève (1978) expansion, a stochastic process ${\mathcal{X}}_t\left(u\right)$ can be expressed as

$$\begin{array}{cc}\hfill {\mathcal{X}}_t\left(u\right)& =\mu \left(u\right)+\sum _{k=1}^{\infty }{\beta }_{t,k}{\varphi }_k\left(u\right)\hfill \\ \hfill & =\mu \left(u\right)+\sum _{k=1}^K{\beta }_{t,k}{\varphi }_k\left(u\right)+e\left(u\right),\hfill \end{array}$$

where $\mu \left(u\right)$ denotes the mean function approximated by $\overline{\mathcal{X}}\left(u\right)={\displaystyle \frac{1}{n}}{\sum }_{t=1}^n{\mathcal{X}}_t\left(u\right)$, the kth set of the principal component scores, ${\beta }_{t,k}=\langle {\mathcal{X}}_t\left(u\right)-\mu \left(u\right),{\varphi }_k\left(u\right)\rangle$, is given by the projection of decentered data ${\mathcal{X}}_t\left(u\right)-\mu \left(u\right)$ in the direction of the kth eigenfunction ${\varphi }_k$, $e\left(u\right)$ denotes the model errors with zero mean and finite variance, and $K<n$ denotes the number of retained functional principal components.

There exist various methods for determining the number of retained components:

(a)

Scree plots or the fraction of variance explained by the first few functional principal components (Chiou 2012);

(b)

Pseudo-versions of Akaike information criterion and Bayesian information criterion (Yao et al. 2005);

(c)

Predictive cross-validation leaving out one or more curves (Ramsay and Silverman 2006);

(d)

Bootstrap methods (Hall and Vial 2006);

(e)

Eigenvalue ratio criterion (EVR) (Lam and Yao 2012).

As pointed out by Hyndman et al. (2013), the forecast accuracy is robust to over-estimating K, but under-estimating K can be detrimental. We consider $K=6$ as in Hyndman et al. (2013). In contrast to the Lee and Carter (1992) method, the higher-order terms capture additional non-random patterns.

Table 1. Sensitivity for choice of K.

Gender/K	K = 1	K = 2	K = 6
Male	0.9385	0.9984	0.9999
Female	0.9629	0.9994	1.0000

gives the total variance explained by different choices of K for comparison.

Conditional on the estimated mean function and functional principal components, the forecast curves can be obtained via the forecast principal component scores. Differing from the random walk with drift used in Lee and Carter (1992), a broader range of univariate time-series forecasting methods can be applied. While Hyndman and Shang (2009) proposed univariate time-series forecasting methods, Aue et al. (2015) considered multivariate time-series forecasting methods. As pointed out by Peña and Sánchez (2007), the multivariate time-series method can capture contemporaneous and instantaneous correlations among the principal component scores. However, it cannot handle a large number of parameters in the model. The univariate time-series method can not model cross-correlation between two sets of the principal component scores. However, it can handle the non-stationarity exhibited in any set of principal component scores. Between the univariate and multivariate time-series forecasting methods, we apply an exponential smoothing state space model of Hyndman et al. (2008).

Considering the typical retirement age in Australia, we set age 65 as the starting point. For instance, when determining the annuity price for an individual aged 65 to 110, a prediction spanning 45 years is required. Based on the Australian female and male data from 1921 to 2020, we forecast age- and sex-specific survival functions from 2021 to 2065 as Figure 3.

3.2. Comparison of the Point Forecast Accuracy

Backtesting is performed by splitting the Australian female and male survival functions into two sets: the training dataset from 1921 to 2000 and the testing dataset from 2001 to 2020. The forecast results from each method are then compared with the testing dataset to evaluate the method’s finite-sample performance. As an illustration, Figure 4 illustrates the forecast results of the survival function in 2001, 2010, and 2020. For both males and females in each sample year, the figure demonstrates that the FPCR curve aligns more closely with the actual data compared to the LC curve.

To evaluate and compare the forecast accuracy, we consider mean absolute forecast error (MAFE) and mean absolute percentage error (MAPE). The MAFE quantifies the average absolute deviation between predicted values and observed values. The MAPE resembles the MAFE but is represented as a percentage of the actual value. With our focus on the post-retirement stage, we observe individuals ranging from 65 to 110 years old, encompassing a total of 46 ages in h-step-ahead survival function forecast. The MAFE and MAPE are then defined as follows:

$$\begin{array}{cc}\hfill \mathrm{MAFE}\left(h\right)& ={\displaystyle \frac{1}{46\times (21-h)}}\sum _{x=65}^{110}\sum _{\zeta =1}^{21-h}\left|S(x,\zeta )-\widehat{S}(x,\zeta )\right|\hfill \\ \hfill \mathrm{MAPE}\left(h\right)& ={\displaystyle \frac{100}{46\times (21-h)}}\sum _{x=65}^{110}\sum _{\zeta =1}^{21-h}{\displaystyle \frac{\left|S(x,\zeta )-\widehat{S}(x,\zeta )\right|}{S(x,\zeta )}}\hfill \end{array}$$

where $S(x,\zeta )$ denotes the actual survival function for the $x\mathrm{th}$ age and $\zeta \mathrm{th}$ forecasting year, while $\widehat{S}(x,\zeta )$ denotes the point forecasts for survival function.

Figure 5 displays the 1-to-20-step-ahead point forecast errors averaged over 20 forecast horizons. The errors become larger for both methods in longer forecast horizons. The FPCR produces lower MAFE and MAPE values than the ones from the LC for the male and female data.

The FPCR produces more accurate survival function forecasts than the ones obtained from the LC method. Another finding is that our model fits better to Australian female data than the male data. A summary of the averaged test results can be found in

Table 2. Point forecast accuracy by average across all forecast horizons.

	Female		Male
Method	MAFE	MAPE	MAFE	MAPE
LC	0.0049	2.6965	0.0073	3.7589
FPCR	0.0032	1.6431	0.0032	1.9565

.

3.3. Comparison of the Interval Forecast Accuracy

Prediction intervals are crucial tools for assessing the probabilistic uncertainty associated with point forecasts. This uncertainty can arise from systematic factors like parameter or model uncertainty, as well as random fluctuations such as error terms in specific models. As emphasized by Chatfield (1993, 2000), interval forecasts are as vital as point forecasts for several reasons:

(1)

Quantify the levels of future uncertainty;

(2)

Help decision makers to understand the range of possible outcomes before making more informed choices;

(3)

Assess the quality of different forecasting methods more thoroughly;

(4)

Take varying assumptions to explore diverse scenarios.

To measure forecast uncertainty, we construct the pointwise prediction intervals by using a nonparametric bootstrap method. The forecast.ftsm function from the ftsa package (Hyndman and Shang 2024) in version 6.4 provides a solution to apply this method. The nonparametric bootstrap method can be applied through a simple procedure as below:

(1)

Sample Resampling: The nonparametric model can be utilized to estimate the h-step-ahead forecasts for the kth principal component scores ${\widehat{\beta }}_{n+h|n,k}$ and model fit errors ${\widehat{e}}_{n+h}$. That is, we can generate bootstrap samples of principal component scores and model residual terms by randomly sampling with replacements from the estimated principal component scores and residual terms.

(2)

Statistic Calculation: After obtaining the bootstrapped principal component scores and errors in the model residuals, one can then obtain the bootstrapped survival function $\widehat{S}\left(x\right)$. We take the first six principal components as suggested by Hyndman and Booth (2008) and express

$${\widehat{S}}_{n+h|n}^{\left(b\right)}\left(x\right)=\widehat{\mu }\left(x\right)+\sum _{k=1}^6{\widehat{\beta }}_{n+h|n,k}^{\left(b\right)}{\varphi }_k\left(x\right)+{\widehat{e}}_{n+h}^{\left(b\right)}\left(x\right),b=1,2,\dots B,$$

where $\widehat{\mu }\left(x\right)$ is the estimated mean function, ${\varphi }_k\left(x\right)$ denotes the kth functional principal component, ${\widehat{\beta }}_{n+h|n,k}^{\left(b\right)}$ denotes the bootstrapped principal component scores, ${\widehat{e}}_{n+h}^{\left(b\right)}\left(x\right)$ represents the bootstrapped model residuals, and $B=1000$ denotes the number of replications.

(3)

Repeat Resampling: Repeat processes 1 and 2 above numerous times to construct a bootstrap distribution. The nonparametric bootstrap method closely resembles a Monte Carlo simulation experiment. We choose to run it B times to ensure that calculations based on the bootstrap distribution remain unaffected by significant Monte Carlo errors. Consequently, we obtain 1000 values of the bootstrapped survival functions $\{{\widehat{S}}_{n+h|n}^{\left(1\right)}\left(x\right),{\widehat{S}}_{n+h|n}^{\left(2\right)}\left(x\right),\dots ,{\widehat{S}}_{n+h|n}^{\left(1000\right)}\left(x\right)\}$, which we use to compute the 95% pointwise prediction intervals of the forecast results.

4. Single-Premium Temporary Immediate Annuity Pricing

Annuities guarantee the purchaser regular future payments, typically monthly and often extending for life. Under this overarching concept, various annuity types cater to distinct financial objectives. The primary categories include fixed or variable annuities, temporary or lifetime annuities, and immediate or deferred annuities.

A fixed annuity ensures a consistent payment throughout the agreement’s term, remaining constant without decrease or increase. In contrast, a variable annuity’s value fluctuates in line with the performance of the underlying mutual funds it is invested in, allowing for potential increase or decrease. Temporary annuities provide payments for a specified period, such as 5, 10, or 20 years, while lifetime annuities cover the remainder of the annuitant’s life. An immediate annuity initiates payouts immediately after the annuitant makes a lump-sum payment to the insurer. In contrast, a deferred annuity commences payments at a future date that is predetermined by the annuitant.

Annuitants commonly purchase annuities to complement their existing retirement income sources, such as pensions and superannuation. An annuity offers an assured income stream to provide a sense of security to annuitants. Even if their other financial resources diminish, they will still receive an additional steady income stream from the annuity.

Guaranteed lifetime annuities can serve as a suitable option for individuals seeking a steady income to complement their superannuation, pensions, or other investments. However, lifetime annuities might be perceived as expensive due to several reasons, such as longevity risk, guaranteed payouts, and administrative costs.

4.1. Annuity Pricing Based on Survival Functions

We study term annuity pricing by considering maturity-specific contracts that make the annuity almost lifetime for those individuals aged over 65 (see also Shang and Haberman 2020). The available maturity dates range from 5 to 45, increasing in 5-year increments, to encompass ages up to 110. These temporary annuities guarantee a fixed income level that exceeds the income provided by a lifetime annuity for a similar premium. They serve as an alternative to lifetime annuities and offer the buyer the flexibility to explore the option of purchasing a deferred annuity at a later date.

The cost of an annuity, which has a maturity term of T up to 45 years, is a stochastic variable since it depends on the zero-coupon bond price and future mortality rates. The price of the annuity for an individual aged x years with a benefit of one Australian dollar per year under a constant interest rate (i%) environment can be expressed as follows:

$$a_x^T=\sum _{\tau =1}^T{\hspace{0.17em}}_{\tau }{p}_x{e}^{-i\tau }$$

(3)

where ${\hspace{0.17em}}_{\tau }{p}_x$ represents the survival probability that can be derived from methods discussed in Section 3. A cohort approach is adopted in calculating the survival probabilities. The probability of survival for a person at time $t=0$ (or year 2021) for $\tau$ years, given their current age x, can be obtained from the formula below:

$${\hspace{0.17em}}_{\tau }{p}_x={\displaystyle \frac{S(x+\tau ,\tau )}{S(x,0)}}$$

(4)

where $S(x+\tau ,\tau )$ represents the survival function for age $x+\tau$ at time $\tau$ (e.g., $\tau$ years after year 2021) and $S(x,0)$ is the survival function for an x-year-old person at time 0.

By considering annuity purchase ages from 65 to 110, 45-step-ahead point and interval forecasts of survival functions are constructed to measure forecast uncertainty. The maturity term of T years is selected carefully to meet the condition that sum of age and maturity is less than 110. By combining Equations (3) and (4), the annuity price for each purchase age and maturity date can be expressed as follows:

$$a_x^T=\sum _{\tau =1}^T{\displaystyle \frac{S(x+\tau ,\tau )}{S(x,0)}}{e}^{-i\tau }.$$

(5)

If we add the overhead cost loading L into the calculation, Equation (5) leads to

$$a_x^T=(1+L)\times \sum _{\tau =1}^T{\displaystyle \frac{S(x+\tau ,\tau )}{S(x,0)}}{e}^{-i\tau }.$$

Taking survival function forecast results from Section 3, we can calculate the best estimate of the annuity prices with different ages and maturity dates for Australian female and male policyholders in

Table 3. Australian female and male annuity prices with different ages and maturities (T).

		Maturity (T)
Sex	Age	5	10	15	20	25	30	35	40	45
Female	65	4.57	8.46	11.71	14.30	16.20	17.37	17.88	18.00	18.01
	70	4.52	8.28	11.24	13.37	14.64	15.17	15.28	15.29
	75	4.44	7.90	10.35	11.75	12.30	12.42	12.42
	80	4.24	7.17	8.78	9.39	9.50	9.51
	85	3.87	5.93	6.65	6.78	6.79
	90	3.22	4.29	4.47	4.48
	95	2.40	2.79	2.82
	100	1.63	1.74
	105	1.12
Male	65	4.53	8.33	11.47	13.99	15.91	17.20	17.81	17.94	17.95
	70	4.46	8.10	10.95	13.05	14.38	14.96	15.08	15.09
	75	4.35	7.68	10.03	11.43	12.00	12.11	12.11
	80	4.14	6.92	8.48	9.06	9.16	9.17
	85	3.75	5.68	6.35	6.47	6.48
	90	3.10	4.09	4.25	4.26
	95	2.33	2.68	2.71
	100	1.61	1.72
	105	1.14

.

We exclude overhead cost loading in our calculations, i.e., $L=0$. The reason is that overhead cost loading varies between companies and can be easily incorporated by multiplying the final results by $1+L$. The constant interest rate is set at 2.5%, based on Australia’s long-term inflation target of 2% to 3%.

The 95% pointwise prediction intervals for the annuity pricing results can also be computed for the Australian female and male residents in

Table 4. The 95% pointwise prediction intervals of age-specific survival functions, from which we calculate annuity prices with different ages and maturities (T) for Australian females and males, based on the assumption of a constant interest rate of 2.5%.

		Maturity (T)
Sex	Age	5	10	15	20	25	30
Female	65	(4.42, 4.63)	(7.89, 8.7)	(10.41, 12.26)	(12.01, 15.35)	(12.79, 18.00)	(13.06, 20.20)
	70	(4.31, 4.62)	(7.4, 8.67)	(9.38, 12.2)	(10.35, 15.22)	(10.69, 17.71)	(10.77, 19.68)
	75	(4.10, 4.62)	(6.65, 8.63)	(7.93, 12.07)	(8.38, 14.92)	(8.49, 17.16)	(8.50, 18.86)
	80	(3.78, 4.59)	(5.58, 8.52)	(6.24, 11.77)	(6.40, 14.35)	(6.42, 16.28)	(6.42, 17.59)
	85	(3.29, 4.54)	(4.40, 8.30)	(4.67, 11.27)	(4.70, 13.57)	(4.71, 15.11)
	90	(2.78, 4.42)	(3.40, 7.92)	(3.49, 10.56)	(3.49, 12.48)
	95	(2.08, 4.31)	(2.30, 7.63)	(2.31, 9.87)
	100	(1.54, 4.24)	(1.61, 7.19)
	105	(1.10, 4.05)
Male	65	(4.61, 4.63)	(8.37, 8.7)	(11.42, 12.27)	(13.58, 15.4)	(14.94, 18.15)	(15.69, 20.54)
	70	(4.50, 4.63)	(8.00, 8.69)	(10.64, 12.25)	(12.30, 15.38)	(13.17, 18.09)	(13.52, 20.36)
	75	(4.39, 4.63)	(7.54, 8.68)	(9.69, 12.23)	(10.79, 15.32)	(11.23, 17.9)	(11.30, 19.75)
	80	(4.21, 4.62)	(6.90, 8.67)	(8.43, 12.18)	(8.99, 15.11)	(9.09, 17.23)	(9.09, 18.10)
	85	(3.97, 4.62)	(6.04, 8.63)	(6.85, 11.98)	(6.99, 14.32)	(7.00, 15.40)
	90	(3.50, 4.59)	(4.87, 8.42)	(5.07, 11.15)	(5.08, 12.33)
	95	(2.79, 4.47)	(3.25, 7.57)	(3.26, 8.92)
	100	(1.61, 3.94)	(1.69, 5.56)
	105	(0.91, 2.77)

.

The assumption of a constant interest rate greatly affects the calculation of annuity prices. To gauge this sensitivity, we perform a sensitivity test using interest rate assumptions of 1% and 5%. Detailed results for both the Australian male and female residents are provided in the Appendix A.

4.2. Compare Results with Annuity Price on Sales

Australia’s annuity landscape would greatly benefit from increased competition and diversity among providers. As of the end of 2023, Challenger Life stands as Australia’s largest and dominant annuity provider, furnishing steady income to thousands of clients. According to Challenger (2024), the Challenger Lifetime Annuity (Liquid Lifetime) offers three options, so customers can tailor it to meet their needs. These three options can be differentiated by payment delivery options, namely, immediate payments, deferred payments, and market-linked payments. An immediate payment annuity allows individuals to convert a lump sum of money into a stream of income payments right away. A deferred payment annuity enables individuals to make either a lump-sum payment or periodic premiums to an insurance company and then receive regular income payments at a future date in return. Market-linked payments, as the name suggests, are tied to market indexes, influencing the payments customers receive.

We take the prices for immediate payment annuity from Challenger (2024) to compare with the annuity prices obtained from Section 4.1 based on the survival function forecast. The provided rates serve as estimates and are intended for illustration purposes only. Challenger typically updates these rates on a weekly basis. Prices are segregated by gender and are based on an initial investment of $100,000 to obtain annuity payments at different starting ages, as detailed in

Table 5. Challenger female and male annual payments—Liquid Lifetime (immediate payments).

		AMOUNT OF INFLATION PROTECTION
Sex	Age	Withdrawal Period	Full	Partial	Nil	RBA Cash Linked
Female	65	22 years	$4924	$5965	$6807	$7063
	70	18 years	$5583	$6593	$7366	$7586
	75	14 years	$6520	$7500	$8203	$8385
	80	10 years	$7986	$8947	$9586	$9761
Male	65	19 years	$5222	$6256	$7074	$7315
	70	15 years	$5997	$7006	$7760	$7966
	75	12 years	$7026	$8001	$8684	$8862
	80	9 years	$8709	$9666	$10,284	$10,473

.

The Challenger lifetime annuity permits withdrawals up to the specified withdrawal period for each starting age. The withdrawal amount is determined by the remaining product value on the withdrawal date, and may be subject to deduction of administration fees and other charges. Following the end of the withdrawal period, the annuity product is considered to have no remaining withdrawal value, and, therefore, cannot be withdrawn or have any purpose for withdrawal. Payment rates vary depending on the chosen inflation protection options, including full protection, partial protection, or no protection. The payment amount can also be linked to the Reserve Bank of Australia (RBA) cash rate rather than the consumer price index.

The prices linked to the longest maturity dates in our survival function model can be considered as lifetime benefits, as Human Mortality Database (2024) assumes no individuals are alive beyond age 110. To build up a fair comparison, we backtrack to determine the prices for the Challenger product, calculating the initial investment needed to obtain a $1 annual benefit for a lifetime for each payment option. We opt to compare the non-inflation-protected annuity payment options Challenger offers with the outcome of our survival function models for each starting age and gender in

Table 6. The survival function annuity prices were compared to those in Challenger (2024). The highlighted values represent payment rates without inflation protection, providing a basis for comparison with the payment rates derived from survival function results.

		Inflation Protection				Survival Function
Sex	Age	Full	Partial	Nil	RBA	1%	2.5%	5%
Female	65	$20.31	$16.76	$14.69	$14.16	$22.00	$18.01	$13.39
	70	$17.91	$15.17	$13.58	$13.18	$18.12	$15.29	$11.84
	75	$15.34	$13.33	$12.19	$11.93	$14.29	$12.42	$10.03
	80	$12.52	$11.18	$10.43	$10.24	$10.63	$9.51	$8.00
Male	65	$19.15	$15.98	$14.14	$13.67	$22.03	$17.95	$13.27
	70	$16.68	$14.27	$12.89	$12.55	$17.92	$15.09	$11.65
	75	$14.23	$12.50	$11.52	$11.28	$13.94	$12.11	$9.77
	80	$11.48	$10.35	$9.72	$9.55	$10.25	$9.17	$7.72

.

The comparison yields intriguing insights. For individuals commencing at ages 65 and 70, Challenger (2024)’s non-inflation-protected annuity prices align neatly between the 2.5% and 5% interest sensitivity test outcomes for survival function, which is logical given the RBA cash rate of 4.35% as of 17 March 2024. However, this seems not the case when we look into the results for starting ages of 75 and 80. At age 75, Challenger’s prices go closer to our model’s 2.5% interest results, but they converge towards the 1% interest outcomes at age 80. This suggests a potential divergence in mortality curve conservatism or added loadings reflecting the greater uncertainty in outcomes and wider prediction intervals for older age groups in Challenger (2024)’s pricing model. These discoveries carry significant weight, suggesting that insurers might be able to fine-tune annuity pricing to foster competition and boost business activity, provided our mortality forecast models attain adequate credibility.

5. Conclusions

Survival probabilities are less prone to random fluctuations due to the cumulative nature of the indicator (from birth to age x), suggesting that their use could provide more robust forecasts. We begin by introducing two distinct mortality forecasting methods employed for projecting survival functions: the Lee–Carter method and the functional principal component regression method. For the FPCR, we choose Hyndman and Ullah (2007)’s model to forecast a time series of age- and sex-specific survival functions. The LC method and the FPCR are subsequently subjected to a backtesting process, utilizing both Australian male and female mortality data spanning from 1921 to 2000 from the Human Mortality Database (2024).

Our objective is to validate the performance of these methods by comparing their projections for the period 2001 to 2020 with the actual data from the same dataset. After the backtesting phase, we recommend the use of the FPCR for forecasting based on multiple evaluation metrics, namely MAFE and MAPE. Following the selection of the FPCR, this article proceeds to apply this method to forecast survival functions for the Australian population from 2021 to 2065. The focus of the forecast is on age groups ranging from 65 to 110, known as the retirement stage.

The survival function forecasts for the Australian population are utilized in a single-premium temporary immediate annuity pricing model. This model facilitates the calculation of annuity prices for Australian policyholders across various starting ages and maturity dates. The integration of these forecast outcomes into the annuity pricing model offers insights into the pricing dynamics for single-premium annuities for policyholders in Australia. The annuity pricing results also come with 95% pointwise prediction intervals, offering a measure of forecast uncertainty.

We compare the annuity prices derived from survival functions with the Challenger immediate payment annuity prices for each starting age. It appears that Challenger sets higher annuity prices for older age groups, which may be due to the inherent uncertainty in long-term mortality forecasts. To ensure reproducibility, the codes for generating point and interval forecasts using our method are available at https://www.github.com/wondersky8086/survival-function-forecast.git (accessed on 19 June 2024).

There are several ways in which the methodology presented here can be further extended, and we briefly mention two:

(1)

We consider an independent functional time-series forecasting method to separately model and forecast female and male survival functions. One could consider a joint modeling technique, such as multivariate and multilevel functional time-series methods of Shang and Kearney (2022), which captures the possible correlation between the series.

(2)

Our data consist of all available years from the Human Mortality Database (2024), but it may be subject to structural change. One possible idea is to apply a change point detection method from Shang and Xu (2022) to identify potential change points, and then fit our forecasting method to a truncated sample period.