Derivation of Closed-Form Gompertz-Makeham Life Expectancy

Question

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)$$

Answer 1

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.