ABSTRACT (original article): Equity-linked insurance products have gained popularity in recent years as retirement products with investment benefits. However, the design and pricing of these products for couples have been overlooked despite empirical evidence showing positive dependence in a couple’s lifetimes. In this paper, we propose some suitable products for couples where the benefits depend on the death times of both lives, and perform valuation using the discounted density approach while allowing the lifetimes to be dependent. By modeling the lifetimes with a bivariate mixed Erlang distribution, closed-form pricing formulas are developed for a variety of benefit types such as income protection for the last survivor (possibly with roll-up guarantee or benefit indexation) and dynamic fund protection/withdrawals. Fitting of bivariate lifetime data is also discussed in relation to bivariate Laguerre series. The impact of dependence on the prices of these products is demonstrated via numerical examples. In particular, our results suggest that incorrectly assuming independence between lifetimes would overprice these products compared to the actual situation of positive dependence, thereby making the products less attractive. CITATION (original article): Tang, K., Cheung, E. C. K., & Woo, J.-K. (2025). Designing and valuing new equity-linked insurance products for couples. Insurance: Mathematics and Economics, 121, 111-132. https://doi.org/10.1016/j.insmatheco.2025.01.003
This notebook is based on the paper above. We propose a variety of equity-linked products where the benefits depend on the times of death of both lives. One such product is a reversionary annuity benefit paid to the last survivor, and we develop closed-form pricing formulas by modelling the joint lifetimes using a bivariate mixed Erlang distribution. This notebook presents the numerical example from Section 5.1 of the paper for the special case of a pure scale mixture of Erlangs.
Bivariate mixed Erlang distribution
Bivariate mixed Erlang distribution
Consider a couple consisting of two persons aged and with remaining lifetimes and respectively. We assume that, for some positive integers , , and , the joint distribution of can be expressed as a bivariate mixed Erlang distribution (shape and scale mixture) with density defined in Equation (2.1). The key is in the mixing weights , which captures the dependence information in the lifetimes.
z
1
z
2
T
z
1
T
z
2
r
1
r
2
l
1
l
2
(,)
T
z
1
T
z
2
q
i
1
i
2
j
1
j
2
For the special case ==1, we have a pure scale bivariate mixture of Erlangs, which reduces to a bivariate mixture of exponentials. We present a hypothetical example to analyse the impact of dependence between the couple’s lifetimes on the prices of various equity-linked products. In Section 4.2, we demonstrate that if we put together two univariate combinations of exponentials with a Farlie-Gumbel-Morgenstern copula, the resulting joint density is a bivariate combination of exponentials. To illustrate, we start by calculating the Erlang mixing weights from the exponential rates using the example in Section 5.1.
r
1
r
2
In[]:=
(*Mixedexponentialrates*)λ1={7/100,9/100,11/100,13/100};λ2={5/100,7/100,10/100,12/100};m1=Length[λ1];m2=Length[λ2];(*a1anda2arecalculatedsuchthatandarehypoexponentialdistributions*)a1=TableIfj≠i,,1,{i,m1};a2=TableIfj≠i,,1,{i,m2};
T
z
1
T
z
2
m1
∏
j=1
λ1[[j]]
λ1[[j]]-λ1[[i]]
m2
∏
j=1
λ2[[j]]
λ2[[j]]-λ2[[i]]
(*CalculatemixedErlangweights*)i1[j1_]:=i1[j1]=Ceiling[j1/m1];i2[j2_]:=i2[j2]=Ceiling[j2/m2];k1[j1_]:=k1[j1]=j1-m1(i1[j1]-1);k2[j2_]:=k2[j2]=j2-m2(i2[j2]-1);
q
q
a1[[i1[j1]]]a2[[i2[j2]]]a1[[k1[j1]]]a2[[k2[j2]]](1+θ) | 1≤j1≤ 2 m1 2 m2 |
a1[[i1[j1]]]a2[[i2[j2]-m2]]a1[[k1[j1]]]a2[[k2[j2]]] -2θλ2[[i2[j2]-m2]] λ2[[i2[j2]-m2]]+λ2[[k2[j2]]] | 1≤j1≤ 2 m1 2 m2 2 m2 |
a1[[i1[j1]-m1]]a2[[i2[j2]]]a1[[k1[j1]]]a2[[k2[j2]]] -2θλ1[[i1[j1]-m1]] λ1[[i1[j1]-m1]]+λ1[[k1[j1]]] | 2 m1 2 m1 2 m2 |
a1[[i1[j1]-m1]]a2[[i2[j2]-m2]]a1[[k1[j1]]]a2[[k2[j2]]] 4θλ1[[i1[j1]-m1]]λ2[[i2[j2]-m2]] (λ1[[i1[j1]-m1]]+λ1[[k1[j1]]])(λ2[[i2[j2]-m2]]+λ2[[k2[j2]]]) | 2 m1 2 m1 2 m2 2 m2 |
q[j1_,j2_,θ_]:=q[j1,j2,θ]=If1≤j1≤2&&2+1≤j2≤2+2,[j1,j2-2,θ],0
2
m1
2
m1
2
m1
2
m2
q
2
m1
In[]:=
(*CalculatemixedErlangscaleparameters*)[j1_]:=[j1]=
[j2_]:=[j2]=
ℓ1=2;ℓ2=2+2;Doγ[j]=[j],{j,1,2}Doγ[j]=[j-2],{j,2+1,2+2}
γ
1
γ
1
λ1[[i1[j1]]] | 1≤j1≤ 2 m1 |
λ1[[i1[j1]-m1]]+λ1[[k1[j1]]] | 2 m1 2 m1 |
γ
2
γ
2
λ2[[i2[j2]]] | 1≤j2≤ 2 m2 |
λ2[[i2[j2]-m2]]+λ2[[k2[j2]]] | 2 m2 2 m2 |
2
m1
2
m1
2
m2
γ
1
2
m1
γ
2
2
m1
2
m1
2
m1
2
m2
Income protection: basic products
Income protection: basic products
The products of interest have a payout that is influenced by the death times of both lives. We present a motivating example where the last survivor receives a guaranteed minimum income benefit (GMIB) starting from the time of the first death until his/her own death, and the benefit characteristics are linked to the performance of a stock or index fund. The time-0 value of such an income protection product is given in Equation (2.3). The value of this product can also be split into two expectations based on a specific individual dying first. Couples may be interested in insuring only one of these cases, for example, when the main income earner dies first. By symmetry, it is sufficient to evaluate one of these expectations in Equation (2.5). We assume that the underlying equity index follows a geometric Brownian motion given by , where . To appropriately discount the cash flows, let be the risk-free force of interest and is set for risk-neutral pricing.
{S(t)}
t>=0
S(t)=S(0)
X(t)
X(t)=μt+σW(t)
δ
μ
In[]:=
(*Stockdynamics*)δ=3/100;σ=15/100;μ=δ-;
2
σ
2
At the moment, the value of the product is an expectation over the couple’s uncertain times of death and an uncertain index price. The key to our valuation is through the use of the discounted density method, which calculates the density function of the index process stopped at a random time and is appropriately discounted to present time (hence discounted density function). As described in Proposition 1, the value of the income protection product can be expressed as sums of integrals over only the index prices, thus simplifying the analysis. A detailed method is provided in Section 3.1. Since the discounted density of a geometric Brownian motion stopped at an Erlang time is known to have analytic form, the resulting integrals over index prices can be solved.
GMIB based on index value at the first death
GMIB based on index value at the first death
In Section 3.1.1, we demonstrate how the required integral over the discounted index prices can be evaluated analytically for a specific GMIB product that pays the higher of the index value and some predetermined guaranteed minimum amount starting from the time of the first death. When the initial index value exceeds the guaranteed amount , the benefit integral is given by Equation (3.7), and when , we get Equation (3.8).
S(0)
K
S(0)<K
In[]:=
(*Benefitintegralfor3.1.1withm=1*)Benefit311[γ_,δ_,μ_,σ_,S0_,K_]:=Benefit311[γ,δ,μ,σ,S0,K]=Blockα=(1/)-μ-+2**(γ+δ),β=(1/)-μ++2**(γ+δ),κ:=(γ/(γ+δ))*((-α*β)/(β-α)),IfK⩾S0,κ*K(-α)+Kβ+*(1/(β-1)-1/β),κ*S0(β-1)+IfK>0,S0(1-α)+(/)*(1/(-α)-1/(1-α)),0(*Valuewhenlife1diesfirst*)V311x[δ_,μ_,σ_,θ_,S0_,K_]:=V311x[δ,μ,σ,θ,S0,K]=**Benefit311[γ[j1]+γ[j2],δ,μ,σ,S0,K](*Valuewhenlife2diesfirst*)V311y[δ_,μ_,σ_,θ_,S0_,K_]:=V311y[δ,μ,σ,θ,S0,K]=**Benefit311[γ[j1]+γ[j2],δ,μ,σ,S0,K]
2
σ
2
μ
2
σ
2
σ
2
μ
2
σ
β
S0
β-1
K
α
S0
α-1
K
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j2]+δ
γ[j1]
γ[j1]+γ[j2]
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j1]+δ
γ[j2]
γ[j1]+γ[j2]
GMIB based on running maximum of index value at the first death
GMIB based on running maximum of index value at the first death
Section 3.1.2 proposes a similar annuity product where the payment rate is based on the running maximum of the index process at the time of the first death. Such features are common guarantees to lock in the optimal price of the index, minimising policyholder regret if the index decreases prior to the first death.
In[]:=
(*Benefitintegralfor3.1.2withm=1*)Benefit312[γ_,δ_,μ_,σ_,S0_,K_]:=Benefit312[γ,δ,μ,σ,S0,K]=Blockα=(1/)-μ-+2**(γ+δ),β=(1/)-μ++2**(γ+δ),κ:=(γ/(γ+δ))*((-α*β)/(β-α)),IfK⩾S0,κ*(-α+β)(-α)*Kβ+*(1/(β-1)-1/β),κ*(-α+β)(-α)(S0/(β-1))(*Valuewhenlife1diesfirst*)V312x[δ_,μ_,σ_,θ_,S0_,K_]:=V312x[δ,μ,σ,θ,S0,K]=**Benefit312[γ[j1]+γ[j2],δ,μ,σ,S0,K](*Valuewhenlife2diesfirst*)V312y[δ_,μ_,σ_,θ_,S0_,K_]:=V312y[δ,μ,σ,θ,S0,K]=**Benefit312[γ[j1]+γ[j2],δ,μ,σ,S0,K]
2
σ
2
μ
2
σ
2
σ
2
μ
2
σ
β
S0
β-1
K
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j2]+δ
γ[j1]
γ[j1]+γ[j2]
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j1]+δ
γ[j2]
γ[j1]+γ[j2]
GMIB based on index value and its running maximum at the first death
GMIB based on index value and its running maximum at the first death
Section 3.1.3 generalises the benefit function, where the payment rate is based on both the index value and its running maximum at the time of the first death. Rather than a fixed guarantee, this benefit can be considered as a floating guarantee that increases as a fraction of the running maximum.
In[]:=
(*Benefitintegralfor3.1.3withm=1*)Benefit313[γ_,δ_,μ_,σ_,S0_,ω_]:=Benefit313[γ,δ,μ,σ,S0,ω]=Blockα=(1/)-μ-+2**(γ+δ),β=(1/)-μ++2**(γ+δ),κ:=(γ/(γ+δ))*((-α*β)/(β-α)),S0*κ*(-α+β)*(β-1)*(1/(-α)-1/(1-α))+1((1-α)(β-1))(*Valuewhenlife1diesfirst*)V313x[δ_,μ_,σ_,θ_,S0_,ω_]:=V313x[δ,μ,σ,θ,S0,ω]=**Benefit313[γ[j1]+γ[j2],δ,μ,σ,S0,ω](*Valuewhenlife2diesfirst*)V313y[δ_,μ_,σ_,θ_,S0_,ω_]:=V313y[δ,μ,σ,θ,S0,ω]=**Benefit313[γ[j1]+γ[j2],δ,μ,σ,S0,ω]
2
σ
2
μ
2
σ
2
σ
2
μ
2
σ
1-α
ω
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j2]+δ
γ[j1]
γ[j1]+γ[j2]
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j1]+δ
γ[j2]
γ[j1]+γ[j2]
GMIB based on index value and its running minimum at the first death
GMIB based on index value and its running minimum at the first death
Section 3.1.4 considers a GMIB based on the index value and its running minimum at the time of the first death. This benefit improves when the difference between the index value and its running minimum is large, minimising policyholder regret from purchasing the contract when the index price is initially high. The final result is expressed in terms of integrals that can be evaluated within Mathematica.
In[]:=
(*Benefitfunctionfor3.1.4withm=1*)Benefit314[γ_,δ_,μ_,σ_,S0_,K_]:=Benefit314[γ,δ,μ,σ,S0,K]=Blockα=(1/)-μ-+2**(γ+δ),β=(1/)-μ++2**(γ+δ),κ:=(γ/(γ+δ))*((-α*β)/(β-α)),IfK⩾S0,κ*(-α+β)*Kβ*1(-α)-NIntegrate,{y,0,Infinity},PrecisionGoal->12+2*S0(β-1)*NIntegrate,{y,0,Infinity},PrecisionGoal->12-S0β*NIntegrate,{y,0,Infinity},PrecisionGoal->12,Null(*Valuewhenlife1diesfirst*)V314x[δ_,μ_,σ_,θ_,S0_,K_]:=V314x[δ,μ,σ,θ,S0,K]=**Benefit314[γ[j1]+γ[j2],δ,μ,σ,S0,K](*Valuewhenlife2diesfirst*)V314y[δ_,μ_,σ_,θ_,S0_,K_]:=V314y[δ,μ,σ,θ,S0,K]=**Benefit314[γ[j1]+γ[j2],δ,μ,σ,S0,K]
2
σ
2
μ
2
σ
2
σ
2
μ
2
σ
α*y-β*Log[1/2*(1+K/S0*)]
y
-(1-α)*y-(β-1)*Log[1/2*(1+K/S0*)]
y
-(1-α)*y-β*Log[1/2*(1+K/S0*)]
y
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j2]+δ
γ[j1]
γ[j1]+γ[j2]
ℓ1
∑
j1=1
ℓ2
∑
j2=1
q[j1,j2,θ]
γ[j1]+δ
γ[j2]
γ[j1]+γ[j2]
Income protection: products with non-level guarantee or payment
Income protection: products with non-level guarantee or payment
A common feature of life insurance contracts is that the coverage or benefits increase by a predetermined percentage to keep up with the cost of living, especially since these contracts can be in effect for long periods.
Guarantee: roll-up
Guarantee: roll-up
In Section 3.2.1, we modify the previous income protection product in Section 3.1.1 where the guaranteed rate increases by a ‘roll-up’ rate from policy inception until the time of the first death. Note that the form of the benefit integral is the same but with the parameters changed.
Payment: benefit indexation
Payment: benefit indexation
In the previous roll-up guarantee, the guaranteed rate increases until the time of the first death, but the subsequent payments to the last survivor stays level. Instead, Section 3.2.2 has the indexation applied once the payments begin. We apply benefit indexation to the product in Section 3.1.1, noting that all modifications appear outside the Section 3.1.1 benefit integral.
Payment: increasing annuity
Payment: increasing annuity
Instead of benefit indexation that increases exponentially, our methodology is also applicable to a linearly increasing payment rate. This may be useful when the surviving spouse does not require a high income rate immediately after the first death. The result from Section 3.2.3 when an increasing annuity is applied to the product in Section 3.1.1 is expressed in terms of the Section 3.1.1 benefit integral with modified parameters.
Death benefits: products with bequest
Death benefits: products with bequest
Dynamic fund protection
Dynamic fund protection
Dynamic fund withdrawals
Dynamic fund withdrawals
Tables in the paper
Tables in the paper
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Designing and valuing new equity-linked insurance products for couples
by Kelvin Tang, Eric C.K. Cheung and Jae-Kyung Woo
Wolfram Community, STAFF PICKS, May 5, 2026
https://community.wolfram.com/groups/-/m/t/3712289
by Kelvin Tang, Eric C.K. Cheung and Jae-Kyung Woo
Wolfram Community, STAFF PICKS, May 5, 2026
https://community.wolfram.com/groups/-/m/t/3712289