2. Dirichlet processes and Dirichlet process mixtures

The Dirichlet process (DP) was first introduced by Ferguson (Reference Ferguson1973) to specify a probability model for the random masses of a discrete distribution G. We thus say that G follows a Dirichlet process on the measurable space $\left(\Gamma,\mathcal{G} \right)$ if for any measurable partition $\left(W_1,\ldots, W_q \right)$ of $\Gamma$ , then

\begin{align*} \left(G\left(W_1 \right), \ldots, G\left(W_q \right)\right) \sim \text{Dirichlet}\left(\phi G_0\left(W_1 \right), \ldots, \phi G_0\left(W_q \right) \right)\end{align*}

which we write $G \sim \text{DP}\left(G_0,\phi\right)$ . Hence, a random draw from the DP yields an almost sure discrete probability distribution over a countably infinite number of points drawn independently from a base distribution denoted as $G_0$ , which specifies the DP together with the concentration parameter $\phi$ (see Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013, Chapter 23 and De Iorio et al., Reference De Iorio, Müller, Rosner and MacEachern2004). The $\phi$ parameter captures the degree of shrinkage of G toward $G_0$ , or in other words, the strength of the prior assumption $G_0$ over G, analogous to the prior assumption about the parameters of a probability distribution in Bayesian statistics.

Sethuraman (Reference Sethuraman1994) outlines a construction of G as a mixture distribution with a countably infinite number of components, indexed by k:

(2.1) \begin{align}G= \displaystyle\sum_{k=1}^{\infty} \pi_k \delta_{\gamma_k},\end{align}

where $\gamma_k \overset{\mathrm{i.i.d.}}{\sim} G_0$ for $k=1,2,\ldots$ , and $\delta_{\gamma_k}$ denotes the Dirac measure that assigns unitary mass if $\gamma=\gamma_k$ and zero otherwise. The mixture weights $\pi_k$ are randomly generated using the so called stick breaking procedure (SBP), which rescales a set of i.i.d. random variables $\psi_k \overset{\mathrm{i.i.d.}}{\sim} \text{Beta}\left(1,\phi\right)$ as follows:

(2.2) \begin{align} \pi_k \equiv \pi_k \left(\psi_{1:k} \right) = \psi_k \displaystyle\prod_{j=1}^{k-1} \left(1 - \psi_j \right).\end{align}

where $\psi_{1:k} = \left(\psi_{1},\ldots,\psi_k \right)$ . The “stick breaking” definition refers to the decreasing size of the mixture weights as the index k increases. From this characterization, we can observe how a random draw from a Dirichlet process yields a discrete distribution over a countably infinite number of atoms from $G_0$ .

In this paper, we are interested in a flexible model for the probability distribution of a random variable $T_i$ , eventually multivariate, for the ith unit, whose density $f\left(\cdot;\;\beta,\gamma_i\right)$ is indexed by a global parameter vector $\beta$ , common to all units, and by a unit-specific parameter vector $\gamma_i$ . Unit-specific parameters, commonly referred as random effects in biostatistics, are introduced in a probability model to characterize the heterogeneity among the units due to unobservable variables.

To achieve this goal, we assume that $\gamma_i$ are drawn from a discrete distribution G, which follows a Dirichlet process. Therefore, convolving the density of $T_i$ with $G \sim \text{DP}\left(G_0,\phi \right)$ yields a Dirichlet process mixture (DPM) model (Lo, Reference Lo1984):

(2.3) \begin{align} f \left(t \mid \beta, G \right) = \int_{\Omega{\gamma}} f \left(t \mid \beta, \gamma \right) \mathrm{d}G\left(\gamma\right) = \displaystyle\sum_{k=1}^{\infty} \pi_k f \left(t \mid \beta,\gamma_k \right)\end{align}

where $\Omega{\gamma}$ denotes the sample space of $\gamma$ .

The DPM model hereby defined allows for a more flexible distribution of $T_i$ , which can capture complex features in the data, such as fat tails and multimodality, as opposed to the case where we specify a parametric model for $\gamma_i$ , such as the normal distribution.

This DP construction is flexible since despite its apparent nonparametric nature, it regularizes the distribution of G toward a simple parametric form $G_0$ through the concentration parameter $\phi$ .

From another perspective, a discrete distribution for $\gamma_i\sim G$ is akin to the creation of ties among the units, which can configure clusters of observations with the same values of $\gamma$ . Let $s_i=k$ to indicate that the ith unit belongs to the kth cluster characterized by the parameter $\gamma^*_k$ . We use the superscript $^*$ to denote the common values of $\gamma_i$ across the n units in the sample. The clustering procedure, tuned by the parameter $\phi$ , allows to sequentially group the observations through a sampling process known in the literature as the Chinese restaurant process (see Heinz, Reference Heinz2014 for an illustration). In words, as we observe the units in a sequence, these are more likely to be in a certain class with a probability that depends on the number of units already therein, or to belong to a newly created class with a probability which depends on $\phi$ (see Blackwell and MacQueen, Reference Blackwell and MacQueen1973 for a characterization in terms of the Pólya urn distribution).

The use of Dirichlet processes can also be seen as a method to infer the number of components of a mixture distribution, as opposed to strategies based on appropriate model selection criteria (see Ungolo and van den Heuvel, Reference Ungolo and van den Heuvel2022 for a discussion). Their advantage is to avoid the specification and the estimation of several models with different numbers of mixture components, which can be time consuming.

Summing up, the DPM model assumes the following data generating process for a sample $t_1,\ldots,t_n$ :

(2.4) \begin{equation}\begin{aligned}& \pi_k \mid \phi \sim \text{SBP}\left(\phi\right), \quad \text{ }k=1,2,\ldots; && \\[5pt] & \gamma^*_k \mid G_0 \overset{\mathrm{i.i.d.}}{\sim} G_0, \quad \text{ }k=1,2,\ldots; && \\[5pt] & s_i \mid \pi_1, \pi_2,\ldots \overset{\mathrm{i.i.d.}}{\sim} \text{Discrete}\left(\pi_1, \pi_2,\ldots \right), \quad \text{ }i=1,\ldots,n; &&\\[5pt] & t_i \mid \beta,\gamma_{s_i}^* \overset{\mathrm{i.i.d.}}{\sim} f\left(t_i \mid \beta,\gamma_{s_i}^* \right) \quad \text{ } i=1,\ldots,n.\end{aligned}\end{equation}

3. The Augmented Variable DPM model

This section introduces the augmented variable Dirichlet process mixture (AVDPM) model for the analysis of nonexchangeable joint dependent lifetimes within a group, where individual as well as group-specific covariates are available. For example, husband and wife lifetimes are likely to be positively associated (Denuit et al., Reference Denuit, Dhaene, Le Bailly De Tilleghem and Teghem2001), since they share the same living conditions (e.g. diet, socioeconomic factors), are exposed to similar risks (e.g. during a catastrophic event they are likely to be in the same place), or eventually subject to the broken-heart syndrome (Parkes et al., Reference Parkes, Benjamin and Fitzgerald1969). Other examples include the joint lifetimes of the primary and secondary head of an insurance policy, families with husband, wife and one child, and so on. We describe the framework in the context of a model for the joint lifetime of male–female couples. In Section 8, we discuss how to extend the framework to more general cases of groups with different numbers of exchangeable lifetimes.

Let $T_i=\left(T_{i,1},\ T_{i,2}\right)$ , where $T_{i,1}$ and $T_{i,2}$ denote the random future lifetime of husband and wife, respectively, for the ith couple, with individual-specific vector of characteristics $x_{i,1}$ and $x_{i,2}$ (such as age, medical status, and so on) fixed and observable, and couple-specific covariates vector $Z_i$ , which can include for example household income and geo-demographic profile (an indicator of the socioeconomic status, see Ungolo et al., Reference Ungolo, Christiansen, Kleinow and MacDonald2019), which we treat as a random variable. Therefore, $Z_i$ includes any features that can explain the statistical association between the lifetimes of the husband and of the wife. For example, the household income can explain the heterogeneity in the longevity profile of each couple member since a wealthy couple can access better healthcare services all else being equal, compared to a deprived one. Similarly, the geo-demographic profile, capturing the effect of the area where the couple lives (e.g. urban or rural), can also be a proxy for the socioeconomic characteristics of a couple.

For the ith couple, we assume that conditional on a couple-specific bivariate random effect $\gamma_i=\left(\gamma_{i,1},\gamma_{i,2}\right)$ , then $T_{i,1}$ and $T_{i,2}$ are independently distributed. This is, for example, the approach followed in Ungolo and van den Heuvel (Reference Ungolo and van den Heuvel2022) when analyzing a joint model for the time to event and the time to censoring and by Ungolo and van den Heuvel (Reference Ungolo and van den Heuvel2024) where they develop a joint model for the time to competing risk events.

The AVDPM model “augments” the joint probability distribution of $\left(T_{i,1},T_{i,2}\right)$ by specifying a joint probability model for $\left(T_{i,1},T_{i,2}, Z_i\right)$ :

(3.1) \begin{align} f & \left(t_{i,1},t_{i,2}, z_i \mid x_{i,1},x_{i,2}; \beta_1,\beta_2,G\right)\nonumber \\[5pt] &= \int_{\Omega_{\gamma,\zeta}} \left[\displaystyle\prod_{j=1}^2 f\left(t_{i,j} \mid x_{i,j}; \beta_j, \gamma_{i,j}\right) \right] f \left(z_i \mid \zeta_i\right) \mathrm{d}G\left(\gamma_{i,1},\gamma_{i,2},\zeta_i \right)\end{align}

where $\beta_j$ is the parameter specific to the future lifetime density of the husband ( $j=1$ ) or of the wife ( $j=2$ ), $\zeta_i$ is the couple-specific parameter vector indexing the distribution of $Z_i$ (which can also include global parameters, i.e. not couple-specific), and $\Omega_{\gamma,\zeta}$ denotes the sample space of $\left(\gamma_1,\gamma_2,\zeta\right)$ . We generically denote by f the probability density function of continuous variables and the mass function of the discrete ones.

The joint distribution of Equation (3.1) allows to capture the dependence between $T_{i,1}$ and $T_{i,2}$ and between $\left(T_{i,1},T_{i,2} \right)$ and $Z_i$ through the joint distribution of $\left(\gamma_{i,1},\gamma_{i,2},\zeta_i \right)$ . The individual-specific fixed covariate vector $x_{i,j}$ directly enters the regression function of $T_{i,j}$ . On the other hand, the ith couple-specific covariate vector $Z_i$ is assumed to affect the individual lifetime distribution only indirectly, through $\left(\gamma_{i,1},\gamma_{i,2},\zeta_i \right)$ . Indeed, it is reasonable to assume that individual variables like medical status and age affect the individual mortality directly, while couple-specific variables like household income can have an indirect effect on mortality through the affordability of good quality healthcare services. This model differs from treating $Z_i$ as an exogenous variable, for example by including this observable variable among the set of fixed covariates within a regression function, such as in the paper of Deresa et al., (Reference Deresa, Van Keilegom and Antonio2022). Both approaches allows to estimate the effect of $Z_i$ on the distribution of $T_{i,j}$ , as shown in Equation (3.3) below. In addition, the density of Equation (3.1) allows to have an understanding of the differing strength of the statistical association between $T_{i,1}$ and $T_{i,2}$ as we show in the empirical analysis of Section 6.

We complete the model specification by assuming that the multivariate discrete distribution G is a random draw from a Dirichlet process, and write $G\sim \text{DP}\left(G_0,\phi \right)$ . For simplicity, we assume $G_0=G_{0,\gamma}\times G_{0,\zeta}$ , that is $\gamma$ and $\zeta$ are independently distributed in the base distribution. This specification still induces a dependence between $\gamma$ and $\zeta$ through the concentration parameter $\phi$ . A further simplification can be to specify $G_{0,\gamma}=G_{0,\gamma_1}\times G_{0,\gamma_2}$ , leaving the distribution of $\left(\gamma_{i,1},\gamma_{i,2},\zeta_i \right)$ fully tuned by $\phi$ . However, the former approach allows to account for the prior information from the researcher about the joint distribution of the parameters, including the statistical association between $\gamma_1$ and $\gamma_2$ , as more reasonable for the analysis of male–female couples.

An additional feature of this factorization, which enhances the flexibility of this joint model is that we can specify different parametric models for $T_{i,1}$ and $T_{i,2}$ , and that it allows for both positive and negative dependencies between lifetimes, similarly to the model of Lu (Reference Lu2017).

As in Equation (2.3), we can rewrite the density in (3.1) as a mixture distribution with an infinite number of components, obtaining the augmented variable DPM (AVDPM)model:

(3.2) \begin{align} f & \left(t_{i,1},t_{i,2}, z_i \mid x_{i,1},x_{i,2}; \beta_1,\beta_2,G\right) \nonumber \\[5pt] &= \displaystyle\sum_{k=1}^{\infty} \pi_k \left[\displaystyle\prod_{j=1}^2 f\left(t_{i,j} \mid x_{i,j}; \beta_j, \gamma_{k,j}^*\right) \right] f \left(z_i \mid \zeta_k^* \right)\end{align}

where $\pi_k$ ( $k=1,2,\ldots$ ) is the mixture weight characterized by the stick-breaking procedure described in Section 2, and the superscript $^*$ denotes the unique values of $\gamma_{i,j}$ and $\zeta_i$ .

The flexibility of this factorization allows to understand the impact of $Z_i$ on the distribution of $T_{i,j}$ (or of $\left(T_{i,1},T_{i,2}\right)$ ) by using standard probability calculus:

(3.3) \begin{align} f\left(t_{i,j} \mid x_{i,j},z_i;\;\beta_j, \gamma_{\cdot,j}^*, \zeta_{\cdot}^* \right) = \frac{\displaystyle\sum_{k=1}^{\infty} \pi_k f\left(t_{i,j} \mid x_{i,j};\;\beta_j, \gamma_{k,j}^* \right) f\left(z_i \mid \zeta_{k}^* \right)}{\displaystyle\sum_{k=1}^{\infty} \pi_k f\left(z_i \mid \zeta_{k}^* \right)} \,.\end{align}

where $\gamma_{\cdot,j}^* = (\gamma_{1,j}^*,\gamma_{2,j}^*,\ldots )$ and $\zeta_{\cdot}^* = (\zeta_{1}^*,\zeta_{2}^*,\ldots )$ .

In a similar fashion, we can derive the probability distribution of the time to death of the last survivor $T_{\overline{1,2}}=\text{max}\left(T_1,T_2 \right)$ . We omit the subscript i for notational convenience. The corresponding survivor function, simplistically denoted as $S_{\overline{x_1,x_2}}\left(t \mid z \right)$ , is useful especially in actuarial calculations, as we illustrate in Section 7:

(3.4) \begin{align} S_{\overline{x_1,x_2}}\left(t \mid z \right) := \mathbb{P} \left(T_{\overline{1,2}} \gt t \mid x_1,x_2,z \right) = 1 - \mathbb{P} \left(T_{1} \lt t, \enspace T_{2} \lt t \mid x_1,x_2,z \right),\end{align}

where

\begin{align}\nonumber \mathbb{P} \left(T_{1} \lt t, \enspace T_{2} \lt t \mid x_1,x_2,z \right) &= \frac{\displaystyle\sum_{k=1}^{\infty} \pi_k f\left(z \mid \zeta_{k}^* \right) \left[ \displaystyle\prod_{j=1}^2 \int_0^t f\left(u \mid x_{j};\;\beta_j, \gamma_{k,j}^* \right)\mathrm{d}u \right]} {\displaystyle\sum_{k=1}^{\infty} \pi_k f\left(z \mid \zeta_{k}^* \right)} \,. \end{align}

Analogous formula can be used for the joint life survival probability, denoted as $S_{x_1,x_2}\left( t \mid z\right)$ , characterizing the random variable $T_{1,2}=\text{min}\left(T_1,T_2\right)$ .

The DPM naturally clusters each couple into different classes. If we say that W is the random class allocation for a couple, we can calculate the probability of belonging to a certain class k conditional on the value of Z:

(3.5) \begin{align} \mathbb{P}\left(W=k \mid Z=z \right) = \frac{\pi_k f\left(z \mid \zeta_k^*\right) }{\displaystyle\sum_{h=1}^{\infty}\pi_h f\left(z \mid \zeta_h^*\right)} \,.\end{align}

In this way, we can understand how the various couples belonging to the kth class can share the same mortality profile and husband–wife mortality dependence relationships.

Compared to earlier literature in the field, the AVDPM model extends the mixture model analyzed by Ungolo and van den Heuvel (Reference Ungolo and van den Heuvel2022) since the authors consider a discrete random component (independent of Z) with unknown number of levels, chosen by means of a model selection procedure. The limitation of this approach is the need of estimating several models, which can be particularly time consuming for larger datasets as mentioned in Section 2. Ungolo and van den Heuvel (Reference Ungolo and van den Heuvel2024) overcome this issue by assuming that the random component is drawn from a discrete distribution following a Dirichlet process as in this paper. However, their joint lifetime model does not account for the statistical association among competing risks due to common factors.

To the best of our knowledge, the aforementioned papers (Section 1) of Youn and Shemyakin (Reference Youn and Shemyakin1999) and Dufresne et al. (Reference Dufresne, Hashorva, Ratovomirija and Toukourou2018) are the sole contributions accounting for the common variable Z (the absolute value of the age difference between the spouses and gender of the oldest partner) within the copula dependence parameter.

The factorization implied by the AVDPM model can also find application in the analysis of the time to competing causes of death, where Z accounts for those genetic factors affecting the dependence between the causes. Alternatively, the AVDPM model can be used to jointly model dependent frequency and/or severity of claims in nonlife and health insurance by type of event or line of business, without necessarily specifying a strong parametric assumption on the dependence, as can be the case with copula models.

4. Data and parametric model

We showcase the approach outlined in Section 3 to the analysis of the Canadian life insurance dataset initially studied by Frees et al. (Reference Frees, Carriere and Valdez1996) and then by Deresa et al. (Reference Deresa, Van Keilegom and Antonio2022) and Dufresne et al. (Reference Dufresne, Hashorva, Ratovomirija and Toukourou2018) among others. After some data processing operations, briefly described in Section 1 of the Supplementary Material, we have information about 13,482 joint and last-survivor annuity contracts in force between December 29, 1988 and December 31, 1993 (the observation period).

Specifically, we focus on a joint model for the lifetime distribution of individual members of male–female couples, for which we observe the starting date of the contract, the date of birth (thus their age at the start of the contract), and the date of death if any couple member dies within the observation period. This dataset contains a large number of censored units: indeed, we only observe 1,424 deaths among males and 500 deaths among females, while the remaining units are all right censored. In addition, most couples’ lifetime data are subject to left truncation. This means that we are able to observe the annuity contract only if both couple members are alive at the start of the observational period. Right censoring and left truncation must then be taken into account when deriving the likelihood function of the observations, as we have shown in Section 5.

In this dataset, the age at the start of the contract is the only individual-specific covariate (average of 65.54 for males and 62.64 for females), which we denote as $x_{i,j}$ for the jth member of the ith couple. As in Dufresne et al. (Reference Dufresne, Hashorva, Ratovomirija and Toukourou2018) and Youn and Shemyakin (Reference Youn and Shemyakin1999), we consider also the age difference (mean of 2.91 years) through the random variable $Z_i^A = \ln \left( |X_{i,1} - X_{i,2} | \right)$ and the indicator variable $Z_i^M$ , which is equal to 1 if the male is older than the female and 0 otherwise (the male is the oldest policyholder in 76.9% of the couples in the dataset). According to their analysis, a model including these two variables captures some additional features of the association between the lifetime of the husband and of the wife. More precisely, they claim that the larger the age difference, the lower the lifetime dependence. In addition, Dufresne et al. (Reference Dufresne, Hashorva, Ratovomirija and Toukourou2018) observe how the $Z_i^M$ covariate has an influence in the relationship between husband and wife, consequently affecting their lifetime dependence. Therefore, we will consider $Z_i^A$ and $Z_i^M$ as the two elements of the couple-specific covariates, $Z_i=\left(Z^A_i,Z^M_i\right)$ .

We randomly split the dataset into a training set, corresponding to 75% of the policyholders within the dataset (10,112 units), and use the remaining 25% to assess the out of sample performance of the model (3,370 units).

Following Frees et al. (Reference Frees, Carriere and Valdez1996) and Carriere (Reference Carriere2000), we specify a Gompertz-type hazard function to characterize the probability distribution of the male ( $j=1$ ) and the female ( $j=2$ ) lifetime within the ith policy:

(4.1) \begin{align} \mu \left(t \mid x_{i,j};\;\alpha_j, \beta_j, \gamma_{i,j}\right) = \exp \left(\alpha_j + \beta_j \left(x_{i,j} - \overline{x} + t\right) + \gamma_{i,j} \right),\end{align}

where $\overline{x}=70$ is set in order to decrease the posterior correlation between $\alpha_j$ and $\beta_j$ , which improves the convergence of the MCMC sampler. The choice of $\overline{x}$ is based on an indicative average value of x as calculated across males and females. For example, Ungolo et al. (Reference Ungolo, Kleinow and Macdonald2020) set $\overline{x}=77.5$ since they analyzed a pension scheme dataset with older scheme members.

We consider the male–female lifetime dependence by assuming a couple-specific frailty term $\exp\left(\gamma_{i,j}\right)$ and a joint model for $\left(\gamma_{i,1}, \gamma_{i,2}\right)$ . The parameter $\alpha_j$ denotes the log-baseline hazard function and $\beta_j$ measures the log-linear increase in the hazard function due to the individual age. The probability density function of $T_{i,j}$ can be written as:

(4.2) \begin{align} f \left(t \mid x_{i,j}; \alpha_j, \beta_j, \gamma_{i,j}\right) = \exp \biggl[ -\int_0^t \mu \left(s \mid x_{i,j};\;\alpha_j, \beta_j, \gamma_{i,j}\right) \mathrm{d}s \biggr] \mu \left(t \mid x_{i,j};\;\alpha_j, \beta_j, \gamma_{i,j}\right) .\end{align}

For the couple-specific covariates, we assume that $Z_i^A \sim \text{N}\left( \zeta_{i}^A,\sigma^2_A \right)$ and $Z_i^M \sim \text{Bernoulli}\left(\zeta^M_i\right)$ . The density of $T_{i,j}$ as a function of $x_{i,j}$ and $z_i$ follows from Equation (3.3). In this way, we can easily observe how the resulting joint model for $\left(T_{i,1},T_{i,2},Z_i\right)$ has an additional layer of flexibility since the effect of $Z_i$ on the hazard function (through dependent $\gamma_{i,j}$ , $\zeta^A_i$ and $\zeta^M_i$ ) is not necessarily proportional nor monotone.

The ith couple-specific parameters are drawn from a multivariate random discrete distribution G, $\left(\gamma_{i,1},\gamma_{i,2},\zeta_{i}^A,\zeta_{i}^M\right)\sim G$ , where G is a draw from a Dirichlet process with concentration parameter $\phi$ and base measure $G_0$ :

(4.3) \begin{align}G_0=\underbrace{g\left(\gamma_{i,1}, \gamma_{i,2} \mid \textbf{0}, \Sigma_{\gamma} \right)}_{\text{MVN}\left( \textbf{0}, \enspace \Sigma_{\gamma} \right)} \underbrace{g\left(\zeta_{i}^A \mid m_{A}, s_A^2 \right)}_{\text{N}\left( m_{A}, \enspace s_{A}^2 \right)}\underbrace{ g \left(\zeta_{i}^M \mid 13.31,4.44 \right)}_{\text{Beta} \left( 13.31,\ 4.44 \right)} \,.\end{align}

This model for the base distribution is motivated by the need to carry out a computationally efficient Bayesian inference by exploiting the conditional conjugacy of the parameters where possible, whilst keeping the model simple and flexible enough to capture the complex features of the data.

The parameters of the Beta distribution for $\zeta_{i}^M$ follow from an elicitation process where we specify a weakly informative distribution. Indeed, we expect that for most contracts the male couple member is the oldest in the couple; therefore, we set $\zeta_i^M$ to have mean 0.75 and standard deviation equal to 0.1 in the base distribution. In our analysis, we have also tried a prior distribution with a standard deviation of 0.05, which did not significantly impact the final results in terms of the convergence of the chain towards a stationary distribution, posterior distribution of the other parameters, and cluster allocation of the couples.

For all the other parameters of the base distribution, we assume that these are random to enhance the robustness of the inference, as from the prior elicitation process and the sensitivity analysis described in Section 5.1.

5. Inference

First of all, we approximate the mixture distribution of Equation (3.2) by setting an upper bound $K=25$ to the number of mixture components as in Ungolo and van den Heuvel (Reference Ungolo and van den Heuvel2024), which results in the truncated SBP of Ishwaran and James (Reference Ishwaran and James2001). Indeed, as motivated in Section 5.1, this is a conservative choice since only few components are sufficient to model the additional heterogeneity in the data, while all other components turn out to have a mixture weight nearly equal to zero. This simplifies the implementation of the Markov Chain Monte Carlo (MCMC) sampler compared to the use of the bound-free slice samplers of Walker (Reference Walker2007) and Kalli et al. (Reference Kalli, Griffin and Walker2011), or the retrospective sampler of Papaspiliopoulos and Roberts (Reference Papaspiliopoulos and Roberts2008).

5.2. Likelihood

Let $d_{i,j}$ denote an indicator variable that is equal to 1 if the male ( $j=1$ ) or the female couple member ( $j=2$ ) is observed to die throughout the observational study and 0 otherwise (hence, the lifetime variable is right censored). We assume that the censoring mechanism can be ignored since the censoring event for a couple member is caused by the end of the observation period (Type I censoring). Furthermore, we assume that the individual times to event for each couple member are independently distributed, conditional on the covariates and the multivariate random effect.

As discussed in Section 4, data are subject to left truncation since each contract is observable upon survival of both members at the start of the study. The level of left truncation for each couple is denoted by the variable $a_i$ , denoting the time (in years) from the start of the contract to the start of the observational study. This means that we need to work with the lifetime density function conditional on both couple members being alive at the beginning of the observation period. Note that $a_i$ is equal to zero if the contract starts during the observation period.

Let $\textbf{t}=\left(t_{1,1},t_{1,2},\ldots,t_{n,1},t_{n,2} \right)$ , $\textbf{x}=\left(x_{i,1}, x_{i,2},\ldots,x_{n,1}, x_{n,2} \right)$ , $\textbf{d}=\left(d_{1,1}, d_{1,2},\ldots,d_{n,1}, d_{n,2} \right)$ , $\textbf{z}^{\textbf{A}}=\left(z^A_1,\ldots,z^A_n \right)$ , $\textbf{z}^{\textbf{M}}=\left(z^M_1,\ldots,z^M_n \right)$ , $\textbf{a}=\left(a_1,\ldots,a_n \right)$ , $\alpha=\left(\alpha_1,\alpha_2 \right)$ , $\beta=\left(\beta_1,\beta_2 \right)$ , $\psi=\left(\psi_1,\ldots,\psi_{K-1} \right)$ , $\gamma^*=\left(\gamma^*_{1,1},\gamma^*_{1,2},\ldots,\gamma^*_{K,2}\right)$ , $\zeta^{A*}=\left(\zeta^{A*}_1,\ldots,\zeta^{A*}_K\right)$ and $\zeta^{M*}=\left(\zeta^{M*}_1,\ldots,\zeta^{M*}_K\right)$ . We introduce the latent indicator variable $s_{i,k}$ , which is equal to 1 if the ith couple belongs to the kth class, and zero otherwise. This auxiliary variable facilitates an efficient computation of the posterior distribution (Müller et al., Reference Müller, Erkanli and West1996).

The likelihood function of the parameters conditional on $\textbf{t}$ , $\textbf{x}$ , $\textbf{a}$ , $\textbf{d}$ , $\textbf{z}^{\textbf{A}}$ , $\textbf{z}^{\textbf{M}}$ and the latent allocation variable $\textbf{s}=\left(s_{1,1},\ldots,s_{n,K}\right)$ is given by

(5.1) \begin{align} & L \left(\alpha, \beta,\gamma^*, \zeta^{A*},\zeta^{M*}, \sigma^2_{A}, \psi \mid \textbf{t}, \textbf{x},\textbf{d},\textbf{a}, \textbf{z}^{\textbf{A}}, \textbf{z}^{\textbf{M}}, \textbf{s} \right)\nonumber \\[5pt] &\propto \displaystyle\prod_{i=1}^n \displaystyle\prod_{k=1}^K \pi_k^{s_{i,k}} \Bigg\lbrace \left[\displaystyle\prod_{j=1}^2 f \left(t_{i,j} \mid x_{i,j}, d_{i,j}, a_i;\;\alpha_j,\beta_j,\gamma^*_{k,j} \right)\right] f\left(z^A_i \mid \zeta^{A*}_{k}, \sigma^2_{A} \right) \left({\zeta^{M*}_k}\right)^{z^M_i} \left(1 - {\zeta^{M*}_k}\right)^{1 - z^M_i} \Bigg\rbrace^{s_{i,k}}\end{align}

where $f\left(\cdot \mid \zeta^{A*}_{k}, \sigma^2_{A} \right)$ denotes the probability density function of the normal distribution with mean $\zeta^{A*}_{k}$ and variance $\sigma^2_{A}$ , and $\left({\zeta^{M*}_k}\right)^{z^M_i} \left(1 - {\zeta^{M*}_k}\right)^{1 - z^M_i}$ is the probability mass function of the Bernoulli-distributed random variable $Z^M_i$ with parameter ${\zeta^{M*}_k}$ . The jth couple member likelihood contribution $f \left(t_{i,j} \mid x_{i,j}, d_{i,j}, a_i;\;\alpha_j,\beta_j,\gamma^*_{k,j} \right)$ is derived to account for right censoring and left truncation by simple algebra as follows:

(5.2) \begin{align}\nonumber f &\left(t_{i,j} \mid x_{i,j},d_{i,j},a_{i};\;\alpha_j, \beta_j,\gamma^*_{k,j} \right) &&\\[5pt] &= \frac{\exp\left[ -\int_0^{t_{i,j}+a_i} \mu \left(s \mid x_{i,j}; \alpha_j, \beta_j,\gamma_{k,j}^*\right) \mathrm{d}s \right] \mu \left(t_{i,j}+a_i \mid x_{i,j};\;\alpha_j, \beta_j,\gamma_{k,j}^*\right)^{d_{i,j}}}{\exp\left[ -\int_0^{a_i} \mu \left(q \mid x_{i,j};\;\alpha_j, \beta_j,\gamma_{k,j}^*\right) \mathrm{d}q \right]}\end{align}

where the numerator is the likelihood contribution for a noninformative right-censored observation, which is divided by the probability of being alive between policy inception (time 0) until the start of the observation period. Given the specification of the hazard function in Equation (4.1), the logarithm of the likelihood contribution is given by

(5.3) \begin{align}\nonumber \ln\, & f \left(t_{i,j} \mid x_{i,j},d_{i,j},a_{i};\;\alpha_j, \beta_j,\gamma^*_{k,j} \right) &&\\[5pt] \nonumber &= -\frac{\exp\left[ \beta_j \left(t_{i,j} + a_i \right)\right] - \exp\left(\beta_j a_i \right)}{\beta_j} \exp \left(\alpha_j + \beta_j \left(x_{i,j}-\bar{x}\right) + \gamma^*_{k,j} \right)&&\\[5pt] & \quad + d_{i,j} \left(\alpha_j + \beta_j \left(x_{i,j} - \bar{x} + t_{i,j} + a_i\right) + \gamma^*_{k,j} \right) \,.\end{align}

5.3 Posterior distribution

The (unnormalized) posterior distribution follows as the product of the likelihood and the prior distribution. In this latter, we assume that all parameters are pairwise independently distributed, with densities generically denoted by p:

(5.4)

In order to efficiently learn the posterior distribution of Equation (5.4), we propose to first data augment the dataset of the missing value of the latent class $s_{i,k}$ and then use a blocked Gibbs sampler scheme (Ishwaran and James, Reference Ishwaran and James2001), consisting of a sequential draws of the parameters (exploiting their conditional conjugacy where possible). The steps of this data augmentation-blocked MCMC sampler are detailed in Appendix A. We implement this sampler in R (R Core Team, 2013) in order to have a full control over the MCMC sampling process. The code implementing the sampler is available at the GitHub repository https://github.com/ungolof/AVDPM.

7. Actuarial illustration of AVDPM

Let $Y_{x_1,x_2}\left(z\right)$ denote the present value of a cash flow of $1 paid continuously to a couple with characteristics z, where the male is aged $x_1$ and the female $x_2$ , as long as both are alive (joint status). Conversely, let $Y_{\overline{x_1,x_2}}\left(z\right)$ denote the present value of a $1 cash flow paid continuously until the death of the last survivor of a similar couple.

Figure 4. Last survivor annuity factor calculated using the AVDPM model (solid line), the PH model (dashed line), and the BG model (dotted line) for different values of age difference ( $\exp \left( Z^A \right)$ , x-axis) and $Z^M$ when the oldest member is aged 60 (top panel) and 70 (bottom panel), for $\iota=1\%$ (left panel) and $\iota=5\%$ (right panel).

Assuming a force of interest $\iota$ , we can obtain the annuity factor of these cash flows as the expected value of $Y_{x_1,x_2}\left(z\right)$ and $Y_{\overline{x_1,x_2}}\left(z\right)$ (Dickson et al., Reference Dickson, Hardy and Waters2013):

\begin{align}\nonumber \overline{a}_{x_1,x_2}\left(z\right) &= \mathbb{E}\left[Y_{x_1,x_2}\left(z\right) \right] = \int_0^{\infty} \exp \left(- \iota t \right) S_{x_1,x_2} \left(t \mid z \right) \mathrm{d}t; &&\\[5pt] \nonumber \overline{a}_{\overline{x_1,x_2}}\left(z\right) &= \mathbb{E}\left[Y_{\overline{x_1,x_2}}\left(z\right) \right]=\int_0^{\infty} \exp \left(- \iota t \right) S_{\overline{x_1,x_2}} \left(t \mid z \right) \mathrm{d}t,\end{align}

where the last survivor function $S_{\overline{x_1,x_2}} \left(t \mid z \right)$ was defined in Section 3, and $S_{x_1,x_2} \left(t \mid z \right)$ denotes the corresponding joint survivor function.

We analyze the effect of the covariates by comparing the annuity value under the AVDPM with the annuity value obtainable using the BG model in order to assess the effect of the covariates (as well as the dependence), and with the value obtained using the PH model to assess the effect of the dependence.

Figure 4 shows the value of the last survivor annuity factor obtainable under the AVDPM model and the other two competing models PH and BG described in Section 6. The annuity factors are evaluated at different male ages (60 and 70), different values of $Z=\left(Z^A,Z^M\right),$ and two different forces of interest $\iota=(0.01, 0.05)$ . The same plot for the joint life annuity is shown in Section 4 of the Supplementary Material (Figure 4.1).

The plot shows that the AVDPM annuity value is always higher than the value obtained using the PH model. When the male is the oldest couple member, the difference between these two values is highest when the age difference is within five years and narrows afterward. Such annuity price difference is also negligible for an age difference of one year. This means that neglecting the effect of the dependence has the effect of underpricing the last survivor annuity, especially when the age difference is between two and five years. A similar evidence is obtained when the female is the oldest couple member, whenever the male is lesser than ten years younger. This evidence is in contrast with earlier findings in the literature which use copula models to account for lifetime dependence. Using a Frank copula, Frees et al. (Reference Frees, Carriere and Valdez1996, Figure 3) show that the price of a last survivor annuity is higher under a model that assumes independent lifetimes, especially for couples where males and females have a similar age. Similarly, Dufresne et al. (Reference Dufresne, Hashorva, Ratovomirija and Toukourou2018, Figure 5.2) observe that the difference between independent and dependent lifetime (using a Clayton and a Joe copula) in the last-survivor life expectancy is highest for an age difference around 0, regardless the value of $Z^M$ .

When comparing the annuity factor using the AVDPM with the value obtained using a Base Gompertz model which assumes independent lifetimes and does not account for the effect of the covariates, we note that using latter in annuity pricing yields a lower value compared to the former when the female is the older couple member. Conversely, when $Z^M=1$ , the BG model underprices the last survivor annuity for an age difference lower than eight years. These evidences are consistent with the findings in Frees et al. (Reference Frees, Carriere and Valdez1996, Figure 4) and Deresa et al. (Reference Deresa, Van Keilegom and Antonio2022), where the latter analyze only the case of a couple where the male is aged 65 years and the female is two years younger.

A higher interest rate, which decreases the value of the annuity all else being equal, yields a reduction of both the effect of the dependence and the effect of the price difference due to the covariates. On the other hand, an increase in the age of the two heads (Figure 4 bottom panel) yields a higher percentage difference between the annuity value priced using the AVDPM model and the two models which assume independence. Therefore, in this latter case, the dependence turns out to have a higher impact in the pricing of the last survivor annuity.

For the joint life annuity, we note that when the male is the oldest couple member, the joint life annuity factor under the AVDPM is always higher compared to the value obtainable under the two models assuming independent lifetimes for an age difference higher than a year, and the same is obtained in case $Z^M=0$ and the age difference is between 1 and 15 years. This is obtained for both male ages and for the two interest rates hereby analyzed. Similarly, the copula-based dependence model analyzed by Deresa et al. (Reference Deresa, Van Keilegom and Antonio2022) yields a higher joint life annuity value compared to the independence case, although they focus on the sole case where males and females are again aged 65 and 63 years, respectively.

8. Extension of AVDPM to the analysis of the joint lifetimes of nonexchangeable units

So far, we have illustrated the method for the case of male–female couples. If we want to allow for groups with differing number of exchangeable members, as can be the case of collective insurance policies, we can extend the framework through a hierarchical model with additional layers.

Suppose the ith group includes $J_i$ members, with common set of variables $z_i$ , and individual (within group) characteristics $x_{i,j}$ ( $j=1,\ldots,J_i$ ). A possible solution is to model the group-specific joint distribution of the lifetimes and Z as:

(8.1) \begin{align} f\left(t_{i,1},\ldots,t_{i,J_i},z_i\mid x_i;\;\beta,G \right) = \int_{\Omega_{\gamma,\eta}} \left[\displaystyle\prod_{j=1}^{J_i} f\left(t_{i,j} \mid x_{i,j} ;\;\beta,\gamma_{i,j} \right) \right] f \left(z_i\mid\zeta_i \right) \mathrm{d}Q \left(\gamma_{i,j}\mid \eta_i \right) \mathrm{d}G \left(\eta_i, \zeta_i \right)\end{align}

where Q denotes a suitable distribution function for $\gamma_{i,j}$ , indexed by the group-specific parameter $\eta_i$ . Again, Q can also be a random draw from a Dirichlet process, although we would opt for a simpler known parametric form for computational reasons and also because we are indexing Q with a group-specific parameter $\eta_i$ whose distribution is drawn from a DP.