2.1 Types of Dependence

Before describing our model, we formalize the notion of short-term dependence, as in Spreeuw (Reference Spreeuw2006).

We consider two lives (x) and (y), who are respectively aged x and y at duration 0. The complete remaining lifetimes of (x) and (y) are denoted by T_x and T_y, respectively. We assume that T_x and T_y are continuously distributed, with upper bounds ω_x−x and ω_y−y, respectively. The variables ω_x and ω_y denote the limiting ages of (x) and (y). For $$$ t\: \in \:\left[ {0,{{\omega }_x}{\rm{ - }}x} \right) $$$ and $$$ s\: \in \:\left[ {0,{{\omega }_y}{\rm{ - }}y} \right) $$$, we define μ ₁(x + t) and μ ₂(y + s) as the forces of mortality relating to T_x and T_y.

Further, define $$$ {{\mu }_1}\left( {x\: + \:t\left| {{{T}_y}\: = \:{{t}_y}} \right.} \right) $$$ as the conditional force of mortality of (x) at duration t (age x + t) given that (y) has died at duration t_y (age y + t_y) with $$$ {{t}_y}\: \in \:\left[ {0,t} \right) $$$. Likewise, we define $$$ {{\mu }_2}\left( {y\: + \:s\left| {{{T}_x}\: = \:{{t}_x}} \right.} \right) $$$ as the conditional force of mortality of (y) at duration s (age x + s) given that (x) has died at duration t_x (age x + t_x) with $$$ {{t}_x}\: \in \:\left[ {0,s} \right) $$$.

We can now specify the notions of long-term and short-term dependence, using the definition in Spreeuw (Reference Spreeuw2006).

2.2 Model A: Independence of Lifetimes

The first model that we consider is the standard model in multiple life contingencies where the remaining lifetimes of joint lives are independent. We refer to this subsequently as Model A. It is easiest to describe in terms of the four-state continuous-time Markov model of Norberg (Reference Norberg1989) and Wolthuis (Reference Wolthuis2003), for the mortality status of a couple consisting of a man aged x and a woman aged y. This is depicted in Figure 1. In this model, μ ₀₁(·) and μ ₂₃(·) are the force of mortality functions for a man whose spouse is still alive and for a man whose spouse has died, respectively. Likewise, μ ₀₂(·) and μ ₁₃(·) represent the respective force of mortality functions for a woman whose spouse is still alive and a woman whose spouse has died.

Figure 1 Markov model as in Norberg (Reference Norberg1989) and Wolthuis (Reference Wolthuis2003)

Independence of remaining lifetimes, in this model, implies that the force of mortality functions μ ₀₁(·) and μ ₂₃(·) for the male are identical and, likewise for the female, μ ₀₂(·) and μ ₁₃(·) are identical functions (Norberg, Reference Norberg1989).

2.3 Model B: Long-term Dependence of Lifetimes

The second model that we consider is that of Denuit & Cornet (Reference Denuit and Cornet1999), which is itself an adaptation of the model of Norberg (Reference Norberg1989) and Wolthuis (Reference Wolthuis2003) as illustrated in Figure 1. Denuit & Cornet (Reference Denuit and Cornet1999) specify the following conditional forces of mortality as functions of the marginal forces of mortality:

$$\mu _{{01}}^{ \ast } \left( t \right)\: = \:{{\mu }_1}\left( {x\: + \:t\left| {{{T}_y}\: \gt \:t} \right.} \right)\: = \:\left( {1{\rm{ - }}\alpha _{{01}}^{ \ast } } \right)\:{{\mu }_1}\left( {x\: + \:t} \right)$$

$$ \mu _{{02}}^{ \ast } \left( t \right)\: = \:{{\mu }_2}\left( {y\: + \:t\left| {{{T}_x}\: \gt \:t} \right.} \right)\: = \:\left( {1{\rm{ - }}\alpha _{{02}}^{ \ast } } \right)\:{{\mu }_2}\left( {y\: + \:t} \right) $$

$$\mu _{{13}}^{ \ast } \left( t \right)\: = \:{{\mu }_2}\left( {y\: + \:t\left| {{{T}_x}\:\leq \:t} \right.} \right)\: = \:\left( {1\: + \:\alpha _{{13}}^{ \ast } } \right)\:{{\mu }_2}\left( {y\: + \:t} \right)$$

$$ \mu _{{23}}^{ \ast } \left( t \right)\: = \:{{\mu }_1}\left( {x\: + \:t\left| {{{T}_y}\:\leq \:t} \right.} \right)\: = \:\left( {1\: + \:\alpha _{{23}}^{ \ast } } \right)\:{{\mu }_1}\left( {x\: + \:t} \right) $$

where $$$ \alpha _{{01}}^{ \ast }, \,\alpha _{{02}}^{ \ast }, \,\alpha _{{13}}^{ \ast }, \,\alpha _{{23}}^{ \ast } \:\geq \:0 $$$. (We reserve the non-starred version of these symbols for our main model, which is Model C below.)

In this model, death of the man leads to a constant increase of the woman's mortality by $$$ \frac{{1\: + \:\alpha _{{13}}^{ \ast } }}{{1{\rm{ - }}\alpha _{{02}}^{ \ast } }} $$$, whereas a man's mortality goes up by a factor of $$$ \frac{{1\: + \:\alpha _{{23}}^{ \ast } }}{{1{\rm{ - }}\alpha _{{01}}^{ \ast } }} $$$ whenever his spouse dies. Note that this model is a special case of long-term dependence, since the mortality of one life only depends on whether the spouse has died or not, and not when (s)he died.

2.4 Model C: Short-term Dependence of Lifetimes

We may now introduce our model, which we label as Model C. We extend the Markov model of Denuit & Cornet (Reference Denuit and Cornet1999) by allowing the mortality of a remaining life to depend on the time elapsed since spouse's death. Upon the death of his spouse, every widower enters an initial bereaved state. He may leave this initial bereaved state either by transition to the death state at any time or by transition to an ultimate widowed state after a fixed period. Similar states of initial bereavement and ultimate widowhood exist for widows. The crucial point is that the survivor's mortality in the ultimate widowed state can differ from that in the initial bereaved state. This model therefore allows explicitly for short-term dependence as per Definition 1.

More precisely, we augment the state space of the four-state model in Figure 1 by splitting each of the widowed states into two further states. This leads to a six-state model as shown in Figure 2. A woman becoming widowed will enter state 1 ((y) alive, (x) died less than t ₁ years ago) in which she will stay for at most t ₁ years, after which, if still alive, she makes the transition to state 2 ((y) alive, (x) died more than t ₁ years ago) from which only the transition to state 5 (both dead) is possible. Likewise, a man losing his spouse will first enter state 3 ((x) alive, (y) died less than t ₂ years ago) where he will stay for at most t ₂ years, after which he will automatically make the transition to state 4 and stay there while alive. Note that t ₁ and t ₂ are not necessarily equal. This allows for a different time-scale of broken-heart effect for males and females.

Figure 2 Six-state model

Using the extended model in a proportional hazards setting requires additional parameters than those in equations (1a)–(1d). The modified specification is as follows:

$${{\mu }_{01}}\left( t \right)\: = \:{{\mu }_1}\left( {x\: + \:t\left| {{{T}_y}\: \gt \:t} \right.} \right)\: = \:\left( {1{\rm{ - }}{{\alpha }_{01}}} \right)\:{{\mu }_1}\left( {x\: + \:t} \right)$$

$${{\mu }_{03}}\left( t \right)\: = \:{{\mu }_2}\left( {y\: + \:t\left| {{{T}_x}\: \gt \:t} \right.} \right)\: = \:\left( {1{\rm{ - }}{{\alpha }_{03}}} \right)\:{{\mu }_2}\left( {y\: + \:t} \right)$$

$$ {{\mu }_{15}}\left( t \right)\: = \:{{\mu }_2}\left( {y\: + \:t\left| {0\:\leq \:t{\rm{ - }}{{T}_x}\: \lt \:{{t}_1}} \right.} \right)\: = \:\left( {1\: + \:{{\alpha }_{15}}} \right)\:{{\mu }_2}\left( {y\: + \:t} \right) $$

$$ {{\mu }_{25}}\left( t \right)\: = \:{{\mu }_2}\left( {y\: + \:t\left| {t{\rm{ - }}{{T}_x}\:\geq \:{{t}_1}} \right.} \right)\: = \:\left( {1\: + \:{{\alpha }_{25}}} \right)\:{{\mu }_2}\left( {y\: + \:t} \right) $$

$$ {{\mu }_{35}}\left( t \right)\: = \:{{\mu }_1}\left( {x\: + \:t\left| {0\:\leq \:t{\rm{ - }}{{T}_y}\: \lt \:{{t}_2}} \right.} \right)\: = \:\left( {1\: + \:{{\alpha }_{35}}} \right)\:{{\mu }_1}\left( {x\: + \:t} \right) $$

$$ {{\mu }_{45}}\left( t \right)\, = \:{{\mu }_1}\left( {x\: + \:t\left| {t{\rm{ - }}{{T}_y}\,\geq \,{{t}_2}} \right.} \right)\: = \:\left( {1\: + \:{{\alpha }_{45}}} \right)\:{{\mu }_1}\left( {x\: + \:t} \right) $$

where $$$ \alpha _{{01}}^{{}}, \,\alpha _{{03}}^{{}}, \,\alpha _{{15}}^{{}}, \,\alpha _{{25}}^{{}}, \,\alpha _{{35}}^{{}}, \,{{\alpha }_{45}}\:\geq \:0 $$$.

Note that the earlier four-state Model B is a special case of our augmented six-state Model C, with α ₁₅ = α ₂₅ = $$$\alpha _{{13}}^{ \ast } $$$ and α ₃₅ = α ₄₅ = $$$\alpha _{{23}}^{ \ast } $$$.

Note also that it is conceivable for short-term dependence to be negative, in the sense that α ₁₅ or α ₃₅ could be negative. For example, the strain of caring for a sick partner could be relieved upon the partner's death. However, this effect is not reported in the literature and, overall, widow(er)hood appears to increase mortality initially. The data set that we describe in the next section also confirms this and gives no evidence of negative short-term dependence.

4.1 Model Identification

Our aim in this section is to estimate our augmented Markov model, as described in section 2.4, using the data set described in section 3. The first step is to give the precise specification of Model C, by identifying the cut-off point between states 1 and 2, and between states 3 and 4.

In particular, we choose t ₁ by testing, for different values of t ₁, whether there is a significant difference between the observed mortality rates of recently bereaved widows (where spouse's death occurs within t ₁ years) and the observed mortality rates of the remaining widows. Our test hypotheses are:

$${{H}_0}:\; \; S_{y}^{{15}} \left( t \right)\: = \:S_{y}^{{25}} \left( t \right),\; \; {\rm{for}}\,{\rm{all}}\,\,t\:{\rm{ \gt }}\:{\rm{0,}} $$

$$ {{H}_1}:\; \; S_{y}^{{15}} \left( t \right)\: \ne \:S_{y}^{{25}} \left( t \right),\; \; {\rm{for}}\,{\rm{at}}\,{\rm{least}}\,{\rm{one}}\,\,t\,{\rm{ \gt }}\,{\rm{0,}} $$

where $$$ S_{y}^{{i5}} \left( t \right) $$$ is the probability that a widow, aged y in state i where i ∈ {1, 2}, does not enter state 5 within t years. If H ₀ is true, then the observed mortality rates in states 1 and 2 are two samples from the same survival function, otherwise they are governed by different survival functions.

To compare the mortality of recent widows with longer-term widows, we perform a two-sample Kolmogorov-Smirnov test. In the usual version of this test, the distributions of two samples are compared by computing the maximal absolute deviation $$$ {{D}_{m,n}}\: = \:{\sup }_t \left| {{{{\mathop\^{F}\limits}}_m}(t){\rm{ - }}{{{\mathop\^{G}\limits}}_n}(t)} \right| $$$, where m and n are the sizes of the two samples and $$$ {{\^{F} }_m}(t) $$$ and $$$ {{\^{G} }_n}(t) $$$ are their respective empirical distribution functions. The null hypothesis that $$$ {{\^{F} }_m}(t)\, = \,{{\^{G} }_n}(t) $$$ for t > 0 is rejected if the test statistic $$$ {{\left( {\frac{{mn}}{{m\: + \:n}}} \right)}^{\frac{1}{2}}} {{D}_{m,n}} $$$ is greater than the critical value of the Kolmogorov distribution for a certain significance level. We use the censored-data version of the Kolmogorov-Smirnov test which is based on empirical survival functions. More precisely, the test employs a statistic that is the largest absolute value of the weighted sum of the differences in the Nelson-Aalen estimates of cumulative hazard rates in the two samples. This test is readily implemented in various statistical software packages and is described by Fleming et al. (Reference Fleming, O'Fallon, O'Brien and Harrington1980), with the asymptotic theory discussed by Andersen et al. (Reference Andersen, Borgan, Gill and Keiding1993, p. 391).

An alternative to the Kolmogorov-Smirnov test is the logrank test (see e.g. Machin et al., Reference Machin, Cheung and Parmar2006, p. 51). However, this does not employ Smirnov-type maximal absolute deviation statistics. It may therefore fail to capture temporary differences in survival distributions when, for example, a large local deviation is offset by a small opposite overall bias (Fleming et al., Reference Fleming, O'Fallon, O'Brien and Harrington1980; Klein & Moeschberger, Reference Klein and Moeschberger1997, p. 209). It is, of course, temporary differences in mortality between recently bereaved and longer-term widow(er)s that we seek to identify.

The results of the Kolmogorov-Smirnov test comparing the mortality of widows aged y = 60 in states 1 and 2, for different values of t ₁, are shown in Table 4. We note that the higher t ₁ is, the higher the p-value. For t ₁ = 0.5 or t ₁ = 1, the test is significant at 5% level, but not for higher values of t ₁. We conclude that both t ₁ = 0.5 and t ₁ = 1 are suitable as cut-off points between the states 1 and 2. For computational convenience, we choose t ₁ = 1 in the numerical applications that follow. A similar test can be undertaken for widower's mortality, of course. We obtain similar results and again t ₂ = 1 is a convenient specification. This is also consistent with the results of Parkes et al. (Reference Parkes, Benjamin and Fitzgerald1969) and Jagger & Sutton (Reference Jagger and Sutton1991) for widowers.

Table 4 Outcomes of Kolmogorov-Smirnov test

Before proceeding to parameterize the six-state model, we make two further remarks. First, the Kolmogorov-Smirnov test is a nonparametric test and the choice of t ₁ = 1 and t ₂ = 1 does not depend on any assumed underlying mortality law. Secondly, the test is used here to identify the six-state model and provide indicative values for t ₁ and t ₂. The model is parameterized and then validated using further tests later.

4.2 Estimation

In their copula models of joint mortality (briefly described in section 1), Frees et al. (Reference Frees, Carriere and Valdez1996) and Carriere (Reference Carriere2000) use marginal distributions that are Gompertz. Carriere (Reference Carriere2000) fits the marginal survival function prior to considering copulas and demonstrates that the Gompertz model performs significantly better than other mortality models. In the case of independent lifetimes, the copula models and the Markov models are exactly the same, of course. Ji et al. (Reference Ji, Hardy and Li2011) also cite the parsimony of the Gompertz mortality function, its smoothness, and the ease of extrapolation to extreme ages as justification for using Gompertz mortality in their semi-Markov model of joint mortality. For these reasons, we also specify the marginal force of mortality functions to be Gompertz, implying that

$$ {{\mu }_1}(x)\: = \:\frac{1}{{{{\sigma }_m}}}\exp \left( {\frac{{x{\rm{ - }}{{m}_m}}}{{{{\sigma }_m}}}} \right),\; \; \; \; {{\mu }_2}(y)\: = \:\frac{1}{{{{\sigma }_f}}}\exp \left( {\frac{{y{\rm{ - }}{{m}_f}}}{{{{\sigma }_f}}}} \right). $$

We follow the two-step estimation procedure of Denuit & Cornet (Reference Denuit and Cornet1999). First, the parameters of the base Gompertz distributions are estimated using all male and female mortality data. Based on this, the α family of parameters defined in equations (2a)–(2f) are then estimated. This is the same approach as in Denuit & Cornet (Reference Denuit and Cornet1999), except that they parameterize a Makeham distribution.

The log-likelihood in terms of the parameters of the male base Gompertz mortality function is:

$${{\ell }_1}\: = \:{{e}^{{\rm{ - }}{{m}_m}/{{\sigma }_m}}} \mathop{\sum}\limits_{i:{\rm{all}}\,{\rm{males}}} \left( {{{e}^{{{u}_i}/{{\sigma }_m}}} {\rm{ - }}{{e}^{{{v}_i}/{{\sigma }_m}}} } \right)\: + \:\frac{1}{{{{\sigma }_m}}}\mathop{\sum}\limits_{{\scriptstyle i:{\rm{male}} \atop \scriptstyle {\rm{deaths}}}} {{v}_i}\:{\rm{ - }}\:\frac{{{{m}_m}{{d}_m}}}{{{{\sigma }_m}}}\;{\rm{ - }}\;{{d}_m}\log {{\sigma }_m},$$

where d_m denotes the total number of male deaths, u_i is the entry age of (male) life i in the investigation, and v_i is his exit age on death or censoring. A similar function applies to the log-likelihood pertaining to the parameters of the female base Gompertz mortality function.

Our maximum likelihood estimates, with standard errors in brackets, of the parameters of the base Gompertz mortality functions in equation (3) are as follows: m_m = 86.37 (0.247), σ_m = 9.76 (0.343), m_f = 92.07 (0.336), and σ_f = 8.06 (0.217). There are slight discrepancies in our estimates of σ_m and σ_f as compared to Frees et al. (Reference Frees, Carriere and Valdez1996) and Carriere (Reference Carriere2000), possibly because our data set is slightly different from the one used by these authors.

An alternative to this two-step estimation procedure would be to estimate all parameters, i.e. Gompertz parameters and α-parameters, in a single maximum-likelihood exercise. This could result in somewhat better estimates of the Gompertz parameters but it would also lead to more complicated likelihood functions. In any event, the estimates of m_m, σ_m, m_f and σ_f given above have small standard errors. The two-step procedure also has the advantage that formulas for the estimators of the α parameters are simple and easy to interpret.

Next, we estimate the α family of parameters in equations (2a)–(2f) by partial maximum likelihood. The partial log-likelihood function for α ₀₁ is given by

$$ \ell _{1}^{p} ({{\alpha }_{01}})\: = \:{\rm{ - }}\left( {1{\rm{ - }}{{\alpha }_{01}}} \right)\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{males}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{0}}}} {\int}_{{{{\!\!\!\!\!u}_{0,i}}}}^{{{{v}_{0,i}}}} {{\mu }_{{{x}_i}}}\:d{{x}_i}\: + \:{{d}_{01}}\log \left( {1{\rm{ - }}{{\alpha }_{01}}} \right)\: + \:\mathop{\sum}\limits_{{\scriptstyle i:{\rm{male}}\,{\rm{deaths}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{0}}}} {{\mu }_{{{v}_{0,i}}}} $$

with d ₀₁ denoting the observed number of male deaths in state 0, u _0,i denoting the entry age of life i in state 0, and v _0,i denoting his age at exit from state 0. A similar partial log-likelihood function $$$ \ell _{2}^{p} ({{\alpha }_{02}}) $$$ for α ₀₂ may be found. For lives in state 0, this leads to estimates

$$ {{\^{\alpha } }_{01}}\: = \:1\:{\rm{ - }}\:{{d}_{01}}\:{{\left( {\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{males}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{0}}}} {\int}_{{{{\!\!\!\!u}_{0,i}}}}^{{{{v}_{0,i}}}} {{\mu }_{{{x}_i}}}d{{x}_i}} \right)}^{{\rm{ - }}1}}, $$

$$ {{\^{\alpha } }_{02}}\: = \:1\:{\rm{ - }}\:{{d}_{02}}\:{{\left( {\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{females}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{0}}}} {\int}_{{{{\!\!\!\!u}_{0,i}}}}^{{{{v}_{0,i}}}} {{\mu }_{{{y}_i}}}d{{y}_i}} \right)}^{{\rm{ - }}1}}, $$

where d ₀₁ is as above and d ₀₂ is the observed number of female deaths in state 0. We see from equations (6) above that, for each life under risk, we take the force of mortality with parameters estimated as above, integrated over the age range in which the life was under observation. The observation in state 0 ends on death of the life, death of the spouse, or at censoring.

Given the partial maximum likelihood, the parameters α ₀₁, α ₀₂, α ₁₅, α ₂₅, α ₃₅ and α ₄₅ are estimated independently of each other. For instance, if $$$ {{\tilde{\alpha }}_{01}} $$$ represents the maximum likelihood estimator of α ₀₁ and D ₀₁ represents the random number of male deaths in state 0, then $$$ {{\tilde{\alpha }}_{01}} $$$ is asymptotically normally distributed with mean α ₀₁ and variance given by the Cramér-Rao lower bound. This can be simplified to

$$ {\rm{ - }}{{\left( {{\Bbb E}\left[ {\frac{{{{\partial }^2} \ell _{1}^{p} }}{{\partial \alpha _{{01}}^{2} }}} \right]} \right)}^{{\rm{ - }}1}} \: = \:{\rm{ - }}{{\left( {{\Bbb E}\left[ {\frac{{{\rm{ - }}{{D}_{01}}}}{{{{{(1{\rm{ - }}{{\alpha }_{01}})}}^2} }}} \right]} \right)}^{{\rm{ - }}1}} \: = \:\frac{{{{{(1{\rm{ - }}{{\alpha }_{01}})}}^2} }}{{{\Bbb E}{{D}_{01}}}}. $$

The standard error of $$$ {{\^{\alpha } }_{01}} $$$ is therefore estimated as $$$ (1{\rm{ - }}{{\^{\alpha } }_{01}})/\sqrt {{{d}_{01}}} $$$, with a similar expression for the standard error of $$$ {{\^{\alpha } }_{02}} $$$.

The remaining parameters α_{j 5}, for j ∈ {1, 2, 3, 4}, are also estimated as above. For example, the partial log-likelihood function for α ₁₅ is

$$\ell _{3}^{p} ({{\alpha }_{15}})\: = \:{\rm{ - }}\left( {1\: + \:{{\alpha }_{15}}} \right)\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{females}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{1}}}} {\int}_{{{{\!\!\!\!u}_{1,i}}}}^{{{{v}_{1,i}}}} {{\mu }_{{{y}_i}}}\:d{{y}_i}\: + \:{{d}_{15}}\log \left( {1\: + \:{{\alpha }_{15}}} \right)\: + \:\mathop{\sum}\limits_{{\scriptstyle i:{\rm{female}}\,{\rm{deaths}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{1}}}} {{\mu }_{{{v}_{1,i}}}}.$$

The notation is self-explanatory, e.g. d_{j 5} represents the observed number of female deaths in state j, u_j,i is the entry age of life i in state j, and υ_j,i is the exit age of life i from state j. Thus, for lives in states 1, …, 4, we get as estimates:

$$ {{\^{\alpha } }_{j5}}\: = \:{{d}_{j5}}\:{{\left( {\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{females}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,j}} {\int}_{{{{\!\!\!\!u}_{j,i}}}}^{{{{v}_{j,i}}}} {{\mu }_{{{y}_i}}}d{{y}_i}} \right)}^{{\rm{ - }}1}} \:{\rm{ - }}\:1,\; \; \; \; {\rm{for}}\,\,j\: \in \:\left\{ {{\rm{1,2}}} \right\}{\rm{,}} $$

$$ {{\^{\alpha } }_{k5}}\: = \:{{d}_{k5}}\:{{\left( {\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{males}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,k}} {\int}_{{{{\!\!\!\!u}_{k,i}}}}^{{{{v}_{k,i}}}} {{\mu }_{{{x}_i}}}d{{x}_i}} \right)}^{{\rm{ - }}1}} \:{\rm{ - }}\:1,\; \; \; \; {\rm{for}}\,k\: \in \:\left\{ {{\rm{3,4}}} \right\}{\rm{.}} $$

The standard errors of the estimates in equations (9) are given by

$$ s.e.\left[ {{{{\^{\alpha } }}_{j5}}} \right]\, = \,\frac{{1\: + \:{{{\^{\alpha } }}_{j5}}}}{{\sqrt {{{d}_{j5}}} }}\; \; \; \; {\rm{for}}\,j\, \in \:\left\{ {{\rm{1,}} \ldots {\rm{,4}}} \right\}{\rm{.}} $$

The estimates turn out to be very sensitive to the choice of age range. When estimating coefficients $$$ {{\^{\alpha } }_{j5}} $$$, j∈{1, …, 4}, we use as ages of entry the intervals [65, 85] for males and [60, 80] for females, as these intervals contain the largest proportion of widowed lives. Our findings also suggest that the dependence between the two lifetimes varies with the ages of individuals.

Tables 5 and 6 display the results of our parameter estimation procedure, assuming different values of t ₁ and t ₂. The last two rows of each table display the parameters pertaining to Model B, i.e. the standard four-state Markov model as in Figure 1. Note that t ₁ = ∞ and t ₂ = ∞ imply that $$$ {{\alpha }_{15}}\: = \:\alpha _{{13}}^{ \ast } $$$ and $$$ {{\alpha }_{35}}\: = \:\alpha _{{23}}^{ \ast } $$$, and therefore the parameters α ₂₅ and α ₄₅ are irrelevant. The tables clearly show that α ₀₁ is smaller than both α ₁₅ and α ₂₅, and similarly that α ₀₃ is smaller than both α ₃₅ and α ₄₅. This implies that the lifetimes are dependent on each other. Tables 5 and 6 also show that α ₁₅ > α ₂₅ and α ₃₅ > α ₄₅. This suggests the presence of short-term dependence for both genders, with such dependence being stronger for widowers than for widows.

Table 5 Estimates of parameters of dependence for bereaved females. Numbers in brackets are relevant number of deaths in data and standard error respectively

Table 6 Estimates of parameters of dependence for bereaved males. Numbers in brackets are relevant number of deaths in data and standard error respectively

However, the estimates of the coefficients α ₁₅, … α ₄₅ all have fairly large standard errors, so the differences between α ₁₅ and α ₂₅, and also between α ₃₅ and α ₄₅, could be due to chance. We therefore perform a formal test for short-term dependence in the next section.

4.3 Model Testing

We perform two validation tests on the model. First, we test whether widows’ mortality depends significantly on the time elapsed since death of the spouse. We use the likelihood ratio test to test the null hypothesis

$${{H}_{\rm{0}}}{\rm{:}}\; \; {{\alpha }_{{\rm{15}}}}\,{\rm{ = }}\,{{\alpha }_{{\rm{25}}}}$$

against the alternative hypothesis

$${{H}_{\rm{1}}}{\rm{:}}\; \; {{\alpha }_{{\rm{15}}}}\, \ne \,{{\alpha }_{{\rm{25}}}}.$$

Let $$$ L_{3}^{p} ({{\alpha }_{15}}) $$$ be the partial likelihood function for α ₁₅, corresponding to the log-likelihood $$$ \ell _{3}^{p} ({{\alpha }_{15}}) $$$ given in equation (8). That is, $$$ \ell _{3}^{p} ({{\alpha }_{15}})\: = \:\log \left( {L_{3}^{p} ({{\alpha }_{15}})} \right) $$$. The likelihood ratio is therefore

$$\openup 6pt \displaylines{ \rLambda \, = \,\frac{{L_{3}^{p} ({{{\^{\alpha } }}_{25}})}}{{{{{\sup }}_{{{\alpha }_{15}}}}L_{3}^{p} ({{\alpha }_{15}})}}\, = \: & \frac{{L_{3}^{p} ({{{\^{\alpha } }}_{25}})}}{{L_{3}^{p} ({{{\^{\alpha } }}_{15}})}}\, = \,\exp \left( {\ell _{3}^{p} ({{{\^{\alpha } }}_{25}}){\rm{ - }}\ell _{3}^{p} ({{{\^{\alpha } }}_{15}})} \right) \cr = & \:{{\left( {\frac{{1\: + \:{{{\^{\alpha } }}_{25}}}}{{1\: + \:{{{\^{\alpha } }}_{15}}}}} \right)}^{{{d}_{15}}}} \exp \left( {({{{\^{\alpha } }}_{15}}{\rm{ - }}{{{\^{\alpha } }}_{25}})\mathop{\sum}\limits_{{\scriptstyle i:{\rm{all}}\,{\rm{females}} \atop \scriptstyle {\rm{in}}\,{\rm{state}}\,{\rm{1}}}} {\int}_{{{{\!\!\!\!u}_{1,i}}}}^{{{{v}_{1,i}}}} {{\mu }_{{{y}_i}}}dy} \right), \cr} $$

where d ₁₅, u _1,i and v _1,i are as defined earlier.

Table 7 exhibits the results of the tests for different values of t ₁. The critical value of the χ ²-distribution with 1 d.f. at 5% significance level is 3.84, and we find that H ₀ is rejected for t ₁ = 0.5, 1, 1.5 and 2. This validates our choice of t ₁ = 1 in section 4.1, although we note that t ₁ = 0.5 yields an even smaller p-value than t ₁ = 1. Similar tests for widowers’ mortality confirm that t ₂ = 1 is acceptable.

Table 7 Results of likelihood ratio tests for dependence of widows’ mortality on time since partner's death

The second validation test that we consider is a test of goodness-of-fit to the mortality data of the base Gompertz mortality distributions in equations (3), with the parameter estimates given in section 4.2. We consider males and females separately and carry out a χ ²-test, losing 2 degrees of freedom because of the Gompertz parameters. For males aged [65,90] and females aged [47,98], the null hypothesis of Gompertz mortality is not rejected at 5% significance level (p-values of 0.062 and 0.065 for males and females respectively). When we include younger and extreme older ages for males, the paucity of data and consequent small risk exposure lead to a poor fit for the Gompertz. We discuss the use of Gompertz mortality further in the next section.

4.4 Data and Modelling Issues

Before concluding this section on estimation, a number of issues concerning the data and the model that we have used are highlighted.

First, even though we have a large data set, the data is sparse at extreme ages. There is a small number of widows (906) and an even smaller number of widowers (337), which exacerbates estimation errors. This can be seen in Tables 5 and 6. For example, the standard error for α ₁₅ is larger than for α ₀₃, the former pertaining to the mortality of recently bereaved widows and the latter to females whose partners are still alive. The standard error for α ₁₅ is however smaller than for α ₃₅, where α ₃₅ concerns the mortality of recently bereaved widowers.

Secondly, we recognise that the Gompertz mortality law is idealized and is chosen for its tractability and, although the test in section 4.3 shows adequate goodness-of-fit, it is unlikely to be the best-fitting mortality model. Indeed, Ji et al. (Reference Ji, Hardy and Li2011) make a similar argument and justify their use of Gompertz mortality on similar grounds. Nevertheless, it is worth noting that Carriere (Reference Carriere2000) finds that alternative distributions, such as the Weibull, have a worse fit than the Gompertz to the same data set. This is also consistent with the findings of Ji et al. (Reference Ji, Hardy and Li2011) on a four-state model that corresponds to Model B in section 2.3 but also includes a constant-intensity transition for the simultaneous death of both partners, in the event of a common shock such as an accident. They calibrate the Gompertz distribution to each of the 4 states separately (whereas we followed Denuit & Cornet (Reference Denuit and Cornet1999) by fitting the Gompertz distribution to all male and all female mortality data) and they find that Gompertz’ law gives an adequate fit to the same data set, although Makeham's law may provide an improved fit.

Thirdly, it is worth remarking on a specific drawback of our 6-state Markov model (Model C), which is the jump in mortality experienced upon transition from state 1 to 2, and from state 3 to 4. It is unrealistic that there would be such an abrupt end to the broken-heart effect, with a step change in mortality after an initial bereavement period. This very issue is addressed by Ji et al. (Reference Ji, Hardy and Li2011) in work that was done concurrently with ours and on the same data set. They use a semi-Markov model, which allows the mortality of the widow(er) to be a continuous function of time elapsed since death of the spouse.

This 4-state semi-Markov model is an attractive alternative to our 6-state Markov model, but it also suffers from a number of shortcomings. For example, Ji et al. (Reference Ji, Hardy and Li2011) use a transition intensity that is exponentially declining in time since the partner's death. It does not seem very realistic to assume that the broken-heart effect declines at its fastest rate at the very time when the partner dies. The broken-heart syndrome may arguably intensify for a period shortly after the death of a spouse as the effects of loneliness and other factors cumulate with adverse psychological consequences (Holden et al., Reference Holden, Kim and Novak2010). It may be more appropriate to consider a reverse sigmoidal function, representing stagnation or very slow initial decline followed by a faster amelioration in mortality rates as recovery from bereavement progresses.

An advantage of our model is that the jump in mortality is like a step function that may indeed approximate the reverse sigmoidal shape that a post-bereavement excess mortality function might take. The broken-heart effect is likely to taper off over a much shorter time-scale than the duration of long-term life insurance products. Indeed, the analysis of Ji et al. (Reference Ji, Hardy and Li2011, Fig. 5) seems to indicate that, after one year, the excess post-bereavement mortality has almost halved for widows and completely disappeared for widowers, which is consistent with our choice of t ₁ = 1 year in the statistical test described in section 4.3. Consequently, the smoothness of such a rapid decline may have little bearing on the pricing of long-term annuities and assurances.

A major constraint of semi-Markov models is that parameterizing them often requires considerably more data than is available. For example, in the model of Ji et al. (Reference Ji, Hardy and Li2011), a parameter is required in the exponential decay function to govern the speed of reduction of the bereavement effect. As is acknowledged by Ji et al. (Reference Ji, Hardy and Li2011), the small number of deaths among widows and widowers leads to high parameter uncertainty. Ignoring catastrophic events such as accidents that lead to simultaneous deaths, there is an implicit assumption in their model that the longer-term mortality of a widowed individual is the same as it would be if the individual's partner were alive. (That is, the exponential decay function specified by Ji et al. (Reference Ji, Hardy and Li2011) tends to 1.) This is not necessarily the case in Model C, where it is possible for μ ₂₅ ≠ μ ₀₃, for example. One may make allowance for this in the semi-Markov model, but only at the expense of an additional parameter for each gender, thereby compounding parameter estimation error. Another advantage of our model is therefore that it has a minimum number of parameters, allowing for both an immediate bereavement effect (α ₁₅, α ₃₅) and a possible long-term bereavement effect after the initial period (α ₂₅, α ₄₅).

Finally, the jump in transition intensities is in fact a common feature of all Markov models that involve splitting of states. It occurs for example in demographic models for disability insurance, which are useful for estimating rates of recovery from disability (Haberman & Pitacco, Reference Haberman and Pitacco1999; Gregorius, Reference Gregorius1993), and in Markov chain models of select mortality (Norberg, Reference Norberg1988; Möller, Reference Möller1990).

5.1 Model Implementation

In this section, we use the models and parameter estimates previously obtained to analyze the impact of the type of dependence on the pricing and valuation of contracts which may still be in force when one of the lives dies. We investigate the following whole-life contracts:

• • Contingent assurance contracts: a benefit of 1 is payable immediately on the death of (y), provided this happens after the death of (x). We consider three distinct premium payment arrangements: (a) Single Premium, (b) Level Premium I, which refers to a level premium payment annually in advance while both alive, and (c) Level Premium II, which refers to a level premium payment annually in advance while (y) is alive. We assume that, at the issue of the contract, ages are x = 55 and y = 50.

• • Joint life and survivor annuity: a benefit of 1 p.a. is payable in arrears until either (x) or (y) dies, reducing to 0.6 p.a. after the first death and continuing while the survivor is alive. We consider a single premium payment only. We assume that, at the issue of the contract, both lives are aged 65.

The data set, that we describe in section 3 and that we use to parameterize our model in section 4, is based on annuities. Annuitant mortality is typically lighter than assurance mortality. Whilst it is fine to investigate the joint life and survivor annuity above, it is not ideal to model a contingent assurance contract based ultimately on annuitant mortality data. Nevertheless, this is the only data set that is available at present. We believe that it is valuable to consider at the very least the qualitative effects of short-term dependence in joint mortality on assurances as well as on annuities.

For both types of contracts, we calculate:

1. 1. the premium payable, under different premium payment arrangements,

2. 2. the state-dependent net premium provisions at several durations.

The term provision, commonly used in UK actuarial practice, may be interpreted synonymously with the term reserve. Note that when calculating provisions in state 0, i.e. when both lives are alive, we calculate the value at 0, 1, 5, 10, 15, 20 years since the start of the contract. When calculating provisions in states 1, 2, 3 or 4, i.e. when one spouse has died and the other is still alive, the time of death of the first spouse can take the values 15, 19.25, 19.5, 19.75 and 20 years (all from the start of the contract) while the provision is calculated at duration 20.

We use all three models, as described in section 2, to perform the above calculations. The models and their associated parameter values are summarized below.

• • Model A is the Markov model with independence of the remaining lifetimes T_x and T_y. We use the Gompertz mortality functions of equation (3) in section 4.2 with estimated parameters given therein.

• • Model B is the four-state model of Denuit & Cornet (Reference Denuit and Cornet1999), as displayed in Figure 1. The mortality of widow(er)s does not depend on time elapsed since spouse's death. This implies that $$$ \alpha _{{13}}^{ \ast } \: = \:{{\alpha }_{15}}\: = \:{{\alpha }_{25}} $$$ and $$$ \alpha _{{23}}^{ \ast } \: = \:{{\alpha }_{35}}\: = \:{{\alpha }_{45}} $$$. Using the data set gives parameter estimates $$$ \^{\alpha } _{{01}}^{ \ast } \: = \:0.06,\,\^{\alpha }\limits_{{02}}^{ \ast } \: = \:0.14,\,\^{\alpha }_{{13}}^{ \ast } \: = \:2.01,\,\^{\alpha } _{{23}}^{ \ast } \: = \:2.93 $$$ (from Table 5 with t ₁ = ∞, and Table 6 with t ₂ = ∞).

• • Model C is our extended Markov model, displayed in Figure 2. It allows the mortality of widow(er)s to depend on time elapsed since death of the spouse. We use t ₁ = t ₂ = 1 and the parameter values $$$ {{\^{\alpha } }_{01}}\: = \:0.06,\,{{\^{\alpha } }_{03}}\: = \:0.14 $$$ (which are obviously the same as in Model B for lives with partner alive) as well as the parameter values $$$ {{\^{\alpha } }_{15}}\: = \:3.40,\,{{\^{\alpha } }_{25}}\: = \:1.15,\,{{\^{\alpha } }_{35}}\: = \:7.19\ {\rm{and}}\ {{\^{\alpha } }_{45}}\: = \:0.41 $$$ (from Table 5 with t ₁ = 1, and Table 6 with t ₂ = 1).

For both sets of policies, we assume interest at 5% per annum.

5.2 Contingent Assurance Contracts

5.2.1 Premiums

Table 8 shows the premiums calculated for the contingent assurance contract under the three different models and for the three different premium payment patterns. Model B results in the highest premiums, followed by Model C and Model A. Dependence causes premiums for contingent assurances to be higher which is why Model A yields the lowest premiums.

Table 8 Premiums for the contingent assurance contracts

Thus, ignoring dependence results in under-pricing by around 20% (Model A relative to either Model B or Model C). But assuming dependence that is persistent rather than short-term results in over-pricing by around 6% (Model B relative to Model C).

The lower premium for Model C relative to Model B is due to the lower mortality experienced by widows after a year in Model C. This effect appears to outweigh the impact of higher mortality in the first year upon entering widowhood.

5.2.2 Provisions

Table 9 gives the provisions in state 0 for contingent assurances. Note again that Model B gives the highest values, and Model A the lowest values. Our modelling shows that an insurance company which assumes independence of joint lifetimes, when they are in fact dependent in the short term, will under-provision for a contingent assurance by about 20% (Model A compared to Model C). And if it assumes long-term dependence, when dependence is in fact short-term, it will over-provision by around 4–6%. This remains true irrespective of the premium payment pattern.

Table 9 Provisions in state 0 for the contingent assurance contracts

In Table 10, we can view the provisions to be held at duration 20 when (x) has died and (y) is still alive. Obviously, death of the male causes the payment of the sum assured to be made with certainty. This is why the provisions in Table 10 are significantly higher than the values in Table 9. Again, the independence assumption and long-term dependence assumption lead to significant under-provisioning and over-provisioning respectively, compared to short-term dependence.

Table 10 Provisions at duration 20 in state 1 (Models A and B) and in state 1 or 2 (Model C) for the contingent assurance contracts

The provision at duration 20 has been tabulated according to different periods elapsed since death of the spouse, i.e. spouse died 1 year ago (duration 19), 9 months ago (duration 19.25), half a year ago (duration 19.5), 3 months ago (duration 19.75), or has just died (duration 20). Model B takes no account of the time elapsed since spouse's death and the provision under Model B is therefore fixed at 0.5784. On the other hand, Model C does allow for time since death of spouse and this leads to the highest provision at duration 20 (because of the short-term nature of the dependence).

5.3 Joint Life and Survivor Annuity Contract

5.3.1 Premiums

The premiums calculated for the joint life and survivor annuity contract under Models A, B and C appear in Table 11. Recall that we only consider a single premium payment for this type of policy. If independence of coupled lifetimes is assumed instead of short-term dependence, a typical joint life and survivor annuity contract is over-priced by about 2.3% (Model A in comparison to Model C). On the other hand, if coupled lifetimes are assumed to be dependent in a persistent rather than transient way, this contract is under-priced by about 1.8% (Model B in comparison to Model C).

Table 11 Premiums for the joint life and survivor annuity contract

Model C results in a premium that is intermediate between the premiums under Models A and B, which is consistent with the results for the contingent assurance product in section 5.2.1. However, it is now Model A that yields the highest premium, followed by Model C, then Model B. This conforms with intuition. The higher the degree of dependence, the more likely the event that the deaths of the two are separated by a short spell only, and therefore the shorter the expected period during which annuity instalments are payable. The higher premium for Model C relative to Model B follows from the lower mortality experienced by widows after a year in Model C. It appears again that this outweighs the impact of higher mortality in the first year after becoming widowed.

Reassuringly, these results mirror the results reported by Ji et al. (Reference Ji, Hardy and Li2011) when they price joint life and survivor annuities using (a) their semi-Markov model of joint lifetimes (analogous to our Model C above), (b) a Markov model of long-term dependence only (Model B), and (c) a model of independent lifetimes (Model A). As in our model annuity contract described in section 5.1, they also assume that the partners to whom annuities are sold are of the same age. However, they do consider a range of different ages, whereas we price annuities only for partners aged 65 (a typical retirement age). Ji et al. (Reference Ji, Hardy and Li2011) find that the model that allows for short-term dependence generates higher annuity values for younger couples, and lower values for older couples, as compared with the model with long-term dependence only. This agrees with our intuition above: older couples do not enjoy a long enough period of lower mortality after the initial bereavement period of high mortality, as predicted by the model with short-term dependence. For them, the higher mortality soon after bereavement outweighs the lower mortality experienced thereafter and their annuities are therefore cheaper than for younger couples.

5.3.2 Provisions

Table 12 gives the provisions in state 0 for the joint life and survivor annuity. Note again that Model A gives the highest values, and Model B the lowest values. Assuming independence, when short-term dependence is in fact prevalent, leads to an excess of 2–12% in the provisions whereas assuming long-term dependence leads to a shortfall of 1–2.5% in the provisions (Models A and B respectively relative to Model C).

Table 12 Provisions in state 0 for the joint life and survivor annuity contract

In Table 13, we can view the provisions to be held at duration 20 when (x) has died and (y) is still alive. Death of one partner triggers reduced annuity payments on the life of the survivor. The provisions in Table 13 are therefore lower than the values in Table 12. As before for the contingent assurance, the provision at duration 20 has been tabulated according to various periods elapsed since death of the spouse. The independence assumption (in Model A) and long-term dependence assumption (in Model B) lead to significant over-provisioning (by 53–74%) and under-provisioning (by 10–21%) respectively, compared to short-term dependence (Model C).

Table 13 Provisions at duration 20 in state 1 (Model A or B) and in state 1 or 2 (Model C) for the joint life and survivor annuity contract

Table 14 contains the provisions arising at duration 20 when the female dies before the male. As anticipated, the provisions under Model C are intermediate between those of Models A and B. Model A indicates over-provisioning by 69–110% and Model B indicates under-provisioning by 26–41%, compared to Model C. Since male survivor mortality is higher than female survivor mortality, provisions in Table 14 are lower than those in Table 13.

Table 14 Provisions at duration 20 in state 2 (Model A or B) and in state 3 or 4 (Model C) for the joint life and survivor annuity contract

As for the contingent assurance contract, the provision after the first death varies with time since the first death in Model C (in both Tables 13 and 14). Ignoring dependence (in Model A) or assuming the wrong type of dependence (long-term, in Model B) means not only that an insurer commits the wrong amount of capital to reserves, but in addition precludes it from managing the provision dynamically.