# Derivation of Closed-Form Gompertz-Makeham Life Expectancy

I am trying to replicate the derivation of the life expectancy formulas from Missov & Lenart (2013). They propose that the following:

Theorem 2. Life expectancy of the Gompertz–Makeham mortality distribution at age $$x$$ is: $$e_{GM}(x) = \frac{1}{b} e^\frac{a}{b} \left( \frac{a}{b} \right)^\frac{c}{b} \Gamma \left( -\frac{c}{b} , \frac{a}{b}e^{bx}\right) \qquad a,b > 0, \; c\geq 0$$ where $$\Gamma(s,z) = \int_z^\infty t^{s-1} e^{-t} dt$$ denotes the upper incomplete gamma function

Their proof goes as follows:

Proof. The individual force of mortality and the survival function for a Gompertz–Makeham homogeneous population is given by (1) and (2) for z = 1, i.e $$\mu_{GM}(x) = a ^{bx} + c \qquad \qquad S_{GM}(x) = \exp \left[ -\frac{a}{b} (e^{bx} - 1) -cx \right]$$ Remaining life expectancy at age $$x$$ can be calculated by $$e_{GM}(x) = \int_x^\infty S_{GM}(t, y) dt \tag{B1}$$

So far I can follow, except for the added $$y$$ argument to the survival function. Nevertheless, they continue:

Substituting $$s = e^{bt}$$ and denoting $$a = \frac{a}{b}$$ reduces (B.1) to $$e_{GM}(x) = \frac{e^a}{b} \int_{e^{bx}}^\infty e^{-as} q^{-\frac{c}{b}-1} dq$$

Here I have a problem: we seem to subsitute using $$s$$, but a variable $$q$$ appears. My intuition is that thefirst step here is to plug in the survival function as $$\exp \left[ -\frac{a}{b} (e^{bt} - 1) -ct \right]$$ and substitute for $$a$$, which allows us to draw $$e^a$$ out of the integral. But the we need to substitute for $$t$$, so what happens to $$ct$$? And how is $$q$$ defined?

To be complete, they the proof continues:

Finally we substitute $$u = as$$ and get $$e_{GM}(x) = \frac{e^a}{b}[a]^\frac{c}{b} \int_{ae^{bx}}^\infty e^{-u} u^{-\frac{c}{b}-1} du = \frac{1}{b} e^\frac{a}{b} \left( \frac{a}{b} \right)^\frac{c}{b} \Gamma \left( -\frac{c}{b} , \frac{a}{b}e^{bx}\right)$$

As happened previously, formulating the question and thinking about it helped to find a solution. Here it is:

We start from \begin{aligned} e(x) &= \int_x^\infty \exp \left[ -ct - \frac{a}{b} (e^{bt}-1) \right] \mathsf{d}t \\ &= \int_x^\infty \exp \left[ -ct - \frac{a}{b} e^{bt} + \frac{a}{b} \right] \mathsf{d}t \\ &= \exp\left( \frac{a}{b} \right) \int_x^\infty \exp \left[ -ct - \frac{a}{b} e^{bt} \right] \mathsf{d}t \end{aligned}

Now, define $$q(t) = e^{bt}$$, so that we can apply integration by substitution:

$$\frac{\mathsf{d} q}{\mathsf{d} t} = q'(t) = b e^{bt} \qquad \Leftrightarrow \qquad \mathsf{d} q = \mathsf{d} t \cdot q'(t)$$

This might be somewhat loose on notation, but it works (I'm not a mathematician). Nevertheless, we can continue: \begin{aligned} \exp\left( \frac{a}{b} \right) \int_x^\infty \exp \left[ -ct - \frac{a}{b} e^{bt} \right] \mathsf{d}t &= \exp\left( \frac{a}{b} \right) \int_x^\infty \exp \left[ -ct - \frac{a}{b} e^{bt} \right] \frac{q'(t)}{q'(t)}\mathsf{d}t \\ &= \exp\left( \frac{a}{b} \right) \int_x^\infty e^{-ct} \exp \left[-\frac{a}{b} e^{bt} \right] \frac{1}{b}e^{-bt} \;\; \mathsf{d} t \cdot q'(t) \\ &= \frac{1}{b}\exp\left( \frac{a}{b} \right) \int_x^\infty \exp \left[-\frac{a}{b} e^{bt} \right] \left(e^{bt} \right)^{-\frac{c}{b} - 1} \;\; \mathsf{d} t \cdot q'(t) \end{aligned}

Where we used the fact that $$\left(e^{bt} \right)^{-\frac{c}{b} - 1} = e^{-bt - ct}$$. We then get to the first step in the proof by substituting in $$A = a/b$$ (less confusing with a capital A) and the definition of $$q(t)$$, while also setting the lower bound to $$q(x)$$: \begin{aligned} \frac{e^A }{b} \int_{e^{bx}}^\infty e^{-Aq} q^{-\frac{c}{b} - 1} \;\; \mathsf{d} q \end{aligned} So $$s$$ seems to have been $$q$$ all along. We can then continue, by defining $$u(q) = Aq$$ so that: $$\frac{\mathsf{d} u}{\mathsf{d} q} = u'(q) = A \qquad \Leftrightarrow \qquad \mathsf{d} u = \mathsf{d} q \cdot u'(q)$$

and perform a second round of integration by substitution.

\begin{aligned} \frac{e^A }{b} \int_{e^{bx}}^\infty e^{-Aq} q^{-\frac{c}{b} - 1} \;\; \mathsf{d} q &= \frac{e^A }{b} \int_{e^{bx}}^\infty e^{-Aq} q^{-\frac{c}{b} - 1} \frac{A^{\frac{c}{b}}}{A^{\frac{c}{b}}}\;\; \frac{u'(q)}{u'(q)}\mathsf{d} q \\ &= \frac{e^A }{b} A^{\frac{c}{b}} \int_{e^{bx}}^\infty e^{-Aq} q^{-\frac{c}{b} - 1} A^{-\frac{c}{b}} A^{-1}\;\; u'(q)\cdot\mathsf{d} q \\ &= \frac{e^A }{b} A^{\frac{c}{b}} \int_{e^{bx}}^\infty e^{-Aq} q^{-\frac{c}{b} - 1} A^{-\frac{c}{b} -1}\;\; u'(q)\cdot\mathsf{d} q \\ &= \frac{e^A }{b} A^{\frac{c}{b}} \int_{e^{bx}}^\infty e^{-Aq} u^{-\frac{c}{b} - 1} \;\; u'(q) \cdot\mathsf{d} q \\ &= \frac{e^A }{b} A^{\frac{c}{b}} \int_{Ae^{bx}}^\infty e^{-u} u^{-\frac{c}{b} - 1} \;\; \mathsf{d} u \end{aligned} The last step to the proof is then simply to apply the definition of the upper-incomplete gamma function.