Age-Period Modeling of Mortality Gaps: The Cases of Cancer and Circulatory Diseases

Abstract

Understanding and modeling mortality patterns, especially differences in mortality rates between populations, is vital for demographic analysis and public health planning. We compare three statistical models within the age-period framework to examine differences in death counts. The models are based on the double Poisson, bivariate Poisson, and Skellam distributions, each of which provides unique strengths in capturing underlying mortality trends. Focusing on mortality data from 1960 to 2015, we analyze the two leading causes of death in Italy, which exhibit significant temporal and age-related variations. Our results reveal that the Skellam distribution offers superior accuracy and simplicity in capturing mortality differentials. These findings highlight the potential of the Skellam distribution for analyzing mortality gaps effectively.

1. Introduction

Since the early nineteenth century, life expectancy in developed countries has shown consistent growth, particularly since the end of the Second World War (Oeppen and Vaupel 2002). Researchers, governments, and related organizations have often noted a faster-than-expected increase in longevity, posing financial challenges for life insurers, pension plans, and social security systems. A shift in longevity trends can substantially affect the number of elderly individuals and centenarians, thereby increasing demand for social services. Consequently, official forecasts that guide future pensions, healthcare, and other social needs can shape individual decisions on savings and retirement timing. In this context, methodological advancements in longevity studies have become increasingly important for anticipating future changes (see, e.g., Pascariu et al. 2018; Raftery et al. 2013; Torri and Vaupel 2012), aiming to reduce the gap between actual outcomes and projections.

Various approaches have been utilized to model the mortality surface, capturing how mortality rates evolve over time. Prior to the 1980s, mortality models were relatively simple and often depended on expert judgment (see Pollard 1987, for a detailed review). However, with the increased availability of mortality data and the advancement of statistical methods, more sophisticated and complex mortality models have since been developed.

According to Booth and Tickle (2008), three main paradigms in demographic modeling can be identified. The first paradigm, explanation, is based on structural or epidemiological models that focus on specific causes of death, such as a well-established link between lung cancer and tobacco smoking. The second paradigm, expectation, relies on expert opinions and may range from informal assessments to more formalized approaches. The third paradigm, extrapolation, utilizes regular patterns observed in age-specific mortality rates and temporal trends. This approach includes complex stochastic models, such as the Lee and Carter’s (1992) model and the generalized age-period-cohort model. Although the Lee–Carter model is regarded as one of the most influential in mortality modeling and forecasting, various extensions and modifications have gained considerable attention in actuarial and demographic literature (see, e.g., Brouhns et al. 2002; Cairns et al. 2006, 2009; Haberman and Renshaw 2011; Marino et al. 2023; Nigri et al. 2019; Shang and Haberman 2020; Shang et al. 2023).

Despite the widespread use of the Lee–Carter model and its variants, which serve as benchmarks for many new methodologies, these models have several limitations. Cairns et al. (2008) addressed the challenge of optimal mortality estimation by proposing criteria for a robust mortality model, emphasizing consistency with historical data and biologically plausible long-term dynamics. Apart from modeling age-specific mortality rates, research has shifted toward developing statistical frameworks that model the age-at-death distribution, allowing for parameter forecasting using econometric models (see, e.g., Aliverti et al. 2022; Basellini and Camarda 2019; Cardillo et al. 2023; Pascariu et al. 2019). However, researchers and practitioners have noted that separate forecasts for different populations or sub-populations can lead to unreliable scenarios, as backtesting often reveals diverging behaviors. Such inconsistencies in forecasts are implausible and may undermine decision-making processes.

Building on this discussion, we propose a novel approach for modeling and forecasting the gap in mortality counts between two populations. Specifically, we leverage recent developments in cause-of-death analysis (see, e.g., Arnold and Sherris 2015; Kjærgaard et al. 2019; Wilmoth 1995) by treating the two target populations as distinct causes of death within the same country. These two causes of death have implications for public health, particularly in identifying trends and temporal changes across different age groups. Understanding the trend associated with any specific cause of death enables policymakers to better allocate and target healthcare resources.

Our focus, therefore, is on the distribution of deaths, highlighting the disparities between two population counts, which represent the two leading causes of death, namely cancer and circulatory diseases, in the Italian population. These causes are major contributors to overall mortality. It is noteworthy that the two populations could also represent other categories, such as different genders, countries, or sub-populations within the same country. We aim to provide a more reliable and coherent forecast by addressing the inconsistencies often seen in separate forecasts for different populations or sub-populations. We develop a statistical toolkit for those working with mortality data, allowing for the joint modeling of mortality evolution between two populations while considering their differences. Alongside traditional alternatives like the double and bivariate Poisson distributions, we introduce the Skellam distribution, commonly used in fields involving count data, such as sport analysis.

Recently, Kung et al. (2021) introduced the Skellam distribution within the Lee-Carter framework to model the log-difference of mortality rates, focusing on age-dependent improvements in mortality and their implications for pricing longevity derivatives. While their work provides valuable insights into the dynamics of mortality improvement, our approach differs fundamentally in both scope and methodology. Specifically, we employ the Skellam distribution within the age-period framework to model differences in mortality counts across causes of death. In this context, it is standard practice to fix the age profile parameters, as the primary focus lies on capturing temporal (period) effects. This distinction highlights the originality of our contribution, which centers on analyzing mortality gaps with a straightforward and coherent structure tailored to demographic and public health applications.

Mortality data are commonly reported as count data within specific age-period domains for various populations. The simplest models for modeling mortality patterns use age and period as the only explanatory variables, since they capture all the necessary information intrinsic to the data. This approach serves as a foundational basis for exploring and understanding the complex dynamics of mortality across different populations.

The remainder of the paper is structured as follows: Section 2 introduces the motivating dataset and demographic measures. Sections 3 and 4 detail the methods used in our study. In Section 5, we present the implementation and evaluation of the distributions considered. Section 6 concludes the paper and offers some insights on potential extensions of the presented methodology.

2. Data and Measures

Let $\mathrm{x}$ represent age, $\{{\mathrm{t}}_1,{\mathrm{t}}_2,\dots ,\mathrm{T}\}$ denote a set of calendar years, and $\{\mathrm{A},\mathrm{B}\}$ represent the set of populations we consider. The death count at age $\mathrm{x}$ and time $\mathrm{t}$ for population $\mathrm{i}$ is denoted by ${\mathrm{D}}_{\mathrm{xt}}^{\mathrm{i}}$ .

We utilize period life tables by sex and single year of age (0–85+) from the Human Mortality Database (2024). We extract death counts by cause of death from the World Health Organization (2022) for the total population. This allows us to compute the proportion of deaths by cause, age group (in five-year intervals), and sex for each year from 1960 to 2015.

The data, in their original form, were classified using the seventh, eighth, ninth, and tenth revisions of the International Classification of Diseases (ICD). To analyze factors influencing Italian longevity, we group deaths into two major categories up to age 85+ (see Table 1): cancers and circulatory system diseases. These categories represent significant causes of death, capturing a substantial proportion of overall mortality. For details on ICD codes and the classification system, refer to Table 1.

Table 1.

ICD Codes and Our Classification into Two Major Groups.

Cause	ICD7	ICD8	ICD9	ICD10
Cancer	A044-A059, B018-B019	A045-A060	B08-B09, B100,B101, B109, B11-B17	C
Circulatory	A070, A079-A086,B022, B024-B029	A080-A088	B25-B30	I00-I99

Source. Authors' own elaborations on data from the World Health Organization (2022).

We define the gap in deaths between two populations, in our case the two causes of death $\mathrm{i}$ , at age $\mathrm{x}$ and year $\mathrm{t}$ as:

$${\mathrm{G}}_{\mathrm{xt}}={\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}-{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}.$$

Figure 1 provides a detailed depiction of the relationship between age and the temporal evolution of dynamics for a specific cause of death. It also illustrates the measure we are modeling: the difference in counts over time.

Figure 1.

Trends in cancer and circulatory death counts, as well as their differences, by age and year in Italy (both sexes combined) from 1960 to 2018.

3. Models and Distributions

3.1. Age-Period Model

Statistical modeling in longevity analysis is essential for understanding its evolution and variability. The literature highlights three key components—age, period, and cohort effects—that provide valuable insights into non-observable causal mechanisms.

Age effect is associated with a specific population’s biological and social processes, linked to its age structure rather than the period. In contrast, the period effect results from external factors that impact all age groups equally at a particular calendar time. These effects may arise from various environmental, social, or economic factors; the Spanish flu, the First World War, and COVID-19 are prime examples. Lastly, cohort effect refers to variations that stem from a group of individuals’ unique experiences or exposures as they progress through time.

The model that integrates all these components is known as the age-period-cohort approach, which presents technical challenges that can complicate the interpretation and clarity of estimates, especially when working with causes of death (see, e.g., Mason et al. 1973, for more details).

The demographic literature centering on longevity risks has primarily proposed age-period models. Continuing this line of research, we focus on the age-period framework. A canonical regression model, based on the ordinary least squares method, is specified as follows:

$${\mathrm{Y}}_{\mathrm{xt}}=\mu +{\beta }_{\mathrm{x}}^{\mathrm{A}}+{\beta }_{\mathrm{t}}^{\mathrm{P}}+{\mathrm{\epsilon }}_{\mathrm{xt}},$$

(1)

where ${\mathrm{Y}}_{\mathrm{xt}}$ represents the outcome variable to be explained, $\mu$ is the intercept, ${\beta }_{\mathrm{x}}^{\mathrm{A}}$ denotes the $\mathrm{x}$ th age effect, ${\beta }_{\mathrm{t}}^{\mathrm{P}}$ represents the ${\mathrm{t}}^{\mathrm{th}}$ period effect, and ${\epsilon }_{\mathrm{xt}}$ is the error term. In both age-period and age-period-cohort models, age and period are treated as categorical variables to capture nonlinear effects, utilizing a generalized linear model approach.

Given the mortality data for two populations, A and B, specific to ages and years, our goal is to determine the average period trend and the age-specific profile of the differences in mortality. We leverage the fundamental assumption in mortality modeling that death counts for each population, $({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})$ , are generated by a Poisson process (Basellini and Camarda 2019; Brillinger 1986) with distinct intensity parameters ${\lambda }_{\mathrm{xt}}^{\mathrm{A}}$ and ${\lambda }_{\mathrm{xt}}^{\mathrm{B}}$ :

$${\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}~\mathrm{Pois}\left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right),{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}~\mathrm{Pois}\left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right).$$

3.2. Poisson Distribution

3.2.1. Double Poisson

The simplest model employs two independent Poisson variables, where the parameters are constructed as the sum of age effects ( ${\beta }_{\mathrm{x}}^{\alpha }$ ) and period effects $({\beta }_{\mathrm{t}}^{\pi }$ ). We refer to this model as the double Poisson. It involves two distinct models for the parameters of the two Poisson processes:

$$\begin{array}{c}\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right)={\mu }^{\mathrm{A}}+{\beta }_{\mathrm{x}}^{\mathrm{A},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{A},\pi }\hfill \\ \multicolumn{1}{c}{\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right)={\mu }^{\mathrm{B}}+{\beta }_{\mathrm{x}}^{\mathrm{B},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{B},\pi }}\end{array}$$

If the goal is to estimate the differences in deaths between the two populations by year and age, one can calculate:

$${\hat{\mathrm{G}}}_{\mathrm{xt}}={\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{A}}-{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{B}},$$

where ${\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{A}}$ and ${\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{B}}$ are the parameter estimates.

The joint probability of the observations $({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})$ is:

$$\mathrm{f}({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}};{\lambda }_{\mathrm{xt}}^{\mathrm{A}},{\lambda }_{\mathrm{xt}}^{\mathrm{B}})=\frac{{\mathrm{e}}^{-{\lambda }_{\mathrm{xt}}^{\mathrm{A}}}{\left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right)}^{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}}}{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}!}·\frac{{\mathrm{e}}^{-{\lambda }_{\mathrm{xt}}^{\mathrm{B}}}{\left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right)}^{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}}}{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}!}.$$

Given observations $({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})$ for all age and period combinations $(x,t)$ , the log-likelihood function denoted by $\ell$ is:

$$\ell =\sum _{\mathrm{x},\mathrm{t}}[-{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}\log ({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}!)-{\lambda }_{\mathrm{xt}}^{\mathrm{A}}+{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right)-{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}\log ({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}!)-{\lambda }_{\mathrm{xt}}^{\mathrm{B}}+{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right)].$$

3.2.2. Bivariate Poisson

The double Poisson model, while straightforward, makes a somewhat unrealistic assumption by treating the processes as generated by two independent Poisson distributions. In cases where we consider two subpopulations within a single population (such as two sexes or two causes of death) or two populations experiencing similar environmental and external factors (such as neighboring countries), the correlation between deaths in the two groups can be quite high.

An alternative to the double Poisson model is the bivariate Poisson distribution, which accounts for the dependence between the two outcome variables. In this model, while the marginal distributions remain Poisson, the random variables are allowed to be dependent (Karlis and Ntzoufras 2005).

Consider random variables ${\mathrm{X}}_{\mathrm{k}}$ , $\mathrm{k}=1,2,3,$ which follow independent Poisson distributions with parameters ${\lambda }_{\mathrm{k}}>0$ . Thus, the random variables $\mathrm{X}={\mathrm{X}}_1+{\mathrm{X}}_3$ and $\mathrm{Y}={\mathrm{X}}_2+{\mathrm{X}}_3$ follow jointly a bivariate Poisson distribution:

$$(\mathrm{X},\mathrm{Y})~\mathrm{BP}({\lambda }_1,{\lambda }_2,{\lambda }_3)$$

Marginally, each random variable follows a Poisson distribution:

$$\begin{array}{c}\mathrm{X}~\mathrm{Pois}({\lambda }_1+{\lambda }_3),\mathrm{Y}~\mathrm{Pois}({\lambda }_2+{\lambda }_3)\hfill \end{array}$$

Moreover, $\mathrm{Cov}(\mathrm{X},\mathrm{Y})={\lambda }_3$ , indicating that ${\lambda }_3$ measures the dependence between the two random variables. If ${\lambda }_3=0$ , the bivariate Poisson distribution reduces to the double Poisson model.

In the case of the death count of two populations, we have:

$$({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})~\mathrm{BP}({\lambda }_{\mathrm{xt}}^{\mathrm{A}},{\lambda }_{\mathrm{xt}}^{\mathrm{B}},{\lambda }_{\mathrm{xt}}^{\mathrm{AB}})$$

Therefore, the joint probability mass function is given by:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{f}({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}};{\lambda }_{\mathrm{xt}}^{\mathrm{A}},{\lambda }_{\mathrm{xt}}^{\mathrm{B}},{\lambda }_{\mathrm{xt}}^{\mathrm{AB}})={\mathrm{e}}^{-({\lambda }_{\mathrm{xt}}^{\mathrm{A}}+{\lambda }_{\mathrm{xt}}^{\mathrm{B}}+{\lambda }_{\mathrm{xt}}^{\mathrm{AB}})}\frac{{\left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right)}^{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}}}{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}}\frac{{\left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right)}^{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}}}{{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}}}\\ \multicolumn{1}{c}{\sum _{\mathrm{k}=0}^{\min (\underset{\mathrm{xt}}{\stackrel{\mathrm{A}}{\mathrm{D}}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})}\left(\begin{array}{c}{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}\hfill \\ \multicolumn{1}{c}{\mathrm{k}}\end{array}\right)\left(\begin{array}{c}{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}\hfill \\ \multicolumn{1}{c}{\mathrm{k}}\end{array}\right)\mathrm{k}!{\left(\frac{{\lambda }_{\mathrm{xt}}^{\mathrm{AB}}}{{\lambda }_{\mathrm{xt}}^{\mathrm{A}}{\lambda }_{\mathrm{xt}}^{\mathrm{B}}}\right)}^{\mathrm{k}}.}\end{array}$$

Given observations $({\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})$ for all age and period combinations $(x,t)$ , the log-likelihood function is:

$$\ell =\sum _{\mathrm{x},\mathrm{t}}[-({\lambda }_{\mathrm{xt}}^{\mathrm{A}}+{\lambda }_{\mathrm{xt}}^{\mathrm{B}}+{\lambda }_{\mathrm{xt}}^{\mathrm{AB}})+\log (\sum _{\mathrm{k}=0}^{\min ({{\mathrm{D}}^{\mathrm{A}}}_{\mathrm{xt}},{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}})}\left(\begin{array}{c}{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{A}}\hfill \\ \multicolumn{1}{c}{\mathrm{k}}\end{array}\right)\left(\begin{array}{c}{\mathrm{D}}_{\mathrm{xt}}^{\mathrm{B}}\hfill \\ \multicolumn{1}{c}{\mathrm{k}}\end{array}\right)\mathrm{k}!{\left(\frac{{\lambda }_{\mathrm{xt}}^{\mathrm{AB}}}{{\lambda }_{\mathrm{xt}}^{\mathrm{A}}{\lambda }_{\mathrm{xt}}^{\mathrm{B}}}\right)}^{\mathrm{k}})].$$

In our modeling using the bivariate Poisson distribution, we explore three distinct models, each concentrating on a different parameter.

$$\begin{array}{c}\multicolumn{1}{c}{\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{A}}\right)={\mu }^{\mathrm{A}}+{\beta }_{\mathrm{x}}^{\mathrm{A},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{A},\pi }}\\ \multicolumn{1}{c}{\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{B}}\right)={\mu }^{\mathrm{B}}+{\beta }_{\mathrm{x}}^{\mathrm{B},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{B},\pi }}\\ \multicolumn{1}{c}{\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{AB}}\right)={\mu }^{\mathrm{AB}}}\end{array}$$

Following the work of Karlis and Ntzoufras (2005), we consider models with a constant ${\lambda }_{\mathrm{xt}}^{\mathrm{AB}}$ (without covariates on ${\lambda }_{\mathrm{xt}}^{\mathrm{AB}}$ ) for ease of interpretation.

In the estimation of the mortality gap, we have simply:

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathrm{G}}}_{\mathrm{xt}}={\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{A}}+{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{AB}}-{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{B}}-{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{AB}}}\\ \multicolumn{1}{c}{{\hat{\mathrm{G}}}_{\mathrm{xt}}={\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{A}}-{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{B}}.}\end{array}$$

3.3. Skellam Distribution

By employing the bivariate Poisson distribution to estimate the gap between deaths, the correlation component is eliminated. We propose using the Skellam (or Poisson Difference) distribution to directly model the gap rather than jointly modeling the number of deaths for the two populations.

The Poisson Difference distribution, introduced by Skellam (1946), describes the difference $\mathrm{Z}=\mathrm{X}-\mathrm{Y}$ between two discrete random variables $\mathrm{X}$ and $\mathrm{Y}$ . The probability function of $\mathrm{Z}$ is a discrete distribution defined over the set of integers $\mathbb{Z}$ , encompassing both positive and negative values. This distribution applies to any pair of random variables $(\mathrm{X},\mathrm{Y})$ that can be expressed as $\mathrm{X}={\mathrm{W}}_1+{\mathrm{W}}_3$ and $\mathrm{Y}={\mathrm{W}}_2+{\mathrm{W}}_3$ , where ${\mathrm{W}}_1~\mathrm{Pois}\left({\lambda }_1\right)$ and ${\mathrm{W}}_2~\mathrm{Pois}\left({\lambda }_2\right)$ , with ${\mathrm{W}}_3$ following any distribution. The joint distribution of $\mathrm{X}$ and $\mathrm{Y}$ forms a bivariate distribution, with correlation arising from the shared stochastic component ${\mathrm{W}}_3$ (Karlis and Ntzoufras 2009). The mean of a Skellam distribution is given by $\mathbb{E}[\mathrm{Z}]={\lambda }_1-{\lambda }_2$ .

The probability mass function of $Z$ is defined as $\mathrm{Pr}(\mathrm{Z}=\mathrm{z})=\mathrm{f}(\mathrm{z};{\lambda }_1,{\lambda }_2)$ with,

$$\mathrm{f}(\mathrm{z};{\lambda }_1,{\lambda }_2)={\mathrm{e}}^{-{\lambda }_1-{\lambda }_2}{\left(\frac{{\lambda }_1}{{\lambda }_2}\right)}^{\frac{\mathrm{z}}{2}}{\mathrm{I}}_{\mathrm{z}}\left(2\sqrt{{\lambda }_1{\lambda }_2}\right),$$

(2)

where ${\mathrm{I}}_{\mathrm{z}}\left(\nu \right)={\left(\frac{\nu }{2}\right)}^{\mathrm{z}}\sum _{\mathrm{k}=0}^{\infty }\frac{{\left(\frac{{\nu }^2}{4}\right)}^{\mathrm{k}}}{\mathrm{k}!·(\mathrm{z}+\mathrm{k})!}$ is the modified Bessel function of the first kind with order $z$ (see Abramowitz and Stegun 1968).

Therefore, we assume that ${\mathrm{G}}_{\mathrm{xt}}~\mathrm{PD}({\lambda }_{\mathrm{xt}}^{\mathrm{C}},{\lambda }_{\mathrm{xt}}^{\mathrm{D}})$ , and, as before, we model it using age and period as covariates:

$$\begin{array}{c}\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{C}}\right)={\mu }^{\mathrm{C}}+{\beta }_{\mathrm{x}}^{\mathrm{C},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{C},\pi },\hfill \\ \multicolumn{1}{c}{\log \left({\lambda }_{\mathrm{xt}}^{\mathrm{D}}\right)={\mu }^{\mathrm{D}}+{\beta }_{\mathrm{x}}^{\mathrm{D},\alpha }+{\beta }_{\mathrm{t}}^{\mathrm{D},\pi },}\end{array}$$

We use the superscripts ^$\mathrm{C}$ and ^$\mathrm{D}$ to denote the two parameters because Skellam’s parameters ${\lambda }_{\mathrm{xt}}^{\mathrm{C}}$ and ${\lambda }_{\mathrm{xt}}^{\mathrm{D}}$ do not represent the intensities of deaths for populations A and B directly. Instead, they represent two intensities that together approximate the differences in deaths.

The fitted values from the regression directly provide estimates of the death gap:

$${\hat{\mathrm{G}}}_{\mathrm{xt}}={\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{C}}-{\hat{\mathrm{\lambda }}}_{\mathrm{xt}}^{\mathrm{D}}.$$

The Skellam distribution requires only the series of differences in deaths between populations A and B, rather than the two death series individually.

We estimate the model parameters by maximizing the log-likelihood function of Skellam distribution, which is

$$\begin{array}{c}\ell (\mu ,{\beta }_{\mathrm{x}},{\beta }_{\mathrm{t}};{\mathrm{G}}_{\mathrm{xt}})=\sum _{\mathrm{x},\mathrm{t}}[-({\lambda }_{\mathrm{xt}}^{\mathrm{C}}+{\lambda }_{\mathrm{xt}}^{\mathrm{D}})+\frac{{\mathrm{G}}_{\mathrm{xt}}}{2}\log \left(\frac{{\lambda }_{\mathrm{xt}}^{\mathrm{C}}}{{\lambda }_{\mathrm{xt}}^{\mathrm{D}}}\right)+\log \left({\mathrm{I}}_{{\mathrm{G}}_{\mathrm{xt}}}\left(2\sqrt{{\lambda }_{\mathrm{xt}}^{\mathrm{C}}{\lambda }_{\mathrm{xt}}^{\mathrm{D}}}\right)\right)].\hfill \end{array}$$

3.4. Estimation and Goodness of Fit

With regard to the estimation procedure, Ho and Singer (2001) proposes a generalized least squares method, whereas Kocherlakota and Kocherlakota (2001) outlines a Newton–Raphson approach for likelihood maximization in the context of bivariate Poisson regression models. In this application we leverage the approach by Karlis and Ntzoufras (2003) using the EM algorithm to obtain maximum likelihood estimates for both, double and bivariate Poisson. The EM algorithm proceeds by estimating the unobserved data via their conditional expectations at the E-step and then it maximizes the complete-data likelihood at the M-step (Karlis and Ntzoufras 2005).

The estimation of the parameters for the Skellam regression is performed by minimizing the negative log-likelihood using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, using the optim function in .

The parameter vector $\mathbf{\theta }=\{{\mu }^{\mathrm{C}},{\mu }^{\mathrm{D}},{\beta }_{\mathrm{x}}^{\mathrm{C}},{\beta }_{\mathrm{t}}^{\mathrm{D}},{\beta }_{\mathrm{x}}^{\mathrm{D}},{\beta }_{\mathrm{t}}^{\mathrm{D}}\}$ includes the intercepts ( ${\mu }^{\mathrm{C}},{\mu }^{\mathrm{D}}$ ), age-specific effects ( ${\beta }_{\mathrm{x}}^{\mathrm{C}},{\beta }_{\mathrm{x}}^{\mathrm{D}}$ ), and period-specific effects ( ${\beta }_{\mathrm{t}}^{\mathrm{C}},{\beta }_{\mathrm{t}}^{\mathrm{D}}$ ). The optimization is carried out iteratively, where at each step the parameters are updated according to the BFGS algorithm:

$${\mathbf{\theta }}^{(\mathrm{k}+1)}={\mathbf{\theta }}^{\left(\mathrm{k}\right)}-{\mathrm{H}}^{-1}\nabla (-\ell ),$$

where ${\mathrm{H}}^{-1}$ is the approximated inverse Hessian matrix, and $\nabla (-\ell )$ represents the gradient of the negative log-likelihood with respect to the parameters.

The optimization terminates when the norm of the gradient satisfies a convergence criterion, such as $\Vert \nabla (-\ell )\Vert <ϵ$ , where $ϵ$ is a predefined tolerance. This approach ensures accurate estimation of parameters while leveraging the full structure of the Skellam distribution.

In order to measure the goodness-of-fit in terms of the log-likelihood, Akaike information criterion (AIC), Akaike information criterion with small-sample correction ( ${\mathrm{AIC}}_{\mathrm{c}}$ ), and Bayesian information criterion (BIC) are used.

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{AIC}=2\mathrm{p}-2\ell }\\ \multicolumn{1}{c}{{\mathrm{AIC}}_{\mathrm{c}}=\mathrm{AIC}+\frac{2{\mathrm{p}}^2+2\mathrm{p}}{\mathrm{n}-\mathrm{p}-1}}\end{array}$$

and

$$\mathrm{BIC}=\mathrm{p}\log \mathrm{n}-2\ell ,$$

where p is the number of effective parameters, and n is the number of total observations. All the results are presented in Tables 2 and 3 of Section 5.

Table 2.

Accuracy Evaluation of the Three Models in the In-Sample and Out-of-Sample Fitting over Three Baseline and Forecasting Periods for All the Ages. Bold and Underlined Numbers Represent the Lowest and Second Lowest Values Among the Three Models for Each Fitting Period, Respectively.

Model	In-sample						Out-of-sample
	BIC	AIC	AICc	RMSE	MAE	MAPE	RMSE	MAE	MAPE
Baseline 1961–2000, Forecast 2001–2015
Skellam	9,884.40	9,377.92	9,422.57	466.79	235.49	0.89	1,624.31	778.76	1.26
Double Poisson	299,456.57	298,950.10	298,994.74	1,557.18	813.07	0.98	2,627.33	1,700.54	2.95
Bivariate Poisson	292,884.12	292,294.79	292,340.31	1,539.02	791.50	0.85	2,593.86	1,666.44	2.49
Baseline 1971–2000, Forecast 2001–2015
Skellam	7,305.39	6,915.83	6,956.86	495.35	238.58	0.53	1,475.54	803.09	0.92
Double Poisson	185,354.15	184,964.59	185,005.62	1,380.38	748.31	0.76	2,194.91	1,413.12	2.05
Bivariate Poisson	184,891.00	184,432.74	184,474.77	1,375.31	740.18	0.71	2,182.50	1,398.93	1.87
Baseline 1981–2000, Forecast 2001–2015
Skellam	4,765.96	4,490.27	4,529.64	520.92	233.36	0.29	1,406.80	865.64	0.42
Double Poisson	91,332.70	91,057.01	91,096.38	1,124.21	589.59	0.47	1,691.55	1,098.81	1.18
Bivariate Poisson	91,316.54	90,986.42	91,027.04	1,122.08	585.60	0.44	1,687.16	1,090.11	1.04

Table 3.

Accuracy Evaluation of the Three Models, Both in-Sample and Out-of-Sample, Across Three Baseline and Forecasting Periods for Ages 40 and Older. Bold and Underlined Numbers, Respectively, Represent the Lowest and Second Lowest Values Among the Three Models for Each Fitting Period.

Model	In-sample							Out-of-sample
	BIC	AIC	AICc	RMSE	MAE	MAPE	RMSE	MAE	MAPE
Baseline 1961–2000,Forecast 2001–2015, Age $\ge 40$
Skellam	9,452.07	9,079.01	9,149.82	601.12	372.36	0.18	2,391.84	1,561.41	0.27
Double Poisson	258,270.67	257,897.60	257,968.42	2,085.50	1,379.95	0.57	3,534.51	2,950.76	0.87
Bivariate Poisson	229,753.76	229,309.58	229,382.14	1,952.41	1,186.06	0.40	3,109.66	2,474.85	0.61
Baseline 1971–2000,Forecast 2001–2015, Age $\ge 40$
Skellam	4,439.31	4,165.83	4,226.47	633.52	380.14	0.17	1,350.22	890.35	0.20
Double Poisson	169,287.45	169,013.97	169,074.61	1,864.26	1,280.84	0.42	2,956.41	2,469.65	0.72
Bivariate Poisson	156,421.77	156,091.32	156,153.88	1,783.47	1,126.53	0.34	2,634.37	2,114.51	0.53
Baseline 1981–2000,Forecast 2001–2015, Age $\ge 40$
Skellam	4,089.25	3,910.44	3,962.35	677.45	408.04	0.12	1,156.49	931.80	0.21
Double Poisson	85,926.79	85,747.99	85,799.89	1,529.94	1,050.06	0.30	2,327.94	1,973.18	0.54
Bivariate Poisson	81,892.53	81,671.02	81,725.22	1,485.71	964.05	0.25	2,183.59	1,771.90	0.41

4. Forecasting

Introducing a novel forecasting model for the differences in mortality between two distinct populations is crucial for several reasons. It allows for generating coherent forecasts for both populations, ensuring consistency in predictive analysis. By leveraging the Skellam distribution, which naturally models the difference between two Poisson-distributed variables, we can more accurately capture the inherent variability and interdependence between the causes of death. This approach enhances the reliability of forecasts and provides a deeper understanding of the underlying mortality dynamics, ultimately supporting more informed public health decisions and better resource allocation to address these gaps.

Given the age and period effects estimated during the baseline period using the regression models discussed in the previous subsections, the proposed forecasting method involves fixing the observed age profile ${\beta }_{\mathrm{x}}^{\alpha }$ and the intercept $\mu$ for the two populations. We then forecast the period effects $\mathrm{h}$ years ahead ${\beta }_{\mathrm{T}+\mathrm{h}}^{\pi }$ . This forecasting approach is applied uniformly across all models presented in Section 3.

We aim to predict the differences in mortality between the two populations $h$ years into the future:

$${\hat{\mathrm{G}}}_{\mathrm{x},\mathrm{T}+\mathrm{h}}={\hat{\mathrm{\lambda }}}_{\mathrm{x},\mathrm{T}+\mathrm{h}}^{\mathrm{A}}-{\hat{\mathrm{\lambda }}}_{\mathrm{x},\mathrm{T}+\mathrm{h}}^{\mathrm{B}},$$

where

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathrm{\lambda }}}_{\mathrm{x},\mathrm{T}+\mathrm{h}}^{\mathrm{A}}=\exp ({\mu }^{\mathrm{A}}+{\beta }_{\mathrm{x}}^{\mathrm{A},\mathrm{\alpha }}+{\hat{\mathrm{\beta }}}_{\mathrm{T}+\mathrm{h}}^{\mathrm{A},\pi })}\\ \multicolumn{1}{c}{{\hat{\mathrm{\lambda }}}_{\mathrm{x},\mathrm{T}+\mathrm{h}}^{\mathrm{B}}=\exp ({\mu }^{\mathrm{B}}+{\beta }_{\mathrm{x}}^{\mathrm{B},\alpha }+{\hat{\mathrm{\beta }}}_{\mathrm{T}+\mathrm{h}}^{\mathrm{B},\pi })}\end{array}$$

We propose using a multivariate random walk model to project the historical period effects for the two populations. This approach is selected because it effectively captures the joint dynamics and correlations between the two causes of death over time. The multivariate random walk model assumes that future values evolve as a stochastic process with a deterministic drift and random fluctuations. This framework offers both flexibility and robustness in modeling temporal dependencies and inherent randomness in mortality data, resulting in more accurate and realistic forecasts.

We denote ${\beta }_{\mathrm{t}}^{\left(\mathrm{i}\right)}$ as the observed value of the period parameter (with superscript $\pi$ omitted for simplicity) for population $\mathrm{i}\in \{\mathrm{A},\mathrm{B}\}$ in year $\mathrm{t}\in \{{\mathrm{t}}_1,{\mathrm{t}}_2,\dots ,\mathrm{T}\}$ . Accordingly, we define the matrix of period parameters as:

$${\beta }_{\mathrm{t}}^{\left(\mathrm{i}\right)}=\left[\begin{array}{c}{\beta }_{{\mathrm{t}}_1}^{\left(\mathrm{A}\right)}{\beta }_{{\mathrm{t}}_2}^{\left(\mathrm{A}\right)}\dots {\beta }_{\mathrm{T}}^{\left(\mathrm{A}\right),}\hfill \\ \multicolumn{1}{c}{{\beta }_{{\mathrm{t}}_1}^{\left(\mathrm{B}\right)}{\beta }_{{\mathrm{t}}_2}^{\left(\mathrm{B}\right)}\dots {\beta }_{\mathrm{T}}^{\left(\mathrm{B}\right)}}\end{array}\right],$$

where the row indicates the specific population population and the column indicates the time.

Therefore, $({\beta }_t^{\left(i\right)})t\in \{t_1,\dots ,T\}$ represents the bivariate time series of period parameters estimated using the models in Section 3. We simultaneously model this time series using a multivariate random walk with drift:

$${\mathbf{\beta }}_{\mathbf{t}}^{\left(\mathbf{i}\right)}={\mathbf{\beta }}_{\mathbf{t}-1}^{\left(\mathbf{i}\right)}+{\mathit{\delta }}^{\left(\mathbf{i}\right)}+{ϵ}_{\mathbf{t}}^{\left(\mathbf{i}\right)},{ϵ}_{\mathbf{t}}^{\left(\mathbf{i}\right)}~\mathrm{N}(0,{\mathbf{\Sigma }}^{\left(\mathrm{i}\right)})$$

(3)

where ${\mathit{\delta }}^{\left(\mathrm{i}\right)}$ is the vector of drift parameters that drive the dynamics for population $i$ , and ${\mathit{\Sigma }}^{\left(i\right)}$ is the $({\mathrm{n}}_{\mathrm{i}}\times {\mathrm{n}}_{\mathrm{i}})$ variance-covariance matrix of the multivariate white noise ${ϵ}_{\mathbf{t}}^{\left(\mathbf{i}\right)}$ .

Using the in-sample data, we compute $h$ -step-ahead point forecasts and assess out-of-sample accuracy by comparing these forecasts with the actual out-of-sample data. Let $\mathrm{T}$ denote the forecast origin and $h$ the forecast horizon. The $h$ -step-ahead forecast of the period parameter is given by:

$${\widehat{\beta }}_T^{\left(i\right)}\left(h\right)=E[{\beta }_{T+h}^{\left(i\right)}|{\beta }_T^{\left(i\right)},{\beta }_{t_{\tau }-1}^{\left(i\right)},\dots ,{\beta }_1^{\left(i\right)}.],$$

The $\mathrm{h}$ -step-ahead forecast for the multivariate random walk with drift model, as specified in Equation (3), can be expressed as:

$${\hat{\beta }}_{\mathbf{T}}^{\left(\mathbf{i}\right)}\left(\mathrm{h}\right)={\mathbf{\beta }}_{\mathbf{T}}^{\left(\mathbf{i}\right)}+{\mathit{\delta }}^{\left(\mathbf{i}\right)}\mathrm{h},$$

We evaluate the model’s performance using the mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE). Let $\mathrm{T}$ represent the final year of the forecast, and ${\mathrm{t}}_{\mathrm{\tau }}$ (where ${\mathrm{t}}_1<{\mathrm{t}}_{\mathrm{\tau }}<\mathrm{T}$ ) denote the calendar year corresponding to the last observation in the in-sample set.

We calculate the MAE, RMSE, and MAPE for the differences in mortality between the two populations $i$ as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{MAE}=\frac{\sum _{\mathrm{t}={\mathrm{t}}_{\mathrm{\tau }}+\mathrm{h}}^{\mathrm{T}}\mid {\mathrm{G}}_{\mathrm{t}}-{\hat{\mathrm{G}}}_{\mathrm{t}}\mid }{(\mathrm{T}-{\mathrm{t}}_{\mathrm{\tau }}-1)},}\\ \multicolumn{1}{c}{\mathrm{RMSE}=\sqrt{\frac{\sum _{\mathrm{t}={\mathrm{t}}_{\mathrm{\tau }}+\mathrm{h}}^{\mathrm{T}}{({\mathrm{G}}_{\mathrm{t}}-{\hat{\mathrm{G}}}_{\mathrm{t}})}^2}{(\mathrm{T}-{\mathrm{t}}_{\mathrm{\tau }}-1)}},}\\ \multicolumn{1}{c}{\mathrm{MAPE}=\frac{1}{(\mathrm{T}-{\mathrm{t}}_{\mathrm{\tau }}-1)}\sum _{\mathrm{t}={\mathrm{t}}_{\mathrm{\tau }}+\mathrm{h}}^{\mathrm{T}}\frac{|{\mathrm{G}}_{\mathrm{t}}-{\hat{\mathrm{G}}}_{\mathrm{t}}|}{{\mathrm{G}}_{\mathrm{t}}}\times 100,}\end{array}$$

where ${\mathrm{G}}_{\mathrm{t}}$ denotes the out-of-sample holdout data, and ${\hat{\mathrm{G}}}_{\mathrm{t}}$ denotes its estimates.

5. Application to Cancer and Circulatory Diseases

We applied the three models, treating age and years as categorical variables with the first category as the reference. This approach yielded two sets of age and period effects for each model, corresponding to the $\beta$ parameters introduced in Section 3.

We compare the accuracy of the three proposed models in estimating the differences in death counts ( ${\hat{\mathrm{G}}}_{\mathrm{xt}}$ ), both in-sample and out-of-sample. Specifically, we use RMSE, MAE, and MAPE for both in-sample and out-of-sample evaluations, and further assess goodness-of-fit using BIC, AIC, and AICc. To compare the predictive accuracy of the competing forecasts, we employ the Diebold-Mariano test.

Estimating and forecasting mortality is inherently challenging, particularly when introducing a novel approach. This complexity arises from the evolving nature of mortality dynamics, which has changed significantly over the past thirty to forty years due to both endogenous and exogenous factors, such as improvements in social and health conditions (Rau et al. 2008). Given this, each statistical model must be tested and calibrated to determine the period that best reflects the estimated parameters. In other words, it is crucial to identify the optimal period length and age range for feeding the model. These factors impact the estimation of model parameters and, consequently, the model’s reliability, as measured by backtesting errors.

This aspect is frequently underestimated in mortality studies. It is important to note that this issue is not limited to our model but is relevant to all mortality models, which should be subjected to a sensitivity analysis regarding period and age structure variations. To address this, we conduct a sensitivity analysis to evaluate how performance varies with different estimation period lengths and age ranges. For our analysis, we consider two distinct age ranges: (i) complete age range, from 0 to 80+ years; and (ii) a specific focus on the adult age range, from forty years onwards. For both age ranges, we examine three estimation periods: forty years (1961–2000), thirty years (1971–2000), and twenty years (1981–2000). All results are evaluated and compared over the same forecasting horizon, spanning the period 2001 to 2015, for consistency in the out-of-sample exercise.

The accuracy results, and the goodness-of-fit metrics related to the in-sample, are presented in Tables 2 and 3 for the full age range and ages 40+, respectively. The tables reveal several insights, the age-period model based on the Skellam distribution consistently provides superior in-sample estimates, demonstrating higher accuracy and superior goodness-of-fit across all periods and both age ranges.

In terms of out-of-sample estimation, the sensitivity analysis indicates that the model based on the Skellam distribution still outperforms the other two models. All the models have better performances on the more recent windows, which is not surprising because this is the period with the most regular longevity dynamics (Nigri et al. 2022).

To offer a more comprehensive analysis and enhance the understanding of the models’ behavior, we present a two-phase graphical investigation. In Figure 2, we provide heatmaps illustrating the out-of-sample RMSE across different age groups: young ages (0–24), adult ages (25–59), and elderly ages (60+). Additionally, scatterplots are presented to demonstrate the model fit, both in-sample and out-of-sample, over time for selected ages. The results for the period 1961 to 2000 are detailed below, with corresponding results for other periods available in Appendix.

Figure 2.

Root Mean Square Error (RMSE) of Forecasting Differences in Death Counts: Baseline 1961 to 2000; Ages 0 to 80+. Dark colors represent inferior forecasting accuracy.

Figure 2 shows that the models achieve better accuracy at younger ages, although even in this age group, the Skellam regression provides slightly more accurate forecasts. The accuracy advantage of the Skellam-based model becomes more pronounced at adult ages, as also confirmed in Figure 3. However, among older individuals, especially as the forecasting horizon approaches the last observations, the accuracy of the bivariate Poisson and double Poisson models surpasses that of the Skellam model. The two Poisson-based models exhibit fairly similar trends, though the bivariate Poisson model displays slightly lighter colors at adult and older ages.

Figure 3.

Root Mean Square Error of forecasting estimation for differences in death counts. Estimations are based on the baseline period 1961 to 2000 for ages 40 to 80+. Dark colors represent inferior forecasting accuracy.

Furthermore, we can observe from Figures 4 and 5, both in-sample and out-of-sample results demonstrate that the Skellam-based model provides an exceptional fit to the data compared to the other models, and importantly, it maintains this consistency over time. When the model is trained on a shorter baseline period (see Appendix), the accuracy decreases but remains outstanding for some ages, mainly in the in-sample fitting.

Figure 4.

In-sample and out-of-sample analysis for specific ages, with estimations based on the baseline period 1961 to 2000; ages 0 to 80+.

Figure 5.

In-sample and out-of-sample exercise for specific ages. Estimations based on the baseline period 1961 to 2000; ages 40 to 80+.

5.1. Diebold-Mariano (DM) Test

To compare the predictive accuracy of competing forecasts, we employ the DM test (Diebold and Mariano 1995). In essence, the forecast error ${\mathrm{e}}_{\mathrm{it}}$ is defined as the difference between the values predicted by model $\mathrm{i}$ and the actual observed values at time $\mathrm{t}$ . The loss associated with the forecast from model $i$ is assumed to be $\mathrm{g}\left({\mathrm{e}}_{\mathrm{it}}\right)={\mathrm{e}}_{\mathrm{it}}^2$ . The two forecasts are considered to have equal accuracy if and only if the expected value of the loss differential, ${\mathrm{d}}_{\mathrm{t}}=\mathrm{g}\left({\mathrm{e}}_{1\mathrm{t}}\right)-\mathrm{g}\left({\mathrm{e}}_{2\mathrm{t}}\right)$ , is zero for all $\mathrm{t}$ . The DM test statistic is defined as follows:

$$\mathrm{DM}=\frac{\overline{\mathrm{d}}}{\sqrt{\frac{\mathrm{s}}{\mathrm{N}}}},$$

(4)

where $\overline{\mathrm{d}}$ is the sample mean of the loss differential, $\mathrm{s}$ is the sample variance of the loss differential, and $\mathrm{N}$ represents the sample size. The null hypothesis of the test states that the competing models have equal forecast accuracy, that is, ${\mathrm{H}}_0:\mathrm{E}[{\mathrm{d}}_{\mathrm{t}}]=0,\forall \mathrm{t}$ , while the alternative hypothesis is ${\mathrm{H}}_1:\mathrm{E}[{\mathrm{d}}_{\mathrm{t}}]\ne 0,\forall \mathrm{t}$ . Under ${\mathrm{H}}_0$ , the DM statistic is asymptotically distributed as a standard normal distribution with a mean of 0 and a standard deviation of 1.

We apply the DM test to compare the forecast accuracy of the Skellam Model against the double Poisson and bivariate Poisson models. The results of the DM test are reported in Tables 4 and 5. When the null hypothesis is rejected at a certain level of significance, it indicates that the differences between the two models are statistically significant, and the forecasting accuracy of the first model (Skellam, in our study) is superior to that of the second model. From Tables, we observe that the forecasts produced by the Skellam approach consistently exhibit higher accuracy compared to those from the other two models.

Table 4.

The DM Test Between Skellam-Based Forecast and Double Poisson and Bivariate Poisson Models, Respectively. Evaluation Period 2001 to 2015 for Age Groups from 0 to 80+. Asterisks ***, **, and * Indicate Two-Side Significance at 0.1%, 1%, and 5% level, respectively.

Age group	Baseline 1961–2000		Baseline 1971–2000		Baseline 1981–2000
	Forecast 2001–2015		Forecast 2001–2015		Forecast 2001–2015
	Double Poisson	Bivariate Poisson	Double Poisson	Bivariate Poisson	Double Poisson	Bivariate Poisson
0–4	−12.014***	−10.614***	−9.687***	−8.694***	−6.527***	−6.09***
5–9	−14.374***	−14.457***	−13.877***	−14.092***	−9.473***	−9.251***
10–14	−12.071***	−11.882***	$-12.472$ ***	$-12.607$ ***	$-9.468$ ***	$-9.29$ ***
15–19	$-10.396$ ***	$-9.462$ ***	$-10.404$ ***	$-10.096$ ***	$-8.017$ ***	$-7.605$ ***
20–24	$-12.004$ ***	$-11.274$ ***	$-12.609$ ***	$-12.38$ ***	$-10.182$ ***	$-9.71$ ***
25–29	$-10.572$ ***	$-9.715$ ***	$-9.752$ ***	$-9.224$ ***	$-6.643$ ***	$-5.955$ ***
30–34	$-7.308$ ***	$-6.681$ ***	$-6.37$ ***	$-6.076$ ***	$-4.349$ **	$-4.08$ **
35–39	$-7.818$ ***	$-7.427$ ***	$-6.586$ ***	$-6.376$ ***	$-4.01$ **	$-3.76$ ***
40–44	$-7.525$ ***	$-7.289$ ***	$-6.701$ ***	$-6.602$ ***	$-4.381$ **	$-4.271$ **
45–49	$-8.928$ ***	$-8.692$ ***	$-9.031$ ***	$-8.953$ ***	$-5.694$ ***	$-5.613$ ***
50–54	$-7.728$ ***	$-7.591$ ***	$-8.17$ ***	$-8.125$ ***	$-8.028$ . ***	$-8.004$ ***
55–59	$-5.366$ ***	$-5.354$ ***	$-5.229$ ***	$-5.232$ ***	$-4.771$ ***	$-4.777$ ***
60–64	$-3.154$ **	$-3.253$ **	$-2.984$ *	$-3.026$ **	– $2.443$ *	– $2.479$ *
65–69	$-0.575$	$-0.637$	$2.185$ *	$1.95$	$5.742$ ***	$5.559$ ***
70–74	$-8.478$ ***	$-7.14$ ***	$-3.183$ **	$-2.876$ *	$2.194$ *	$2.281$ *
75–79	$-9.641$ ***	$-9.185$ ***	$-7.294$ ***	$-7.065$ ***	$-3.891$ **	$-3.773$ **
80+	$2.449$ *	$2.405$ *	$2.155$ *	$2.121$	$0.477$	$0.419$

Table 5.

The DM Test Between Skellam-Based Forecast and Double Poisson and Bivariate Poisson Models, Respectively. Evaluation Period 2001 to 2015 for Age Groups from 40 to 80+. Asterisks ***, **, and * Indicate Two-Side Significance at 0.1%, 1%, and 5% Level, Respectively.

Age group	Baseline 1961–2000		Baseline 1971–2000		Baseline 1981–2000
	Forecast 2001–2015		Forecast 2001–2015		Forecast 2001–2015
	Double Poisson	Bivariate Poisson	Double Poisson	Bivariate Poisson	Double Poisson	Bivariate Poisson
40–44	$-7.931$ ***	$-3.03$ **	$-7.209$ ***	$-1.056$	$-4.687$ ***	$4.388$ **
45–49	$-9.848$ ***	$-7.663$ ***	$-9.217$ ***	$-6.603$ ***	$-6.003$ ***	$0.118$
50–54	$-8.959$ ***	$-7.367$ ***	$-8.179$ ***	$-7.118$ ***	$-7.972$ ***	$-7.325$ ***
55–59	$-7.52$ ***	$-6.215$ ***	$-5.752$ ***	$-5.276$ ***	$-5.41$ ***	$-5.083$ ***
60–64	$-6.247$ ***	$-5.582$ ***	$-3.964$ **	$-3.938$ **	$-3.689$ **	$-3.661$ **
65–69	$2.697$ *	$3.057$ **	$-2.187$ *	$-2.405$ *	$-0.782$	$-1.334$
70–74	$-2.575$ *	$-0.63$	$-9.03$ ***	$-6.402$ ***	$-3.997$ **	$-1.568$
75–79	$-8.342$ ***	$-7.642$ ***	$-8.558$ ***	$-7.034$ ***	$-4.874$ ***	$-3.621$ **
80+	$2.576$ *	$2.144$	$1.733$	$0.938$	$-0.399$	$-1.824$

6. Conclusion

We propose three distinct modeling approaches to analyze death counts and their differences, each grounded in the fundamental assumption that deaths at each age and year follow a Poisson distribution. Our approaches include three age-period models: the classic double Poisson, the more advanced bivariate Poisson, and the novel Skellam distribution. These models represent a progression from traditional methods to more sophisticated and innovative techniques in studying mortality gaps.

We applied the three modeling approaches to analyze differences in death counts for two major causes of death in Italy: cancer and diseases of the circulatory system. By focusing on causes of death, we bypass the need to account for risk exposures, as the exposure populations are consistent for both causes, allowing us to directly compare mortality rates and death counts.

Our comparison of the three models reveals that the Skellam distribution model, an innovative approach in this context, provides the best fit to the data and requires the least amount of input information.

Given the complexity of mortality, which is influenced by numerous factors, no single statistical model can fully capture its dynamics. To address this, our study includes a sensitivity analysis to evaluate how performance varies with different estimation periods and age ranges.

Leveraging model selection criteria such as AICc and BIC in the in-sample study, and also exploiting both metric errors and hypothesis tests (DM test) for the out-of-sample exercise, we found that the Skellam distribution is valuable for jointly and consistently modeling mortality differences. Its versatility also allows for generalization to other contexts, such as comparing deaths between different sexes, countries, or economic conditions. Utilizing the Skellam distribution enables the development of coherent forecasting models for these scenarios.

All codes, data, and supplementary results are available on the online Repository in GitHub.

There are at least two ways in which our work can be further extended, and we briefly outline two: (1) Leveraging the Skellam-based framework to improve the forecasting of the death count of one of the two populations, such as one of the two causes. Combining forecasts for one of the two causes using Poisson methods with gap forecasts made using the more accurate Skellam method can also enhance the forecasting accuracy for the other cause of death. (2) Including other covariates in the models to account for the geospatial dimension in the difference in the death counts of the two populations.