## 1 Introduction

Copula models are useful tools for the analysis of multivariate data since by using the well-known Sklar’s theorem, any multivariate joint distribution can be decomposed into its univariate marginal distributions and a copula function, which allows for capturing the arbitrary dependence structure between several random variables. As a result, copulas have been widely used in the field of finance, insurance, system reliability, etc. among many other application areas. See, e.g.,

Jaworski et al. (2010), Joe (2014) and Nelsen (2007) for more details about copulas and their applications.Given a random vector

with joint cumulative distribution function (CDF)

and continuous marginal CDFs , by Sklar’s theorem (Sklar (1959)), the CDF can be expressed uniquely as , where denotes the copula function. A copula is itself the joint CDF of a random vectorhaving its marginals as uniform distributions on

, henceforth denoted by . It is to be noted that the original results in Sklar (1959) are also applicable to discrete-valued random variables. However, the focus of this paper is modeling continuous-valued multivariate random vectors. Thus, for the rest of the paper, we assume the marginal CDFs are absolutely continuous.As copula plays an important role in capturing the general dependence structure between multiple variables, it is critical to estimate copula in an accurate way, especially in higher dimensions where the dependence structure becomes much more complicated, often illusive and may even be supported on a lower-dimensional manifold. One of the primary goals of this paper is to estimate a smooth copula function from a random sample of independent identically distributed (iid) observations for .

Many parametric families have been proposed for modeling multivariate copulas and there has been previous work addressing the corresponding parametric estimation methods. For detailed discussions, see, e.g., Joe (2005), Joe (2014), McNeil et al. (2009), Nelsen (2007), Smith (2011) and Žežula (2009)

, etc. However, no matter how sophisticated and flexible the parametric models that we may use, they might still lead to biased copula estimates when the parametric model is misspecified and thus may not be able to capture complex dependence structures required in practice. Compared to standard multivariate copulas, vine copula models allow for more flexibility in capturing complex dependency structures using appropriate vine tree structures by choosing bivariate copula families for each node of pair copulas from a vast array of parametric bivariate copulas. But it is often challenging to obtain estimates of multivariate dependence measures which involve high-dimensional integrals that are often algebraically intractable by using vine copulas.

Thus, recognizing some of the above-mentioned limitations of parametric copula models, a variety of nonparametric estimators have been proposed for multivariate copula estimation. Most of the available nonparametric estimators rely on the empirical methods, e.g., the empirical copula and its multilinear extension, the empirical multilinear copula (Deheuvels (1979); Fermanian et al. (2004); Genest et al. (2017)), or kernel-based methods such as local linear estimator (Chen and Huang (2007)), mirror reflection estimator (Gijbels and Mielniczuk (1990)) and improvements of these two estimators (Omelka et al. (2009)). See alsoRémillard and Scaillet (2009) and Scaillet and Fermanian (2002)

for other nonparametric copula estimators. However, except for the empirical multilinear copula, most of these estimators are valid copulas only asymptotically, meaning that they are not necessarily genuine copulas for finite samples. Moreover, multivariate dependent measures (e.g., Spearman’s rho, Kendall’s tau, etc.) based on such estimated copula could take values outside of their natural range, thus making them unattractive in practice. On the other hand, there has been recent work on Bayesian nonparametric methods for estimating general d-dimensional copula and among many others, a noteworthy Bayesian nonparametric model is based on an infinite mixture of multivariate Gaussian or the skew-normal copulas proposed by

Wu et al. (2014). The infinite mixture models provide a lot of flexibility in modeling various dependence structures but those typically lack simple (analytic) expressions of dependence measures making them harder to compute in practice.The primary focus of this paper is the nonparametric estimation of multivariate copulas for any arbitrary dimensions that are genuine copulas for any finite sample size and are uniformly consistent as the sample size gets large. We consider an extension of the Bernstein copula (Sancetta and Satchell (2004)), which is a family of copulas defined in terms of multivariate Bernstein polynomials. One of the primary advantages of the Bernstein copula is that it provides a class of nested models that are able to uniformly approximate any multivariate copula with minimal regularity conditions. A simple case of the Bernstein copula is the empirical Bernstein copula, which is a nonparametric copula estimator proposed by Sancetta and Satchell (2004). The asymptotic properties of the empirical Bernstein copula are well studied in Janssen et al. (2012) and its application in testing independence is described in Belalia et al. (2017). The application of the Bernstein copula to the modeling of dependence structures of non-life insurance risks is provided in Diers et al. (2012), among many other applications.

However, the empirical Bernstein copula has two main drawbacks that could prevent us from obtaining accurate copula estimation for small samples: (i) the empirical Bernstein copula is not necessarily a valid copula itself, which is a common disadvantage for most nonparametric copula estimators; and (ii) the degrees of Bernstein polynomials are often set to be equal to an integer across different dimensions, which limits the flexibility of the Bernstein copula and thus might not be appropriate for large dimensions.

In order to address the above-described problem (i), Segers et al. (2017)

show that the empirical Bernstein copula is a genuine copula if and only if all the polynomial degrees are divisors of the sample size, and further propose a new copula estimator called the empirical beta copula, which can be seen as a special case of the empirical Bernstein copula when the degrees of Bernstein polynomials are all set equal to the sample size. The empirical beta copula is a valid copula itself and has been shown to outperform some classical copula estimators in terms of bias and variance. But it always has a larger variance compared to the empirical Bernstein copula with smaller polynomial degrees. It is surprising that much less attention has been given to the problem (ii), and even for equally set degrees, there has been limited work on the data-dependent choice of degrees in the literature.

Janssen et al. (2012) recommend an optimal choice of the equal degrees in the bivariate case by minimizing the asymptotic mean squared error. Nevertheless, such a choice requires the knowledge of the first- and second-order partial derivatives, which might not be easy to estimate in practice. Burda and Prokhorov (2014) put priors on the polynomial degrees, however, their priors don’t rely on data or sample size and they use multivariate Bernstein density instead of Bernstein copula density. The Dirichlet process assigned as the prior for the copula doesn’t guarantee the copula estimate to be a valid copula itself. Besides, the number of weights grows exponentially as the dimension increases, leading to computational inefficiency of MCMC methods for larger dimensions. To the best of our knowledge, Lu and Ghosh (2020) first develop a data-dependent grid search algorithm for the selection of polynomial degrees which has shown superior empirical estimation properties for small to moderate sized samples. But the methodology is limited to bivariate cases and extension to larger dimensions remained challenging.For the purpose of addressing the two problems as described above, we introduce a new nonparametric smooth estimator for multivariate copula that we call the empirical checkerboard Bernstein copula (ECBC), which is constructed by extending the Bernstein copula allowing for varying degrees of the polynomials. It is shown to be a genuine smooth copula for any number of degrees and any finite sample size. Furthermore, we develop an empirical Bayesian method that takes the data into account to automatically choose the degrees of the proposed estimator using its posterior distribution thereby accounting for the uncertainty of such tuning parameter selection. As shown in Segers et al. (2017), larger degrees of the Bernstein copula lead to a larger variance of the estimation, so a choice of degrees that is relatively small compared to the sample size but sufficient for a good copula estimation is desirable. The degrees are allowed to be dimension-varying within the Bayesian model, which provides much more flexibility and accuracy, especially in higher dimensions.

It is especially noteworthy that while the focus of the paper is to estimate the copula function, it is straightforward to get a closed-form estimate of the corresponding copula density by taking derivatives of the ECBC. However, direct estimation of a closed-form copula function has many advantages compared to first estimating a copula density, e.g., it is often easier to differentiate than to integrate for higher dimensions. Besides, for those copulas which are not absolutely continuous such as Marshall-Olkin copulas (Embrechts et al. (2001)) having support on a possibly lower-dimensional manifold, the direct estimation of the copula density could be difficult. Owing to the closed form of the estimated copula function and its density, it can be shown that the proposed ECBC allows for straightforward estimation of various dependence measures.

The rest of the paper is organized as follows: in Section 2 we present an empirical Bayes nonparametric copula model. In Section 2.1, we derive the closed-form expression of estimates of popular multivariate dependence measures based on the novel methodology of multivariate copula estimation. We then illustrate the performance of the proposed methodology in Section 3. 3.1 shows the finite-sample performance for bivariate cases. The accuracy of the estimation of multivariate dependence measures is investigated in Section 3.2. Section 3.3 illustrates the estimation of tuning parameters of the proposed ECBC copula estimator and the comparison with the empirical Bernstein copulas is provided in Section 3.4. Section 4 provides an application to portfolio risk management. Finally, we make some general comments in Section 5.

## 2 An Empirical Bayes Nonparametric Copula Model

Suppose we have i.i.d. samples , where is a cumulative distribution function and is the absolutely continuous marginal CDF of the -th component. By Sklar ’s theorem (Sklar (1959)), there exists a unique copula such that

and

The Bernstein copula is a family of copulas defined in terms of Bernstein polynomials and it was first introduced by Sancetta and Satchell (2004). It is a flexible model that can be used to uniformly approximate any copula. The Bernstein polynomial with degrees of a function is defined as

(1) |

and is called the Bernstein copula when C is a copula.

A general estimation for the Bernstein copula of an unknown copula C is the empirical Bernstein copula (Sancetta and Satchell (2004)) , where is the rank-based empirical copula. We denote the empirical Bernstein copula as

where

where denotes indicator function and slightly modified empirical marginal distribution functions are defined as

where the modification instead of modifies the standard empirical marginal distribution to be away from 1 in order to reduce potential problems at boundaries.

However, the empirical Bernstein copula is not guaranteed to be a valid copula for finite samples as the empirical copula is not necessarily a genuine copula. Segers et al. (2017) show that the empirical Bernstein copula is a copula if and only if all the degrees are divisors of . In order to obtain a valid copula estimation for any degrees, we replace the empirical copula with the empirical checkerboard copula , which is a simple multilinear extension of the empirical copula defined as

where is the rank of among , see, e.g., Carley and Taylor (2002) and Li et al. (1997) for more details. Notice that the main difference between the empirical copula and the empirical checkerboard copula is that is a genuine copula, so we can obtain a valid copula estimation taking the form

(2) |

where

(3) |

and we call the proposed empirical checkerboard Bernstein copula (ECBC).

Unlike the empirical Bernstein copula, the ECBC is a genuine copula for any degrees and any fixed sample size . It is known that Bernstein polynomials with smaller values of degrees ’s may lead to biased estimates while unnecessary larger degrees of Bernstein polynomials will necessarily lead to larger variances. So it is critical to choose the proper degrees of the ECBC based on a given sample. In order to do that, we develop an empirical Bayes method for choosing ‘optimal’ degrees , where ’s are allowed to be different for different and also depend on the random sample of observations.

As illustrated in Sancetta and Satchell (2004), using partial derivatives of (2) with respect to each and rearranging, we can obtain the density corresponding to ECBC as follows:

(4) | ||||

where

(5) |

Clearly, the Bernstein copula is a mixture of independent Beta distributions leading to a tensor product form. For notational convenience, let us denote

Following the work by Gijbels et al. (2010), the pseudo-observations can be treated as samples from . We then use this approximation to build an empirical Bayesian hierarchical model:

where is the ‘floor’ function denoting the largest integer not exceeding the value and

, i.e., are samples from the empirical checkerboard copula . It then follows that

Based on the proposition 1 in Genest et al. (2017), can be drawn using the following hierarchical scheme:

where and denotes the discrete uniform distribution, i.e., for . Assuming that there are no ties in the pseudo samples (owing to absolute continuity of marginal distributions or breaking it by random assignment in practice), we can equivalently represent the ’s more conveniently as

Next, to account for the uncertainty in the estimation of the degrees ’s, we propose to introduce a sample-size dependent empirical prior distribution on the degrees

and obtain posterior estimates by Markov Chain Monte Carlo (MCMC) methods. This would not only allow for the almost automatic adaptive estimation of the degrees (based on the observed data) but also allow for quantifying the uncertainty of this crucial tuning parameter vector. Notice that the idea of putting priors on the polynomial degrees has also been adopted by

Burda and Prokhorov (2014). However, their priors don’t rely on data or sample size and they use multivariate Bernstein density instead of Bernstein copula density, i.e., the weights are belonged to a simplex without any more constraints. A Dirichlet process with a baseline of uniform distribution on is assigned as the prior for the copula in (1), which doesn’t guarantee to be a valid copula. In order to avoid the construction of priors under constraints, we use the empirical estimates for the coefficients of the Bernstein copula instead of assigning priors to them.Motivated by the asymptotic theory of the empirical Bernstein estimator, e.g., as in Janssen et al. (2014)

, we propose the hierarchical shifted Poisson distributions as the prior distribution for

:(6) |

and

(7) |

The following theorem provides the large-sample consistency of the ECBC using the same set of assumptions as required for the large-sample consistency of the empirical checkerboard copula.

###### Theorem 1.

###### Proof:.

We denote the ECBC as for simplicity. Also, the empirical Bernstein copula and the Bernstein copula are denoted as and , respectively. Let denote the supremum norm of a function defined on d-dimensional square . Using the familiar triangle inequality we have

First, under the assumption that the marginal CDFs are continuous, it follows from the Remark 2 in Genest et al. (2017) that

Next, notice that

In above the second inequality follows from the fact that since

are binomial probabilities,

for any and for any .Next by using the Lemma 1 in Janssen et al. (2012) and equation (3) in Kiriliouk et al. (2019), we obtain

Hence, it now follows that

Also, by using Lemma 3.2 in Segers et al. (2017) we have

Thus, combining the above inequalities that are satisfied almost surely (a.s.) for every fixed ’s we obtain

Next we consider the proposed empirical priors on the degrees to be

We make use of the following simple Lemma:

###### Lemma 1.

Suppose , then

Proof of Lemma 1: By Jensen’s inequality for the square-root function, it follows that

.

Notice that as , and conditioning on , we then have by the above Lemma, as . Thus, taking expectation with respect to the prior distribution, it follows that

∎

Notice that in the above result the a.s. convergence is with respect to the empirical marginal distribution of the data integrating out the conditional empirical distribution of data (given the ’s) weighted by the empirical prior distribution of the tuning parameters ’s. This is not the usual notion of posterior consistency but rather the notion can be viewed as using integrated likelihood approach (Berger et al. (1999)) with respect to the empirical marginal distribution obtained by integrating the priors given by equations (6) and (7).

It is to be noted that the joint posterior distribution of () may not preserve necessarily an exchangeable structure as the above prior. Using the empirical Bayes hierarchical structure of the above-proposed model it can be shown that efficient MCMC methods can be utilized to draw approximate samples from the path of a geometrically ergodic Markov Chain with posterior distribution as its stationary distribution. By generating a sufficiently large number of MCMC samples we can estimate the marginal posterior mode of the discrete-valued parameter ’s as final estimates. Let denote MCMC samples of and for each , let denote the distinct values among these MCMC samples. Then the (marginal) posterior mode of is estimated by

The final estimate of the smooth copula based on the proposed ECBC is then given by

where

It is to be noted that other posterior estimates (e.g., posterior mean when it exists or coordinate-wise posterior median or some version of multivariate posterior median) can also be used but for simplicity (and the requirement that these posterior estimates of ’s be necessarily integer-valued) we have chosen to use posterior mode based on the marginal posterior distributions of ’s. Through many numerical illustrations we show the easy applicability of this choice in various examples.

### 2.1 Multivariate Dependence Estimation

In higher dimensions, it is often of interest to evaluate the strength of dependence among variables. This is often done using copulas since most dependence measures can be expressed as a function of copulas. Spearman’s rank correlation coefficient (Spearman’s rho) is one of the most widely used dependence measures. For a bivariate copula , Spearman’s rho can be written as

A multivariate extension of Spearman’s rho given in Nelsen (1996) takes the form

(8) |

Compared to vine copulas that rely on pair copulas and complex tree structures, one of the advantages of our copula estimator is that it is straightforward to obtain an estimate of multivariate Spearman’s rho as

(9) |

where is the beta function.

It can be shown that the multivariate Spearman’s rho is bounded by

where the lower bound approaches to zero as dimension increases. Since our copula estimator is a genuine copula, the estimate of multivariate d-dimensional Spearman’s rho can avoid taking values out of the parameter space, which might be an issue for estimates built on other nonparametric copula estimators, e.g., the empirical copula (see Pérez and Prieto-Alaiz (2016)).

Similar to Spearman’s rho, Kendall’s tau is another common dependence measure and has its multivariate version as well, which is given by Nelsen (1996) as

(10) |

By applying (4) and (LABEL:eq:density2), it is also easy to obtain an estimate of multivariate Kendall’s tau based on our copula estimator as

Thus, using our proposed ECBC copula not only are we able to obtain a fully non-parametric estimate of any copula function in closed form (once the tuning parameters are estimated by their posterior modes) but also we are able to derive the closed-form expression of estimates of the popular multivariate measures of dependence for any arbitrary dimension .

Moreover, although we only illustrate the use of multivariate extensions of Kendall’s tau and Spearman’s rho as possible measures of multivariate dependence, any other multivariate notion of dependence measures that are suitable functionals of the underlying copula can also be computed using our closed-form expression of the ECBC estimator. This is particularly advantageous compared to even some of the flexible yet complicated parametric copula family (e.g., Archimedian, multivariate Gaussian or t, etc.) for which it may require high-dimensional numerical integration to compute multivariate versions of Spearman’s rho as given in (8) and/or Kendall’s tau given in (10). For vine copulas, it is particularly challenging to obtain estimates of these multivariate measures of dependence as such high-dimensional integrals are often algebraically and even numerically intractable say for dimension , whereas for ECBC even when a new measure of dependence is created as a functional of the copula that may be more complicated than those defined in eq.s (8) and (10), we can easily obtain a large number of Monte Carlo (MC) samples from the ECBC and use MC based approximation to estimate such new measures of multivariate dependence (we illustrate such a case in our real case study involving portfolio risk optimization in Section 4).

## 3 Numerical Illustrations using Simulated Data

### 3.1 Finite-Sample Performance for Bivariate Cases

We investigate the finite-sample performance of the ECBC through a Monte Carlo simulation study. Samples from the true copula are generated using the package in (Hofert et al. (2014)). In order to visualize the results using contour plots, we first restrict our illustration to bivariate copulas. Three copula families with various parameters and an asymmetric copula are considered. The first four examples are the bivariate Frank copulas

with parameter equal to and , which reflects a wide range of dependence from negative to positive. The next two examples are the Clayton copula with parameter

and the Gumbel copula with parameter

The value of Kendall’s tau is 0.33 for Clayton copula with parameter and 0.5 for Gumbel copula with parameter . Both cases have a moderate positive dependence.

The next example is the independent copula . Finally, we consider an asymmetric copula

In the simulation study, samples are drawn from the true copula for each replicate (of size ) and there are replicates. Degrees of the ECBC are estimated by posterior modes by obtaining 5000 MCMC samples following 2000 burn-in samples of two chains for each of the replicated data sets generated from a chosen true copula model. It is to be noted that it takes about 2, 11, and 60 minutes to run 7000 iterations for two chains by using MacOS with 16GB of RAM for , , and , respectively, when data are drawn from multivariate Frank copula. Convergence of MCMC runs was monitored based on preliminary runs using standard diagnostics available in R packages rjags and coda. We show the results for eight copulas in Fig. 1 and 2. The contour plot of the true underlying copula and the empirical MC average of copula estimates are given for comparison and represents the MC mean of the posterior modes of degree parameters.

We can see from the contour plots that the average of the estimated copula is extremely close to the underlying true copula across all different dependence structures irrespective of the assumed parametric models. This illustrates that the proposed ECBC has a robust performance in estimating various true copulas.

### 3.2 Accuracy of Multivariate Dependence Measures ()

To assess the finite sample performance of the estimate of multivariate Spearman’s rho , we conduct Monte Carlo simulations for copulas. We consider the independent copula and Clayton copulas with parameter value , respectively. For each copula model, Monte Carlo replicates are generated with size . For each replicate, we compute the proposed estimator in (9) and the estimator based on empirical copula

(11) | ||||

Finally, for each estimator we compute the mean, bias, variance and mean square error (MSE) over all replicates.

Table 1 shows the results on the estimation of multivariate Spearman’s rho for four different copulas. An approximated value of the true multivariate Spearman’s rho can be obtained by numerical integration since there is no analytical expression as a function of the parameter (see Pérez and Prieto-Alaiz (2016)), which is also given in Table 1. The corresponding boxplots for the two estimates based on replicates along with a horizontal line for true multivariate Spearman’s rho are shown in Fig. 3.

From the results we can see our estimator outperforms with respect to variance and MSE. In terms of bias, tends to underestimate and have a larger bias as strength of dependence increases, but there is not a clear superiority of one estimator over the other. As shown in Fig. 3 (c) and (d) where there is a moderate or strong dependence in trivariate cases, can take values out of parameter space ( and of are taking values greater than in (c) and (d), respectively), which can be problematic in measuring dependence.

Copula | Estimator | Mean | Bias | Variance | MSE | |
---|---|---|---|---|---|---|

Independent | 0 | 0.007 | 0.007 | 0.008 | 0.009 | |

-0.015 | -0.015 | 0.014 | 0.015 | |||

Clayton | 0.308 | 0.294 | -0.014 | 0.007 | 0.008 | |

0.287 | -0.021 | 0.028 | 0.029 | |||

Clayton | 0.504 | 0.477 | -0.027 | 0.007 | 0.008 | |

0.519 | 0.015 | 0.039 | 0.040 | |||

Clayton | 0.717 | 0.680 | -0.037 | 0.004 | 0.006 | |

0.732 | 0.015 | 0.050 | 0.051 | |||

### 3.3 Estimation of Tuning parameters of ECBC

We now illustrate one of the primary advantages of our proposed empirical Bayes estimate of the ECBC that allows for data-dependent automatic selection of dimension-varying degree parameters ’s. We further explore the special case of choosing equal degrees by using the following prior distribution

(12) | |||

(13) |

For our empirical illustration, we consider three true bivariate copulas and one trivariate copula to explore the comparative performance of choosing dimension-varying degrees compared to setting them equal across all dimensions. The first example is the Farlie-Gumbel-Morgenstern (FGM) copula with parameter . The next two choices are the independent copula (e.g., FGM with ) and the Gaussian copula with positive dependence (correlation

). The last one is the trivariate t-copula with degrees of freedom

and pairwise dependence and . Again samples of size are obtained for each four cases and repeated times for MC evaluation. For each sample, we choose the degrees of the ECBC by computing the posterior modes using our proposed empirical Bayesian method.Fig. 4 presents the scatterplot of estimated values of or for each chosen true copula model. From the plots we can observe that for bivariate copulas, posterior estimates of and

are significantly different in most cases without any prior restrictions of equality. In fact, posterior probability of choosing equal

is only about and for FGM copula, for independent copula, and for Gaussian copula, respectively indicating against forcing as is popularly done in the literature. For the trivariate t-copula case, the posterior probability of is 0, decisively suggesting that equality assumption is sub-optimal in general and particular as the dimension increases. We have conducted further studies with dimensions (not shown here for limitation of space) and the conclusions remain very similar.In order to further compare the performances of copula estimators with flexible degrees vs. equal degrees, we fit our Bayesian models with two different settings of priors: (i) (flexible) the original priors given in (6) and (7) where degrees are allowed to vary in different dimensions; (ii) (equal) modified priors given in (12) and (13) where degrees are set to be equal; using the same data set generated from the three bivariate copulas. Following Segers et al. (2017), we consider three global performance measures: the integrated squared bias, the integrated variance, and the integrated mean squared error. Given a copula estimator