Some points for integration

The survivor function from age \(x\) to age \(x+t\), denoted \({}_tp_x\) by actuaries, is a useful tool in mortality work. As mentioned in one of our earliest blogs, a basic feature is that the expected time lived is the area under the survival curve, i.e. the integral of \({}_tp_x\). This is easy to express in visual terms, but it often requires numerical integration if there is no closed-form expression for the integral of the survival curve. In this article we look at some of the options available to actuaries who need to integrate numerically.

A general result is that the survivor function has the following form:

\[{}_tp_x = e^{-H_x(t)}\]

where \(H_x(t)\) is the integrated hazard function:

\[H_x(t) = \int_0^t \mu_{x+s}ds\]

and \(\mu_x\) is the force of mortality at age \(x\). Richards (2014) tabulates closed-form expressions for \(\mu_x\) and \(H_x(t)\) for a wide variety of mortality models. However, to calculate the expected time lived over an age interval \([x, u)\), we need to perform a second integral:

\[\int_0^{u-x} {}_tp_x dt\]

If we let \(u\to\infty\) then we would get the complete life expectancy from age \(x\). Now, how could we go about approximating this integral between the ages of 65 and 100 (say)? One approach is simple brute force: we approximate the survivor function by a series of steps of width \(h\) and sum the area of the rectangles. Mathematically this would be:

\[\int_0^{u-x}{}_tp_x \approx h\sum_{i=1}^{nh} {}_tp_x\qquad(1)\]

where we have split the domain of integration into \(n\) equally sized intervals of width \(h\). If we can evaluate \({}_tp_x\) for any value of \(t\), we can select an arbitrarily small value of \(h\) to get an approximation of the integral using equation (1). This is illustrated in Table 1 for the Gompertz law of mortality using various values of \(h\).

Table 1. Expected time lived between ages 65 and 100 under Gompertz law of mortality, \(\mu_x=e^{\alpha+\beta x}\), with \(\alpha=-12\) and \(\beta=0.12\).

The brute-force approach in Table 1 looks like it is converging to 15.33. However, this is at the cost of a fair amount of computing resource: for \(h=0.0001\) the survival function has had to be computed 350,000 times. A smarter approach is to use something like adaptive quadrature, which requires far fewer evaluations of the survival function. R offers this in the \(\tt integrate()\) function for single-variable calculations: for our Gompertz example this function returns the value 15.3315 with an absolute error less than 8.1e-8.

The brute-force and adaptive-quadrature methods both rely on being able to calculate \({}_tp_x\) at arbitrary points. However, in actuarial work it is common to have to calculate a life expectancy — or an annuity factor — where \({}_tp_x\) is only available for integer values of \(t\). How can we accurately integrate under a curve where the values of the curve are only available at fixed points? One approach is to use some results from the history of mathematics:

For two points separated by one year we use the Trapezoidal Rule:

\[\int_a^{a+1} f(x) \approx \frac{1}{2}\left[f(a)+f(a+1)\right]\]

For three points spaced one year apart we use Simpson's Rule:

\[\int_a^{a+2} f(x) \approx \frac{1}{3}\left[f(a) + 4f\left(a+1\right) + f(a+2)\right]\]

For four points spaced one year apart we use Simpson's 3/8 Rule:

\[\int_a^{a+3} f(x) \approx \frac{3}{8}\left[f(a) + 3f(a+1) + 3f(a+2) + f(a+3)\right]\]

To integrate over \(n\) points at yearly intervals we would first apply Simpson's 3/8 Rule as many times as possible, then Simpson's Rule and then the Trapezoidal Rule for any remaining points at the highest ages. How accurate would this simple approach be? If we apply it to our Gompertz example with \({}_tp_{65}\) evaluated for \(t = 0, 1, 2,\ldots,35\) we get 15.3315, the same value returned by R's \(\tt integrate()\) function. This shows that a remarkable degree of accuracy can be obtained just from using values at yearly-spaced intervals.