# A non-parametric estimator for Archimedean copulas under flexible censoring
scenarios and an application to claims reserving ††thanks: Citation:
Michaelides, M., Cossette, H. & Pigeon, M., A non-parametric estimator for
Archimedean copulas under flexible censoring scenarios and an application to
claims reserving (2024). Preprint.

Marie Michaelides
Chaire Co-operators en analyse des risques actuariels
Département de mathématiques
Université du Québec à Montréal
Montréal
<EMAIL_ADDRESS>
&Mathieu Pigeon
Chaire Co-operators en analyse des risques actuariels
Département de mathématiques
Université du Québec à Montréal
Montréal
&Hélène Cossette
École d’actuariat
Université Laval
Québec

###### Abstract

With insurers benefiting from ever-larger amounts of data of increasing
complexity, we explore a data-driven method to model dependence within
multilevel claims in this paper. More specifically, we start from a non-
parametric estimator for Archimedean copula generators introduced by Genest
and Rivest (1993), and we extend it to diverse flexible censoring scenarios
using techniques derived from survival analysis. We implement a graphical
selection procedure for copulas that we validate using goodness-of-fit methods
applied to complete, single-censored, and double-censored bivariate data. We
illustrate the performance of our model with multiple simulation studies. We
then apply our methodology to a recent Canadian automobile insurance dataset
where we seek to model the dependence between the activation delays of
correlated coverages. We show that our model performs quite well in selecting
the best-fitted copula for the data at hand, especially when the dataset is
large, and that the results can then be used as part of a larger claims
reserving methodology.

##  1 Introduction

With the increasing ease in data collection and greater availability of
computational resources for financial institutions, actuarial practitioners
and academics face everyday larger quantities of complex data. All this data
grants insurers access to detailed information about their insureds and the
claims that are typically collected over long periods of time, allowing for a
better understanding of their portfolio and, among other, more accurate loss
modelling.

This increase in quantity, however, also leads to an increase in data
complexity. Actuaries must sort through large datasets to retrieve relevant
information and derive predictions for premiums or the development of claims,
using data that is more often than before unstructured, missing or untimely.
One particular complexity frequently encountered in insurance is the issue of
incomplete or censored data, arising when part of the information is missing
for a specific observation. This situation can arise, for example, in a model
for coupled lifetimes when death occurs for one of the individuals in the
study, as illustrated in [1]. Claims severities can be censored as well if the
insurer benefits from a reinsurance treaty that covers a predefined portion of
the losses. We can also consider the case where a given prescription period on
the payments made by the insurer has passed, thereby stopping further payments
towards the claimant even if the claim is still open. Multiple factors can
also cause censoring. In the example with the coupled lifetimes, observations
can be censored either by the death of one individual or if they simply revoke
their contract with the insurance company. Insurers can thus observe a
multitude of censoring scenarios in which one or more variables can be
censored or where multiple factors can cause censoring.

Additionally, a second complexity often encountered in insurance data is the
correlation between different - potentially censored - variables of interest
to the insurers. This dependence can take the form, for example, of temporal
dependence between the occurrence of a specific event and a claim’s
development delay as in [2], [3] and [4]. [5], [6] or [7] also illustrate
hierarchical dependence within multi-features claims. We can even find cross-
sectional dependence between losses and expenses, as in [8], or between a
claim’s lifetime and its final amount, as in [9].

To account for this, many authors in the actuarial literature have used the
framework of copulas, i.e., multivariate distribution functions whose
marginals are uniformly distributed on the unit interval. Copulas allow a
representation of the joint distribution function of multivariate data by
separating their marginal distributions from the inherent dependence
structure. Examples of their application abound in the statistical literature,
and since their introduction by [10], many authors have extensively studied
and analyzed their properties.

Among the different classes of copulas, [11] introduced the Archimedean
family, widely applied today in many fields beyond insurance. Contrarily to
other copula families, Archimedean copulas are not derived from Sklar’s
theorem, as highlighted in [12]. They are characterized instead by a function
$\phi(.)$, defined with one or more so-called dependence parameters, called
the generator of the copula. Thus, estimating $\phi(.)$ directly allows us to
identify the copula model most appropriate for a given set of multivariate
observations. For this purpose, [13] propose a non-parametric approach relying
on a study of the bivariate probability integral transformation for
Archimedean copulas. Their estimator is, however, only applicable to complete
data, seldom encountered in insurance, as mentioned earlier.

Although some authors have worked with [13]’s non-parametric estimator or
extensions of it, these contributions are still rare, especially in the
actuarial literature. Considering the large amounts of data available to
insurers and the highly diverse censoring schemes that these datasets can
contain, we aim in this paper to fill this gap by suggesting the use of a non-
parametric estimator of Archimedean copula generators applicable to various
censoring scenarios. For this purpose, we base ourselves on the approach
proposed in [13] and slightly modify it, using results from [14] and an
estimator for the join distribution proposed by [15] in survival analysis. To
the best of our knowledge, this estimator has never been used in an actuarial
setting nor in the application of a non-parametric estimator for Archimedean
copulas.

The motivation of our work is three-fold. First, we use a non-parametric
estimator that readily applies to large amounts of data. The approach we
follow in this paper allows insurers to benefit from the large datasets at
their disposition and let the data speak for itself to select the most
appropriate model. Second, we account for the complexity of this data by using
a methodology that encapsulates its multivariate features and any possible
censoring schemes, including censoring caused by multiple factors. Third, we
gather and provide methods to validate and assess the goodness-of-fit of our
estimation, thereby strengthening it.

We structure our paper as follows: first, we provide a brief literature review
of the main papers we rely on in our approach in Section 2. In Section 3, we
describe the statistical model for the non-parametric estimator using methods
borrowed from survival theory. Section 4 presents three model validation
approaches under different censoring scenarios. In Section 5, we test our
model within these scenarios using extensive simulation studies. Section 6
presents an application to a Canadian automobile claims dataset in which
policyholders benefit from multiple insurance coverages, and that provides a
good illustration of a more complex censoring scheme. To build a micro-level
claims reserving model, we use copulas to account for the underlying
dependence between coverages. We conclude our work in Section 7.

##  2 Literature review

In recent years, attention has been paid to Archimedean copulas, particularly
in finance and insurance. Allowing for very flexible dependence structures,
this class of copulas has become popular thanks to their simple form and
statistical tractability, among other interesting properties. Notably, many of
them allow high-dimensional dependence modeling using only one parameter
governing the strength of dependence. Although initially more widely used in
life insurance where they arise from frailty models, Archimedean copulas have
grown popular in non-life insurance as well.

Thanks to the large quantities of data more easily available, benefits can be
reaped from using a non-parametric estimator such as the one proposed by [13]
to directly retrieve the generator of Archimedean copula models. We discuss
two interesting examples of the application of this estimator in the case of
censored data. First, the work of [16] combines an estimator for the generator
from [14] to an estimator of the joint distribution proposed by [17] to model
the dependence between the allocated loss adjustment expenses and right-
censored losses. Their methodology is, however, only applicable to cases where
only one variable is subject to censoring, thereby limiting possible
replications to more flexible scenarios.

To the best of our knowledge, the only non-parametric approach for more
flexible censoring scenarios is the work proposed in [18]. In their paper, the
authors propose a discrete estimator and two smooth estimators for Archimedean
copulas and any type of copula applicable to various censoring scenarios.
Their discrete estimator extends the empirical copula estimator proposed in
[19] by using random weights in the empirical joint distribution to account
for the bias caused by censored observations. For the two smooth copula
estimators, [18] introduce kernel functions in addition to the weights of the
empirical joint distribution estimator.

In the present paper, we use an estimator situated between those of [16] and
[18]. As in [16], we focus only on Archimedean copulas and use a non-
parametric estimator specifically for Archimedean generators. We, however,
seek to allow for the same level of flexibility in the choice of censoring
scenarios as [18]. As such, we extend the work of [16] to the case where any
or all variables can be subject to censoring by replacing the estimator for
the joint distribution with that proposed in [15].

In addition, we want to use robust approaches to validate the results obtained
from such a non-parametric estimator. Goodness-of-fit tests for copulas have
been discussed by some authors in the literature, notably [20] who present a
review of available tests, with a focus on so-called blanket tests, i.e., that
are not constrained to one specific class of copulas or that do not require
excessive data manipulation, based on the empirical copula and Kendall
distribution. The authors also propose a new testing procedure based on
Rosenblatt’s transform. More recently, [12] reviewed all the latest
developments in copula theory, including estimation and testing methods.
However, all these tests only apply to complete data. The literature on
goodness-of-fit tests for copulas with censored data is much scarcer. For the
specific case of Archimedean copulas, [21] proposes a test based on a Cramer-
von-Mises type distance and Kendall distribution applicable to Archimedean
copulas. [22] proposes a test based on a multiple-imputation procedure of the
censored data for that same class of copulas, and [23] use a cross-ratio
function. Moving away from Archimedean copulas, [24] derive likelihood and
pseudo-likelihood tests, embedding the copulas under the null and alternative
hypotheses. [25] propose a general smooth test that assumes a specific form
for the alternative copulas. Others, such as [26], describe a goodness-of-fit
test easily applicable to high-dimensional copulas, and that allows for any
parametric copula specification. They propose a penalized pseudo-likelihood
ratio statistic to compare various semi-parametric multivariate survival
functions subject to copula misspecification and general censorship. Recently,
[27] and [28] presented general tests based on information matrices applicable
to various copula families. While [27] remains focused on right-censored data,
[28] allow for greater flexibility by including interval-censoring and
recurrent events in their test. Finally, in addition to their three non-
parametric copula estimators, [18] also propose a test similar to that of [21]
that can be used with different kinds of distances such as Kolmogorov-Smirnov,
Cramer-von-Mises or the square-root of a quadratic integrated distance.

In this paper, we use three approaches to validate the non-parametric
estimation of the Archimedean copula generator. More specifically, we use an
omnibus procedure based on the likelihood function, a method based on the
$L^{2}$-norm distance, inspired by [21] and [18], and the test proposed by
[22].

##  3 Statistical Model

In this section, we describe an approach to retrieve the Archimedean copula
generator $\phi(.)$ in the case of flexible censoring. More specifically, we
present a non-parametric estimator building on the work of [15] in survival
analysis. The advantage of using such a non-parametric estimator lies in the
possibility of directly tailoring the model to the data, even in the presence
of complex censoring. This contrasts with the more naive copula modelling
approach that usually consists of fitting a certain number of copulas and
comparing their performance using information criteria such as the Akaike
Information Criterion (AIC) or Bayesian Information Criterion (BIC).

We begin this section by defining the notation before describing the
estimator. We end by detailing the graphical comparison method, thanks to
which we identify the copula model most appropriate to a given dataset. For
clarity, we focus solely on the bivariate case, but the notation and model
described hereafter can easily be extended to a multivariate framework.

###  3.1 Notation

We define the following model components:

  * •

$\bm{T}=(T_{1},T_{2})$ with $T_{i}\in\mathbb{R}^{+}$ for $i=1,2$ is the
initial vector of variables of interest. They can represent times-to-event,
losses, etc. To avoid identifiability issues with the copula, we assume these
variables to be continuous;

  * •

$\bm{X}=(X_{1},X_{2})$ with $X_{i}\in\mathbb{R}^{+}$ for $i=1,2$ is the vector
of censoring variables, such that $\bm{T}$ and $\bm{X}$ are independent;

  * •

$\bm{\omega}=(\omega_{1},\omega_{2})$, with $\omega_{i}>0$ for $i=1,2$ are
constants acting as additional censoring factors on $\bm{T}$. We call these
the limits.

  * •

$\bm{Y}=(Y_{1},Y_{2})$ where $Y_{i}=\min(T_{i},X_{i},\omega_{i})$ for $i=1,2$
is the observed bivariate vector, taking censoring of $\bm{T}$ into account;

  * •

$\bm{\Delta}=(\mathds{1}_{[Y_{1}=T_{1}]},\mathds{1}_{[Y_{2}=T_{2}]})$ is the
vector of censoring indicators.

To illustrate the difference between the censoring variable $X$ and the limit
$\omega$, consider a study of the time to recovery after administration of a
treatment of patients suffering from a medical condition. The time to
recovery, denoted by $T$, is censored if the patient dies, with the time to
death being the censoring variable $X$. In addition, health professionals
might consider the treatment ineffective if the patient does not recover after
a given amount of time, even if he is still alive. Further treatments will
then need to be considered. In this example, the time elapsed until the
treatment is deemed to have failed is the limit $\omega$. Even if the patient
is still alive, i.e., $T>X$, the imposition of $\omega$ will censor his time
to recovery. The observed delay $Y=\min(T,X,\omega)$ is of interest to medical
practitioners. $Y$ represents the delay until the first event, being either
full recovery, death, or need for additional treatments.

Further applications can be found in finance, insurance, and other fields.
Section 6 focuses on an application to bivariate claims delays.

###  3.2 Non-parametric estimator of the generator

We now describe the non-parametric approach leading to the identification of
the copula linking the pairs of censored variables $(Y_{1},Y_{2})$.

Assuming flexible censoring scenarios where either one or both variables can
be subject to random censoring, we follow the methodology laid out in [15] and
work with the following extension of [29]’s estimators for the conditional
distributions of $Y_{1}$ (resp. $Y_{2}$) given $Y_{2}=y_{2}$ (resp.
$Y_{1}=y_{1}$):

$\displaystyle\hat{F}_{1|2}(y_{1}|y_{2})=1-\prod_{Y_{i1}\leq
y_{1},\Delta_{i1}=1}\Bigg{(}1-\frac{W_{ni2}(y_{2};h_{n})}{\sum_{j=1}^{n}W_{nj2}(y_{2};h_{n})\mathds{1}_{Y_{j1}\geq
Y_{i1}}}\Bigg{)},$ (1)

where $y_{2}$ must be uncensored and

$\displaystyle
W_{ni}(y;h_{n})=\begin{cases}\frac{k\Big{(}\frac{y-Y_{i2}}{h_{n}}\Big{)}}{\sum_{\Delta_{j1}=1}k\Big{(}\frac{y-Y_{j2}}{h_{n}}\Big{)}},&\hskip
8.5359pt\text{if}\hskip 5.69046pt\Delta_{i2}=1\\\ 0,&\hskip
8.5359pt\text{if}\hskip 5.69046pt\Delta_{i2}=0.\end{cases}$

Here, $k(.)$ is a known kernel function and $\\{h_{n}\\}$ is the bandwidth: a
sequence of positive constants such that $h_{n}\rightarrow 0$ as
$n\rightarrow\infty$. We work with a similar estimator for
$\operatorname{P}[Y_{2}\leq y_{2}|Y_{1}=y_{1}]=\hat{F}_{2|1}(y_{2}|y_{1})$.

The estimator of the joint distribution proposed in this context by [15] is
then given by:

$\displaystyle\hat{F}(\bm{y})=w(\bm{y})\int_{0}^{y_{2}}\hat{F}_{1|2}(y_{1}|z_{2})d\tilde{F}_{2}(z_{2})+(1-w(\bm{y}))\int_{0}^{y_{1}}\hat{F}_{2|1}(y_{2}|z_{1})d\tilde{F}_{1}(z_{1}),$
(2)

where $\tilde{F}_{1}$ and $\tilde{F}_{2}$ are the marginal estimators of
Kaplan and Meier (1958):

$\displaystyle\tilde{F}_{j}(y_{i})=1-\prod_{Y_{i,j}\leq
y_{i},\Delta_{ij}=1}\Big{(}1-\frac{1}{n-i+1}\Big{)},\hskip 8.5359ptj=1,2.$

The weights $w(\bm{y})$ minimise the mean-squared error of $\hat{F}(\bm{y})$.
Equation (2) is the estimator that, combined with [13]’s approach, allows us
to retrieve the generator of the copulas in the presence of different
censoring scenarios. We start from [13]’ proposition laid out just below.

Proposition 1.1. [[13]] Let $S_{1}(t_{1})$ and $S_{2}(t_{2})$ be marginal
survival functions whose dependence function $C(.)$ is of the form

$\displaystyle
C(S_{1}(t_{1}),S_{2}(t_{2}))=\phi^{-1}\Big{\\{}\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))\Big{\\}},$

where $\phi(.)\in(0,1]$ is a convex decreasing function such that $\phi(1)=0$.
Let

$\displaystyle
U=\frac{\phi(S_{1}(t_{1}))}{\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))}\hskip
8.5359pt\text{and}\hskip
8.5359ptV=\phi^{-1}\Big{\\{}\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))\Big{\\}}=C(S_{1}(t_{1}),S_{2}(t_{2})).$

Then,

  1. 1.

$U$ is uniformly distributed on $(0,1)$;

  2. 2.

$V$ follows the Kendall distribution defined as
$K(\nu)=\nu-\frac{\phi(\nu)}{\phi^{(1)}(\nu)}$, for $0<\nu\leq 1$; and

  3. 3.

$U$ and $V$ are independent.

[13] suggest the following estimator for the Kendall distribution in the
presence of complete data:

$\displaystyle\hat{K}_{n}(\nu)=\frac{1}{n}\\#\\{i|\nu_{i}\leq\nu\\},$

where $\\#$ denotes the cardinality of a set.

We cannot use this estimator in the presence of censored data. Instead, we
choose to apply the extension proposed by [14] to estimate the Kendall
distribution by

$\displaystyle\hat{K}_{n}(\nu)=\int_{0}^{\infty}\int_{0}^{\infty}\mathds{1}_{[\hat{F}(\bm{y})\leq\nu]}d\hat{F}(\bm{y}),$
(3)

in which we insert (2) as the estimator of the joint distribution, allowing
for all possible censoring scenarios.

Finally, we insert $\hat{K}_{n}(\nu)$ in the estimator for the Archimedean
generator proposed by [13]:

$\displaystyle\hat{\phi}_{n}(\nu)=\exp\Bigg{\\{}\int_{\nu_{0}}^{\nu}\frac{1}{t-\hat{K}_{n}(t)}dt\Bigg{\\}},$
(4)

with $0<\nu_{0}<1$ an arbitrarily chosen constant. Thanks to the use of
Equation (2), we can retrieve the generator non-parametrically, regardless of
the level of censoring present in the data.

###  3.3 Graphical comparison

Starting from the empirical estimator $\hat{K}_{n}(\nu)$, we illustrate the
graphical procedure introduced by [13] to select the best fitting model among
competing copulas for a given dataset. We define the univariate function
$\lambda(\nu)$ from which [13] recover the generator in Equation (4) as

$\displaystyle\lambda(\nu)=\frac{\phi(\nu)}{\phi^{(1)}(\nu)}.$

An estimator of $\lambda(\nu)$ is easily retrieved by

$\displaystyle\hat{\lambda}_{n}(z)=z-\hat{K}_{n}(z).$ (5)

The idea is now to graphically compare the empirical estimator
$\hat{\lambda}_{n}(\nu)$ to the equivalent functions
$\lambda_{\hat{\alpha}}(\nu)$ of each competing copula model under
consideration, where $\hat{\alpha}$ is the dependence parameter estimated with
(6). Table 1 presents some of the most popular Archimedean copulas and the
one-to-one relationship between their dependence parameter and Kendall’s tau.
The estimation of $\hat{\tau}$ is straightforward:

$\displaystyle\hat{\tau}=4\int_{0}^{1}\hat{\lambda}_{n}(\nu)d\nu=3-4\int_{0}^{1}\hat{K}_{n}(\nu)d\nu.$
(6)

Table 1: $\lambda_{\alpha}(\nu)$ for some popular Archimedean copulas.

Copula | $\tau$ |  $\lambda_{\alpha}(\nu$)
---|---|---
Clayton | $\alpha/(\alpha+2)$ | $\nu(\nu^{\alpha}-1)/\alpha$
Frank | $1+\frac{4}{\alpha}\Big{(}\int_{0}^{\alpha}\frac{\xi}{\alpha(e^{\xi}-1)}d\xi-1\Big{)}$ | $-\ln{\Big{\\{}(e^{-\alpha\nu}-1)/(e^{-\alpha}-1)}\Big{\\}}/(-\alpha/(e^{\alpha\nu}-1))$
Gumbel-Hougaard | $1-1/\alpha$ | $\nu\ln({\nu})/\alpha$
Joe | - | $-\ln{\Big{\\{}1-(1-\nu)^{\alpha}}\Big{\\}}/-\Big{(}(\alpha(1-\nu)^{(\alpha-1)})/(1-(1-\nu)^{\alpha})\Big{)}$

To illustrate the methodology laid out in this section, we consider the Loss-
ALAE dataset used in [8], [16], and [18]. We are in a single-censoring
scenario where only the loss variable is subject to right-censoring. We
estimate $\hat{K}_{n}(\nu)$ and $\hat{\lambda}_{n}(\nu)$ using, respectively,
Equation (3) and Equation (5). We then plot $\hat{\lambda}_{n}(\nu)$ and
$\lambda_{\hat{\alpha}}(\nu)$ for the Clayton, Frank, Gumbel, and Joe copulas.
We estimate $\hat{\alpha}$ for each using the one-to-one relationship between
Kendall’s tau, retrieved from Equation (6) using $\hat{\lambda}_{n}(\nu)$ and
the dependence parameters, as illustrated in Table 1.

Figure 1 compares the empirical estimate and the copulas under consideration.
We observe that when using the approach laid out in Section 3.2 with [15]’s
estimator in (2) for the joint distributions, we obtain very similar results
to those presented in [16] who used the estimator from [17]. More
specifically, we clearly see in Figure 1 that $\hat{\lambda}_{n}(\nu)$ closely
follows the curve $\lambda_{\hat{\alpha}}(\nu)$ of the Gumbel copula. This
indicates that the Gumbel copula best fits the Loss-ALAE data. We also show in
Table 2 the estimated Kendall’s tau and dependence parameters obtained in,
respectively, [16] using [17]’s estimator in the second column (single
censoring model), compared to those in the third column obtained using the
methodology from Section 3.2 (flexible censoring model). We observe that all
estimated parameters are very similar. This illustrates the adequacy of the
approach that we use in this paper.

Figure 1: Graphical comparison of different copulas for the Loss-ALAE dataset. Table 2: Comparison of parameters estimates. | Single censoring model | Flexible censoring model
---|---|---
Copula | $\tilde{\tau}=0.3567$ | $\hat{\tau}=0.3507$
Clayton | 1.1088 | 1.0803
Frank | 3.5919 | 3.5177
Gumbel-Hougaard | 1.5544 | 1.5403
Joe | 2.0074 | 1.9805

##  4 Results validation

While practical and easy to implement, the estimator described in Section 3.2
is non-parametric, and the choice of the copula model is driven at this stage
by a graphical analysis only, as illustrated in Figure 1. Goodness-of-fit
tests for copulas are rare in the literature and even more so in the presence
of censoring. To validate our choice and increase the level of confidence in
our approach, we gather three different approaches in this section: a pseudo-
maximum likelihood estimation procedure of the dependence parameters, a
validation based on the $L^{2}$-norm distance, and a statistical test for
Archimedean copulas with censored data introduced by [22].

###  4.1 Omnibus estimation procedure

The omnibus estimation procedure, or maximum pseudo-likelihood procedure, is a
semi-parametric optimization procedure that substitutes empirical versions of
the marginal distributions in the (parametric) likelihood function of the
copula model, $L(.)$.

We use rescaled Kaplan-Meier estimators for $T_{1}$ and $T_{2}$, multiplying
the original Kaplan-Meier versions by $n/(n+1)$, where $n$ is the total number
of observations in the sample. As mentioned in [30] as well as [16], this
rescaling prevents issues related to the potential unboundedness of the copula
log-density as $u_{1}$ and $u_{2}$ tend to one.

Once we have the marginal distribution functions, we use them in the
estimation step, i.e., the optimization of the model likelihood to obtain the
pseudo-likelihood estimator $\hat{\alpha}^{*}$ of the dependence parameter:

$\displaystyle\hat{\alpha}^{*}=\operatorname*{arg\,max}L(u_{1},u_{2},$
$\displaystyle\delta_{1},\delta_{2};\alpha),$

where the likelihood function under flexible censoring scenarios is given by:

$\begin{split}L(u_{1},u_{2};&\delta_{1},\delta_{2};\alpha)=\prod_{i=1}^{n}c_{\alpha}(u_{1i},u_{2i};\alpha)^{\delta_{1i}\delta_{2i}}+\bigg{(}\frac{\partial
C_{\alpha}(u_{1i},u_{2i};\alpha)}{\partial
u_{1}}\bigg{)}^{\delta_{1i}(1-\delta_{2i})}\\\ &+\bigg{(}\frac{\partial
C_{\alpha}(u_{1i},u_{2i};\alpha)}{\partial
u_{2}}\bigg{)}^{(1-\delta_{1i})\delta_{2i}}+C_{\alpha}(u_{1i},u_{2i};\alpha)^{(1-\delta_{1i})(1-\delta_{2i})}.\end{split}$
(7)

$C_{\alpha}$ is the candidate Archimedean copula under consideration,
$c_{\alpha}$ its density and $\delta_{ji}=1$ if $T_{j}<X_{j}$ for $j=1,2$ for
observation $i$ and equals zero otherwise (i.e., for a censored observation).
The partial derivatives with respect to either $u_{1}$ or $u_{2}$ can be found
in Table 3.

Table 3: Partial derivatives for some popular Archimedean copulas with $\tilde{u}=-\ln u$ and $\bar{u}=1-u$. Copula | $\partial C_{\alpha}(u_{1}.u_{2})/\partial u_{1}$
---|---
Clayton | $[1+u_{1}^{\alpha}(u_{2}^{-\alpha}-1)]^{-1-1/\alpha}$
Frank | $[e^{-\alpha u_{1}}-e^{-\alpha(u_{1}+u_{2})}]\times[(1-e^{-\alpha})-(1-e^{-\alpha u_{2}})(1-e^{-\alpha u_{1}})]^{-1}$
Gumbel | $u_{1}^{-1}e^{\\{-(\tilde{u}_{2}^{\alpha}+\tilde{u}_{1}^{\alpha})^{1/\alpha}\\}}[1+(\tilde{u}_{2}/\tilde{u}_{1})^{\alpha}]^{-1+1/\alpha}$
Joe | $(1-\bar{u}_{2}^{\alpha})(1-\bar{u}_{2}^{\alpha}+\bar{u}_{2}^{\alpha}\bar{u}_{1}^{-\alpha})^{-1+1/\alpha}$

In order to use the omnibus procedure to validate the model choice resulting
from the graphical comparison, we retrieve the pseudo-maximum likelihood
estimator $\hat{\alpha}^{*}$ for the chosen candidate copula models and
compare them to the estimators $\hat{\alpha}$ found with the non-parametric
model using (6) and Table 1. The copula with the smallest difference between
these two estimators is the most appropriate for the data.

###  4.2 $L^{2}$-norm

The second approach to validate the choice of copula model is based on the
$L^{2}$-norm. We follow the methodology laid out in [14] that introduce a
goodness-of-fit statistic specifically designed for Archimedean copulas with
censored data under the assumption that a consistent non-parametric estimator
of the bivariate joint distribution function is available. Note that with
estimated parameters present in the specification of the null hypothesis for
this test, the approach we show here is not an actual statistical test but
rather a way to validate our model. However, we will use the word test in this
section for simplicity. As mentioned in the introduction, more formal and
general methods have been proposed in the literature. We discuss one in
Section 4.3.

[14]’s test relies on the $L^{2}$-norm distance between the empirical
estimator $\hat{K}_{n}(\nu)$ and the corresponding $K_{\hat{\alpha}}(\nu)$ of
the hypothesized parametric copula model,

$\displaystyle
D(\hat{\alpha})=\int_{\xi}^{1}\big{(}\hat{K}_{n}(\nu)-K_{\hat{\alpha}}(\nu)\big{)}^{2}d\nu.$

We use Riemann sum approximations to simplify numerical analysis:

$\displaystyle
D(\hat{\alpha})=\sum_{i=1}^{n}\big{(}\hat{K}_{n}(\nu_{(i)})-K_{\hat{\alpha}}(\nu_{(i)})\big{)}^{2}(\nu_{(i)}-\nu_{(i-1)}),$

where $\nu_{(i)}$ is the $i^{\text{th}}$ ordered value of
$\\{\nu_{1},\nu_{2},...,\nu_{n}\\}$. However, the asymptotic distribution of
$D(\hat{\alpha})$ is difficult to retrieve analytically. We propose below a
parametric bootstrap procedure inspired by [21] and [18] to select the most
appropriate copula model based on a pseudo $p$-value:

Step 1. Compute $\hat{K}_{n}(\nu)$ for the data at hand and, for the $M$
candidate models under consideration, estimate $\hat{\alpha}_{m}$, for
$m=1,...,M$.
Step 2. For each candidate model $C_{\hat{\alpha}_{m}}$, generate $B$ censored
samples of size $n$.
Step 3. For each of the $B$ samples, estimate $\hat{\alpha}_{m,b}$ and
$D(\hat{\alpha}_{m,b})$.
Step 4. For each model $m$, obtain the pseudo $p$-value as

$\displaystyle
p_{m}=\frac{1}{B}\sum_{b=1}^{B}\mathds{1}[\min_{l}D(\hat{\alpha}_{b,l})>D(\hat{\alpha}_{b,m})]$

with $m\neq l$. The copula model $m$ best fitted to the data is the one for
which the pseudo $p$-value is the smallest. In this case, we can interpret
$p_{m}$ as the number of times among $B$ replications that copula model $m$ is
selected as most fitted for the data at hand.

###  4.3 Goodness-of-fit test for censored dependent data

The omnibus procedure and $L^{2}$-norm validation approach from Sections 4.1
and 4.2 are easy to implement and give a good idea of whether a copula model
is appropriate for the data at hand. However, neither are formal statistical
tests. To further increase the confidence in our choice of copula model, we
present a more formal approach for Archimedean copulas, proposed by [22].

This goodness-of-fit test applies to uncensored and right-censored data thanks
to a multiple imputation procedure. It bases its premise on Proposition 1.1.
from [13].

For a bivariate vector of dependent observations $(T_{1},T_{2})$ whose
marginal survival functions $S_{1}(.)$ and $S_{2}(.)$ have a dependence
function $C(.)$ of the form

$\displaystyle
C(S_{1}(t_{1}),S_{2}(t_{2}))=\phi^{-1}\Big{\\{}\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))\Big{\\}},$

then, the random variables

$\displaystyle
U=\frac{\phi(S_{1}(t_{1}))}{\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))}\hskip
8.5359pt\text{and}\hskip
8.5359ptV=\phi^{-1}\Big{\\{}\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))\Big{\\}}=C(S_{1}(t_{1}),S_{2}(t_{2})).$
(8)

are independent. In other words, under the null hypothesis that an Archimedean
copula can model the dependence between $(T_{1},T_{2})$ with generator
$\phi(.)$, the correlation coefficient between $U$ and $V$ is null:

$\displaystyle H_{0}:\rho=0\hskip 8.5359pt\text{vs}\hskip
8.5359ptH_{1}:\rho\neq 0.$

Let

$\displaystyle
r_{n}=\frac{\sum_{i=1}^{n}(\hat{U}_{i}-\bar{\hat{U}})(\hat{V}_{i}-\bar{\hat{V}})}{\sqrt{\sum_{i=1}^{n}(\hat{U}_{i}-\bar{\hat{U}})^{2}\sum_{i=1}^{n}(\hat{V}_{i}-\bar{\hat{V}})^{2}}}$

where $\hat{U}$ and $\hat{V}$ are consistent estimators of $U$ and $V$ and
where $\bar{\hat{U}}$ and $\bar{\hat{V}}$ are the sample means of
$\hat{U}_{i}$ and $\hat{V}_{i}$. The test statistic is defined as

$\displaystyle Z_{n}=\frac{1}{2}\log\Big{[}\frac{1+r_{n}}{1-r_{n}}\Big{]}$

and we have that $\sqrt{n}Z_{n}\rightarrow N(0,1)$ in distribution.

For censored data, we can not directly retrieve $\hat{U}_{i}$ and
$\hat{V}_{i}$ from the estimators proposed by [13] based on Equation (8).
Instead, [22] shows that, depending on the censoring pattern, the data can be
simulated using a multiple imputation procedure from one of the distributions
presented in A. All proofs and further details are in [22].

##  5 Simulation study

Using simulation studies, we now illustrate the approach described in Section
3. We provide in this section a few examples to showcase the graphical model
selection procedure and the validation approaches introduced in Section 4.

###  5.1 Independence case

The first way to validate the results from the estimator for the generator in
Section 3 is to analyze how it handles independent data. Setting Kendall’s tau
equal to 0 in Table 1, each of the four copulas should converge to the
independence copula. This is the case if $\alpha=0$ for both the Clayton and
Frank copulas and $\alpha=1$ for the Gumbel and Joe copulas.

We simulate samples of size $n=1000$ from each of the four copulas with
$\tau=0$. We then use the methodology laid out in Section 3 to estimate
$\hat{K}_{n}(\nu)$, $\hat{\tau}$ and, finally, $\hat{\alpha}$ for each copula.
Table 2 presents the results of 1000 simulations where we plot the estimated
dependence parameters. We also present the average estimates over the 1000
simulations in Table 4.

Figure 2 shows that the Clayton and Frank copulas simulations are clearly
centered around $0$, with the Frank copula presenting slightly larger
variations. These are, however, simply due to the wider range of possible
values that $\alpha$ can take for this specific copula. The simulations are
very close to 1 for the Gumbel and Joe copulas. Table 4 confirms these
results. Our approach using [15]’s estimator for the joint distribution
provides adequate estimates of the dependence parameters when $\tau=0$.

Figure 2: Estimated dependence parameters for different copulas in case of independence ($\tau=0$). Table 4: Average estimates of $\hat{\alpha}$ for $1000$ simulations when $\tau=0$. | Clayton | Frank | Gumbel | Joe
---|---|---|---|---
$\hat{\alpha}$ | 0.0249 | 0.0352 | 1.0163 | 1.0150

###  5.2 Graphical comparison

We now move away from the independence scenario and focus on the non-
parametric estimator of $\hat{\lambda}_{n}(\nu)$ used for the graphical
comparison.

We simulate bivariate samples $(T_{1},T_{2})$ of $n=1000$ observations from
the Clayton, Gumbel, Frank, and Joe copulas. In each case, we select the
dependence parameter such that Kendall’s tau equals $0.4$, and we work with
unit-exponential distributions for the marginals. We simulate the vector of
censoring variables $(X_{1},X_{2})$ from exponential distributions with their
parameters set such that $20\%$ of observations have at least one censored
component. We are in a double-censoring scenario where one or both variables
can be censored.

Figure 3 displays the empirical estimator $\hat{\lambda}_{n}(\nu)$ and the
curves of $\lambda_{\hat{\alpha}}(\nu)$ for the four copulas under
consideration, for each of the four simulated samples. In each plot, the
continuous curve depicts the non-parametric estimator presented in Equation
(5). For all, $\hat{\lambda}_{n}(\nu)$ appears to be closest to the curve of
the correct copula, i.e., the copula from which the sample was simulated. In
addition to the one from Figure 1, these graphical comparisons suggest that
the approach proposed in Section 3 performs quite well under flexible
censoring scenarios. In the remainder of this section, we propose further
validating these results using three different approaches.

Figure 3: Graphical comparison of the estimated $\hat{\lambda}(\nu)$ functions
for the four simulated samples. We simulate $n=1000$ bivariate observations
from a Clayton (top left), Frank (top right), Gumbel (bottom left) and Joe
(bottom right) copulas.

###  5.3 Omnibus procedure

We begin validating the results from Figure 3 with the omnibus procedure. In
addition to the double-censoring scenario considered here, we add two other
censoring schemes: one with complete data (no censoring) and one where we
allow for only one censored variable. We work again for these two additional
scenarios with simulated samples of size $n=1000$, Kendall’s tau equal to
$0.4$, and unit-exponential margins. For the single-censoring scenario, we
simulate the censoring variable following an exponential distribution such
that around $20\%$ of observations have a censored component, similar to the
double-censoring scheme.

For each of the three scenarios, Table 5 displays the results of the omnibus
procedure. For each simulated sample, we compare the dependence parameters
$\hat{\alpha}$ retrieved from the estimation of Kendall’s tau in Equation (6)
to the maximum likelihood estimator $\hat{\alpha}^{*}$. We observe that in
most cases, the distance between $\hat{\alpha}$ and $\hat{\alpha}^{*}$ is
minimal for the correct copula, i.e., the omnibus procedure leads to selecting
the copula model from which the data was originally simulated. We count three
occurrences where the procedure selects an alternative model. In the single-
censoring scenario, the Gumbel model is chosen for the Clayton sample, and the
Joe copula is selected for the Gumbel sample. The Gumbel copula is selected
for the Frank sample in the double-censoring scenario. However, we observe
that the correct copulas present in each case are the second smallest distance
between $\hat{\alpha}$ and $\hat{\alpha}^{*}$ and that the difference with the
third and fourth distances is quite marked. For the Gumbel sample in the
single censoring scenario for example, the smallest distance between the
parameters is equal to $0.0428$ for the Joe copula, $0.0992$ for the Gumbel
copula and then, $0.3988$ and $0.4971$ for, respectively, the Clayton and
Frank copulas.

To increase the confidence in model selection, we show in Table 6 the results
of performing the omnibus comparison between $\hat{\alpha}$ and
$\hat{\alpha}^{*}$ not once but one thousand times. Each value in this table
represents the percentage of simulations in which the competing models were
not considered best-fitted for the sample in terms of the difference between
the estimated dependence parameters. We observe, for example, that when we
simulate a sample from the Clayton copula with $\tau=0.4$ in a double-
censoring scenario, only $276$ simulations out of a thousand reject the
Clayton model as best-fitted. Among these $276$ simulations, $188$ opt for a
Joe model, $71$ for a Gumbel model, and $17$ for a Frank model.

The results displayed in Tables 5 and 6 illustrate that the omnibus procedure
is a good first indicator of the appropriateness of the copula model selected
via the graphical procedure.

Table 5: Omnibus procedure for different censoring scenarios.

True copula | Candidate model |  | No censoring |  | single-censoring |  | double-censoring
---|---|---|---|---|---|---|---
| $\hat{\alpha}$ | $\hat{\alpha}^{*}$ |  | $\hat{\alpha}$ | $\hat{\alpha}^{*}$ |  | $\hat{\alpha}$ | $\hat{\alpha}^{*}$
Clayton ($\alpha=1.3332$) | Clayton |  | 1.7694 | 1.8284 |  | 1.1789 | 1.3263 |  | 0.9240 | 0.8541
Frank |  | 5.2063 | 5.7575 |  | 3.7724 | 4.1154 |  | 2.1323 | 2.7655
Gumbel |  | 1.8846 | 1.7294 |  | 1.5895 | 1.4926 |  | 1.4621 | 1.3070
Joe |  | 2.6347 | 1.7373 |  | 2.0732 | 1.4713 |  | 1.5231 | 1.2792
Frank ($\alpha=4.1611$) | Clayton |  | 1.3202 | 0.8049 |  | 1.2055 | 0.7428 |  | 1.0034 | 0.5326
Frank |  | 4.1284 | 4.2527 |  | 3.8403 | 3.9827 |  | 3.3149 | 2.9126
Gumbel |  | 1.6602 | 1.5274 |  | 1.6029 | 1.4533 |  | 1.5018 | 1.1694
Joe |  | 2.2066 | 1.6870 |  | 2.0983 | 1.6260 |  | 1.9084 | 1.1759
Gumbel ($\alpha=1.6667$) | Clayton |  | 1.1377 | 0.6557 |  | 1.1380 | 0.7392 |  | 0.9464 | 0.6493
Frank |  | 4.1641 | 4.3241 |  | 3.6674 | 4.1645 |  | 2.8121 | 3.6805
Gumbel |  | 1.5690 | 1.5934 |  | 1.5691 | 1.6683 |  | 1.4734 | 1.5283
Joe |  | 2.0345 | 1.8735 |  | 2.0347 | 1.9919 |  | 1.8671 | 1.7554
Joe ($\alpha=2.2191$) | Clayton |  | 1.2609 | 0.6067 |  | 1.1494 | 0.5514 |  | 0.9344 | 0.5170
Frank |  | 3.9802 | 4.2203 |  | 3.6968 | 3.8131 |  | 3.1298 | 3.4508
Gumbel |  | 1.6305 | 1.7284 |  | 1.5748 | 1.6356 |  | 1.4674 | 1.5638
Joe |  | 2.151 | 2.1997 |  | 2.0454 | 2.0431 |  | 1.8441 | 1.9360

Table 6: Percentage of simulations in which different candidate copula models are selected with the omnibus procedure. Scenario |  | True copula |  | Candidate model
---|---|---|---|---
|  | Clayton | Frank | Gumbel | Joe
No censoring |  | Frank |  | 0.876 | 0.322 | 0.853 | 0.949
| Joe |  | 0.897 | 0.680 | 0.845 | 0.578
Single-censoring |  | Frank |  | 0.880 | 0.274 | 0.884 | 0.962
| Joe |  | 0.892 | 0.844 | 0.868 | 0.396
Double-censoring |  | Clayton |  | 0.276 | 0.983 | 0.929 | 0.822
| Frank |  | 0.865 | 0.392 | 0.947 | 0.796
| Gumbel |  | 0.906 | 0.992 | 0.370 | 0.732
| Joe |  | 0.936 | 0.942 | 0.962 | 0.160

###  5.4 $L^{2}$-norm

Next, we perform a second study in which we simulate bivariate samples of
$n=500$ observations from various copulas with Kendall’s tau equal to $0.2$,
$0.4$, and $0.6$. We use the same exponential distributions of parameters
equal to one for the marginal distributions of $(T_{1},T_{2})$. We again
consider three censoring schemes: one with no censoring, one where we allow
only one variable to be censored, and one where both variables are subject to
censoring. For the double-censoring scheme, we also analyze the independence
case.

When working with complete data, we use [13]’s estimator $\hat{K}_{n}(\nu)$ to
compute $D(\hat{\alpha})$:

$\displaystyle\hat{K}_{n}(\nu)=\frac{1}{n}\\#\\{i|\nu_{i}\leq\nu\\},$

with

$\displaystyle\nu_{i}=\frac{1}{n-1}\\#\\{(t_{1,(j)},t_{2,(j)})|t_{1,(j)}<t_{1,(i)},t_{2,(j)}<t_{2,(i)}\\}.$

For the single-censoring scenario, we use a [14]’s estimator
$\hat{K}_{n}(\nu)$ similar to Section 3.2, but with a different estimator for
the joint distribution $\hat{F}(\bm{y})$. Namely, we compute the estimator
proposed by [17] for single-censoring scenarios. Assuming that $T_{1}$ is the
variable that can be censored by $X_{1}$, we have:

$\displaystyle\hat{F}(y_{1},t_{2})$
$\displaystyle=\int_{0}^{t_{2}}\hat{F}_{Y_{1}|T_{2}}(y_{1}|z)d\hat{F}_{T_{2}}(z)$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}\mathds{1}[0\leq t_{2,k}\leq
t_{2}]\hat{F}_{Y_{1}|T_{2}}(y_{1}|t_{2,k}),$

where $\hat{F}_{Y_{1}|T_{2}}(y_{1}|z)$ is [29]’s estimator from Equation (1).

For the two censoring scenarios, we simulate the censoring variables using
exponential distributions with parameters such that at least $20\%$ of
observations have one (or two) censored components. For each sample, we then
perform the bootstrap procedure described in Section 4.2 to obtain estimates
of the $L^{2}$-norm distance between $\hat{K}_{n}(\nu)$ estimated non-
parametrically for the data and the corresponding $K_{\hat{\alpha}}(\nu)$ for
different candidate models.

Table 7 presents the results of $1000$ simulations for each sample. Each value
in this table corresponds to the pseudo $p$-value from Section 4.2 and should
be understood as the percentage of bootstrap simulations in which the
candidate model was not selected as the best model for the simulated sample.
That is the percentage of simulations for which the $L^{2}$-norm distance
between $\hat{K}_{n}(\nu)$ and $K_{\hat{\alpha}}(\nu)$ for the candidate model
was not the smallest one compared to the other models under consideration. We
use samples from both Frank and Joe copulas in the no-censoring and single-
censoring schemes. At the same time, in the double-censoring scenario, we
analyze samples from the Clayton, Frank, Gumbel, and Joe copulas.

We observe that in all three scenarios for each copula, the pseudo $p$-value
is the smallest for the candidate model corresponding to the true copula. For
example, in the double-censoring scenario, when we simulate the data from a
Frank copula with Kendall’s tau equal to $0.4$, we observe that the Frank
copula is rejected as an appropriate candidate model in just $1.5\%$ of the
bootstrapped simulations performed. The Gumbel model is rejected in $98.5\%$
of simulations, and the Clayton and Joe copulas are rejected in all
simulations. In addition, we observe that as the value of Kendall’s tau
increases, the percentage of rejection for the correct model decreases. For
the Frank sample in the double-censoring scenario, $8.7\%$ of the bootstrap
simulations reject the Frank model when $\tau=0.2$. This percentage decreases
to $1.5\%$ when $\tau=0.4$ and to $0.9\%$ when $\tau=0.6$. We observe similar
results and trends in all scenarios for each copula considered. This showcases
the strength of our test.

In addition, for the doubly-censored scenario, we simulate samples from each
of the four copulas with Kendall’s tau equal to $0$. This is equivalent to
simulating four bivariate samples from the independence copula. We observe
that the results of the bootstrap simulations are much more random. For the
Joe sample, for example, $45\%$ of simulations reject the Clayton model, all
reject the Frank model, $80\%$ reject the Gumbel model, and $75\%$ reject the
Joe model. This is to be expected since, when $\tau=0$, all four copulas
converge towards the same independence copula. The estimates that we thereby
get for $K_{\hat{\alpha}}(\nu)$ and $\lambda_{\hat{\alpha}}(\nu)$ all converge
towards the same functions. The choice of the smallest distance between
$\hat{K}_{n}(\nu)$ and $K_{\hat{\alpha}}(\nu)$ for the different candidate
models thus becomes much more random. This is well illustrated in Figure 4
that shows the plot of $\hat{\lambda}(\nu)$ and those of
$\lambda_{\hat{\alpha}}(\nu)$ for the four candidate models considered—the
curves for the Clayton, Frank, Gumbel, and Joe copula overlap.

Figure 4: Graphical comparison of the estimated $\hat{\lambda}(\nu)$ functions for the independent samples.  Table 7: Pseudo $p$-values for different candidate copula models under different censoring scenarios. Scenario |  | True copula | $\tau$ | Candidate model
---|---|---|---|---
| Clayton | Frank | Gumbel | Joe
No censoring |  | Frank | 0.2 | 1.000 | 0.091 | 0.909 | 1.000
| 0.4 | 1.000 | 0.067 | 0.948 | 0.985
| 0.6 | 1.000 | 0.001 | 0.999 | 1.000
| Joe | 0.2 | 1.000 | 0.999 | 0.924 | 0.077
| 0.4 | 1.000 | 1.000 | 0.998 | 0.002
| 0.6 | 1.000 | 1.000 | 0.999 | 0.001
single-censoring |  | Frank | 0.2 | 1.000 | 0.078 | 0.923 | 0.999
| 0.4 | 1.000 | 0.054 | 0.946 | 1.000
| 0.6 | 1.000 | 0.012 | 0.999 | 0.989
| Joe | 0.2 | 1.000 | 1.000 | 0.928 | 0.072
| 0.4 | 1.000 | 1.000 | 0.987 | 0.013
| 0.6 | 1.000 | 1.000 | 0.999 | 0.001
double-censoring |  | Clayton | 0.0 | 0.470 | 0.890 | 0.790 | 0.850
| 0.2 | 0.170 | 0.831 | 0.999 | 1.000
| 0.4 | 0.079 | 0.921 | 1.000 | 1.000
| 0.6 | 0.000 | 1.000 | 1.000 | 1.000
| Frank | 0.0 | 0.230 | 1.000 | 0.770 | 1.000
| 0.2 | 0.986 | 0.087 | 0.928 | 0.999
| 0.4 | 1.000 | 0.015 | 0.985 | 1.000
| 0.6 | 1.000 | 0.009 | 0.991 | 1.000
| Gumbel | 0.0 | 0.620 | 1.000 | 1.000 | 0.038
| 0.2 | 1.000 | 0.970 | 0.168 | 0.862
| 0.4 | 1.000 | 0.975 | 0.135 | 0.890
| 0.6 | 1.000 | 0.986 | 0.057 | 0.957
| Joe | 0.0 | 0.450 | 1.000 | 0.800 | 0.750
| 0.2 | 1.000 | 1.000 | 0.843 | 0.157
| 0.4 | 1.000 | 1.000 | 0.929 | 0.071
|  | 0.6 | 1.000 | 1.000 | 0.949 | 0.051

###  5.5 Goodness-of-fit test

Using the same censoring scenarios, marginal distributions, and copulas as in
Section 5.4 with Kendall’s tau equal to $0.2$, $0.4$, and $0.6$, we now apply
[22]’s goodness-of-fit test described in Section 4.3. We present the results
of 1000 simulations for simulated samples of size $n=200$ in Table 8. For the
single-censoring and double-censoring scenarios, we use $M=5$ imputed
datasets.

We observe similar results to those obtained by [22] and to those from Section
5.4. When dependence is weak, the test does not always accurately predict the
best model for the data. In particular, we see that in the single-censoring
scenario when the sample comes from the Frank copula with $\tau=0.2$, even
though the lowest percentage of rejection of the null hypothesis is correctly
observed when the candidate model is the Frank copula, this percentage is also
relatively low for the three alternative models. As already observed in [22],
the power of the test, however, increases when we increase the sample size. To
illustrate this, we present in Figure 5 what happens to the $p$-values in the
no-censoring scenario when the true model is the Frank copula with $\tau=0.4$
when we progressively increase the sample size from $n=200$ to $n=2000$. We
work with increasing steps of $200$ observations; in each case, we perform
1000 simulations. The strength of the test increases quite steeply with the
sample size. For the Frank candidate model, the resulting $p$-value starts at
a low level of approximately $5\%$, as was reported in Table 8, and quickly
tends to $0$ as we increase the sample size. The Gumbel copula displays a low
$p$-value of less than $10\%$ for the smallest sample size. However, it
greatly increases such that for a sample of $2000$ observations, it is
correctly rejected by over $85\%$ of simulations. The Joe and Clayton
candidate models display similar increasing trends, even though they started
with a higher $p$-value for small sample sizes than the Gumbel copula.

Similarly to the results of the $L^{2}$-norm in Table 7, we also observe in
Table 8 that the strength of the test increases with the level of dependence.
In all scenarios and for all copula models, the percentage of rejection of the
null hypothesis increases with $\tau$ when the null hypothesis is incorrectly
specified and decreases with $\tau$ when the true and null copula coincide.

Figure 5: Evolution of the $p$-values of competing models based on $1000$ simulations for varying sample sizes simulated from the Frank copula with Kendall’s tau equal to $0.4$ in a no censoring scenario. Table 8: Percentage of rejection of the null hypothesis for different copulas $n=200$ Scenario |  | True copula | $\tau$ | Copula under $H_{0}$
---|---|---|---|---
| Clayton | Frank | Gumbel | Joe
No censoring |  | Frank | 0.2 | 0.351 | 0.082 | 0.030 | 0.074
| 0.4 | 0.728 | 0.058 | 0.086 | 0.679
| 0.6 | 0.859 | 0.022 | 0.196 | 0.956
| Joe | 0.2 | 0.814 | 0.206 | 0.248 | 0.022
| 0.4 | 0.922 | 0.776 | 0.724 | 0.018
| 0.6 | 0.944 | 0.970 | 0.942 | 0.003
single-censoring |  | Frank | 0.2 | 0.142 | 0.032 | 0.229 | 0.142
| 0.4 | 0.712 | 0.084 | 0.800 | 1.000
| 0.6 | 0.760 | 0.000 | 0.968 | 0.990
| Joe | 0.2 | 0.720 | 0.030 | 0.052 | 0.016
| 0.4 | 0.970 | 0.232 | 0.081 | 0.007
| 0.6 | 0.980 | 0.400 | 0.131 | 0.004
double-censoring |  | Clayton | 0.2 | 0.097 | 0.031 | 0.512 | 0.239
| 0.4 | 0.021 | 0.203 | 0.835 | 0.945
| 0.6 | 0.016 | 0.645 | 0.968 | 0.958
| Frank | 0.2 | 0.391 | 0.022 | 0.057 | 0.029
| 0.4 | 0.810 | 0.012 | 0.850 | 0.999
| 0.6 | 0.944 | 0.001 | 0.975 | 1.000
| Gumbel | 0.2 | 0.456 | 0.082 | 0.010 | 0.012
| 0.4 | 0.848 | 0.172 | 0.030 | 0.051
| 0.6 | 0.962 | 0.458 | 0.025 | 0.775
| Joe | 0.2 | 0.990 | 0.468 | 0.535 | 0.081
| 0.4 | 1.000 | 0.530 | 0.727 | 0.002
| 0.6 | 1.000 | 0.880 | 0.843 | 0.000

###  5.6 Impact of the limit

In Section 3.1, we introduced a limit $\omega_{i}$ for $i=1,2$ acting as an
additional censoring layer. Depending on the context, these limits can be
imposed by practitioners, authorities, etc. They can be the maximum time
before treatment is considered ineffective, even if the patient is still
alive, or the maximum payment delay authorized for a claim for example.

In this section, we conduct a simulation study to assess the impact of
$\omega$ on the strength and form of dependence between two censored
variables.

Consider a bivariate vector of time-to-events $\bm{T}=(T_{1},T_{2})$ with log-
normal distributions such that $\mu_{T_{1}}=8$, $\mu_{T_{2}}=7$,
$\sigma_{T_{1}}=1$ and $\sigma_{T_{2}}=3$. We set the correlation $\rho$
between $T_{1}$ and $T_{2}$ at $0.35$. Let $\bm{X}=(X_{1},X_{2})$ be the
censoring vector with $X_{1}$ independent from $X_{2}$. Both $X_{1}$ and
$X_{2}$ are log-normally distributed, with parameters $\mu_{X_{1}}$ and
$\mu_{X_{2}}$ chosen to allow for different levels of censorship. We also set
$\sigma_{X_{1}}=1$ and $\sigma_{X_{2}}=1$. We perform three simulation
studies: one with low levels of censorship for both variables (approx. $2\%$
and $1\%$ for, respectively, $T_{1}$ and $T_{2}$), one with medium levels of
censoring (approx. $25\%$ and $40\%$) and one with higher levels of censorship
(approx. $65\%$ and $75\%$ for, respectively, $T_{1}$ and $T_{2}$). Finally,
consider the limits $\omega_{1}$ and $\omega_{2}$ such that we observe the
bivariate vector
$\bm{Y}=(Y_{1},Y_{2})=\big{(}\min(T_{1},X_{1},\omega_{1}),\min(T_{2},X_{2},\omega_{2})\big{)}$.

We use the strategy laid out in Section 3 to derive the joint distribution (2)
from Beran’s estimators. As in [15], we choose the weights $w(\bm{y})=0.5$ for
the sake of simplicity although better results could be obtained by optimizing
them. Using the estimated joint distribution, we then derive Kendall’s tau and
estimate $\hat{\lambda}_{n}(\nu)$ for a sample of $n=500$ observations. We use
different quantiles from the distributions of $Y_{1}$ and $Y_{2}$ to set the
values of the limits $\omega_{1}$ and $\omega_{2}$.

Figure 6 displays the results of decreasing the limit as well as the level of
censoring on the form of dependence via the function $\hat{\lambda}_{n}(\nu)$
and the strength of dependence via Kendall’s tau. We observe that neither the
limit nor the level of censoring seem to have an impact on the form of
dependence: in all cases, the shape of the generator function remains similar.
As we increase the level of censoring and the limit, we notice that the shape
of dependence becomes even more pronounced. This is in line with the results
from Sections 5.4 and 5.5 where we noted that, even though both tests were
able to correctly choose the most appropriate copula for a given sample and
low dependence levels, their strength increased with Kendall’s tau.

Figure 6: Plots of $\lambda(\nu)$ and Kendall’s tau for the joint distribution
using different limits and different levels of censoring. Low, medium and high
levels of censoring correspond to, approximately, $5\%$, $30\%$ and $75\%$.

##  6 Application to automobile insurance claims

In this section, we apply our proposed non-parametric approach to a recent
Canadian automobile insurance dataset that includes over $600,000$ claims that
occurred between January 2015 and June 2021. Each policy in force in this
dataset provides four different coverages to its holder, and each claim can
fall under one or more of these coverages.

We focus on the Accident Benefits and Bodily Injury coverages, which are the
most important cost-wise in the portfolio. Our goal is to model the dependence
between them using their activation delays. These are defined as the time
elapsed between the reporting date of the claim and the date at which the
insurer first triggers the coverage in the claims management system. This is
the date the insurer records that the claim falls under that specific
coverage. Focusing on the bivariate case, i.e., a policy providing two
coverages, we denote the vector of activation delays by
$\bm{T}=(T_{1},T_{2})$. If the claim settles before one of the coverages is
triggered, then the delay for this coverage is censored. The censoring
variable is the claim settlement delay we denote by $\bm{X}=(X_{1},X_{2})$.
Note that we have $X_{1}=X_{2}$ in this specific application since we assume
that a claim only settles once. In addition, policyholders may be unable to
receive compensation for a specific coverage if a certain amount of time has
passed since reporting. This limit on the delays can depend on the nature of
the claims, local or governmental regulations, or even on company-level rules
of the insurer. Both internal and external agents can set it. This is our
limit $\omega_{i}$, for $i=1,2$ introduced in Section 3.1. Suppose a coverage
has not been triggered yet when this limit passes. In that case, its
activation delay becomes censored, even if the claim remains open.

Figure 7 illustrates the typical development of a claim and the notation
introduced above. In this example, a claim occurred on June $1^{\text{st}}$
and was reported the following day. On June $5^{\text{th}}$, the insurer
records that the claim has impacted coverage 1. Recording the delays in days,
the activation delay for this coverage is then $T_{1}=3$. Assuming a limit of
$\omega_{1}=\omega_{2}=30$ days on the activation delays, coverage 2 becomes
censored on July $2^{\text{nd}}$, even if the claim is still open. It settles
two days later, leading to $X_{1}=X_{2}=32$. The observed delays are then
given by
$(Y_{1},Y_{2})=\big{(}\min(T_{1},X_{1},\omega_{1}),\min(T_{2},X_{2},\omega_{2})\big{)}=\big{(}\min(3,32,30),\min(\infty,32,30)\big{)}=(3,30)$.

Figure 7: Typical claim development and illustration of the notation used.

Considering the observed activation delays shown (in periods of six months) in
Table 9 for the Accident Benefits and Bodily Injury coverages in the portfolio
at hand, we set the limit at $\omega_{1}=\omega_{2}=730$ days, equivalent to
two calendar years. Almost $94\%$ of all Accident Benefits claims trigger that
coverage within six months of reporting, and most of the remaining claims will
be labelled as Accident Benefits no later than a year after the reporting
date. The activation delays are longer for Bodily Injury claims, with $1.69\%$
of them taking two years or more after reporting to be labelled as such.

Table 9: Percentage of claims with different activation delays, shown in periods of six months. | Activation delays
---|---
Coverage | No delay | 1 period | 2 periods | 3 periods |  $\geq 4$ periods
Accident Benefits | 93.84 | 5.73 | 0.29 | 0.08 | 0.06
Bodily Injury | 85.86 | 9.86 | 1.46 | 1.13 | 1.69

We first apply the approach described in Section 3 to graphically find the
copula best fitted to the dataset at hand. We present the plots of the non-
parametric estimators $\hat{K}_{n}(\nu)$ and $\hat{\lambda}_{n}(\nu)$ in
Figure 8, along with the corresponding fitted curves for four candidate
Archimedean models. In both plots, the non-parametric estimator is closest to
the curves of the Joe copula. With an estimated $\hat{\tau}=0.2705$, the Joe
copula with $\hat{\alpha}=1.6652$ best fits our data.

The omnibus procedure first confirms this result. We estimate the dependence
parameters for the competing models using the likelihood in Equation (7). The
results shown in the third and fourth columns of Table 10 indicate that the
Joe copula is most appropriate since it presents the smallest difference
between $\hat{\alpha}$ and $\hat{\alpha}^{*}$. Similarly, the results of 1000
bootstrapped simulations for the $L^{2}$-norm validation approach shown in the
fifth column of Table 10 show that $D(\hat{\alpha})$ is the smallest for the
Joe copula. We finally apply [22]’s goodness-of-fit test and find again in the
last column of Table 10 that the Joe model presents the smallest percentage of
rejection when we perform 1000 simulations for each competing model.

Figure 8: $K(\nu)$ and $\lambda(\nu)$ for the different copulas and for the Canadian automobile insurance dataset with a limit at 730 days. Table 10: Results validation, based on 1000 bootstrapped simulations for the $L^{2}$-norm and 1000 simulations for [22]’s test. Copula | $\hat{\tau}$ | Omnibus procedure |  $L^{2}$-norm | [22]
---|---|---|---|---
$\hat{\alpha}^{*}$ | $\hat{\alpha}$ | $D(\hat{\alpha})$ |  $\%$ rejection $H_{0}$
Clayton | 0.2705 | 0.2432 | 0.7417 | 0.00399 | 0.9994
Frank | 1.1166 | 2.5612 | 0.00180 | 0.9265
Gumbel | 1.0554 | 1.3708 | 0.00084 | 0.8742
Joe | 1.3821 | 1.6652 | 0.00033 | 0.0476

Taking the resulting Joe copula as best fitted to our data, we now show the
results of some simulations for the activation delays, taking the
$1^{\text{st}}$ of January 2019 as the valuation date. Figure 9 shows the
simulated densities for the delays of both coverages. The dotted lines in both
plots represent the average of the simulations, while the continuous line
shows the true mean delays. We observe that the simulations provide a good
match for the observed data. The simulated mean activation delay for the
Accident Benefits coverage is 134.79 days compared to an observed average
delay of 131.52. The average simulated delay for the Bodily Injury claims is
257.26 days, and the observed average is slightly higher at 271.08 days. This
larger difference in the predictions for the Bodily Injury claims stems from
the greater variability typically observed for this coverage.

When a claim occurs, we can thus predict the activation delays for both
coverages with a rather good level of accuracy. This can then be incorporated
into a larger claims reserving model and enhance the predictions for the total
portfolio reserve.

Figure 9: Densities of the simulated activation delays for the Accident
Benefits and Bodily Injury coverages, taking the $1^{\text{st}}$ of January
2019 as valuation date.

##  7 Conclusion

In this paper, we present a non-parametric estimator for the generator
function of the Archimedean family of copulas and extend it to flexible
censoring scenarios. Following the initial idea of [13], we propose a simple
graphical selection procedure that works well under various levels of
censoring. We then present three approaches to validate the results, all
applicable to flexible censoring schemes for the Archimedean family.

We illustrate both the non-parametric approach and validation methods using
various simulation studies. We discuss the impact of a limit we introduce as a
particular type of censoring and show that it does not affect the model
selection. We then apply this methodology to a real automobile insurance
dataset using activation delays for two insurance coverages provided within
each policy in force. We show good results when using real, large datasets and
can easily retrieve the most appropriate copula model. We further show that
the copula model resulting from the non-parametric approach provides accurate
estimates of the activation delays for both coverages.

The approach presented in this paper could be used as part, for example, of a
claims reserving model such as that proposed in [31] or [32] where the
estimation of the activation delays can replace, respectively, the activation
patterns and the multinomial estimation of the claim type. A similar
methodology could also be used in different models such as that proposed in
[33] where, for example, the next event’s exact time for reported but not
settled claims could be extended to a multivariate framework by considering
simultaneously, as in Section 6, different business lines or insurance
coverages. Alternative dependence models could also be considered,
particularly other families of copulas, as in [18].

## References

  * [1] N.W. Deresa, I. Van Keilegom, and K. Antonio.  Copula-based inference for bivariate survival data with left truncation and dependent censoring.  Insurance: Mathematics and Economics, 107:1–21, 2022.
  * [2] X. Zhou and X.B. Zhao.  Applying copula models to individual claim loss reserving methods.  Insurance: Mathematics and Economics, 46(2):10, 2010.
  * [3] P. Shi, X. Feng, and J.-P. Boucher.  Multilevel modeling of insurance claims using copulas.  The Annals of Applied Statistics, 10(2):834–863, 2016.
  * [4] P. Shi and L. Yang.  Pair copula constructions for insurance experience rating.  Journal of the American Statistical Association, 113:122–133, 2018.
  * [5] E. Frees and E. Valdez.  Hierarchical insurance claims modeling.  Journal of the American Statistical Association, 103(484):1457–1469, 2008.
  * [6] E. Frees, P. Shi, and E. Valdez.  Actuarial applications of a hierarchical insurance claims model.  ASTIN Bulletin, 39(1):165–197, 2009.
  * [7] L. Yang and P. Shi.  Multiperil rate making for property insurance using longitudinal data.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 182(2):647–668, 2019.
  * [8] E.W. Frees and E. A. Valdez.  Understanding relationships using copulas.  North American Actuarial Journal, 2, 1998.
  * [9] Olivier Lopez.  A censored copula model for micro-level claim reserving.  Insurance: Mathematics and Economics, 87(4):1–14, 2019.
  * [10] Abe Sklar.  Fonctions de répartition à $n$ dimensions et leurs marges.  Publ. Inst. Statist. Univ. Paris 8, pages 229–231, 1959.
  * [11] B. Schweizer and A. Sklar.  Probabilistic Metric Spaces. North-Holland Series in Probability and Applied Mathematics.  North-Holland Publishing Co., New York, 1983.
  * [12] J. Größer and O. Okhrin.  Copulae: An overview and recent developments.  Wiley Interdisciplinary Reviews: Computational Statistics, 2021.
  * [13] C. Genest and L. Rivest.  Statistical inference procedures for bivariate archimedean copulas.  Journal of the American Statistical Association, 88:1034–1043, 1993.
  * [14] W. Wang and M.T. Wells.  Model selection and semiparametric inference for bivariate failure-time data.  Journal of the American Statistical Association, 95:62–72, 2000.
  * [15] M.G. Akritas and I. Van Keilegom.  Estimation of bivariate and marginal distributions with censored data.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 65:457–471, 2003.
  * [16] M. Denuit, O. Purcaru, and I. Van Keilegom.  Bivariate archimedean copula models for censored data in non-life insurance.  Journal of Actuarial Practice 1993-2006, 5, 2006.
  * [17] M.G. Akritas.  Nearest neighbor estimation of a bivariate distribution under random censoring.  Annals of Statistics, 22:1299–1327, 1994.
  * [18] S. Gribkova and O. Lopez.  Non-parametric copula estimation under bivariate censoring.  Scandinavian Journal of Statistics, 42, 2015.
  * [19] P. Deheuvels.  La fonction de dépendance empirique et ses propriétés. un test non paramétrique d’indépendance.  Acad. Roy. Belg. Bull. Cl. Sci. (5), 65(6):274–292, 1997.
  * [20] C. Genest, B. Rémillard, and D. Beaudoin.  Goodness-of-fit tests for copulas: A review and a power study.  Insurance: Mathematics and Economics, 4:199–213, 2009.
  * [21] M.L. Lakhal-Chaieb.  Copula inference under censoring.  Biometrika, 97:505–512, 2010.
  * [22] A. Wang.  Goodness-of-fit tests for archimedean copula models.  Statistica Sinica, 20:441–453, 2010.
  * [23] T. Emura, C.-W. Lin, and W. Wang.  A goodness-of-fit test for archimedean copula models in the presence of right-censoring.  Computational statistics and data analysis, 54(12):3033–3043, 2010.
  * [24] Y.E. Yilmaz and J.F. Lawless.  Likelihood ratio procedures and tests of fit in parametric and semiparametric copula models with censored data.  Lifetime Data Analysis, 17(3):386–408, 2011.
  * [25] J. Lin and X. Wu.  A diagnostic test for specification of copulas under censorship.  Econometric Reviews, pages 1–17, 2020.
  * [26] X. Chen, Y. Fan, D. Pouzo, and Z. Ying.  Estimation and model selection of semiparametric multivariate survival functions under general censorship.  Journal of Econometrics, 157:129–142, 2008.
  * [27] Q.M. Zhou.  Information matrix equivalence in the presence of censoring: A goodness-of-fit test for semiparametric copula models with multivariate survival data.  arXiv:2109.09782, 2021.
  * [28] T. Sun, Y. Cheng, and Y. Ding.  An information ratio-based goodness-of-fit test for copula models on censored data.  Biometrics, 79(3):1713–1725, 2023.
  * [29] R. Beran.  Nonparametric regression with randomly censored survival data, 1981.
  * [30] C. Genest, K. Ghoudi, and L.P. Rivest.  A semiparametric estimation procedure of dependence parameters in multivariate families of distributions.  Biometrika, 82:543–552, 1995.
  * [31] M. Michaelides, M. Pigeon, and H. Cossette.  Individual claims reserving using activation patterns.  European Actuarial Journal, 2023.
  * [32] M-P. Côté, C. Genest, and D.A. Stephens.  A bayesian approach to modeling multivariate multilevel insurance claims in the presence of unsettled claims.  Bayesian Anal., 17(1):67–93, 2022.
  * [33] K. Antonio and R. Plat.  Micro-level stochastic loss reserving for general insurance.  Scandinavian Actuarial Journal, 2014(7):649–669, 2014.

##  Appendix A Distributions of $U$ and $V$ proposed in [22] for diverse
censoring patterns

The following results are based on Corollary 1., Corollary 2. and Theorem 4.
from [22].

Let $(T_{1},T_{2})$ be a bivariate vector submitted to censoring by the
bivariate continuous vector $(X_{1},X_{2})$, as described in Section 3.1. Let
$S_{1}(t_{1})$ and $S_{2}(t_{2})$ be the marginal survival functions of,
respectively, $T_{1}$ and $T_{2}$, whose dependence function can be modelled
by an absolutely continuous Archimedean copula such that

$\displaystyle
C(S_{1}(t_{1}),S_{2}(t_{2}))=S(t_{1},t_{2})=\phi^{-1}\Big{\\{}\phi(S_{1}(t_{1}))+\phi(S_{2}(t_{2}))\Big{\\}}.$

Then, we have that:

  1. 1.

if both observations are censored, the distribution function of
$(V|T_{1}>x_{1},T_{2}>x_{2})$ is

$\displaystyle
F_{1}(v,x_{1},x_{2})=\frac{1}{S(x_{1},x_{2})}\Bigg{[}v-\frac{\phi(v)-\phi\\{S(x_{1},x_{2})\\}}{\phi^{\prime}(v)}\Bigg{]},\hskip
14.22636pt0\leq v\leq S(x_{1},x_{2})$

and the distribution function of $(U|T_{1}>x_{1},T_{2}>x_{2})$ is uniformly
distributed on the interval

$\displaystyle\Bigg{[}\frac{\phi\\{S_{1}(x_{1})\\}}{\phi(v)},1-\frac{\phi\\{S_{2}(x_{2})\\}}{\phi(v)}\Bigg{]};$

  2. 2.

if only the second observation is censored, the distribution function of
$(V|T_{1}=t_{1},T_{2}>x_{2})$ is

$\displaystyle
F_{2}(v,t_{1},c_{2})=\frac{p^{\prime}(\phi(v))}{p^{\prime}(\phi(S(t_{1},x_{2})))},\hskip
14.22636pt0\leq v\leq S(t_{1},x_{2})$

with $p(.)=\phi^{-1}(.)$ and the distribution function of
$(U|T_{1}=t_{1},T_{2}>x_{2})$ is

$\displaystyle
G_{2}(y,t_{1},x_{2})=\frac{p^{\prime}\\{\phi(S_{1}(t_{1}))/u\\}}{p^{\prime}\\{\phi(S(t_{1},x_{2}))\\}},\hskip
14.22636pt0\leq u\leq\frac{\phi\\{S_{1}(t_{1})\\}}{\phi\\{S(t_{1},x_{2})\\}};$

  3. 3.

if only the first observation if censored, the distribution function of
$(V|T_{1}>x_{1},T_{2}=x_{2})$ is

$\displaystyle
F_{3}(v,x_{1},t_{2})=\frac{p^{\prime}(\phi(v))}{p^{\prime}(\phi(S(x_{1},t_{2})))},\hskip
14.22636pt0\leq v\leq S(x_{1},t_{2})$

and the distribution function of $(U|T_{1}>x_{1},T_{2}=t_{2})$ is

$\displaystyle
G_{3}(y,x_{1},t_{2})=\frac{p^{\prime}\\{\phi(S_{2}(t_{2}))/(1-u)\\}}{p^{\prime}\\{\phi(S(x_{1},t_{2}))\\}},\hskip
14.22636pt\frac{\phi\\{S_{1}(x_{1})\\}}{\phi\\{S(x_{1},t_{2})\\}}\leq u\leq
1;$

  4. 4.

if both observations are uncensored, the distribution function of
$(V|T_{1}=t_{1},T_{2}=t_{2})$ is the same as in Proposition 1.1. from [13] and
the conditional distribution of $(U|V=v)$ is just a uniform distribution on
$[0,1]$ independent of $V$.