Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Modelling the age distribution of longevity leaders

Csaba Kiss Department of Stochastics, Institute of Mathematics, Budapest University of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary. László Németh Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstraße 39, D-10117 Berlin, Germany. Max Planck Institute for Demographic Research, Konrad-Zuse-Str. 1., D-18057 Rostock, Germany. nemeth@wias-berlin.de Bálint Vető Department of Stochastics, Institute of Mathematics, Budapest University of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary. ELKH–BME Stochastics Research Group, Műegyetem rkp. 3., H-1111 Budapest, Hungary.
Abstract

Human longevity leaders with remarkably long lifespan play a crucial role in the advancement of longevity research. In this paper, we propose a stochastic model to describe the evolution of the age of the oldest person in the world by a Markov process, in which we assume that the births of the individuals follow a Poisson process with increasing intensity, lifespans of individuals are independent and can be characterized by a gamma–Gompertz distribution with time-dependent parameters. We utilize a dataset of the world’s oldest person title holders since 1955, and we compute the maximum likelihood estimate for the parameters iteratively by numerical integration. Based on our preliminary estimates, the model provides a good fit to the data and shows that the age of the oldest person alive increases over time in the future. The estimated parameters enable us to describe the distribution of the age of the record holder process at a future time point.

keywords:
human longevity, world’s oldest person, stochastic model, gamma–Gompertz lifespan distribution

Introduction

Exceptionally long human lifespans are one of the cornerstones of demography and mortality research. Studying the group of record holders may reveal not only the underlying mortality mechanism of a population but also potentially shed some light on the future developments of human longevity. As life expectancy increases [1] with deaths shifting to older ages, the distribution of deaths at the oldest-old ages [2, 3] gains more interest of demographers, actuaries and decision makers of numerous disciplines.

The pattern of adult human mortality has already been described [4] but there is still a debate about the exact distribution of deaths at adult ages. More details on the necessity of a correctly specified model for the underlying mortality process and its impact on further research are discussed in [5, 6]. Numerous publications review the ongoing discussion on the existence of a mortality plateau [7, 8, 9, 10, 11, 12], and the levelling-off of adult human death rates at the oldest ages is supported by the findings in [13, 14, 15] while others cast some doubt on this observation [16, 17, 18]. If there is a mortality plateau then the distribution of deaths at the oldest-old ages must be gamma-Gompertz and human lifespan can increase further without any maximum [19]. Further studies discuss the existence of a limit to human lifespan with more focus on the extreme value distribution aspect of the deaths at the oldest-old ages for various populations and different lifetime distributions [20, 21, 22, 23, 24]. These models can be helpful in determining the plausibility of longevity leaders as well.

We contribute to this discussion by proposing a stochastic model to describe the evolution of the age of the world’s oldest person. Based on our estimates the model provides a good fit to the titleholder data since 1955, collected by the Gerontology Research Group [25]. With the model results, it is possible to predict the age of the oldest person in the world in the future. When should we expect to see the next Jeanne Calment, the supercentenarian with the longest human lifespan ever documented? Will her record ever be surpassed? Our results provide a prediction for the age distribution of the record holder in the coming decades to answer these questions.

Results

Our model describes the evolution of the age of the oldest living person under the following assumptions. We assume that the births of individuals follow a Poisson process with time-dependent intensity [26]. The lifespans of individuals in the population are independent and their distribution may depend on the date of birth. Then the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities. As the first main result of this paper, we explicitly compute the distribution of the age of the record holder for any given birth rate parameter and lifespan distribution. The detailed mathematical description and the properties of the general model are described in Section 1.

We apply our general result to the case which approximates the human birth rate over the world and the human lifespan. We specify the intensity function of the Poisson process of births to have an exponential growth in time. The underlying force of mortality is chosen so that it follows an extension of the Gompertz mortality model [4] and the lifespan distribution of individuals is given by the gamma–Gompertz (ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G) distribution with time-dependent parameters. This distribution adequately captures the slowing down of senescence mortality at the oldest old ages. Given the growth parameters of the birth rate, we fit the model parameters to the statistics of the oldest person titleholder data using maximum likelihood method. The optimal parameters of the model fit well to the data. It shows in particular that the age of the oldest person alive increases over time, and it will most likely increase further in the future. We compute the expected value and a confidence interval for the age of the world’s oldest person using the fitted model parameters for each year between 1955 and 2019 shown by the green curves on Figure 1. The detailed discussion of the model specification, likelihood calculations as well as the parameter fitting are given in Section 2. Section 3 contains calculations related to the gamma–Gompertz–Makeham generalization of the gamma–Gompertz distribution.

Refer to caption
Figure 1: Blue curve: age of the oldest person in the world since 1955 with two outliers indicated (Jeanne Calment: JC and Sarah Knauss: SK). Green curves: the model-estimated mean age of the oldest person (solid line) and the standard deviation (distance from the dashed lines) with the estimation based on the full dataset. Red curves: mean age and standard deviation of the oldest person estimated using the data between 1955 and 1988.

Our results enable us to predict the age distribution of the world’s oldest person at future time points. We compute the probability density of the age of the world’s oldest person in different years not only in the past but also in the future. These densities are shown in Figure 2. When comparing the age distribution of the oldest person in the world in different years to the age of Jeanne Calment at her death, we find that on Jan 1st 2060 we can expect that the age of the world’s oldest person will exceed her age with probability around 0.50.50.50.5. This also means that with high probability her age record will already be broken by that time.

Refer to caption
Figure 2: The probability densities of the age of the oldest person in the world in the following years: 1960 (red), 1980 (green), 2000 (blue), 2020 (pink), 2040 (orange), 2060 (magenta), 2080 (yellow) and 2100 (black).

In Figure 1 two extreme outliers with unexpectedly long lifetimes can be observed. Jeanne Calment died at the age of 122.45 years in 1997, and Sarah Knauss died at the age 119.27 in 2000. In our model, the probability of observing an age greater than or equal to their actual age at the time of their death is 0.0002860.0002860.0002860.000286 for Calment and 0.01160.01160.01160.0116 for Knauss. See details in Subsection 2.4.

The fact that Calment and Knauss are outliers among the oldest old in the world became even more evident when we performed a backtesting of our model. We estimated the parameters based on the data on the world’s oldest person between 1955 and 1988 where the ending date is the time when Calment became the world’s oldest person. The model-estimated mean and confidence interval of the world’s oldest person using the full data and the partial dataset (before Calment) are shown in Figure 1 by green and red, respectively. The estimate using the data until 1988 is less reliable after 2000 which is shown by the fact that the observed data is out of the confidence interval in the majority of the time after 2000. When we compare the two confidence intervals, we can conclude that, based on the data before 1988, Calment and Knauss already had extremely high ages at their deaths. Adding the remaining data set, the estimated mean age of the world’s oldest person becomes lower. Hence we can conclude that now we can consider the ages of the two outliers between 1988 and 2000 at their death to be more extreme than based on the information available in 1988.

The other important observable is the reign length of a record holder. The numerical value of the expected reign length with our estimated model parameters is 1.1951.1951.1951.195 in 1955 and it is 1.1881.1881.1881.188 in 2019. The empirical value of the reign length is 1.0081.0081.0081.008 which is not much less than the model-based estimate.

Our approach to studying the age of the oldest old is completely novel because it takes into account jointly the age and the time of birth of individuals. Although the age and the reign length of the world’s oldest person depend in a complex non-linear way on the total lifespan and the time of birth of supercentenarians, we compute explicitly the probability distribution of the age of the oldest person. Hence, the performance of our predictions cannot be directly compared to previous results in the literature because in the usual approach, the oldest person in each cohort is considered separately, and it is not relevant whether this person was ever the oldest in the population, see e.g. the extreme value method in [20]. Our model contributes to the mathematical understanding of the evolution of the oldest individual, which is the extra benefit compared to a prediction using the trend in the data, e.g. in linear regression. In this way, we not only observe but also prove mathematically that the dynamics of the birth process and that of the lifespan distribution which we consider in this paper necessarily imply the increase of the expected age of the world’s oldest person.

1 Mathematical model for the age of the oldest person

In this section, we provide the mathematical definition of a general model for the age of the world’s oldest person, where the births of individuals follow a Poisson process and their lifespans are independent. Under the assumptions of the general model, the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities. In Subsections 1.11.2, we describe the exact distribution of the age of the record holder in this generality for any given birth rate parameter and lifespan distribution using the two-dimensional representation of the age process of the oldest person. In the time-homogeneous case with constant birth rate and identical lifespan distributions the reign length distribution of a record holder is computed in Subsections 1.31.5. We explain the role of the entry age parameter in Subsection 1.6.

1.1 Model description and two-dimensional representation

The model is formally defined as follows. Let λ(t)𝜆𝑡\lambda(t)italic_λ ( italic_t ) be the birth rate parameter which depends on time and let Ftsubscript𝐹𝑡F_{t}italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT and ftsubscript𝑓𝑡f_{t}italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT be a family of cumulative distribution functions and density functions corresponding to non-negative random variables which are also time-dependent. We assume that individuals are born according to a Poisson point process at rate λ(t)𝜆𝑡\lambda(t)italic_λ ( italic_t ) and that the lifespan of an individual born at time t𝑡titalic_t is given by Ftsubscript𝐹𝑡F_{t}italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT so that lifespans are independent for different individuals.

Let Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT denote the age of the oldest person in the population at time t𝑡titalic_t. The process (Yt:t):subscript𝑌𝑡𝑡(Y_{t}:t\in\mathbb{R})( italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT : italic_t ∈ blackboard_R ) is Markovian. The Markov property holds because at any time t𝑡titalic_t the history of the process (Ys:st):subscript𝑌𝑠𝑠𝑡(Y_{s}:s\leq t)( italic_Y start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT : italic_s ≤ italic_t ) provides information about the lifetime of individuals born before the current record holder while any transition of (Ys:st):subscript𝑌𝑠𝑠𝑡(Y_{s}:s\geq t)( italic_Y start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT : italic_s ≥ italic_t ) depends only on the lifetime of the current record holder and of those born after them.

The evolution of the Markov process Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT is the following. It has a deterministic linear growth with slope 1111 due to the ageing of the current record holder. This happens until the death of the record holder. Additionally, given that Yt=limstYs=ysubscript𝑌limit-from𝑡subscript𝑠𝑡subscript𝑌𝑠𝑦Y_{t-}=\lim_{s\uparrow t}Y_{s}=yitalic_Y start_POSTSUBSCRIPT italic_t - end_POSTSUBSCRIPT = roman_lim start_POSTSUBSCRIPT italic_s ↑ italic_t end_POSTSUBSCRIPT italic_Y start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = italic_y for some t𝑡titalic_t with y>0𝑦0y>0italic_y > 0, the process has a downward jump at time t𝑡titalic_t at rate ft(y)/(1Ft(y))subscript𝑓𝑡𝑦1subscript𝐹𝑡𝑦f_{t}(y)/(1-F_{t}(y))italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_y ) / ( 1 - italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_y ) ) which is the hazard rate of the distribution Ftsubscript𝐹𝑡F_{t}italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT at y𝑦yitalic_y. This corresponds to the possibility that the record holder dies at time t𝑡titalic_t which happens at rate ft(y)/(1Ft(y))subscript𝑓𝑡𝑦1subscript𝐹𝑡𝑦f_{t}(y)/(1-F_{t}(y))italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_y ) / ( 1 - italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_y ) ). The conditional distribution of the jump is given by

𝐏(Yt<x|Yt=y,Yt<y)=exp(xyλ(tu)(1Ftu(u))du)𝐏formulae-sequencesubscript𝑌𝑡bra𝑥subscript𝑌limit-from𝑡𝑦subscript𝑌𝑡𝑦superscriptsubscript𝑥𝑦𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢\mathbf{P}(Y_{t}<x\,|\,Y_{t-}=y,Y_{t}<y)=\exp\left(-\int_{x}^{y}\lambda(t-u)(1% -F_{t-u}(u))\,\mathrm{d}u\right)bold_P ( italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x | italic_Y start_POSTSUBSCRIPT italic_t - end_POSTSUBSCRIPT = italic_y , italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_y ) = roman_exp ( - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) (1)

for all x>0𝑥0x>0italic_x > 0. The jump distribution in (1) has an absolutely continuous part supported on [0,y]0𝑦[0,y][ 0 , italic_y ] with density

jy,t(x)=exp(xyλ(tu)(1Ftu(u))du)λ(tx)(1Ftx(x))subscript𝑗𝑦𝑡𝑥superscriptsubscript𝑥𝑦𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢𝜆𝑡𝑥1subscript𝐹𝑡𝑥𝑥j_{y,t}(x)=\exp\left(-\int_{x}^{y}\lambda(t-u)(1-F_{t-u}(u))\,\mathrm{d}u% \right)\lambda(t-x)\left(1-F_{t-x}(x)\right)italic_j start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT ( italic_x ) = roman_exp ( - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) italic_λ ( italic_t - italic_x ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_x end_POSTSUBSCRIPT ( italic_x ) ) (2)

and a point mass at 00 with probability

ay,t=exp(0yλ(tu)(1Ftu(u))du).subscript𝑎𝑦𝑡superscriptsubscript0𝑦𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢a_{y,t}=\exp\left(-\int_{0}^{y}\lambda(t-u)(1-F_{t-u}(u))\,\mathrm{d}u\right).italic_a start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT = roman_exp ( - ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) . (3)

As we shall see in the relevant parameter regime the probability ay,tsubscript𝑎𝑦𝑡a_{y,t}italic_a start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT of the point mass at 00 is negligible. The transition formula in (1) can be proven using the description below.

We introduce a two-dimensional representation of the process Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT as follows. Let Λ={(ti,xi):iI}Λconditional-setsubscript𝑡𝑖subscript𝑥𝑖𝑖𝐼\Lambda=\{(t_{i},x_{i}):i\in I\}roman_Λ = { ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) : italic_i ∈ italic_I } be a marked Poisson process in ×+subscript\mathbb{R}\times\mathbb{R}_{+}blackboard_R × blackboard_R start_POSTSUBSCRIPT + end_POSTSUBSCRIPT where {ti:iI}conditional-setsubscript𝑡𝑖𝑖𝐼\{t_{i}:i\in I\}{ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT : italic_i ∈ italic_I } forms a Poisson point process on \mathbb{R}blackboard_R with intensity λ(t)𝜆𝑡\lambda(t)italic_λ ( italic_t ) and xi0subscript𝑥𝑖0x_{i}\geq 0italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ≥ 0 is sampled independently for each iI𝑖𝐼i\in Iitalic_i ∈ italic_I according to the distribution Ftisubscript𝐹subscript𝑡𝑖F_{t_{i}}italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT. The point (ti,xi)subscript𝑡𝑖subscript𝑥𝑖(t_{i},x_{i})( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) represents an individual born at time tisubscript𝑡𝑖t_{i}italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT with lifespan xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT for all iI𝑖𝐼i\in Iitalic_i ∈ italic_I, that is, the individual i𝑖iitalic_i is alive in the time interval [ti,ti+xi)subscript𝑡𝑖subscript𝑡𝑖subscript𝑥𝑖[t_{i},t_{i}+x_{i})[ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) and their age at time t𝑡titalic_t is tti𝑡subscript𝑡𝑖t-t_{i}italic_t - italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT if t[ti,ti+xi)𝑡subscript𝑡𝑖subscript𝑡𝑖subscript𝑥𝑖t\in[t_{i},t_{i}+x_{i})italic_t ∈ [ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ). Hence the marked Poisson process ΛΛ\Lambdaroman_Λ contains all relevant information about the age statistics of the population at any time. In particular the age of the oldest person Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT can be expressed in terms of ΛΛ\Lambdaroman_Λ as

Yt=max{(tti)𝟙{t[ti,ti+xi)}:iI}.subscript𝑌𝑡:𝑡subscript𝑡𝑖subscript1𝑡subscript𝑡𝑖subscript𝑡𝑖subscript𝑥𝑖𝑖𝐼Y_{t}=\max\left\{(t-t_{i})\mathbbm{1}_{\{t\in[t_{i},t_{i}+x_{i})\}}:i\in I% \right\}.italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = roman_max { ( italic_t - italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) blackboard_1 start_POSTSUBSCRIPT { italic_t ∈ [ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) } end_POSTSUBSCRIPT : italic_i ∈ italic_I } . (4)

where the indicator 𝟙{t[ti,ti+xi)}subscript1𝑡subscript𝑡𝑖subscript𝑡𝑖subscript𝑥𝑖\mathbbm{1}_{\{t\in[t_{i},t_{i}+x_{i})\}}blackboard_1 start_POSTSUBSCRIPT { italic_t ∈ [ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) } end_POSTSUBSCRIPT is 1111 exactly if the i𝑖iitalic_ith person is alive at time t𝑡titalic_t.

The transition distribution formula (1) can be seen using the two-dimensional representation as follows. Given that the current record holder dies at time t𝑡titalic_t at age y𝑦yitalic_y the event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x } means that nobody with age between x𝑥xitalic_x and y𝑦yitalic_y can be alive at time t𝑡titalic_t. This event can be equivalently characterized in terms of the Poisson process of birth at rate λ()𝜆\lambda(\cdot)italic_λ ( ⋅ ) thinned by the probability that the person is still alive at time t𝑡titalic_t. Indeed the event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x } can be expressed as a Poisson process of intensity at time tu𝑡𝑢t-uitalic_t - italic_u given by λ(tu)(1Ftu(u))𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢\lambda(t-u)(1-F_{t-u}(u))italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) for u[x,y]𝑢𝑥𝑦u\in[x,y]italic_u ∈ [ italic_x , italic_y ] not having any point in the time interval [ty,tx]𝑡𝑦𝑡𝑥[t-y,t-x][ italic_t - italic_y , italic_t - italic_x ]. This probability appears exactly on the right-hand side of (1). In other words, for any u[x,y]𝑢𝑥𝑦u\in[x,y]italic_u ∈ [ italic_x , italic_y ] people are born at time tu𝑡𝑢t-uitalic_t - italic_u at rate λ(tu)𝜆𝑡𝑢\lambda(t-u)italic_λ ( italic_t - italic_u ). On the other hand, the probability for a person born at time tu𝑡𝑢t-uitalic_t - italic_u to be alive at time t𝑡titalic_t (that is, at age u𝑢uitalic_u) is 1Ftu(u)1subscript𝐹𝑡𝑢𝑢1-F_{t-u}(u)1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ). See Figure 3 for illustration.

Refer to caption
Figure 3: Two-dimensional representation of the event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x } conditionally given {Yt=y,Yt<y}formulae-sequencesubscript𝑌limit-from𝑡𝑦subscript𝑌𝑡𝑦\{Y_{t-}=y,Y_{t}<y\}{ italic_Y start_POSTSUBSCRIPT italic_t - end_POSTSUBSCRIPT = italic_y , italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_y }: the marked Poisson process does not have any point in the grey region, see (1).

1.2 Exact distribution of the oldest person’s age process

We assume that all birth events are already sampled on (,t]𝑡(-\infty,t]( - ∞ , italic_t ] together with the corresponding lifespans. Then the distribution of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT can be computed explicitly for all t𝑡t\in\mathbb{R}italic_t ∈ blackboard_R using the two-dimensional representation. For all t𝑡t\in\mathbb{R}italic_t ∈ blackboard_R the density

ht(x)=exp(xλ(tu)(1Ftu(u))du)λ(tx)(1Ftx(x))subscript𝑡𝑥superscriptsubscript𝑥𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢𝜆𝑡𝑥1subscript𝐹𝑡𝑥𝑥h_{t}(x)=\exp\left(-\int_{x}^{\infty}\lambda(t-u)(1-F_{t-u}(u))\,\mathrm{d}u% \right)\lambda(t-x)(1-F_{t-x}(x))italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_x ) = roman_exp ( - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) italic_λ ( italic_t - italic_x ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_x end_POSTSUBSCRIPT ( italic_x ) ) (5)

for all x>0𝑥0x>0italic_x > 0 and the point mass at 00

mt=exp(0λ(tu)(1Ftu(u))du)subscript𝑚𝑡superscriptsubscript0𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢m_{t}=\exp\left(-\int_{0}^{\infty}\lambda(t-u)(1-F_{t-u}(u))\,\mathrm{d}u\right)italic_m start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = roman_exp ( - ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) (6)

characterize the distribution of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT which can be seen as follows. We mention that the point mass at 00 is negligible in the application.

Similarly to the proof of the transition formula in (1) the event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x } for any x>0𝑥0x>0italic_x > 0 is the same as the event that nobody with age at least x𝑥xitalic_x is alive at time t𝑡titalic_t. We express this event in terms of the Poisson process of birth at rate λ()𝜆\lambda(\cdot)italic_λ ( ⋅ ) thinned by the probability that the person is still alive at time t𝑡titalic_t. The event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x } means that a Poisson process of intensity at time tu𝑡𝑢t-uitalic_t - italic_u given by λ(tu)(1Ftu(u))𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢\lambda(t-u)(1-F_{t-u}(u))italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) for ux𝑢𝑥u\geq xitalic_u ≥ italic_x does not have any point in (,tx]𝑡𝑥(-\infty,t-x]( - ∞ , italic_t - italic_x ] yielding

𝐏(Yt<x)=exp(xλ(tu)(1Ftu(u))du).𝐏subscript𝑌𝑡𝑥superscriptsubscript𝑥𝜆𝑡𝑢1subscript𝐹𝑡𝑢𝑢differential-d𝑢\mathbf{P}(Y_{t}<x)=\exp\left(-\int_{x}^{\infty}\lambda(t-u)(1-F_{t-u}(u))\,% \mathrm{d}u\right).bold_P ( italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x ) = roman_exp ( - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u ) . (7)

In other words for any ux𝑢𝑥u\geq xitalic_u ≥ italic_x individuals are born at time tu𝑡𝑢t-uitalic_t - italic_u at rate λ(tu)𝜆𝑡𝑢\lambda(t-u)italic_λ ( italic_t - italic_u ). A person born at time tu𝑡𝑢t-uitalic_t - italic_u is alive at time t𝑡titalic_t at age u𝑢uitalic_u with probability 1Ftu(u)1subscript𝐹𝑡𝑢𝑢1-F_{t-u}(u)1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT ( italic_u ). (5)–(6) follow by differentiation in (7) and by taking the x0𝑥0x\to 0italic_x → 0 limit. See Figure 4 for illustration.

Refer to caption
Figure 4: Two-dimensional representation of the event {Yt<x}subscript𝑌𝑡𝑥\{Y_{t}<x\}{ italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT < italic_x }: the marked Poisson process does not have any point in the grey region, see (7).

1.3 Homogeneous model

The exact computation of the reign length distribution (see Subsection 1.5) can only be performed in a special case of our general model described in Subsection 1.1. We introduce this special case as the homogeneous model where individuals are born at the times of a Poisson process of constant rate λ=λ(t)>0𝜆𝜆𝑡0\lambda=\lambda(t)>0italic_λ = italic_λ ( italic_t ) > 0 for all t𝑡titalic_t. The lifespan of individuals are independent and identically distributed with a fixed density f=ft𝑓subscript𝑓𝑡f=f_{t}italic_f = italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT and cumulative distribution function F=Ft𝐹subscript𝐹𝑡F=F_{t}italic_F = italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT for all t𝑡titalic_t which does not depend on time.

In the homogeneous model the jump distribution given in (2)–(3) simplifies to

jy(x)subscript𝑗𝑦𝑥\displaystyle j_{y}(x)italic_j start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_x ) =jy,t(x)=exp(λxy(1F(u))du)λ(1F(x)),absentsubscript𝑗𝑦𝑡𝑥𝜆superscriptsubscript𝑥𝑦1𝐹𝑢differential-d𝑢𝜆1𝐹𝑥\displaystyle=j_{y,t}(x)=\exp\left(-\lambda\int_{x}^{y}(1-F(u))\,\mathrm{d}u% \right)\lambda(1-F(x)),= italic_j start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT ( italic_x ) = roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) italic_λ ( 1 - italic_F ( italic_x ) ) , (8)
aysubscript𝑎𝑦\displaystyle a_{y}italic_a start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT =ay,t=exp(λ0y(1F(u))du).absentsubscript𝑎𝑦𝑡𝜆superscriptsubscript0𝑦1𝐹𝑢differential-d𝑢\displaystyle=a_{y,t}=\exp\left(-\lambda\int_{0}^{y}(1-F(u))\,\mathrm{d}u% \right).= italic_a start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT = roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) .

The distribution of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT does not depend on time in this case hence it is a stationary distribution as well. The formulas for the density of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT and point mass at 00 reduce in the homogeneous case to

h(x)𝑥\displaystyle h(x)italic_h ( italic_x ) =exp(λx(1F(u))du)λ(1F(x)),absent𝜆superscriptsubscript𝑥1𝐹𝑢differential-d𝑢𝜆1𝐹𝑥\displaystyle=\exp\left(-\lambda\int_{x}^{\infty}(1-F(u))\,\mathrm{d}u\right)% \lambda(1-F(x)),= roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) italic_λ ( 1 - italic_F ( italic_x ) ) , (9)
m𝑚\displaystyle mitalic_m =exp(λ0(1F(u))du)absent𝜆superscriptsubscript01𝐹𝑢differential-d𝑢\displaystyle=\exp\left(-\lambda\int_{0}^{\infty}(1-F(u))\,\mathrm{d}u\right)= roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u )

where the integral 0(1F(u))dusuperscriptsubscript01𝐹𝑢differential-d𝑢\int_{0}^{\infty}(1-F(u))\,\mathrm{d}u∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u is equal to the expected lifespan.

The equilibrium condition for the homogeneous density hhitalic_h can be written as

h(x)+h(x)f(x)1F(x)xh(y)f(y)1F(y)jy(x)dy=0.superscript𝑥𝑥𝑓𝑥1𝐹𝑥superscriptsubscript𝑥𝑦𝑓𝑦1𝐹𝑦subscript𝑗𝑦𝑥differential-d𝑦0h^{\prime}(x)+h(x)\frac{f(x)}{1-F(x)}-\int_{x}^{\infty}h(y)\frac{f(y)}{1-F(y)}% j_{y}(x)\,\mathrm{d}y=0.italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x ) + italic_h ( italic_x ) divide start_ARG italic_f ( italic_x ) end_ARG start_ARG 1 - italic_F ( italic_x ) end_ARG - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) divide start_ARG italic_f ( italic_y ) end_ARG start_ARG 1 - italic_F ( italic_y ) end_ARG italic_j start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_x ) roman_d italic_y = 0 . (10)

After differentiation and using the fact that

ddxjy(x)=jy(x)(λ(1F(x))f(x)1F(x))dd𝑥subscript𝑗𝑦𝑥subscript𝑗𝑦𝑥𝜆1𝐹𝑥𝑓𝑥1𝐹𝑥\frac{\mathrm{d}}{\mathrm{d}x}j_{y}(x)=j_{y}(x)\left(\lambda(1-F(x))-\frac{f(x% )}{1-F(x)}\right)divide start_ARG roman_d end_ARG start_ARG roman_d italic_x end_ARG italic_j start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_x ) = italic_j start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_x ) ( italic_λ ( 1 - italic_F ( italic_x ) ) - divide start_ARG italic_f ( italic_x ) end_ARG start_ARG 1 - italic_F ( italic_x ) end_ARG ) (11)

one can derive from (10) the second order differential equation

h′′(x)+(2f(x)1F(x)λ(1F(x)))h(x)+(2f(x)2(1F(x))2+f(x)1F(x))h(x)=0.superscript′′𝑥2𝑓𝑥1𝐹𝑥𝜆1𝐹𝑥superscript𝑥2𝑓superscript𝑥2superscript1𝐹𝑥2superscript𝑓𝑥1𝐹𝑥𝑥0h^{\prime\prime}(x)+\left(\frac{2f(x)}{1-F(x)}-\lambda(1-F(x))\right)h^{\prime% }(x)+\left(\frac{2f(x)^{2}}{(1-F(x))^{2}}+\frac{f^{\prime}(x)}{1-F(x)}\right)h% (x)=0.italic_h start_POSTSUPERSCRIPT ′ ′ end_POSTSUPERSCRIPT ( italic_x ) + ( divide start_ARG 2 italic_f ( italic_x ) end_ARG start_ARG 1 - italic_F ( italic_x ) end_ARG - italic_λ ( 1 - italic_F ( italic_x ) ) ) italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x ) + ( divide start_ARG 2 italic_f ( italic_x ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 - italic_F ( italic_x ) ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG italic_f start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x ) end_ARG start_ARG 1 - italic_F ( italic_x ) end_ARG ) italic_h ( italic_x ) = 0 . (12)

The point mass m𝑚mitalic_m at 00 satisfies

mλ=0h(y)f(y)1F(y)jy(0)dy.𝑚𝜆superscriptsubscript0𝑦𝑓𝑦1𝐹𝑦subscript𝑗𝑦0differential-d𝑦\ m\lambda=\int_{0}^{\infty}h(y)\frac{f(y)}{1-F(y)}j_{y}(0)\,\mathrm{d}y.italic_m italic_λ = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) divide start_ARG italic_f ( italic_y ) end_ARG start_ARG 1 - italic_F ( italic_y ) end_ARG italic_j start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( 0 ) roman_d italic_y . (13)

1.4 The peaks process

In the homogeneous model the sequence of peaks in Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT forms a discrete time Markov chain. By peak we mean a local maximum of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT with value being equal to the lifespan of the last record holder. Each time the oldest person dies the process Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT has a peak with a downward jump following it. Let Znsubscript𝑍𝑛Z_{n}italic_Z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT denote the age of record holders at which they die which are the values of the peaks of the process Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT. The sequence Znsubscript𝑍𝑛Z_{n}italic_Z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT forms a discrete time Markov chain. The Markov property follows by the fact that ages at death of previous record holders only give information on people born before the current record holder but transitions depend on the lifespan of the current record holder and that of people born after them.

The stationary density of Znsubscript𝑍𝑛Z_{n}italic_Z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT is given by

z(x)=f(x)exp(λx(1F(u))du)0f(y)exp(λy(1F(u))du)dy.𝑧𝑥𝑓𝑥𝜆superscriptsubscript𝑥1𝐹𝑢differential-d𝑢superscriptsubscript0𝑓𝑦𝜆superscriptsubscript𝑦1𝐹𝑢differential-d𝑢differential-d𝑦z(x)=\frac{f(x)\exp\left(-\lambda\int_{x}^{\infty}(1-F(u))\,\mathrm{d}u\right)% }{\int_{0}^{\infty}f(y)\exp\left(-\lambda\int_{y}^{\infty}(1-F(u))\,\mathrm{d}% u\right)\mathrm{d}y}.italic_z ( italic_x ) = divide start_ARG italic_f ( italic_x ) roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) end_ARG start_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_f ( italic_y ) roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) roman_d italic_y end_ARG . (14)

The formula can be seen as follows. To have a record holder who dies at age x𝑥xitalic_x there has to be a person who has lifespan x𝑥xitalic_x which gives the factor f(x)𝑓𝑥f(x)italic_f ( italic_x ) in the numerator on the right-hand side of (14). The exponential factor is by the two-dimensional representation equal to the probability that no people born before the record holder who just died can be alive at the time the record holder dies. The denominator on the right-hand side of (14) makes z(x)𝑧𝑥z(x)italic_z ( italic_x ) a probability density function.

The density of Znsubscript𝑍𝑛Z_{n}italic_Z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT can also be characterized by the following description. It satisfies the integral equation

z(x)=0z(w)0min(x,y)jw(y)f(x)1F(y)dydw+0z(w)awf(x)dw𝑧𝑥superscriptsubscript0𝑧𝑤superscriptsubscript0𝑥𝑦subscript𝑗𝑤𝑦𝑓𝑥1𝐹𝑦differential-d𝑦differential-d𝑤superscriptsubscript0𝑧𝑤subscript𝑎𝑤𝑓𝑥differential-d𝑤z(x)=\int_{0}^{\infty}z(w)\int_{0}^{\min(x,y)}j_{w}(y)\frac{f(x)}{1-F(y)}\,% \mathrm{d}y\,\mathrm{d}w+\int_{0}^{\infty}z(w)a_{w}f(x)\,\mathrm{d}witalic_z ( italic_x ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_z ( italic_w ) ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_min ( italic_x , italic_y ) end_POSTSUPERSCRIPT italic_j start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ( italic_y ) divide start_ARG italic_f ( italic_x ) end_ARG start_ARG 1 - italic_F ( italic_y ) end_ARG roman_d italic_y roman_d italic_w + ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_z ( italic_w ) italic_a start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT italic_f ( italic_x ) roman_d italic_w (15)

which comes from the possible transitions of the peak process as follows. If the previous record holder had a total lifetime w[0,)𝑤0w\in[0,\infty)italic_w ∈ [ 0 , ∞ ) then at the death the process Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT jumps down to some value y𝑦yitalic_y at rate jw(y)subscript𝑗𝑤𝑦j_{w}(y)italic_j start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ( italic_y ) or to 00 with probability awsubscript𝑎𝑤a_{w}italic_a start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT. The density of the age at which a person dies who becomes a record holder at age y𝑦yitalic_y is f(x)/(1F(y))𝑓𝑥1𝐹𝑦f(x)/(1-F(y))italic_f ( italic_x ) / ( 1 - italic_F ( italic_y ) ). From the integral equation in (15) one can derive the second order differential equation for the function g(x)=z(s)/f(x)𝑔𝑥𝑧𝑠𝑓𝑥g(x)=z(s)/f(x)italic_g ( italic_x ) = italic_z ( italic_s ) / italic_f ( italic_x ) given by

g′′(x)λ(1F(x))g(x)+λf(x)g(x)=0superscript𝑔′′𝑥𝜆1𝐹𝑥superscript𝑔𝑥𝜆𝑓𝑥𝑔𝑥0g^{\prime\prime}(x)-\lambda(1-F(x))g^{\prime}(x)+\lambda f(x)g(x)=0italic_g start_POSTSUPERSCRIPT ′ ′ end_POSTSUPERSCRIPT ( italic_x ) - italic_λ ( 1 - italic_F ( italic_x ) ) italic_g start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x ) + italic_λ italic_f ( italic_x ) italic_g ( italic_x ) = 0 (16)

which is satisfied by g(x)=cexp(λx(1F(u))du)𝑔𝑥𝑐𝜆superscriptsubscript𝑥1𝐹𝑢differential-d𝑢g(x)=c\exp\left(-\lambda\int_{x}^{\infty}(1-F(u))\,\mathrm{d}u\right)italic_g ( italic_x ) = italic_c roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) in accordance with (14).

1.5 Reign length distribution

In the homogeneous model, let Wnsubscript𝑊𝑛W_{n}italic_W start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT denote the reign length of the n𝑛nitalic_nth record holder, that is, the time length for which this person is the oldest person of the population. The density of the random reign length is given by

r(w)=0h(y)f(y+w)dy+m0wf(z)λeλ(wz)f(z)dz0h(y)(1F(y))dy+m.𝑟𝑤superscriptsubscript0𝑦𝑓𝑦𝑤differential-d𝑦𝑚superscriptsubscript0𝑤𝑓𝑧𝜆superscript𝑒𝜆𝑤𝑧𝑓𝑧differential-d𝑧superscriptsubscript0𝑦1𝐹𝑦differential-d𝑦𝑚r(w)=\frac{\int_{0}^{\infty}h(y)f(y+w)\,\mathrm{d}y+m\int_{0}^{w}f(z)\lambda e% ^{-\lambda(w-z)}f(z)\,\mathrm{d}z}{\int_{0}^{\infty}h(y)(1-F(y))\,\mathrm{d}y+% m}.italic_r ( italic_w ) = divide start_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) italic_f ( italic_y + italic_w ) roman_d italic_y + italic_m ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_w end_POSTSUPERSCRIPT italic_f ( italic_z ) italic_λ italic_e start_POSTSUPERSCRIPT - italic_λ ( italic_w - italic_z ) end_POSTSUPERSCRIPT italic_f ( italic_z ) roman_d italic_z end_ARG start_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) ( 1 - italic_F ( italic_y ) ) roman_d italic_y + italic_m end_ARG . (17)

The density formula in (17) can be derived based on the stationary density of the peaks process given by (14) as follows.

It holds for the density of the reign length that

r(w)=0z(x)(0xjx(y)f(y+w)1F(y)dy+ax0wf(z)λeλ(wz)dz)dx𝑟𝑤superscriptsubscript0𝑧𝑥superscriptsubscript0𝑥subscript𝑗𝑥𝑦𝑓𝑦𝑤1𝐹𝑦differential-d𝑦subscript𝑎𝑥superscriptsubscript0𝑤𝑓𝑧𝜆superscript𝑒𝜆𝑤𝑧differential-d𝑧differential-d𝑥r(w)=\int_{0}^{\infty}z(x)\left(\int_{0}^{x}j_{x}(y)\frac{f(y+w)}{1-F(y)}\,% \mathrm{d}y+a_{x}\int_{0}^{w}f(z)\lambda e^{-\lambda(w-z)}\,\mathrm{d}z\right)% \mathrm{d}xitalic_r ( italic_w ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_z ( italic_x ) ( ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT italic_j start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ( italic_y ) divide start_ARG italic_f ( italic_y + italic_w ) end_ARG start_ARG 1 - italic_F ( italic_y ) end_ARG roman_d italic_y + italic_a start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_w end_POSTSUPERSCRIPT italic_f ( italic_z ) italic_λ italic_e start_POSTSUPERSCRIPT - italic_λ ( italic_w - italic_z ) end_POSTSUPERSCRIPT roman_d italic_z ) roman_d italic_x (18)

based on the decomposition with respect to the previous value of the peaks process Znsubscript𝑍𝑛Z_{n}italic_Z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. The integral 0wf(z)λeλ(wz)dzsuperscriptsubscript0𝑤𝑓𝑧𝜆superscript𝑒𝜆𝑤𝑧differential-d𝑧\int_{0}^{w}f(z)\lambda e^{-\lambda(w-z)}\,\mathrm{d}z∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_w end_POSTSUPERSCRIPT italic_f ( italic_z ) italic_λ italic_e start_POSTSUPERSCRIPT - italic_λ ( italic_w - italic_z ) end_POSTSUPERSCRIPT roman_d italic_z is the density of the convolution of the density f𝑓fitalic_f with an independent exponential distribution of parameter λ𝜆\lambdaitalic_λ. On the right-hand side of (18), one can use the definitions of the density z𝑧zitalic_z given by (14) and the homogeneous jump distribution jxsubscript𝑗𝑥j_{x}italic_j start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT and axsubscript𝑎𝑥a_{x}italic_a start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT given by (8). Then in the numerator after the exchange of the order of integrations in the first term and by using the formula for the stationary distribution given by (9) one gets the numerator of (17). In the denominator one can use the equality

0f(y)exp(λy(1F(u))du)dy=0h(y)(1F(y))dy+m.superscriptsubscript0𝑓𝑦𝜆superscriptsubscript𝑦1𝐹𝑢differential-d𝑢differential-d𝑦superscriptsubscript0𝑦1𝐹𝑦differential-d𝑦𝑚\int_{0}^{\infty}f(y)\exp\left(-\lambda\int_{y}^{\infty}(1-F(u))\,\mathrm{d}u% \right)\mathrm{d}y=\int_{0}^{\infty}h(y)(1-F(y))\,\mathrm{d}y+m.∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_f ( italic_y ) roman_exp ( - italic_λ ∫ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F ( italic_u ) ) roman_d italic_u ) roman_d italic_y = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) ( 1 - italic_F ( italic_y ) ) roman_d italic_y + italic_m . (19)

which follows by integration by parts.

Note also that the density is not equal to the remaining reign length of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT under the stationary distribution because it would involve the integral 0h(y)f(y+w)/(1F(y))dysuperscriptsubscript0𝑦𝑓𝑦𝑤1𝐹𝑦differential-d𝑦\int_{0}^{\infty}h(y)f(y+w)/(1-F(y))\,\mathrm{d}y∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_h ( italic_y ) italic_f ( italic_y + italic_w ) / ( 1 - italic_F ( italic_y ) ) roman_d italic_y in place of the first term in the numerator on the right-hand side of (17).

1.6 Entry age parameter

Next we introduce another parameter which we call the entry age and we denote it by E𝐸Eitalic_E. As opposed to our original model we consider individuals as being born at age E𝐸Eitalic_E at a modified birth rate and with a modified lifespan distribution. As a result we obtain a model to the age of the world’s oldest person all values of the entry age parameter E0𝐸0E\geq 0italic_E ≥ 0 and we can fit the model parameters with different values of the entry age.

For any value E0𝐸0E\geq 0italic_E ≥ 0 of the entry age we denote by λE(t)subscript𝜆𝐸𝑡\lambda_{E}(t)italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) the rate at which people reach the age E𝐸Eitalic_E at time t𝑡titalic_t, that is,

λE(t)=λ(tE)(1FtE(E))subscript𝜆𝐸𝑡𝜆𝑡𝐸1subscript𝐹𝑡𝐸𝐸\lambda_{E}(t)=\lambda(t-E)(1-F_{t-E}(E))italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) = italic_λ ( italic_t - italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_E ) ) (20)

because the new birth process of rate λE(t)subscript𝜆𝐸𝑡\lambda_{E}(t)italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) is obtained by an inhomogeneous thinning of the original Poisson process of the birth events. The lifespan distribution of those born at time t𝑡titalic_t with age E𝐸Eitalic_E becomes the remaining lifetime distribution at age E𝐸Eitalic_E. The modified cumulative distribution function and density are given by

FtE(x)=FtE(x+E)FtE(E)1FtE(E),ftE(x)=ftE(x+E)1FtE(E).formulae-sequencesuperscriptsubscript𝐹𝑡𝐸𝑥subscript𝐹𝑡𝐸𝑥𝐸subscript𝐹𝑡𝐸𝐸1subscript𝐹𝑡𝐸𝐸superscriptsubscript𝑓𝑡𝐸𝑥subscript𝑓𝑡𝐸𝑥𝐸1subscript𝐹𝑡𝐸𝐸F_{t}^{E}(x)=\frac{F_{t-E}(x+E)-F_{t-E}(E)}{1-F_{t-E}(E)},\qquad f_{t}^{E}(x)=% \frac{f_{t-E}(x+E)}{1-F_{t-E}(E)}.italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_x ) = divide start_ARG italic_F start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_x + italic_E ) - italic_F start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_E ) end_ARG start_ARG 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_E ) end_ARG , italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_x ) = divide start_ARG italic_f start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_x + italic_E ) end_ARG start_ARG 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_E end_POSTSUBSCRIPT ( italic_E ) end_ARG . (21)

2 Model specification and parameter fitting

In Subsection 2.1 we specify the general model introduced and discussed in Section 1, that is, we assume that the birth rate parameter increases exponentially in time and that the lifespan distribution is given by the gamma–Gompertz–Makeham distribution with time-dependent parameters. We provide the details of the computation of the likelihood as a function of the model parameters in Subsection 2.2. We show the way to maximize the likelihood and how the optimal parameters can be found using the Nelder–Mead method in Subsection 2.3. With these values of the parameters, the age of the world’s oldest person and the reign length of the record holder can be computed as described in Subsection 2.4.

2.1 Model specification: birth rate parameter and lifespan distribution

For the rest of the paper we specifiy our general model described in Section 1 to the following choice of the birth rate parameter and of the lifespan distribution. We choose the value of the entry age to be E=0,30,60𝐸03060E=0,30,60italic_E = 0 , 30 , 60 and we fit the model parameters for all three values of E𝐸Eitalic_E separately.

First we specify the intensity function of the Poisson process of births with an exponential growth in time. For any of the three values of the entry age E𝐸Eitalic_E we assume that the birth rate at age E𝐸Eitalic_E is given by

λE(t)=CEeκEtsubscript𝜆𝐸𝑡subscript𝐶𝐸superscript𝑒subscript𝜅𝐸𝑡\lambda_{E}(t)=C_{E}e^{\kappa_{E}t}italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) = italic_C start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_κ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT italic_t end_POSTSUPERSCRIPT (22)

where the numerical values of the parameters CEsubscript𝐶𝐸C_{E}italic_C start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT and κEsubscript𝜅𝐸\kappa_{E}italic_κ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT are obtained by linear regression of the logarithmic data of newborns, people at age 30303030 and 60606060 published by the United Nations since 1950. We extrapolate the linear regression backwards in time and we use the numerical values shown in Table 1.

We assume that the underlying force of mortality is chosen so that the lifespan distribution of individuals follows the gamma–Gompertz (ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G) distribution with cumulative distribution function and density

Fa,b,γΓG(x)subscriptsuperscript𝐹ΓG𝑎𝑏𝛾𝑥\displaystyle F^{\operatorname{\Gamma G}}_{a,b,\gamma}(x)italic_F start_POSTSUPERSCRIPT roman_Γ roman_G end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_γ end_POSTSUBSCRIPT ( italic_x ) =1(1+aγb(ebx1))1/γ,absent1superscript1𝑎𝛾𝑏superscript𝑒𝑏𝑥11𝛾\displaystyle=1-\left(1+\frac{a\gamma}{b}(e^{bx}-1)\right)^{-1/\gamma},= 1 - ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ( italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT - 1 ) ) start_POSTSUPERSCRIPT - 1 / italic_γ end_POSTSUPERSCRIPT , (23)
fa,b,γΓG(x)subscriptsuperscript𝑓ΓG𝑎𝑏𝛾𝑥\displaystyle f^{\operatorname{\Gamma G}}_{a,b,\gamma}(x)italic_f start_POSTSUPERSCRIPT roman_Γ roman_G end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_γ end_POSTSUBSCRIPT ( italic_x ) =aebx(1+aγb(ebx1))11γabsent𝑎superscript𝑒𝑏𝑥superscript1𝑎𝛾𝑏superscript𝑒𝑏𝑥111𝛾\displaystyle=ae^{bx}\left(1+\frac{a\gamma}{b}(e^{bx}-1)\right)^{-1-\frac{1}{% \gamma}}= italic_a italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ( italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT - 1 ) ) start_POSTSUPERSCRIPT - 1 - divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG end_POSTSUPERSCRIPT

for x0𝑥0x\geq 0italic_x ≥ 0 where a,b,γ𝑎𝑏𝛾a,b,\gammaitalic_a , italic_b , italic_γ are positive parameters. We mention that the gamma–Gompertz–Makeham (ΓGMΓGM\operatorname{\Gamma GM}roman_Γ roman_GM) distribution differs from the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G distribution by the presence of a non-negative extrinsic mortality parameter c𝑐citalic_c which appears as an additive term in the force of mortality. See (29) for the definition of the ΓGMΓGM\operatorname{\Gamma GM}roman_Γ roman_GM distribution. In our model, we exclude the extrinsic mortality for the following two reasons. Since the extrinsic mortality becomes irrelevant at high ages and we aim to model the front-end of the death distribution at the oldest-old ages, we do not expect to obtain a reliable estimate on the extrinsic mortality using the data about the world’s oldest person. On the other hand, as explained later in Subsection 2.3, the likelihood maximization provides unrealistic lifespan distributions even for the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G model if one tries to optimize in all the parameters at the same time. In order for the algorithm to result in a distribution close to the actual human lifespan distribution, the number of model parameters had to be decreased.

For our model we suppose that in the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G lifespan distribution, parameters b=bE𝑏subscript𝑏𝐸b=b_{E}italic_b = italic_b start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT, the rate of aging and γ=γE𝛾subscript𝛾𝐸\gamma=\gamma_{E}italic_γ = italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT, the magnitude of heterogeneity are constants over time and that they only depend on the value of the entry age parameter E𝐸Eitalic_E. The parameter a𝑎aitalic_a, the initial level of mortality at the entry age for individuals born at time t𝑡titalic_t, depends on time given by the exponentially decreasing function

aE(t)=KEeαE(t2000)subscript𝑎𝐸𝑡subscript𝐾𝐸superscript𝑒subscript𝛼𝐸𝑡2000a_{E}(t)=K_{E}e^{-\alpha_{E}(t-2000)}italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) = italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT - italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t - 2000 ) end_POSTSUPERSCRIPT (24)

where the exponent αEsubscript𝛼𝐸\alpha_{E}italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT and the constant KEsubscript𝐾𝐸K_{E}italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT only depends on the entry age E𝐸Eitalic_E. The reason for subtracting 2000200020002000 in (24) is only technical, the numerical values of the parameters do not become tiny with this definition.

In the model with entry age E𝐸Eitalic_E, we assume that the birth rate λE(t)subscript𝜆𝐸𝑡\lambda_{E}(t)italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) is given by (22) and we fit the gamma–Gompetz distribution with parameters bE,γEsubscript𝑏𝐸subscript𝛾𝐸b_{E},\gamma_{E}italic_b start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT and aE(t)subscript𝑎𝐸𝑡a_{E}(t)italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t ) given by (24) for the modified distribution function FtE(x)superscriptsubscript𝐹𝑡𝐸𝑥F_{t}^{E}(x)italic_F start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_x ) and density ftE(x)superscriptsubscript𝑓𝑡𝐸𝑥f_{t}^{E}(x)italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_x ) in (21). This means that we search for the best fitting values of the parameters αE,KE,bE,γEsubscript𝛼𝐸subscript𝐾𝐸subscript𝑏𝐸subscript𝛾𝐸\alpha_{E},K_{E},b_{E},\gamma_{E}italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , italic_b start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT which results in an approximation of the remaining lifetime distribution at the age E𝐸Eitalic_E.

2.2 Likelihood calculations

The aim of the maximum likelihood method is to give an estimate to the parameters αEsubscript𝛼𝐸\alpha_{E}italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT, KEsubscript𝐾𝐸K_{E}italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT, bEsubscript𝑏𝐸b_{E}italic_b start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT and γEsubscript𝛾𝐸\gamma_{E}italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT for E=0,30,60𝐸03060E=0,30,60italic_E = 0 , 30 , 60 by finding those values for which the likelihood of the full sample is the largest. The sample is obtained from the the historical data on the world’s oldest person available in [25]. We transform this information into a list of triples of the form (ti,yi,zi)subscript𝑡𝑖subscript𝑦𝑖subscript𝑧𝑖(t_{i},y_{i},z_{i})( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) for i=1,,n𝑖1𝑛i=1,\dots,nitalic_i = 1 , … , italic_n where tisubscript𝑡𝑖t_{i}italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the i𝑖iitalic_ith time in the sample when the oldest person dies at age yisubscript𝑦𝑖y_{i}italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and the new record holder has age zisubscript𝑧𝑖z_{i}italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT at time tisubscript𝑡𝑖t_{i}italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT. Then the data has to satisfy the consistency relation tizi=ti+1yi+1subscript𝑡𝑖subscript𝑧𝑖subscript𝑡𝑖1subscript𝑦𝑖1t_{i}-z_{i}=t_{i+1}-y_{i+1}italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_t start_POSTSUBSCRIPT italic_i + 1 end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i + 1 end_POSTSUBSCRIPT since the two sides express the date of birth of the same person.

In the model with entry age E𝐸Eitalic_E, the likelihood of the i𝑖iitalic_ith data point (ti,yi,zi)subscript𝑡𝑖subscript𝑦𝑖subscript𝑧𝑖(t_{i},y_{i},z_{i})( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) given the previous data point is equal to

ftiyi+EE(yiE)1Ftiyi+EE(zi1E)jyi,tiE(zi)=ftiyi+EE(yiE)1Ftiyi+EE(zi1E)exp(ziyiλE(tiu+E)(1Ftiu+EE(uE))du)λE(tizi+E)(1Ftizi+EE(ziE))superscriptsubscript𝑓subscript𝑡𝑖subscript𝑦𝑖𝐸𝐸subscript𝑦𝑖𝐸1superscriptsubscript𝐹subscript𝑡𝑖subscript𝑦𝑖𝐸𝐸subscript𝑧𝑖1𝐸subscriptsuperscript𝑗𝐸subscript𝑦𝑖subscript𝑡𝑖subscript𝑧𝑖superscriptsubscript𝑓subscript𝑡𝑖subscript𝑦𝑖𝐸𝐸subscript𝑦𝑖𝐸1superscriptsubscript𝐹subscript𝑡𝑖subscript𝑦𝑖𝐸𝐸subscript𝑧𝑖1𝐸superscriptsubscriptsubscript𝑧𝑖subscript𝑦𝑖subscript𝜆𝐸subscript𝑡𝑖𝑢𝐸1superscriptsubscript𝐹subscript𝑡𝑖𝑢𝐸𝐸𝑢𝐸differential-d𝑢subscript𝜆𝐸subscript𝑡𝑖subscript𝑧𝑖𝐸1superscriptsubscript𝐹subscript𝑡𝑖subscript𝑧𝑖𝐸𝐸subscript𝑧𝑖𝐸\frac{f_{t_{i}-y_{i}+E}^{E}(y_{i}-E)}{1-F_{t_{i}-y_{i}+E}^{E}(z_{i-1}-E)}j^{E}% _{y_{i},t_{i}}(z_{i})\\ =\frac{f_{t_{i}-y_{i}+E}^{E}(y_{i}-E)}{1-F_{t_{i}-y_{i}+E}^{E}(z_{i-1}-E)}\exp% \left(-\int_{z_{i}}^{y_{i}}\lambda_{E}(t_{i}-u+E)(1-F_{t_{i}-u+E}^{E}(u-E))\,% \mathrm{d}u\right)\lambda_{E}(t_{i}-z_{i}+E)(1-F_{t_{i}-z_{i}+E}^{E}(z_{i}-E))start_ROW start_CELL divide start_ARG italic_f start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) end_ARG start_ARG 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i - 1 end_POSTSUBSCRIPT - italic_E ) end_ARG italic_j start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_CELL end_ROW start_ROW start_CELL = divide start_ARG italic_f start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) end_ARG start_ARG 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i - 1 end_POSTSUBSCRIPT - italic_E ) end_ARG roman_exp ( - ∫ start_POSTSUBSCRIPT italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_u - italic_E ) ) roman_d italic_u ) italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) ) end_CELL end_ROW (25)

for all i=2,3,,n𝑖23𝑛i=2,3,\dots,nitalic_i = 2 , 3 , … , italic_n except for i=1𝑖1i=1italic_i = 1 in which case the 1Fy1t1+EE(z0E)1superscriptsubscript𝐹subscript𝑦1subscript𝑡1𝐸𝐸subscript𝑧0𝐸1-F_{y_{1}-t_{1}+E}^{E}(z_{0}-E)1 - italic_F start_POSTSUBSCRIPT italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_z start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_E ) factor in the denominator is missing. In (25) above we use the transition probabilities of the model with entry age E𝐸Eitalic_E given by

jy,tE(x)=exp(xyλE(tu+E)(1Ftu+EE(uE))du)λE(tx+E)(1Ftx+EE(xE))subscriptsuperscript𝑗𝐸𝑦𝑡𝑥superscriptsubscript𝑥𝑦subscript𝜆𝐸𝑡𝑢𝐸1superscriptsubscript𝐹𝑡𝑢𝐸𝐸𝑢𝐸differential-d𝑢subscript𝜆𝐸𝑡𝑥𝐸1superscriptsubscript𝐹𝑡𝑥𝐸𝐸𝑥𝐸j^{E}_{y,t}(x)=\exp\left(-\int_{x}^{y}\lambda_{E}(t-u+E)(1-F_{t-u+E}^{E}(u-E))% \,\mathrm{d}u\right)\lambda_{E}(t-x+E)(1-F_{t-x+E}^{E}(x-E))italic_j start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_y , italic_t end_POSTSUBSCRIPT ( italic_x ) = roman_exp ( - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y end_POSTSUPERSCRIPT italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t - italic_u + italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_u + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_u - italic_E ) ) roman_d italic_u ) italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t - italic_x + italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_x + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_x - italic_E ) ) (26)

as a generalization of (2). The explanation of the left-hand side of (25) is that the person died at time tisubscript𝑡𝑖t_{i}italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT at age yisubscript𝑦𝑖y_{i}italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT had age E𝐸Eitalic_E at time tiyi+Esubscript𝑡𝑖subscript𝑦𝑖𝐸t_{i}-y_{i}+Eitalic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E. The previous data point ensures that this person has already reached age zi1subscript𝑧𝑖1z_{i-1}italic_z start_POSTSUBSCRIPT italic_i - 1 end_POSTSUBSCRIPT hence we condition their lifetime distribution on this fact. The transition probabilities in (26) are obtained similarly to (2) with the difference that a person at age u𝑢uitalic_u with u[x,y]𝑢𝑥𝑦u\in[x,y]italic_u ∈ [ italic_x , italic_y ] at time t𝑡titalic_t had age E𝐸Eitalic_E at time tu+E𝑡𝑢𝐸t-u+Eitalic_t - italic_u + italic_E.

Note that when computing the likelihood of the full data by multiplying the right-hand side of (25) for different values of i𝑖iitalic_i the consistency relation of the data implies that the factor 1Ftizi+E(ziE)1subscript𝐹subscript𝑡𝑖subscript𝑧𝑖𝐸subscript𝑧𝑖𝐸1-F_{t_{i}-z_{i}+E}(z_{i}-E)1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) of the i𝑖iitalic_ith term cancels with the factor 1Fti+1yi+1+E(ziE)1subscript𝐹subscript𝑡𝑖1subscript𝑦𝑖1𝐸subscript𝑧𝑖𝐸1-F_{t_{i+1}-y_{i+1}+E}(z_{i}-E)1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i + 1 end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i + 1 end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) coming from the (i+1)𝑖1(i+1)( italic_i + 1 )st term. Hence the log-likelihood of the full sample is given by

l(α,K,b,γ)𝑙𝛼𝐾𝑏𝛾\displaystyle l(\alpha,K,b,\gamma)italic_l ( italic_α , italic_K , italic_b , italic_γ ) (27)
=i=1n(logftiyi+EE(yiE)ziyiλE(tiu+E)(1Ftiu+EE(uE))du+logλE(tizi+E))+log(1FtnznEE(zn+E))absentsuperscriptsubscript𝑖1𝑛superscriptsubscript𝑓subscript𝑡𝑖subscript𝑦𝑖𝐸𝐸subscript𝑦𝑖𝐸superscriptsubscriptsubscript𝑧𝑖subscript𝑦𝑖subscript𝜆𝐸subscript𝑡𝑖𝑢𝐸1superscriptsubscript𝐹subscript𝑡𝑖𝑢𝐸𝐸𝑢𝐸differential-d𝑢subscript𝜆𝐸subscript𝑡𝑖subscript𝑧𝑖𝐸1superscriptsubscript𝐹subscript𝑡𝑛subscript𝑧𝑛𝐸𝐸subscript𝑧𝑛𝐸\displaystyle=\sum_{i=1}^{n}\left(\log f_{t_{i}-y_{i}+E}^{E}(y_{i}-E)-\int_{z_% {i}}^{y_{i}}\lambda_{E}(t_{i}-u+E)(1-F_{t_{i}-u+E}^{E}(u-E))\,\mathrm{d}u+\log% \lambda_{E}(t_{i}-z_{i}+E)\right)+\log(1-F_{t_{n}-z_{n}-E}^{E}(z_{n}+E))= ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ( roman_log italic_f start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) - ∫ start_POSTSUBSCRIPT italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_u - italic_E ) ) roman_d italic_u + roman_log italic_λ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E ) ) + roman_log ( 1 - italic_F start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT - italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT + italic_E ) )
=i=1n(logfKeα(tiyi+E2000),b,γΓG(yiE)ziyiCeκ(tiu+E)(1FKeα(tiu+E2000),b,γΓG(uE))du)absentsuperscriptsubscript𝑖1𝑛subscriptsuperscript𝑓ΓG𝐾superscript𝑒𝛼subscript𝑡𝑖subscript𝑦𝑖𝐸2000𝑏𝛾subscript𝑦𝑖𝐸superscriptsubscriptsubscript𝑧𝑖subscript𝑦𝑖𝐶superscript𝑒𝜅subscript𝑡𝑖𝑢𝐸1subscriptsuperscript𝐹ΓG𝐾superscript𝑒𝛼subscript𝑡𝑖𝑢𝐸2000𝑏𝛾𝑢𝐸differential-d𝑢\displaystyle=\sum_{i=1}^{n}\left(\log f^{\operatorname{\Gamma G}}_{Ke^{-% \alpha(t_{i}-y_{i}+E-2000)},b,\gamma}(y_{i}-E)-\int_{z_{i}}^{y_{i}}Ce^{\kappa(% t_{i}-u+E)}\left(1-F^{\operatorname{\Gamma G}}_{Ke^{-\alpha(t_{i}-u+E-2000)},b% ,\gamma}(u-E)\right)\mathrm{d}u\right)= ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ( roman_log italic_f start_POSTSUPERSCRIPT roman_Γ roman_G end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_K italic_e start_POSTSUPERSCRIPT - italic_α ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E - 2000 ) end_POSTSUPERSCRIPT , italic_b , italic_γ end_POSTSUBSCRIPT ( italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_E ) - ∫ start_POSTSUBSCRIPT italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_C italic_e start_POSTSUPERSCRIPT italic_κ ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E ) end_POSTSUPERSCRIPT ( 1 - italic_F start_POSTSUPERSCRIPT roman_Γ roman_G end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_K italic_e start_POSTSUPERSCRIPT - italic_α ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_u + italic_E - 2000 ) end_POSTSUPERSCRIPT , italic_b , italic_γ end_POSTSUBSCRIPT ( italic_u - italic_E ) ) roman_d italic_u )
+log(1FKeα(tnzn+E2000),b,γΓG(znE))+nlogC+i=1nκ(tizi+E).1subscriptsuperscript𝐹ΓG𝐾superscript𝑒𝛼subscript𝑡𝑛subscript𝑧𝑛𝐸2000𝑏𝛾subscript𝑧𝑛𝐸𝑛𝐶superscriptsubscript𝑖1𝑛𝜅subscript𝑡𝑖subscript𝑧𝑖𝐸\displaystyle\qquad+\log\left(1-F^{\operatorname{\Gamma G}}_{Ke^{-\alpha(t_{n}% -z_{n}+E-2000)},b,\gamma}(z_{n}-E)\right)+n\log C+\sum_{i=1}^{n}\kappa(t_{i}-z% _{i}+E).+ roman_log ( 1 - italic_F start_POSTSUPERSCRIPT roman_Γ roman_G end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_K italic_e start_POSTSUPERSCRIPT - italic_α ( italic_t start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT + italic_E - 2000 ) end_POSTSUPERSCRIPT , italic_b , italic_γ end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT - italic_E ) ) + italic_n roman_log italic_C + ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT italic_κ ( italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + italic_E ) .

where we suppress the dependence of the parameters α,K,b,γ𝛼𝐾𝑏𝛾\alpha,K,b,\gammaitalic_α , italic_K , italic_b , italic_γ on the entry age. Note that the last two terms do not depend on the parameters α,K,b,γ𝛼𝐾𝑏𝛾\alpha,K,b,\gammaitalic_α , italic_K , italic_b , italic_γ hence we can omit these terms in the maximization of the log-likelihood.

2.3 Likelihood maximization

We implemented the calculation of the log-likelihood function l(α,K,b,γ)𝑙𝛼𝐾𝑏𝛾l(\alpha,K,b,\gamma)italic_l ( italic_α , italic_K , italic_b , italic_γ ) given by (27) in Python. We used numerical integration to obtain the integrals on the right-hand side of (27). We mention that the general integral formula in (30) could not be used because the parameter a𝑎aitalic_a of the gamma–Gompertz distribution in the integrand depends on the integration variable on the right-hand side of (27).

In order to maximize the value of the log-likelihood function l(α,K,b,γ)𝑙𝛼𝐾𝑏𝛾l(\alpha,K,b,\gamma)italic_l ( italic_α , italic_K , italic_b , italic_γ ) we applied the Nelder–Mead method [27] which is already implemented in Python. We mention that initially we used the gamma–Gompertz–Makeham distribution as lifespan distribution, see (29) for the definition, which contains the extra parameter c𝑐citalic_c to be fitted but it turned out that the number of model parameters has to be reduced. The behaviour of the optimization algorithm in the five parameters α,K,b,c,γ𝛼𝐾𝑏𝑐𝛾\alpha,K,b,c,\gammaitalic_α , italic_K , italic_b , italic_c , italic_γ using the gamma–Gompretz–Makeham model was very similar the case of four parameters α,K,c,γ𝛼𝐾𝑐𝛾\alpha,K,c,\gammaitalic_α , italic_K , italic_c , italic_γ in the gamma–Gompertz model. Running the optimization in the full set of parameters (α,K,b,c,γ𝛼𝐾𝑏𝑐𝛾\alpha,K,b,c,\gammaitalic_α , italic_K , italic_b , italic_c , italic_γ in the gamma–Gompertz–Makeham model or α,K,b,γ𝛼𝐾𝑏𝛾\alpha,K,b,\gammaitalic_α , italic_K , italic_b , italic_γ in the gamma–Gompertz model), it turned out that after a few rounds the parameter K𝐾Kitalic_K started to decrease dramatically and reached values below 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT. The resulting lifespan distribution seemed very unrealistic with almost no mortality before the age of 100. This happened for all values of the entry age E=0,30,60𝐸03060E=0,30,60italic_E = 0 , 30 , 60.

We explain this phenomenon by the fact that historical data about the oldest person in the world only gives information about the behaviour of the lifespan distribution between the ages 107 and 123. The simple optimization in the four parameters α,K,b,γ𝛼𝐾𝑏𝛾\alpha,K,b,\gammaitalic_α , italic_K , italic_b , italic_γ simultaneously yields an excellent fit for the tail decay of the lifespan distribution with the historical data but the result may be very far from the actual human lifespan. This would limit the practical relevance of our results.

The mathematical reason for the fact that the four-parameter optimization does not result in a satisfactory approximation to the human lifespan distribution is the following. In these cases, the optimization procedure diverges to those regimes of the parameter space +4superscriptsubscript4\mathbb{R}_{+}^{4}blackboard_R start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT where the corresponding gamma–Gompertz distribution is degenerate. One can prevent reaching these unrealistic combinations of parameters by reducing the amount of freedom in the optimization. Hence we specify some of the parameters a priori and we perform the optimization in the remaining ones so that it provides a good fit to the data on the age of the oldest old as well as a realistic lifespan distribution.

We believe that the most robust of the four parameters of the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G model is b𝑏bitalic_b which is the exponent in the time dependence of the mortality rate. By setting the rate of aging b=0.09𝑏0.09b=0.09italic_b = 0.09 the algorithm gives the optimal triple α,K,γ𝛼𝐾𝛾\alpha,K,\gammaitalic_α , italic_K , italic_γ with the best likelihood which is very stable under changing the initial values of these parameters. The running time is also very short.

The Nelder–Mead algorithm, being a numerical maximization method, heavily relies on the tolerance parameter, which determines the minimal improvement required for the algorithm to continue running. If this parameter is set too high, the algorithm might stop before reaching the optimum. Conversely, if set too low, the algorithm might take excessively long to converge. To address this, we drew inspiration from dynamic learning rate algorithms used in neural network training and developed the following meta-algorithm.

First, we run the Nelder–Mead optimization. Based on the improvement from the starting point, we dynamically adjust the tolerance factor, similar to how learning rates are modified during neural network training. We then run the optimization again, recalibrating the tolerance factor based on the observed improvement, and repeat the process. This iterative adjustment allows us to get closer to the optimum, a hypothesis supported by our practical experience with this meta-algorithm. Following this meta-algorithm, only a few calls of the Nelder–Mead method is enough to reach the optimum. The Python codes for the likelihood calculations as well as the Nelder–Mead optimization implemented to this problem are available in [28].

The numerical values of the resulting parameters for the three choices of the entry age are shown in Table 2. The survival probability functions with the parameters given in Table 2 for individuals born in 2000 corresponding to the entry age E=0,30,60𝐸03060E=0,30,60italic_E = 0 , 30 , 60 are shown on Figure 5 as a function of the age. We also computed the optimal values of the parameters α,K,γ𝛼𝐾𝛾\alpha,K,\gammaitalic_α , italic_K , italic_γ for other values of the rate of aging b𝑏bitalic_b as a sensitivity analysis. The resulting parameter values for the choices b=0.11𝑏0.11b=0.11italic_b = 0.11, b=0.13𝑏0.13b=0.13italic_b = 0.13 and b=0.15𝑏0.15b=0.15italic_b = 0.15 are shown in Table 3.

Refer to caption
Figure 5: Survival probabilities as a function of the age for different values of the entry age: E=0 in red, E=30 in black, E=60 in blue.

We mention as an alternative approach that scaling the parameters could enhance the optimization process, but this requires prior knowledge of the range within which the parameters vary. This range could be determined through our iterative application of the Nelder–Mead algorithm.

2.4 Computation of the oldest person’s age and of the reign length

We observe that the model with the parameters in Table 2 fits well to the titleholder data. We focus on two statistics of the process in order to support this observation about the comparison: the age of the world’s oldest person and the reign length of the record holder. In the case of both statistics exact formulas are only available for the homogeneous model introduced in Subsection 1.3 where the birth rate is constant as well as the lifespan distribution does not depend on time. Hence we apply an approximation where the error is negligible compared to the difference from the statistics computed using the data.

In the general model the distribution of the age of the world’s oldest person at time t𝑡titalic_t is given by the density ht(x)subscript𝑡𝑥h_{t}(x)italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_x ) in (5) and by the point mass mtsubscript𝑚𝑡m_{t}italic_m start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT at 00 in (6). For the numerical computations, we ignore the point mass mtsubscript𝑚𝑡m_{t}italic_m start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT which is below the round-off error in the numerical results. The difficulty in computing the mean age of the oldest person at time t𝑡titalic_t is that parameter a𝑎aitalic_a of the gamma–Gompertz–Makeham distribution function Ftusubscript𝐹𝑡𝑢F_{t-u}italic_F start_POSTSUBSCRIPT italic_t - italic_u end_POSTSUBSCRIPT in the exponent of (5) also depends on the integration variable u𝑢uitalic_u.

In our approximation we fix the value of the parameter a𝑎aitalic_a of the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G distribution in ht(x)subscript𝑡𝑥h_{t}(x)italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_x ) in (5) to a value which is equal to a0(td)subscript𝑎0𝑡𝑑a_{0}(t-d)italic_a start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( italic_t - italic_d ) in (24) with some delay d𝑑ditalic_d. The delay d𝑑ditalic_d is chosen so that the mean age of the oldest person computed using a0(td)subscript𝑎0𝑡𝑑a_{0}(t-d)italic_a start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( italic_t - italic_d ) as parameter a𝑎aitalic_a for all times in the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G distribution function in (5) is equal to the same value d𝑑ditalic_d. For a given t𝑡titalic_t, this value of d𝑑ditalic_d can be obtained as the fixed point of the contraction map

d0xexλ(tu)(1Ftd(u))duλ(tx)(1Ftd(x))dxmaps-to𝑑superscriptsubscript0𝑥superscript𝑒superscriptsubscript𝑥𝜆𝑡𝑢1subscript𝐹𝑡𝑑𝑢differential-d𝑢𝜆𝑡𝑥1subscript𝐹𝑡𝑑𝑥differential-d𝑥d\mapsto\int_{0}^{\infty}xe^{-\int_{x}^{\infty}\lambda(t-u)(1-F_{t-d}(u))\,% \mathrm{d}u}\lambda(t-x)(1-F_{t-d}(x))\,\mathrm{d}xitalic_d ↦ ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_x italic_e start_POSTSUPERSCRIPT - ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_u ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_d end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u end_POSTSUPERSCRIPT italic_λ ( italic_t - italic_x ) ( 1 - italic_F start_POSTSUBSCRIPT italic_t - italic_d end_POSTSUBSCRIPT ( italic_x ) ) roman_d italic_x (28)

which provides a reasonable approximation for the mean age of the world’s oldest person. For the comparison with the data and for the prediction, we use the model with entry age E=0𝐸0E=0italic_E = 0. Hence in (28), the function λ𝜆\lambdaitalic_λ is given in (22) with E=0𝐸0E=0italic_E = 0 and Ftdsubscript𝐹𝑡𝑑F_{t-d}italic_F start_POSTSUBSCRIPT italic_t - italic_d end_POSTSUBSCRIPT is the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G distribution function with parameters given by the E=0𝐸0E=0italic_E = 0 values in Table 2 and with a=a0(td)𝑎subscript𝑎0𝑡𝑑a=a_{0}(t-d)italic_a = italic_a start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( italic_t - italic_d ) in (24). This approximation is reasonable because the distribution of the age of the oldest person is highly concentrated.

The fixed point of the map in (28) as the expected age of the world’s oldest person can be found in a few steps of iterations. We show the result on Figure 1. We applied 10 iterations using the fitted model parameters for each year between 1955 and 2019. By computing the standard deviation of the age of the oldest person as well we obtain the mean and a confidence interval for the age.

The predictions for the age distribution of the world’s oldest person in the future shown on Figure 2. We obtained them by using the exact age distribution formula in Subsection 1.2 along with the numerical values of the parameters C,κ,α,K,γ𝐶𝜅𝛼𝐾𝛾C,\kappa,\alpha,K,\gammaitalic_C , italic_κ , italic_α , italic_K , italic_γ given in Tables 1 and 2 for entry age E=0𝐸0E=0italic_E = 0.

In our model, the distribution function of the age of the world’s oldest person is given in (7) where the ΓGΓG\operatorname{\Gamma G}roman_Γ roman_G distribution function can be substituted with the estimated parameter values at any time. In this way, the probability of observing an age greater than or equal to Calment’s or Knauss’ actual age at the time of their death can be computed exactly. The numerical values are 0.0002860.0002860.0002860.000286 for Calment and it is 0.01160.01160.01160.0116 for Knauss.

The backtesting mentioned in the Results section is performed as follows. We estimated the best parameter values with entry age 00 based on the reduced data on the world’s oldest person between 1955 and 1988 where the ending date is the time when Calment became the world’s oldest person. The resulting parameters α=0.01516,K=0.00002064,γ=0.08413formulae-sequence𝛼0.01516formulae-sequence𝐾0.00002064𝛾0.08413\alpha=0.01516,K=0.00002064,\gamma=0.08413italic_α = 0.01516 , italic_K = 0.00002064 , italic_γ = 0.08413 are numerically not very far from the optimal parameters in Table 2 but the difference is more visible on Figure 1. The figure shows the model-based mean age and confidence interval for the age of the world’s oldest person computed using the full data as well as the data until 1988.

For the reign length of record holders, we again used the expected age at a given time obtained as the fixed point of the iteration in (28). The numerical value of the expected reign length obtained from the iteration is 1.1951.1951.1951.195 in 1955 and it is 1.1881.1881.1881.188 in 2019. The empirical value of the reign length is 1.0081.0081.0081.008 computed from the data by dividing the total length of the time interval between 1955 and 2019 by the number of record holders.

3 Methods

In this section, we provide supplementary information related to the main result of this paper. We perform explicit computations with the gamma–Gompertz–Makeham model and we express the integral of the survival function in terms of a hypergeometric function.

The cumulative distribution function and the density of the gamma–Gompertz–Makeham (ΓGMΓGM\operatorname{\Gamma GM}roman_Γ roman_GM) distribution are given by

Fa,b,c,γΓGM(x)subscriptsuperscript𝐹ΓGM𝑎𝑏𝑐𝛾𝑥\displaystyle F^{\operatorname{\Gamma GM}}_{a,b,c,\gamma}(x)italic_F start_POSTSUPERSCRIPT roman_Γ roman_GM end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_c , italic_γ end_POSTSUBSCRIPT ( italic_x ) =1ecx(1+aγb(ebx1))1/γ,absent1superscript𝑒𝑐𝑥superscript1𝑎𝛾𝑏superscript𝑒𝑏𝑥11𝛾\displaystyle=1-\frac{e^{-cx}}{\left(1+\frac{a\gamma}{b}(e^{bx}-1)\right)^{1/% \gamma}},= 1 - divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_c italic_x end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ( italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT - 1 ) ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT end_ARG , (29)
fa,b,c,γΓGM(x)subscriptsuperscript𝑓ΓGM𝑎𝑏𝑐𝛾𝑥\displaystyle f^{\operatorname{\Gamma GM}}_{a,b,c,\gamma}(x)italic_f start_POSTSUPERSCRIPT roman_Γ roman_GM end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_c , italic_γ end_POSTSUBSCRIPT ( italic_x ) =ecx(1+aγb(ebx1))1+1γc(baγ)+a(b+cγ)ebxbabsentsuperscript𝑒𝑐𝑥superscript1𝑎𝛾𝑏superscript𝑒𝑏𝑥111𝛾𝑐𝑏𝑎𝛾𝑎𝑏𝑐𝛾superscript𝑒𝑏𝑥𝑏\displaystyle=\frac{e^{-cx}}{\left(1+\frac{a\gamma}{b}(e^{bx}-1)\right)^{1+% \frac{1}{\gamma}}}\frac{c(b-a\gamma)+a(b+c\gamma)e^{bx}}{b}= divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_c italic_x end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ( italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT - 1 ) ) start_POSTSUPERSCRIPT 1 + divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG end_POSTSUPERSCRIPT end_ARG divide start_ARG italic_c ( italic_b - italic_a italic_γ ) + italic_a ( italic_b + italic_c italic_γ ) italic_e start_POSTSUPERSCRIPT italic_b italic_x end_POSTSUPERSCRIPT end_ARG start_ARG italic_b end_ARG

for x0𝑥0x\geq 0italic_x ≥ 0 where a,b,c,γ𝑎𝑏𝑐𝛾a,b,c,\gammaitalic_a , italic_b , italic_c , italic_γ are positive parameters. The positivity of parameters implies the finiteness of all moments and, in particular, the convergence of the integral of the survival function x(1Fa,b,c,γΓGM(u))dusuperscriptsubscript𝑥1subscriptsuperscript𝐹ΓGM𝑎𝑏𝑐𝛾𝑢differential-d𝑢\int_{x}^{\infty}(1-F^{\operatorname{\Gamma GM}}_{a,b,c,\gamma}(u))\,\mathrm{d}u∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F start_POSTSUPERSCRIPT roman_Γ roman_GM end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_c , italic_γ end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u. In the homogeneous model, the integral of the survival function appears in the density of the distribution of Ytsubscript𝑌𝑡Y_{t}italic_Y start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT in (9) and in the stationary density of the peaks process in (14). We show below that in the gamma–Gompertz–Makeham model the integral of the survival function can be computed explicitly and it is given by

x(1Fa,b,c,γΓGM(u))du=(baγ)1/γe(c+bγ)xbγ+cF12(1γ,1γ+cb;1+1γ+cb;aγbaγebx)superscriptsubscript𝑥1subscriptsuperscript𝐹ΓGM𝑎𝑏𝑐𝛾𝑢differential-d𝑢superscript𝑏𝑎𝛾1𝛾superscript𝑒𝑐𝑏𝛾𝑥𝑏𝛾𝑐subscriptsubscript𝐹121𝛾1𝛾𝑐𝑏11𝛾𝑐𝑏𝑎𝛾𝑏𝑎𝛾superscript𝑒𝑏𝑥\int_{x}^{\infty}\left(1-F^{\operatorname{\Gamma GM}}_{a,b,c,\gamma}(u)\right)% \mathrm{d}u=\left(\frac{b}{a\gamma}\right)^{1/\gamma}\frac{e^{-(c+\frac{b}{% \gamma})x}}{\frac{b}{\gamma}+c}{}_{2}F_{1}\left(\frac{1}{\gamma},\frac{1}{% \gamma}+\frac{c}{b};1+\frac{1}{\gamma}+\frac{c}{b};\frac{a\gamma-b}{a\gamma}e^% {-bx}\right)∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F start_POSTSUPERSCRIPT roman_Γ roman_GM end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_c , italic_γ end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u = ( divide start_ARG italic_b end_ARG start_ARG italic_a italic_γ end_ARG ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT divide start_ARG italic_e start_POSTSUPERSCRIPT - ( italic_c + divide start_ARG italic_b end_ARG start_ARG italic_γ end_ARG ) italic_x end_POSTSUPERSCRIPT end_ARG start_ARG divide start_ARG italic_b end_ARG start_ARG italic_γ end_ARG + italic_c end_ARG start_FLOATSUBSCRIPT 2 end_FLOATSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG , divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ; 1 + divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ; divide start_ARG italic_a italic_γ - italic_b end_ARG start_ARG italic_a italic_γ end_ARG italic_e start_POSTSUPERSCRIPT - italic_b italic_x end_POSTSUPERSCRIPT ) (30)

where F12(a,b;c;z)subscriptsubscript𝐹12𝑎𝑏𝑐𝑧{}_{2}F_{1}(a,b;c;z)start_FLOATSUBSCRIPT 2 end_FLOATSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_a , italic_b ; italic_c ; italic_z ) is the hypergeometric function. See 15.1.1 in [29] for the definition and properties.

We prove (30) based on the following integral representation 15.3.1 in [29] of the hypergeometric function

2F1(a,b,c,z)=Γ(c)Γ(b)Γ(cb)01tb1(1t)cb1(1tz)adt_{2}F_{1}(a,b,c,z)=\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}t^{b-1}(1% -t)^{c-b-1}(1-tz)^{-a}\,\mathrm{d}tstart_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_a , italic_b , italic_c , italic_z ) = divide start_ARG roman_Γ ( italic_c ) end_ARG start_ARG roman_Γ ( italic_b ) roman_Γ ( italic_c - italic_b ) end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT italic_t start_POSTSUPERSCRIPT italic_b - 1 end_POSTSUPERSCRIPT ( 1 - italic_t ) start_POSTSUPERSCRIPT italic_c - italic_b - 1 end_POSTSUPERSCRIPT ( 1 - italic_t italic_z ) start_POSTSUPERSCRIPT - italic_a end_POSTSUPERSCRIPT roman_d italic_t (31)

which holds whenever Re(c)>Re(b)>0Re𝑐Re𝑏0\operatorname{Re}(c)>\operatorname{Re}(b)>0roman_Re ( italic_c ) > roman_Re ( italic_b ) > 0. First we prove an identity for complex parameters α,β,δ𝛼𝛽𝛿\alpha,\beta,\deltaitalic_α , italic_β , italic_δ which satisfy Re(α+β)>0Re𝛼𝛽0\operatorname{Re}(\alpha+\beta)>0roman_Re ( italic_α + italic_β ) > 0 and we compute

xeβu(1+δeu)αdusuperscriptsubscript𝑥superscript𝑒𝛽𝑢superscript1𝛿superscript𝑒𝑢𝛼differential-d𝑢\displaystyle\int_{x}^{\infty}\frac{e^{-\beta u}}{(1+\delta e^{u})^{\alpha}}\,% \mathrm{d}u∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_β italic_u end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + italic_δ italic_e start_POSTSUPERSCRIPT italic_u end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG roman_d italic_u =1δαxe(α+β)u(1+eu/δ)αduabsent1superscript𝛿𝛼superscriptsubscript𝑥superscript𝑒𝛼𝛽𝑢superscript1superscript𝑒𝑢𝛿𝛼differential-d𝑢\displaystyle=\frac{1}{\delta^{\alpha}}\int_{x}^{\infty}\frac{e^{-(\alpha+% \beta)u}}{(1+e^{-u}/\delta)^{\alpha}}\,\mathrm{d}u= divide start_ARG 1 end_ARG start_ARG italic_δ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT divide start_ARG italic_e start_POSTSUPERSCRIPT - ( italic_α + italic_β ) italic_u end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + italic_e start_POSTSUPERSCRIPT - italic_u end_POSTSUPERSCRIPT / italic_δ ) start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG roman_d italic_u (32)
=e(α+β)xδα01yα+β1(1+exy/δ)αdyabsentsuperscript𝑒𝛼𝛽𝑥superscript𝛿𝛼superscriptsubscript01superscript𝑦𝛼𝛽1superscript1superscript𝑒𝑥𝑦𝛿𝛼differential-d𝑦\displaystyle=\frac{e^{-(\alpha+\beta)x}}{\delta^{\alpha}}\int_{0}^{1}y^{% \alpha+\beta-1}\left(1+e^{-x}y/\delta\right)^{-\alpha}\,\mathrm{d}y= divide start_ARG italic_e start_POSTSUPERSCRIPT - ( italic_α + italic_β ) italic_x end_POSTSUPERSCRIPT end_ARG start_ARG italic_δ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT italic_y start_POSTSUPERSCRIPT italic_α + italic_β - 1 end_POSTSUPERSCRIPT ( 1 + italic_e start_POSTSUPERSCRIPT - italic_x end_POSTSUPERSCRIPT italic_y / italic_δ ) start_POSTSUPERSCRIPT - italic_α end_POSTSUPERSCRIPT roman_d italic_y
=e(α+β)x(α+β)δαF12(α,α+β;1+α+β;exδ)absentsuperscript𝑒𝛼𝛽𝑥𝛼𝛽superscript𝛿𝛼subscriptsubscript𝐹12𝛼𝛼𝛽1𝛼𝛽superscript𝑒𝑥𝛿\displaystyle=\frac{e^{-(\alpha+\beta)x}}{(\alpha+\beta)\delta^{\alpha}}\cdot{% }_{2}F_{1}\left(\alpha,\alpha+\beta;1+\alpha+\beta;-\frac{e^{-x}}{\delta}\right)= divide start_ARG italic_e start_POSTSUPERSCRIPT - ( italic_α + italic_β ) italic_x end_POSTSUPERSCRIPT end_ARG start_ARG ( italic_α + italic_β ) italic_δ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG ⋅ start_FLOATSUBSCRIPT 2 end_FLOATSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_α , italic_α + italic_β ; 1 + italic_α + italic_β ; - divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_x end_POSTSUPERSCRIPT end_ARG start_ARG italic_δ end_ARG )

where we applied a change of variables y=exu𝑦superscript𝑒𝑥𝑢y=e^{x-u}italic_y = italic_e start_POSTSUPERSCRIPT italic_x - italic_u end_POSTSUPERSCRIPT in the second equality above and we applied the hypergeometric identity (31) in the last equality with a=α𝑎𝛼a=\alphaitalic_a = italic_α, b=α+β𝑏𝛼𝛽b=\alpha+\betaitalic_b = italic_α + italic_β, c=1+α+β𝑐1𝛼𝛽c=1+\alpha+\betaitalic_c = 1 + italic_α + italic_β, z=ex/δ𝑧superscript𝑒𝑥𝛿z=-e^{-x}/\deltaitalic_z = - italic_e start_POSTSUPERSCRIPT - italic_x end_POSTSUPERSCRIPT / italic_δ together with the observation that with these values of the parameters the prefactor of the integral on the right-hand side of (31) simplifies to α+β𝛼𝛽\alpha+\betaitalic_α + italic_β. Note that the condition Re(c)>Re(b)>0Re𝑐Re𝑏0\operatorname{Re}(c)>\operatorname{Re}(b)>0roman_Re ( italic_c ) > roman_Re ( italic_b ) > 0 for (31) to hold is satisfied by our assumption Re(α+β)>0Re𝛼𝛽0\operatorname{Re}(\alpha+\beta)>0roman_Re ( italic_α + italic_β ) > 0 which also makes the integrals in (32) convergent.

Next we show (30) using (32) as follows. We write

x(1Fa,b,c,γΓGM(u))dusuperscriptsubscript𝑥1subscriptsuperscript𝐹ΓGM𝑎𝑏𝑐𝛾𝑢differential-d𝑢\displaystyle\int_{x}^{\infty}\left(1-F^{\operatorname{\Gamma GM}}_{a,b,c,% \gamma}(u)\right)\,\mathrm{d}u∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ( 1 - italic_F start_POSTSUPERSCRIPT roman_Γ roman_GM end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a , italic_b , italic_c , italic_γ end_POSTSUBSCRIPT ( italic_u ) ) roman_d italic_u =1(1aγb)1/γxecu(1+aγbaγebu)1/γduabsent1superscript1𝑎𝛾𝑏1𝛾superscriptsubscript𝑥superscript𝑒𝑐𝑢superscript1𝑎𝛾𝑏𝑎𝛾superscript𝑒𝑏𝑢1𝛾differential-d𝑢\displaystyle=\frac{1}{\left(1-\frac{a\gamma}{b}\right)^{1/\gamma}}\int_{x}^{% \infty}\frac{e^{-cu}}{\left(1+\frac{a\gamma}{b-a\gamma}e^{bu}\right)^{1/\gamma% }}\,\mathrm{d}u= divide start_ARG 1 end_ARG start_ARG ( 1 - divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT end_ARG ∫ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_c italic_u end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b - italic_a italic_γ end_ARG italic_e start_POSTSUPERSCRIPT italic_b italic_u end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT end_ARG roman_d italic_u (33)
=1(1aγb)1/γbbxecv/b(1+aγbaγev)1/γdvabsent1superscript1𝑎𝛾𝑏1𝛾𝑏superscriptsubscript𝑏𝑥superscript𝑒𝑐𝑣𝑏superscript1𝑎𝛾𝑏𝑎𝛾superscript𝑒𝑣1𝛾differential-d𝑣\displaystyle=\frac{1}{\left(1-\frac{a\gamma}{b}\right)^{1/\gamma}b}\int_{bx}^% {\infty}\frac{e^{-cv/b}}{\left(1+\frac{a\gamma}{b-a\gamma}e^{v}\right)^{1/% \gamma}}\,\mathrm{d}v= divide start_ARG 1 end_ARG start_ARG ( 1 - divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT italic_b end_ARG ∫ start_POSTSUBSCRIPT italic_b italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT divide start_ARG italic_e start_POSTSUPERSCRIPT - italic_c italic_v / italic_b end_POSTSUPERSCRIPT end_ARG start_ARG ( 1 + divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b - italic_a italic_γ end_ARG italic_e start_POSTSUPERSCRIPT italic_v end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT end_ARG roman_d italic_v
=1(1aγb)1/γbe(1γ+cb)bx(1γ+cb)(aγbaγ)1/γF12(1γ,1γ+cb;1+1γ+cb;aγbaγebx)absent1superscript1𝑎𝛾𝑏1𝛾𝑏superscript𝑒1𝛾𝑐𝑏𝑏𝑥1𝛾𝑐𝑏superscript𝑎𝛾𝑏𝑎𝛾1𝛾subscriptsubscript𝐹121𝛾1𝛾𝑐𝑏11𝛾𝑐𝑏𝑎𝛾𝑏𝑎𝛾superscript𝑒𝑏𝑥\displaystyle=\frac{1}{\left(1-\frac{a\gamma}{b}\right)^{1/\gamma}b}\frac{e^{-% \left(\frac{1}{\gamma}+\frac{c}{b}\right)bx}}{\left(\frac{1}{\gamma}+\frac{c}{% b}\right)\left(\frac{a\gamma}{b-a\gamma}\right)^{1/\gamma}}{}_{2}F_{1}\left(% \frac{1}{\gamma},\frac{1}{\gamma}+\frac{c}{b};1+\frac{1}{\gamma}+\frac{c}{b};% \frac{a\gamma-b}{a\gamma}e^{-bx}\right)= divide start_ARG 1 end_ARG start_ARG ( 1 - divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b end_ARG ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT italic_b end_ARG divide start_ARG italic_e start_POSTSUPERSCRIPT - ( divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ) italic_b italic_x end_POSTSUPERSCRIPT end_ARG start_ARG ( divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ) ( divide start_ARG italic_a italic_γ end_ARG start_ARG italic_b - italic_a italic_γ end_ARG ) start_POSTSUPERSCRIPT 1 / italic_γ end_POSTSUPERSCRIPT end_ARG start_FLOATSUBSCRIPT 2 end_FLOATSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG , divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ; 1 + divide start_ARG 1 end_ARG start_ARG italic_γ end_ARG + divide start_ARG italic_c end_ARG start_ARG italic_b end_ARG ; divide start_ARG italic_a italic_γ - italic_b end_ARG start_ARG italic_a italic_γ end_ARG italic_e start_POSTSUPERSCRIPT - italic_b italic_x end_POSTSUPERSCRIPT )

where we applied the change of variables v=bu𝑣𝑏𝑢v=buitalic_v = italic_b italic_u in the second equality above and we used (32) with α=1/γ𝛼1𝛾\alpha=1/\gammaitalic_α = 1 / italic_γ, β=c/b𝛽𝑐𝑏\beta=c/bitalic_β = italic_c / italic_b, δ=aγ/(baγ)𝛿𝑎𝛾𝑏𝑎𝛾\delta=a\gamma/(b-a\gamma)italic_δ = italic_a italic_γ / ( italic_b - italic_a italic_γ ) and with x𝑥xitalic_x replaced by bx𝑏𝑥bxitalic_b italic_x. The right-hand side of (33) simplifies to that of (30).

References

  • [1] Oeppen, J. & Vaupel, J. W. Broken limits to life expectancy. \JournalTitleScience 296, 1029–1031 (2002).
  • [2] Canudas-Romo, V. Three measures of longevity: Time trends and record values. \JournalTitleDemography 47, 299–312 (2010).
  • [3] Vaupel, J. W., Villavicencio, F. & Bergeron-Boucher, M.-P. Demographic perspectives on the rise of longevity. \JournalTitleProceedings of the National Academy of Sciences 118, e2019536118 (2021).
  • [4] Gompertz, B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. \JournalTitlePhilosophical Transactions of the Royal Society of London 115, 513–583 (1825).
  • [5] Németh, L. & Missov, T. I. Adequate life-expectancy reconstruction for adult human mortality data. \JournalTitlePlos one 13, e0198485 (2018).
  • [6] Missov, T. I., Németh, L. & Dańko, M. J. How much can we trust life tables? sensitivity of mortality measures to right-censoring treatment. \JournalTitlePalgrave Communications 2, 1–10 (2016).
  • [7] Maier, H., Jeune, B. & Vaupel, J. W. Exceptional lifespans (Springer Nature, 2021).
  • [8] Dang, L. H. K. et al. The question of the human mortality plateau. \JournalTitleDemographic Research 48, 321–338 (2023).
  • [9] Vijg, J. & Le Bourg, E. Aging and the Inevitable Limit to Human Life Span. \JournalTitleGerontology 63, 432–434 (2017).
  • [10] Alvarez, J.-A., Villavicencio, F., Strozza, C. & Camarda, C. G. Regularities in human mortality after age 105. \JournalTitlePloS one 16, e0253940 (2021).
  • [11] Belzile, L. R., Davison, A. C., Gampe, J., Rootzén, H. & Zholud, D. Is there a cap on longevity? A statistical review. \JournalTitleAnnual Review of Statistics and Its Application 9, 21–45 (2022).
  • [12] Rootzén, H. & Zholud, D. Human life is unlimited–but short. \JournalTitleExtremes 20, 713–728 (2017).
  • [13] Barbi, E., Lagona, F., Marsili, M., Vaupel, J. W. & Wachter, K. W. The plateau of human mortality: Demography of longevity pioneers. \JournalTitleScience 360, 1459–1461 (2018).
  • [14] Modig, K., Andersson, T., Vaupel, J., Rau, R. & Ahlbom, A. How long do centenarians survive? Life expectancy and maximum lifespan. \JournalTitleJournal of internal medicine 282, 156–163 (2017).
  • [15] Wilmoth, J. R. & Robine, J.-M. The world trend in maximum life span. \JournalTitlePopulation and Development Review 29, 239–257 (2003).
  • [16] Gavrilov, L. & Gavrilova, N. Mortality measurement at advanced ages: A study of the social security administration death master file. \JournalTitleNorth American Actuarial Journal 15, 442–447 (2011).
  • [17] Newman, S. J. Errors as a primary cause of late-life mortality deceleration and plateaus. \JournalTitlePLoS biology 16, e2006776 (2018).
  • [18] Camarda, C. G. The curse of the plateau. measuring confidence in human mortality estimates at extreme ages. \JournalTitleTheoretical Population Biology 144, 24–36 (2022).
  • [19] Missov, T. I. & Vaupel, J. W. Mortality implications of mortality plateaus. \JournalTitlesiam REVIEW 57, 61–70 (2015).
  • [20] Gbari, S., Poulain, M., Dal, L. & Denuit, M. Extreme value analysis of mortality at the oldest ages: a case study based on individual ages at death. \JournalTitleNorth American Actuarial Journal 21, 397–416 (2017).
  • [21] Hanayama, N. & Sibuya, M. Estimating the upper limit of lifetime probability distribution, based on data of japanese centenarians. \JournalTitleJournals of Gerontology Series A: Biomedical Sciences and Medical Sciences 71, 1014–1021 (2016).
  • [22] Einmahl, J. J., Einmahl, J. H. & de Haan, L. Limits to human life span through extreme value theory. \JournalTitleJournal of the American Statistical Association 114, 1075–1080 (2019).
  • [23] Li, J. & Liu, J. A modified extreme value perspective on best-performance life expectancy. \JournalTitleJournal of Population Research 37, 345–375 (2020).
  • [24] Milholland, B. Jeanne Calment, actuarial paradoxography and the limit to human lifespan. \JournalTitleRejuvenation Research 23, 17–18 (2020).
  • [25] Group, G. R. World’s oldest person titleholders since 1955. https://grg.org/Adams/C.HTM (2018). [Online; accessed 14-06-2023].
  • [26] Brillinger, D. R. A biometrics invited paper with discussion: the natural variability of vital rates and associated statistics. \JournalTitleBiometrics 42, 693–734 (1986).
  • [27] Nelder, J. A. & Mead, R. A Simplex Method for Function Minimization. \JournalTitleThe Computer Journal 7, 308–313 (1965).
  • [28] Kiss, C. Modelling the age of the oldest person in the world. https://github.com/csabi0312/modelling-the-age-of-the-oldest-person-in-the-world (2024).
  • [29] Abramowitz, M. & Stegun, I. A. Pocketbook of Mathematical Functions (Verlag Harri Deutsch, Thun-Frankfurt am Main, 1984).

Acknowledgements

We thank Katalin Kovács for some useful advice at an early stage of the project which led to this collaboration. The work of Cs. Kiss and B. Vető was supported by the NKFI (National Research, Development and Innovation Office) grant FK142124. B. Vető is also grateful for the support of the NKFI grant KKP144059 “Fractal geometry and applications” and for the Bolyai Research Scholarship of the Hungarian Academy of Sciences. L. Németh was supported by MaRDI, funded by the Deutsche Forschungsgemeinschaft (DFG), project number 460135501, NFDI 29/1 “MaRDI – Mathematische Forschungsdateninitiative.

Author contributions

B. V. initiated the research. Cs. K. and B. V. derived the model, prepared the scripts and figures, and carried out the estimations. Cs. K., L. N. and B. V. analyzed the results and wrote the manuscript.

Data availability

The titleholder data are freely available at https://grg.org/Adams/C.HTM

Additional information

The authors declare no competing interests.

Tables

E=0E=30E=60CE62704.6801093.2491011κE0.0049870.018760.02085missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression𝐸0𝐸30𝐸60missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝐶𝐸62704.680superscript1093.249superscript1011subscript𝜅𝐸0.0049870.018760.02085\begin{array}[]{|c|ccc|}\hline\cr&E=0&E=30&E=60\\ \hline\cr C_{E}&6270&4.680\cdot 10^{-9}&3.249\cdot 10^{-11}\\ \kappa_{E}&0.004987&0.01876&0.02085\\ \hline\cr\end{array}start_ARRAY start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_E = 0 end_CELL start_CELL italic_E = 30 end_CELL start_CELL italic_E = 60 end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_C start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 6270 end_CELL start_CELL 4.680 ⋅ 10 start_POSTSUPERSCRIPT - 9 end_POSTSUPERSCRIPT end_CELL start_CELL 3.249 ⋅ 10 start_POSTSUPERSCRIPT - 11 end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL italic_κ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.004987 end_CELL start_CELL 0.01876 end_CELL start_CELL 0.02085 end_CELL end_ROW end_ARRAY
Table 1: Birth rate parameters in (22) for different values of the entry age E𝐸Eitalic_E.
E=0E=30E=60αE0.012770.011240.01110KE0.000029510.00059500.01208b0.090.090.09γE0.085960.080610.08026l(αE,KE,0.09,0,γE)119.64117.68117.27missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression𝐸0𝐸30𝐸60missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝛼𝐸0.012770.011240.01110subscript𝐾𝐸0.000029510.00059500.01208𝑏0.090.090.09subscript𝛾𝐸0.085960.080610.08026𝑙subscript𝛼𝐸subscript𝐾𝐸0.090subscript𝛾𝐸119.64117.68117.27\begin{array}[]{|c|ccc|}\hline\cr&E=0&E=30&E=60\\ \hline\cr\alpha_{E}&0.01277&0.01124&0.01110\\ K_{E}&0.00002951&0.0005950&0.01208\\ b&0.09&0.09&0.09\\ \gamma_{E}&0.08596&0.08061&0.08026\\ l(\alpha_{E},K_{E},0.09,0,\gamma_{E})&-119.64&-117.68&-117.27\\ \hline\cr\end{array}start_ARRAY start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_E = 0 end_CELL start_CELL italic_E = 30 end_CELL start_CELL italic_E = 60 end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.01277 end_CELL start_CELL 0.01124 end_CELL start_CELL 0.01110 end_CELL end_ROW start_ROW start_CELL italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.00002951 end_CELL start_CELL 0.0005950 end_CELL start_CELL 0.01208 end_CELL end_ROW start_ROW start_CELL italic_b end_CELL start_CELL 0.09 end_CELL start_CELL 0.09 end_CELL start_CELL 0.09 end_CELL end_ROW start_ROW start_CELL italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.08596 end_CELL start_CELL 0.08061 end_CELL start_CELL 0.08026 end_CELL end_ROW start_ROW start_CELL italic_l ( italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT , 0.09 , 0 , italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ) end_CELL start_CELL - 119.64 end_CELL start_CELL - 117.68 end_CELL start_CELL - 117.27 end_CELL end_ROW end_ARRAY

Table 2: The optimal parameters obtained for b=0.09𝑏0.09b=0.09italic_b = 0.09 and for various values of the entry age E=0,30,60𝐸03060E=0,30,60italic_E = 0 , 30 , 60.
E=0E=30E=60αE0.015610.013760.01352KE3.7281060.00014950.005901b0.110.110.11γE0.11600.11170.1110l119.48117.51117.09E=0E=30E=60αE0.018450.016280.01596KE4.5301073.6131050.002784b0.130.130.13γE0.14410.140650.1399l119.37117.40116.98missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression𝐸0𝐸30𝐸60missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝛼𝐸0.015610.013760.01352subscript𝐾𝐸3.728superscript1060.00014950.005901𝑏0.110.110.11subscript𝛾𝐸0.11600.11170.1110𝑙119.48117.51117.09missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression𝐸0𝐸30𝐸60missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝛼𝐸0.018450.016280.01596subscript𝐾𝐸4.530superscript1073.613superscript1050.002784𝑏0.130.130.13subscript𝛾𝐸0.14410.140650.1399𝑙119.37117.40116.98\begin{array}[]{|c|ccc|}\hline\cr&E=0&E=30&E=60\\ \hline\cr\alpha_{E}&0.01561&0.01376&0.01352\\ K_{E}&3.728\cdot 10^{-6}&0.0001495&0.005901\\ b&0.11&0.11&0.11\\ \gamma_{E}&0.1160&0.1117&0.1110\\ l&-119.48&-117.51&-117.09\\ \hline\cr\end{array}\quad\begin{array}[]{|c|ccc|}\hline\cr&E=0&E=30&E=60\\ \hline\cr\alpha_{E}&0.01845&0.01628&0.01596\\ K_{E}&4.530\cdot 10^{-7}&3.613\cdot 10^{-5}&0.002784\\ b&0.13&0.13&0.13\\ \gamma_{E}&0.1441&0.14065&0.1399\\ l&-119.37&-117.40&-116.98\\ \hline\cr\end{array}start_ARRAY start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_E = 0 end_CELL start_CELL italic_E = 30 end_CELL start_CELL italic_E = 60 end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.01561 end_CELL start_CELL 0.01376 end_CELL start_CELL 0.01352 end_CELL end_ROW start_ROW start_CELL italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 3.728 ⋅ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT end_CELL start_CELL 0.0001495 end_CELL start_CELL 0.005901 end_CELL end_ROW start_ROW start_CELL italic_b end_CELL start_CELL 0.11 end_CELL start_CELL 0.11 end_CELL start_CELL 0.11 end_CELL end_ROW start_ROW start_CELL italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.1160 end_CELL start_CELL 0.1117 end_CELL start_CELL 0.1110 end_CELL end_ROW start_ROW start_CELL italic_l end_CELL start_CELL - 119.48 end_CELL start_CELL - 117.51 end_CELL start_CELL - 117.09 end_CELL end_ROW end_ARRAY start_ARRAY start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_E = 0 end_CELL start_CELL italic_E = 30 end_CELL start_CELL italic_E = 60 end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.01845 end_CELL start_CELL 0.01628 end_CELL start_CELL 0.01596 end_CELL end_ROW start_ROW start_CELL italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 4.530 ⋅ 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT end_CELL start_CELL 3.613 ⋅ 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT end_CELL start_CELL 0.002784 end_CELL end_ROW start_ROW start_CELL italic_b end_CELL start_CELL 0.13 end_CELL start_CELL 0.13 end_CELL start_CELL 0.13 end_CELL end_ROW start_ROW start_CELL italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.1441 end_CELL start_CELL 0.14065 end_CELL start_CELL 0.1399 end_CELL end_ROW start_ROW start_CELL italic_l end_CELL start_CELL - 119.37 end_CELL start_CELL - 117.40 end_CELL start_CELL - 116.98 end_CELL end_ROW end_ARRAY
E=0E=30E=60αE0.021290.018800.01842KE5.3621088.5001060.001281b0.150.150.15γE0.17080.16810.1674l119.31117.34116.91missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression𝐸0𝐸30𝐸60missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝛼𝐸0.021290.018800.01842subscript𝐾𝐸5.362superscript1088.500superscript1060.001281𝑏0.150.150.15subscript𝛾𝐸0.17080.16810.1674𝑙119.31117.34116.91\begin{array}[]{|c|ccc|}\hline\cr&E=0&E=30&E=60\\ \hline\cr\alpha_{E}&0.02129&0.01880&0.01842\\ K_{E}&5.362\cdot 10^{-8}&8.500\cdot 10^{-6}&0.001281\\ b&0.15&0.15&0.15\\ \gamma_{E}&0.1708&0.1681&0.1674\\ l&-119.31&-117.34&-116.91\\ \hline\cr\end{array}start_ARRAY start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_E = 0 end_CELL start_CELL italic_E = 30 end_CELL start_CELL italic_E = 60 end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_α start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.02129 end_CELL start_CELL 0.01880 end_CELL start_CELL 0.01842 end_CELL end_ROW start_ROW start_CELL italic_K start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 5.362 ⋅ 10 start_POSTSUPERSCRIPT - 8 end_POSTSUPERSCRIPT end_CELL start_CELL 8.500 ⋅ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT end_CELL start_CELL 0.001281 end_CELL end_ROW start_ROW start_CELL italic_b end_CELL start_CELL 0.15 end_CELL start_CELL 0.15 end_CELL start_CELL 0.15 end_CELL end_ROW start_ROW start_CELL italic_γ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_CELL start_CELL 0.1708 end_CELL start_CELL 0.1681 end_CELL start_CELL 0.1674 end_CELL end_ROW start_ROW start_CELL italic_l end_CELL start_CELL - 119.31 end_CELL start_CELL - 117.34 end_CELL start_CELL - 116.91 end_CELL end_ROW end_ARRAY
Table 3: The optimal parameters for b=0.11𝑏0.11b=0.11italic_b = 0.11, b=0.13𝑏0.13b=0.13italic_b = 0.13 and b=0.15𝑏0.15b=0.15italic_b = 0.15.