1 Introduction
Over the past century, human life expectancy has trended up persistently. For instance, in Australia, (period) life expectancy rose from 60.97 years in 1921 to 83.47 years in 2021, and in New Zealand it increased from 68.81 years in 1948 to 82.30 years in 2021. These sustained gains, driven by continual improvements in public health, medical care, and living conditions, represent one of the most remarkable demographic shifts in modern history. However, longer lifespans also have substantial socioeconomic implications, including population ageing and increasing pressure on pension and healthcare systems. Consequently, the accurate modelling and forecasting of mortality and longevity are of utmost importance for policymakers, demographers, and actuaries in assessing the economic and social impact of longevity.
Since the influential contributions of Lee and Carter (1992) and Cairns et al. (2006), mortality modelling and forecasting have developed into a substantial and active field of research. The Lee-Carter (LC) family of models generally incorporates an age-specific component for capturing the age schedule and a time-varying component reflecting overall mortality improvements over time. By contrast, the Cairns-Blake-Dowd (CBD) family focuses on older ages and incorporates period effects linked to age in a more parsimonious structure. Their subsequent extensions integrate additional terms like cohort effects across birth cohorts. A wide range of these factor-based models has been proposed, with comprehensive reviews conducted by Cairns et al. (2009) and Haberman and Renshaw (2011). More recently, Li and Lu (2017) introduced the spatial-temporal vector autoregressive (STAR) model, imposing structured sparsity to achieve parsimony. Li et al. (2024) provided a review of its subsequent extensions, which further expand the methodological landscape of mortality forecasting.
In this paper, we introduce the demofit R package (Li (2026)), which provides a unified methodological and computational framework for mortality modelling by integrating parametric mortality laws and stochastic mortality models within a coherent interface. The package includes a collection of practical mortality curves with both constrained and unconstrained parameterisations, where the former offer greater interpretability and the latter provide increased flexibility in fitting complex mortality patterns. To enhance robustness, model estimation is performed using multiple optimisation routines, with the final selection based on minimising a least squares criterion. Moreover, the demofit package incorporates three major families of stochastic mortality models, namely the LC, CBD, and vector autoregressive (VAR) frameworks, and also ensemble and multi-population features, together with their estimation and forecasting procedures. A key feature of the package is that the forecasts produced from these stochastic models can further be guided by a selected mortality curve, which helps produce more realistic age patterns, particularly for small or noisy populations encountered in practical applications such as insurance pricing. The overall design emphasises a unified workflow and ease of use for empirical research, actuarial practice, and teaching purposes.
In the next section, we present the proposed framework and outline the list of parametric mortality laws and stochastic mortality models included. We then introduce the demofit package and compare it with a few related existing R packages. Afterwards, we apply the demofit package to specific datasets for illustration. Finally, we conclude with other remarks and a discussion of possible future developments.
2 Mortality modelling
This section provides an overview of the parametric mortality laws and stochastic mortality models included in the demofit package. For the parametric mortality curves, a unified methodological and computational framework is adopted for model fitting, enabling consistent estimation and comparison across different specifications. These mortality curves are chosen to cover a broad range of modelling needs, including different age ranges, levels of flexibility, and practical applications, thereby providing users with a versatile set of tools for mortality analysis.
Let \(m_x\) be the central death rate at age \(x\). Users can choose one of the following fourteen mortality curves, with their indicative age ranges (or subsets thereof) stated in the following brackets (see Gompertz (1825), Makeham (1860), Oppermann (1870), Thiele (1871), Wittstein and Bumsted (1883), Perks (1932), Weibull (1951)1, Beard (1971), Heligman and Pollard (1980), Rogers and Planck (1983), Siler (1983), Martinelle (1987), Thatcher (1999), and Tabeau (2001)2): \[ \begin{array}{l l l} \text{Gompertz 1825} & m_x = B e^{Cx} & (30\text{--}90) \\ \text{Makeham 1860} & m_x = A + B e^{Cx} & (20\text{--}90) \\ \text{Oppermann 1870} & m_x = \frac{A}{\sqrt{x+1}} + B + C\sqrt{x+1} & (0\text{--}20) \\ \text{Thiele 1871} & m_x = A_1 e^{-B_1 x} + A_2 e^{-0.5 B_2 (x-C)^2} + A_3 e^{B_3 x} & (0\text{--}90) \\ \text{Wittstein and Bumsted 1883} & m_x = \frac{A^{-(Bx)^N}}{B} + A^{-(M-x)^N} & (0\text{--}90) \\ \text{Perks 1932} & m_x = \frac{A+BC^x}{1+DC^x} & (20\text{--}100+) \\ \text{Weibull 1939} & m_x = Bx^C & (50\text{--}100+) \\ \text{Van der Maen 1943} & m_x = A + Bx + Cx^2 + \frac{I}{N-x} & (80\text{--}100+) \\ \text{Beard 1971} & m_x = \frac{A e^{Cx}}{1+B e^{Cx}} & (50\text{--}100+) \\ \text{Heligman and Pollard 1980} & m_x = A^{(x+B)^C} + D e^{-E(\ln(x)-\ln(F))^2} + \frac{GH^x}{1+GH^x} & (0\text{--}100+) \\ \text{Rogers and Planck 1983} & m_x = A_0 + A_1 e^{-Ax} + A_2 e^{-B(x-U)-e^{-C(x-U)}} + A_3 e^{Dx} & (0\text{--}100+) \\ \text{Siler 1983} & m_x = A_1 e^{-B_1 x} + A_2 + A_3 e^{B_3 x} & (0\text{--}90) \\ \text{Martinelle 1987} & m_x = \frac{A+B e^{Cx}}{1+D e^{Cx}} + E e^{Cx} & (20\text{--}100+) \\ \text{Thatcher 1999} & m_x = A + \frac{B e^{Cx}}{1+D e^{Cx}} & (20\text{--}100+) \\ \end{array} \]
For consistency within the unified framework, all the mortality curves are expressed in terms of the central death rate \(m_x\). Under the common assumption that the force of mortality is constant within each age interval (e.g. Li (2013)), the force of mortality \(\mu_x\) coincides with the central death rate \(m_x\), allowing those mortality curves originally formulated in terms of \(\mu_x\) to be represented equivalently. From a practical perspective, \(m_x\) can be readily estimated from observed annual data as the ratio of deaths to mid-year exposure (e.g. HMD (2026)), making it a natural choice for empirical implementation. In addition, empirically, for most ages, the one-year death probability \(q_x\) satisfies \(q_x \approx 1-e^{-\mu_x} = 1-e^{-m_x} \approx m_x\), which justifies expressing the mortality curves originally defined in terms of \(q_x\) using \(m_x\) instead. At very old ages like 100+, however, there is conflicting empirical evidence on how \(\mu_x\) behaves with age (increasing, decreasing, or flattening), and at such extremely old ages, the numerical values of \(q_x\) and \(m_x\) can become quite different. For instance, at around age 100, \(m_x\) often reaches 0.5 for many populations, for which the corresponding \(q_x\) is only around 0.4. Consequently, the exact formulas \(q_x = \int_0^1 {}_tp_x\,\mu_{x+t}\,dt\) and \(m_x=\frac{q_x}{\int_0^1 {}_tp_x\,dt}\) need to be used. Nevertheless, those mortality curves originally formulated in terms of \(q_x\) are generally much less suitable for capturing these centenarian mortality patterns anyway, and alternative formulations defined directly on \(\mu_x\) or \(m_x\), such as the logistic curve in Thatcher (1999), are more appropriate.
For the stochastic mortality models, a unified computational structure is similarly adopted to ensure consistency in model specification, estimation, and forecasting. These stochastic mortality models are selected to capture a wide range of dynamic mortality patterns, including different age structures, temporal behaviours, and cohort patterns, offering users a flexible and coherent set of tools for mortality forecasting. This coherent approach facilitates systematic comparison across different forecasting models and enables their integration with the parametric mortality curves within a consistent modelling framework.
Now let \(m_{x,t}\) be the central death rate at age \(x\) in year \(t\), extended with the time dimension. Users can select one of the seven forecasting models below: \[ \begin{array}{l l} \text{Lee-Carter (LC)} & \ln(m_{x,t}) = \alpha_x + \beta_x \kappa_t + \epsilon_{x,t} \\ \text{Renshaw-Haberman (RH)} & \ln(m_{x,t}) = \alpha_x + \beta_x \kappa_t + \gamma_{t-x} + \epsilon_{x,t} \\ \text{Age-Period-Cohort (APC)} & \ln(m_{x,t}) = \alpha_x + \kappa_t + \gamma_{t-x} + \epsilon_{x,t} \\ \text{CBD (M5)} & \ln(m_{x,t}) = \kappa_{1,t} + \kappa_{2,t} (x-\bar{x}) + \epsilon_{x,t} \\ \text{CBD with Cohort (M6)} & \ln(m_{x,t}) = \kappa_{1,t} + \kappa_{2,t} (x-\bar{x}) + \gamma_{t-x} + \epsilon_{x,t} \\ \text{CBD with Curvature and Cohort (M7)} & \ln(m_{x,t}) = \kappa_{1,t} + \kappa_{2,t} (x-\bar{x}) + \kappa_{3,t} ((x-\bar{x})^2-\sigma^2) + \gamma_{t-x} + \epsilon_{x,t} \\ \text{Spatial-Temporal VAR (STAR)} & \ln(M_{t}) = \mu + R \ln(M_{t-1}) + E_t \\ \end{array} \]
For consistency and similar reasons as above, all the forecasting models are also expressed in terms of the central death rate \(m_{x,t}\), where \(M_t\) is the vector of central death rates in year \(t\). For the first three models, \(\alpha_x\) refers to the age effect, \(\kappa_t\) allows for the period effect with sensitivity \(\beta_x\), and \(\gamma_{t-x}\) represents the cohort effect. The constraints imposed include \(\sum_x \beta_x = 1\), \(\sum_t \kappa_t = 0\), and \(\sum_c \gamma_c = 0\). For the CBD models, \(\kappa_{1,t}\) is the intercept parameter, \(\kappa_{2,t}\) is the slope parameter, \(\kappa_{3,t}\) is the curvature parameter, \(\bar{x}\) is the average age of the age range, \(\sigma^2\) is the average value of \((x-\bar{x})^2\), and \(\gamma_{t-x}\) is the cohort effect with \(\sum_c c^i \gamma_c = 0\) for \(i\) = 0, 1, 2. For the STAR model, \(R\) is a banded matrix with non-zero entries in the main diagonal and the two subdiagonals below that. The error terms are \(\epsilon_{x,t}\) and \(E_t\). Note that the LC, RH, APC, and STAR models can be applied to the whole age range or subsets of it, while the three CBD models are designed mainly for ages 50-90. More details of these forecasting models can be found in Lee and Carter (1992), Bray (2002), Cairns et al. (2006), Cairns et al. (2009), Haberman and Renshaw (2011), Li and Lu (2017), and Li et al. (2024).
Within our integrated framework, the forecasts generated from the stochastic mortality models can be further regularised flexibly through a selected parametric mortality curve. Although the stochastic models effectively capture temporal dynamics, the resulting forecasted age profiles may exhibit irregular, non-smooth, or demographically unexplainable patterns across ages, particularly for small or noisy populations, in the presence of cohort parameters, or when the forecasts are anchored to the observed values in the final year of the data. Such features can pose challenges for practical implementation, such as insurance pricing and communication with stakeholders. By further fitting a parametric curve to the forecasted rates for each projection year, the age patterns become smoother and more structured, while remaining consistent with the underlying temporal evolution. This pragmatic approach facilitates interpretation and implementation by expressing the forecasts in a functional form familiar to insurance practitioners and demographers, while maintaining coherence with the stochastic modelling framework. Note that the CBD (M5) model does not have this issue though, as it is by nature a simple linear curve fitted to each year. In this aspect, McNown and Rogers (1989) proposed a more natural approach of directly projecting the parameters of a mortality curve. As commented in Li and Liu (2019), however, estimating and forecasting too many parameters of a mortality curve (e.g. eight parameters in the Heligman-Pollard curve) can lead to a biologically unreasonable age schedule, despite the mathematical smoothness of the forecasts.
3 Package demofit
The demofit package provides a unified interface for fitting mortality curves and forecasting mortality rates within a coherent workflow. For mortality curve fitting, the core function MC() (MC denoting “mortality curve”) fits a selected mortality curve to observed mortality rates by picking the solution with the smallest weighted sum of squares value on the log death rates. The package includes both constrained and unconstrained parameterisations of the mortality curves. The former is in line with the traditional designs and provides greater interpretability, while the latter offers increased flexibility for capturing more complex age patterns. To enhance robustness and numerical stability, multiple optimisation strategies are employed simultaneously, including the PORT routines, Nelder-Mead method, and Levenberg-Marquardt algorithm, with the final estimates chosen based on the best fit among the three methods. The fitted curve is returned as an S3 object with associated methods for obtaining parameter estimates (coef), fitted values (fitted), interpolated / extrapolated values (predict), weighted sum of squares (deviance), residuals (residuals), and graphical display (plot).
This multi-strategy approach differs from the existing MortalityLaws package (Pascariu (2025)), where optimisation is carried out using a single selected routine or loss function. By contrast, demofit applies several optimisation methods in parallel and selects the best-performing solution according to a common fitting criterion, thereby improving robustness to starting values, mitigating convergence issues, and producing more stable results for data of different calendar years and populations over a wide range of mortality curves. (See the next section for examples.) In addition, the function CC() is provided for evaluating and constructing a mortality curve given a set of parameters, allowing users to directly generate the fitted mortality schedule.
For stochastic mortality modelling, the demofit package provides the function FCS() (FC standing for “forecasting” and S for “smoothed”) to fit a range of forecasting models within the same framework. All the models are formulated in terms of the log death rate, and estimation is generally carried out by minimising the sum of squared residuals (i.e. Gaussian likelihood). For computational efficiency, the LC model is estimated via singular value decomposition (SVD)3 (e.g. Renshaw and Haberman (2006)), and the CBD (M5) and STAR models are fitted using regression-based methods (via OLS, SQP). For the other more complex specifications, parameter estimation is performed using an iterative Newton updating scheme (e.g. Renshaw and Haberman (2006), Li (2013)). After model fitting, the period and cohort parameters (\(\kappa_t\), \(\kappa_{1,t}\), \(\kappa_{2,t}\), \(\kappa_{3,t}\), \(\gamma_c\)) are forecasted using the function auto.arima() in the forecast package (Hyndman et al. (2026)). Moreover, there is a jumpoff argument where users can make the forecasts from the first six models anchored to the observed values in the final year of the data and so by nature more in line with the STAR forecasts. (As noted in Lee and Miller (2001), an adjustment to the correct jump-off level can give a more accurate life expectancy forecast and remove the bias from simply using the fitted values.) These settings ensure consistency in model fitting and facilitate comparison between different forecasts. The fitted model is returned as an S3 object with associated methods for extracting parameter estimates (coef), forecasts (forecast::forecast), standardised residuals (residuals), and graphical plots (plot).
In comparison, the StMoMo package (Villegas et al. (2025)) adopts a likelihood-based framework in which its mortality models are fitted assuming Poisson or Binomial distributions for deaths, with estimation carried out via generalised nonlinear models implemented in the gnm package (Turner et al. (2025)). The demofit package instead adopts a least squares formulation4 directly on the log death rates, as a default, providing a unified, numerically stable, and computationally efficient estimation approach across different model specifications. This facilitates systematic comparison and the integration of mortality curves with stochastic mortality models. Moreover, two features further extend the scope of the demofit package. First, the inclusion of the STAR model, the ensemble and multi-population features, and the option to regularise the forecasts as a post-processing step using different mortality curves expand the range of modelling approaches available within a coherent computational framework. (See the next section for illustrations.) Second, the demofit package provides tailored functions, including FCS() (FC for “forecasting” and S for “smoothed”), LCS() (LC model), RHS() (RH model), APCS() (APC model), CBDS() (CBD model), CBDCS() (CBD with cohort model), CBDQCS() (CBD with curvature and cohort model), STARS() (STAR model), ENS() (ensemble), CFMS() (common factor model), and CAES() (common age effect model), which are designed for ease of use in practical applications without the need for extensive configuration.
Overall, the demofit package provides a useful addition to the relatively small list of packages for mortality modelling in R. While the two packages MortalityLaws and StMoMo offer well-established implementations for specific modelling approaches, the range of available R tools remains rather limited compared to other areas of statistical modelling. By integrating parametric mortality laws and stochastic mortality models within a coherent and user-friendly framework, demofit complements the existing packages and broadens the set of practical tools available for actuarial and demographic applications.
Besides these frequentist approaches, recent Bayesian packages such as StanMoMo (Barigou and Goffard (2023)) and BayesMortalityPlus (Figueiredo et al. (2025)) provide complementary frameworks for mortality modelling. The StanMoMo package implements Bayesian versions of several stochastic mortality models, including the LC, RH, APC, CBD, and M6 models, together with Bayesian model averaging and uncertainty quantification. The BayesMortalityPlus package focuses on Bayesian mortality graduation (smoothing) using the Heligman-Pollard and dynamic linear models, together with Bayesian LC (only) mortality forecasting. It also provides mortality closure methods for advanced ages. In contrast, the demofit package is designed as a unified methodological and computational framework spanning both parametric mortality curves and stochastic mortality models under a common interface, with an emphasis on practical implementation, systematic model comparison, and computational efficiency instead of Bayesian inference. It additionally provides multi-population modelling and forecast regularisation using mortality curves, thereby complementing the existing Bayesian packages with a wider methodological scope for practical applications.
4 Examples
In this section, we first illustrate the use of the function MC() and compare it numerically with the corresponding fitting function in MortalityLaws. The fitted curve under MC() is stored as an S3 object with associated methods including coef() for parameter estimates, fitted() for fitted values, predict() for interpolation / extrapolation, deviance() for weighted sum of squares, residuals() for residuals, and plot() for graphical display. We obtain Australian and New Zealand mortality data from the Human Mortality Database (HMD (2026)) for the period 1950-2021 and ages 0-100. We follow Lee and Miller (2001) and base the analysis on data since 1950 so as to avoid the major structural changes in the distant past. As an example, we use the following code to load the packages, fit the Rogers-Planck curve (with the usual parameter constraints) to the death rates of a particular calendar year, and compute the fitting errors. The same procedure applies to the other mortality curves in both packages.
library(demofit)
library(MortalityLaws)
x <- 0:100
fit1 <- MC(x=x, m=m, curve="rogersplanck2")
err1 <- deviance(fit1)
fit2 <- MortalityLaw(x=x, mx=m, law="rogersplanck")
err2 <- sum((log(m) - log(fitted(fit2)))^2)Table 1 presents the average reductions (in %) in the errors of fitting thirteen mortality curves to the four populations across the data period. We do not include the Thatcher logistic curve (for centenarian mortality) as it is available in demofit only. Note that the values shown in parentheses correspond to the reductions obtained using unconstrained parameterisations instead. There are a number of interesting observations. First, for the mortality curves with only two or three parameters, both packages basically give rise to the same fitting errors. But for the other mortality curves with more than three parameters, demofit generates noticeable reductions in the fitting errors. These improvements are primarily driven by the multi-optimisation strategy, which reduces sensitivity to starting values and helps avoid local minima for more complex mortality curves, hence providing more robust estimation. (The only exception is the Heligman-Pollard curve, where MortalityLaws may have adopted special starting values or a tailored solution.) Second, using unconstrained parameterisations generally provides further reductions in the fitting errors, reflecting their greater flexibility in capturing detailed age patterns. Figure 1 illustrates the fitted Rogers-Planck curves for Australian female death rates in 2020 from different methods. It can be seen that demofit’s unconstrained parameterisation provides the best fit to the whole age range. The next is its constrained parameterisation version, which appears to miss some of the age patterns from around age 70 onward. It is then followed by the constrained version from MortalityLaws, where the age patterns both below age 20 and above age 70 are not adequately captured by comparison. (Their fitting errors are 7.97, 9.71, and 14.42, respectively.)
| Curve | Australia F | Australia M | New Zealand F | New Zealand M |
|---|---|---|---|---|
| Gompertz | 0.00 (0.00) | 0.00 (0.00) | 0.00 (0.00) | 0.00 (0.00) |
| Makeham | -0.00 (-0.00) | -0.00 (-0.00) | -0.00 (-0.00) | -0.00 (-0.00) |
| Oppermann | 0.00 (-0.00) | -0.00 (-0.00) | -0.00 (-0.00) | -0.00 (-0.00) |
| Thiele | 21.15 (24.17) | 8.91 (8.99) | 18.63 (19.66) | 12.91 (12.91) |
| Wittstein-Bumsted | 2.84 (2.84) | 0.42 (0.42) | 0.61 (0.61) | 0.24 (0.24) |
| Perks | -22.33 (8.35) | 0.01 (4.71) | -2.81 (3.07) | 0.20 (0.92) |
| Weibull | 0.00 (0.00) | 0.00 (0.00) | 0.00 (0.00) | 0.00 (0.00) |
| Van der Maen | 1.48 (11.57) | -21.45 (15.30) | 14.77 (17.27) | 17.98 (21.63) |
| Beard | -0.23 (45.69) | 0.00 (6.58) | 0.00 (9.86) | -0.00 (4.16) |
| Heligman-Pollard | 0.00 (0.00) | 0.00 (0.00) | 0.21 (0.23) | 0.00 (0.00) |
| Rogers-Planck | 24.46 (36.54) | 62.75 (64.16) | 21.19 (22.96) | 50.12 (49.95) |
| Siler | 0.64 (0.64) | 0.52 (0.58) | 0.29 (0.29) | -0.00 (0.09) |
| Martinelle | 2.42 (4.84) | 32.30 (54.44) | 2.59 (7.36) | 23.55 (35.84) |
Figure 1: Fitted mortality curves (Rogers-Planck) from different methods for Australian female death rates in 2020. (Crosses are actual observations.)
Furthermore, we demonstrate how the Thatcher curve (available in demofit only) is more suitable for centenarian mortality than the other traditional mortality curves. We collect high-quality old-age data of Sweden from the Kannisto-Thatcher Database (KTD (2018)). Sweden lies in the best data quality group in that database as it has a very long history of national population registers. Since the data at extreme ages are sparse, we aggregate the deaths and exposures over 1991-2010 and calculate the corresponding death rates at ages 80-105. Figure 2 displays the fitted Thatcher, Weibull, and Gompertz curves (constrained) for Swedish female death rates. It is visually clear that the Thatcher curve is able to capture the flattening of the log death rates at extreme ages, while the traditional Weibull and Gompertz curves, having power and exponential growth, fail to do so. (Their fitting errors are 0.00, 0.17, and 0.30, respectively.) In practical applications, when old-age data are scarce, one may fit the Thatcher curve to the age range where data are relatively more reliable and then use the fitted curve for mortality closure by extrapolating the death rates to advanced ages. Figure 2 also illustrates the extrapolated death rates up to age 120 (via predict()). The flattening of the Thatcher curve here is consistent with the broader evidence reported by Thatcher (1999), where the conclusions were drawn from pooled mortality data across 13 industrialised countries, analysed both collectively and for different periods and birth cohorts, thereby reducing the uncertainty associated with sparse observations at the oldest ages.
Figure 2: Fitted and extrapolated mortality curves (Thatcher, Weibull, Gompertz) for Swedish female death rates in 1991-2010. (Dots are actual observations.)
We then demonstrate the use of the functions LCS(), RHS(), APCS(), CBDS(), CBDCS(), CBDQCS(), and STARS()5 and examine their estimated parameters and forecasts. The fitted models are returned as S3 objects with associated methods including coef() for parameter estimates, forecast::forecast() for smoothed forecasts, residuals() for standardised residuals, and plot() for graphical plots. We continue to use the HMD data for the more recent period 1980-2021 and ages 0-90. As an illustration, we apply the code below to fit the LC, RH, APC, and STAR models to the death rates of a particular population, where the forecasts are regularised by the Heligman-Pollard curve (constrained), which is suitable for the whole age range. (Note that the raw forecasts, without adjustment, can still be obtained via, for instance, f1$forecast.) The same process applies to the CBD (M5, M6, M7) models, but the age range has to be restricted to ages 50-90 only. Figure 3 shows the forecasted age profiles in 2051 under the four models for New Zealand female death rates. First, the models differ in the projected shape of mortality at young adult ages, with the “accident hump” around age 20 being more pronounced under the LC and APC models. Second, noticeable differences emerge in the projected old-age mortality improvements, where the four models imply different rates of decline and hence different future age profiles. Third, while all four models preserve the overall age pattern of mortality, the flexibility introduced by the cohort parameters in the RH and APC models and the age-varying improvement structure in the STAR model leads to visibly different forecast trajectories compared with the classical LC model. Figure 4 compares the raw and smoothed forecasted death rates in 2051 under the APC model. For a small population like New Zealand, there is a relatively high level of noise in the data, which can be propagated into the forecasts, especially when the forecasts are anchored to the last observed year. The irregular and biologically unreasonable age patterns appear to be substantially alleviated by the selected Heligman-Pollard curve, which is familiar to actuaries and practitioners and is helpful for insurance product pricing and communication purposes. Note that the number of policyholders in a typical insurance portfolio is generally much smaller than the population size of a country and so the issue of forecasting irregular mortality schedules can be more severe in practical valuation work. More specifically, mortality schedules are often used directly in insurance pricing and valuation, where irregular age-to-age fluctuations can translate into unusual patterns in insurance premiums or other quantities derived from mortality rates. Such patterns may be difficult to justify from a demographic perspective and can also complicate the communication of insurance products to prospective customers. The regularised forecasts can therefore provide a more sensible and familiar mortality schedule for subsequent actuarial applications (Atkinson (2018)), without requiring the underlying stochastic mortality model to be replaced.
Figure 3: Mortality forecasts in 2051 (top left panel – LC; top right – RH; bottom left – APC; bottom right – STAR) using New Zealand female data of years 1980-2021 and ages 0-90.
Figure 4: Smoothed and unadjusted forecasted death rates in 2051 for New Zealand females under APC model.
We now reduce the age range to 50-90 so as to fit all the seven forecasting models to Australian male data using both demofit and StMoMo (with the binomial assumption therein). Figure 5 plots the parameter estimates with respect to their corresponding dimensions under different models. The parameter estimates under the LC and APC models (first and third panels) are generally close between the two packages. The differences under the CBD (M5, M6, M7) models (last three panels) are relatively more apparent, though the overall age / temporal / cohort patterns are still broadly in line between the two packages. These differences can be attributed mainly to the distributional assumptions and estimation algorithms of the two packages. The major differences lie in the RH model outputs (second panel), where the classical linearity in overall mortality decline (i.e. \(\kappa_t\)) is preserved by demofit even after the inclusion of the cohort effect, while this overall temporal effect appears to be “absorbed” into the cohort parameters produced by StMoMo. It is well known that separating age, period, and cohort effects is inherently challenging due to identification and computational issues, and different modelling and estimation choices may lead to alternative but statistically equivalent representations. In this context, demofit maintains the cohort effect as a secondary component while preserving the classical Lee-Carter-style temporal structure6, which has important implications for the stability, interpretability, and extrapolation behaviour of mortality forecasts. From a forecasting perspective, this decomposition avoids conflating long-term mortality improvement with cohort-specific effects, allowing the projected period trend to remain more stable and less sensitive to cohort fluctuations.
In addition, compared to the above six forecasting models, the STAR model (available in demofit only) has a distinctive feature in ensuring the age coherence property, whereby the ratio of forecasted death rates between any two ages converges to a constant over time. Figure 6 compares the ratios of forecasted death rates from 2022 to 2051 between the STAR and LC models. Under the STAR model, the forecasted ratios gradually converge, implying that the projected death rates across ages remain proportionally stable in the long term. By contrast, under the LC model, the forecasted ratios continue to increase or decrease over time, leading to persistent and growing divergence between ages. Such behaviour would be biologically unreasonable and inconsistent with the general expectation of stable relative mortality patterns across ages in the long term. Note that this long-run age coherence property is implied by the fitted STAR model specification, while the short- and medium-term forecast trajectories remain data-driven through the estimated model parameters and the latest trends in the data.
Moreover, we present the function CFMS(), which fits the common factor model proposed by Li and Lee (2005). The model specification is \(\ln(m_{x,t,i}) = \alpha_{x,i} + B_x K_t + \beta_{x,i} \kappa_{t,i} + \epsilon_{x,t,i}\), where \(B_x K_t\) is the common factor for two populations and \(\beta_{x,i} \kappa_{t,i}\) is the specific factor for population \(i\). The existing R packages for mortality modelling are primarily designed for single-population analysis and generally do not provide built-in support for multi-population modelling. We apply CFMS() to New Zealand female and male data (treated as population 1 and population 2, respectively) of years 1980-2021 and ages 50-90. As noted in Li (2013), under suitable time series assumptions, the ratio of forecasted death rates at each age between two related populations converges to a constant over time. Figure 7 shows how the ratios of forecasted death rates for both sexes at different ages converge under the common factor model. If the two populations are modelled and forecasted separately instead, this ratio is likely to increase or decrease indefinitely. The common factor model setting here ensures that there is a biologically sensible long-term relationship between female and male mortality. Similarly, this long-term convergence is a structural property of the common factor model under the assumed time series dynamics, whereas the estimated population-specific factors and the latest trends in the data determine the short-term forecasting behaviour.
f9 <- CFMS(x=x, M1=M1, M2=M2, curve="makeham2", h=30, jumpoff=2)Finally, we briefly introduce the function ENS(), which performs weighted averaging over the forecasts from the seven models. Some recent mortality studies (e.g. Li (2023)) suggest that this kind of ensemble or stacking approach can lead to enhanced forecasting accuracy, particularly in applications where forecasting accuracy is prioritised over model interpretability (i.e. age schedule, overall trend, cohort pattern, etc). As a backtesting exercise for illustration, we apply this function to New Zealand female and male data of years 1980-2011 and obtain both the unadjusted and smoothed forecasts for years 2012-2021. Table 2 reports the mean absolute errors of the forecasts under different forecasting models and the ensemble approach (with equal weights). There are a number of interesting observations. First, for both populations, the ensemble approach yields the lowest errors when the forecasts are unadjusted, followed by the STAR model. But when the forecasts are smoothed, the improvement from the ensemble approach is much reduced. Second, the smoothed forecasts produce lower errors than the unadjusted forecasts for both populations, suggesting that the smoothing step does not necessarily jeopardise forecasting accuracy. Table 3 presents the mean absolute errors of the 10-year forecasts for a different data period of years 1975-2006. While the ensemble approach appears to perform even better in this case, the smoothed forecasts continue to have lower errors than the unadjusted forecasts for each model tested. It is also worth noting that allowing the ensemble weights to be determined based on in-sample validation may further enhance the forecasting performance of the ensemble approach.
f8 <- ENS(x=x, M=M, curve="makeham2", jumpoff=2)
Figure 5: Parameter estimates (top to bottom panel – LC, RH, APC, M5, M6, M7) from demofit (solid lines) and StMoMo (dashed lines) using Australian male data of years 1980-2021 and ages 50-90.
Figure 6: Ratios of death rates (top left – 50 vs 60; top right – 60 vs 70; bottom left – 70 vs 80; bottom right – 80 vs 90) from STAR (dashed lines) and LC (dotted lines) models using Australian male data of years 1980-2021 and ages 50-90.
Figure 7: Female-to-male ratios of death rates at ages 50, 60, 70, and 80 from common factor model using New Zealand data of years 1980-2021 and ages 50-90.
| Model | NZ F (unadjusted) | NZ M (unadjusted) | NZ F (smoothed) | NZ M (smoothed) |
|---|---|---|---|---|
| LC | 0.1050 | 0.1070 | 0.0901 | 0.0921 |
| RH | 0.1183 | 0.1241 | 0.1039 | 0.1091 |
| APC | 0.1074 | 0.1203 | 0.0963 | 0.1067 |
| M5 | 0.1007 | 0.1042 | 0.0926 | 0.0926 |
| M6 | 0.1956 | 0.1039 | 0.1438 | 0.0902 |
| M7 | 0.1315 | 0.1066 | 0.1206 | 0.0918 |
| STAR | 0.1001 | 0.1038 | 0.0990 | 0.1005 |
| Ensemble | 0.0995 | 0.1010 | 0.0933 | 0.0929 |
| Model | NZ F (unadjusted) | NZ M (unadjusted) | NZ F (smoothed) | NZ M (smoothed) |
|---|---|---|---|---|
| LC | 0.0896 | 0.0796 | 0.0745 | 0.0690 |
| RH | 0.1081 | 0.0844 | 0.0946 | 0.0683 |
| APC | 0.0957 | 0.0876 | 0.0765 | 0.0696 |
| M5 | 0.0895 | 0.0802 | 0.0765 | 0.0714 |
| M6 | 0.1091 | 0.0885 | 0.0973 | 0.0691 |
| M7 | 0.1183 | 0.0869 | 0.1065 | 0.0716 |
| STAR | 0.0897 | 0.0851 | 0.0883 | 0.0839 |
| Ensemble | 0.0860 | 0.0730 | 0.0748 | 0.0641 |
As a final note regarding the differences between the raw and smoothed forecasts, we fit the LC, RH, APC, and STAR models to New Zealand data of years 1980-2011 and ages 0-90, where the forecasts are regularised by the Heligman-Pollard curve (constrained), being suitable for the entire age range. Compared to the results above, where the age range is 50-90 only, the age profile is more complex here and the results in Table 4 are mixed between the raw and smoothed forecasts. For the first three models, the smoothed forecasts yield lower errors, but the situation is reversed for the STAR model. Figure 8 compares the trajectories between the raw and smoothed STAR forecasts at ages 0 (birth), 20 (accident hump), and 65 (retirement). Despite their higher errors under the STAR model in this example, the smoothed forecasts still largely preserve the declining trends of the unadjusted forecasts. It appears that while the smoothing step improves age-profile coherence, its effect on forecasting accuracy generally depends on the data, model, and age range involved.
| Model | NZ F (unadjusted) | NZ M (unadjusted) | NZ F (smoothed) | NZ M (smoothed) |
|---|---|---|---|---|
| LC | 0.3380 | 0.2673 | 0.3174 | 0.2535 |
| RH | 0.4055 | 0.3046 | 0.3702 | 0.2881 |
| APC | 0.3198 | 0.2544 | 0.3118 | 0.2371 |
| STAR | 0.2853 | 0.2144 | 0.3111 | 0.2377 |
Figure 8: Observed log death rates (solid) in 1980-2011 and smoothed (dotted) and unadjusted (dashed) forecasted log death rates in 2012-2021 for New Zealand under STAR model.
5 Concluding remarks
In this paper, we have introduced the R package demofit, which provides a unified methodological and computational framework for mortality modelling by integrating parametric mortality laws and stochastic mortality models within a coherent interface. The package offers a comprehensive collection of mortality curves alongside widely used forecasting models. Parameter estimation for the mortality curves is carried out using multiple optimisation routines to enhance robustness, while computation remains efficient for typical mortality datasets. In addition, the ability to regularise mortality forecasts using parametric curves provides a practical mechanism for obtaining smooth and interpretable age profiles, particularly in insurance applications involving small or noisy policyholder populations. While similar functionality may be obtained by combining the few existing packages in a certain way, this requires substantial manual integration, inconsistent model specifications, and separate estimation frameworks. The demofit package instead provides a unified workflow with consistent estimation, model comparison, and forecast regularisation within a single interface. The package also includes ensemble forecasting functionality, allowing users to combine multiple models to improve predictive performance. In addition, it implements the common factor model for multi-population mortality modelling, enabling coherent forecasting across related populations.
We have described the underlying methodological and computational framework and illustrated the main functionalities of the package using real mortality data. Through a series of examples, we have demonstrated its use in mortality curve fitting, model comparison, forecasting, and the analysis of age-specific and time-varying mortality patterns. The results highlight the flexibility of the package in accommodating different modelling approaches within a single workflow, making it suitable for empirical research, actuarial applications, demographic studies, and teaching.
These examples also highlight that forecast regularisation / smoothing / graduation should be viewed as an optional post-processing step7. The additional backtesting results show that smoothing can reduce forecast errors in some settings, while producing higher errors in others. While this procedure does not necessarily alter the underlying temporal mortality trends, its effect depends on the selected mortality curve and the characteristics of the stochastic model. Users should therefore assess the effect of regularisation for their particular application by comparing the raw and smoothed forecasts and performing diagnostic checks8. The choice of mortality curve may also introduce systematic bias if its age profile is inappropriate for the population under study, and users are advised to fit and compare several mortality curves. Moreover, since the smoothed forecasts are obtained by transforming the original forecasts, the uncertainty measures generated by the underlying stochastic model should be interpreted with care rather than automatically regarded as uncertainty intervals for the smoothed forecasts.
There are several directions for future development. These include extending the range of models available within the package, such as incorporating other mortality and fertility curves and forecasting models, embedding life expectancy calculations, developing bootstrapping tools for uncertainty quantification, and building more multi-population models. In addition, further enhancements may include automated comparison and selection of models to assist users in identifying suitable specifications for different applications. We will continue to expand the functionality of the package to support a broader range of demographic and actuarial analyses.
6 Acknowledgements
The author thanks the editors and reviewers for their valuable comments and suggestions, which have greatly enhanced the presentation of the paper. The author also thanks Dr Emi Tanaka and Professor Rob Hyndman for their inspiration and insightful discussions at the initial stage of this R project. The author acknowledges support from the SMU Internal Research Grant (24-SOE-SMU-117).
6.1 CRAN packages used
demofit, MortalityLaws, forecast, StMoMo, StanMoMo, BayesMortalityPlus
6.2 CRAN Task Views implied by cited packages
ActuarialScience, Econometrics, Environmetrics, Finance, MissingData, TimeSeries
As noted in Renshaw and Haberman (2006), for the LC model, the maximum likelihood estimates from the iterative Newton updating scheme under the Gaussian error structure are numerically the same as the SVD estimates.↩︎
For applications where exposure information is available and likelihood-based estimation is preferred, demofit also provides Poisson estimation for selected models.↩︎
The last letter “S” refers to the feature that the mortality forecasts can be “smoothed”.↩︎
demofit initialises the RH model from the classical LC decomposition before optimisation, which empirically tends to retain the principal period trend while allowing the cohort effect to capture residual cohort-specific variation.↩︎
In the latest demofit version 0.1.5, if the
curveargument is missing, no smoothing will be performed.↩︎Note that the post-model regularisation adopted here is deliberately a pragmatic (and rather subjective) actuarial device rather than a replacement for statistically integrated smoothing approaches. In insurance practice, graduation may also be applied to modelled mortality rates when the resulting schedules require further smoothing, with the choice of graduation method depending on the data quality and intended use of the resulting mortality table (Atkinson 2018). Overall, this package function offers a simple and flexible practical option that can be evaluated empirically for a particular application. Further investigation of its statistical properties is left for future research.↩︎