Modeling Best Practice Life Expectancy Using Gumbel Autoregressive Models

Abstract

Best practice life expectancy has recently been modeled using extreme value theory. In this paper we present the Gumbel autoregressive model of order one—Gumbel AR(1)—as an option for modeling best practice life expectancy. This class of model represents a neat and coherent framework for modeling time series extremes. The Gumbel distribution accounts for the extreme nature of best practice life expectancy, while the AR structure accounts for the temporal dependence in the time series. Model diagnostics and simulation results indicate that these models present a viable alternative to Gaussian AR(1) models when dealing with time series of extremes and merit further exploration.

1. Introduction

Best practice life expectancy (BPLE) (Oeppen and Vaupel 2002; Vallin and Meslé 2009) is the maximum life expectancy from among national populations during a given year, at a particular age. Since BPLE is just an annual maximum, Medford (2017) proposed to model it using extreme value theory (EVT). Since then, there have been interesting contributions surrounding the modeling and potential applications of BPLE from Liu and Li (2019) and Li and Liu (2020). Liu and Li (2019) used BPLE to approximate and extrapolate lower and upper bounds in life expectancy. Medford (2017) assumed that the dependencies between successive annual BPLEs were captured in time effects and fitted a time-varying generalized extreme value distribution (GEV), in effect a separate GEV for each year. However, Li and Liu (2020) explored a more sophisticated approach, using Archimedean copulas to account for the dependence based on the work of Wüthrich (2004). For detailed discussions of the issues around using EVT to fit BPLEs in practice, readers may refer to Medford (2017), which includes detailed exposition for the assumptions.

Gaussian models are the basic models for linear time series (Hamilton 1994; Harvey 1993). The underlying error process driving the series is assumed to have no autocorrelation, with a mean of zero and a constant variance. It is also not uncommon to additionally assume that the error process is itself Gaussian. Hence both the assumed innovations (error) and marginal distributions are often assumed to be normally distributed. However, in situations where data record extreme events, Gaussian time series are unsuitable.

In the context of Autoregressive integrated moving average (ARIMA) models, assuming GEV errors may be a reasonable assumption when modeling extremes: maxima or minima. Hughes et al. (2007) used the GEV to model innovations of time series of extreme Antarctic temperatures. However, there was no discourse around the stationary marginal distribution or any attempt to set it out explicitly. Toulemonde et al. (2010) used autoregressive models with a stationary marginal Gumbel extreme value distribution to model atmospheric methane and nitrous oxide levels. More recently, Balakrishna and Shiji (2014) applied a similar model to daily maxima of the Bombay Stock Exchange and the Standard and Poor’s Index.

The models employed by Balakrishna and Shiji (2014) and Toulemonde et al. (2010) were obtained by mixing extreme value distributions over a positive stable distribution: adding a Gumbel distributed variable to a positive $\alpha$-stable variable results in a Gumbel distributed variable too. This was previously highlighted by Crowder (1989); Hougaard (1986); Watson and Smith (1985) in the context of survival analysis. Tawn (1990) presented such models in the modeling of multivariate extremes, while Fougeres et al. (2009) did so in a mixture context. Other interesting applications are in Crowder (1998). Gumbel models were previously fit to BPLE (Liu and Li 2019; Medford 2017) but not with the modeling framework that we present here. The Gumbel distribution accounts for the extreme nature of BPLE, while the autoregressive (AR) structure accounts for the temporal dependence in the time series of BPLE. Theoretically, extreme value innovations are more suitable for time series of maxima than Gaussian innovations. We remain in the familiar ARIMA model framework but with some added complexity to reflect the extreme value marginal distributions from the innovations.

The aim of this paper is to present Gumbel autoregressive model of order one (Gumbel AR(1)) as an option for modeling BPLE, since they are not particularly well known, and much less widely used. They offer an alternative to modeling short-term dependencies among maxima coming from light-tailed distributions. We do not attempt to introduce new theory or methodology, as the focus is on demonstrating how this family of models can be straightforwardly applied and their advantages over traditional Gaussian AR models. We believe that these models merit further study and development and could be used in a range of contexts. This paper is structured as follows. Section 2 presents some background on the Gumbel distribution and stable distributions. Section 3 provides detail on the statistical basis of the Gumbel AR(1) model. Section 4 steps through the application of the model and Section 5 concludes.

2. Some Preliminaries

2.1. The Gumbel Distribution

The GEV distribution function is given by

$$F\left(x\right)=\exp \left\{-{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]}^{\frac{-1}{\xi }}\right\},$$

(1)

defined on $\{x:1+\xi (z-\mu )/\sigma >0\}$. The model is described by three parameters: $\mu (-\infty <\mu <\infty )$, $\sigma (\sigma >0)$ and $\xi (-\infty <\xi <\infty )$ referred to as the location, scale and shape parameters, respectively.

The shape parameter, $\xi$, determines the heaviness of the right tail, and this leads to three types of distributions. When $\xi <0$, the distribution has finite support and is short-tailed, leading to the Weibull distribution. When $\xi >0$, there is polynomial tail decay, leading to heavy tails, and the GEV is of the Fréchet type. The case where $\xi =0$ is taken to be the limit of Equation (1), as $\xi \to 0$ and there is exponential tail decay leading to light tails, and the GEV is of the Gumbel type with distribution function

$$F\left(x\right)=\exp \left\{-\exp \left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\right\}.$$

(2)

The characteristic function (cf) of the Gumbel distribution is given by

$${\varphi }_x\left(t\right)=E\left(\exp \left(itX\right)\right)=\Gamma (1-i\sigma t)e^{i\mu t}$$

(3)

where $\Gamma (.)$ is the gamma integral and $i=\sqrt{-1}$ is the complex number.

2.2. Positive $\alpha$-Stable Distributions

A random variable S is said to be $\alpha$-stable if for all non-negative real numbers $c_1$ and $c_2$, there exist positive real numbers a and b such that $c_1{S}_1$ + $c_2{S}_2$ is equal in distribution to $aS+b$ where $S_1$ and $S_2$ are independent, identically distributed (iid) copies of S. There are three cases where one can write down closed form expressions for the density and verify directly that they are stable. These are the normal, Cauchy and Lévy distributions. For other distributions, there are no closed form solutions.

The Gumbel AR model ultimately arises from the additive relationship between the Gumbel distribution and positive $\alpha$-stable variables.

Proposition 1.

Let S be an α-stable variable defined by its Laplace transform:

$$\mathbb{E}\left(exp(-uS)\right)=exp(-u^{\alpha }),$$

(4)

for all $u\ge 0$ and $\alpha \in (0,1)$. Let X be Gumbel distributed with parameters μ and σ and be independent of S. Then the sum $X+\sigma logS$ is also Gumbel distributed with parameters μ and $\sigma /\mu$.

Proof. 

The moment generating function of X is $e^{\mu t}\Gamma (1-\sigma t)$. Let $\sigma \log S=Y$. Y is obviously exponential-stable and has moment generating function (mgf) $\Gamma (1-t\sigma /\alpha )/\Gamma (1-\sigma t)$. Therefore, the mgf of $X+Y$ is

$$\mathbb{E}{e}^{tX}\mathbb{E}{e}^{tY}=e^{\mu t}\Gamma (1-\sigma t){\displaystyle \frac{\Gamma (1-t\sigma /\alpha )}{\Gamma (1-\sigma t)}}=e^{\mu t}\Gamma (1-t\sigma /\alpha )$$

which is the mgf of a Gumbel distribution with parameters $\mu$ and $\sigma /\alpha )$. □

3. Gumbel AR(1) Model

We begin with the framework of the simple stationary AR(1) model. Let a sequence of independent, identically distributed random variables be $\left\{{ϵ}_t\right\}$ and define a stationary AR(1) sequence $\left\{X_t\right\}$ by

$$X_t=\alpha X_{t-1}+{ϵ}_t,0\le \alpha \le 1,=1,2,\dots$$

(5)

with $X_0$ and ${ϵ}_1$ being independent. Suppose that ${ϵ}_t=\alpha \log S_t$ where S is a positive stable random variable with the Laplace transform given in (4). Our goal is to have the marginal value of $\left\{X_t\right\}$ be an extreme value distribution of Gumbel type, as in Equation (2).

If F is the marginal distribution of $\left\{X_t\right\}$ in (6), a proper distribution for the innovation ${ϵ}_t$ exists if and only if F is self-decomposable. The random variable X is self-decomposable if the ratio ${\varphi }_X\left(t\right)/{\varphi }_X\left(\alpha t\right)$ is also well-defined for all $\alpha$ ($0<\alpha <1$). For the Gumbel AR(1) model this ratio is given by

$${\varphi }_{ϵ}\left(t\right)={\displaystyle \frac{{\varphi }_X\left(t\right)}{{\varphi }_X\left(\alpha t\right)}}{e}^{it(1-\alpha )\mu )}{\displaystyle \frac{\Gamma (1-i\sigma t)}{\Gamma (1-i\sigma \alpha t)}}.$$

(6)

Based on results from Brockwell and Brown (1978); Nolan (2020); Zolotarev (1986), it can be shown that ${ϵ}_t$ has distribution $(1-\alpha )\mu +\sigma Z$ where Z has the distribution −log$\left(S^{-\alpha }\right)$. The mean and variance of the innovation random variable are given by

$$E\left({ϵ}_t\right)=(1-\alpha )(\mu +\sigma \gamma )\mathrm{and}Var\left({ϵ}_t\right)={\displaystyle \frac{(1-{\alpha }^2)\left({\pi }^2{\sigma }^2\right)}{6}}$$

(7)

where $\gamma \approx 0.5772$ is Euler’s constant. It follows then that

$$E\left(X_t\right)=\mu +\sigma \gamma \mathrm{and}Var\left(X_t\right)={\displaystyle \frac{{\pi }^2{\sigma }^2}{6}}$$

(8)

Various estimation methods are possible (Balakrishna and Shiji 2014), but we adopt a simple method of moments approach, similar in spirit to Toulemonde et al. (2010) in order to obtain estimates of the three unknown parameters of our Gumbel AR(1) model: $\mu$, $\sigma$ and $\alpha$. Therefore, we derive the following equations which we solve for $\mu$ and $\sigma$. Since $X_t$ are strictly stationary (Balakrishna and Shiji 2014; Toulemonde et al. 2010), we can write

$$\overline{X}=\mu +\sigma \gamma \mathrm{and}{s}^2={\displaystyle \frac{{\pi }^2{\sigma }^2}{6}}$$

(9)

leading to the method of moment estimators:

$$\widehat{\mu }=\overline{X}-\gamma {\displaystyle \frac{\sqrt{6}}{\pi }}s\mathrm{and}\widehat{\sigma }={\displaystyle \frac{\sqrt{6}}{\pi }}s$$

where $\overline{X}={\sum }_{t=1}^n{X}_t/n$ and $s^2={\sum }_{t=1}^n(X_t-\overline{X})/n-1$. We estimate $\sigma$ using a Yule–Walker type estimate, since the method of moments cannot generate an estimator for it. Hence, in line with Toulemonde et al. (2010)

$$\widehat{\alpha }={\displaystyle \frac{{\sum }_{t=1}^{n-1}(X_t-\overline{X})(X_{t+1}-\overline{X})}{{\sum }_{t=1}^n{(X_t-\overline{X})}^2}}.$$

Note that $\sqrt{n}(\widehat{\alpha }-\alpha )\stackrel{d}{\to }N(0,1-{\alpha }^2)asn\to \infty$ and allows asymptotic standard errors to be estimated. The standard errors for $\widehat{\mu }$ and $\widehat{\sigma }$ are rather more complicated, but details of them can be found in (Toulemonde et al. 2010).

4. Illustration

The data used in model implementation and testing came from two sources. First, the Human Mortality Database (HMD, 2020) (Human Mortality Database 2020) covers the low-mortality countries that have the best data and the highest life expectancies. It contains life tables for 41 countries plus all the raw data used in constructing those tables. The specific data used were life expectancy at birth for both males and females for the period 1965 to 2017. The second source of data was the United Nations world population prospects (United Nations Population Division 2019). This was used to supplement HMD data where they were not yet available for the most recent years (New Zealand, Taiwan, Canada and Israel) and to obtain the life expectancy for all the other countries and territories not found in the HMD. The choice of fitting period (1965 to 2017) was somewhat arbitrary, but was chosen in an attempt to strike a balance between using the most recent data and trends in BPLE (Liu and Li 2019; Vallin and Meslé 2009), and having sufficient data for reasonable parameter estimation.

To illustrate, we fit the AR(1) model with a Gumbel marginal distribution to time series of male and female best practice life expectancy at birth. The BPLEs which have been extracted from the data are presented in Figure 1.

We modeled the time series using standard ARIMA time series approaches. Since the time series trend upward, we first made it stationary. As the series evolves almost linearly and it was previously argued that the trend in BPLE is deterministic (Medford and Vaupel 2020), we achieved stationarity by subtracting the linear trend found by fitting a simple linear regression model. BPLE increased at about 0.21 years per annum for males and 0.22 years per annum for females in this data window.

We needed to check the order of the fitted series and confirm that it has order one. We did this via the partial autocorrelation function (PACF). That is shown in Figure 2. In particular, the PACF cuts off at lag 1, suggesting that the autoregressive model of order 1 assumption is reasonable. An AR(1) model was also found to be the optimal model based on the AICc (corrected AIC) criterion.

The parameters $\widehat{\alpha }$, $\widehat{\mu }$ and $\widehat{\sigma }$ were estimated using the method outlined in the previous section. These estimates, along with their standard errors, are presented in

Table 1. Estimated parameter values of a fitted Gumbel autoregressive model of order one—Gumbel AR(1) with estimated asymptotic standard errors.

	Male	Female
$\widehat{\alpha }$	0.30 (0.13)	0.49 (0.12)
$\widehat{\mu }$	70.5 (0.05)	75.9 (0.06)
$\widehat{\sigma }$	0.27 (0.04)	0.24 (0.04)

. Furthermore, the confidence interval for $\widehat{\alpha }$ was (0.04, 0.56) for males and (0.24, 0.72) for females. Zero does not fall within these intervals, indicating that $\alpha$ is a significant parameter and that the AR model may be appropriate.

Appendix A presents diagnostics for the fitted model. The diagnostic tests confirm that an AR(1) model fits the time series well and that the Gumbel distribution is an acceptable choice for the marginal distribution. Therefore, we conclude that the Gumbel AR(1) model is adequate for the data.

A Comparison with Gaussian AR(1)

Given the extra effort involved in Gumbel AR(1) modeling, is it more accurate than using a Gaussian AR(1) model? The predictive distribution of $[X_{t+1}\mid X_t=x]$ in Equation (1) is a log positive $\alpha$-stable random variable (Toulemonde et al. 2010). Using the estimated parameters in Table 1, we simulated 10,000 observations from the fitted Gumbel AR(1) model and from a counter-factual Gaussian AR(1) model. We used a visual check of the histograms of the simulated observations from the two models which were then superimposed with the true Gumbel density. These comparisons are presented in Figure 3 and Figure 4. It is evident that the Gaussian distribution, because it is symmetric, fits the tails poorly and is unable to capture the asymmetry (short left tail, long right tail) in the data.

5. Conclusions

In this paper we presented a first-order Gumbel autoregressive model and fit it to a time series of best practice life expectancies. Gumbel AR(1) models can be used to model short term temporal dependence among extreme values which come from distributions that have reasonably light tails. The Gumbel AR(1) is rather limited, however, and it would be more appropriate to have greater flexiblity in the model. Greater flexibility includes, for example, being able to handle data with heavier tails and being able to fit different types of time series. A more general model should be able to fit the other extreme value distributions, the Frechet and Weibull distributions. It should also be able to accommodate the general class of ARIMA time series, not just AR(1) models. This opens up a potentially interesting area for future research. Forecasting using these types of models might also be an avenue that merits further exploration. In our view the usefulness and applicability of these models have not been fully explored and appear to be ripe for further development