Oskar Laverny
email: oskar.laverny@univ-lyon1.fr

Esterina Masiello

Véronique Maume-Deschamps

Didier Rullière

Dependence structure estimation using Copula Recursive Trees

Keywords: Bagging, CORT, Density estimation trees, Nonparametric estimation, Patchwork copula, Piecewise linear copula, Quadratic program 2010 MSC: 62E17, 62H10, 62H20, 62G30

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Although the estimation of copula [START_REF] Durante | Principles of Copula Theory[END_REF][START_REF] Joe | Dependence Modeling with Copulas[END_REF][START_REF] Nelsen | An Introduction to Copulas[END_REF] is a wide-treated subject, most efficient estimators available in the literature are based on restricted, parametric estimation. Vine copulas [START_REF] Müller | Selection of sparse vine copulas in high dimensions with the Lasso[END_REF][START_REF] Müller | Selection of Sparse Vine Copulas in Ultra High Dimensions[END_REF][START_REF] Nagler | Nonparametric Estimation in Simplified Vine Copula Models[END_REF][START_REF] Nagler | Evading the curse of dimensionality in nonparametric density estimation with simplified vine copulas[END_REF][START_REF] Okhrin | Copulae in High Dimensions: An Introduction[END_REF], although useful in high dimensions, often use parametric models, such as Archimedean copulas, as base building blocks. On the other hand, graphical models [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Li | PANDA: AdaPtive Noisy Data Augmentation for Regularization of Undirected Graphical Models[END_REF] assume a Gaussian dependence structure and therefore are fast but under restrictive assumptions. Classical nonparametric density estimators such as kernels [START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF][START_REF] Siloko | Boosting and bagging in kernel density estimation[END_REF][START_REF] Siloko | New gradient methods for bandwidth selection in bivariate kernel density estimation[END_REF][START_REF] Siloko | On asymptotic mean integrated squared error's reduction techniques in kernel density estimation[END_REF] or wavelets [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications to semi supervised learning[END_REF][START_REF] Genest | Copula Density Estimation By Wavelet Methods[END_REF][START_REF] Morettin | Wavelet-smoothed empirical copula estimators[END_REF] are not suited to satisfy constraints such as the uniformity of margins (one counter-example may be found in [START_REF] Chen | Nonparametric estimation of copula functions for dependence modelling[END_REF][START_REF] Fermanian | Weak convergence of empirical copula processes[END_REF]). We explore here a specific class of non-parametric copula density estimators with tree-structured piecewise constant densities, and design an estimator that lies in this class, the CORT estimator.

The CORT estimator is based on the density estimation tree from [START_REF] Ram | Density estimation trees[END_REF], a tree-structured non-parametric density estimator, and on the framework of patchwork copulas from [START_REF] De Amo | Extensions of subcopulas[END_REF][START_REF] Durante | Convergence results for patchwork copulas[END_REF][START_REF] Durante | Multivariate patchwork copulas: A unified approach with applications to partial comonotonicity[END_REF]. There already exist several other piecewise constant density estimators: the cascaded histograms of [START_REF] Goh | Cascaded High Dimensional Histograms: A Generative Approach to Density Estimation[END_REF], the Dirichlet-based Polya tree [START_REF] Ning | A nonparametric Bayesian approach to copula estimation[END_REF], the distribution element trees by [START_REF] Meyer | Density Estimation with Distribution Element Trees[END_REF], the adaptative sparse grids of [START_REF] Peherstorfer | Density Estimation with Adaptive Sparse Grids for Large Data Sets[END_REF], the framework of Gaussian mixtures by [START_REF] Criminisi | Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning[END_REF], the Bayesian sequential partitioning techniques by [START_REF] Li | Density Estimation via Discrepancy Based Adaptive Sequential Partition[END_REF][START_REF] Lin | Estimation of high-dimensional graphical models using regularized score matching[END_REF] with their interesting asymptotic consistency results, and the Wasserstein compression techniques provided by [START_REF] Luini | Density estimation of multivariate samples using Wasserstein distance[END_REF] are all worth noting in the field of non-parametric piecewise density estimation. But these models are built to estimate densities without taking into account uniformity of margins, and they do not always lead to proper copulas when applied on pseudo-observations or true copula samples.

The CORT estimator has the particularity of being tree-shaped which ensures on one hand that the computation of the estimated density and the distribution function on new data points is fast, and on the other hand that the storage of the model is efficient. Thus, it could be used for tasks such as re-sampling a dataset outside the already existing points, or for compression purposes, when dealing with big-data dependencies. Finally, under mild conditions, the estimator is a proper copula, where classical non-parametric estimators, such as Deheuvel's empirical copula, are not. This paper is organized as follows. Section 2 describes the class of piecewise linear copulas and gives some of their properties. In Section 3, we propose an estimation procedure, allowing localized dimension reduction, and we establish a convergence result for this procedure. Section 4 deals with ensemble models based on the CORT estimator: bagging techniques and out-of-bag generalization statistics are developed in the field of copula density estimation, and applied to the CORT estimator. Finally, Section 5 investigates the performance of the model by applications on some simulated examples, and Section 6 concludes.

The piecewise linear copula

Let X = (X 1 , . . . , X d ) be a multivariate random vector of dimension d. We are interested in the dependence structure of X. The concept of copula, whose formalization is due to [START_REF] Sklar | Fonctions de repartition à n dimension et leurs marges[END_REF], allows to study this dependence separately from the marginal distributions. Consider the distribution function (d.f.) F of the random vector X with marginal d.f.s F 1 , . . . , F d , and define the copula C as:

C(u) = F F -1 1 (u 1 ), . . . , F -1 d (u d ) , ∀u ∈ [0, 1] d ,
where F -1 i : y → inf x ∈ R : F i (x) ≥ y is the generalized inverse of F i . Then the distribution F of X is characterized by F 1 , . . . , F d and C. Sklar's Theorem [START_REF] Sklar | Fonctions de repartition à n dimension et leurs marges[END_REF] states that a copula C satisfying the previous equality always exists, and that it is uniquely defined provided that the marginal random variables are absolutely continuous. In particular, it is unique if the random vector is continuous. Note that C is the distribution function of a d-dimensional random vector with uniform margins.

The estimation of the full distribution can then be split into the estimation of one-dimensional margins, a widely treated subject, and the estimation of the copula for which we propose here to work with a piecewise linear distribution function. In the following, we define the piecewise linear copulas and then present some of their properties.

Definition

Let I = [0, 1] d be the unit hypercube, the domain of any copula, and let 1 A be the indicator function of an event A.

Definition 1 (Piecewise linear copula). Let L be a finite partition of I into subsets called leaves, usually denoted by . The piecewise linear copula with partition L and weighs p is defined by its distribution function:

∀u ∈ I, C p,L (u) = ∈L p λ (u), (1)
where λ denotes the Lebesgue measure of a set, λ (u) = λ( ) -1 λ([0, u] ∩ ), and p is a vector of non-negative weights summing to one. The corresponding copula density, which is piecewise constant, is given by:

c p,L (u) = ∈L p λ( ) 1 u∈ . (2)
This type of histogram, where leaves might not all have same shape and size, has already been used in the case of density estimation with different construction schemes and leaves shapes, e.g. with a Voronoï diagram [START_REF] Cervellera | Voronoi tree models for distribution-preserving sampling and generation[END_REF][START_REF] Geiß | Optimally solving a transportation problem using Voronoi diagrams[END_REF] or a Delaunay tessellation [START_REF] Boissonnat | Randomized incremental construction of Delaunay triangulations of nice point sets[END_REF][START_REF] Hornus | An efficient implementation of Delaunay triangulations in medium dimensions[END_REF][START_REF] Watson | Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[END_REF] as partition, or more trivially with simple sets of hyper-boxes [START_REF] Anderlini | Density Estimation Trees as fast non-parametric modelling tools[END_REF][START_REF] Lin | Estimation of high-dimensional graphical models using regularized score matching[END_REF][START_REF] Luini | Density estimation of multivariate samples using Wasserstein distance[END_REF][START_REF] Ram | Density estimation trees[END_REF].

Remark 1 (Existence). Depending on the choice of the partition L and the weights p, the distribution function C p,L is not always a copula. However, if ∀ ∈ L, p = λ( ), then C p,L is the independence copula. Therefore, for any partition, there exists at least one set of weights making the model a proper copula. On the other hand, for a partition that is too complex, this is frequently the only solution: the assumption of uniform leaves and the marginals uniformity constraints restrict the shape of bins. For polytopic leaves that are not hyper-rectangles, we do not know if weighting them efficiently can be achieved inside the copula constraints.

Before looking more precisely at the copula constraints on these piecewise constant distribution functions, we therefore restrict ourselves to the case of hyper-rectangular leaves, leading to the following definition: We call suitable a partition where every leave is a hyper-rectangle with strictly positive Lebesgue measure.

Remark 1 partly drives the definition of a suitable partition. It is also why we chose to extend the density estimation trees from [START_REF] Ram | Density estimation trees[END_REF] instead of another piecewise constant density estimator: it produces a suitable partition. If not specified, we consider by default that partitions we are dealing with are suitable. In the next subsection, we propose some properties of the dependence structure induced by such a copula.

Properties

With the above formulation of a piecewise linear copula, we can easily obtain closed-form expressions for classical quantities of interest in copula modeling. We recall some of those quantities and then derive their expression for piecewise linear copulas. The Kendall τ and Spearman ρ (see [START_REF] Nelsen | An Introduction to Copulas[END_REF]) are common dependence measures that can be computed from a copula. They are respectively defined for a copula C and its density c as:

τ = 4 C(u) c(u) du -1, ρ = 12 C(u)du -3.
Note that both τ and ρ are always in [-1, 1]. The piecewise constant expression of the density in the piecewise linear class allows for simple computation of τ and ρ, although the expressions can be somewhat cumbersome.

Proposition 1 (Common dependence measures). Let C p,L be a piecewise linear copula. Its Kendall τ and Spearman ρ are given in closed form by:

τ = -1 + 4 ∈L k∈L p p k λ( )λ(k) d i=1 (b ∧ d) 2 -(a ∨ c) 2 2 + c ((a ∨ c) -(b ∧ d)) + (d -c) (b -(a ∨ d)) , ρ = -3 + 6 ∈L p d i=1 (2 -b i -a i ) ,
where we denote = (a, b] and k = (c, d], and ∧, ∨ denote respectively the minimum and maximum operator.

The proof is postponed to the appendix. Matrices of bivariate dependence measures can be obtained by projection of the partition on couples of dimensions and using the same formula on the projected models. The closed form expression for piecewise linear copulas facilitates greatly these kinds of computations. This property will be exploited later, when introducing penalization techniques. Furthermore, these closed form expressions will be used to assess the performance of the fitting procedure that we describe in the next section.

Estimation

Suppose that we have a dataset u i, j n×d of (pseudo-)observations from an unknown copula C. We seek parameters ( p, L) of c (n) p,L , an approximation of c in the piecewise linear copula class based on these n observations. To find the optimal parameters (p * , L * ), we will adopt a two stage mechanism, considering first that the partition L is known. From now on, we denote by f 2 2 the squared L 2 norm of a function f , given by f 2 2 = f (x) 2 dx.

3.1. Optimal weights p * L knowing the partition L Suppose that a partition L is already constructed, and that we want to construct weights p to complete the approximation. As was done by [START_REF] Ram | Density estimation trees[END_REF] (see also [START_REF] Alquier | Density estimation with quadratic loss: A confidence intervals method[END_REF][START_REF] Birgé | Model selection for density estimation with L2-loss[END_REF]), we will use an Integrated Square Error (ISE) loss to build the weights. Given a copula density c, the ISE of an estimator ĉ is the squared L 2 distance to c:

I(ĉ) = ĉ -c 2 2 = (ĉ(u) -c(u)) 2 du.
To approximate the true copula density c by a piecewise linear copula, we are looking to solve:

arg min p,L I(c (n) p,L ) = arg min p,L c (n) p,L -c 2 2 = arg min p,L (c (n) p,L (u) -c(u)) 2 du = arg min p,L c (n) p,L (u) 2 -2c(u)c (n) p,L (u)du = arg min p,L c (n) p,L 2 2 -2 c (n) p,L , c , since c 2 2 is irrelevant to the minimization. Remark that c (n) p,L , c = E c (n) p,L (
U) , so that, with a slight abuse of notation, we write c (n) p,L , c even if C does not admit a density c. Furthermore, since the true copula C and its density c are unknowns, using empirical observations (u i, j ) from the copula, E c (n) p,L (U) can be approximated by n -1 n i=1 c (n) p,L (u i ). We then define our empirical loss.

Definition 3 (Empirical Integrated Square Error). Given observations (u i, j ) from a copula, define the Empirical Integrated Square Error (EISE) of an estimator ĉ of the copula density as:

Î(ĉ) = ĉ 2 2 - 2 n n i=1 ĉ(u i ). (3)
To minimize this loss, we find the weights p * L knowing the partition L by solving the following problem:

arg min p Î(c p,L ), s.t. c p,L is a copula. ( 4
)
The copula constraints in Eq. ( 4) are classically expressed in the literature as in Definition 4:

Definition 4 (Copula constraints). Let C be the distribution function of a signed measure µ C . C is a copula if the following three conditions are satisfied:

∀u ∈ I, ∀i ∈ {1, . . . , d}, C(u 1 , . . . , u i-1 , 0, u i+1 , . . . , u d ) = 0, (Ground constraint), ∀u ∈ I, ∀i ∈ {1, . . . , d}, C(1, . . . , 1, u i , 1, . . . , 1) = u i , (Marginal uniformity), ∀a ≤ b ∈ I, µ C ([a, b]) ≥ 0,
(d-increasingness). that, under these constraints, our optimization problem is in fact a quadratic program, as Proposition 2 shows.

Proposition 2 (Quadratic program). Let L be a suitable partition. Then the weights minimizing the empirical integrated square error (3) under the copula constraints from Definition 4 are the unique solution of the quadratic program:

arg min p p A L p -2 p A L f L , s.t. B L p = g L and p ≥ 0,
where the matrices A L , B L and the vectors f L , g L are given by:

A L = λ( ) -1 1 =k ∈L,k∈L , (size |L| × |L|), B 1 = λ i (u i ) (i,u)∈{1,...,d}×M L , ∈L , (size nd × |L|), B L = (B 1 , B 2 ), (size (nd + 1) × |L|), f L =        1 n n i=1 1 u i ∈        ∈L , (size |L|), B 2 = (1) ∈L , (size 1 × |L|), g 1 = (u i ) (i,u)∈{1,...,d}×M L , (size nd), g L = ( g 1 , 1), (size (nd + 1)),
where we denoted |L| the cardinal of L, and M L the set of middle-points of leaves.

Proof This proof is in three parts. First, we show that the objective formulation is correct, then we discuss the constraints formulation and finally we prove existence and uniqueness of a solution. Introduce first a new scalar product on R |L| :

x, y L =

∈L

x y λ( ) .

We denote by . 2 L its associated square norm, and by d L its associated distance. Note that the associated bilinear symmetric operator has matrix A L (defined above), a diagonal and positive definite matrix.

Using the definition of A L , f L and ., . L , the objective function in Eq. (4) rewrites:

Î(c p,L ) = p 2 L -2 p, f L L = p A L p -2 p A L f L . (5)
Then, the unconstrained version of the convex optimization problem from Eq. (4) corresponds to the projection of f L onto [0, 1] |L| (convex, closed), with respect to the norm . 2 L . Note that, if we set the weights to be equal to the empirical frequencies f L , this result yields, for the optimization of L, a loss of the form:

-f L 2 L = - ∈L f 2 λ( ) , (6)
which is the loss that was directly used by [START_REF] Ram | Density estimation trees[END_REF].

We now discuss the constraints. Denote by C L the subset of [0, 1] |L| containing vectors p such that C p,L satisfies the set of constraints from Definition 4, for a given L. The claim is that:

C L = p ∈ R |L| : B L p = g L and p ≥ 0 .
The first constraint of Definition 4 is trivially satisfied by our model. We show that the second constraint can be evaluated on only one point per leaf. Remember that the piecewise linear copula is uniformly distributed on each leaf.

Hence, for all i ∈ {1, . . . , d}, if on some point u = (u 1 , . . . , u d ) , ∈L p λ i (u i ) = u i , then defining x such that u + x is in the same leaf as u will give us that ∈L p λ i (u i + x i ) = u i + x i since only one leaf is active in the sum. Therefore, the marginal uniformity constraint may be evaluated on only one point per leaf, and hence is equivalent to:

∀ ∈ L, ∃ u ∈ , ∀i ∈ {1, . . . , d}, ∈L p λ i (u i ) = u i , where for = [a, b] we still denote i = [a i , b i ] its marginals and λ i (u i ) = λ([0,u i ]∩ i ) λ( i )
. Then, if we choose evaluation points as middle-points of leaves to put the constraints in matrix-vector form, we have the expression B 1 p = g 1 for these constraints. Furthermore, we need to force the sum of weights p to be equal to 1 (so that the total marginal mass is 1), giving a last line of ones to B L and a last value of one to g L .

The third constraint states that the measure associated with the copula C has to be positive on any hyper-rectangle

[a, b]. Recall that: ∀a ≤ b, µ C ([a, b]) = [a,b] c(u)du = ∈L p λ ([a, b])
If all weights are positive, then ∀a ≤ b, µ C ([a, b]) ≥ 0. On the other hand, if one weight is negative, then taking a, b inside the corresponding leaf would make µ C ([a, b]) negative. The last constraint is therefore equivalent to positivity of weights, which gives the last part of the wanted expression for C L . Finally, existence and uniqueness of a solution are trivial since the objective function is a quadratic function, C L is a closed and convex set, and C L is non-empty by Remark 1.

Denoting P L,S (x) the orthogonal projection of a vector x onto a set S regarding the distance d L , the quadratic program from Proposition 2 gives the optimal weights knowing the partition as:

p * L = P L,C L (f L ).
The empirical frequencies f L , which are the unconstrained solution, can then be used as a good starting point for a solver. We concentrate now on the construction of the partition L.

Locally optimal splitting point x * ,D

Suppose that we have already a suitable partition L and associated weights p L such that C p L ,L is a copula. For a given leaf ∈ L, denote L = 1 , . . . , k a partition of into k new leaves such that L 2 = L \ { } ∪ L is a new suitable partition. Then we have, as in [START_REF] Ram | Density estimation trees[END_REF], the following property: Proposition 3 (Independence of surrogate loss). Define the surrogate loss associated to the additional split from

( p * L , L) to ( p * L 2 , L 2 ) as the difference of empirical integrated squared errors Î(c p * L 2 ,L 2 ) -Î(c p * L ,L ).
Then the surrogate loss depends only on the partition L and data inside it.

Proof Index the objects of Proposition 2 by the partitions L and L 2 , and the part corresponding to L\{ } cancels.

This locality of the loss allowed [START_REF] Ram | Density estimation trees[END_REF] to use a recursive partitioning algorithm. We then only perform simple splits.

Definition 5 (Simple split and splitting dimensions). Denote x a given breakpoint in the leaf and D ⊆ {1, . . . , d} the set of splitting dimensions. Then the simple split of on x with dimensions D is defined as the partition L( , x, D) given by:

L((a, b], x, D) = × j∈ D (a j , x j ], (x j , b j ] × j∈{1,...,d}\ D (a j , b j ] .
The full partition of I obtained after the split is denoted L x, D = L \ { } ∪ L( , x, D). When D = {1, . . . , d}, we might omit it in the subscripts.

Remark 2 (Degrees of freedom). In a simple split on a set of dimension D, the weighting of the new leaves is a quadratic program with 2 |D| parameters responding to | D| + 1 linear constraints. Hence, there exists a trade-off between complexity and expressiveness of the model in the dimension | D| of the breakpoints. We will exploit this characteristic of the recursive procedure for sparsity purposes later on.

Remark 3 (No one-dimensional splits). Note that the copula constraints will not allow for estimation with only onedimensional splits (| D| = 1), as in [START_REF] Ram | Density estimation trees[END_REF], since there would be no degrees of freedom in the weights. This represents a huge problem as multivariate splits often imply bigger computational burden. Furthermore, the constraints themselves are not localized, but on the global scale, hence including them forbids a parallel implementation. We will see later that this issue can be avoided by delaying the constraint problem to a later stage of the optimization process.

Note that, neglecting the constraints, knowing D, we are going to choose the splitting point x * , D as:

x * ,D = arg min x∈ - k∈L( ,x, D) f 2 k λ(k) (7)
In the next subsection, before giving a complete description of the fitting algorithm, we talk about the localized dimension reduction procedures that are possible with these simple splits, and describe the construction of splitting dimensions D.

Optimal splitting dimensions D *

Suppose that we found x * to split a leaf in all the dimensions {1, . . . , d}. Before effectively splitting the leaf, we will check if, by chance, some dimensions can be removed from the splitting dimensions. Indeed, we will choose the splitting dimensions D based on a statistical test whose hypothesis writes:

Hypothesis 1 (Sparsity hypothesis H j ). Denoting U ∼ C, define for a given dimension j ∈ {1, . . . , d} the sparsity hypothesis as:

H j : U j ⊥ ⊥ U -j | U ∈ and U j | U ∈ ∼ U( j ),
where we denoted U j the j th marginal of the random vector U and U -j the vector of all other marginals.

The hypothesis H j literally says that, inside the leaf , the dimension j of the data is uniform and independent of the others. When H j is accepted in a leaf , we will reduce the dimension of the breakpoint x * in this leaf and in all child leaves accordingly, by removing j from the set of splitting dimensions D * . This will have the result that, inside the leaf , the final model will consider the dimension j to be uniform and independent of others.

To check this hypothesis using the integrated square error, as analyzed by [START_REF] Bowman | Density based tests for goodness-of-fit[END_REF], suppose without loss of generality that is rescaled to I, containing n observations of the random vector U ∼ F, for F the restriction of C to , rescaled to I (note that F is not a copula). This removes the conditioning in the hypothesis. Then, define the test statistic as follows.

Definition 6 (Test statistic). Denote by f (n) f,L the piecewise constant density that will be estimated on data U ∼ F, using the surrogate loss [START_REF] Cervellera | Voronoi tree models for distribution-preserving sampling and generation[END_REF], and by e j,n (x) = E f (n) f,L (x)|H j the expectation of this estimate under H j . The test statistic is given by:

I j = e j,n -f (n) f,L 2 2 ,
where L, e j,n and f (n) f,L are stochastic objects, depending on the independent and identically distributed (i.i.d.) random vectors U 1 , . . . , U n . This test statistic does not test hypothesis H j per se, but rather tests H j under the hypothesis of piecewise constant density. This is a weaker assertion, but it is enough to decide if local dimension reduction is possible or not given the current estimation stage and the data. More classical tests for independence and/or uniformity can be founded in [START_REF] Wet | Cramér-von Mises tests for independence[END_REF][START_REF] Genest | Testing for independence in arbitrary distributions[END_REF][START_REF] Hofert | A Graphical Goodness-of-Fit Test for Dependence Models in Higher Dimensions[END_REF][START_REF] Song | Testing conditional independence via Rosenblatt transforms[END_REF][START_REF] Yang | Multivariate tests of uniformity[END_REF][START_REF] Yao | Testing mutual independence in high dimension via distance covariance[END_REF].

The statistic I j has the nice property that it is constructed as a square distance, and hence is always non-negative and is 0 only under H j . On the other hand, it requires that we compute the full patchwork distribution in the two cases (under H j or not), which can be costly. Instead, we can only compute the next split, which will reduce the computation and simplify the statistic. The drawback is that the test is weakened. The next property gives an empirical form of this statistic, using this simplification.

Proposition 4 (Empirical form of the statistic). On a sample of data (u i, j ) n×d , we can approximate the statistic I j by:

I j = k∈L x * ,{1,...,d}\{ j}              f 2 k λ(k) + ∈L x * ,{1,..,d} ⊂k       f 2 λ( ) -2 f k f λ(k)                    . (8)
Proof The estimator will cut on the same breakpoints on dimensions other than j whether we work under H j or not. This gives the definition of the partitions, and hence the expression for e j,n and f (n) f,L , giving the expression of I j .

The law of the statistic (8) under H j cannot be computed explicitly. We use a Monte-Carlo simulation to compute the p-value of the test. To that purpose, simulate n uniform random variables to replace u ., j , and recompute the statistic (8), T times. Indeed, note that under the null, the values of u .,-j can be held fixed.

The full localized dimension reduction procedure is formalized in Algorithm 1. Note that the threshold α is a probability threshold for each individual test, and not a global threshold. To avoid multiple testing issues, it should be treated as a hyper-parameter and not a type 1 error. We now formalize in the next subsection the complete estimation procedure.

Full estimation procedure

Suppose that we start from the independence copula, which writes C {1},{I} in our framework. Then, if the sample of observations belonging to the sole leaf is too far from independence, i.e., if H j does not hold for all j, we construct the first breakpoint by greedily minimizing the loss (7) over the splitting point x. Rescaling the new leaves to I allows starting over and split again, until a proper stopping condition is reached: either there are no points anymore in the leaf or the leaf passes the uniformity tests. A third stopping condition is that the loss is no more reduced by splitting.

To fasten the computation, we decided to ignore the copula constraints while splitting, and enforce them at the end on the constructed partition to correct the empirical weights. Later properties of convergence back up this decision, and this furthermore allows to parallelize the splitting process. Experiments showed that the algorithm that enforces the constraints at each split is much slower (since for the optimization of the breakpoint, sub-optimization corresponding to the quadratic program of Proposition 2 must be run for each evaluation) and does not provide much better results.

More formally, Algorithm 2 below states the complete estimation procedure. The resulting estimator of the copula density, denoted by c (n) p * L ,L , is called CORT for Copula Recursive Tree. The current implementation of this algorithm as an R package cort is available on CRAN. Note that the conditions for a leaf to be splitable can vary: by default, we consider that a leaf becomes non-splitable when it contains less than two observations. Options exist in the implementation to randomize the splitting dimensions (instead of optimizing them), or to constraint the maximum number of splitting dimensions to be inferior to some threshold. We do not provide the analysis of these parametrizations. f L ,L . This additional outcome of the fitting procedure will be used in the next sections for performance analysis.

Remark 5. If we restrict the breakpoint possibility to all points with coordinates in i m+1 i∈{1,...,m} where m divides the number of observations n, this gives us only m d candidates. This corresponds to a form of checkerboard copula, see [START_REF] Cuberos | Copulas checker-type approximations: Application to quantiles estimation of sums of dependent random variables[END_REF] for more details. Since the ISE loss we use is tractable enough, the breakpoint can be chosen by directly minimizing the criterion over the continuous space I if the dimension of the problem is not too big.

After discussing the consistency of this estimation procedure, Section 4 discusses potential extensions, mainly through bagging principles.

Consistency

We will show the consistency of the CORT estimator in the L 2 almost-sure sense. Recall from [START_REF] Ram | Density estimation trees[END_REF]Theorem 1] the following result about the unconstrained estimator.

Proposition 5 (Consistency of f (n)

f L ,L ). Assuming the maximum diameter of leaves goes to 0 as n goes to ∞, we have:

Pr lim n →+∞ f (n) f L ,L -c 2 2 = 0 = 1
A detailed proof, based on a generalization by Vapnik-Chervonenkis [START_REF] Vapnik | On the uniform convergence of relative frequencies of events to their probabilities[END_REF] of the Glivenko-Cantelli Theorem, can be found in [START_REF] Ram | Density estimation trees[END_REF]. Denote now by q L the volumes given by the true copula on the leaves:

∀ ∈ L, q = c(u)du.
Then, one can easily check that q L ∈ C L , and Proposition 5 leads to the following useful corollary:

Corollary 1. d L (f L , q L ) 2 → 0, a.s.
Indeed, replacing c by a piecewise constant density with partition L and weights q L does not change the value of 2 . Definition 7 (Integrated constraint influence). Define the integrated constraint influence as the following squared norm:

f (n) f L ,L -c 2 2 and hence f (n) f L ,L -c 2 2 = d L (f L , q L )
c (n) p * L ,L -f (n) f L ,L 2 2 = d L (p * L , f L ) 2 (9)
In the simulation studies in Section 5, this quantity will be monitored via burn-in techniques, to see how it behaves as n grows. Proposition 6 below gives the corresponding theoretical result.

Proposition 6 (Asymptotic effect of constraints). As n → ∞, the integrated constraint influence goes to 0.

Proof Recall from Proposition 2 that p L * is the orthogonal projection of f L into C L . Since q L ∈ C L , we have that d L (f L , p * L ) 2 ≤ d L ( f L , q L ) 2
, which concludes the argument by Corollary 1.

The consistency of the estimator is now easy to obtain.

Proposition 7 (Consistency). For c the density of the true copula, assuming the diameter of the leaves goes to 0 as n goes to ∞, the estimator c (n) p * L ,L is consistent, i.e.:

Pr lim n →+∞ c (n) p * L ,L -c 2 2 = 0 = 1 Proof Remark that c (n) p * L ,L -c 2 2 = d L ( p * L , q L ) 2 .
Then, p * L being the orthogonal projection of f L into C L , and q L being in C L , we have d L ( p * L , q L ) 2 ≤ d L (f L , q L ) 2 and we conclude by Corollary 1. 10

Bagging of density estimators

While using the CORT algorithm, overfitting is probable since we do not use any kind of pruning. Although an implementation of pruning might be possible in our setting (remove one leaf, compute new weights through the adequate quadratic program, and compare the loss in the objective function to the final loss, then design a criterion to know where to stop), we preferred to augment the 'complexity' of the model by mixing (bagging) rather than reduce it by pruning, since the model is already quite simple. Furthermore, with a minimum node size (minimum number of points in a leaf to initiate a split) set to one, the model produces an average number of leaves growing exponentially with the dimension, a lot of them being empty. The bagging procedure we propose in this section should overcome these two issues by mixing different trees fitted on different resamples of the dataset.

Two estimates ĝ1 and ĝ2 of a function g can be bagged into an estimate ĝ1 +ĝ 2 2 , candidate for the estimation of g if the estimates are close to be uncorrelated. This principle gave rise to the bagging algorithm in regression, developed by Breiman in [START_REF] Breiman | Out-of-bag estimation[END_REF]. In density estimation, bagging can also be exploited: if the estimator has a high variance and a small bias, then bagging it might yield a better result, regarding the bias/variance trade-off.

Recall that, while bootstrapping over n observations, the chance that an observation does not appear in the bootstrap sample is given by 1 -

1 n n → 1
e . Hence, asymptotically 36.8 percent of the dataset will end up out-of-bag. These samples can be used to check for accuracy of the model, and even set hyper-parameters when some are needed. Following the work of [START_REF] Breiman | Out-of-bag estimation[END_REF] in regression models, [START_REF] Sain | Cross-validation of multivariate densities[END_REF] formalized the cross-validation and bagging process for density estimation. Note that usually, the leave-one-out method is used in kernel density estimation to select hyper-parameters (mainly the bandwidth), see, e.g., [START_REF] Miecznikowski | Bootstrap MISE estimators to obtain bandwidth for kernel density estimation[END_REF], but the more involved out-of-bag procedure we propose is inspired by [START_REF] Wu | Density estimation via the random forest method[END_REF].

In the following, we denote by Ĉ the empirical copula of the whole dataset u, a n × d matrix.

Definition 8 (CORT Forest). Define u (1) , . . . , u (N) as N bootstrap resamples of the same size n, and ĉ(1) , ..., ĉ(N) (resp. Ĉ(1) , ..., Ĉ(N) ) the densities (resp. d.f.) of the CORT estimators on these resamples fitted by Algorithm 2. Define the CORT forest with weights ω = (ω 1 , . . . , ω N ) as the mixture distribution with density:

ĉω (v) = N j=1 ω j ĉ( j) (v) ∀v ∈ I.
For each observation u i in the original training set, we recall the out-of-bag density estimate as:

ĉoob ω (u i ) = N j=1 ω j ĉ( j) (u i )1 u i u ( j) N j=1 ω j 1 u i u ( j)
.

We define Ĉω and Ĉoob ω accordingly. Note that ĉoob ω is not a proper density as it may fail to sum to 1 and it is only defined on the observation points. Because trees were fitted independently, this is however, on observation points, a good approximation of the forest density itself. From ĉoob ω , based on [START_REF] Miecznikowski | Bootstrap MISE estimators to obtain bandwidth for kernel density estimation[END_REF], [START_REF] Wu | Density estimation via the random forest method[END_REF] defines out-of-bag version of common fitting statistics as follows.

Definition 9 (Out-of-bag statistics). Define respectively the out-of-bag empirical integrated square error, the out-ofbag Kullback-Leibler divergence and two out-of-bag Cramer-Von-Mises distances associated to the forest as:

Ĵ(ĉ ω ) = ĉω 2 2 - 2 n n i=1 ĉoob ω (u i ), K(ĉ ω ) = - 1 n n i=1 ln ĉoob ω (u i ) M( Ĉω ) = 1 n n i=1 Ĉoob ω (u i ) -Ĉ(u i ) 2 , N( Ĉω ) = 1 n n i=1 Ĉω (u i ) 2 -2 Ĉoob ω (u i ) Ĉ(u i ).
Note that K(ĉ ω ) is obtained as an empirical version of ln c(u) ĉoob ω (u) dC(u). M(ĉ ω ) estimates the Cramer-Von-Mises distance between the out-of-bag d.f. of the forest and the empirical copula. On the other hand, N(ĉ ω ) keeps the true norm of the model, in the same spirit as Ĵ(ĉ ω ), see [START_REF] Miecznikowski | Bootstrap MISE estimators to obtain bandwidth for kernel density estimation[END_REF][START_REF] Wu | Density estimation via the random forest method[END_REF] for more details about these cross-validation tools.

We denote:

Ω =          ω ∈ R N , N j=0 ω j = 1 and ω j ≥ 0 ∀ j ∈ {1, . . . , N}         
the set of possible weights. Indeed, we want ĉω to be a density, which implies that weights sum to one and are positives. We call optimal the forest with the weights minimizing Ĵ(ĉ ω ). We use an optimization program to find weights given the resampling and the trees:

ω * = arg min ω∈Ω Ĵ(ĉ ω ).
Note that this method makes out-of-bag observations contribute to the final estimation through weights.

The forest estimation is studied in the next section. Note that the constructions of this section are the same if you use another base estimator than the CORT estimator: in the next section, we will use these tools to compare bagging of the CORT estimator with bagging of other copula estimators.

Investigation of performance

In this section, we investigate the performance of the proposed estimation procedure on several simulated datasets.

To ensure reproducibility, we provide all the datasets in the R package cort, with the code and parameters needed to re-simulate them. We will compare our results with several other models:

• Deheuvel's empirical copula, hereafter denoted by "Empirical";

• The empirical beta copula [START_REF] Segers | The Empirical Beta Copula[END_REF], hereafter denoted "Beta";

• A Checkerboard copula [START_REF] Cuberos | Copulas checker-type approximations: Application to quantiles estimation of sums of dependent random variables[END_REF] with m = 10, denoted by "Cb(m=10)";

• Another less-precise Checkerboard copula with m = 5, denoted by "Cb(m=5)".

Recall from [START_REF] Cuberos | Copulas checker-type approximations: Application to quantiles estimation of sums of dependent random variables[END_REF] that a checkerboard copula with a parameter m is a piecewise linear copula with a partition L composed of a regular grid of hypercubes with side length m -1 . See the reference for details about the empirical beta copula. They are all non-parametric or semi-parametric models, with a straightforward estimation procedure. The two checkerboard copulas are provided by the R package cort, and the empirical beta copula, as well as the empirical copula, are from the R package copula.

We will compare results of different models in terms of dependence measures, Kendall tau and Spearman rho first. Then, we will look at predictive performance metrics defined in the previous section, Ĵ, K, M and N, computed via a weighted bagging of each model.

Dataset 1 (Piecewise linear copula). This dataset is a simple test: we simulate random samples from a density inside the piecewise copula class, and test whether the estimator can recover it. For that, we will use a 2-dimensional sample with 500 observations, uniform on I, and apply the following function:

h 1 (u) =       u 1 , u 2 + 1 u 1 ≤ 1 4 + 21 u 1 ≤ 1 2 + 1 3 4 ≤u 1 4       .
To illustrate the behavior of the algorithm, we propose in Fig. 1 to look at a running example of the first few splits. Fig. 2 (a) shows a simulation from the final (fully fitted) CORT estimator on Dataset 1, along with results from the CORT forest. observe that all models perform correctly regarding dependence measures on this dataset, although the Checkerboard with m = 5 (which has a pretty rough partition, not including the 1 4 multiples) has a Spearman ρ a little too high. Furthermore, bagging the CORT model gives the worst results.

Performing a standard weighted bagging procedure, we obtain fit statistics K, Ĵ, M and N, displayed in Table 2. The experiment fits every model 500 times on resamples of the dataset and weights the resulting models to minimize Ĵ. Note that the Kullback-Leibler out-of-bag divergence K is infinite for the empirical copula, since it does not assign weights to points it did not see. We observe in bagging results that the predictive performance of the checkerboards is quite poor, that the forest based on empirical Beta copulas is a lot better and that the CORT forest generalized even better regarding the density-based measures Ĵ and K. Note that the bagging of empirical beta copula is a very powerful model.

Dataset 2 (Another piecewise linear copula). As for Dataset 1, we simulate from a density inside the piecewise linear copula class, by applying the function:

h 2 (u) = u 1 , u 2 2 + 1 2 1 u 1 [ 1 3 , 2 3 )
to a 200 × 2 uniform sample, and taking ranks.

This second dataset is also in the piecewise linear class, but it splits the space in a ternary way, which the recursive splitting procedure of the CORT estimator cannot reproduce. In Fig. 3 (a), you can observe simulation from the CORT copula and the CORT forest fitted on Dataset 2. Remember that the algorithm splits recursively on only one breakpoint. If the data is split in a ternary way, as in Dataset 2, it will not succeed. Looking at details from the fitting procedure, we see that the constraints forced the algorithm to return exactly the independence copula (by setting weights of each leaf equal to its volume). Indeed, while splitting, two splits in adjacent leaves are not synchronized: since the optimization routines are independent of each other, it is unlikely that they return the same breakpoint value for a given dimension, meaning that weights will not be transferable between the two zones: in our case, the constraints forced the weights back to independence. We see that the forest tried to correct this behavior, but the result is quite bad. The dependence measures (τ and ρ) are here structurally 0, and every model respected this correctly. Table 3 shows the same results as for the previous model. Note that the out-of-bag ISE did not decrease during the fit of the forest, so that more complexity was added without any gain. We also observe that many trees in the forest have zero weights, they do not fit the data enough. Finally, Table 3 shows that the predictive performance of the CORT algorithm is quite bad.

We now turn ourselves to the third dataset.

Dataset 3 (Modified Clayton). This dataset is a simulation of 200 points from a 3-dimensional Clayton copula [START_REF] Joe | Dependence Modeling with Copulas[END_REF] with θ = 7 (hence highly dependent), for the first, third and fourth marginals. The second marginal is added as independent uniform draws. Lastly, the third marginal is flipped, inducing a negative dependence structure.

Dataset 3 is based on the Clayton copula, a commonly used dependence structure in many fields of application. The estimator developed in Algorithm 2 has several options: the most important one is the inclusion, or not, of the localized dimension reduction through Algorithm 1. Since here we have a completely independent dimension, this option is worth it: it reduces by a factor of 2 the number of leaves, and hence the complexity of the model, by setting the same second edge [0, 1] to each leaf. Fig. 4 gives a representation of the tree and the statistics of the forest. The comparison of Table 4 shows that, even if the CORT algorithm performs correctly, the forest tends to be biased in the dependence measures, toward more independence. On the other hand, the predictive performance of the model from Table 5 is really high on the density-based estimates ( Ĵ and K), and is less good on the distribution function based versions ( M and N).

The last dataset was produced based on a function h 3 defined by:

h 3 (u) = u 1 , sin(2πu 1 ) - u 2 π , 1 + u 3 π 2 u 3 2 1 1 4 ≥u 1 -sin(π u 1 )1 1 4 <u 1 .
Dataset 4 (Simulated functional). We chose to produce a voluntarily hard to estimate dependence structure, by applying h 3 to uniformly drawn 3-dimensional random vectors. The dataset is the ranks of (h 3 (u i ))) i∈1,...,500 .

Fig. 5 shows the CORT estimator and the bagging statistics on this dataset. Finally, the predictive performance from Table 7 is still two-sided: the density-based results are quite good, but the distribution-function based ones are not very good.

More details and experiments are available in the supplementary section.

Conclusion

From a simple density estimation procedure designed by [START_REF] Ram | Density estimation trees[END_REF], we constructed a piecewise constant, tree-shaped, recursive copula density estimator. We computed several closed-form expressions for this estimator, and we gave an asymptotic result.

If, intuitively, constraining the space of potential weighting solutions will help the convergence of optimization routines, the copula constraints forced us to split the space in more than one dimension, making the resulting estimation procedure complex with the increasing dimension. The localized dimension reduction procedure helps to reduce the complexity.

The CORT estimator has good generalization performance and is straightforward to use since it does not have restrictive hypothesis on the true dependence structure. Although the implementation we provide is very fast, a balance between computation time and precision is available in the number of trees used in the bagging procedure. However, more work needs to be done to correct defaults of the splitting procedure, which is not able to understand certain kinds of dependence structures. On Fig. 6, note that smaller values of P (d.f. based) and Q (density based) mean a better model. We observe that, although the bagging procedure is not worth it, the CORT estimator is very good on this example, both for density and d.f. estimation. Unfortunately, it is not always the case, as shown by the experiments we did on Dataset 2. The box plot of P error on Dataset 2 confirms the analysis we already had: in this case, the empirical beta copula performs a lot better.

For Dataset 3, however, the box plot of P and Q in Fig. 8 confirms the performance of the estimator: To understand more precisely what happened with Kendall's tau and Spearman's rho on Dataset 3, we ran a burn-in experiment: we fitted trees on subsamples of increasing size. We can then observe the burn in Kendall taus and Spearman rhos, represented on Figure 9. We see that the fitted values, on the right of Fig. 9, are convergent but biased compared to the empirical observed values of the dependence measures directly computed on resamples of the dataset, on the left. We also observe the high variance of the estimator, which is one of the good reasons to use a bagging procedure.

Finally, the box plots on Dataset 4 are given in Fig. 10.

Definition 2 (

 2 Hyper-rectangles and suitable partitions). Let a and b both be in I. Then, if a ≤ b (component-wise), we define the hyper-rectangle (a, b] as (a, b] = (a 1 , b 1 ] × ... × (a d , b d ].

Algorithm 1 :

 1 Localized dimension reduction. Data: T ∈ N, a leaf , observations u 1 , .., u n ∈ , and a threshold probability α Result: The splitting dimensions D * 1 Obtain x * and L x * ,{1,..,d} by greedily minimizing the loss in (7) on . 2 foreach j ∈ 1, . . . , d do 3Denote by s the statistic I j given by[START_REF] Chen | Nonparametric estimation of copula functions for dependence modelling[END_REF].

4 foreach i ∈ 1 , 6

 416 . . . , T do 5 Simulate n uniform random variables on j and replace the jth coordinate of the data by this simulation; Denote s i the value of the statistic from (8) on this new dataset; return D * = j ∈ 1, . . . , d : p j > 1 -α

Algorithm 2 : 3 foreach ∈ L that is splitable do 4 Run 5 if D * ∅ then 6 7 Set 11

 23456711 CORT estimation. Data: Observed ranks u 1 , .., u n ∈ I Result: Parameters p and L of the estimated piecewise linear copula 1 Initialize the tree by L = {I} and p L = {1}. 2 while there exist leaves ∈ L that are still splitable do Algorithm 1 in to find D * . Find x * , D * minimizing the surrogate loss (7). L = L x * ,D * , D * . Compute p * L = P L,C L (f L ) from Proposition 2. 12 return ( p * L , L).

Remark 4 .

 4 Before step 11 of Algorithm 2, the built model is simply the density estimate f(n)

Fig. 1 :Fig. 2 :

 12 Fig.1:(Dataset 1) Running example. Data points are in red, and darker zones mean boxes with higher weights (black is the maximum possible in each leaf). On the left, the model after the first split. In the middle, the model after the next round of splits, and on the right after the third round.

Fig. 3 :

 3 Fig. 3: (Dataset 2) (a) The estimated tree: in black, lower left, the input data. In red, upper-right, a simulation from the estimated tree. (b) The estimated forest: in black, lower left, the input data. In red, upper-right, a simulation from the estimated forest.

Fig. 4 :

 4 Fig. 4: (Dataset 3) (a) Representation from the tree: in black, lower left, the input data. In red, upper-right, a simulation from the estimated tree. (b) Forest Statistics: on the left, K and Ĵ in function of the number of trees. On the right, the Integrated Constraint Influence and square norm of each tree against the weight of the tree in the forest.

Fig. 5 :

 5 Fig. 5: (Dataset 4) (a) Representation from the tree: in black, lower left, the input data. In red, upper-right, a simulation from the estimated tree. (b) Forest Statistics: On the left, K and Ĵ in function of the number of trees. On the right, the Integrated Constraint Influence and square norm of each tree against the weight of the tree in the forest.

( a )Fig. 6 :

 a6 Fig. 6: (Dataset 1) Box plots of resulting errors for 20 resamples for each model (the lower the better). P is focused on the distribution function and Q on the density.

( a )Fig. 7 :

 a7 Fig. 7: (Dataset 2) Box plots of resulting errors for 20 resamples for each model (the lower the better). P is focused on the distribution function and Q on the density.

( a )Fig. 8 :

 a8 Fig. 8: (Dataset 3) Box plots of resulting errors for 20 resamples for each model (the lower the better). P is focused on the distribution function and Q on the density.

Fig. 9 :

 9 Fig. 9: (Dataset 3) Burn in of dependence measures: on the top row, Kendall's tau; on the bottom row, Spearman's rho; on the left column, empirical values from subsamples of increasing sizes; on the right column, values obtained by the fitted CORT estimators on the same subsamples. The size of the subsamples is in abscissa. Each line type corresponds to a couple of variables.

( a )Fig. 10 : 2 - 2 -

 a1022 Fig. 10: (Dataset 4) Box plots of resulting errors for 20 resamples for each model (the lower the better). P is focused on the distribution function and Q on the density.

Table 2 :

 2 Results of the bagging of each model on Dataset 1. Each row represents a different performance metric: in all cases, lower is better.

		Empirical	Cb(m=10)	Cb(m=5)	Beta	CORT
	Ĵ(ĉ ω )	0.002	-3.16	-2.37	-2.98	-4.81
	K(ĉ ω )	Inf	-1.06	-0.743	-1.18	-0.837
	M(ĉ ω )	8.72e-06	4.13e-05	0.000281	1.58e-05	0.000633
	N(ĉ ω )	-0.0277	-0.0277	-0.0275	-0.0277	-0.0262

Table 3 :

 3 Results of the bagging of each model on Dataset 2. Each row represents a different performance metric: in all cases, lower is better.

		Empirical	Cb(m=10)	Cb(m=5)	Beta	CORT
	Ĵ(ĉ ω )	0.00501	-1.48	-1.39	-1.21	-2.55
	K(ĉ ω )	Inf	-0.346	-0.271	-0.426	-0.311
	M(ĉ ω )	3.03e-05	5.16e-05	0.00016	7.38e-05	0.000907
	N(ĉ ω )	-0.134	-0.134	-0.134	-0.134	-0.134

Table 5 :

 5 Results of the bagging of each model on Dataset 3. Each row represents a different performance metric: in all cases, lower is better.

		Empirical	Cb(m=10)	Cb(m=5)	Beta	CORT
	Ĵ(ĉ ω )	0.00501	-9.27	-7.69	120	-54.6
	K(ĉ ω )	Inf	Inf	Inf	-0.582	-1.97
	M(ĉ ω )	2.84e-05	2.51e-05	5.25e-05	3.38e-05	9.41e-05
	N(ĉ ω )	-0.000666	-0.000669	-0.000641	-0.000657	-0.000639

Acknowledgments

We are particularly grateful to the two referees, the associate editor and the editor who all made useful comments and suggestions, improving the manuscript. We would like to thank SCOR SE for funding this work through a CIFRE grant. Any errors are ours.

We observe on Fig. 2 (a) that the algorithm splits the space as requested. The few simulated points outside the four main boxes are there because the algorithm did not split exactly on ( 14 , 3 4 ), ( 1 2 , 1 2 ) and ( 3 4 , 1 4 ), and the constraints forced him to put some weight on some leaves that do not contain points.

Bagging the CORT algorithm on this dataset, we obtain statistics given by Fig. 2 (b). We observe that the out-of-bag statistics are decreasing in the number of trees fitted, although the weighting of the forest did select less than 10 trees over 500. Altogether, the algorithm succeeded into finding the right breakpoints. A comparison of the fit in terms of dependence measure to other models is available in Table 1. 1, we display pairwise dependence measures (Kendall τ and Spearman ρ) of the obtained fits. To read these measures, consider the first column, corresponding to the empirical copula, as the goal for other models. We

On the left of Fig. 4 (b), the convex, decreasing shape of the Kullback-Leibler divergence with respect to the number of trees shows that the generalization error of the forest decreases with the number of trees. The decreasing trend of the constraint influence and the square norm of trees with respect to the assigned weights by the forest, on the right of Fig. 4 (b) shows how the weighting procedure selected trees. Table 4 shows dependence measures obtained from the different models. Although the estimation on the second and third marginals is not very good, the statistics of the forests are good, showing that the estimator get better and better while adding trees. Indeed, on Fig. 5 (b), left, we observe decreasing out-of-bag errors, but we also observe on Fig. 5 (b), right, a very skewed and fat-tailed density for the constraint influence, meaning that certain trees did not produce very good partitions. 6 shows that the bivariate projections were not all treated as well as others by the CORT algorithm: the values of τ 1,2 and ρ 1,2 are surprisingly quite good compared to τ 1,3 , τ 2,3 , ρ 1,3 , ρ 2,3 , for the CORT estimator. Hopefully, the bagging corrects this bias quite correctly.

Supplementary material

On each dataset from Section 5, we ran an additional experiment. To observe the predictive performance of each model, we designed a cross-validation procedure: on 20 resamples of each dataset, we computed the Cramer-Von-Mises and ISE errors on test samples, given respectively by:

where T is the test dataset used to obtain D, a given copula estimator with density d. Note that P is focused on the distribution function and Q on the density. Below are the resulting box plots for each of the datasets.

Fig. 6 gives a box plot of P Cramer-Von-Mises errors on the first dataset.