2. Archimedean Copulas

Archimedean copulas are a group of copulas that share certain fundamental characteristics. There are a variety of options for modelling the dependency structure. For the bivariate case, considering two associated random variables ${X_1}$ and ${X_2}$ with distribution functions ${F_1}$ and ${F_2}$ , they have the basic form (e.g. Nelsen, Reference Nelsen1999)

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

in which ${u_1} = {F_1}({x_1})$ , ${u_2} = {F_2}({x_2})$ , and $\varphi $ is a function called generator with inverse ${\varphi ^{ - 1}}$ . Each type of Archimedean copula has a unique form of $\varphi $ . When $\varphi $ is continuous, convex, and strictly decreasing with $\varphi (0) = \infty $ and $\varphi (1) = 0$ , $C$ is a strict Archimedean copula. Table 1 presents a list of 13 strict one-parameter Archimedean copulas that are considered in this paper. In particular, copulas 1, 3, 4, 5, 6, and 9 belong to the Clayton, Ali–Mikhail–Haq, Gumbel–Hougaard, Frank, Joe, and Gumbel–Barnett families, respectively. Table 2 provides the formulae of Kendall’s tauFootnote ¹ , upper tail dependence ( ${\lambda _U}$ )Footnote ² , and lower tail dependence ( ${\lambda _L}$ )Footnote ³ that we have derived for the copulas. The parameter $\theta $ can be estimated from the sample tau or maximum likelihood, the latter of which is used here (see Section 6 for more estimation methods). The copula density function can be derived as

\begin{align*}c({u_1},{u_2}) = \frac{{{\partial ^2}}}{{\partial {u_1}\partial {u_2}}}C({u_1},{u_2})\end{align*}

based on which the log-likelihood function can be evaluated.

Table 1. Bivariate Archimedean copulas.

Note: Copulas 1, 3, 4, 5, 6, and 9 are the Clayton, Ali–Mikhail–Haq, Gumbel–Hougaard, Frank, Joe, and Gumbel–Barnett copulas, respectively. Copulas 1, 3, 5, 9, 10, and 20 are equivalent to the product copula when $\theta = 0$ . Copulas 4, 6, and 13 produce the product copula when $\theta = 1$ . Copula 17 becomes the product copula when $\theta = - 1$ . The copula numbers are those used by Nelsen (Reference Nelsen1999).

Table 2. Kendall’s tau, upper tail dependence, and lower tail dependence.

Note: The copula numbers are those used by Nelsen (Reference Nelsen1999).

It is interesting to notice the different features that these copulas possess. Seven of them have a positive (or zero for some) tau, two of them have a negative (or zero) tau, and the rest can have a positive or negative (or zero) tau. Copulas 1, 19, and 20 have lower tail dependence only, copulas 4 and 6 have upper tail dependence only, copulas 12 and 14 have both upper and lower tail dependence, while the others have neither. This variety of combinations offers a lot of flexibility for modelling dependence in insurance and financial risks.

To identify the most optimal copula for modelling the dependency structure of a particular data set of size $n$ , we use the test statistic (Guégan & Ladoucette, Reference Guégan and Ladoucette2004) below:

$$\sum\limits_{j = 1}^n {{{(\tilde{\text{Pr}}({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}}) - \hat C({u_{1,j}},{u_{2,j}}))}^2}} $$

which is the total squared differencesFootnote ⁴ between the empirical bivariate distribution function (denoted as ~) and the fitted bivariate distribution function $\hat C$ . The copula which gives the smallest test statistic value would be chosen over the other candidates. Note that this overall test covers the entire plane of the two random variables. To supplement the analysis and put more focus on tail risks, we also compute the likelihood ratio test statistic (Kupiec, Reference Kupiec1995):

$$ 2(\ln ({\tilde p^m}{(1 - \tilde p)^{n - m}}) - \ln ({\hat p^m}{(1 - \hat p)^{n - m}}))$$

in which $ p = \Pr ({X_1} \gt F_1^{ - 1}(1 - \gamma ),{X_2} \gt F_2^{ - 1}(1 - \gamma ))$ has empirical estimate $ \tilde p = m/n $ , where $ m $ is the observed number of data points falling in the $ 100\gamma $ % upper-right-quadrant region, and $\hat p$ is calculated from the fitted copula. A smaller test statistic value is preferred, and the (one-sided) critical value can be taken as the 95th percentile of the $\chi _1^2$ distribution (3.84). Surprisingly, testing of the tail performance, which is a major strength of copula modelling, has typically been omitted in previous actuarial applications. Strictly speaking, tail dependence is a limiting property which cannot be tested for a finite data set. But the likelihood ratio test provides a practical and convenient way to assess how close a copula models the association in the tail. Note also that when the null hypothesis is not rejected, the fitted copula can be regarded as giving a good description of the shape and extent of association, but it does not mean that the copula is necessarily a true representation of the underlying relationship. For other overall goodness-of-fit tests for copulas, see Genest et al. (Reference Genest, Rémillard and Beaudoin2009).

3. Modelling Bivariate Claims in General Insurance

The first data set contains 1,500 general liability claims, which were previously studied by Frees & Valdez (Reference Frees and Valdez1998), Klugman & Parsa (Reference Klugman and Parsa1999), and Denuit et al. (Reference Denuit, Purcaru and Van Keilegom2006). Each claim has an indemnity payment and an allocated loss adjustment expense. The policy limits are recorded for the majority of the claims, in which 34 claim payments are censored from above. For the 148 claims with no recording of the policy limit, it is simply assumed that the policy limit does not exist. Figure 1 plots the log payments against the log expenses, and also the corresponding copula data from their empirical distribution functions. There is a clear positive relationship (left graph) between the two quantities ( $\hat \tau = 0.31$ ), especially in the upper-right-quadrant tail, which can be modelled by a suitable copula structure. Note that for those censored claims, the contribution to the likelihood function has to be adjusted as $1 - {\textstyle{\partial \over {\partial {u_2}}}}C({u_1},{u_2})$ (e.g. Frees & Valdez, Reference Frees and Valdez1998), in which ${u_1}$ refers to the censored payments. It also appears that (right graph) the two variables are exchangeable, that is, the dependency structure is symmetric between the two arguments in the copula function. A Cramér–von Mises type test (Genest et al., Reference Genest, Nešlehová and Quessy2012) shows that the null hypothesis of symmetry cannot be rejected at 5% significance level. Hence, Archimedean copulas, which are symmetric by nature, are suitable candidates for modelling the losses.

Figure 1 Payments versus expenses in logarithmic scale (left) and copula data (right).

Under Sklar’s (Reference Sklar1959) theorem, if the marginal distribution functions are continuous, there exists a unique copula that links the marginal distribution functions to form the corresponding joint distribution function. It implicitly allows a separate consideration between selecting the marginal distributions and choosing the dependency structure. Accordingly, we first examine the marginal distributions. Figure 2 shows two histograms of the claim payments and associated expenses. Both distributions appear to be heavily skewed to the right. We then test the suitability of the gamma, Weibull, lognormal, Pareto, and Burr distributions on each of the two components. The lognormal and Burr distributions are selected, respectively, for the payments and expenses, based on the results from the chi-square test.

Table 3. Fitting results of bivariate Archimedean copulas on censored general liability claims.

Figure 2 Histograms of claim payments and associated expenses.

Using the two tests in Section 2 together can give a better insight into how well a copula can describe the overall association as well as the level of dependence in the right tails between the two related quantities. As the data are censored, however, the tests cannot be conducted precisely. For the 34 censored claims, 22 of their claim payments (i.e. policy limits) are close to or more than the estimated 90th percentile, while the others have their associated expenses below the corresponding 90th percentile. As an approximation, we assume that the underlying losses of all these 22 cases are above the 90th percentile and apply the likelihood ratio test at $\gamma = 0.1$ (noted as + in Table 3). For comparison, we also exclude the 34 censored claims and conduct the two tests (noted as ^*). Furthermore, we invert all the copulas by 180° (i.e. $c(1 - {u_1},1 - {u_2})$ , the so-called survivor copula) and examine their performance as well. Any lower (upper) tail dependence property will then become upper (lower) tail dependence after rotation.

The fitting results of the general liability claims are shown in Table 3. Under each test, the rankings are stated in parentheses, and the best three options are made bold. It is clear that copulas 4 and 6, which both have upper tail dependence, outperform all the others in terms of both the whole data range (2nd column) and the upper-right-quadrant tail (3rd and 4th columns). In fact, the null hypothesis is rejected at 5% significance level for all the other copulas under the likelihood ratio test. The results of the 180° rotated versions (5th to 7th columns) are relatively more dispersed. It appears that copulas 1, 3, and 20 produce a reasonable fit both overall and in the right tails, and their performance is largely comparable to those of the two non-rotated copulas selected above. It is interesting to note that while copulas 1 and 20 have upper tail dependence after rotation, copula 3 does not have tail dependence, which is reflected in its lower performance in the right tails when compared to the other chosen models. Note also that although copulas 12, 14, and 19 fit the uncensored claims reasonably well overall, they are not the best choices for modelling the association in the largest claims, as shown by the likelihood ratio test statistics. Furthermore, considering the possible ranges of tau, we also test the 90° and 270° rotated versions (i.e. $c(1 - {u_1},{u_2})$ and $c({u_1},1 - {u_2})$ _, respectively) of copulas 9 and 10, but we find that their fitting performance (not shown here) turns out to be worse than many of those in Table 3.

In general insurance pricing, valuation, and capital assessment, it is very important to measure the magnitude of tail events adequately. We use each of the fitted copulas selected above (and also the Gaussian copula for comparison) to simulate one million scenarios, and compute the expected shortfall (conditional value-at-risk) as the sample mean of those simulated total claim amounts which exceed a certain percentile. Table 4 shows that our estimates of expected shortfall are very close between the selected copulas at the 95th and lower percentiles. The differences become more obvious starting from the 99th percentile. The rotated copulas 1 and 20 give the largest estimates, followed by copulas 4 and 6, in which all these four copulas possess upper tail dependence. By contrast, the rotated copula 3 and the Gaussian copula, without tail dependence, produce smaller estimates, the significance of which increases sharply with the percentile level. It is clear that the tail dependence property plays a significant role in assessing extreme events and should be a major consideration in copula fitting. Using copulas without any tail dependence runs the risk of a serious underestimation of the capital requirement for an insurer. Regarding the final choice between the four copulas (4, 6, rotated 1, and rotated 20), the decision may not be straightforward under the data constraints (e.g. censoring issues, no recording of the policy limit for some claims). Additional data fields may need to be collected in order to make a better differentiation; or else an actuarial judgement has to be made on how conservative the estimation should be, depending on the purpose of the analysis.

Table 4. Estimated expected shortfall of total claim amounts under selected copulas.

We now turn to our second data set which comprises 2,167 Danish fire insurance claims (collected from R package “fitdistrplus”), which was used earlier by Haug et al. (Reference Haug, Klüppelberg and Peng2011). Each claim has a loss amount to building and a loss amount to contents. The claims are recorded only for those cases with the sum of losses to building and contents and profits being at least 1 million Danish Kroner (DK). To avoid the potential negative dependence between the items induced by this restriction, we model only the 301 claims which have both the losses to building and contents each greater than or equal to 1 million DK. Figure 3 displays these losses to building and contents on a logarithmic scale. Some positive dependence (left graph) between the two amounts can be identified ( $\hat \tau = 0.21$ ). It also shows that (right graph) the two variables are exchangeable, in which the null hypothesis of symmetry is not rejected at 5% significance level. Figure 4 illustrates two histograms of the losses to building and contents, which show heavy right skewness. Since the selected data are truncated from below, we deduce the conditional marginal densities as ${f_1}({x_1})/(1 - {F_1}(1))$ and ${f_2}({x_2})/(1 - {F_2}(1))$ , and the conditional joint distribution function as

\begin{align*}\Pr ({X_1} \le {x_1},{X_2} \le {x_2}|{X_1},{X_2} \ge 1) = \frac{{F({x_1},{x_2}) - F({x_1},1) - F(1,{x_2}) + F(1,1)}}{{1 - {F_1}(1) - {F_2}(1) + F(1,1)}}\end{align*}

in which ${f_1}$ and ${f_2}$ are the (unconditional) marginal densities and $F$ is the (unconditional) bivariate distribution function for the losses in millions of DK. The copula density functions and the test statistics are adjusted accordingly. Based on the chi-square test, the lognormal and Pareto distributions are chosen for the losses to building and contents, respectively.

Figure 3 Losses to building versus losses to contents in logarithmic scale (left) and copula data (right).

Figure 4 Histograms of losses to building and contents.

Table 5 provides the fitting results of the fire insurance claims. To examine the goodness of fit in the right tails, three cases of $\gamma = 0.01$ , $\gamma = 0.1$ , and $\gamma = 0.25$ are considered. It can be seen that the results are less clear-cut than those in Table 3, probably due to the much smaller sized data set after truncation. Under the overall test (2nd column), copulas 9 and 10 lead to the best fit, but their performances in the tails (3rd to 5th columns) are poor and many of the cases are rejected. In fact, their parameter estimates are close to zero, which refer to the product copula with independent marginals. This implication is clearly not supported by the plots in Figure 3. By contrast, copulas 4, 6, 12, and 14, which possess upper tail dependence, seem to provide a good description of both the overall data range and the tails. In particular, copulas 4 and 6 perform better at $\gamma = 0.25$ , while copulas 12 and 14 perform better in the extreme tails. For the 180° rotated versions (6th to 9th columns), copulas 1 and 20, which have upper tail dependence after rotation, seem to produce a good fit. Interestingly, the rotated copulas 4 and 13 also give some comparable results, but they do not have upper tail dependence.

Table 5. Fitting results of bivariate Archimedean copulas on truncated fire insurance claims.

Table 6 presents our estimates of expected shortfall for the total losses. Again, the differences between the selected copulas are more obvious at the 99th and higher percentiles. Copula 4 clearly generates the largest estimates. The next largest ones are produced by copulas 6, 14, and the rotated copulas 1 and 20, the results of which are very similar. All these copulas have upper tail dependence, though, interestingly, copula 12 yields smaller estimates. On the other hand, the rotated copulas 4 and 13, and the Gaussian copula, having no upper tail dependence, produce much smaller estimates. These results highlight once more the importance of tail dependence in the allowance for tail risks. Note that the small size of the truncated data here limits the scope of our investigation. Nevertheless, it is quite reassuring to realise that as long as the focus is on tail events and the chosen copulas have upper tail dependence, the resulting tail estimates tend to be fairly consistent with one another.

Table 6. Estimated expected shortfall of total losses under selected copulas.

4. Modelling Mortality Dependence and Pricing Mortality Bonds

We have obtained the mortality data of England and Wales, Netherlands, Norway, and Sweden from the Human Mortality Database (HMD, 2018) for the period of 1900 to 2014. Figure 5 plots the log mortality indices ( $\ln q_t^{(i)}$ ) of these four regions. There are obviously some extreme mortality co-movements caused by epidemics and wars before 1950. Since a mortality or catastrophe bond typically has its payments dependent on a weighted mortality index between pre-specified populations (e.g. Swiss Re’s catastrophe bonds), an allowance for such mortality co-movements in extreme events should be made adequately. Otherwise, the potential impact would be understated, and the bond’s premium would then be underestimated. Suppose $q_t^{(i)}$ is the mortality rate or index in year $t$ of population $i$ for the entire age range. Figure 6 illustrates the sample autocorrelations of $\Delta \ln q_t^{(i)}$ and ${(\Delta \ln q_t^{(i)})^2}$ . The significant patterns shown in the graphs call for time series modelling with a suitably chosen ARIMA-GARCH process.

Figure 5 Log mortality indices.

Figure 6 Sample autocorrelations of $\Delta \ln q_t^{(i)}$ and ${(\Delta \ln q_t^{(i)})^2}$ .

As in the previous section, we consider the marginal distributions and the dependency structure separately. We first use the ARIMA-GARCH process to remove the autocorrelations and conditional heteroskedasticity from each region’s mortality trend over time, and then model their innovations by various Archimedean copulas. The ARIMA( $p$ , 1, $q$ )-GARCH( $m$ , $s$ ) process for $\ln q_t^{(i)}$ is taken as (e.g. Tsay, Reference Tsay2002)

\begin{align*}\Delta \ln q_t^{(i)} &= {\phi _0} + \sum\limits_{j = 1}^p {{\phi _j}\Delta \ln q_{t - j}^{(i)}} + \sum\limits_{j = 1}^q {{\varphi _j}{a_{t - j}}} + {a_t}\\ \hspace*{-60pt} a_t=\sigma_t\varepsilon _t\qquad\qquad\qquad\qquad\qquad\sigma _t^2 &= {\alpha _0} +\sum\limits_{j = 1}^m {{\alpha _j}a_{t - j}^2} + \sum\limits_{j = 1}^s {{\beta _j}\sigma _{t - j}^2}\end{align*}

where ${\phi _j}$ ’s are the autoregressive parameters, ${\varphi _j}$ ’s are the moving-average parameters, ${a_t}$ is the innovation, ${\sigma _t}$ is the time-varying volatility with its parameters ${\alpha _j}$ and ${\beta _j}$ , and ${\varepsilon _t}$ ’s are independent and identically distributed across time with mean zero and variance equal to one.

Table 7 lists the selected orders for each population based on the Akaike information criterion (AIC), Bayesian information criterion (BIC), and parameter significanceFootnote ⁵ . Figure 7 shows the sample autocorrelations of the standardised residuals from the fitted processes. The previous significant residual patterns are largely removed. Table 7 also gives the selected distributions based on the chi-square test. We consider three probability distributions including the normal, Student $t$ , and skewed Student distributions (e.g. Fernandez & Steel, Reference Fernandez and Steel1998). The normal distribution is the preferred option for England and Wales and Netherlands, and the skewed Student distribution is the preferred one for Norway and Sweden.

Table 7. Selected time series processes and probability distributions.

Figure 7 Sample autocorrelations of standardised residuals.

Table 8 presents the results of fitting the 13 Archimedean copulas to two pairs of the populations under study ( ${\hat \tau _{{\rm{E \& W,NLD}}}} = 0.44$ and ${\hat \tau _{{\rm{NOR,SWE}}}} = 0.44$ ). The relative performances between the copulas are more consistent than those in Table 5. For the likelihood ratio test, the rankings are mostly in line between $\gamma = 0.05$ and $\gamma = 0.1$ . However, there are some small differences between these rankings and those under the total squared differences test, which suggest that the best overall-fitting copula is not necessarily the best tail-fitting copula. Between England and Wales and Netherlands, the best overall-fitting one is copula 5, but its rankings are 5th and 2nd in the likelihood ratio tests. Similarly, between Norway and Sweden, the best overall-fitting one is copula 14, with its rankings being only 9th and 8th in the likelihood ratio tests. The final choice depends on the purpose of the modelling – since we will use the fitted models to price a mortality bond and incorporate extreme risks, we would need to strike a balance between different aspects of the fitting performance. In general, for England and Wales and Netherlands, the optimal choice appears to be copula 14, since its performance is excellent in the tails (2nd and 1st), while providing a reasonable overall fit (4th). For Norway and Sweden, the optimal choice appears to be copula 17 (followed probably by copulas 5 and 13), which shows good performance (2nd) in all the tests. Note that copula 14 has upper tail dependence while copula 17 does not. These choices reflect the severity of historical simultaneous mortality jumps of the two pairs of populations. By contrast, copulas 9, 10, and 20 are clearly not suitable for modelling the mortality data here. Particularly, copulas 9 and 10 are rejected in six of the eight likelihood ratio tests. Overall, it is obvious that different copulas with varying characteristics would be suitable for dealing with different countries’ data. Though being commonly used in the actuarial literature, the Clayton, Gumbel–Hougaard, and Frank copulas are indeed not the only feasible choices for mortality dependence modelling. (Note that the null hypothesis of symmetry is not rejected at 5% significance level for both pairs of countries).

Table 8. Fitting results of bivariate Archimedean copulas for neighbouring countries.

We now adopt the selected models from above and consider a mortality bond structure in line with those issued by the Vita programme of Swiss Re. Suppose the combined mortality index (CMI) of the bond is specified as ${q_t} = (q_t^{(i)} + q_t^{(k)})/2$ for two regions $i$ and $k$ with an equal weight of 0.5. Let $a$ be the attachment point and $b$ be the detachment point of the bond, with a term of $T$ years. The principal loss (in proportion) during year $t$ can be expressed as

\begin{align*}{L_t} = \min \left(1,\max \left(0,\frac{{{q_t} - a}}{{b - a}}\right)\right)\end{align*}

in which any losses will accumulate over the term and the whole principal will be depleted when the total principal losses exceed one. The coupon rate is set as the risk-free rate plus a spread (i.e. a premium for taking mortality risks). We simulate 100,000 scenarios of future mortality indices with an equal real-world probability of 0.00001 and adopt the Wang transform (Wang, Reference Wang2000) to calculate the risk-neutral probabilities and the spread. Based on Wang (Reference Wang2004), the market price of risk is assumed to be −0.45. It is also assumed that the risk-free rate is 1% p.a., consistent with the current low interest rate environment.

We include two sets of attachment and detachment points in our example, in which the base level is taken as the CMI in year 2014. The term of the bond is supposed to be 5 years. Table 9 provides the estimated spreads for each pair of populations. For comparison, we also obtain the results from two other copulas with similar tail dependence properties. As expected, the spread required is lower for the higher attachment and detachment points, since the chance for the mortality levels to hit the higher range is lower. Moreover, the spreads estimated for England and Wales and Netherlands are higher than the spreads computed for Norway and Sweden, the results of which are in line with the upper tail dependence properties of the selected copulas. Lastly, the spreads calculated from the selected copulas 14 and 17 are very close to those from copulas 12 and 5. These similarities can be explained by their common tail dependence properties, despite the overall copula differences.

Table 9. Estimated spreads for a mortality bond based on two neighbouring countries.

5. Simulation Exercise

In this section, we perform a large-scale simulation study to investigate further the potential problems in copula selection. We simulate random samples of varying data sizes ( $n$ = 100, 500, 1,000) from each copula using different tau values ( $\tau $ = 0.25, 0.5, 0.75) and fit all the copulas in the list (excluding copulas 3, 9, and 10, due to their ranges of tau) to these simulated samples. Embrechts et al. (Reference Embrechts, Lindskog and McNeil2001) provided an algorithm for generating bivariate random samples from an Archimedean copula. First, two random numbers, $S$ and $Q$ , are generated from the Uniform(0, 1) distribution. Then two associated Uniform(0, 1) samples, ${U_1}$ and ${U_2}$ , are computed from ${\varphi ^{ - 1}}(S\varphi (T))$ and ${\varphi ^{ - 1}}((1 - S)\varphi (T))$ _, respectively, in which $T = K_C^{ - 1}(Q)$ and ${K_C}(t) = t - \varphi (t)/\varphi '(t)$ . Lastly, samples of two associated variables of interest are produced by using ${U_1}$ and ${U_2}$ in the inverse transform method. Based on 1,000 simulated scenarios, we examine how the data size and the level of dependence impact the results of statistical tests and the accuracy of final copula choice.

As in the previous sections, the fitting performances are ranked under each test in every simulated scenario. Tables 10, 11 and 12 give the number of times (out of 1,000 scenarios) that the “correct” copula is selected as the best choice or one of the best three choices, for each copula in Table 1 taken in turn as the simulation model. Some interesting observations are noted below:

1. 1. As expected, the frequency of picking the underlying copula assumption correctly increases with the data size generally. The increase tends to be more noticeable for the likelihood ratio test at the extreme tails (when there is tail dependence), because the accuracy in capturing the tails is very low when the sample size is small.

2. 2. The overall test results are more stable and accurate than those of the likelihood ratio test when the data size is smaller. But when the size is large, the accuracy of the likelihood ratio test at the 5th and 95th percentiles becomes more comparable to (and sometimes better than) that of the overall test.

3. 3. Even when the sample size is 1,000, the accuracy level of the likelihood ratio test at the 1st and 99th percentiles is still quite low in many cases. Regardless of tail dependence, the performance at the 5th and 95th percentiles is often better than at the 1st and 99th percentiles.

4. 4. While the accuracy levels differ between the various copula assumptions, the accuracy seems to have a tendency to be higher when the tau value is 0.5. Comparatively, for a lower tau of 0.25, the samples are more randomly scattered, and the specific copula properties may become less important, making the selection more ambiguous. For a higher tau of 0.75, the simulated samples are more concentrated around the diagonal line. Those other copulas with similar overall and tail dependence properties to the simulation assumption would then be more competitive against the simulating copula and achieve high rankings (not shown here). For instance, copulas 4 and 6 (sometimes joined by copulas 12 and 14) would compete with one another, and the same happens between copulas 1, 19, and 20, and between copulas 5 and 17.

5. 5. The performance of copula 17 is hampered by the dominance of copula 5 with a similar shape. Although copula 13 has no tail dependence, it does have a lower tail shape, which draws some competition from those copulas with lower tail dependence. Moreover, copulas 12 and 14 have both upper and lower tail dependence, which cause them subject to competition from more copulas.

Table 10. Simulated frequencies based on 1,000 scenarios of selecting the correct copulas as the best choice (left figure) or one of the best three choices (right figure), for n = 100 (top), n = 500 (middle), and n = 1000 (bottom), when τ = 0.25.

Note: When copulas 12, 14, and 19 are used, it is assumed that τ = 0.35 instead.

Table 11. Simulated frequencies based on 1,000 scenarios of selecting the correct copulas as the best choice (left figure) or one of the best three choices (right figure), for n = 100 (top), n = 500 (middle), and n = 1000 (bottom), when τ = 0.5.

Table 12. Simulated frequencies based on 1,000 scenarios of selecting the correct copulas as the best choice (left figure) or one of the best three choices (right figure), for n = 100 (top), n = 500 (middle), and n = 1000 (bottom), when τ = 0.75.

There are some major implications from this simulation study. First, the overall test and the tests on the tails can supplement each other well, especially when the data size is large. The simultaneous use of different tests would give a much better picture of the model fitting. Second, it appears that a sample size of at least a few thousand is needed if one wants to really identify the copula very accurately. When the sample size is 1,000, using the overall test, the average frequency of choosing the correct copula as the best choice is only about 45%, though the average frequency of including the correct copula in the best three is around 85%. The data sizes in the previous sections are not large enough in this sense, and unavoidably, there is some degree of model uncertainty, which requires certain judgement to make the final choice. Despite the difficulty to differentiate between the copula choices at times, it is comforting to see from the previous sections that the estimates of tail measures are fairly consistent between using copulas with similar tail dependence properties.

6. Other Estimation Methods

In this section, we first consider five other estimation methods that have not been tested in previous actuarial applications. These methods include matching theoretical and sample values of Blomqvist’s beta, maximum likelihood using the diagonal of a copula, and three minimum distance estimators based on the total squared differences, total absolute differences, and maximum absolute difference, respectively, between the empirical joint distribution function and the fitted joint distribution function. Moreover, we experiment with a non-parametric estimator for Archimedean copulas, as well as two recently proposed empirical copulas called the empirical beta copula and empirical checkerboard copula, and compare their resulting estimates of tail measures.

Blomqvist’s beta (e.g. Nelsen, Reference Nelsen1999) is defined as $\Pr (({X_1} - F_1^{ - 1}(0.5))({X_2} - F_2^{ - 1}(0.5)) \gt 0) - $ $\Pr (({X_1} - F_1^{ - 1}(0.5))({X_2} - F_2^{ - 1}(0.5)) \lt 0)$ . It can be expressed as $4C(0.5,0.5) - 1$ , which makes it convenient to equate the theoretical and samples values to find the copula parameter. This method-of-moments approach is similar to that for Kendall’s tau.

The diagonal section of a copula is defined as $\delta (u) = C(u,u)$ (e.g. Nelsen, Reference Nelsen1999). If $C$ is the copula between ${U_1} \sim {\rm{Uniform}}\ (0,1)$ and ${U_2} \sim {\rm{Uniform}}\ (0,1)$ , $\delta (u)$ is the distribution function of $\max ({U_1},{U_2})$ . Accordingly, maximum likelihood estimation (MLE) can be performed using the density function ${\textstyle{\partial \over {\partial u}}}\delta (u)$ and the observed values $\max ({u_{1,j}},{u_{2,j}})$ . This method provides an alternative to the usual maximum likelihood approach involving the copula density function.

The three minimum distance estimators being explored here are based on the Cramér–von Mises and Kolmogorov–Smirnov test statistics (e.g. Genest et al., Reference Genest, Rémillard and Beaudoin2009). The following distances can be minimised to obtain the copula parameter:

\begin{align*}\sum\limits_{j = 1}^n {{{(\tilde Pr({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}}) - \tilde C({u_{1,j}},{u_{2,j}}))}^2}} \end{align*}

\begin{align*}\sum\limits_{j = 1}^n {|\tilde Pr({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}}) - \tilde C({u_{1,j}},{u_{2,j}})|} \end{align*}

\begin{align*}\mathop {\max }\limits_j |\tilde Pr({X_1} \le {x_{1,j}},{X_2} \le {x_{2,j}}) - \tilde C({u_{1,j}},{u_{2,j}})|\end{align*}

We apply the above-mentioned estimation methods to fit copulas 1, 4, 12, and 14 to the general liability claims data in Section 3. Table 13 lists the corresponding copula parameter estimates for each copula. It is interesting to see that the parameter estimates from different methods are very close to one another for the best-performing copula (4), while they vary a lot for the worst-performing copula (1). Besides providing initial values for more tedious estimation methods, applying different methods to the same set of data can provide a rough indication on whether the model being fitted is reasonable and robust.

Table 13. Copula parameter estimates calculated from different estimators.

We then conduct a simulation study to compare the performances of different estimation methods. We simulate 1,000 samples from each of copulas 1, 4, 12, and 14 using a tau value of 0.5 and then fit the same copulas to the simulated samples. Table 14 presents the estimated bias and standard error of the copula parameter under each estimation method based on 1,000 simulated scenarios. The usual MLE outperforms the other estimation methods in terms of both the bias and standard error. It is followed by the two method-of-moments approaches of matching Kendall’s tau and Blomqvist’s beta, in which their performances are close to that of the MLE. These two estimators have the advantage of being straightforward to compute, especially when Blomqvist’s beta is expressed directly in terms of the copula. Comparatively, applying the MLE to Archimedean copulas can be challenging when dealing with the required derivatives, particularly for multiple dimensions. The diagonal MLE comes next, the performance of which is not too far from those of the three methods above. For more tedious cases, the method-of-moments approaches and the diagonal MLE can be used to generate the initial values for the MLE. Lastly, the three minimum distance estimators clearly have larger bias and standard error than the other methods. In particular, the one based on the Kolmogorov–Smirnov test statistic performs the worst amongst the last three choices.

Table 14. Estimated bias (left) and standard error (right) of different estimators.

Hitherto, we have used fully parametric copulas and margins. In fact, there are some other semi-parametric and non-parametric approaches for modelling the dependency structure. One semi-parametric approach is to use a parametric copula, while estimating the margins non-parametrically. The use of the Kendall distribution function in Genest & Rivest (Reference Genest and Rivest1993) can be considered as semi-parametric as well. The Kendall distribution function is specified as $K(t) = \Pr (C({U_1},{U_2}) \le t)$ , which can be deduced as $K(t) = t - \varphi (t)/{\textstyle{\partial \over {\partial t}}}\varphi (t)$ for an Archimedean copula. Genest & Rivest (Reference Genest and Rivest1993) proposed a non-parametric estimator of $K(t)$ , which then provides an estimator of $C$ within the class of Archimedean copulas. The estimator of $K(t)$ can be seen as a decomposition of Kendall’s tau and is computed as ${\textstyle{1 \over n}}\sum\nolimits_{j = 1}^n {\lambda (t - {V_j})} $ , where $\lambda $ is the distribution function of a point mass at the origin and ${V_j} = {\textstyle{1 \over {n - 1}}}\# \{ ({U_{1,h}},{U_{2,h}})\,:\,{U_{1,h}} \lt {U_{1,j}},{U_{2,h}} \lt {U_{2,j}}\} $ . Using the general liability claims, Figure 8 compares the non-parametric estimate of $K(t)$ with the corresponding estimates of $K(t)$ of four Archimedean copulas. It can be seen that the non-parametric estimate highly resembles the function under copula 4 but deviates significantly from those of copulas 1 and 12. These observations are in line with the rankings in Table 3, which again suggest that upper tail dependence is an important feature of the data that cannot be overlooked.

Figure 8 Non-parametric estimate of Kendall distribution function versus the closest functions under copulas 1, 4, 5, and 12 (calculated from minimising the squared differences between the two curves).

A non-parametric approach takes both the copula and margins as parameter-free in order to provide the highest generality. One natural option for constructing an empirical copula would be the empirical bivariate distribution function. Recently, Segers et al. (Reference Segers, Sibuya and Tsukahara2017) introduced the empirical beta copula and empirical checkerboard copula which are two different “smoothed” versions and compared them with the traditional (unsmoothed) empirical copula. The empirical beta copula is given as ${\textstyle{1 \over n}}\sum\nolimits_{j = 1}^n {{F_{n,{\rm{Rank}}({x_{1,j}})}}({u_{1,j}}){F_{n,{\rm{Rank}}({x_{2,j}})}}({u_{2,j}})} $ , in which ${F_{n,r}}(u)$ $ = \sum\nolimits_{s = r}^n {\left( {\begin{array}{c}n\\s\end{array}} \right){u^s}{{(1 - u)}^{n - s}}} $ . On the other hand, the empirical checkerboard copula is specified as ${\textstyle{1 \over n}}\sum\nolimits_{j = 1}^n {\prod\nolimits_{i = 1}^2 {\min (\max (n{u_{i,j}} - {\rm{Rank}}({x_{i,j}}) + 1,{\rm{ }}0),{\rm{ }}1)} } $ . Under certain necessary and sufficient conditions, both are genuine copulas. These authors found that both copulas outperformed the traditional empirical copula in their simulation study.

In the following, we use the general liability claims data again and compare the estimates of expected shortfall under the three approaches: fully parametric (copula 4 or 6 with parametric margins), semi-parametric (estimator of $K(t)$ with parametric margins) and non-parametric (empirical or empirical beta or empirical checkerboard copula with empirical margins). Table 15 gives the expected shortfall estimates from taking each approach in turn. It can be seen that the fully parametric estimation produces the largest tail estimates. The semi-parametric approach based on Genest & Rivest (Reference Genest and Rivest1993) leads to similar results, which suggest that the estimate of $K(t)$ indicates an Archimedean copula with tail properties like those of copulas 4 and 6. Comparatively, the non-parametric estimation yields smaller estimates, the significance of which increases with the percentile level. The major implication of these results is that although a non-parametric approach can provide the highest generality for the data and is not restricted by the mechanics of a model, the resulting structure may lead to an underrepresentation of the potential impact of tail events and so a serious underestimation of tail measures. Recent regulatory developments have put an increasing focus on the assessment of extreme events. From an actuarial perspective, it is very important to have an adequate allowance for the tails, which have a significant impact on risk management and capital measurement.

Table 15. Estimated expected shortfall of total claim amounts under parametric, semi-parametric, and non-parametric approaches.

Note: Full I: copula 4 + parametric margins; Full II: copula 6 + parametric margins; Semi: estimator of $K(t)$ + parametric margins; Non I: empirical copula + empirical margins; Non II: empirical beta copula + empirical margins; Non III: empirical checkerboard copula + empirical margins.

7. Concluding Remarks

In this paper, we apply an extensive list of Archimedean copulas to some general and life insurance bivariate modelling problems and deal with both censoring and truncation issues. We cover not only the few copula choices that are commonly used in the literature, but also several others, as well as their rotated versions, which have not been fully tested in earlier applications. Our analysis clearly suggests that exploiting the rich diversity of Archimedean copula structures can provide a lot of flexibility for coping with different shapes of overall and tail dependence in different data sets. In particular, an adequate allowance for extreme tail events is of utmost importance for an insurer’s capital requirement. Archimedean copulas can serve as a very useful tool for calculating such capital allowance. We have shown that the copula assumption has a significant impact on the estimation of the expected shortfall. Moreover, we also highlight the usefulness of testing both the overall fitting and the tail fitting, in order to make a more careful copula selection. Furthermore, through a thorough simulation exercise, we notice some potential issues in copula selection and examine the practical implications. Tail dependence should be a major consideration when adopting any copula function in insurance modelling. Although it is a limiting property which is hard to identify precisely for real data, conducting a suitable statistical test on the tail association would help us choose a copula with more appropriate tail dependence properties, which affect the calculation of tail measures in pricing and capital management. Finally, we experiment with a number of other estimation methods including matching Blomqvist’s beta, maximum likelihood based on the diagonal section of a copula, minimum distance estimators, a non-parametric estimator for Archimedean copulas, and empirical copulas. It appears that some of them can serve as useful alternatives to the usual approaches adopted in previous applications.

Figure 9 A partially nested Archimedean copula and a D-vine copula optimally fitted to the mortality of England and Wales, Netherlands, Norway, and Sweden.

Despite the potential usefulness of Archimedean copulas, there have been some issues when extending the copulas to more than two dimensions. First, a multi-dimensional Archimedean copula is permutation-symmetric and the single parameter does not allow for different dependence levels for different pairs. Its exchangeability limits the use in multivariate modelling. Second, the traditional conditional simulation algorithm for multivariate Archimedean copulas is cumbersome and is not efficient enough for practical use. Some recent developments which attempt to address these problems can be tested for insurance applications in future research. For instance, hierarchical Archimedean copulas (e.g. Joe, Reference Joe2014; Okhrin & Ristig, Reference Okhrin and Ristig2014) and vine copulas (e.g. Brechmann & Schepsmeier, Reference Brechmann and Schepsmeier2013; Czado, Reference Czado2019) can readily be built from the less commonly used copulas in Table 1 for handling multivariate random variables. One feasible approach to construct hierarchical Archimedean copulas is to fit a bivariate copula to every possible pair, choose the pair with the highest dependence and convert it into a pseudo variable, and repeat the process iteratively with the remaining variables. Similarly, the sequential method can be used to construct vine copulas, in which the optimal structure in each tree is determined by finding the so-called maximum spanning tree that maximises cumulative pairwise dependencies. As an example, Figure 9 presents the optimal structures of a partially nested Archimedean copula and a D-vine copula for the mortality of the four regions studied in Section 4. The former adheres closely to the geographical locations of the four regions and is highly interpretable, while the latter provides some flexibility in using different bivariate copulas within the structure. In both cases, besides the traditional use of the overall tests on dependence (e.g. tau, AIC, BIC), a more specific testing on the tails, as demonstrated earlier, can also be incorporated into the construction process.

Moreover, considering the impact of the data size, a Bayesian approach, rather than the usual maximum likelihood, can be taken to fit the copulas. One way is to incorporate extra information via the prior distributions, such that additional references can be deduced from certain relevant data. Another possible modification is to specify a prior distribution for the selection amongst different copula candidates. The resulting posterior distribution would then be a hybrid between the copulas under consideration and blend their desirable properties into one structure. Lastly, other interesting developments such as quasi-copulas and semi-copulas (e.g. Durante & Sempi, Reference Durante and Sempi2015) and partition-of-unity copulas (e.g. Pfeifer et al., Reference Pfeifer, Mändle, Ragulina and Girschig2019) may also be explored for their potential use in insurance modelling.