This article values a continuous whole-life insurance benefit under a Gompertz–Makeham mortality law. For a life aged 40, a death benefit of \(\$500{,}000\), and a 4% force of interest, the expected present value is approximately \(\$94{,}035\). Its standard deviation is approximately \(\$76{,}584\), reflecting substantial uncertainty in the timing of death. The calculation is verified using two numerically independent integral representations, then extended to interest-rate and mortality-parameter sensitivity tests.
Consider a continuous whole-life insurance policy issued to a person aged 40. The policy pays a death benefit of \(\$500{,}000\) immediately on death.
The force of mortality is Gompertz–Makeham:
\[ \mu(x)=A+Bc^x, \]
where
\[ A=0.0005,\qquad B=0.000075,\qquad c=1.075. \]
The force of interest is constant:
\[ \delta=0.04. \]
Key results
These are model-based values for the death benefit only; they are not a commercial premium quotation.
Let \(T_{40}\) denote the future lifetime of the insured. The probability that the insured survives at least another \(t\) years is
\[ {}_tp_{40}=\Pr(T_{40}>t) =\exp\left[-\int_0^t\mu(40+s)\,ds\right]. \]
Under the specified mortality law,
\[ \begin{aligned} \int_0^t\mu(40+s)\,ds &=\int_0^t\left(A+Bc^{40+s}\right)ds\\ &=At+Bc^{40}\int_0^t c^s\,ds\\ &=At+\frac{Bc^{40}}{\log c}(c^t-1). \end{aligned} \]
Therefore,
\[ \boxed{ {}_tp_{40} = \exp\left[ -At-\frac{Bc^{40}}{\log c}(c^t-1) \right]. } \]
For the given parameters,
\[ {}_tp_{40} = \exp\left[ -0.0005t -\frac{0.000075(1.075)^{40}}{\log(1.075)} \left(1.075^t-1\right) \right]. \]
The future-lifetime density is
\[ \boxed{ f_{T_{40}}(t) =\mu(40+t){}_tp_{40} =\left(A+Bc^{40+t}\right){}_tp_{40}. } \]
mu <- function(t, A_par = A, B_par = B, c_par = c_base) {
A_par + B_par * c_par^(x + t)
}
survival <- function(t, A_par = A, B_par = B, c_par = c_base) {
exp(
-A_par * t -
(B_par * c_par^x / log(c_par)) * (c_par^t - 1)
)
}
density <- function(t, A_par = A, B_par = B, c_par = c_base) {
mu(t, A_par, B_par, c_par) *
survival(t, A_par, B_par, c_par)
}
mu_40 <- mu(0)
gompertz_integral_coefficient <- B * c_base^x / log(c_base)
mu_40## [1] 0.001853318
## [1] 0.01871274
The force of mortality at age 40 is approximately 0.00185332, or about 0.185% per year on an instantaneous basis.
The complete expectation of future lifetime can be found from
\[ E[T_{40}]=\int_0^\infty {}_tp_{40}\,dt. \]
The median solves \({}_tp_{40}=0.5\). The mode maximizes the future-lifetime density.
expected_lifetime <- integrate(
survival,
lower = 0,
upper = 150,
rel.tol = 1e-11,
subdivisions = 2000
)$value
quantile_lifetime <- function(probability) {
uniroot(
function(t) survival(t) - (1 - probability),
interval = c(0, 150),
tol = 1e-11
)$root
}
median_lifetime <- quantile_lifetime(0.50)
mode_lifetime <- optimize(
function(t) -density(t),
interval = c(0, 150),
tol = 1e-11
)$minimum
survival_60 <- survival(60)
lifetime_summary <- data.frame(
Quantity = c(
"Force of mortality at age 40",
"Expected future lifetime",
"Median future lifetime",
"Modal future lifetime",
"Probability of survival for 60 years",
"Probability of death within 60 years"
),
Value = c(
mu_40,
expected_lifetime,
median_lifetime,
mode_lifetime,
survival_60,
1 - survival_60
),
Units = c(
"per year",
"years",
"years",
"years",
"probability",
"probability"
)
)
knitr::kable(
lifetime_summary,
digits = 6,
caption = "Selected lifetime results for the insured."
)| Quantity | Value | Units |
|---|---|---|
| Force of mortality at age 40 | 0.001853 | per year |
| Expected future lifetime | 47.546354 | years |
| Median future lifetime | 49.820566 | years |
| Modal future lifetime | 54.819432 | years |
| Probability of survival for 60 years | 0.235604 | probability |
| Probability of death within 60 years | 0.764396 | probability |
Because the benefit is paid immediately at death, its random present value is
\[ Z=500{,}000e^{-\delta T_{40}}. \]
For a unit benefit, define
\[ \overline A_{40}(\delta) =E[e^{-\delta T_{40}}] =\int_0^\infty e^{-\delta t}f_{T_{40}}(t)\,dt. \]
Thus,
\[ E[Z] =500{,}000\,\overline A_{40}(0.04). \]
discounted_moment_density <- function(
discount_force,
A_par = A,
B_par = B,
c_par = c_base) {
if (discount_force == 0) {
return(1)
}
integrand <- function(t) {
exp(-discount_force * t) *
density(t, A_par, B_par, c_par)
}
# For c > 1, survival beyond 150 years is numerically zero. When a
# sensitivity scenario gives c <= 1, use a longer range for the Makeham tail.
integration_upper <- if (c_par > 1) 150 else 5000
integrate(
integrand,
lower = 0,
upper = integration_upper,
rel.tol = 1e-11,
subdivisions = 3000
)$value
}
unit_epv <- discounted_moment_density(delta)
benefit_epv <- benefit * unit_epv
unit_epv## [1] 0.1880691
## [1] 94034.56
The expected present value is therefore
\[ \boxed{ E[Z] =500{,}000(0.188069116556) =\text{\$94,034.56}. } \]
This amount is the net single premium for the death benefit alone. It does not include expenses, profit, capital costs, or margins for adverse deviation.
The calculation above integrates the discounted future-lifetime density. A second representation follows from integration by parts:
\[ \boxed{ \overline A_{40}(\delta) =1-\delta\int_0^\infty e^{-\delta t}{}_tp_{40}\,dt. } \]
This second method integrates the survival function rather than the density.
discounted_moment_survival <- function(
discount_force,
A_par = A,
B_par = B,
c_par = c_base) {
if (discount_force == 0) {
return(1)
}
integrand <- function(t) {
exp(-discount_force * t) *
survival(t, A_par, B_par, c_par)
}
integration_upper <- if (c_par > 1) 150 else 5000
annuity_value <- integrate(
integrand,
lower = 0,
upper = integration_upper,
rel.tol = 1e-11,
subdivisions = 3000
)$value
1 - discount_force * annuity_value
}
unit_epv_check <- discounted_moment_survival(delta)
verification_epv <- data.frame(
Method = c(
"Discounted density integral",
"Integration-by-parts survival integral"
),
Result = c(unit_epv, unit_epv_check),
Difference_from_first = c(0, unit_epv_check - unit_epv)
)
knitr::kable(
verification_epv,
digits = 15,
caption = "Independent verification of the unit-benefit EPV."
)| Method | Result | Difference_from_first |
|---|---|---|
| Discounted density integral | 0.1880691 | 0 |
| Integration-by-parts survival integral | 0.1880691 | 0 |
Since
\[ Z^2=(500{,}000)^2e^{-2\delta T_{40}}, \]
the second moment is
\[ E[Z^2] =(500{,}000)^2 \int_0^\infty e^{-2\delta t}f_{T_{40}}(t)\,dt. \]
It follows that
\[ \operatorname{Var}(Z) =E[Z^2]-E[Z]^2. \]
unit_second_moment <- discounted_moment_density(2 * delta)
unit_second_moment_check <- discounted_moment_survival(2 * delta)
benefit_second_moment <- benefit^2 * unit_second_moment
benefit_variance <- benefit^2 * (unit_second_moment - unit_epv^2)
benefit_sd <- sqrt(benefit_variance)
risk_results <- data.frame(
Quantity = c(
"Unit-benefit first moment",
"Unit-benefit second moment",
"Benefit EPV",
"Benefit second moment",
"Benefit variance",
"Benefit standard deviation"
),
Value = c(
unit_epv,
unit_second_moment,
benefit_epv,
benefit_second_moment,
benefit_variance,
benefit_sd
),
Units = c(
"dimensionless",
"dimensionless",
"dollars",
"dollars squared",
"dollars squared",
"dollars"
)
)
knitr::kable(
risk_results,
digits = 6,
caption = "Present-value moments and risk measures."
)| Quantity | Value | Units |
|---|---|---|
| Unit-benefit first moment | 0.188069 | dimensionless |
| Unit-benefit second moment | 0.058830 | dimensionless |
| Benefit EPV | 94034.558278 | dollars |
| Benefit second moment | 14707549736.738316 | dollars squared |
| Benefit variance | 5865051586.227781 | dollars squared |
| Benefit standard deviation | 76583.624792 | dollars |
The principal monetary results are
\[ \boxed{ E[Z]=\text{\$94,034.56} } \]
\[ \boxed{ \operatorname{Var}(Z) =5.8651e+09 \ \text{dollars}^2 } \]
\[ \boxed{ \operatorname{SD}(Z) =\text{\$76,583.62} } \]
The standard deviation is large because the time of payment is highly uncertain. An early death produces a present value close to the full \(\$500{,}000\), whereas a payment many decades later is heavily discounted.
The second-moment verification is:
verification_second_moment <- data.frame(
Method = c(
"Discounted density integral",
"Integration-by-parts survival integral"
),
Result = c(unit_second_moment, unit_second_moment_check),
Difference_from_first = c(
0,
unit_second_moment_check - unit_second_moment
)
)
knitr::kable(
verification_second_moment,
digits = 15,
caption = "Independent verification of the unit-benefit second moment."
)| Method | Result | Difference_from_first |
|---|---|---|
| Discounted density integral | 0.0588302 | 0 |
| Integration-by-parts survival integral | 0.0588302 | 0 |
t_survival <- seq(0, 60, length.out = 1000)
plot(
t_survival,
survival(t_survival),
type = "l",
lwd = 3,
col = "#276FBF",
ylim = c(0, 1),
xlab = "Years after age 40",
ylab = "Survival probability",
main = "Survival Probability for a Life Aged 40"
)
grid()Modelled probability that a life aged 40 survives for at least another t years.
The model gives
\[ \Pr(T_{40}>60)={}_{60}p_{40} =0.235604. \]
Thus the probability of surviving from age 40 to age 100 is approximately 23.56%.
t_density <- seq(0, 80, length.out = 1200)
density_values <- density(t_density)
plot(
t_density,
density_values,
type = "l",
lwd = 3,
col = "#238636",
xlab = "Future lifetime t (years)",
ylab = expression(f[T[40]](t)),
main = "Probability Density of Future Lifetime"
)
polygon(
c(t_density, rev(t_density)),
c(density_values, rep(0, length(density_values))),
col = adjustcolor("#54AEFF", alpha.f = 0.30),
border = NA
)
lines(t_density, density_values, lwd = 3, col = "#238636")
grid()Probability density of the future lifetime of a life aged 40.
The density peaks at approximately 54.82 future years, corresponding to an age of approximately 94.82.
For a general force of interest \(\delta\), define
\[ V(\delta) =500{,}000E[e^{-\delta T_{40}}]. \]
delta_grid <- seq(0, 0.10, by = 0.01)
interest_results <- data.frame(
Force = delta_grid
)
interest_results$Unit_EPV <- vapply(
interest_results$Force,
discounted_moment_density,
numeric(1)
)
interest_results$Benefit_EPV <- benefit * interest_results$Unit_EPV
interest_table <- transform(
interest_results,
Force = paste0(sprintf("%.0f", 100 * Force), "%"),
Unit_EPV = round(Unit_EPV, 6),
Benefit_EPV = paste0(
"$",
format(round(Benefit_EPV), big.mark = ",", scientific = FALSE)
)
)
knitr::kable(
interest_table,
col.names = c("Force of interest", "Unit-benefit EPV", "Benefit EPV"),
caption = "Expected present value as the force of interest varies."
)| Force of interest | Unit-benefit EPV | Benefit EPV |
|---|---|---|
| 0% | 1.000000 | $500,000 |
| 1% | 0.629939 | $314,970 |
| 2% | 0.408194 | $204,097 |
| 3% | 0.272669 | \(136,335 | |4% | 0.188069|\) 94,035 |
| 5% | 0.134056 | $ 67,028 |
| 6% | 0.098745 | $ 49,372 |
| 7% | 0.075086 | $ 37,543 |
| 8% | 0.058830 | $ 29,415 |
| 9% | 0.047375 | $ 23,687 |
| 10% | 0.039098 | $ 19,549 |
delta_plot <- seq(0, 0.10, length.out = 151)
epv_plot <- benefit * vapply(
delta_plot,
discounted_moment_density,
numeric(1)
)
plot(
100 * delta_plot,
epv_plot,
type = "l",
lwd = 3,
col = "#CF222E",
xlab = "Force of interest (%)",
ylab = "Expected present value ($)",
main = "Expected Present Value versus Force of Interest",
yaxt = "n"
)
axis(
2,
at = pretty(epv_plot),
labels = paste0("$", format(pretty(epv_plot), big.mark = ","))
)
grid()
points(
100 * delta,
benefit_epv,
pch = 19,
cex = 1.3,
col = "#0969DA"
)
text(
100 * delta,
benefit_epv,
labels = paste0(" 4%: $", format(round(benefit_epv), big.mark = ",")),
pos = 4,
col = "#0969DA"
)Expected present value of the death benefit as a function of the force of interest. The point marks the baseline 4% assumption.
These rates are forces of interest. A force \(\delta\) corresponds to an annual effective rate
\[ i=e^\delta-1. \]
For example, a 10% force is equivalent to an annual effective rate of approximately 10.52%.
Differentiating the expected present value with respect to \(\delta\),
\[ \boxed{ V'(\delta) =-500{,}000E[T_{40}e^{-\delta T_{40}}]<0. } \]
Since \(T_{40}>0\), the derivative is negative. Increasing the interest rate therefore always decreases the expected present value.
The second derivative is
\[ \boxed{ V''(\delta) =500{,}000E[T_{40}^2e^{-\delta T_{40}}]>0. } \]
The EPV is therefore a decreasing but convex function of the force of interest.
Economically, when the insurer can earn a higher investment return, it needs less money today to finance a benefit expected to be paid many years in the future.
Each mortality parameter is now varied separately by \(\pm10\%\), while the other two parameters remain fixed.
parameter_scenarios <- data.frame(
Parameter = rep(c("A", "B", "c"), each = 2),
Change = rep(c(-0.10, 0.10), times = 3),
stringsAsFactors = FALSE
)
parameter_scenarios$New_parameter <- mapply(
function(parameter, change) {
baseline <- switch(
parameter,
A = A,
B = B,
c = c_base
)
baseline * (1 + change)
},
parameter_scenarios$Parameter,
parameter_scenarios$Change
)
scenario_epv <- function(parameter, change) {
A_scenario <- A
B_scenario <- B
c_scenario <- c_base
if (parameter == "A") A_scenario <- A * (1 + change)
if (parameter == "B") B_scenario <- B * (1 + change)
if (parameter == "c") c_scenario <- c_base * (1 + change)
benefit * discounted_moment_density(
delta,
A_par = A_scenario,
B_par = B_scenario,
c_par = c_scenario
)
}
parameter_scenarios$EPV <- mapply(
scenario_epv,
parameter_scenarios$Parameter,
parameter_scenarios$Change
)
parameter_scenarios$Percent_change_in_EPV <- 100 * (
parameter_scenarios$EPV / benefit_epv - 1
)
sensitivity_table <- transform(
parameter_scenarios,
Change = paste0(ifelse(Change > 0, "+", ""), sprintf("%.0f%%", 100 * Change)),
New_parameter = format(New_parameter, digits = 8, scientific = FALSE),
EPV = paste0(
"$",
format(round(EPV, 2), big.mark = ",", nsmall = 2)
),
Percent_change_in_EPV = paste0(
ifelse(Percent_change_in_EPV > 0, "+", ""),
sprintf("%.3f%%", Percent_change_in_EPV)
)
)
knitr::kable(
sensitivity_table,
col.names = c(
"Parameter",
"Parameter change",
"New parameter",
"Benefit EPV",
"Change in EPV"
),
caption = "One-at-a-time mortality-parameter sensitivity."
)| Parameter | Parameter change | New parameter | Benefit EPV | Change in EPV |
|---|---|---|---|---|
| A | -10% | 0.0004500 | $ 93,693.81 | -0.362% |
| A | +10% | 0.0005500 | $ 94,374.83 | +0.362% |
| B | -10% | 0.0000675 | $ 89,741.68 | -4.565% |
| B | +10% | 0.0000825 | $ 98,075.42 | +4.297% |
| c | -10% | 0.9675000 | $ 6,307.14 | -93.293% |
| c | +10% | 1.1825000 | $389,191.35 | +313.881% |
scenario_labels <- paste0(
parameter_scenarios$Parameter,
ifelse(parameter_scenarios$Change > 0, " +10%", " -10%")
)
scenario_colours <- ifelse(
parameter_scenarios$Percent_change_in_EPV >= 0,
"#238636",
"#CF222E"
)
bar_positions <- barplot(
parameter_scenarios$Percent_change_in_EPV,
names.arg = scenario_labels,
las = 2,
col = scenario_colours,
border = NA,
ylab = "Change in expected present value (%)",
main = "Mortality-Parameter Sensitivity",
ylim = c(
min(parameter_scenarios$Percent_change_in_EPV) - 30,
max(parameter_scenarios$Percent_change_in_EPV) + 45
)
)
abline(h = 0, col = "#24292F", lwd = 1)
text(
bar_positions,
parameter_scenarios$Percent_change_in_EPV,
labels = sprintf("%+.1f%%", parameter_scenarios$Percent_change_in_EPV),
pos = ifelse(parameter_scenarios$Percent_change_in_EPV >= 0, 3, 1),
cex = 0.8
)Percentage change in EPV under one-at-a-time mortality-parameter stresses. Note the much wider scale for c.
Increasing \(A\), \(B\), or \(c\) generally brings deaths forward. Earlier payments are discounted for less time, so their expected present value is higher.
The value is only mildly sensitive to \(A\), because its baseline level is small. The value is more sensitive to \(B\), which controls the scale of the age-dependent mortality component.
The parameter \(c\) has an especially large effect because it is raised to attained age:
\[ \mu(40+t)=A+Bc^{40+t}. \]
In particular,
\[ \frac{\partial(Bc^x)}{\partial c}=Bx c^{x-1} \]
and
\[ \boxed{ \frac{\partial\log(Bc^x)}{\partial\log c}=x. } \]
At age 40, a 1% local change in \(c\) therefore produces approximately a 40% local change in the Gompertz component of mortality.
There is an important modelling warning. Reducing \(c=1.075\) by 10% gives \(c=0.9675<1\). The Gompertz component then decreases with age rather than increases. The calculation is a valid mechanical stress test, but it is no longer a conventional adult human-mortality model.
Because \(c\) acts exponentially, a useful supplementary test is to vary \(\log(c)\) by \(\pm10\%\). This keeps \(c>1\), so mortality continues to increase with age.
log_c_changes <- c(-0.10, 0.10)
local_c_values <- exp(log(c_base) * (1 + log_c_changes))
local_c_epv <- vapply(
local_c_values,
function(c_value) {
benefit * discounted_moment_density(delta, c_par = c_value)
},
numeric(1)
)
local_c_table <- data.frame(
Change_in_log_c = paste0(
ifelse(log_c_changes > 0, "+", ""),
sprintf("%.0f%%", 100 * log_c_changes)
),
New_c = local_c_values,
Benefit_EPV = local_c_epv,
Change_in_EPV = 100 * (local_c_epv / benefit_epv - 1)
)
knitr::kable(
local_c_table,
digits = c(0, 6, 2, 2),
col.names = c(
"Change in log(c)",
"New c",
"Benefit EPV ($)",
"Change in EPV (%)"
),
caption = "Supplementary sensitivity that preserves increasing mortality."
)| Change in log(c) | New c | Benefit EPV ($) | Change in EPV (%) |
|---|---|---|---|
| -10% | 1.067254 | 74108.36 | -21.19 |
| +10% | 1.082803 | 115738.04 | 23.08 |
The primary calculation uses adaptive numerical quadrature. The mathematical integral is over \([0,\infty)\), while the baseline numerical calculation truncates at 150 future years, by which point the survival probability is numerically indistinguishable from zero:
\[ M(k) =\int_0^\infty e^{-k\delta t} \left(A+Bc^{40+t}\right) \exp\left[ -At-\frac{Bc^{40}}{\log c}(c^t-1) \right]dt. \]
The moments are then
\[ E[Z]=500{,}000M(1) \]
and
\[ E[Z^2]=(500{,}000)^2M(2). \]
As a second numerical method, integration by parts gives
\[ M(k) =1-k\delta\int_0^\infty e^{-k\delta t}{}_tp_{40}\,dt. \]
The two methods integrate different functions but produce effectively identical results. This provides a check against numerical truncation and integration error.
This policy promises \(\$500{,}000\), but the insurer does not expect to pay that amount immediately. Under the mortality model, the insured’s expected remaining lifetime is about 47.5 years. Discounting the distant and uncertain payment at a 4% force reduces its expected value today to approximately \(\$94,035\).
The present value is uncertain because the insured could die next year or many decades from now. That timing uncertainty produces a standard deviation of approximately \(\$76,584\).
Higher interest rates reduce the policy value because money invested today has more time to grow before the claim is paid. Higher mortality normally raises the value because the claim is expected to be paid sooner and is therefore discounted for less time.
Finally, the sensitivity analysis illustrates an important actuarial lesson: parameters inside an exponential expression must be stress-tested carefully. Apparently modest parameter changes can produce very large valuation changes.
This is a deliberately stylised educational valuation. In practice, an actuary would also consider:
The results should therefore be interpreted as a model-based net present value, not as a quoted premium or prediction for a particular individual.
All numerical results and figures in this article are generated when the document is rendered. The R session information is included below.
## R version 4.6.0 (2026-04-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Sequoia 15.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Australia/Melbourne
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0 xfun_0.57
## [5] cachem_1.1.0 knitr_1.51 htmltools_0.5.9 rmarkdown_2.31
## [9] lifecycle_1.0.5 cli_3.6.6 sass_0.4.10 jquerylib_0.1.4
## [13] compiler_4.6.0 rstudioapi_0.18.0 tools_4.6.0 evaluate_1.0.5
## [17] bslib_0.11.0 yaml_2.3.12 rlang_1.3.0 jsonlite_2.0.0