Elena Di Bernardino
email: elena.di_bernardino@cnam.fr

Didier Rullière
email: didier.rulliere@univ-lyon1.fr

Estimation of multivariate critical layers: Applications to rainfall data

Keywords:

Calculating return periods and critical layers (i.e. multivariate quantile curves) in a multivariate environment is a difficult problem. A possible consistent theoretical framework for the calculation of the return period, in a multi-dimensional environment, is essentially based on the notion of copula and level sets of the multivariate probability distribution. In this paper we propose a fast and parametric methodology to estimate the multivariate critical layers of a distribution and its associated return periods. The model is based on transformations of the marginal distributions and transformations of the dependence structure within the class of Archimedean copulas. The model has a tunable number of parameters, and we show that it is possible to get a competitive estimation without any global optimum research. We also get parametric expressions for the critical layers and return periods. The methodology is illustrated on rainfall 5-dimensional real data. On this real data-set we obtain a good quality of estimation and we compare the obtained results with some classical parametric competitors. Finally we provide a simulation study.

Introduction

Return Periods

The notion of Return Period (RP) is frequently used in environmental sciences for the identification of dangerous events, and provides a means for rational decision making and risk assessment. Roughly speaking, the RP can be considered as an analogue of the "Value-at-Risk" in Economics and Finance, since it is used to quantify and assess the risk (see, e.g., [START_REF] Nappo | Kendall distributions and level sets in bivariate exchangeable survival models[END_REF]). In engineering practice, finance, insurance and environmental science the choice of the RP depends on the impact/magnitude of the considered event and the consequences of its realisation.

Equally important is the related concept of design quantile, usually defined as "the value of the variable characterizing the event associated with a given RP". In the univariate case the design quantile is usually identified without ambiguity. Conversely in the multivariate setting different definitions are possible (see [START_REF] Serfling | Quantile functions for multivariate analysis: approaches and applications[END_REF]). For this reason, the identification problem of design events in a multivariate context has recently attracted the attention of many researches. The interested reader is referred for example to [START_REF] Embrechts | Bounds for functions of multivariate risks[END_REF], [START_REF] Belzunce | Quantile curves and dependence structure for bivariate distributions[END_REF], [START_REF] Nappo | Kendall distributions and level sets in bivariate exchangeable survival models[END_REF], in the economics and finance context; to [START_REF] Chebana | Index flood-based multivariate regional frequency analysis[END_REF], [START_REF] Chebana | Multivariate quantiles in hydrological frequency analysis[END_REF] (and references therein) in the hydrological context.

During the last years, researchers in environmental fields joined efforts to properly answer the following crucial question: "How is it possible to calculate the critical design event(s) in the multivariate case?" (see for instance [START_REF] Salvadori | Extremes in Nature: An Approach Using Copulas[END_REF]). In this sense, a possible consistent theoretical framework for the calculation of the design event(s) and the associated return period(s) in a multi-dimensional environment, is proposed, e.g., by [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF], [START_REF] Salvadori | Semi-parametric approximation of Kendall's distribution function and multivariate Return Periods[END_REF], [START_REF] Gräler | Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation[END_REF]. In particular the authors define the multivariate return period using the notion of upper and lower level sets of multivariate probability distribution F and of the associated Kendall's measure.

In the following, we will consider a sequence X = {X 1 , X 2 , . . .} of independent and identically distributed d-dimensional random vectors, with d > 1. Thus each X k , k ∈ N, has the same multivariate distribution F X : R d + → [0, 1] as the nonnegative real-valued random vector X ∼ F X = C(F X1 , . . . , F X d ) describing the hydrological phenomenon under investigation. The function C is the d-dimensional copula associated to F (see [START_REF] Nelsen | An introduction to copulas[END_REF]). We write I = {1, . . . , d} the set of indexes of the considered random variables and of their associated cumulative distribution functions, i.e., F Xi (x i ) = P (X i ≤ x i ), for i ∈ I.

In the following we will consider multivariate distribution functions F X satisfying these regularity conditions for all u ∈ [0, 1], the diagonal of the copula C, i.e. C(u, . . . , u), is a strictly increasing function of u;

for any i ∈ I, the marginal F Xi is continuous and strictly monotonic distribution function.

In this setting, we introduce the notion of critical layer (see, e.g., [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF], [START_REF] Salvadori | Semi-parametric approximation of Kendall's distribution function and multivariate Return Periods[END_REF], [START_REF] Gräler | Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation[END_REF]).

Definition 1.1 (Critical layer). The critical layer ∂L(α) associated to the multivariate distribution function F X of level α ∈ (0, 1) is defined as

∂L(α) = {x ∈ R d : F X (x) = α}.
Then ∂L(α) is the iso-hyper-surface (with dimension d -1) where F equals the constant value α. Thus, ∂L(α) is a (iso)line for bivariate distributions, a (iso)surface for trivariate ones, and so on. The critical layer ∂L(α) partitions R d into three non-overlapping and exhaustive regions:

   L < (α) = {x ∈ R d : F X (x) < α}, ∂L(α) = the critical layer itself, L > (α) = {x ∈ R d : F X (x) > α}.
Practically, at any occurrence of the hydrological considered phenomenon, only three mutually exclusive events may happen: either a realization of the considered hydrological event lies in one of these 3 Borel sets L < (α), ∂L(α), or L > (α).

In the applications, usually, the event of interest is of the type {X ∈ A}, where A is a non-empty Borel set in R d collecting all the values judged to be "dangerous" according to some suitable criterion. A natural choice for A is the set L > (α) (see [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF], [START_REF] Gräler | Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation[END_REF]). The first random index N where X k reaches the set L > (α) is N = min k∈N {X k ∈ L > (α)}. Assuming P [X ∈ L > (α)] ∈ (0, 1), one easily shows that N is a geometric random variable. The Return Period is defined as the average time required for reaching the set L > (α), that is:

RP > (α) = ∆ t • E [N ] = ∆ t P [X ∈ L > (α)] , (1)
where ∆ t > 0 is the (deterministic) average time elapsing between X k and X k+1 , k ∈ N. The probability that a realization of this vector belongs to L < (α) is given by the Kendall's function, which only depends on the copula C of this random vector, i.e., K C (α) = P X ∈ L < (α) = P [C(U 1 , . . . , U d ) ≤ α] , for α ∈ (0, 1).

(2)

Then, the considered Return Period can be expressed using Kendall's function in (2), RP > (α) = ∆ t • 1 1-K C (α) . Obviously, Return Periods can naturally be associated to other sets than L > (α), the interested reader is referred for example to [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF].

This paper aims at:

• giving a parametric representation of the multivariate distribution F of a random vector X, here representing rain measurements (for applications see Section 6),

• giving direct estimation procedure for this representation,

• giving closed parametric expressions, both for critical layers in Definition 1.1 and Return Periods in (1),

• adapting this methodology to some asymmetric dependencies (as, for instance, non-exchangeable random vectors; for a possible investigation in this sense see Section 7).

In the next section, we introduce the model used to answer the issues introduced above.

The model

We consider the following model, which is detailed in Di Bernardino and Rullière (2013a),

F (x 1 , . . . , x d ) = T • C 0 (T -1 1 • F 1 (x 1 ), . . . , T -1 d • F d (x d )), (3)
or equivalently      F (x 1 , . . . , x d ) = C( F 1 (x 1 ), . . . , F d (x d )), with C(u 1 , . . . , u d ) = T • C 0 (T -1 (u 1 ), . . . , T -1 (u 1 ))

F i (x) = T • T -1 i • F i (x), for i ∈ I, (4)
where F 1 , . . . , F d are given parametric initial marginal cumulative distribution functions, and where C 0 is a given initial copula. Hence the distribution F (x 1 , . . . , x d ) is built from transformed marginals F i , i ∈ I and from a transformed copula C, under regularity conditions. Transformation T permits to transform the initial dependence structure C 0 . For a given T , transformations T i permit to transform marginals, i ∈ I. All these transformations are described hereafter.

As we will see in the following, the initial copula C 0 in (3) is not estimated but it is chosen at the beginning of the estimation procedure. So it can represent some kind of a priori belief on dependence structure of the data or on the considered problem, that will be transformed in order to improve the fit.

Furthermore, as in Di Bernardino and Rullière (2013b), we will assume in the following that C 0 is an Archimedean copula. This means that in this paper, we mainly consider copulas that can be written as

C φ (u 1 , . . . , u d ) = φ(φ -1 (u 1 ) + . . . + φ -1 (u d )),
where the function φ is called the generator of the Archimedean copula C φ . Some conditions like d-monotony are given in [START_REF] Mcneil | Multivariate archimedean copulas, d-monotone functions and l 1 -norm symmetric distributions[END_REF]. Here we choose strict generator, i.e., φ(t) > 0, ∀ t ≥ 0 and lim

t→+∞ φ(t) = 0, with proper inverse φ -1 such that φ • φ -1 (t) = t.
The function T : [0, 1] → [0, 1] is a continuous and increasing function on the interval [0, 1], with T (0) = 0, T (1) = 1, with supplementary assumptions that will be chosen to guarantee that C is also a copula (detailed hereafter). Internal transformations T i : [0, 1] → [0, 1] are continuous non-decreasing functions, such that T i (0) = 0, T i (1) = 1, for i ∈ I. Conditions on transformations such that C is a copula are discussed for example in [START_REF] Durante | Distorted copulas: Constructions and Tail Dependence[END_REF], Di Bernardino and Rullière (2013a), Di Bernardino and Rullière (2013b).

Remark that among problems generated by transformations of Archimedean copulas, one can point out in particular the problem of uniqueness. Transformations of a given initial copula leading to a given target copula are not unique. This raises some problems for the analysis of the convergence of estimators of the transformation. This also causes problems to compare transformations and to understand their impact on the dependence structure. A further analysis shows that also a generator of an Archimedean copula is not unique, causing the same kind of problems. Then in Di Bernardino and Rullière (2013b), the definition of equivalence classes for both transformations and generators is provided to select some standardized forms for practical use, for the comparison and the interpretation of obtained distribution functions. Firstly equivalent classes for transformations can be characterized. Furthermore one can ensures the uniqueness of the transformation T by passing through the point (x 0 , y 0 ), among the invariant class for transformations (see Lemma 2.2 and Corollary 2.1 of the aforementioned paper). Obviously, in an iterative procedure of estimation, the uniqueness of the transformation T is essential in order to permit the convergence of the procedure and the identifiability of the considered transformation model.

In the following we detail the proposed semi-parametric estimation procedure in order to easily fit the model in (3). We show in particular how to estimate the transformations T and T i , i ∈ I.

Organization of the paper

The paper is organized as follows. In Section 2 we developed the estimation procedure for the chosen model in (3). In Section 2.1 we focus on a central tool in our estimation procedure: the estimation of the diagonal section of a copula. This estimated diagonal section is used in the non-parametric estimation of the external transformation T and the internal transformations T i (see Section 2.2). The parametric estimation using composite hyperbolic transformations is detailed in Section 3. Then in Section 4 we propose the parametric form for the desired quantities: the multivariate distribution function F , its critical layers ∂L(α) and the associated Kendall's function K C (α), for α ∈ (0, 1). Finally, Section 6 is devoted to a detailed study of a 5-dimensional rainfall data-set. Furthermore, a nested model is proposed for this 5-dimensional data-set, using the whole estimation procedure presented in the paper (see Section 7). Finally we provide a simulation study in Section 8. Directions for future research and the conclusion are in Section 9.

2 Nonparametric estimation

Initial diagonal section estimation

In order to give a non-parametric estimation of the external transformation T , we will propose a non-parametric estimator of the copula within the Archimedean class of copulas. Several non-parametric estimators of the generator of an Archimedean copula are available. One can cite for example those of [START_REF] Dimitrova | GeD spline estimation of multivariate Archimedean copulas[END_REF] or [START_REF] Genest | Inference in multivariate Archimedean copula models[END_REF] both based on empirical Kendall's function. Here, proposed estimations will rely on an initial estimation of the diagonal section of the empirical copula. Firstly, we propose some (classical) estimators of the diagonal section δ 1 of a copula,

δ 1 (u) = C(u, . . . , u) , u ∈ [0, 1].
Remark that if (U 1 , U 2 , . . . , U d ) has as distribution function the copula C then

P[U 1 ≤ u, . . . , U d ≤ u] = C(u, . . . , u) = P[max{U 1 , U 2 , . . . , U d } ≤ u] = δ 1 (u).
As a consequence, estimators based on the diagonal section of a copula are relying on the distribution of the maximum of U 1 , . . . , U d , as explained in [START_REF] Sungur | Diagonal copulas of Archimedean class[END_REF]. As pointed out in Di Bernardino and Rullière (2013b), the diagonal of an Archimedean copula, under some suitable conditions, is essential to describe the copula. So, in the following we recall important assumptions (which are fulfilled for many Archimedean copulas, including the independent copula) for the unique determination of an Archimedean copulas starting from its diagonal section (see, for instance, [START_REF] Erdely | Frank's condition for multivariate Archimedean copulas[END_REF] and references therein). Some constructions of copulas starting from the diagonal section are given for example in [START_REF] Nelsen | On the construction of copulas and quasi-copulas with given diagonal sections[END_REF] and [START_REF] Wysocki | Constructing Archimedean copulas from diagonal sections[END_REF].

Proposition 2.1 (Identity of Archimedean copulas, Theorem 3.5 by [START_REF] Erdely | Frank's condition for multivariate Archimedean copulas[END_REF]). Let C a d-dimensional Archimedean copula whose diagonal section δ 1 satisfies δ 1 (1 -) = d. Then C is uniquely determined by its diagonal.

Condition in Proposition 2.1 is referred to as Frank's condition in [START_REF] Erdely | Frank's condition for multivariate Archimedean copulas[END_REF] (see their Theorems 1.2 and 3.5). Note that if |φ (0)| < +∞ then the condition on the diagonal in Proposition 2.1 is automatically satisfied. In this case, the function φ can be reconstructed from the diagonal δ (see also [START_REF] Segers | Diagonal sections of bivariate Archimedean copulas. Discussion of "Inference in multivariate Archimedean copula models[END_REF]). As pointed out by [START_REF] Embrechts | Comments on: Inference in multivariate Archimedean copula models[END_REF] a possible limitation is that if φ has finite right-hand derivative at zero, the Archimedean copula generated by φ has upper tail independent structure. To show that the situation of many Archimedean copulas having the same diagonal is far from exceptional, a recipe to construct further examples is given in [START_REF] Segers | Diagonal sections of bivariate Archimedean copulas. Discussion of "Inference in multivariate Archimedean copula models[END_REF]. For further details the interested reader is referred to Section 3 in Di Bernardino and Rullière (2013b). It has been shown however that if one aims at fitting both lower and upper tail dependence, then some specific dependence structures can be suited to this purpose, see Di Bernardino and Rullière (2014). Here we will not focus on very extreme tail behavior, but on the main part of the distribution. Using the necessary few data in the tails would require some specific estimators of tail dependence coefficients, as detailed in above mentioned article.

In the following we present different estimation of the diagonal section. Some empirical copula estimators for δ 1 are given in [START_REF] Deheuvels | La fonction de dépendance empirique et ses propriétés[END_REF]. A comparison between several more recent estimators is presented in [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF].

Empirical diagonal

We present here an estimator that is detailed in [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF], and directly inspired by the one of [START_REF] Deheuvels | La fonction de dépendance empirique et ses propriétés[END_REF]. Consider observations

{X (k) = (X (k) 1 , . . . , X (k)
d )} k∈{1,...,n} of the random vector X. Define pseudo-observations

U (k) = (U (k) 1 , . . . , U (k) d ) by setting every component i for observation number k to U (k) i = 1 n + 1 n j=1 1 X (j) i ≤X (k) i
, with i ∈ I, k ∈ {1, . . . , n}. One can check that for any i ∈ I, k ∈ {1, . . . , n}, U (k) i ∈ (0, 1).

The empirical estimator C of the copula of vector X, is

C(u 1 , . . . , u d ) = 1 n n k=1 1 U (k) 1 ≤u1,...,U (k)
d ≤u d .
We thus obtain for any u ∈ [0, 1], the empirical estimation for the diagonal of C and its inverse, i.e.,

δ emp 1 (u) = C(u, . . . , u) , δ emp -1 (u) = arginf x ∈ [0, 1]; δ emp 1 (x) ≥ u .

Smooth empirical diagonal

We do not present here all possible smooth estimators of a copula and the associated smoothed estimates of the diagonal section. We will restrict ourselves in estimators based on transformations, in the same spirit as the other parts of this paper. These estimators perform reasonably well (see [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF]), especially considering Cramér-von Mises distance.

Using a smooth estimation of the empirical cumulative distribution of U could create some bias since the distribution support must be [0, 1], this is the reason why we consider a transformation of pseudo-observations. For given smoothing coefficients h 1 , . . . , h d , we can define

L (k) = (L (k) 1 , . . . , L (k)
d ) where

L (k) i = G -1 (U (k) i ),
with i ∈ I, k ∈ {1, . . . , n}, and where G -1 is the inverse of a cumulative distribution function, for example G -1 (x) = logit(x), and one can check that ((G) (x)) 2 /G(x) is bounded, which is a required condition detailed in [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF].

A smooth version of C can be defined by:

C(u 1 , . . . , u d ) = 1 n n k=1 i∈I K G -1 (u i ) -L (k) i h i ,
where K is a suited kernel function (we took here a multiplicative multivariate kernel).

We thus obtain for any u ∈ [0, 1], a smooth estimator δ of the diagonal section of C, and a smooth estimator δ -1 of the inverse function of this diagonal section:

δ 1 (u) = C(u, . . . , u) δ -1 (u) = arginf x ∈ [0, 1]; δ 1 (x) ≥ u . (5)
Note that initial pseudo-data can also be smoothed, by defining

U (k) i = 1 n + 1 j∈{1,...,n} Φ X (k) i -X (j) i h i .
Discussions on this last choices can be found in [START_REF] Chen | Nonparametric estimation of copula functions for dependence modelling[END_REF] and [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF].

A discussion on the possible choices for h 1 , . . . , h d is given in [START_REF] Chiu | A comparative review of bandwidth selection for kernel density estimation[END_REF] (for d = 1), [START_REF] Wand | Comparison of smoothing parameterizations in bivariate kernel density estimation[END_REF], [START_REF] Wand | Multivariate plug-in bandwidth selection[END_REF], [START_REF] Zhang | A Bayesian approach to bandwidth selection for multivariate kernel density estimation[END_REF] (for d ≥ 2) and references therein. However, according to [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF], "a good bandwidth selection rule is missing, for the moment, and is subject of further research."

A summary of the needed input parameter and the estimation of δ 1 and δ -1 1 is given in Algorithm 1.

Algorithm 1 Smooth initial diagonal estimation

Input parameters

Choose G -1 inverse of a cdf, e.g., G -1 (x) = logit(x) Choose h i , i ∈ I bandwidth sizes, e.g., the ones proposed by the Silverman's Rule of Thumb

Estimation

For any x ∈ [0, 1], get δ 1 (x) and δ -1 (x) by Equation ( 5)

2.2 Nonparametric estimation of transformations T and T i , i ∈ I

A non parametric estimator of the external transformation T in (3) is given in Di Bernardino and Rullière (2013b), for x ∈ (0, 1), by

T (x) = δ r(x) (y 0 ), (6)
with r(x) such that δ 0 r(x) (x 0 ) = x, where δ r refers to the estimator of the self-nested diagonal δ r of the target copula C, and δ 0 r(x) (x 0 ) refers to the self-nested diagonal of the initial copula C 0 . These estimators δ r and self-nested diagonals are defined hereafter.

In the case where the initial copula C 0 is an Archimedean copula of generator φ 0 , then

δ 0 r(x) (x 0 ) = φ 0 d r(x) φ -1 0 (x 0 ) and r(x) = 1 ln d ln φ -1 0 (x) φ -1 0 (x 0 )
.

In particular, if C 0 is the independence copula, with generator φ 0 (x) = exp(-x), and setting for example x 0 = y 0 = exp(-1), then T (x) = δ ln(-ln(x))/ ln d (e -1 ).

Let δ 1 be an estimator of δ 1 , and δ -1 be an estimator of the inverse function δ -1 as in Equation ( 5). At a relative integer order k ∈ Z, the self-nested diagonals estimators are defined as

     δ k (u) = δ 1 • . . . • δ 1 (u), (k times), k ∈ N δ -k (u) = δ -1 • . . . • δ -1 (u), (k times), k ∈ N δ 0 (u) = u. (7)
At any real order r ∈ R, an estimator δ r of the self-nested diagonal is

δ r (x) = z z -1 • δ k (x) 1-α z -1 • δ k+1 (x) α , x ∈ [0, 1]
with α = r -r and k = r , where r denotes the integer part of r, and where z is a function driving the interpolation, ideally (if known) the generator of the considered copula C. Some consistency results about this estimator T are detailed in Di Bernardino and Rullière (2013b).

For invertible marginal transformations T i , one easily shows

T i = F i • F -1 i • T, i ∈ I. (8)
Then, by replacing the transformed marginal distribution F i in (8) by an estimator of the target i-th marginal distribution, denoted by F i (e.g., empirical cdf), and given the external transformation or its consistent estimation T , a non-parametric estimation of T i , for i ∈ I, is provided by

T i (u) = F i • F -1 i • T (u), u ∈ (0, 1), (9)
where F i is the chosen i-th initial marginal. Following Definition 2.1 summaries the expression of non parametric estimators for both T and T i , i ∈ I.

Definition 2.1 (Non-parametric estimators of T and T i ). For a given arbitrary couple (x 0 , y 0 ) ∈ (0, 1) 2 , a non-parametric estimator of T is given by

T (x) = δ r(x) (y 0 ), for all x ∈ (0, 1) with r(x) such that δ 0 r(x) (x 0 ) = x,
where δ 0 r(x) refers to the self nested diagonal of the initial copula C 0 . In particular, if the initial copula C 0 is the independence copula, r(x) = 1 ln d ln -ln x -ln x0 . For any i ∈ I, non-parametric estimators T i are

T i (x) = F i • F -1 i • T (x).
A summary of the steps for the non-parametric estimation of T and T i is given in Algorithm 2.

Algorithm 2 Non-parametric estimation of T and T i

Input parameters

Choose (x 0 , y 0 ) in (0, 1) 2 , initial copula C 0 and initial marginals F i , for i ∈ I Choose k max pre-computation range

Eventual pre-calculations Get δ 1 and δ -1 by Algorithm 1 Calculate and store δ k (y 0 ), k ∈ {-k max , . . . , 0, . . . , k max } by Equation ( 7)

Non-Parametric estimation Get T by Equation ( 6), using previous precomputations Get T i by Equation ( 9)

Subset of points of transformations

In practice, we will propose in Section 3 some parametric estimators for the transformations T and T i , by requiring that these transformations are passing through a finite set of points. The interested reader is also referred to Di Bernardino and Rullière (2013a). We aim here at determining some good set of points to be chosen.

Firstly, we define some sets of points to be chosen, starting from given sets of quantile levels. The chosen expressions for these sets is motivated by several interesting properties (see Proposition 2.2 below).

Definition 2.2 (Set Ω). Let J ⊂ N be a finite set of indexes and Q = q (T ) j j∈J be a given set of targeted percentiles, q

(T ) j ∈ (0, 1), j ∈ J. One defines Ω(Q) = {(α j , β j } j∈J , with α j = q (T ) j β j = T (α j ),
where T is the estimator in Equation (6).

Definition 2.3 (Sets Ω i , i ∈ I). Let J i ⊂ N be finite sets of indexes and Q i = q (i) j j∈Ji
be finite given sets of targeted percentiles, q (i) j

∈ (0, 1), j ∈ J i , i ∈ I. For a given transformation τ , one defines the sets

Ω i (Q i , τ ) = (α i j , β i j ) j∈Ji where α i j = τ -1 (q (i) j ), β i j = F i • F -1 i (q (i) j ), (10)
where F -1 i is an estimator of the target i-th marginal distribution, denoted by F i .

Given Ω and Ω i , for i ∈ I, some set of points as in Definitions 2.2 and 2.3, we firstly derive some properties of some parametric transformations T Ω and T Ωi , i ∈ I that are chosen to pass through these sets Ω and Ω i , for i ∈ I (see Proposition 2.2 below). Theses properties justify the choice of sets Ω and Ω i in our estimation procedure.

Proposition 2.2 (Set of points for T i ). Assume that estimators T and F i are given continuous and invertible functions. Let J i ⊂ N be finite sets of indexes and

Q i = q (i) j j∈Ji
be a finite given set of targeted percentiles,

q (i) j ∈ (0, 1), j ∈ J i , i ∈ I. Then for transformations T i defined in Equation (9), i ∈ I T i is passing through all points of Ω i (Q i , T ),
where Ω i (Q i , T ) is defined as in Definition 2.3. Furthermore, for any invertible transformation T Ωi passing through all points of

Ω i (Q i , T ), i ∈ I F -1 Ωi (q) = F -1 i (q) for all q ∈ Q i , (11)
where

F Ωi = T • T -1 Ωi • F i is the transformed marginal distribution which uses T Ωi . Proof: Let i ∈ I, for any q ∈ Q i , one can define (α, β) ∈ Ω i (Q i , T ) by Equation (10), with α = T -1 (q) and β = F i • F -1 i (q).
We assume here that all estimators are continuous and invertible functions, so that for example

T • T -1 = Id. Firstly, one can check that T i (α) = T i • T -1 (q) and from Equation (9), T i (α) = F i • F -1 i • T • T -1 (q) = F i • F -1 i (q) = β, so that T i (α) =
β and the first result holds. Secondly, assuming all estimators are invertible, one can write

F -1 Ωi = F -1 i • T Ωi • T -1 , so that F -1 Ωi (q) = F -1 i • T Ωi (α). By assumption T Ωi (α) = β since T Ωi is passing through all points of Ω i (Q i , T ), so that finally F -1 Ωi (q) = F -1 i (β) = F -1 i • F i • F -1 i (q) and F -1 Ωi (q) = F -1 i (q)
which gives the second result.

Once chosen the thresholds q (i) j , j ∈ J i , for which we want to identify the transformed margins with targeted margins, one thus get a finite set of passage points for T i . Indeed, Proposition 2.2 shows that, using Ω i from Definition 2.3, one can select a further parametric estimator T Ωi passing trough points of Ω i and such that quantiles of the target F i are identified with quantiles of the transformed margin using T Ωi (see Equation ( 11)), for each quantile level q ∈ Q i and for any i ∈ I. Figure 1 illustrates this identification of quantile levels. Here an univariate data-set and the set Q = {0.1, 0.5, 0.9} are chosen. The black line represents the empirical cdf of data, F , and the red one the cdf F using the procedure in Proposition 2.2. As proved before, the quantiles of F are identified with quantiles of the obtained transformed F for each quantile level q ∈ Q (black points and associated dotted lines).

Such a result would be difficult to establish for the external transformation T . A weaker form is given by following Proposition 2.3.

Proposition 2.3 (Set of points for T ). Let J ⊂ N be a finite set of indexes and Q = q (T ) j j∈J be a given set of targeted percentiles, q

(T ) j ∈ (0, 1), j ∈ J. Define Ω(Q) as in Definition 2.2. Then obviously T is passing through all points of Ω(Q) . Furthermore, if T Ω is another transformation passing through all points of Ω(Q), then max (αj ,βj )∈Ω C(β j , . . . , β j ) -C Ω (β j , . . . , β j ) ≤ sup u∈[0,1] T (u) -T Ω (u)
where

C Ω (u 1 , . . . , u d ) = T Ω • C 0 (T -1 Ω (u 1 ), . . . , T -1 Ω (u d )) is the transformed copula using transformation T Ω , C Ω (u 1 , . . . , u d ) = T • C 0 ( T -1 (u 1 ), . . . , T -1 (u d ))
is the transformed copula using the full non-parametric transformation T and C 0 is the initial copula.

Proof: From Definition 2.2, it is obvious that T is passing through all points of Ω(Q). For any

(α j , β j ) ∈ Ω(Q), C Ω (β j , . . . , β j ) = T Ω •C 0 (T -1 Ω (β j ), . . . , T -1 Ω (β j )) and C Ω (β j , . . . , β j ) = T •C 0 ( T -1 (β j ), . . . , T -1 (β j )). Since both T and T Ω are passing trough all points of Ω(Q), it follows that T -1 Ω (β j ) = T -1 (β j ) = α j . As a consequence, C(β j , . . . , β j ) -C Ω (β j , . . . , β j ) = T (v j ) -T Ω (v j ) , with v j = C 0 (α j , . . . , α j ), v j ∈ [0, 1]. Hence the result.
Proposition 2.3 gives an interpretation of the set Q defined in Definition 2.2. Indeed, Q can be interpreted as a set of percentiles for which we want to minimize the difference between diagonals of the transformed nonparametric copula and a transformed parametric copula using a transformation T Ω instead of T .

Once obtained these finite sets of reaching points for T i and T , one can find parametric estimators without any optimization procedure. A summary of the estimation procedure for obtaining piecewise linear estimators of T and T i , i ∈ I, is given in Algorithm 3.

Algorithm 3 Subsets from non-parametric estimation

Input parameters Choose Q and Q i , i ∈ I, initial set of quantile levels, e.g. Q = Q1 = . . . = Q d = {0.25, 0.5, 0.75}.

Subsets of points Get T and Ti

, i ∈ I, by Algorithm 2 Get Ω(Q) by Definition 2.2 Get Ω i (Q i , T ) by Definition 2.3
We have seen that, starting from given thresholds q (i) j ∈ (0, 1), it was possible to propose some set of points Ω i , i ∈ I and Ω, as in Proposition 2.2 and 2.3, such that each T i is passing trough points of Ω i , and each T is passing trough points of Ω. Starting from these considerations, one can also introduce in the following definition for some piecewise linear estimators of T and T i , i ∈ I.

Definition 2.4 (Piecewise linear estimators of T and T i ). One can define two piecewise linear estimators of the external and internal transformations T and T i , i ∈ I.

• T PL is defined as a piecewise linear function passing trough the points (0, 0), all points of Ω(Q), and (1, 1),

where Ω(Q) is given in Definition 2.2.

• for each i ∈ I, T PL i is defined as a piecewise linear function passing through the points (0, 0), all points of Ω i (Q i , T PL ), and (1, 1), where Ω i (Q i , T PL ), i ∈ I is given in Definition 2.3.

Piecewise linear estimators T PL and T PL i , i ∈ I, are estimators relying on a chosen finite number of parameters. However, these estimators are not differentiable everywhere on (0, 1). In the following, we will try to find differentiable parametric estimators passing through the points of considered subsets Ω and Ω i , i ∈ I.

Parametric estimation

Composite hyperbolic transformations

As initially defined in [START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF] and [START_REF] Bienvenüe | On hyperbolic iterated distortions for the adjustment of survival functions[END_REF] we recall here the definition of hyperbolic composite transformations. Definition 3.1 (Conversion and transformation functions). Let f any bijective increasing function, f : R → R. It is said to be a conversion function. Furthermore the associated transformation function T to f is defined by:

T f : [0, 1] → [0, 1] such that T f (u) =    0 if u = 0, logit -1 (f (logit(u))) if 0 < u < 1, 1 if u = 1.
Remark that T f is a continuous non-decreasing function, such that T f (0) = 0, T f (1) = 1. Furthermore we remark that the transformation functions in Definition 3.1 are chosen in a way to be easily invertible. In particular in a way such that

T f • T g = T f •g , T -1 f = T f -1 .
When f is easily invertible, these readily invertible transformations help sampling transformed distributions (see [START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF].

In this section we consider the following particular class of hyperbolic conversion function (for further details see [START_REF] Bienvenüe | On hyperbolic iterated distortions for the adjustment of survival functions[END_REF]).

Definition 3.2 (A class of hyperbole). The considered hyperbole H is H m,h,ρ1,ρ2,η (x) = m -h + (e ρ1 + e ρ2 ) x -m -h 2 -(e ρ1 -e ρ2 ) x -m -h 2 2 + e η-ρ 1 +ρ 2 2 ,
with m, h, ρ 1 , ρ 2 ∈ R, and one smoothing parameter η ∈ R.

After some calculations, one can check that

H -1 m,h,ρ1,ρ2,η (x) = H m,-h,-ρ1,-ρ2,η (x).
Functions H in Definition 3.2 are thus readily invertible: a simple change of the sign of some parameters leads to the inverse function. As a consequence, transformations T H based on conversion functions H will also be readily invertible. In the following, we consider the generic hyperbolic conversion function defined in Definition 3.2. First remark that when the smoothing parameter η tends to -∞, the hyperbole H tends to the angle function:

A m,h,ρ1,ρ2 (x) = m -h + (x -m -h) e ρ1 1 {x<m+h} + e ρ2 1 {x>m+h} . ( 12
)
As remarked in [START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF], it thus appears that hyperbolic transformations have the advantage of being smooth versions of angle functions. They show in their paper that initial parameters for the estimation are easy to obtain with angle compositions. Another advantage of the consider hyperbolic transformations is the flexibility in the tail parameters estimation. We discuss this point in Remark 1.

Remark 1 (Tail behavior). In Di Bernardino and Rullière (2014), it is shown that for the hyperbolic transformations defined above (see Definition 3.2) and starting from some particular initial Archimedean copulas C 0 , it is possible to produce Archimedean copulas having tunable regular variation properties, and thus to get specific targeted multivariate lower and upper tail coefficients. Indeed using the hyperbolic conversion functions H, the aforementioned paper proposes a generic way to construct families of Archimedean generators presenting a chosen couple of lower and upper tail coefficients. This construction is based on the fact that the conversion function H in Definition 3.2 has an asymptote ax + b at -∞ with a = e ρ1 , and an asymptote αx + β at +∞ with α = e ρ2 . Illustrations proposed in Section 4 in the above mentioned article show that, when fitting some data, it is thus possible to propose a fit that respects some estimated tail dependence coefficients, by deducing parameters ρ 1 and ρ 2 from tail coefficients and by estimating other parameters m, h and η.

For sake of clarity, we recall below some definitions in Di Bernardino and Rullière (2013a). These definitions of composite transformations and suited parameters will be useful in the following.

Definition 3.3 (Composite transformations). Let k ∈ N. Consider η ∈ R and a given parameter vector θ = (m, h, ρ 1 , ρ 2 , a 1 , r 1 , . . . , a k , r k ) if k ≥ 1, or θ = (m, h, ρ 1 , ρ 2 ) if k = 0.
We define the angle composite transformation A θ as:

A θ = T f θ , with f θ = A a k ,0,0,r k • • • • • A a1,0,0,r1 • A m,h,ρ1,ρ2 if k ≥ 1, A m,h,ρ1,ρ2 if k = 0,
and the hyperbolic composite transformation H θ,η as:

H θ,η = T f θ,η , with f θ,η = H a k ,0,0,r k ,η • • • • • H a1,0,0,r1,η • H m,h,ρ1,ρ2,η if k ≥ 1, H m,h,ρ1,ρ2,η if k = 0,
with A m,h,ρ1,ρ2 as in Equation ( 12) and H m,h,ρ1,ρ2,η as in Definition 3.2.

Remark 2 (Inverse composite transformations). Let k ∈ N. Consider η ∈ R and a given parameter vector

θ = (m, h, ρ 1 , ρ 2 , a 1 , r 1 , . . . , a k , r k ) if k ≥ 1, or θ = (m, h, ρ 1 , ρ 2 ) if k = 0. Since T -1 f = T f -1 , the angle composite transformation A -1
θ is such that:

A -1 θ = T f θ , with f θ = A m,-h,-ρ1,-ρ2 • A a1,0,0,-r1 • • • • • A a k ,0,0,-r k if k ≥ 1, A m,-h,-ρ1,-ρ2 if k = 0.
The hyperbolic inverse composite transformation H θ,η is such that:

H -1 θ,η = T f θ,η , with f θ,η = H m,-h,-ρ1,-ρ2,η • H a1,0,0,-r1,η • • • • • H a k ,0,0,-r k ,η if k ≥ 1, H m,-h,-ρ1,-ρ2,η if k = 0.
Definition 3.4 (Suited parameters from Ω). Let k ∈ N. Consider one given set Ω = {ω 1 , . . . , ω 3+k }, ω j ∈ (0, 1) 2 . Denote by u j and v j the two respective components of each ω j in the logit scale, such that ω j = (logit -1 u j , logit -1 v j ), j ∈ {1, . . . , 3 + k}. Assume that u j and v j are increasing sequences of j. We define:

Θ(Ω) = (m, h, ρ 1 , ρ 2 , a 1 , r 1 , . . . , a k , r k ) if k ≥ 1, (m, h, ρ 1 , ρ 2 ) if k = 0, where m = u2+v2 2 , h = u2-v2 2 , ρ 1 = ln v2-v1 u2-u1 , ρ 2 = ln v3-v2 u3-u2 , r k = ln v 3+k -v 2+k u 3+k -u 2+k u 2+k -u 1+k v 2+k -v 1+k , a k = v 2+k , k ≥ 1.
Proposition 3.1 (Suited composite transformations). Let k ∈ N. Consider one given set Ω = {ω 1 , . . . , ω 3+k }, ω j ∈ (0, 1) 2 and a smoothing parameter η ∈ R. Set θ = Θ(Ω), then -the transformation A θ (x) is piecewise linear in the logit scale and will be called logit-piecewise linear. It links point (0, 0), points of Ω, and point (1, 1).

-the transformation H θ,η converges pointwise to A θ as η tends to -∞. It results that the continuous and differentiable transformation H θ,η can fit as precisely as desired the set of points Ω when η tends to -∞.

Proof: The first result is proved in [START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF]. It simply comes from the fact that A Θ(Ω) (u j ) = v j for all j ∈ {1, . . . , 3 + k}, where ω j = (logit -1 u j , logit -1 v j ). The convergence of the hyperbole composite transformation toward the angle composite transformation is straightforward and also evoked in [START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF].

Smoothed parametric estimators

For the estimation of transformations T and T i , we assume that are given:

• Initial estimators of δ and δ -1 (see Section 2.1).

• The sets of quantile levels

Q = q (T ) j j∈J , Q i = q (i) j j∈Ji
, i ∈ I (see Section 2.3).

• Some smoothing parameters η ∈ R and η i ∈ R, i ∈ I (see Section 3.1).

In the following we introduce the smooth parametric estimators T and T i , i ∈ I, for external and internal transformations respectively.

Estimation of T Using estimators of δ and δ -1 (see Section 2.1), one easily gets the resulting set Ω(Q) by Proposition 2.3 and suited parameters by Definition 3.4:

θ = Θ( Ω(Q)). ( 13
)
Then one defines: T = H θ,η .

Estimation of T i , i ∈ I Once T estimated by T , one gets the passage set Ω i , for i ∈ I by Proposition 2.2 and suited associated parameters θ i by Definition 3.4, i.e.,

θ i = Θ( Ω i (Q i , T ))). (14)
Note that once given thresholds sets Q and Q i and smooth parameters η and η i , all estimated parameters are directly and analytically defined, so that we do not need here any inversion or optimization procedure. Then one defines:

T i = H θi,ηi , for i ∈ I.
A summary of the expressions of smooth estimators T and T i , i ∈ I, is given in following definition.

Definition 3.5 (Smooth estimation of T and T i , i ∈ I). Let Q and Q i , i ∈ I be given sets of quantile levels. Let η ∈ R and η i ∈ R, i ∈ I be given smoothing parameters. One defines

T = H θ,η , T i = H θi,ηi , i ∈ I, (15)
where parameters θ and θ i are given by

θ = Θ( Ω(Q))), θ i = Θ( Ω i (Q i , T )), i ∈ I.
The sets Ω and Ω i are defined in Definitions 2.2 and 2.3. The function Θ is defined in Definition 3.4.

Estimation procedure for obtaining smoothed parametric estimators T and T i , i ∈ I is gathered in Algorithm 4.

Algorithm 4 Parametric estimation

Input parameters

Choose Q and Q i , i ∈ I the sets of quantile levels, e.g. Q = Q1 = . . . = Q d = {0.25, 0.5, 0.75}, Choose smoothing parameters η ∈ R and η i ∈ R, i ∈ I,

Parametric estimation

Get Ω(Q)) and Ω i (Q i , T ), i ∈ I by Algorithm 3 Get θ by Equation ( 13) Get θ i by Equation ( 14) Get smooth estimators T and T i , i ∈ I by Equation ( 15).

One can define the complete vector parameter presented in Algorithms 3 and 4:

Θ Θ Θ = ( θ 1 , . . . , θ d , θ, η 1 , . . . , η d , η).
What is noticeable here is that given thresholds Q, Q i and smoothing parameters η, η i , one have direct expressions for the parametric estimators θ and θ i , i ∈ I. The estimation can also change, depending on some other estimation choices (generator among its equivalence class (x 0 , y 0 ), bandwidths h i , i ∈ I or kernel function K), but one aims here at finding good estimators whatever the choice of any reasonable value of these parameters.

In numerical Section 6 we will illustrate this point. Finding a global optimum in high dimension would lead to a curse a dimensionality and numerical problems. As a consequence, it is very important to get correct estimators without jointly optimizing a lot of parameters in high dimension.

Final parametric results

Once all parameters estimated as in Section 3, the previous parametric model allows to get various analytical results for both the multivariate distribution function, its associated critical layers, Kendall's function and multivariate return periods.

Firstly we obtain the parametric expression for the transformed copula C:

C(u 1 , . . . , u d ) = φ( φ -1 (u 1 ) + . . . + φ -1 (u d )), ( 16
)
where φ is the final parametric transformed Archimedean generator (see ( 4)). This generator can be easily written in terms of the external transformation T , i.e.,

φ(t) = T (φ 0 (t)), (17)
where φ 0 is the generator associated to the initial copula C 0 in (4) and T is the smooth estimator of the external transformation presented in Section 3.2.

Using model in (4), the corresponding estimated transformed multivariate distribution is given by:

F Θ Θ Θ (x 1 , . . . , x d ) = H θ,η • C 0 (H -1 θ1,η1 • F 1 (x 1 ), . . . , H -1 θ d ,η d • F d (x d )), (18)
where Θ Θ Θ is the complete estimated vector parameter presented in Algorithms 3 and 4.

Furthermore, using expression in (18) for the estimated transformed multivariate distribution F , the associated parametric α critical-layers are given by

∂ L Θ Θ Θ (α) = {(F -1 1 • H θ1,η1 (u 1 ), . . . , F -1 d • H θ d ,η d (u d )), (u 1 , . . . , u d ) ∈ (0, 1) d , C 0 (u 1 , . . . , u d ) = H -1 θ,η (α)},
where a direct analytic expression

H -1 θ,η
is given by Remark 2 (see also Proposition 2.4 in Di Bernardino and Rullière (2013a)).

As presented in the introduction, we aim at providing a parametric estimation of the multivariate Return Period α) . [START_REF] Genest | On the multivariate probability integral transformation[END_REF] obtained the following explicit expression for the Kendall distribution in the case of multivariate Archimedean copulas with a given generator ϕ, i.e.,

RP > (α) = ∆ t • 1 1-K C (
K C (α) = α + d-1 i=1 1 i! -ϕ -1 (α) i ϕ (i) ϕ -1 (α) , for α ∈ (0, 1),
where the notation f (i) corresponds to the i-th derivatives of a function f .

Consider the case of transformed generator φ(t) = T (φ 0 (t)), in ( 16)-( 17). Then, the analytical formula for the estimated transformed Kendall distribution K C can be easily written as:

K C (α) = α + d-1 i=1 1 i! -φ -1 0 (T -1 (α)) i T • φ 0 (i) φ -1 0 (T -1 (α) , for α ∈ (0, 1). ( 19
)
Remark that the hyperbolic transformation T = H θ,η , in this paper are chosen in a way to be easily invertible.

Then estimated smooth inverse transformation T -1 in ( 19) is straightforwardly obtained by changing the signs of some parameters of the hyperbolic composition (see Remark 2).

Comprehensive scheme of the estimation procedure

The method presented in this paper involves different notations progressively introduced in different sections. Moreover, the final scheme of the procedure relies in many tuning parameters. To improve clarity, a comprehensive scheme of the estimation procedure is presented below (see Table 1). The relationship between the different algorithms presented above is pointed out using arrows and implications. Moreover, the nature of considered parameters is specified, in particular, we distinguish between tuning parameters and estimated/calculated quantities.

As one can see in Table 1 (left column), the algorithm depends on many tuning parameters, we discuss here the sensitivity of the results on these choices. The results expressed below refer to variations of mean absolute errors and are based on our simulation study (see Section 8, case n = 1000, using as a reference

x 0 = y 0 = e -1 , Q = Q 1 = Q 2 = {0
.25, 0.5, 0.75}, η = -1, η 1 = η 2 = -3, with very small bandwidths).

• We always used in our numerical illustrations (see Sections 6, 7 and 8), G -1 = logit(x) in order to simplify inversion procedures and C 0 the independent copula. We believe that these choices might have more important implications for the tails. A deep analysis in this sense is developed in Di Bernardino and Rullière (2014).

• The use of positive bandwidth sizes h i ensures continuity and possible inversion of fitted copula diagonal section and marginals distributions. Using small bandwidth sizes may increase the variance of some non-parametric estimators, but in practice this leads to estimate parametric copula and margins using empirical diagonal and empirical margins instead of smoothed ones. The impact is thus limited on the parametric estimation. As an example in further simulation study (see Section 8), using a bandwidth given by classical rule of thumb or dividing this bandwidth by 100 leads to an absolute variation of 0.6% on resulting mean absolute errors.

• Concerning the choice of the point (x 0 , y 0 ), we observed absolute variations of order 4% in mean absolute errors compared to the choice of x 0 = y 0 = e -1 (which simplifies some calculations).

• Concerning the choice of quantile levels set Q, Q i , for i ∈ I, using regular thresholds

Q = Q 1 = Q 2 = {0
.25, 0.5, 0.75}, without any optimization, leads to mean absolute error of M AE = 3.6% (instead of approx. 1% for optimized thresholds), which is still reasonable for non optimized parameters.

• At last, choices of η and η i , i ∈ I, mostly depend on the desired smoothness of the final resulting distribution. Except for clearly excessive smoothing, the impact on resulting distribution is quite small. As an example on simulated data of Section 8, using η = 0 and η 1 = η 2 = -2 or using η = -2 and η 1 = η 2 = -4 leads to absolute variations of 0.5% on resulting mean absolute errors.

Scheme of the estimation procedure Chosen tuning parameters in Algorithm Calculated quantities

1) G -1 (inverse of a cdf ) 2) h i , i ∈ I bandwidth sizes 3) Point (x 0 , y 0 ) in (0, 1) 2 4) Initial copula C 0 5) Initial marginals F i , for i ∈ I 6) Q and Q i , i ∈ I (initial set of quantile levels) 7) η ∈ R and η i ∈ R, i ∈ I ( smoothing
Ω i (Q i , T ), for i ∈ I g) Using 7
), e) and f ),

• Get θ = Θ( Ω(Q)). • Get θ i = Θ( Ω i (Q i , T )) for i ∈ I , where T = H θ,η .
h) Using 4), 5) and g), finally get

F Θ Θ Θ (x 1 , . . . , x d ) = H θ,η • C 0 (H -1 θ1,η1 • F 1 (x 1 ), . . . , H -1 θd,ηd • F d (x d )), where Θ Θ Θ = ( θ 1 , . . . , θ d , θ, η 1 , . . . , η d , η).
Table 1: Comprehensive scheme of the global estimation procedure. We distinguish between tuning parameters (left column) and estimated/calculated quantities (right column). This quite small sensitivity of these parameters is important. Indeed in high dimension, due to the curse of dimensionality, it is not possible in practice to find the global optimum of a function. It is thus essential that any chosen tuning parameter lead to small average errors, in order that calculated quantities are near a local optimizer of the considered error, and in order to reduce the difference between local optima and global one. Resulting calculated parameters like (m, h, ρ 1 , ρ 2 ) for each transformation can still be optimized once they are close to one optimizer.

In further numerical illustrations (see Sections 6, 7 and 8), we never optimize resulting calculated quantities, but we sometimes choose tuning parameters that permit to get reduced final errors.

Numerical results on the rainfall real data

In the following we illustrate our methodology presented in the previous sections using a 5-dimensional rainfall data-set.

Presentation of the data-set

Data comes from the website CISL Research Data Archive (RDA), http://rda.ucar.edu, and is available for registered users. The user is granted the right to use the Site for non-commercial, non-profit research, or educational purposes only, without any fee or cost, as specified in the terms of use of the website (January 2014). The name of the data in this website is: ds570.0, and the direct url to this data is http://rda.ucar.edu/ datasets/ds570.0/. The whole citation information for this data is: We detail here data selection in order to help researchers to retrieve the same dataset. We have taken all monthly rainfall data available from CSV files for five chosen stations of India and Sri-Lanka. After a first importation step, only numerical or missing values of the field Precip(mm) were kept, and the field Date was considered as a numerical primary key in order to avoid repeated data. Considered stations and corresponding lines number after this step are: In a second step, all five stations were grouped by date, and we have selected dates for which all fields Precip (mm) of the five station were non-missing and non-zero (the suppression of zero precipitation data aims at easing the parametric representation of margins, as detailed in [START_REF] Koning | Are precipitation levels getting higher? Statistical evidence for the Netherlands[END_REF]). At last precipitations have been expressed in decimetres: all fields Precip (mm) have been divided by 100. As a result, the data has 797 lines giving monthly precipitation in decimetres, for some dates in the period from 1893-01 to 2013-06. Excluded dates in this period are those for which at least one field Precip (mm) was missing, zero or non-numerical.

As one can see part of the data is very old and should certainly be interpreted with caution. We have drawn on Figure 2 the autocorrelation functions to show that a small autocorrelation is still remaining in considered data. Indeed as remarked in [START_REF] Reiss | Statistical analysis of extreme values with applications to insurance, finance, hydrology and other fields[END_REF], on the one hand, considering annual data is a way to avoid the problems of serial correlation and seasonal variation. On the other hand, it may represent a lost of information contained in the data. A compromise can be to base the inference on seasonal or monthly maxima (see [START_REF] Reiss | Statistical analysis of extreme values with applications to insurance, finance, hydrology and other fields[END_REF], Section 14.1). This is the approach followed in this paper. Remark that there is some significant seasonality in the considered rainfall real data (see Figure 2). To overcome this type of problem it could be possible to fit the following model : X i,t = a i S t + Y i,t , for i = {1, . . . , 5} where S t represents some cycle (that might be non-parametric) and Y i,t represents the deviation to the cycle. However this type of study is beyond the scope of the present paper. Geographical position of 5 stations and the scatter plot of ranks of data are provided in Figure 3. Remark that the numerical illustrations presented here aims at showing the feasibility of the estimation but do not aim at furnishing a complete hydrological study, which would require more data treatments for handling seasonality, analysis of peak and durations of rainfall, etc. The interested reader is referred to [START_REF] Salvadori | On the return period and design in a multivariate framework[END_REF][START_REF] Salvadori | Semi-parametric approximation of Kendall's distribution function and multivariate Return Periods[END_REF], [START_REF] Gräler | Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation[END_REF].

Estimation results

We consider the model as in Equation ( 3). For the sake of clarity, each transformation is here involving only one hyperbola, thus requiring the choice of three quantiles thresholds per hyperbola for the estimation. We take as initial copula C 0 the independent one, and the initial margins F i (x) = 1 -e -x , i ∈ I. Concerning the impact of the choice of the initial copula C 0 on transformed generator φ, the interested reader is referred to Proposition 3.12 in Di Bernardino and Rullière (2013b).

For the estimation, we have chosen very small smoothing bandwidths h i , i ∈ I (see Equation ( 5)). Diagonal section of the copula and smooth empirical margins are thus very close to empirical ones, but small smoothing allows to get proper invertible functions. As stated in Algorithm 2, the estimation relies on an arbitrary chosen initial point (x 0 , y 0 ), which corresponds to the arbitrary choice of one point of the Archimedean generator among its equivalence class (see Remarks 7 and 8 in Di Bernardino and Rullière (2013b)). The choice of this point has an impact on the estimation, and given an arbitrary abscissa x 0 = 0.5, we have kept the value of y 0 giving the best results, here y 0 = 0.24. The non-parametric estimation relies on the choice of quantile thresholds (see Algorithm 3). Here we have chosen:

Q = {5%, 50%, 95%} and Q i = {20%, 50%, 80%} , i ∈ I = {1, . . . , 5} .
At last, for the parametric estimation (see Algorithm 4), one have to choose smoothing parameters. Selected smoothing parameters are: η = -1 and η i = -3, i ∈ I = {1, . . . , 5} , (see Table 2 below). It would be naturally possible to optimize all these parameters, Q, Q i , η, η i , for i ∈ I, but we have chosen here fixed threshold and smoothing parameters in order to show the feasibility of the estimation.

This also shows that it is possible to get good fits with non optimized parameters, which is important when the dimension is high, since a global optimization procedure would have to face a curse of dimensionality.

As presented in Sections 2 and 3, we obtain the complete estimated vector of parameters Θ Θ Θ given in Table 2.

Then we get the transformed multivariate copula C and distribution function F Θ Θ Θ , obtained by Equations ( 16) and ( 18) respectively.

We perform a goodness-of-fit test based on the empirical process in order to test the quality of the adjustment of copula C on these multivariate data. In the large scale Monte Carlo experiments carried out by [START_REF] Genest | Goodness-of-fit tests for copulas: A review and a power study[END_REF], the statistic S n gave the best results overall (see Section 4 in Ivan [START_REF] Kojadinovic | Modeling multivariate distributions with continuous margins using the copula R package[END_REF]). An approximate p-value for S n can be obtained by means of a parametric bootstrap-based procedure (see Section 4.1 in Ivan Kojadinovic and Jun Yan ( 2010)), and whose validity was recently shown by [START_REF] Genest | Validity of the parametric bootstrap for goodness-of-fit testing in semiparametric models[END_REF]. In order to apply the bootstrap-based procedure, we need to generate random samples from the transformed copula C with generator as in (17). To this aim we use the Marshall and Olkin's algorithm with a numerical inversion of the Laplace Transform of generater φ using the Talbot method. We obtain a p-value= 0.37129. Furthermore, we test different competitor copula families, with maximum likelihood estimated parameters.

Parameters Θ Θ Θ m h ρ 1 ρ 2 η θ for external T -0.
Obtained p-values, with different goodness-of-fit tests, are gathered in Table 3. Therefore, among all the copula families that we have tested, the transformed copula C is the only one that is not rejected at the 5% significance level (see Ivan [START_REF] Kojadinovic | Modeling multivariate distributions with continuous margins using the copula R package[END_REF]).

Copula 2009)) for competitor copula families on the considered 5-dimensional rainfall data, with n = 797. In all cases, the number of Monte Carlo experiments is fixed at N = 1000.

To appreciate the quality of multivariate parametric adjustments we evaluate the Supremum Absolute Error on a lattice G (see ( 20)) and the Mean Absolute Error on the data (see ( 21)), respectively defined as:

SAE = sup (x1,x2,x3,x4,x5)∈G |F (x 1 , x 2 , x 3 , x 4 , x 5 ) -F n (x 1 , x 2 , x 3 , x 4 , x 5 )| , (20)
MAE = 1 n n k=1 F (X (k) 1 , X (k) 2 , X (k) 3 , X (k) 4 , X (k) 5 ) -F n (X (k) 1 , X (k) 2 , X (k) 3 , X (k) 4 , X (k) 5 ) , (21)
with n = 797, F n the 5-dimensional empirical distribution function and F a parametric model on the multivariate rainfall data. The first SAE criterion is a classical Kolmogorov-Smirnov statistic evaluated on a particular lattice G, and the same choice of an absolute value (L 1 -norm) has been done for the second MAE criterion in order to use comparable norms, and to get results in accordance with further Figures 456showing absolute differences. This second criterion can be linked to L 1 -variant Cramér-von-Mises distances like

|F (x) -F n (x)| dF (x)
, where the continuous measure dF (x) has been replaced by the empirical one dF n (x) thus leading to a sum. This way, the MAE criterion can be efficiently computed without integration on a multivariate domain, while the difference between the sum and the integral can be bounded using the difference

|F (x) -F n (x)|.
We consider our transformed model (see ( 18)) and three different parametric models using Frank, Gumbel and Clayton copulas and parametric marginals. The dependence parameters of copulas, fitted by Maximum likelihood, are: θ = 2.45 (Frank copula), θ = 1.33 (Gumbel copula) and θ = 0.47 (Clayton copula). For marginals we have tested 15 different classes of distributions and we have fitted the best model using the Akaike Information Criterion. Following this criterion, the best parametric marginals were Gamma distributed. We recall here that zero valued precipitations have been excluded from this data. Results are gathered in Table 4.

Gamma distributions

X 1 X 2 X 3 X 4 X 5 shape
1.7553 0.6344 0.9600 1.0373 1.0116 rate 0.7994 0.5968 0.7844 0.6724 0.5799 Table 4: Shape and rate MLE parameters of fitted marginal Gamma distributions for X i for i = 1, . . . , 5.

Then the obtained Supremum Absolute Error in (20) and Mean Absolute Errors in ( 21) are given in Table 5. As one can see, our model performs better both on the Supremum Absolute Error (SAE) and the Mean Absolute Error (MAE) criteria. However, quantifying the goodness of fit on a multivariate data is a difficult problem. Some models with small Mean Absolute Error criteria may behave poorly when considering a specific projection of the 5-dimensional space. A control of the performance of the model for the distribution fit of each pair of random variable is highly recommended. 20)-( 21) for the considered parametric 5-dimensional models. Best results are indicated in bold font.

Models

Using parametric multivariate models introduced above, we now consider the fit of the bivariate distributions F (Xi,Xj ) and of the marginals F Xi , for i, j = 1, . . . , 5. They are particular projections of the 5-dimensional distribution. We consider the errors in ( 20)-( 21) for (X i , X j ) data, for i, j = 1, . . . , 5, i.e.,

SAE i,j = sup (x,y)∈G F (Xi,Xj ) (x, y) -F n (x, y) , (22)
MAE i,j = 1 n n k=1 F (Xi,Xj ) (X (k) i , X (k) j ) -F n (X (k) i , X (k) j ) , (23)
with F n the bivariate empirical distribution function and F (Xi,Xj ) the projection of the multivariate parametric model on (X i , X j ). Errors in ( 22)-( 23), evaluated using our transformed model F and classical parametric models on the considered rainfall data, are gathered in Tables 6 and7. Remark that, without any optimization procedure, the transformed model performs better in terms of errors in ( 22)-( 23) for almost all the couples (X i , X j ). In order to make the reading of the Tables 6 and7 easier, we give in Table 8 the associated synthetic statistics. As we can seen the best values (displayed in bold font) are provided by the transformed model F . As remarked above, some graphical illustrations of Tables 6 and7 are provided in Figures 4 for (X 1 , X 4 ), Figure 5 for (X 2 , X 5 ) and Figure 6 for (X 3 , X 4 ). The maximal range for these figures corresponds to the 95th percentile of each random variable, in order to focus on main part of the data and to preserve the readability of each figure.

Errors model F X 1 X 2 X 3 X 4 X 5 X 1 0.0110 0.
Errors Gumbel copula and Gamma marginals 6 and7, where SSAE= 22)). Black dots represent the associated rainfall data (X 1 , X 4 ).

X 1 X 2 X 3 X 4 X 5 X 1 0.0119 0.
sup ij SAE ij , SMAE= sup ij MAE ij , MSAE= mean ij SAE ij and MMAE= mean ij MAE ij .
Furthermore, from (4), we get Fi = T

• T -1 i • F i , for i ∈ I.
Then, using the smooth estimation of external and internal transformations T and T -1 i , for i ∈ I, one can obtain the transformed parametric marginal distributions. Results for some margins are presented in Figure 7 below (results for other margins are completely analogous). Finally we consider another projections of F , i.e. the 5-dimensional diagonal. In Figure 8 (left) we present the parametric estimation of the 5-dimensional diagonal using the transformed model and the classical parametric competitors introduced above.

In Figure 8 (right) we present the Kendall distribution for the considered 5-dimensional rainfall data-set using our transformed generator φ(t) = T • exp(-t). We also display K C (α), for Gumbel, Frank and Clayton copulas. 22)). Black dots represent the associated rainfall data (X 2 , X 5 ). The vertical lines show estimates of 90% and 95% univariate marginal quantiles. Red lines the Gamma distribution with parameters as in Table 4. Green lines are ln(1 -F i ), for i = 1, 4, 5 as in the 5-dimensional transformed model with parameters as in Table 2.

7 A nested model on the rainfall real data

Choice of clusters

In this section we intend to show the flexibility of the proposed model and associated estimation procedure.

In particular, we adapt our methodology in the case of some asymmetric dependencies (as, for instance, nonexchangeable random vectors). The correlation matrix of the considered rainfall data is displayed in Figure 9 (left). As we can see, some pairs of stations present a bigger correlation. To illustrate how our model can be adapted to this situation we have decided to create 2 different clusters. We have grouped together pairs of variables presenting correlation greater than 66%, this leads to a first (tri-variate) cluster composed by stations (X 2 , X 3 , X 5 ).

        X 1 X 2 X 3 X 4 X 5 X 1 1 0.
The remaining second (bivariate) cluster is (X 1 , X 4 ). One can check that all correlations inside each cluster are greater than 50%. Figure 3 gives the geographical position of each station and helps visualizing each cluster. Furthermore an hierarchical cluster analysis on the set of dissimilarities produce by the distance of the X i is developed. We use different types of distance to create the dissimilarities (Euclidian, maximum, Manhattan, Canberra, Binary, Minkowski). In all cases we obtain the result in Figure 9. As one can see, whatever the distance chosen for dissimilarities, the dendrogram gives a justification to chosen clusters of station indexes {2, 3, 5} and {1, 4}. However, chosen nested structure has not been optimized; the choice of an optimal nested structure is analyzed for instance in [START_REF] Segers | Nonparametric estimation of the tree structure of a nested Archimedean copula[END_REF].

Then, following these considerations, we firstly fit a 3-dimensional model for the first group and a 2-dimensional one for the second one. We generate the pseudo-data coming from these two models and finally we construct the joint copula for these bivariate data-set. As in Section 6, in the following we take as initial copula C 0 the independent one, and the initial margins

F i (x) = 1 -e -x , i ∈ A, B.
The multivariate distribution for the cluster A = {2, 3, 5} is assumed to be written:

F A (x 2 , x 3 , x 5 ) = T A • C 0 (T -1 A • F 2 (x 2 ), T -1 A • F 3 (x 3 ), T -1 A • F 5 (x 5 )), (24)
with Fi =

T A • T -1 Ai • F i , for i ∈ A.
The multivariate distribution for the cluster B = {1, 4} is assumed to be written:

F B (x 1 , x 4 ) = T B • C 0 (T -1 B • F 1 (x 1 ), T -1 B • F 4 (x 4 )), (25)
with

Fi = T B • T -1 Bi • F i , for i ∈ B.
The whole 5-dimensional distribution is assumed to be written:

F (x 1 , x 2 , x 3 , x 4 , x 5 ) = T • C 0 T -1 • F A (x 2 , x 3 , x 5 ), T -1 • F B (x 1 , x 4 ) , (26)
where

C(u, v) = T • C 0 T -1 (u), T -1 (v)
is referred as the root copula at point (u, v). It is effectively a proper copula if the transformation T satisfies admissibility conditions that are given in Proposition 2.1 of Di Bernardino and Rullière (2013b), in order to satisfy 2-monotony as detailed in [START_REF] Mcneil | Multivariate archimedean copulas, d-monotone functions and l 1 -norm symmetric distributions[END_REF]. Model in ( 26) corresponds to a Nested Archimedean Copula model, with two nested levels, as described in [START_REF] Hofert | Densities of nested Archimedean copulas[END_REF]. In this article authors give also conditions such that the resulting nested copula is a proper distribution function. Despite it may not be the case in general, we will assume that F with expression as in ( 26) is a proper multivariate distribution. In the following we will check the admissibility for considered transformation (see Figure 12).

As one will see in Figure 13, conditions like the admissibility of used copulas and of the final nested distribution should be carefully checked. Here the copula C is used without uniform margins to create a multivariate distribution when the final nested density remains positive on the whole domain. Other constructions ensuring this admissibility could be investigated, like the use of Hierarchical Kendall Copula, as in Brechmann ( 2014). However such supplementary investigations are beyond the scope of this paper.

Remark that in the model presented here, for the sake of simplicity, the root copula C in (26) has only two arguments (two child copulas). The methodology detailed here is however applicable with more arguments. The hierarchical copula detailed here is also a two-step hierarchical copula, with only one imbrication level, but the methodology can be extended to more levels. The interest reader is referred to [START_REF] Hofert | Densities of nested Archimedean copulas[END_REF] or Brechmann (2014).

Estimation results

The parameters of each cluster copula transformation T A and T B in ( 24) and (25) will be respectively denoted by θ A and θ B . At last, parameter of marginal transformation T Ai will be denoted θ Ai , for i ∈ A, and parameter of marginal transformation T Bi will be denoted θ Bi , for i ∈ B. To estimate the external transformation T of model ( 26) we firstly construct a bivariate pseudo data-set:

Z 1 = F A (X 1 , X 4 ), Z 2 = F B (X 2 , X 3 , X 5 ).
Then we fit on this bivariate data-set a model

F (Z1,Z2) (z 1 , z 2 ) = T • C 0 (T -1 • F1 (z 1 ), T -1 • F2 (z 2 )), (27)
with Fi = T

• T -1 1 • F i , for i = 1, 2.
The parameter of transformation T in (27) will be denoted θ, parameters of transformation T 1 and T 2 will be respectively denoted by Both data-set (Z 1 , Z 2 ) and a graphical illustration of the fit of model ( 27) on this data using the estimated parameters above, are given in Figure 10. We evaluate the Mean Absolute Error in (21) for final model in (26) and we obtain MAE = 0.0088 and SAE = 0.0636. This value is smaller than all values in Table 5. Furthermore, as in Tables 6 and7, we consider the fit of the bivariate distributions F (Xi,Xj ) and of the marginals F Xi , for i, j = 1, . . . , 5. To this aim we evaluate the SAE i,j and MAE i,j errors in ( 22)-( 23) using the nested model in (26). Results are gathered in Table 9.

Remark that errors in Table 9 are smaller than all values previously obtained in Tables 67. Analogously to Figures 456, and for the same couples of random variables, graphical illustrations of Table 9 are provided in Figure 11.

As one can see in Table 9, results with this nested model are very good. In particular, one can see in Figure 11 and in Table 9 that absolute errors in every bi-dimensional projection remain very small. This model is however a simple illustration to feasible improvements of the initial model, when the dimension is greater than 2, without using heavy optimization algorithms.

As remarked above, one have to ensure that the Equations ( 24), ( 25) and ( 26) are proper distribution functions.

We verify here the admissibility for external transformations 6 and7 are indicated in bold font. 26) for (i, j)= (1, 4) (left), (i, j)= (2, 5) (centre panel), (i, j)= (3, 4) (right). Black cross represents the maximum error SAE i,j in the considered lattice (see ( 22)). Black dots represent the associated rainfall data (X i , X j ).

Consider an external transformation τ , as T A , T B or T in respective Equations ( 24), ( 25) and ( 26). Denote by d τ the dimension of this external transformation, and τ (i) the i-th derivative of the function τ , i = 1, . . . , d τ . For any considered external transformation τ , the transformation is admissible if and only if each quantity f τ,i (x) is nonnegative for x ∈ (0, 1) and for each i = 1, . . . , d τ , where :

   f τ,1 (x) = τ (1) (x) f τ,2 (x) = τ (1) (x) + x τ (2) (x), f τ,3 (x) = τ (1) (x) + 3 x τ (2) (x) + x 2 τ (3) (x).
In Figure 12 we have drawn, for x ∈ (0, 1), When functions f τ,i are continuous, i ≤ d τ , their logarithm tends to -∞ before f τ,i becomes negative, and a lower bound of the logarithm ensure that m τ is well-defined for x ∈ (0, 1), so that the quantity m τ helps checking the admissibility of τ , especially for very small values of f τ,i . As one can see in Figure 12, for each transformation T, T B , T A , the quantity in ( 28) is well-defined for each x ∈ (0, 1) and inferiorly bounded. Then, we can deduce the admissibility for considered external transformations T , T B and T A .

m τ = min {ln f τ,i (x), i = 1, . . . , d τ } , (28)
To illustrate the danger to use a single criterion, and the need to check admissibility conditions, we propose in Figure 13 an illustration of situations that can happen with deliberately non admissible transformations. In the left panel of the Figure 13, we have drawn the function m τ in (28) for a typical non-admissible external transformation τ (in this example associated to the cluster A, see ( 24)). What is noticeable is that this non-admissible transformation is however leading to a good value of global MAE and SAE criteria for the final nested model in ( 26), MAE = 0.0057 and SAE = 0.0654. These error values have to be compared with a MAE = 0.0088 and SAE = 0.0636 previously obtained using admissible transformations.

In the right panel of Figure 13, we have drawn errors |F (X3,X4) (x, y)-F n (x, y)|, for (x, y) in a lattice of 100×100 points, where F n is the empirical distribution function and F (X3,X4) is the parametric nested model in ( 26) using this non-admissible external transformation τ . The maximum error is SAE 3,4 = 0.0906 in the considered lattice (at the black cross point). Furthermore MAE 3,4 = 0.0201. These values are larger then associated values in Table 9, i.e., SAE 3,4 = 0.0351 and MAE 3,4 = 0.0100. Despite good global criteria MAE and SAE, admissibility conditions are not fulfilled, and projected criteria MAE 3,4 and SAE 3,4 are disappointing. This shows that projected criteria may behave differently than global criteria, and that the admissibility conditions have to be checked carefully to avoid undesirable behavior of the nested adjustment.

Critical Layers for nested model

Let α ∈ (0, 1) be a targeted level for a critical layer. Let C 0 be the initial copula to be transformed, and assume that C 0 is the independent copula.

Models

Conclusion

We described an estimation procedure for multivariate distribution functions. This methodology provides also parametric expressions of associated quantities as critical layers, Kendall's function, return periods. The considered model is based on transformations of the marginals and of the dependence structure within the class of Archimedean copulas. The proposed estimation is straightforward, it has a tunable number of parameters and it does not rely on any optimization procedure. Furthermore the proposed adjustment is flexible, it can be adapted to different types of data (multimodal distribution or non-exchangeable vectors, see for instance illustration in Section 7). Numerical illustrations are provided using a rainfall real data-set and using simulated data.

Some perspectives for future work are the following ones. Firstly, as we remarked above our procedure does not require any optimization procedure. However an optimization can improve the quality of the estimation. Firstly the choice of thresholds Q i and the smoothing parameters η i can be optimized (in this paper these sets are arbitrarily chosen). Also the parameters linked to the estimation of the Archimedean copula can be optimized, as the choice of a generator among its equivalence class via the point (x 0 , y 0 ), or the kernel for smoothing empirical diagonal of the copula.

Furthermore, the impact on the tail of transformed copulas has to be investigated (see for instance [START_REF] Durante | Distorted copulas: Constructions and Tail Dependence[END_REF]). In particular the relationship between the asymptote of the parametric transformation T and the regular variation of the transformed tails has been recently studied by Di [START_REF] Di Bernardino | On tail dependence coefficients of transformed multivariate Archimedean copulas[END_REF].

A good understanding of the tail behavior is indeed required to estimate the shape of the transformation near 0 and 1, in extreme quantiles where there is a lack of data. More precisely in the aforemention paper is shown that some parameters of hyperbola are linked to the upper and lower multivariate tail dependence coefficients. This implies the possibility to modify the tail dependency of the transformed distribution without changing the global adjustment. Using these results the derivation of a complete estimation procedure both for the center of the distribution and for the tails is an open interesting point.

Figure 1 :

 1 Figure1: Graphical illustration of procedure described in Proposition 2.2, for d = 1. Here an univariate data-set and the set Q = {0.1, 0.5, 0.9} are chosen. The black line represents the empirical cdf of data, F , and the red one the cdf F using the procedure in Proposition 2.2. As proved before, the quantiles of F are identified with quantiles of the obtained transformed F for each quantile level q ∈ Q (black points and associated dotted lines).

  parameters) a) Using 1) and 2), estimate diagonal δ 1 (x) and δ -1 (x), for any x ∈ [0, 1]. b) Using 3), 4) and a), get T c) Estimate of the target i-th marginal distribution, F i for i ∈ I (e.g., empirical cdf ) d) Using 5), b) and c), get T i , for i ∈ I e) Using 6) and b), get Ω(Q) f ) Using 5), 6), b) and c), get

•

  X1: Colombo (station Id: 1722 lines from 1870-01 to 2013-06, 0 lines excluded) • X2: Pamban (station Id: 433630, 1470 lines from 1891-01 to 2013-06, 3 lines excluded) • X3: Puttalam (station Id: 434240, 1734 lines from 1869-01 to 2013-06, 3 lines excluded) • X4: Thiruvananthapuram (station Id: 433710, 1926 lines from 1853-01 to 2013-06, 8 lines excluded) • X5: Trincomalee (station Id: 434180, 1734 lines from 1869-01 to 2013-06, 4 lines excluded)

Figure 2 :

 2 Figure 2: Estimates of the autocorrelation functions for the considered 5-dimensional rainfall data. Dotted horizontal lines give indicative 10% autocorrelation thresholds.

Figure 3 :

 3 Figure 3: Left: Scatter plot of ranks for the considered 797 monthly rainfall measurements (in decimeter) in 5 stations of Sri-Lanka and India between January 1893 and June 2013. Right: Geographical positions of 5 considered stations.

Figure 4 :

 4 Figure4: Errors |F (X1,X4) (x, y) -F n (x, y)|, for (x, y) in a lattice of 100 × 100 points, where F n is the empirical distribution function and F (X1,X4) is parametric model with Gamma marginals and Clayton copula (left), Frank copula (centre panel). Errors of our transformed model F are displayed in the right plot. Black cross represents the maximum error SAE 1,4 in the considered lattice (see (22)). Black dots represent the associated rainfall data (X 1 , X 4 ).

Figure 5 :

 5 Figure5: Errors |F (X2,X5) (x, y) -F n (x, y)|, for (x, y) in a lattice of 100 × 100 points, where F n is the empirical distribution function and F (X2,X5) is parametric model with Gamma marginals and Clayton copula (left), Frank copula (centre panel). Errors of our transformed model F are displayed in the right plot. Black cross represents the maximum error SAE 2,5 in the considered lattice (see (22)). Black dots represent the associated rainfall data (X 2 , X 5 ).

Figure 6 :

 6 Figure6: Errors |F (X3,X4) (x, y) -F n (x, y)|, for (x, y) in a lattice of 100 × 100 points, where F n is the empirical distribution function and F (X3,X4) is parametric model with Gamma marginals and Clayton copula (left), Frank copula (centre panel). Errors of our transformed model F are displayed in the right plot. Black cross represents the maximum error SAE 3,4 in the considered lattice (see (22)). Black dots represent the associated rainfall data (X 3 , X 4 ).

Figure 7 :

 7 Figure 7: Fit of marginals. Black full lines are the empirical marginal tail distribution functions (in log scale).The vertical lines show estimates of 90% and 95% univariate marginal quantiles. Red lines the Gamma distribution with parameters as in Table4. Green lines are ln(1 -F i ), for i = 1, 4, 5 as in the 5-dimensional transformed model with parameters as in Table2.

Figure 8 :

 8 Figure8: Left: Estimation of the 5-dimensional survival diagonal in logarithmic scale. Right: Kendall distribution function for the 5-dimensional rainfall data. Empirical diagonal and empirical Kendall (as in[START_REF] Barbe | On Kendall's process[END_REF]) are presented in black thick line; diagonal of the transformed model F and associated K φ with φ = T • exp(-t) in full red line; diagonal of parametric Gumbel model and K φ Gumbel with φ Gumbel (t) = exp(-t 1/1.33 ) in blue dashed line; diagonal of parametric Frank model and K φ F rank with φ F rank (t) = -(log(exp(-x) (exp(-2.45) -1) + 1))/2.45 in orange thick dashed line; diagonal of parametric Clayton model and K φ Clayton where φ Clayton (t) = (1 + t) -1/0.47 in green dotted line.

Figure 9 :

 9 Figure 9: Left: Correlation matrix of the considered rainfall data. Correlations greater than 60% are indicated in bold font. Right: Dendrogram resulting to the hierarchical cluster analysis on the set of dissimilarities produced by the Euclidian distance on the rainfall data. Red boxes show the two considered clusters.

Figure 10 :

 10 Figure10: Left: Adjustment on the diagonal section for the pseudo-data (Z 1 , Z 2 ); empirical estimation of the diagonal F n,(Z1, Z2) (x, x) (dotted line) versus our parametrical model (full line). Right: Plot of the pseudo-data (Z 1 , Z 2 ).

Figure 11 :

 11 Figure11: Errors |F (Xi,Xj ) (x, y) -F n (x, y)|, for (x, y) in a lattice of 100 × 100 points, where F n is the empirical distribution function and F (Xi,Xj ) is the parametric nested model in (26) for (i, j)= (1, 4) (left), (i, j)= (2, 5) (centre panel), (i, j)= (3, 4) (right). Black cross represents the maximum error SAE i,j in the considered lattice (see (22)). Black dots represent the associated rainfall data (X i , X j ).

  for τ ≡ T in (26) and d τ = 2 (left), for τ ≡ T B in (25) and d τ = 2 (centre panel), for τ ≡ T A in (24) and d τ = 3 (right).

Figure 12 :

 12 Figure 12: Admissibility conditions in (28) for transformations T in (26) (left) and T B in (25) (centre panel) and T A in (24) (right).

Figure 13 :

 13 Figure13: Left: Function m τ in (28) for a typical non-admissible external transformation τ . This non-admissible transformation τ corresponds here to a non-admissible external transformation T A in (24) however leading to a good value of global MAE and SAE criteria for the final nested model in (26), MAE = 0.0057 and SAE = 0.0654. Right: Errors |F (X3,X4) (x, y) -F n (x, y)|, for (x, y) in a lattice of 100 × 100 points, where F n is the empirical distribution function and F (X3,X4) is the parametric nested model in (26) using the non-admissible external transformation τ for the transformation T A . Black cross represents the maximum error SAE 3,4 = 0.0906 in the considered lattice. Black dots represent the associated rainfall data (X 3 , X 4 ).

Figure 16 :

 16 Figure 16: Left: Admissibility conditions in (28) for transformation T for the S 1 with n = 500. Right: Estimation of the 2-dimensional survival diagonal in logarithmic scale for S 1 with n = 500. Empirical diagonal is presented in black thick line; diagonal of the transformed model F in full red line; Gumbel model in blue dashed line; Frank model in orange dashed line; Clayton model in green dotted line; AMH model in violet dashed line; Joe model in black dashed-dotted line; Copula 4.2.12 in dashed dark red line.

  National Climatic Data Center/NESDIS/NOAA/U.S. Department of Commerce, Harvard College Observatory/Harvard University, Meteorology Department/Florida State University, and Climate Analysis Section/Climate and Global Dynamics Division/National Center for Atmospheric Research/University Corporation for Atmospheric Research (1981): World Monthly Surface Station Climatology. Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory. Dataset. http://rda.ucar.edu/datasets/ds570.0. Accessed 08 nov 2013.

Table 2 :

 2 Complete estimated vector of parameters Θ Θ Θ.

		576 0.576 -0.0566 -0.185 -1
	θ for T 1	1.509 -0.089 -0.211 0.0624 -3
	θ for T 2	0.532	0.888	0.216	0.244 -3
	θ for T 3	0.921	0.499 -0.0057 -0.083 -3
	θ for T 4	1.097	0.323	0.067	-0.001 -3
	θ for T 5	1.147	0.274	-0.102	0.116 -3

Table 3 :

 3 The bootstrapped p-values for different goodness-of-fit tests (seeGenest et al. (

	under H 0	S n	S B n	S C n	A n
	Gumbel-Hougaard 0.00331 0.00495 0.00454 0.03465
	Clayton	0.00381 0.00980 0.00704 0.00981
	Frank	0.00617 0.00941 0.00819 0.08416
	t-Student	0.00495 0.00592 0.00498 0.00963
	Normal	0.00980 0.00719 0.00454 0.00205
	Joe	0.00819 0.00495 0.00454 0.00916

Table 5 :

 5 Supremum Absolute Errors (SAE) and Mean Absolute Errors (MAE) as in (

		F	Frank Gumbel Clayton
	SAE	0.0791 0.0941	0.0995	0.1367
	MAE	0.0095 0.0138	0.0169	0.0236

Table 6 :

 6 Left: Supremum Absolute Errors SAE i,j in (22) (first lines) and Mean Absolute Errors MAE i,j in (23) (second lines) using the transformed model F 5 with parameters as in Table2. Right: Errors using parametric model with Frank copula θ = 2.450 and Gamma marginals with parameters as in Table4. Best results are indicated in bold font.

							Errors Frank copula and Gamma marginals
							X 1	X 2	X 3	X 4	X 5
		0052				X 1	0.0119 0.0053
	X 2	0.0306 0.0157 0.0088 0.0078				X 2	0.0410 0.0389 0.0195 0.0225
	X 3	0.0316 0.0677 0.0247 0.0138 0.0327 0.0076			X 3	0.0363 0.0742 0.0351 0.0147 0.0362 0.0175
	X 4	0.0328 0.0643 0.0117 0.0247	0.0536 0.0160	0.0193 0.0060		X 4	0.0396 0.0567 0.0672 0.0137 0.0180 0.0205	0.0282 0.0129
	X 5	0.0437 0.0643 0.0356 0.0155 0.0296 0.0123	0.0867 0.0318	0.0299 0.0115	X 5	0.0493 0.0838 0.0344 0.0882 0.0246 0.0174 0.0356 0.0101 0.0323 0.0055

Table 7 :

 7 Supremum Absolute Errors SAE i,j in (22) (first lines) and Mean Absolute Errors MAE i,j in (23) (second lines) using parametric model using Gumbel copula θ = 1.33 (left), and Clayton copula θ = 0.47 (right) and Gamma marginals with parameters as in Table4. Best results are indicated in bold font.

					Errors Clayton copula and Gamma marginals
					X 1	X 2	X 3	X 4	X 5
		0053		X 1	0.0119 0.0053
	X 2	0.0406 0.0389 0.0121 0.0225		X 2	0.0482 0.0158	0.0389 0.0225
	X 3	0.0401 0.0820 0.0351 0.0168 0.0373 0.0175		X 3	0.0621 0.0229	0.0959 0.0448	0.0351 0.0175
	X 4	0.0462 0.0521 0.0621 0.0282 0.0170 0.0172 0.0195 0.0129		X 4	0.0583 0.0221	0.0403 0.0429 0.0282 0.0119 0.0144 0.0129
	X 5	0.0444 0.0880 0.0352 0.0798 0.0246 0.0161 0.0379 0.0102 0.0309 0.0055	X 5	0.0456 0.0114 0.0437 0.1104	0.0534 0.0187	0.0628 0.0246 0.0237 0.0055
		Model	F	Frank Gumbel Clayton
		SSAE	0.086 0.088	0.088	0.110
		MSAE 0.011 0.012	0.012	0.012
		SMAE 0.032 0.036	0.039	0.045
		MMAE 0.004 0.005	0.005	0.005

Table 8 :

 8 Syntectic statistics associated to Tables

  The obtained results are gathered in tables below.

	Parameters F A in (24)	m	h	ρ 1	ρ 2	η
	θ A	-2.418	-0.168	-0.831	-0.517	-2
	θ A2	-0.262	0.747	-0.136	-0.135	-4
	θ A3	1.711	0.576	-0.389	-0.315	-4
	θ A5	2.061	0.300	-0.401	-0.104	-4
	Parameters F B in (25)	m	h	ρ 1	ρ 2	η
	θ B	0.168	0.776 -0.168 -0.271 -2
	θ B1	1.823	0.094 -0.277 -0.003 -4
	θ B4	-0.449 0.491	0.032	-0.029 -4

  θ 1 and θ 2 . The obtained values are gathered in table below.

	Parameters F (Z1,Z2) in (27)	m	h	ρ 1	ρ 2	η
	θ	-1.225 0.126 -0.243 -0.150 -2
	θ 1	0.401 0.682 -0.055 -0.384 -4
	θ 2	0.478 0.864 0.081 -0.250 -4

Table 9 :

 9 Supremum Absolute Errors SAE i,j (first lines) and Mean Absolute Errors MAE i,j (second lines) using the nested model F in (26) with parameters as in table above. Better results than those in Tables

	Nested model F in (26) with optimized thresholds
		X 1	X 2	X 3	X 4	X 5
	X 1	0.0088 0.0039
	X 2	0.0308 0.0146 0.0080 0.0051
	X 3	0.0305 0.0320 0.0209 0.0114 0.0066 0.0058
	X 4	0.0245 0.0393 0.0351 0.0139 0.0045 0.0100 0.0100 0.0053
	X 5	0.0342 0.0376 0.0342 0.0619 0.0237 0.0071 0.0079 0.0103 0.0204 0.0041

T , T A and T B from Proposition 2.1 in Di Bernardino and Rullière (2013b).

Table 11 :

 11 Mean and Standard deviation of the Supremum Absolute Errors (SAE) and Mean Absolute Errors (MAE) (see Equations (20)-(21) in dimension 2) for the considered 2-dimensional parametric models. Best results are indicated in bold font. Here the sample size is n =500.

		F	Frank	Gumbel	Clayton Ali-Mikhail-Haq Joe Copula Copula 4.2.12
	SAE	0.1002681 0.1174508 0.1174527 0.1174563	0.1799419	0.1174562	0.1174524
	sd(SAE) 0.0122472 0.0087554 0.0087524 0.0087497	0.0259653	0.0087466	0.0087530
	MAE 0.0122261 0.0356336 0.0386076 0.0384325	0.1166529	0.0460649	0.0367786
	sd(MAE) 0.0035945 0.0032964 0.0034439 0.0034859	0.0032843	0.0030135 0.0035753
	Models	F	Frank	Gumbel Clayton Ali-Mikhail-Haq Joe Copula Copula 4.2.12
	SAE	0.0923287 0.1167403 0.1167405 0.1167439	0.1722203	0.1167402	0.1167401
	sd(SAE) 0.0111054 0.0062205 0.0062206 0.0062198	0.0282795	0.0062203 0.0062201
	MAE 0.0101135 0.0346139 0.0378465 0.0372831	0.1154758	0.0454735	0.0359502
	sd(MAE) 0.0017111 0.0025339 0.0026389 0.0025753	0.0024937	0.0022598	0.0027549

Table 12 :

 12 Mean and Standard deviation of the Supremum Absolute Errors (SAE) and Mean Absolute Errors (MAE) (see Equations (20)-(21) in dimension 2) for the considered 2-dimensional parametric models. Best results are indicated in bold font. Here the sample size is n =1000.

Acknowledgements: The authors wish to thank the editor and the two anonymous referees whose comments helped to improve a previous version of this paper. The authors thank Gianfausto Salvadori and Fabrizio Durante for fruitful discussions. This work has been partially supported by the BNP Paribas Cardif Insurance Chair "Management de la modélisation", by the MIRACCLE-GICC project, and by ANR Research Project LOLITA.

The analytical critical layers of the distributions F B and F A are easy to obtain. For F B in (25), we have

Analytical expressions of the inverse of any transformed margins

, are available since inverse transformations are given and since the initial distribution F i is chosen to be readily invertible.

Analogously, we get, for F A in (24)

For the nested distribution F in (26), one can write,

An illustration of critical-layers ∂L A (α) and ∂L B (α) derived above is provided in Figure 14. Black dots represent rainfall data (X 1 , X 4 ) (left) and (X 2 , X 3 , X 5 ) (right).

Simulation study

In order to illustrate the replication of the good performances of the estimation procedure provide in Sections 6-7 in the case of rainfall real data, we develop in the following a simulation study.

Let (X, Y ) be a bivariate vector follows a Copula 4.2.12 in [START_REF] Nelsen | An introduction to copulas[END_REF] with θ = 3. Furthermore X ∼ α Exp(1) + (1 -α) Pareto(3) and X ∼ (1 -α) Exp(1) + α Pareto(3). In the following we consider α = 0.2 and M = 100 Monte Carlo bivariate independent samples S 1 , . . . , S M , with S j = {(X i , Y i )} i∈{1,...,n} with sample size n = 500 and n = 1000 from the distribution presented above. To illustrate the shape of the considered bivariate dependence, the pseudo-observations from the sample S 1 with n = 500 are displayed in Figure 15 (left panel).

For both marginal distributions X and Y , we tested 15 different classes of classical marginals and we fitted the best model using the Akaike Information Criterion. The best fitted marginal distributions are log-normal and the corresponding parameters for the sample S 1 with n = 500 are gathered in Table 10. The adjustment on this considered sample S 1 with n = 500 is illustrated in Figure 15 (centre and right panel).

Log

Second marginal distribution

x Figure 15: Pseudo-observations from the sample S 1 with n = 500 from Copula 4.2.12 in [START_REF] Nelsen | An introduction to copulas[END_REF] with θ = 3 (left panel). The adjustment of marginals on the sample S 1 is displayed in the centre and right panels using the log-normal distributions with parameters in Table 10 (dashed red lines). The marginal empirical distributions are represented using a black line.

As in Sections 6-7, to appreciate the quality of the adjustment we use the Supremum Absolute Error (see Equation ( 20) in dimension 2) on a bivariate lattice G and the Mean Absolute Error on the simulated data (see Equation ( 21) in dimension 2). In particular, in Tables 11-12 we gathered the mean (denoted by SAE and MAE) and the standard deviation (denoted by sd(SAE) and sd(MAE)) of these considered distances on the samples S 1 , . . . , S M for M = 100, with n = 500 (Table 11) and n = 1000 (Table 12).

In Tables 11-12 we tested different parametric models. In particular we choose a specific copula structure (see columns in Tables 1112) and the associated log-normal marginals with maximum likelihood estimated parameters (as discussed above, see Figure 15, centre and right panels) . In the first column we displayed the results for the transformed distribution function F obtained using our estimation procedure (see Section 5). In the last column we gathered the results obtained using the true copula model (Copula 4.2.12 in Nelsen (1999) with θ = 3) and with estimated log-normal marginals. Finally other copula models are illustrated in the remaining columns in order to quantify the misspecification model error.

In Figure 16 (right) we present the parametric estimation of the 2-dimensional diagonal for the sample S 1 with n = 500, using the transformed model and the classical parametric competitors introduced above. Furthermore the admissibility conditions in (28) are verified for the transformation T for the sample S 1 with n = 500 (see Figure 16, left).