Olivier Lopez
email: olivier.lopez0@upmc.fr.

Keywords: Kaplan-Meier estimator, bivariate censoring, left-truncation, survival analysis, lifetime analysis, nonparametric estimation, copula models, bootstrap

A generalization of Kaplan-Meier estimator for analyzing bivariate mortality under right-censoring and left-truncation with applications to model-checking for survival copula models

Olivier Lopez

To cite this version:

Olivier Lopez. A generalization of Kaplan-Meier estimator for analyzing bivariate mortality under right-censoring and left-truncation with applications to model-checking for survival copula models. 2012.

Introduction

Pension contracts with reversion clauses depend on two random lifetimes: the one of the policyholder, and the one of his/her wife/husband. While studying the lifetime of a single individual is a classical problem for insurers, bivariate mortality analysis is quite more delicate. The easiest solution would be to study each lifetime separately, assuming that these two random variables are independent. However, this independence assumption may be far to be true when we are considering two members of the same couple. Several approaches have been proposed to model such lifetimes jointly, see e.g. [START_REF] Frees | Annuity valuation with dependent mortality[END_REF], [START_REF] Carriere | Bivariate survival models for coupled lives[END_REF], [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF]. These approaches are essentially based on parametric or semiparametric models (typically copula models), and therefore rely on rather strong assumptions. On the other hand, in this setting there is a serious lack of nonparametric estimators that may be used, for example, to perform model checking. The main reason for this absence stands in the diculty to handle bivariate truncation and bivariate censoring, which compose the main specicities of lifetime datas. The main aim of this paper is to dene a new nonparametric estimator of the joint probability law of a random vector of lifetimes, under a particular censoring and truncation scheme which corresponds to the typical situation encountered in insurance data-sets. As an example of application, we consider the Canadian data-set initially studied by [START_REF] Frees | Annuity valuation with dependent mortality[END_REF] .

In the univariate case, Kaplan-Meier estimator [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF], generalized to the case of left-truncation by [START_REF] Tsai | A note on the product-limit estimator under right censoring and left truncation[END_REF]) plays this role of a purely nonparametric estimator of the survival function which can be used as a benchmark for any parametric or semiparametric model. However, in the bivariate case, such an estimator of the joint survival function is more dicult to nd. The most classical one can be found in [START_REF] Dabrowska | Kaplan-Meier estimate on the plane[END_REF]. The main drawback of this estimator stands in the fact that it does not dene a true probability denition, making it hard to use when it comes to estimate any other quantity than the survival function itself. In fact, Dabrowska's estimator assigns negative mass to some observations (see Pruitt (1991a)) which is, at least, counterintuitive. Alternative estimators usually present similar drawbacks. [START_REF] Akritas | Estimation of bivariate and marginal distributions with censored data[END_REF] and [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF], under bivariate censoring, proposed dierent kind of estimators that happen to dene true distribution. However, the estimator of [START_REF] Akritas | Estimation of bivariate and marginal distributions with censored data[END_REF] requires a careful choice of a smoothing parameter, while the one dened by [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF] require additional assumptions on the censoring

The author wishes to thank the Society of Actuaries, through the courtesy of Edward J. Frees and Emiliano Valdez, for allowing use of the data in this paper variables which may not be realistic for insurance data-sets. In the present paper, we restrain ourselves to a particular situation which corresponds to the situation observed for the Canadian data-set. The two main assumptions that we make is that censoring occurs at the same time for the two members of the couple, and that the only cause of censoring is the fact that one member (at least) of the couple is still alive at the end of the statistical study. Using these two assumptions allows to dene a simple estimator of the joint distribution. We show that this estimator is uniformly n 1/2 -consistent for estimating the survival function and a large number of related quantities (such as Kendall's τ coecient for instance). Applications to parametric estimation, model checking and bootstrap are considered. The rest of the paper is organized as follows. In section 2, we dene the specic bivariate truncation and censoring scheme under which we work in the following. We use the assumptions of the model to dene our nonparametric estimator of the joint distribution. In Section 3, we provide some theoretical properties of this estimator, giving some uniform convergence rate. We particulary focus on the estimation of quantities which may be used to perform copula modelling and tting of the two lifetimes, and provide a bootstrap procedure that manage to take account for the presence of censoring and truncation. A detailed application of these techniques to the Canadian data-set is performed in section 4.

A nonparametric estimator of the joint distribution of two lifetimes

The aim of this section is to dene a nonparametric estimator for the joint distribution of two lifetimes, under a right-censoring and left-truncation model which is presented in section 2.1. Our estimator is based on a partial derivative equation which is described in section 2.2, while the estimator itself is presented in section 2.3. Comparison with the construction of Kaplan-Meier estimator in the univariate case, and with alternative estimators, is done in section 2.4.

2.1

Bivariate right-censoring and left-truncated model

In the following, we consider two lifetimes (T, U ), and i.i.d. replications (T i , U i ) 1≤i≤n of these random variables. We assume that these random variables are absolutely continuous with respect to Lebesgue's measure on the plane. A bivariate right-censoring model occurs when the random variables (T i , U i ) 1≤i≤n are not directly observed, that is when only the following data is available,

{ Y i = inf(T i , C i ), & δ i = 1 T i ≤C i Z i = inf(U i , D i ), & γ i = 1 U i ≤D i ,
where (C i , D i ) 1≤i≤n consists of i.i.d. replications of a random bivariate censoring vector (C, D), and (δ i , γ i ) 1≤i≤n are indicator functions allowing the distinction between censored and uncensored observations. In addition to right-censored observation, we consider bivariate left-truncation. Introducing two random variables (µ, ν), the random vector

(Y, Z, δ, γ) is observed only if Y ≥ µ and Z ≥ ν. Therefore, observations are made of (Y i , Z i , δ i , γ i , µ i , ν i ) 1≤i≤n where Y i ≥ µ i and Z i ≥ ν i .
In order to be able to estimate the joint distribution of (T, U ), some identiability assumptions are required. In the univariate case, classical assumptions consists of assuming the independence of T from (C, µ), and of C from (T, µ), and of assuming that the probability of observing non-truncated observations is non-zero (otherwise, the set of observations is almost surely empty). Throughout this paper, we will use a multivariate extension of such assumptions, which is summarized in Assumption 1 below.

Assumption 1 Assume that

1. (T, U ) is independent from (C, D, µ, ν), 2. (C, D) is independent from (T, U, µ, ν), 3. P(Y ≥ µ, Z ≥ ν) > 0.
It is important to mention that no additional assumption is made under the random variables (C, D), or (µ, ν), which may have a dependence structure that we do not need to specify.

If we consider the application of such a model in the study of the joint distribution of the lifetimes of two members of the same couple, T (resp. U ) will denote the total lifetime of the husband (resp. his wife), C (resp. D) will denote the age at which the husband (resp. the wife) stops being under observation for any other cause than death. Causes for exiting the statistical study may be various. If one of the members of the couple is not dead at the end of the statistical study, his/her lifetime is censored, but other causes may exist, especially in cases where people can surrender their contract (and therefore stop being under observation). Regarding the truncation variables, µ (resp. ν) will denote the age of the husband (resp. his wife) at the beginning of the statistical study (or at the time when he joins the study after its beginning). Of course, Y and Z, which correspond to the age at which the individuals leave the study for whatever cause, can only be observed if they are both larger than the ages at which the two members of the couple join the study.

Our aim is to provide a nonparametric estimation of the distribution of the random vector (T, U ). In absence of truncation or censoring, this task can be easily performed by using the empirical distribution, which is the distribution dened by the multivari-

ate empirical distribution function, F emp (t, u) = n -1 ∑ n i=1 1 T i ≤t,U i ≤u , which estimates F (t, u) = P(T ≤ t, U ≤ u).
Indeed, this standard procedure consists of assigning the mass n -1 at each observation, since, in an non-censored and non-truncated framework, each observation is expected to bring the same amount of information. In our model, F emp is of course unavailable, since it would require complete observation of (T, U ). We propose to rely on an estimator of the following form,

F (t, u) = n ∑ i=1 W i,n 1 Y i ≤t,Z i ≤u .
(2.1)

An advantage of such a representation of our estimator F is that there is a quick connection between such an estimator and the underlying distribution he denes. Indeed, W i,n can be interpreted as the mass assigned to observation i. The weight W i,n will be designed to compensate the lack of information caused by censoring and truncation in some parts of the distribution. Based on an estimator of the type (2.1), it becomes easy to estimate general quantities such as E[ϕ(T, U )] = ∫ ∫ ϕ(t, u)dF (t, u) for some function

ϕ, by computing ∫ ∫ ϕ(t, u)d F (t, u) = n ∑ i=1 W i,n ϕ(Y i , Z i ). (2.2)
Estimation of the joint distribution or survival function is only a special case of this estimation problem. Our theoretical results will show that we can obtain consistency of integrals of the type (2.2) under appropriate conditions.

To properly dene the weights W i,n , we will proceed in two steps. Section 2.2 provides a key relationship that will only allow to dene an ideal (but unfortunately uncomputable)

estimator. Nevertheless, this relationship will allow us to dene a true estimator of the type (2.1), which achieves some kind of asymptotic equivalence with the ideal one, as it is shown in section 3.1.

2.2

A key relationship for our bivariate lifetime model

In order to state the key result of our approach, let us rst introduce some notation.

Dene

S F (t, u) = P (T > t, U > u) , V (t, u) = P (Y ≤ t, Z ≤ u, δ = γ = 1|Y ≥ µ, Z ≥ ν) , H(t, u) = P (µ ≤ t ≤ Y, ν ≤ u ≤ Z|Y ≥ µ, Z ≥ µ) .
The following proposition establishes a relationship between these three functions and dF (t, u).

Proposition 2.1 Under Assumption 1,

dF (t, u) = S F (t, u)dV (t, u) H(t, u) = {∫ t + ∫ u + dF (t ′ , u ′ ) } dV (t, u) H(t, u) , (2.3)
or, equivalently, for any

ϕ ∈ L 1 , E [ϕ(T, U )] = ∫ ∫ ϕ(y, z)S F (y, z)dV (y, z) H(y, z) = E [δγw(Y, Z)ϕ(Y, Z)|Y ≥ µ, Z ≥ ν] ,
with w(y, z) = S F (y, z) H(y, z) .

(2.4) Before proving Proposition 2.1, let us mention how it allows us to dene an ideal

estimator of F. Dening W * i,n = n -1 δ i γ i w(Y i , Z i ), consider F * (t, u) = n ∑ i=1 W * i,n 1 Y i ≤t,Z i ≤u .
(2.5)

Based on Proposition 2.1, we see that

∫ ∫ ϕ(t, u)dF * (t, u) → ∫ ∫ ϕ(y, z)w(y, z)dV (y, z) = E [ϕ(T, U )] , a.s.,
from the strong law of large numbers, when the sample size increases. Of course, F * can not be computed, since it would require the knowledge of function H and of function

S F (t, u) = ∫ t + ∫ u + dF (t ′ , u ′ ).
Nevertheless, as announced, F * has to be understood as some benchmark estimator, which presents the advantage of depending only on

(Y i , Z i , δ i , γ i ) 1≤i≤n
and not on the unobservable variables (T i , U i ) 1≤i≤n .

Let us mention that F * does not necessarily denes a true probability distribution, since there is no reason than

∑ n i=1 W * i,n = 1. We will denote as W * ∞,∞ = 1 - ∑ n i=1 W i,n
which can be interpreted as the additional mass that has to be put at innity (that is (∞, ∞)) in order to achieve a total mass equal to 1 (let us note that W * ∞,∞ can even be negative, in the case where the total mass of F * is larger than 1). Similar phenomenon occur when we look at Kaplan-Meier estimator in the univariate case. If the largest observation is uncensored, then the total mass of Kaplan-Meier estimator is one, otherwise, an additional mass is missing that can be assumed to be put at innity. We now state the proof of Proposition 2.1.

Proof of Proposition 2.1. Dene

m * (ϕ) = ∫ ∫ ϕ(y, z)S F (y, z)dV (y, z) H(y, z) .
This quantity is equal to

E [δγw(Y, Z, µ, ν)ϕ(Y, Z)|Y ≥ µ, Z ≥ ν] = E [1 T ≤C,U ≤D w(T, U, µ, ν)ϕ(T, U )|Y ≥ µ, Z ≥ ν] .
Next, the right-hand side can be rewritten in the following way,

m * (ϕ) = E [1 T ≤C,U ≤D 1 T ≥µ,U ≥ν w(T, U, µ, ν)ϕ(T, U )] P(Y ≥ µ, Z ≥ ν) = E [E [1 T ≤C,U ≤D 1 T ≥µ,U ≥ν w(T, U, µ, ν)ϕ(T, U )|T, U, µ, ν]] P(Y ≥ µ, Z ≥ ν) = E [E [1 T ≤C,U ≤D |T, U, µ, ν] 1 T ≥µ,U ≥ν w(T, U, µ, ν)ϕ(T, U )] P(Y ≥ µ, Z ≥ ν) .
We now use the fact that (C, D) is independent from (T, U, µ, ν), which leads to

E [1 T ≤C,U ≤D |T, U, µ, ν] = S G (T -, U -), where S G (t, u) = P(C > t, D > u) is the joint survival function of (C, D). Next, we get m * (ϕ) = E [E [1 T ≥µ,U ≥ν |T, U ] S G (T -, U -)w(T, U, µ, ν)ϕ(T, U )] P(Y ≥ µ, Z ≥ ν) ,
and the independence of (µ, ν) from (T, U )

leads to E [1 T ≥µ,U ≥ν |T, U ] = L(T, U ), where L(t, u) = P(µ ≤ t, ν ≤ u)
. The proof will be complete if we show that w dened by (2.4)

is equal to S G (t-, u-) -1 L(t, u) -1 P(Y ≥ µ, Z ≥ ν).
To this aim, observe that, by denition,

H(t, u) = P(µ ≤ t ≤ Y, ν ≤ u ≤ Z|Y ≥ µ, Z ≥ ν) = P (Y ≥ t, Z ≥ u) P (µ ≤ t, ν ≤ u) P (Y ≥ µ, Z ≥ ν) ,
where we used the independence between (Y, Z) and (µ, ν).

Observing that P (Y ≥ t, Z ≥ u) = S G (t-, u-)S F (t-, u-)
from the independence between (T, U ) and (C, D), we nally get

H(t, u) = S G (t-, u-)S F (t-, u-)L(t, u)P(Y ≥ µ, Z ≥ µ)
, and the result follows, since S F (t-, u-) = S F (t, u) for absolutely continuous variables.

A computable estimator

In the previous section, we showed that the distribution function F was solution of the partial derivative equation (2.3). We use the term of partial derivative equation, since, in the case where F corresponds to an absolutely continuous measure, , u). This equation depends on two unknown quantities which are

dF (t, u) = f (t, u)dtdu, where f (t, u) = ∂ t ∂ u F (t
V and H. Nevertheless, V and H can be estimated empirically with n 1/2 -consistent estimators, that is

V (t, u) = 1 n n ∑ i=1 δ i γ i 1 Y i ≤t,Z i ≤t , Ĥ(t, u) = 1 n n ∑ i=1 1 µ i ≤t≤Y i ,ν i ≤u≤Z i .
Therefore, we propose to dene our estimator F as a solution of an empirical version of

equation (2.3), that is F solution of d F (t, u) = { ∫ t + ∫ u + d F (t ′ , u ′ ) } d V (t, u) Ĥ(t, u) . (2.6)
with the condition that the total mass of F is equal to one, in order to dene a true probability distribution. Since V is a sub-distribution function which puts mass only at observations such that δ i = γ i = 1, F is of the form (2.1), with W i,n = 0 when δ i = 0 or γ i = 0. This simple form will allow us to show existence and unicity of this function F , and will also provide an easy way to compute it.

Indeed, in this case, solving (2.6) is equivalent to determine a vector

W = (W 1,n , ..., W n,n , W ∞,∞ )
, where, again, W ∞,∞ represents the eventual additional mass that as to be put at innity. Equation (2.6) can then be rewritten in a matricial form, as

(I n+1 -Â)W = 0, (2.7)
where I n+1 denotes the identity matrix of R n+1 , and  = (â i,j ) 1≤i,j≤n , with ân+1,n+1 = 1, and, for

(i, j) ∈ {1, ..., n} 2 , âi,j = δ j γ j δ i γ i 1 Y j ≥Y i ,Z j ≥Z i n Ĥ(Y i , Z i ) , and âi,n+1 = δ i γ i n -1 Ĥ(Y i , Z i ) -1 , âj,n+1 = 0. Let us observe that there is no problem of denition for coecients âi,j , since Ĥ(Y i , Z i ) > 0, for all i.
Theorem 2.2 ensures existence of solutions of (2.7). Moreover, the proof furnishes a way to solve explicitely this equation.

Theorem 2.2 Under Assumption 1, there exists a unique solution of (2.7) such that:

1. W ∞,∞ > 0, 2. W i,n = 0 if δ i = 0 or γ i = 0, 3. ∑ n+1 i=1 W i,n = 1. Moreover, W i,n > 0 for all i such that δ i = γ i = 1.
Proof. Let m be the number of doubly uncensored observations. Since the weight of all other observations is assumed to be zero, solving (2.7) under the restriction that

W i,n = 0 if δ i = 0 or γ i = 0 is equivalent to solving (I m+1 -Ã) W = 0,
where à is a matrix of size (m + 1) × (m + 1), where we removed all the lines and columns corresponding to the censored observations, and W is the vector of all weights for doubly uncensored observations. Moreover, we can reorder the terms in W with respect to the values of Y (the rst component of W corresponds to the weight that is aected to the doubly uncensored observation with the smallest value of Y, and so on until the m-th component which corresponds to the largest uncensored value of Y ). We will denote by

Y (1) < . . . < Y (m
I m+1 -à = I m+1 -        0 . . . . . . . . . 0 . . . . . . . . . . . . . . . 0 . . . 0 . . . 0 1        =        1 . . . . . . . . . 0 . . . . . . . . . . . . . . . 1 . . . 0 . . . 0 0        , (2.8) since the i-th component of à W corresponds to n -1 Ĥ(Y (i) , Z (i) ) -1 ∑ j Wj 1 j∈R i , where R i = {j : Y (j) > Y (i) , Z (j) > Z (i) } ∪ {∞}. Therefore, we see that R i ⊂ {i + 1, ..., n} ∪ {∞},
which leads to the triangular form of à with zeros on the diagonal, except for the last coecient corresponding to the mass at innity. It follows from (2.8) that the dimension of the kernel of I m+1 -à is exactly one, showing the existence and unicity of a non-zero solution of (2.7) satisfying conditions 2 and 3 of the statement of Theorem 2.2.

We will now prove that the coecients of W are strictly positive if we impose W ∞,∞ > 0. To this aim, it suces to consider a solution of (2.7) with W ′ ∞ = 1, the conclusion being deduced after normalization by the sum of the coecients W ′ i . This proof can be done recursively, and furnishes an explicit solution of (2.7). Indedd,

W ′ m = n -1 Ĥ(Y (m) , Z (m) ) -1 W ′ ∞ is strictly positive since W ′ ∞ is. Since R m-1 ⊂ {m} ∪ {(∞, ∞)}, since Ĥ is a positive function and since W ′ m > 0 and W ′ ∞ > 0, we can conclude that W ′ m-1 > 0.
The same procedure can be applied to show that W ′ m-2 > 0, ..., W ′ 1 > 0. Moreover, we see that these weights can be computed recursively, beginnning with the mass at innity, since the mass at an observation only depends on the mass aected to the observations with a largest value of Y. The proof is complete by taking

W i,n = Wi,n ( ∑ Wj,n ) -1 ,
and

W ∞ = W∞ ( ∑ Wj,n ) -1 .
Remark 2.3 Adding a mass at innity is necessary if we wish to ensure existence of solutions of (2.7). However, as it will appear in the proof of our Theorem 3.1, we will see that this mass will be quite close to zero in practice (of order n -3/2 ). Nevertheless, there exists a particular situation in which we can slightly modify our estimator in order to reallocate the small weight of (∞, ∞) at the "largest" observation. Consider the case where there exists some observation i 0 such that

Y i 0 = max 1≤i≤n (Y i ) and Z i 0 = max 1≤i≤n (Z i ).
In this case (Y i 0 , Z i 0 ) can be considered as the largest observation in a bivariate sense, since all observations are smallest with respect to both Y and Z. If this observation is uncensored, that is δ i 0 = γ i 0 = 1, we can slightly modify equation (2.7) by simply putting the residual mass at point (Y i 0 , Z i 0 ) instead of putting it at (∞, ∞). The new system of n equations then becomes,

W i,n = n -1 Ĥ(Y j , Z j ) -1 n ∑ i=1 W j,n 1 Y j >Y i ,Z j >Z i , f or i ̸ = i 0 , W i 0 ,n = W i 0 ,n ,
and the same proof as for Theorem 2.2 applies, by searching a solution of the system with

W i 0 ,n = 1 (instead of W ∞,∞ = 1).

Comparison with other estimators

Comparison to Kaplan-Meier estimator.

In the univariate case, Kaplan-Meier estimator [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF]) is used to estimate the distribution of right-censored random variable, and [START_REF] Tsai | A note on the product-limit estimator under right censoring and left truncation[END_REF] proposed an adaptation which is used for dealing with left-truncation. This estimator satises a dierential equation which has to be connected with our Proposition 2.1. Indeed, this estimator F (t) of F (t) = P(T ≤ t) can be rewritten as

F (t) = n ∑ i=1 W in 1 Y i ≤t .
Looking at the expression of the weight W in in Sánchez [START_REF] Sánchez Sellero | Uniform representation of product-limit integrals with applications[END_REF], we see that this wait can be expressed as

W in = δ i (1 -F (T i )) Ĥ(T i ) ,
where

H(t) = P(µ ≤ t ≤ T )
, which is the one-dimensional counterpart of (2.4).

Comparison to other existing estimators in the case of bivariate censoring.

Most of the estimators used for studying bivariate censoring only focus on the estimation of the joint survival function. A consequence is that many of them do not dene regular distributions. The estimator proposed by [START_REF] Campbell | Large-sample properties of nonparametric bivariate estimators with censored data[END_REF] is not a proper survival function due to nonmonotonicity. Similarly, Dabrowska's estimator [START_REF] Dabrowska | Kaplan-Meier estimate on the plane[END_REF]) assigns negative mass to some points in the plane, as it is shown in Pruitt (1991a). As a consequence, since Dabrowska's does not dene a probability distribution, in can not be used for perfoming bootstrap. Moreover, assigning negative mass to some observations is, at most, counterintuitive. More generally, dening the mass that some estimator aects to each observation is generally dicult for previously existing estimators, making them dicult to use when it come to estimate integrals with respect to the distribution of (T, U ). See also [START_REF] Prentice | Covariance and survivor function estimation using censored multivariate failure time data[END_REF], [START_REF] Gill | Inecient estimators of the bivariate survival function for three models[END_REF] for modications of these estimators that suer from the same drawbacks.

Another class of estimators relies on nonparametric maximum likelihood techniques (NPMLE). The computation of the estimator is usually delicate, since it is dened implicitely (see [START_REF] Pruitt | Identiability of bivariate survival curves from censored data[END_REF], [START_REF] Pruitt | Strong consistency of self-consistent estimators: general theory and an application to bivariate survival analysis[END_REF] and van der Laan (1996)), while the convergence rate can be lower than n 1/2 . This problem of slow convergence rates appears in [START_REF] Akritas | Estimation of bivariate and marginal distributions with censored data[END_REF], which is one of the rare estimator being a true survival function. The estimator proposed by [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF], achieves this usual convergence rate and denes a true distribution, but relies on a restrictive assumption on the joint law of the censoring, that can be dicult to justify in the framework of insurance contracts.

As we will show in the following, the convergence rate of the estimator that we propose in the present paper is n 1/2 , while its computation is quite simple (and can be done without any approximation). Moreover, the fact that it denes a true distribution allows to use it to perform bootstrap procedure, as it is shown in section 3.3.

3 Theoretical behaviour of the new estimator and bootstrap procedure

The aim of this section is to show that our new estimator of the distribution is n 1/2 -consistent.

In section 3.1, we provide a general result for integrals of the type

∫ ∫ ϕ(t, u)d F (t, u).
Some particular cases are considered in section 3.2. Finally, in section 3.3, we show how a bootstrap procedure can be implemented in order to evaluate estimation errors.

3.1 n 1/2 -consistency of the new estimator

To prove n 1/2 -consistency, we compare the estimator F proposed in the previous section to the ideal estimator F * dened in (2.5). Indeed, estimator F * is very simple to study, since it is a sum of i.i.d. quantities. Therefore, the most delicate result that we need to prove is contained in Theorem 3.1 below, which aims to compare the mass W i,n to the ideal one W * i,n .

Theorem 3.1 Under Assumption 1, assume that there exists a constant K such that

1 H(Y, Z) < K, a.s. (3.1) then sup 1≤i≤n |W i,n -W * i,n | + |W ∞,∞ -W * ∞,∞ | = O P (n -3/2 ).
The proof of this result is postponed to the Appendix section. In this statement, let us discuss (3.1). This assumption has to be linked with the diculty to estimate the tail of the distribution (in a bivariate sense, this means when y and z tend to innity).

A particular case in which this assumption holds is when the censoring has a support included in some compact set [0, t max ]×[0, u min ], strictly included in the support of (T, U ), with P((C, D) = (t max , u min )) ̸ = 0. However, as it will be shown in the simulations we conduct in section 3.2.2, the fact that this assumption does not hold does not seem to have a strong impact in practice. Some similar problems occur when studying the asymptotics of Kaplan-Meier estimator in the univariate case. Many existing results are only valid on compact subsets strictly included in the support of the distribution (see e.g. Van Keilegom

and [START_REF] Van Keilegom | Transfer of tail information in censored regression models[END_REF], [START_REF] Du | Uniform strong representation of the conditional Kaplan-Meier process[END_REF], [START_REF] Lopez | Nonparametric estimation of the multivariate distribution function in a censored regression model with applications[END_REF]), because of similar problems of denominators close to zero that are delicate to handle from a theoretical point of view.

With additional work, these conditions can be replaced by moments conditions (which, unfortunately, can not be checked in practice), such as the one provided by [START_REF] Stute | The central limit theorem under random censorship[END_REF]). However, developing similar arguments in order to relax condition (3.1) is quite delicate in our framework (without signicant practical impact), and is therefore left to future research.

With at hand Theorem 3.1, we can now state our main result on integrals of the form ∫ ∫ ϕ(t, u)d F (t, u).

Theorem 3.2 Let F denote a class of functions such that, for all ϕ ∈ F, ϕ(∞, ∞) = 0,

and |ϕ(y, z)| ≤ Φ(y, z) with E[Φ(Y, Z)] < ∞.
Then, under the assumptions of Theorem 3.1,

∫ ∫ ϕ(t, u)d F (t, u) = ∫ ∫ ϕ(t, u)dF * (t, u) + R n (ϕ), (3.2) with sup ϕ∈F |R n (ϕ)| = O P (n -1/2 ). Moreover, if sup ϕ∈F 1 n n ∑ i=1 δ i γ i S F (Y i , Z i )ϕ(Y i , Z i ) H(Y i , Z i ) -E[ϕ(Y, Z)] = O P (n -1/2 ), (3.3)
we get

sup ϕ∈F ∫ ∫ ϕ(t, z)d F (t, z) - ∫ ∫ ϕ(t, z)dF (t, z) = O P (n -1/2 ).
(3.4) Condition (3.3) will hold as long as F is a Donsker class of functions (see [START_REF] Van Der Vaart | Weak convergence and empirical processes with applications to statistics[END_REF] for a precise dention), that is if F satises an uniform central theorem property. In particular, the class of indicator functions (T, U ) → 1 T ≤t,U ≤u where (t, u) belong to some set T satises condition (3.3), allowing to obtain uniform consistency of F as an estimator of F (see section 3.2.1 for more details).

Proof of Theorem 3.2. Write

∫ ∫ ϕ(t, u)d F (t, u) = ∫ ∫ ϕ(t, z)dF * (t, z) + n ∑ i=1 (W i,n -W * i,n )ϕ(Y i , Z i ).
The second term is bounded by (

n -1 n ∑ i=1 Φ(Y i , Z i ) ) × sup i n|W i,n -W * i,n |. Next, n -1 ∑ n i=1 Φ(Y i , Z i ) = O P (1), while sup i n|W i,n -W * i,n | = O P (n -1/2
), from Theorem 3.1, which shows that (3.2) holds. Finally, (3.4) can be obtained from Assumption (3.3).

3.2

Examples of applications of Theorems 3.1 and 3.2

We now present some applications of the general results of the previous section. The three examples that we give are motivated by survival copula models. Indeed, the bivariate survival function S F (t, u) = P(T > t, U > u) of the random vector (U, V ) admits, by Sklar's Theorem [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF]), a copula representation, that is

S F (t, u) = C(S T (t), S U (u)),
where S T (t) = P(T > t) and S U (u) = P(U > u), and where C is a copula function (see e.g. [START_REF] Nelsen | An introduction to copulas[END_REF]). To understand the dependence between T and U, which is represented by the copula function C, it is natural to search for an estimator C, usually based on a parametric or semiparametric model. Nevertheless, even in this framework, nonparametric inference is required, at least to validate the model.

In section 3.2.1, we show how our estimator can be used to compute a n 1/2 -consistent estimator of S F . Section 3.2.2 illustrates its behaviour in view of estimating Kendall's τ coecient. This nonparametric estimation of τ can be used, in a copula modelling approach, to obtain preliminary estimators that can be used in the maximum likelihood estimation procedure. Finally, section 3.2.3 provides an uniform consistency result for an extension of Kendall's process, that is classically used to assess the validity of Archimedean copula models in absence of right-censoring and left-truncation.

Estimation of the joint distribution function

As a particular example of the use of Theorem 3.2, we provide an uniform consistency result for the estimation of the joint distribution function (T, U ). The survival function S F can be obtained in the same way, since both function are related. Indeed, we can apply Theorem 3.2 to the class of functions F = {(T, U ) :→ 1 T ≤t,U ≤u }. This class is a Donsker class of functions (see [START_REF] Van Der Vaart | Weak convergence and empirical processes with applications to statistics[END_REF], for a complete denition), and therefore satises condition (3.3). Moreover,

F (t, u) = ∫ ∫ 1 t ′ ≤t,u ′ ≤u d F (t ′ , u ′ ), which leads to sup t,u | F (t, u) -F (t, u)| = O P (n -1/2 ). (3.5)
The joint survival function can be studied in the same way by considering 1 T >t,U >u . We see that this convergence rate is of the same order as the one of Dabrowska's estimator [START_REF] Dabrowska | Kaplan-Meier estimate on the plane[END_REF]).

Estimation of Kendall's τ coecient

Kendall's τ coecient can be dened as

τ = P [(T 1 -T 2 )(U 1 -U 2 ) > 0] -P [(T 1 -T 2 )(U 1 -U 2 ) < 0] .
This coecient is a popular measure of dependence between to random variables, see e.g. [START_REF] Nelsen | An introduction to copulas[END_REF] for more details. After some standard calculations, Kendall's τ coecient can be rewritten as

τ = 4 ∫ ∫ F (t, u)dF (t, u) -1.
Therefore, our estimator F can be used to dene a nonparametric estimator of τ in the following way,

τ = 4 ∫ ∫ F (t, u)d F (t, u) -1. (3.6)
Convergence of τ towards τ can be shown combining (3.5) and Theorem 3.2. Indeed,

τ -τ = ∫ ∫ ( F (t, u) -F (t, u) ) dF (t, u) + ∫ ∫ F (t, u)d { F (t, u) -F (t, u) } + ∫ ∫ (F (t, u) -F (t, u))d { F (t, u) -F (t, u) } .
The rst and third term are O P (n -1/2 ) using (3.5), while the second is also O P (n -1/2 ) from Theorem 3.2.

In the following, we will use this estimator τ to dene preliminary estimators for the association parameter of survival copulas models under bivariate censoring and left truncation. Indeed, for some copula families, there exists a closed formula linking the parameter of association to Kendall's τ coecient. Otherwise one can use numerical techniques to determine, for example through simulations, which value of the association parameter corresponds to τ . The preliminary estimators obtained by this method can then be used as starting point for maximum likelihoood optimization. Indeed, these maximization techniques, especially in our incomplete data framework, rely on delicate numerical techniques that can be very sensitive to the initial point of the iterative maximization algorithm that is used. Since our estimator (3.6) is n 1/2 -consistent, it is expected to furnish a good initialization for such procedures.

Illustration.

To investigate the practical behaviour of our procedure, we consider a simulation study in order to assess the convergence of our estimator τ . We consider a couple of lifetimes (T, U ) with the same Weibull distribution for T and U. Weibull distribution is parametrized through a shape parameter α and a scale parameter β, and admits a density

f (t) = α β ( t β ) α-1 exp ( - t α β α
) ,

for t ≥ 0. The two random variables are linked by a Clayton copula (see denition in

Table 2) with parameter θ, so that Kendall's τ coecient is known, since

τ = θ[θ + 2] -1
for a Clayton copula. The censoring variables are assumed to be i. Table 1: Estimation of the mean-squared error and related quantities for the estimation of Kendall's τ coecient.

Kendall's process

Kendall's process is a standard tool to perform model checking for copulas in absence of censoring or truncation, see [START_REF] Genest | Statistical inference procedures for bivariate Archimedean copulas[END_REF], or [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF]. Applications will be shown in section 4.3. Dening K(v) = P(S F (T, U ) ≤ v), empirical versions of K can be computed in absence of censoring using the empirical distribution function.

In our framework, we can extend this denition, dening

K(v) = ∫ ∫ 1 S F (t,u)≤v d F (t, v) = n ∑ i=1 W i,n 1 S F (Y i ,Z i )≤v . (3.7)
The proof of the next result is technical and is postponed to the appendix section.

Proposition 3.3 Under the assumptions of Theorem 3.1, assume that the density of (T, U ) is uniformly bounded by some constant C ≥ 0. Then,

sup v∈[0,1] | K(v) -K(v)| = O P (n -1/2 ).
This theoretical result is a key tool to assess the validity of testing procedures for copula models. More precisions on this are given in section 4.3. [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] proposed a similar estimator based on Dabrowska's estimator, without proving its convergence.

Bootstrap procedure

Theorem 2.2 only provides an uniform convergence rate, without giving an asymptotic law, which could be used to compute condence intervals for estimations of integrals ∫ ∫ ϕ(t, u)dF (t, u). Nevertheless, the most appropriate way to estimate the asymptotic variance of such quantities is to rely on bootstrap approximations (even if we had determined the asymptotic law, bootstrap procedures would be more adapted for nite sample size).

Suppose that we wish to estimate the law of a quantity θ( F ), the particular example we have in mind being V ar

( ∫ ∫ ϕ(t, u)d F (t, u)
)

, but more general quantities may be studied with the similar method. For censored variables in the univariate, there exists two principal strategies for dening resampling plans. [START_REF] Reid | Estimating the median survival time[END_REF]) methodology consists of directly resampling the points (Y i , δ i ) 1≤i≤n , while Efron's methodology [START_REF] Efron | Censored data and the bootstrap[END_REF]) consists of simulating lifetimes and censoring variables to compute (Y i , δ i ). Although Reid's methodology is easier to implement, it can lead to unconsistant results, which is not the case of Efron's methodology, see [START_REF] Akritas | Bootstrapping the Kaplan-Meier estimator[END_REF]. Therefore, we try to extend this second approach to the bivariate case (in the univariate case, Efron's resampling plan has already been extended to the case of left truncation by [START_REF] Wang | Nonparametric estimation from cross-sectional survival data[END_REF]).

In order to simulate censoring variables (C, D), it is required to use an estimator of their joint distribution. Dening G(c, d) = P (C ≤ c, D ≤ d) , we can, reversing the role of the censoring and lifetime variable, compute an estimator Ĝ using the same method that we used for estimating the distribution of the lifetimes. We also need to simulate truncation variables (µ, ν). For this, we use

L(t, u) = ( n ∑ i=1 ŜF (Y i , Z i ) -1 ) -1 ( n ∑ i=1 ŜF (Y i , Z i ) -1 1 µ i ≤t,ν i ≤u ),
as an estimator of their joint distribution.

We here provide a resampling plan, inspired by the methodology of [START_REF] Efron | Censored data and the bootstrap[END_REF] for censored observations in the univariate case.

Resampling plan.

Let B denote the number of bootstrap sample that we whish to generate. For b = 1, ..., B,

1. Simulate (T b i , U b i ) according to F . 2. Simulate (C b i , D b i ) according to Ĝ. 3. Compute (Y b i , Z b i , δ b i , γ b i ). 4. Simulate (τ b i , υ b i ) according to L. 5. If Y b i ≥ τ b i and Z b i ≥ υ b i , keep observation i. Otherwise, resimulate from
Step 1 until this is the case.

6. Obtain an n-sample reiterating this procedure independently, and compute Fb and θ( Fb ).

The vector (θ( Fb )) 1≤b≤B can then be used to estimate dierent characteristics, such as the variance, of θ( F ).

Remark 3.4 Simulating random vector with respect to discrete distributions such as F , Ĝ and L is easy. We here see one of the advantages of having used an estimator F which is linked to a true distribution. If we use estimators such as the one provided by Dabrowska (1984), this resampling plan can not be adapted since Step 1 and 2 are not feasible. Nevertheless, there is a specicity of F and Ĝ that we must take into account: in some cases, these two distributions will present some mass at (∞, ∞). In these cases, some observations which are resimulated should have the value (∞, ∞). The probability occurence of such events should be small (both distributions put a small mass at innity). Such observations will only have the consequence to add some weight at (∞, ∞) for Fb , without leading to failures. However, an alternative would be to remove these observations from the simulated sample, which would not strongly modify the results.

Application to the Canadian data set

In this section, we apply our technique to a portfolio from a Canadian insurer, which has been studied in several papers. In this paper, the authors provided a very detailed study of the joint mortality using Archimedean copula models. The choice of the copula family is made using a technique based on a version of Kendall's process proposed by [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] and based on Dabrowska's estimator. We here compare our results with the new estimator that we propose in order to assess its validity. Moreover, in section 4.3, we provide a test procedure to perform model-checking, and explain how to compute the critical values.

4.1

Description of the data-set

The dataset that we study has been initially introduced by [START_REF] Frees | Annuity valuation with dependent mortality[END_REF], and was studied by [START_REF] Carriere | Bivariate survival models for coupled lives[END_REF], [START_REF] Youn | Statistical aspects of joint life insurance pricing[END_REF], [START_REF] Youn | Pricing practices for joint last survivor insurance[END_REF] and [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF]. We refer to [START_REF] Frees | Annuity valuation with dependent mortality[END_REF] for a more detailed description of this dataset. The dataset concerns 14947 contracts from a large Canadian insurer, observed between Decembre 29th, 1988 andDecembre 31th, 1993. After elimination of same-sex contracts and of couples with more than one policies (for which we only keep one policy), 11454 contracts remain. To model the data, we consider three copula models proposed by [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] to study a particular generation (here, we apply our method to the whole population), and show how our technique allows to perform estimation and goodness-of-t. where ŜT (resp. ŜU ) denotes the Kaplan-Meier estimator of S T (t) = P(T > t) (resp.

S U (u) = P(U > u)), or any other consistent estimator of the survival function. We perform maximization of this quantity using a Newton-Raphson algorithm, with initial value taken as the preliminary estimator of Table 2. The nal estimators are presented in Table 3, along with the corresponding value of Kendall's τ coecient.

4.3

A class of model-checking procedures using Kendall's process

Test-statistics

Archimedean families of copulas present the advantage of being uniquely dened by their generator ϕ, which can be linked to function K through the relation,

K(v) = v - ϕ(v) ϕ ′ (v) .
Genest and Rivest (1993) and [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF] proposed test procedures (in absence or censoring) that can be used for model-checking based on the comparison of

K θ(v) = v-ϕ θ(v)/ϕ ′ θ(v)
, where θ denotes the estimator of the association parameter in a parametric copula model, and a nonparametric version K.

Let C denote the true copula function linking the marginal distributions of T and U.

If one whishes to test H

0 : C ∈ C Θ against H 1 : C / ∈ C Θ ,
where C Θ is a parametric family of copulas, the procedure can be summarized as follows:

1. Estimate θ by θ assuming that model C Θ is true.

2. Compute K, K θ, and d( K, K θ) where d denotes some distance.

3. If d( K, K θ) is larger than some critical value s α , reject H 0 .

Luciano, Vigna and Spreeuw (2008) use a L 2 -distance, that is

d(K, K ′ ) = ∫ 1 0 {K(v) -K ′ (v)} 2 dv, (4.1)
and we will use the same distance for easier comparison. Next, we explain in section 4.3.2 how to compute the critical values.

Based on Proposition 3.3, one can see that this testing approach is consistent. Indeed, if H 0 is not true, our estimator K converges towards some limit K such that

inf θ∈Θ d(K, K θ ) > 0. Hence, since d( K, K θ) > |d(K, K θ) -d( K, K)|, we get d( K, K θ) > inf θ∈Θ d(K, K θ ) + O P (n -1/2 ).
As a consequence, for any α, P

( d( K, K θ) > s α |H 1 ) → n→∞ 1.

Computation of the critical values

We propose to compute the critical values using a procedure which is adapted from [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF] in the uncensored case. The critical values are computed using a bootstrap resampling plan. Let us mention that [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF], in our framework, proposed a very simple resampling plan. However, [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF] showed that this resampling plan was unconsistent and therefore could not be used. Therefore, we propose to adapt our bootstrap procedure dened in section 3.3 in order to extend the procedure of [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF]. We apply this technique for θ( F ) = d( K, K θ), where d is the distance dened in (4.1). However, since we are interesting in knowing the distribution of d( K, K θ) under the null hypothesis, we simply change Step 1 in our bootstrap procedure of section 3.3 by 1'. Simulate (T, U ) according to the distribution C θ( FT (t), ĜU (u)).

The bootstrap vector (d( Kb , K θb )) that we obtain using this procedure can then be used to determine approximated p-values of our test procedure.

Another modication of the bootstrap procedure is required to reduce computation time. Indeed, this simulation scheme suers from the diculty to estimate the distribution in the left-tail (caused by truncation). Even if condition (3.1) holds, the constant K may be large and impact on the performance for nite sample-size (this is actually the case in the univariate situation). In practice, we see that the nonparametric estimator of the distribution of (C, D) is very sensitive to this phenomenon. As an illustration, the smallest censored observations (for both components) has mass 0.5. Moreover, the set of observations such that one at least is less than 20000 days is aected of a mass greater than 0.8, although it represents only 2% of the observations. This is not a problem from a theoretical point of view: if we draw some point (C, D) in this part of the distribution, the truncation phenomenon will frequently reject this observation, and therefore it will not appear in the nal simulated sample. However, in practice, according to our procedure, each time a point is rejected due to truncation, another simulation is performed. If these rejection happen very often, this dramatically increases the computation time. Therefore, we chose to use modify the distribution. We aect to the set of censored points such that C < 20000 or D < 20000 a mass 0.2, and renormalize the rest of the mass in order to obtain a total mass equal to 1.

Goodness-of-t on the Canadian data-set

In Figure 1, we graphically present the function K computed from the Canadian data-set, and compare it with the function Kθ obtained for the three copula models. One can observe that the estimated version using Clayton's model or Nelsen's 4.2.20 model seem to be closer to the nonparametric curve. Nevertheless, all three models seem to be less accurate for high values of v, which may be caused by the presence of censoring.

The application of the test procedure described in the previous section is presented in Table 4. In this table, we compare the value of the distance obtained by our method for the dierent models that we considered, to the one obtained by [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] using the version of K proposed by [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] According to these results, Frank's model corresponds to the smallest p-value, and seems less adapted than the two other models, which corresponds to what was observed on Figure 1. At a classical 5% level, one can see that, for all three models, the null hypothesis is not rejected.

Conclusion

In this paper, we proposed a new nonparametric estimator of the joint distribution of a random vector of lifetimes (T, U ). This estimator is designed to compensate the classical problem that impact lifetime data, such as left truncation and right-censoring. Unlike previously existing estimators in this framework, our estimator denes a true distribution, and does not only focus on the particular case of estimating the joint survival function.

An important advantage of this approach, is that it allows to perform simulations quite easily, as it has been proposed in the bootstrap procedures that we discussed in section 3.3 and 4.3.2. These simulations procedure may be quite valuable in the computation of critical values of goodness-of-t tests, as it is illustrated on the Canadian data-set of [START_REF] Frees | Annuity valuation with dependent mortality[END_REF].

6 Appendix 6.1 Proof of Theorem 3.1

In this proof, we introduce several notations in order to make some simplications. We will denote H i = H(Y i , Z i ), ∆ i = δ i γ i , S i = S F (Y i , Z i ), and

S * i = 1 -F * (Y i , Z i ) = n -1 ∑ n j=1 ∆ j H -1 j 1 Y j >Y i ,Z j >Z i .
The vector W * will represent the weights aected by F * to each observation. Only the doubly uncensored observations (which are the one where some mass is put) are present in this vector of size m + 1 (as in the proof of Theorem 2.2, m will denote the number of doubly uncensored observations) This vector is sorted accordingly to Y (the rst component corresponds to the mass assigned at the doubly uncensored observation with the smallest Y, and so on). W * i will denote the i-th component of W * . We recall some facts that can be easily proven using empirical process theory, see e.g. [START_REF] Van Der Vaart | Weak convergence and empirical processes with applications to statistics[END_REF] for more details:

1 n n ∑ i=1 ∆ i S i H i ,
that is the residual mass that the ideal estimator F * puts at innity. W * ∞,∞ is a sum of i.i.d. quantities, with E[W * ∞,∞ ] = 0 (from Proposition 2.1), and with nite variance. Then, the CLT applies, leading to W * ∞,∞ = O P (n -1/2 ). Next, as in the proof of Theorem 2.2, we remove the components of W which are supposed to be zero, considering the equation (I m+1 -Ã) W = 0, where m denotes the number of doubly uncensored responses. We will also assume that the components of W are sorted with respect to the variable Y. Dene W ′ as a solution of (I m+1 -Ã)W ′ = 0 with

W ′ ∞ = W * ∞ .
Note that the vector W ′ does not necessarily dene a probability measure, since the some of its components is not equal to one. However, W ′ is proportional to W, and denes a positive measure. In this proof, we rst show that W ′ is close to W * , and then show that W is close to W ′ .

We denote

S j = ∑ i≥j |W ′ i -W * i | + |W ′ ∞,∞ -W * ∞,∞ |, with S m+1 = |W ′ ∞,∞ -W * ∞,∞ | = 0.
Let R i be dened as in the proof of Theorem 2.2. We have, for j ≤ m, 0 ≤ S j ≤ S j+1 + |W ′ j -W * j |. Next, observe that

W ′ j = 1 n ∑ i>j W ′ i 1 i∈R j Ĥj = 1 n ∑ i>j W ′ i 1 i∈R j H j + 1 n ∑ i>j W ′ i 1 i∈R j [H j -Ĥj ] Ĥj H j , (6.3) W * j = 1 n S j H j = 1 n ∑ i>j W * i 1 i∈R j H j + 1 n S j -S * j H j .
(6.4) Therefore,

S j ≤ ( 1 + K n ) S j+1 + K 2 n sup j H j | Ĥj -H j | Ĥj + K n sup j |S * j -S j |.
Since S m+1 = 0, one gets

S j ≤ ε n (1 + Kn -1 ) m-j -1 Kn -1 ≤ C(K)ε n ,
where

ε n = K 2 n sup j H j | Ĥj -H j | Ĥj + K n sup j |S * j -S j |,
and where C(K) is a positive constant which does not depend from n or j. Moreover, on

A n we have

sup j H j | Ĥj -H j | Ĥj ≤ sup t,u | Ĥ(t, u) -H(t, u)| × sup t,u
H(t, u) Ĥ(t, u) = O P (n -1/2 ), from (6.1). Indeed, the rst supremum in the right-hand side is exactly the supremum considered in (6.1). To see that sup |S * (t, u) -S(t, u)| = O P (n -1/2 ), from (6.2). This leads to sup j S j = O P (n -1/2 ).

Going back to (6.3) and (6.4), we see that

|W ′ j -W * j | ≤ 1 n sup j S j + C(K) n ε n := η n n ,
with η n = O P (n -1/2 ). From this, we deduce that ∑ n j=1 W ′ j =

∑ n j=1 W * j + r n = 1 + r n , with r n = O P (n -1/2 ). Next, we use the fact that

W j = W ′ j [ ∑ n j=1 W ′ j ] -1 to deduce that |W j -W * j | ≤ |W ′ j -W * j | 1 + r n + W * j r n = θ n n ,
with θ n = O P (n -1/2 ), since W * j ≤ Kn -1 .

  ) the ordered values of uncensored variables Y, and Z (1) < . . . < Z (n) the corresponding values of uncensored variables Z (ordered with respect to Y ). With a slight abuse in notations, we do not change the notation à after reordering. Next, one can observe that

  i.d. exponential with scale parameter equal to one. The truncation variables (µ, ν) are assumed to be i.i.d. with marginal distribution dened as the product B × W of a Bernoulli random variable B with mean value 0.7, and a Weibull random variable W of parameters α = 2 and β = 10 which is independent from B.We considered two dierent situations for the values of the parameters of the marginal laws of T and U. In the rst situation, we consider α = 2 and β = 1.1, which corresponds approximatively to a proportion of 20% of observations which are uncensored for both marginals. Additionally, we considered the case where α = 2 and β = 1.7, corresponding to 10% of observations which are uncensored for both marginals. In each case, we took θ = 2, corresponding to τ = 0.5, and we considered two sample sizes n = 1000, and n = 2000. Through 1000 replications of each setting, we computed an estimation of the mean squared error E [(τ -τ ) 2 ] , and decomposed it into a bias term and a variance term.

Figure 1 :

 1 Figure 1: Comparison of the estimators Kθ to the nonparametric estimator K.

  sup t,z | Ĥ(t, z) -H(t, z)| = O P (n -1/2 ), (6.1) sup t,z |F * (t, z) -F (t, z)| = O P (n -1/2 ).

  u) = O P (1), we can use the uniform consistency of Ĥ coming from (6.1) to see that, for n large enough, on a set with probability tending to one, | Ĥ(t, u) -H(t, u)| ≤ 1/2K, and thereforeH(t, u)/ Ĥ(t, u) ≤ 1/[1 -1/2]. Next, recall that sup j |S * j -S j | ≤ sup t,u

Table 3 :

 3 Final maximum likelihood estimators of the association parameter in the three copula models, with corresponding value of Kendall's τ coecient.

	Model	Final estimator Corresponding value of Kendall's τ
	Clayton	1.3116	0.3960
	Frank	2.91290	0.3012
	Nelsen 4.2.20	0.5098	0.4020

Table 4 :

 4 Goodness-of-t procedure for a the three survival copula models considered(Clayton, Frank, Nelsen 4.2.20).

	instead of our estimator	K

4.2

Copula models and estimation of the association parameter In this section, we consider a semiparametric copula model for describing the distribution of (T, U ). In a semiparametric copula model, the marginal distribution of each variable is estimated nonparametrically (using Kaplan-Meier estimator since we are in presence of censoring and truncation), while the dependence between the two variables is described by a parametric copula family {C θ : θ ∈ Θ} for Θ ⊂ R k , see e.g. [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF].

We consider the three families of copulas that have been used by [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] for studying the joint mortality of a specic generation in the Canadian data-set, that is Clayton, Frank and Nelsen 4.2.20. These three families are of Archimedean type, which means that the copula function is uniquely dened by a decreasing convex function ϕ,

The generators corresponding to the three families are presented in Table 2 below, and we recall the link between τ and the association parameter in these families.

Using our estimator (3.6), we obtain τ = 0.4332. Using the relationship between the association parameter θ and τ, we obtain the following preliminary estimators of θ,

Table 2: Expression of the dierent copula families considered. The column τ presents the expression of Kendall's τ if available, the column preliminary represents the preliminary estimator obtained by calibration with τ coecient.

These preliminary estimators are then used to perform maximum likelihood estimation, using a procedure which has been studied by [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF], [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF], and which is also used in [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF]. In this bivariate censoring framework, the log-likelihood can be decomposed into four terms, that is

)

and the denition of the O P (n -1/2 ) convergence rate, lim x→∞ lim n P(A c n,ε ) = 0. On this set,

We see that, on

3) (it is a Donsker class of bounded functions, see [START_REF] Van Der Vaart | Weak convergence and empirical processes with applications to statistics[END_REF]).

Hence, Theorem 3.2 applies. Therefore, on A n,ε

). Since we assumed that the density of (T, U ) is bounded by some constant C, the probability in the right-hand side of (6.5) is less than Cεn -1/2 .

We can deduce that

Using (6.5), we see that, for x large enough, the rst probability on the right-hand side of (6.6) is less than P(sup ε,v n 1/2 |R n (ε, v)| > x). By denition of the O P (n -1/2 ) rate,