2. Mortality Data and Copula Models

We have acquired the mortality data of several countries from the Human Mortality Database (HMD, 2018) for ages 0–89 and the period of 1900–2014. Based on geographical closeness, we sort them into two distinct groups: Nordic countries (Denmark, Finland, Iceland, Norway, Sweden) and Western Europe (England and Wales, France, Italy, Netherlands, Switzerland). During this data period, there were drastic and massive events which led to mortality jumps, including the influenza pandemic in 1918 and World War II in the 1940s (see Figure 1). These data allow us to take into account the possibility of large-scale casualties when pricing mortality-linked securities. Those aged 90 and above are excluded from the analysis because of data volatility and reliability issues.

Figure 1. Lee–Carter mortality indices (${\kappa_{t,i}}$) of Nordic countries and Western Europe.

Sklar (Reference Sklar1959) proved that given continuous marginal distribution functions, there is a unique copula linking the marginal distribution functions for any joint distribution function. This property allows separate modelling treatments between the marginal distributions and the dependency structure. This flexibility leads to numerous choices for handling multivariate random variables and frees us from the traditional restriction of using multivariate elliptical distributions. Accordingly, for each country $i$, we first assume that the (log) central death rate at age $x$ in year $t$ follows the Lee & Carter (Reference Lee and Carter1992) model structure:

(1) \begin{equation}\ln {m_{x,t,i}} = {\alpha_{x,i}} + {\beta_{x,i}}\,{\kappa_{t,i}},\end{equation}

in which ${\alpha_{x,i}}$ is the age schedule and ${\kappa_{t,i}}$ is the mortality index with ${\beta_{x,i}}$ as the sensitivity measure. In estimating the parameters, we use the Poisson assumption in Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002), that is, ${D_{x,t,i}}~{\textrm{Poisson}}({E_{x,t,i}}\,{m_{x,t,i}})$ with the number of deaths ${D_{x,t,i}}$ and exposure ${E_{x,t,i}}$. Many applications of this model focused on more recent data and used the random walk with drift to model ${\kappa_{t,i}}$. Since the data considered here cover a long period of time, we apply the autoregressive integrated moving-average process, ARIMA($p$,1,$q$), instead (Tsay, Reference Tsay2002):

(2) \begin{equation}\Delta {\kappa_{t,i}} = {\varsigma_{0,i}} + \sum\limits_{j=1}^p {{\varsigma_{j,i}}{\rm{ }}\Delta {\kappa_{t - j,i}}} + \sum\limits_{j=1}^q {{\gamma_{j,i}}{\rm{ }}{\varepsilon_{t - j,i}}} + {\varepsilon_{t,i}},\end{equation}

where ${\varsigma_{j,i}}$’s are the autoregressive parameters, ${\gamma_{j,i}}$’s are the moving-average parameters, and ${\varepsilon_{t,i}}$’s are the innovations, which are independent and identically distributed over time.

In order to allow for the possibility of extreme events, besides the normal distribution, we consider the Student $t$, skew normal, and skew $t$ distributions for the innovations. These three distributions possess skewness and/or heavy tails, which make them possible candidates for incorporating the possibility of extreme events being reflected in the innovations. The skew normal distribution proposed by Azzalini (Reference Azzalini1985) takes the form:

(3) \begin{equation}f(z)=2\phi (z)\Phi (\lambda z),\end{equation}

where $\phi$ and $\Phi$ are the standard normal density and distribution functions, and $\lambda$ is the shape parameter. A linear transformation of $X = \mu + \sigma Z$ can be performed with the location parameter $\mu$ and scale parameter $\sigma$. Azzalini & Capitanio (Reference Azzalini and Capitanio2003) then proposed the skew $t$ distribution by setting $Y = \xi + \omega X/\sqrt V$ with the location parameter $\xi $, scale parameter $\omega$, and another shape parameter $\nu$, in which $V~{\rm{Gamma}}(\nu /2,\nu /2)$ is independent of $Z$. The resulting density function is:

(4) \begin{equation}f(y) = \frac{2}{\omega }{t_\nu }\left(\frac{{y - \xi }}{\omega }\right){T_{\nu +1}}\left(\lambda \frac{{y - \xi }}{\omega }\sqrt {\frac{{\nu +1}}{{\nu + {{(y - \xi )}^2}/{\omega ^2}}}} \right),\end{equation}

with ${t_\nu }$ and ${T_\nu }$ being the Student $t$ density and distribution functions with a degree of freedom $\nu $. These two skewed distributions have some desirable tractability due to their affinity to the regular normal and $t$ distributions. For instance, when $\lambda =0$, the former distributions reduce to the latter. When $\nu $ approaches infinity, the skew $t$ distribution converges to the skew normal distribution.

Figure 1 shows the Lee–Carter mortality indices (${\kappa_{t,i}}$) of the 10 countries under consideration. There are clearly a few spikes in the declining mortality trends. Table 1 presents the chosen ARIMA orders for each country based on the Akaike information criterion (AIC)Footnote ¹, significance of the parameters, and residual patterns. All the selected $p$’s and $q$’s are no greater than two, and the computed autoregressive and moving-average parameters are generally much smaller than one in magnitude. Table 1 also lists the probability distributions selected via the chi-square test. The skew $t$ distribution is the optimal choice for Sweden, France, Italy, the $t$ distribution is optimal for Denmark, Finland, Norway, England and Wales, Netherlands, and Switzerland, while the normal distribution is suitable for Iceland. All the fitted $t$ and skew $t$ distributions (except one) have a degree of freedom of about three or less, reflecting a heavy tail for those very large mortality changes and so a proper allowance for the past extreme events.

Table 1. Selected ARIMA processes and marginal distributions

The next step is to model the relationships between the countries (innovations) with Archimedean copulas. Suppose ${X_1}$ and ${X_2}$ are two associated random variables (ARIMA innovations of two countries) having distributions functions ${F_1}$ and ${F_2}$. For the most basic bivariate case, an Archimedean copula is defined as (Nelsen, Reference Nelsen1999):

(5) \begin{equation}C({u_1},{u_2}) = {\varphi ^{ - 1}}(\varphi ({u_1}) + \varphi ({u_2})),\end{equation}

where ${u_1} = {F_1}({x_1})$, ${u_2} = {F_2}({x_2})$, and the generator $\varphi $ is specific for each type of copula. For mathematical tractability, we only include strict Archimedean copulas, in which $\varphi $ is continuous, convex, and strictly decreasing with $\varphi (0) = \infty $ and $\varphi (1)=0$. Table 2 gives a list of 11 choices, including the Clayton (1), Ali-Mikhail-Haq (3), Gumbel–Hougaard (4), Frank (5), and Joe (6) copulas. For each copula, the main parameter is $\theta $. More than half of the copulas in the list have not been much tested in the actuarial literature. We hence develop our own computer algorithm to perform the fitting and simulations in the next two sections.

In this article, we consider multi-dimensional, fully nested, and partially nested Archimedean copulas in turn. As there are five countries in each group, we use a five-dimensional Archimedean copula (Joe, Reference Joe1997):

(6) \begin{equation}C({u_1},{u_2},{u_3},{u_4},{u_5}) = {\varphi ^{ - 1}}(\varphi ({u_1}) + \varphi ({u_2}) + \varphi ({u_3}) + \varphi ({u_4}) + \varphi ({u_5})).\end{equation}

It gives a valid multivariate distribution if ${\varphi ^{ - 1}}$ is completely monotonicFootnote ². This five-dimensional Archimedean copula has been widely studied in empirical finance (e.g. Savu & Trede, Reference Savu and Trede2008; Gaiduchevici, Reference Gaiduchevici2014). Wang et al. (Reference Wang, Yang and Huang2015) applied a four-dimensional version (with Gumbel-Hougaard, Clayton copulas) to the mortality indices of four countries. It is relatively easy to use and is parsimonious in terms of having only one parameter. The main problem of this approach is that the copula structure is permutation symmetric and the exchangeability feature imposes a serious constraint on setting different dependence levels. It would still work reasonably well for countries with a similar extent of association amongst them (e.g. between Norway, Sweden, Finland), but it cannot sufficiently cater for the situation when there are more countries with varying degrees of dependence.

Table 2. A list of strict Archimedean copulas

Note: The copula identification numbers are as given in Nelsen (Reference Nelsen1999). The copulas marked with * give the product copula when $\theta =0$. Those marked with † give the product copula when $\theta =1$. The one marked with ‡ gives the product copula when $\theta = - 1$. The notation U (L) means that the copula has upper (lower) tail dependence.

The second option is a fully nested Archimedean copula (Joe, Reference Joe1997):

(7) \begin{align}C({u_1},{u_2},{u_3},{u_4},{u_5}) &= \varphi_1^{ - 1}({\varphi_1} \circ \varphi_2^{ - 1}({\varphi_2} \circ \varphi_3^{ - 1}({\varphi_3} \circ \varphi_4^{ - 1}({\varphi_4}({u_1}) + {\varphi_4}({u_2}))\nonumber\\[4pt] & \quad + {\varphi_3}({u_3})) + {\varphi_2}({u_4})) + {\varphi_1}({u_5})).\end{align}

The generators ${\varphi_1}$, ${\varphi_2}$, ${\varphi_3}$, and ${\varphi_4}$ have different parameters ${\theta_1}$, ${\theta_2}$, ${\theta_3}$, and ${\theta_4}$. It produces a proper multivariate distribution if ${\varphi_1}$, ${\varphi_2}$, ${\varphi_3}$, and ${\varphi_4}$ are completely monotonic and ${\varphi_1} \circ \varphi_2^{ - 1}$, ${\varphi_2} \circ \varphi_3^{ - 1}$, and ${\varphi_3} \circ \varphi_4^{ - 1}$ belong to ${{\mathcal{L}}}_\infty ^*$ Footnote ³ (Embrechts et al., Reference Embrechts, Lindskog and McNeil2001), which requires the condition ${\theta_1} < {\theta_2} < {\theta_3} < {\theta_4}$. This copula structure is partially exchangeable and is more flexible than the five-dimensional copula above, while being more parsimonious than traditional multivariate elliptical distributions under which the number of correlation parameters increases dramatically with the dimension. We take a similar binary copula approach to Okhrin et al. (Reference Okhrin, Okhrin and Schmid2013) in determining the sequence of the countries and join the two variables with the strongest dependence (in terms of the sample tauFootnote ⁴ or estimated parameter) at each level, starting from ${\varphi_4}$.

A further alternative is a partially nested Archimedean copula (Joe, Reference Joe1997). Its specification depends on the selected hierarchical structure. For example, if a three-dimensional copula is chosen for ${u_1}$, ${u_2}$, and ${u_3}$ and a bivariate copula for ${u_4}$ and ${u_5}$, and if the two groups are then linked by another bivariate copula, the overall structure can be expressed as:

\begin{equation*}C({u_1},{u_2},{u_3},{u_4},{u_5}) = \varphi_1^{ - 1}({\varphi_1} \circ \varphi_2^{ - 1}({\varphi_2}({u_1}) + {\varphi_2}({u_2}) + {\varphi_2}({u_3})) + {\varphi_1} \circ \varphi_3^{ - 1}({\varphi_3}({u_4}) + {\varphi_3}({u_5}))).\end{equation*}

As in the previous approach, the same conditions of completely monotonic ${\varphi_1}$, ${\varphi_2}$, and ${\varphi_3}$, with ${\varphi_1} \circ \varphi_2^{ - 1} \in {{\mathcal{L}}}_\infty ^*$ and ${\varphi_1} \circ \varphi_3^{ - 1} \in {{\mathcal{L}}}_\infty ^*$ (such that ${\theta_1} < {\theta_2}$ and ${\theta_1} < {\theta_3}$), are required. This approach gives even more flexibility in choosing the structure and represents an extension of the previous recursive way. Since the number of possible combinations is huge, we follow the method used by Okhrin & Ristig (Reference Okhrin and Ristig2014). (The R package HAC therein covers only a few commonly used copulas and distributions.) In the first step, a bivariate copula is fitted to every possible couple within the five variables, and the couple with the strongest dependence is selected and then turned into a pseudo-variable using the fitted copula function. In the next step, the procedure is repeated with the remaining variables and the newly formed pseudo-variable as a new set of four variables. The process is iterated until the last couple is handled.

Okhrin & Ristig (Reference Okhrin and Ristig2014) noted that if the unknown, underlying copula is not binary, this process may lead to a slightly mis-specified structure, but the corresponding difference in the distribution functions is generally immaterial. We also modify the method by a final step of finding the global maximum likelihood estimates of all the copula parameters after confirming the hierarchical structure, since the interdependence between different branches in a copula tree is ignored in the iteration process and the resulting stepwise conditional estimates could be suboptimal.

The dimension of the copulas above can be changed and the same arguments apply. As in Okhrin & Ristig (Reference Okhrin and Ristig2014), within each structure, we apply one generator only. The maximum likelihood estimates are calculated from deriving the copula density function:

(8) \begin{equation}c({u_1},{u_2},{u_3},{u_4},{u_5}) = \frac{{{\partial ^5}}}{{\partial {u_1}\partial {u_2}\partial {u_3}\partial {u_4}\partial {u_5}}}C({u_1},{u_2},{u_3},{u_4},{u_5}).\end{equation}

3. Fitting Performance

Table 3 shows the sample taus between the countries under study. Amongst Nordic countries, Norway and Sweden have the strongest overall dependence, while Iceland is comparatively less associated. Within Western Europe, the highest overall dependence occurs between France and Italy and also between England and Wales and Netherlands. All these observations are broadly in line with the geographic locations of the countries. (Zhu et al., Reference Zhu, Tan and Wang2017 showed similar findings for these countries using their hierarchical clustering analysis.) Figure 2 provides a graphical display of the selected layers and orderings of the countries under each of the three approaches, which do not change with the copula functions in most of the cases. (The copula functions shown are based on the results later in Tables 4–6.) It means that the selection of structures is fairly robust to the copula choice. Interestingly, for Nordic countries, the selected partially nested structure turns out to be just the fully nested structure. For Western Europe, the selected partially nested structures vary slightly under a few copula functions.

Table 3. Sample taus amongst Nordic countries and western Europe

Figure 2. Samples of selected hierarchical structures for Nordic countries and Western Europe.

As noted in Savu & Trede (Reference Savu and Trede2008), there is not much consensus on the recommended method to check the goodness-of-fit of a bivariate copula, not to say a multivariate copula. Here we experiment with three different statistical tests. First, we extend the test used in Guégan & Ladoucette (Reference Guégan and Ladoucette2004) to the multivariate case as:

(9) \begin{equation}\sum\limits_{j=1}^n {{{(\tilde {\rm{P}}{\rm{r}}({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}},{X_3} \le {x_{3,j}},{X_4} \le {x_{4,j}},{X_5} \le {x_{5,j}}) - \hat C({u_{1,j}},{u_{2,j}},{u_{3,j}},{u_{4,j}},{u_{5,j}}))}^2}},\end{equation}

which measures the total squared differences between the empirical joint distribution function and the fitted joint distribution function. The copula that leads to the smallest value of this test statistic would be preferred over the other choices.

In the spirit of the Kolmogorov–Smirnov test, we also devise the following test:

(10) \begin{equation}\mathop {\max }\limits_j |\tilde {\rm{P}}{\rm{r}}({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}},{X_3} \le {x_{3,j}},{X_4} \le {x_{4,j}},{X_5} \le {x_{5,j}}) - \hat C({u_{1,j}},{u_{2,j}},{u_{3,j}},{u_{4,j}},{u_{5,j}})|,\end{equation}

which records the maximum difference between the empirical and fitted functions. Again, the copula with a smaller test statistic value is ranked higher. Finally, the AIC is also included in the comparison of copula fitting.

Table 4 provides the testing results for eleven five-dimensional Archimedean copulas. Based on the three statistical tests, for Nordic countries, copulas 4 (Gumbel-Hougaard), 5 (Frank), and 17 are the top three choices overall. (Note that the three lowest values observed under each test are bolded.) By contrast, for Western Europe, the best ones are copulas 14 and 12. These differences reflect that the nature of dependence is dissimilar between the two groups of countries. On the other hand, using fully nested Archimedean copulas, allowing the use of different parameters between countries, tends to produce better fitting performance, as shown in Table 5. Under this approach, Copula 4 is the most optimal choice (followed by copulas 5 and 6 (Joe)) for Nordic countries, while copulas 14 and 12 are the preferred ones for Western Europe. Finally, Table 6 illustrates that using partially nested structures for Western Europe, copulas 14 and 12 are still the best choices, and they tend to produce lower test values in the total squared differences and maximum difference than those in Table 5. These results suggest that the extra flexibility of partially nested structures can potentially further enhance the fitting performance. (For Nordic countries, the selected partially nested structure turns out to be exactly the fully nested structure, so their testing results are the same in Tables 5 and 6.)

Table 4. Statistical tests results of five-dimensional archimedean copulas (with ranks in parentheses)

Table 5. Statistical tests results of fully nested archimedean copulas (with ranks in parentheses)

Table 6. Statistical tests results of partially nested archimedean copulas (with ranks in parentheses)

Nevertheless, the optimal copula choices are quite similar between the three hierarchical approaches. It appears that the copula function choices are affected mainly by the data rather than the hierarchical structures adopted. Note also that the worst candidate so far is copula 20, which shows a poor performance in almost all the cases. It seems that this copula is not suitable for modelling the dependence in mortality rates between countries. Interestingly, by contrast, Luciano et al. (Reference Luciano, Spreeuw and Vigna2008, Reference Luciano, Spreeuw and Vigna2016) demonstrated a good fitting performance of copula 20 in modelling the two survival times of a couple, using a large Canadian dataset. In summary, it can clearly be seen that different copulas with their own unique features would be suitable for mortality data of different neighbouring countries. Although the Gumbel–Hougaard and Clayton copulas are popular choices in the actuarial literature, they are not necessarily always the best options for modelling mortality dependence.

Overall, it appears that the fully nested structure with copula 4 is the most suitable combination for Nordic countries, while a specific partially nested structure (see Figure 2) with copula 14 is the optimal one for Western Europe. One common feature between copulas 4 and 14 in the bivariate case is their upper tail dependence, calculated as $2 - {2^{1/\theta }}$. It is a measure of the association in the upper-right-quadrant tail. This implied upper tail dependence of the fitted copulas reflects possible mortality co-movements between the countries in adverse incidents. On the other hand, copula 14 also has lower tail dependence but copula 4 does not. It may imply that the recovery from extreme mortality is more prominent and in line amongst the countries in Western Europe than between Nordic countries (see Figure 1).

4. Pricing of Mortality Bonds

In this section, we consider a typical mortality bond structure based on the securities issued under Swiss Re’s Vita program (e.g. Lin & Cox, Reference Lin and Cox2008; Chen & Cox, Reference Chen and Cox2009). Let the combined mortality index (CMI) in year $t$ be ${q_t} = (q_t^{(1)} + q_t^{(2)} + q_t^{(3)} + q_t^{(4)} + q_t^{(5)})/5$ for the five countries in a group with the same weights, $a$ be the attachment point, and $b$ be the detachment point, in which $q_t^{(j)}$ is the population mortality rate of country $j$. Suppose the maturity of the bond is 5 years. The principal loss in percentage in year $t$ can be defined as:

\begin{equation*}{L_t} = \min\!\left(1,\max \!\left(0,\frac{{{q_t} - a}}{{b - a}}\right)\right) \times 100\%, \end{equation*}

where the losses (if any) accrue over time, and the entire principal is exhausted to zero when the aggregated losses exceed 100%. The coupons are payable annually on the remaining principal at a rate of the risk-free rate plus a spread. Three settings of $(a,b)$ are tested: $(1.05,{\rm{ }}1.10) \times {\rm{base level}}$, $(1.10,{\rm{ }}1.15) \times {\rm{base level}}$, and $(1.15,{\rm{ }}1.20) \times {\rm{base level}}$, in which the base level is assumed to be the CMI observed in 2014. The principal loss in a year may be seen as a combination of a long call and another short call option with a higher strike price.

We use the selected models from the previous sections to generate 10,000 paths of future death rates, each with an equal real-world probability. We adopt the following procedure in Mathematica to simulate multivariate random variables from the fitted hierarchical copula structures (e.g. Embrechts et al., Reference Embrechts, Lindskog and McNeil2001):

1. (1) Define ${C_i}({u_1},{u_2},...,{u_i}) = C({u_1},{u_2},...,{u_i},1,...,1)$ for $i=1,{\rm{ }}2,{\rm{ }}3,{\rm{ }}4,{\rm{ }}5$.

2. (2) Derive ${C_i}({u_i}|{u_1},{u_2},...,{u_{i - 1}}) = \frac{{\dfrac{{{\partial ^{i - 1}}{C_i}({u_1},{u_2},...,{u_i})}}{{\partial {u_1}\partial {u_2}...\partial {u_{i - 1}}}}}}{{\dfrac{{{\partial ^{i - 1}}{C_{i - 1}}({u_1},{u_2},...,{u_{i - 1}})}}{{\partial {u_1}\partial {u_2}...\partial {u_{i - 1}}}}}}$ Footnote ⁵ for $i=2,{\rm{ }}3,{\rm{ }}4,{\rm{ }}5$.

3. (3) Simulate ${V_1}~ {\rm{U}}(0,1)$.

4. (4) Simulate ${Q_i}~ {\rm{U}}(0,1)$ and compute ${V_i} = C_i^{ - 1}({Q_i}|{u_1} = {V_1},{u_2} = {V_2},...,{u_{i - 1}} = {V_{i - 1}})$ for $i=2,{\rm{ }}3,{\rm{ }}4,{\rm{ }}5$.

5. (5) Calculate ${X_i} = F_i^{ - 1}({V_i})$ for $i=1,{\rm{ }}2,{\rm{ }}3,{\rm{ }}4,{\rm{ }}5$.

6. (6) Repeat steps (3)–(5) for 10,000 times for each future year.

Then we apply the Wang transform (Wang, Reference Wang2000) to deduce the risk-neutral probabilities and calculate the spread. The Wang transform can be specified as

\begin{equation*}{F^*}(x) = \Phi ({\Phi ^{ - 1}}(F(x)) + \gamma ),\end{equation*}

where $F$ is the real-world distribution function, ${F^*}$ is the risk-neutral distribution function, $\Phi $ is the standard normal distribution function, and $\gamma $ is the market price of risk. We set $\gamma = - 0.45$, which was the average market price of risk found by Wang (Reference Wang2004). The risk-neutral distribution has a heavier right tail than the real-world distribution, and it represents the market perception of large-scale mortality risks. As the current interest rates are very low globally (e.g. 5-year UK Gilt yield was 0.81% p.a. as at 21 March 2019), we assume that the risk-free rate is 1% p.a. flat.

Table 7 displays the estimated spreads for Nordic countries and Western Europe and for the three sets of attachment and detachment points. There are two general observations. First, as the attachment and detachment points increase, the spread decreases. This pattern is driven by the fact that the probability of reaching a higher mortality level is lower than otherwise. Second, the spreads estimated for Western Europe are much higher than those for Nordic countries. It is a reflection of the severity of past extreme events in Western Europe compared to Nordic countries. The computed spreads for Western Europe are broadly in line with those being offered in the market for similar countries (e.g. Chen et al., Reference Chen, MacMinn and Sun2015).

Table 7. Estimated spreads (in basis points) for a mortality bond

For further comparison, Table 7 also includes the results from using only the three most associated countries in each group and from using copula 5 instead. With fewer countries constituting the CMI, there is less diversification, and the resulting spreads are much higher than previously. The impact of diversification appears to be very large in our demonstration. Moreover, copula 5 does not possess upper tail dependence. The calculated spreads are then clearly lower than those from copulas 4 or 14. These results show that upper tail dependence is an important feature to consider when choosing a copula function to model mortality dependence.

Table 8 gives the estimated (real-world) probability of any principal loss (partial or total) in each case. These probabilities are in agreement with the spreads in Table 7. The probability of principal loss is higher for Western Europe, if there are fewer countries in the CMI, and under copulas 4 and 14. Approximately, each additional 1% chance of principal loss would induce an extra risk premium of about 30–40 basis points in our calculations.

Table 8. Estimated probabilities of principal loss for a mortality bond

Similar pricing exercises can be performed for other mortality-linked securities such as q-forward (Coughlan et al., Reference Coughlan, Epstein, Sinha and Honig2007) and mortality option (Cairns et al., Reference Cairns, Blake and Dowd2008) to incorporate mortality or catastrophe risks. Despite that some previous work used copulas to price longevity-linked securities (e.g. Chen et al., Reference Chen, MacMinn and Sun2017; Zhu et al., Reference Zhu, Tan and Wang2017), we deem that copulas by nature are more suitable for modelling temporary, one-off mortality jumps rather than modelling long-term longevity improvements, particularly longevity risk. Empirically, we find that the fitted time series processes for the mortality indices are generally weakly stationary and the resulting projections converge to their long-term trends rather quickly. The potential volatility or deviation allowed for by the use of copulas would tend to be short-term. It is probably better to use trend-shifting models (e.g. Li et al., Reference Li, Chan and Cheung2011) or regime-switching models (e.g. Milidonis et al., Reference Milidonis, Lin and Cox2011) instead to cater for radical permanent changes in the longevity trend.

5. Concluding Remarks

We have demonstrated the use of multi-dimensional, fully nested, and partially nested Archimedean copulas in co-modelling the Lee–Carter mortality indices of multiple countries and pricing mortality bonds. In particular, we have tested the application of 11 types of generators and also the skew normal and skew t distributions. Our fitting results suggest that copulas with upper tail dependence, hierarchical structures with more flexibility, and distributions with heavy tails provide a decent way to model mortality dependence and allow for extreme mortality risks. Moreover, the calculated prices of mortality bonds strongly reflect the mortality history of the countries involved and also the level of diversification between the countries. We argue that hierarchical Archimedean copulas would be more suitable for modelling short-term mortality jumps and pricing mortality-linked securities than for modelling long-term longevity improvements and pricing longevity-linked securities.

In this research, we have used the popular Lee–Carter model to generate the mortality indices. For future research, it would be useful to perform a similar investigation on the well-documented CBD model (Cairns et al., Reference Cairns, Blake and Dowd2006) for comparison. Recently, Zhu et al. (Reference Zhu, Tan and Wang2017) suggested the use of the Lévy subordinated hierarchical Archimedean copulas (with Gumbel-Hougaard, Clayton, Joe copulas) with the CBD model and Gaussian errors for modelling mortality dependence and pricing survivor swaps for a cohort aged 65. It would be interesting to further compare with their approach by extending it with the Lee–Carter model on a wider age range, other less commonly used copula functions, heavy-tailed distributions, and pricing mortality-linked securities. It would also be interesting to compare with the vine copula approach in Zhou (Reference Zhou2019), as the hierarchical Archimedean copulas and the vine copulas are the two main state-of-the-art alternatives in modelling multivariate random variables in the finance and other disciplines.