Modelling the age distribution of longevity leaders

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.


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) 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.
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.5.This also means that with high probability her age record will already be broken by that time.
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.000286 for Calment and 0.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.195 in 1955 and it is 1.188 in 2019.The empirical value of the reign length is 1.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.

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.1-1.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.3-1.5.We explain the role of the entry age parameter in Subsection 1.6.

Model description and two-dimensional representation
The model is formally defined as follows.Let λ (t) be the birth rate parameter which depends on time and let F t and f t 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) and that the lifespan of an individual born at time t is given by F t so that lifespans are independent for different individuals.
Let Y t denote the age of the oldest person in the population at time t.The process (Y t : t ∈ R) is Markovian.The Markov property holds because at any time t the history of the process (Y s : s ≤ t) provides information about the lifetime of individuals born before the current record holder while any transition of (Y s : s ≥ t) depends only on the lifetime of the current record holder and of those born after them.
The evolution of the Markov process Y t is the following.It has a deterministic linear growth with slope 1 due to the ageing of the current record holder.This happens until the death of the record holder.Additionally, given that Y t− = lim s↑t Y s = y for some t with y > 0, the process has a downward jump at time t at rate f t (y)/(1 − F t (y)) which is the hazard rate of the distribution F t at y.This corresponds to the possibility that the record holder dies at time t which happens at rate f t (y)/(1 − F t (y)).The conditional distribution of the jump is given by for all x > 0. The jump distribution in (1) has an absolutely continuous part supported on [0, y] with density 3/14 and a point mass at 0 with probability As we shall see in the relevant parameter regime the probability a y,t of the point mass at 0 is negligible.The transition formula in (1) can be proven using the description below.
We introduce a two-dimensional representation of the process Y t as follows.Let Λ = {(t i , x i ) : i ∈ I} be a marked Poisson process in R × R + where {t i : i ∈ I} forms a Poisson point process on R with intensity λ (t) and x i ≥ 0 is sampled independently for each i ∈ I according to the distribution F t i .The point (t i , x i ) represents an individual born at time t i with lifespan x i for all i ∈ I, that is, the individual i is alive in the time interval [t i ,t i + x i ) and their age at time Hence the marked Poisson process Λ contains all relevant information about the age statistics of the population at any time.In particular the age of the oldest person Y t can be expressed in terms of Λ as where the indicator 1 {t∈[t i ,t i +x i )} is 1 exactly if the ith person is alive at time 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 at age y the event {Y t < x} means that nobody with age between x and y can be alive at time t.This event can be equivalently characterized in terms of the Poisson process of birth at rate λ (•) thinned by the probability that the person is still alive at time t.Indeed the event {Y t < x} can be expressed as a Poisson process of intensity at time t − u given by λ . This probability appears exactly on the right-hand side of (1).In other words, for any u ∈ [x, y] people are born at time t − u at rate λ (t − u).On the other hand, the probability for a person born at time t − u to be alive at time t (that is, at age u) is 1 − F t−u (u).See Figure 3 for illustration.

Exact distribution of the oldest person's age process
We assume that all birth events are already sampled on (−∞,t] together with the corresponding lifespans.Then the distribution of Y t can be computed explicitly for all t ∈ R using the two-dimensional representation.For all t ∈ R the density for all x > 0 and the point mass at 0 characterize the distribution of Y t which can be seen as follows.We mention that the point mass at 0 is negligible in the application.Similarly to the proof of the transition formula in (1) the event {Y t < x} for any x > 0 is the same as the event that nobody with age at least x is alive at time t.We express this event in terms of the Poisson process of birth at rate λ (•) thinned by the probability that the person is still alive at time t.The event {Y t < x} means that a Poisson process of intensity at time t − u given by λ In other words for any u ≥ x individuals are born at time t − u at rate λ (t − u).A person born at time t − u is alive at time t at age u with probability 1 − F t−u (u).( 5)-( 6) follow by differentiation in (7) and by taking the x → 0 limit.See Figure 4 for illustration.

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 for all t.The lifespan of individuals are independent and identically distributed with a fixed density f = f t and cumulative distribution function F = F t for all t which does not depend on time.
In the homogeneous model the jump distribution given in ( 2)-(3) simplifies to The distribution of Y t does not depend on time in this case hence it is a stationary distribution as well.The formulas for the density of Y t and point mass at 0 reduce in the homogeneous case to where the integral ∞ 0 (1 − F(u)) du is equal to the expected lifespan.The equilibrium condition for the homogeneous density h can be written as

5/14
After differentiation and using the fact that one can derive from (10) the second order differential equation The point mass m at 0 satisfies The formula can be seen as follows.To have a record holder who dies at age x there has to be a person who has lifespan x which gives the factor f (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) a probability density function.
The density of Z n can also be characterized by the following description.It satisfies the integral equation which comes from the possible transitions of the peak process as follows.If the previous record holder had a total lifetime w ∈ [0, ∞) then at the death the process Y t jumps down to some value y at rate j w (y) or to 0 with probability a w .The density of the age at which a person dies who becomes a record holder at age y is f (x)/(1 − F(y)).From the integral equation in (15) one can derive the second order differential equation for the function g(x) = z(s)/ f (x) given by which is satisfied by g ) du) in accordance with (14).

Reign length distribution
In the homogeneous model, let W n denote the reign length of the 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 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 6/14 based on the decomposition with respect to the previous value of the peaks process Z n .The integral w 0 f (z)λ e −λ (w−z) dz is the density of the convolution of the density f with an independent exponential distribution of parameter λ .On the right-hand side of (18), one can use the definitions of the density z given by ( 14) and the homogeneous jump distribution j x and a x 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 which follows by integration by parts.Note also that the density is not equal to the remaining reign length of Y t under the stationary distribution because it would involve the integral ∞ 0 h(y) f (y + w)/(1 − F(y)) dy in place of the first term in the numerator on the right-hand side of ( 17).
because the new birth process of rate λ E (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 with age E becomes the remaining lifetime distribution at age E.
The modified cumulative distribution function and density are given by

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.

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 and we fit the model parameters for all three values of 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 we assume that the birth rate at age E is given by where the numerical values of the parameters C E and κ E are obtained by linear regression of the logarithmic data of newborns, people at age 30 and 60 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) distribution with cumulative distribution function and density 7/14 for x ≥ 0 where a, b, γ are positive parameters.We mention that the gamma-Gompertz-Makeham (ΓGM) distribution differs from the ΓG distribution by the presence of a non-negative extrinsic mortality parameter c which appears as an additive term in the force of mortality.See (29) for the definition of the Γ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 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 lifespan distribution, parameters b = b E , the rate of aging and γ = γ E , the magnitude of heterogeneity are constants over time and that they only depend on the value of the entry age parameter E. The parameter a, the initial level of mortality at the entry age for individuals born at time t, depends on time given by the exponentially decreasing function where the exponent α E and the constant K E only depends on the entry age E. The reason for subtracting 2000 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, we assume that the birth rate λ E (t) is given by ( 22) and we fit the gamma-Gompetz distribution with parameters b E , γ E and a E (t) given by ( 24) for the modified distribution function F E t (x) and density f E t (x) in ( 21).This means that we search for the best fitting values of the parameters α E , K E , b E , γ E which results in an approximation of the remaining lifetime distribution at the age E.

Likelihood calculations
The aim of the maximum likelihood method is to give an estimate to the parameters α E , K E , b E and γ E for 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 (t i , y i , z i ) for i = 1, . . ., n where t i is the ith time in the sample when the oldest person dies at age y i and the new record holder has age z i at time t i .Then the data has to satisfy the consistency relation t i − z i = t i+1 − y i+1 since the two sides express the date of birth of the same person.
In the model with entry age E, the likelihood of the ith data point (t i , y i , z i ) given the previous data point is equal to for all i = 2, 3, . . ., n except for i = 1 in which case the 1 − F E y 1 −t 1 +E (z 0 − E) factor in the denominator is missing.In ( 25) above we use the transition probabilities of the model with entry age E given by as a generalization of (2).The explanation of the left-hand side of ( 25) is that the person died at time t i at age y i had age E at time t i − y i + E. The previous data point ensures that this person has already reached age z i−1 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 with u ∈ [x, y] at time t had age E at time t − u + E.
Note that when computing the likelihood of the full data by multiplying the right-hand side of (25) for different values of i the consistency relation of the data implies that the factor 1 − F t i −z i +E (z i − E) of the ith term cancels with the factor 1 − F t i+1 −y i+1 +E (z i − E) coming from the (i + 1)st term.Hence the log-likelihood of the full sample is given by l(α, K, b, γ) where we suppress the dependence of the parameters α, K, b, γ on the entry age.Note that the last two terms do not depend on the parameters α, K, b, γ hence we can omit these terms in the maximization of the log-likelihood.

Likelihood maximization
We implemented the calculation of the log-likelihood function l(α, K, b, γ) 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 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, γ) 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 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, γ using the gamma-Gompretz-Makeham model was very similar the case of four parameters α, K, c, γ in the gamma-Gompertz model.Running the optimization in the full set of parameters (α, K, b, c, γ in the gamma-Gompertz-Makeham model or α, K, b, γ in the gamma-Gompertz model), it turned out that after a few rounds the parameter K started to decrease dramatically and reached values below 10 −10 .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.
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, γ 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 R 4 + 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 model is b which is the exponent in the time dependence of the mortality rate.By setting the rate of aging b = 0.09 the algorithm gives the optimal triple α, K, γ 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   3.
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.

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 is given by the density h t (x) in ( 5) and by the point mass m t at 0 in (6).For the numerical computations, we ignore the point mass m t 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 is that parameter a of the gamma-Gompertz-Makeham distribution function F t−u in the exponent of (5) also depends on the integration variable u.
In our approximation we fix the value of the parameter a of the ΓG distribution in h t (x) in (5) to a value which is equal to a 0 (t − d) in (24) with some delay d.The delay d is chosen so that the mean age of the oldest person computed using a 0 (t − d) as parameter a for all times in the ΓG distribution function in ( 5) is equal to the same value d.For a given t, this value of d can be obtained as the fixed point of the contraction map 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. Hence in (28), the function λ is given in ( 22) with E = 0 and F t−d is the ΓG distribution function with parameters given by the E = 0 values in Table 2 and with a = a 0 (t − 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, γ given in Tables 1 and 2 for entry age E = 0.
In our model, the distribution function of the age of the world's oldest person is given in (7) where the Γ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.000286 for Calment and it is 0.0116 for Knauss.

10/14
The backtesting mentioned in the Results section is performed as follows.We estimated the best parameter values with entry age 0 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.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.195 in 1955 and it is 1.188 in 2019.The empirical value of the reign length is 1.008 computed from the data by dividing the total length of the time interval between 1955 and 2019 by the number of record holders.

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) distribution are given by for x ≥ 0 where a, b, c, γ 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 (1 − F ΓGM a,b,c,γ (u)) du.In the homogeneous model, the integral of the survival function appears in the density of the distribution of Y t 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 where we applied a change of variables y = e x−u in the second equality above and we applied the hypergeometric identity (31) in the last equality with a = α, b = α + β , c = 1 + α + β , z = −e −x /δ 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 α + β .Note that the condition Re(c) > Re(b) > 0 for (31) to hold is satisfied by our assumption Re(α + β ) > 0 which also makes the integrals in (32) convergent.

Figure 1 .
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.

Figure 3 .
Figure 3. Two-dimensional representation of the event {Y t < x} conditionally given {Y t− = y,Y t < y}: the marked Poisson process does not have any point in the grey region, see (1).

Figure 4 .
Figure 4. Two-dimensional representation of the event {Y t < x}: the marked Poisson process does not have any point in the grey region, see (7).

Figure 5 .
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.

) 1 .4 The peaks process
In the homogeneous model the sequence of peaks in Y t forms a discrete time Markov chain.By peak we mean a local maximum of Y t with value being equal to the lifespan of the last record holder.Each time the oldest person dies the process Y t has a peak with a downward jump following it.LetZ n denote the age of record holders at which they die which are the values of the peaks of the process Y t .The sequence Z n 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 Z n is given by

1.6 Entry age parameter
Next we introduce another parameter which we call the entry age and we denote it by E. As opposed to our original model we consider individuals as being born at age 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 E ≥ 0 and we can fit the model parameters with different values of the entry age.For any value E ≥ 0 of the entry age we denote by λ E (t) the rate at which people reach the age E at time t, that is, Table 2 for individuals born in 2000 corresponding to the entry age

Table 2 .
The optimal parameters obtained for b = 0.09 and for various values of the entry age E = 0, 30, 60.