Nonparametric Estimation of Multivariate Copula Using Empirical Bayes Methods

: In the ﬁelds of ﬁnance, insurance, system reliability, etc., it is often of interest to measure the dependence among variables by modeling a multivariate distribution using a copula. The copula models with parametric assumptions are easy to estimate but can be highly biased when such assumptions are false, while the empirical copulas are nonsmooth and often not genuine copulas, making the inference about dependence challenging in practice. As a compromise, the empirical Bernstein copula provides a smooth estimator, but the estimation of tuning parameters remains elusive. The proposed empirical checkerboard copula within a hierarchical empirical Bayes model alleviates the aforementioned issues and provides a smooth estimator based on multivariate Bernstein polynomials that itself is shown to be a genuine copula. Additionally, the proposed copula estimator is shown to provide a more accurate estimate of several multivariate dependence measures. Both theoretical asymptotic properties and ﬁnite-sample performances of the proposed estimator based on simulated data are presented and compared with some nonparametric estimators. An application to portfolio risk management is included based on stock prices data.


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. [17], Joe [19] and Nelsen [26] for more details about copulas and their applications.
Given a random vector (X 1 , . . ., X d ) with joint cumulative distribution function (CDF) F and continuous marginal CDFs F j , j ∈ {1, . . ., d}, by Sklar's theorem (Sklar [37]), the CDF F can be expressed uniquely as F(x 1 , . . ., x d ) = C(F 1 (x 1 ), . . ., F d (x d )), where C(•) denotes the copula function.A copula is itself the joint CDF of a random vector (U 1 = F 1 (X 1 ), . . ., U d = F d (X d )) having its marginals as uniform distributions on [0, 1], henceforth denoted by Uni f [0, 1].It is to be noted that the original results in Sklar [37] 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 F j , j ∈ {1, . . ., d} 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 C from a random sample of n independent identically distributed (iid) observations (X i1 , . . ., X id ) McNeil et al. [24], Nelsen [26], Smith [38] and Žežula [41], 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 [7]; Fermanian et al. [10]; Genest et al. [11]), or kernel-based methods such as local linear estimator (Chen and Huang [6]), mirror reflection estimator (Gijbels and Mielniczuk [12]) and improvements of these two estimators (Omelka et al. [27]).See alsoRémillard and Scaillet [30] and Scaillet and Fermanian [34] 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. [40].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 [33]), 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 [33].The asymptotic properties of the empirical Bernstein copula are well studied in Janssen et al. [15] and its application in testing independence is described in Belalia et al. [2].The application of the Bernstein copula to the modeling of dependence structures of non-life insurance risks is provided in Diers et al. [8], 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. [35] 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. [15] 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 [4] 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 [23] 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. [35], 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. [9]) having support on a possibly lowerdimensional 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.

An Empirical Bayes Nonparametric Copula Model
Suppose we have i.i.d.samples (X i1 , . . ., X id ) ∼ F(x 1 , . . ., x d ), i ∈ {1, . . ., n}, where F is a cumulative distribution function and F j is the absolutely continuous marginal CDF of the j-th component.By Sklar 's theorem (Sklar [37]), there exists a unique copula C(•) such that

and
(F 1 (X 1 ), . . ., The Bernstein copula is a family of copulas defined in terms of Bernstein polynomials and it was first introduced by Sancetta and Satchell [33].It is a flexible model that can be used to uniformly approximate any copula.The Bernstein polynomial with degrees (m 1 , . . ., m d ) of a function C : [0, 1] d → R is defined as and B m (C) 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 [33]) B m (C n ), where C n is the rank-based empirical copula.We denote the empirical Bernstein copula as where I(•) denotes indicator function and slightly modified empirical marginal distribution functions are defined as where the modification 1/(n + 1) instead of 1/n modifies the standard empirical marginal distribution to be away from 1 in order to reduce potential problems at boundaries.However, the empirical Bernstein copula C m,n is not guaranteed to be a valid copula for finite samples as the empirical copula C n is not necessarily a genuine copula.Segers et al. [35] show that the empirical Bernstein copula C m,n is a copula if and only if all the degrees m 1 , ..., m d are divisors of n.In order to obtain a valid copula estimation for any degrees, we replace the empirical copula C n with the empirical checkerboard copula C # n , which is a simple multilinear extension of the empirical copula defined as where R (n) i, j is the rank of X i j among X 1 j , . . ., X n j , see, e.g., Carley and Taylor [5] and Li et al. [22] for more details.Notice that the main difference between the empirical copula C n and the empirical checkerboard copula C # n is that C # n is a genuine copula, so we can obtain a valid copula estimation C # m,n taking the form where and we call the proposed empirical checkerboard Bernstein copula (ECBC).
Unlike the empirical Bernstein copula, the ECBC is a genuine copula for any degrees m 1 , m 2 , . . ., m d ∈ Z + and any fixed sample size n.It is known that Bernstein polynomials with smaller values of degrees m j '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 (m 1 , m 2 , . . ., m d ), where m j 's are allowed to be different for different j ∈ {1, . . ., d} and also depend on the random sample of observations.
As illustrated in Sancetta and Satchell [33], using partial derivatives of (2) with respect to each u j and rearranging, we can obtain the density corresponding to ECBC as follows: where 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. [13], the pseudo-observations (U i1 , . . .U id ), i ∈ {1, . . ., n} can be treated as samples from (F 1 (U i1 ), . . ., F d (U id )) ∼ C. We then use this approximation to build an empirical Bayesian hierarchical model: where a is the 'floor' function denoting the largest integer not exceeding the value a and , i.e., (V i1 , . . ., V id ), i ∈ {1, . . ., n} are samples from the empirical checkerboard copula Based on the proposition 1 in Genest et al. [11], V i j can be drawn using the following hierarchical scheme: ∼ Uni f (0, 1) i ∈ {1, . . ., n}, j ∈ {1, . . ., d}.
where F n, j (x j ) = 1/n n i=1 I(X i j ≤ x j ) and DisUni f denotes the discrete uniform distribution, i.e., Pr[π i = j] = 1/n for j ∈ {1, . . ., d}.Assuming that there are no ties in the pseudo samples U 1 j , . . ., U n j (owing to absolute continuity of marginal distributions or breaking it by random assignment in practice), we can equivalently represent the V i j 's more conveniently as Next, to account for the uncertainty in the estimation of the degrees m j 's, we propose to introduce a samplesize dependent empirical prior distribution on the degrees m 1 , . . ., m d 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 [4].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 [0, 1] d is assigned as the prior for the copula C in (1), which doesn't guarantee C 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. [16], we propose the hierarchical shifted Poisson distributions as the prior distribution for m 1 , . . .m d : and 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.Given the empirical priors distribution of m 1 , . . ., m d as in ( 6) and ( 7), and assuming the regularity conditions for the consistency of the empirical checkerboard copula, the proposed ECBC is consistent in the following sense: where the expectation is taken with respect to the empirical prior distribution.
First, under the assumption that the marginal CDFs are continuous, it follows from the Remark 2 in Genest et al. [11] that Next, notice that In above the second inequality follows from the fact that since and for any j ∈ {1, . . ., d}.Next by using the Lemma 1 in Janssen et al. [15] and equation (3) in Kiriliouk et al. [21], we obtain

Hence, it now follows that
Also, by using Lemma 3.2 in Segers et al. [35] we have Thus, combining the above inequalities that are satisfied almost surely (a.s.) for every fixed m j 's we obtain Next we consider the proposed empirical priors on the degrees m 1 , . . ., m d to be We make use of the following simple Lemma: Proof of Lemma 1: By Jensen's inequality for the square-root function, it follows that .
Notice that as α j ∼ Uni f (1/3, 2/3), Pr(α j > 1/3) = 1 and conditioning on α j , we then have by the above Lemma, E 1/m j α j ≤ (1 − e −n α j )/n α j → 0 as n → ∞.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 m j 's) weighted by the empirical prior distribution of the tuning parameters m j '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. [3]) 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 (m 1 , . . .m d ) 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 m j 's as final estimates.Let m j1 , . . ., m jK denote K MCMC samples of m j , j ∈ {1, . . ., d} and for each j, let m j1 < m j2 < • • • < m jD j denote the distinct values among these MCMC samples.Then the (marginal) posterior mode of m j is estimated by The final estimate of the smooth copula based on the proposed ECBC is then given by 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 m j 's be necessarily integer-valued) we have chosen to use posterior mode based on the marginal posterior distributions of m j 's.Through many numerical illustrations we show the easy applicability of this choice in various examples.

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 C, Spearman's rho can be written as A multivariate extension of Spearman's rho given in Nelsen [25] takes the form 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 where B 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 ρd 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 [28]).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 [25] as By applying ( 4) and ( 5), 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 m j , j ∈ {1, . . ., d} 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 d ≥ 2.
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 highdimensional integrals are often algebraically and even numerically intractable say for dimension d ≥ 5, 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).

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 copula in R (Hofert et al. [14]).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 −2, −1, 1 and 2, which reflects a wide range of dependence from negative to positive.The next two examples are the Clayton copula with parameter 1 and the Gumbel copula with parameter 2 The value of Kendall's tau is 0.33 for Clayton copula with parameter 1 and 0.5 for Gumbel copula with parameter 2. Both cases have a moderate positive dependence.
The next example is the independent copula C(u, v) = uv.Finally, we consider an asymmetric copula In the simulation study, n = 100 samples are drawn from the true copula for each replicate (of size n = 100) and there are N = 100 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 N 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 d = 2, d = 10, and d = 50, 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 N copula estimates are given for comparison and ( m1 , m2 ) 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.

Accuracy of Multivariate Dependence Measures (d = 3)
To assess the finite sample performance of the estimate of multivariate Spearman's rho ρd , we conduct Monte Carlo simulations for d = 3 copulas.We consider the independent copula and Clayton copulas with parameter value {0.5, 1, 2}, respectively.For each copula model, N = 100 Monte Carlo replicates are generated with size n = 100.For each replicate, we compute the proposed estimator ρd in (9) and the estimator based on empirical copula ρd 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 ρ d can be obtained by numerical integration since there is no analytical expression as a function of the parameter (see Pérez and Prieto-Alaiz [28]), which is also given in Table 1.The corresponding boxplots for the two estimates based on N = 100 replicates along with a horizontal line for true multivariate Spearman's rho ρ d are shown in Fig. 3.
From the results we can see our estimator ρd outperforms ρd with respect to variance and MSE.In terms of bias, ρd 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, ρd can take values out of parameter space [−2/3, 1] (3% and 12% of ρd are taking values greater than 1 in (c) and (d), respectively), which can be problematic in measuring dependence.

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 m j 's.We further explore the special case of choosing equal degrees m 1 = . . .= m d = m by using the following prior distribution 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 The next two choices are the independent copula (e.g., FGM with a = 0) and the Gaussian copula with positive dependence (correlation ρ = 0.5).The last one is the trivariate t-copula with degrees of freedom 4 and pairwise dependence ρ 12 = −0.2,ρ 13 = 0.5 and ρ 23 = 0.4.Again samples of size n = 100 are obtained for each four cases and repeated N = 100 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 (m 1 , m 2 ) or (m 1 , m 2 , m 3 ) for each chosen true copula model.From the plots we can observe that for bivariate copulas, posterior estimates of m 1 and m 2 are significantly different in most cases without any prior restrictions of equality.In fact, posterior probability of choosing equal m 1 = m 2 is only about 0.08, 0.08 and 0.13 for FGM copula, for independent copula, and for Gaussian copula, respectively indicating against forcing m 1 = m 2 as is popularly done in the literature.For the trivariate t-copula case, the posterior probability of m 1 = m 2 = m 3 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. [35], we consider three global performance measures: the integrated squared bias, the integrated variance, and the integrated mean squared error.Given a copula estimator Ĉn , the performance measures are defined as integrated squared bias: We compute the performance measures by applying the computation method described in Segers et al. [35], which relies on Monte Carlo simulation to get a Monte Carlo estimate of each performance measure.Table 2 presents the results for the four cases, where the first three are bivariate (d = 2) and the last one is trivariate (d = 3).The standard errors of the Monte Carlo estimates are not reported in the table as they are negligibly small.We can see that the copula estimators with flexible degrees perform better than those with equal degrees in terms of IV and IMSE in all cases.As the difference in IMSE is dominated by the IV term, the choice of flexible degrees leads to smaller uncertainty, and hence smaller IMSE while the biases remain relatively unaffected.It is also interesting to observe that estimated degrees are far smaller than sample sizes indicating the empirical beta copula may not have optimal performance, which we next explore.

Comparison with the Empirical Bernstein Copulas
In this section, we compare the finite-sample performance of the ECBC with other nonparametric copula estimators only as it has been already demonstrated that parametric models based methods lead to biased estimates under model misspecification.First, we consider the empirical beta copula introduced by Segers et al. [35], which is a special case of the empirical Bernstein copula where the degrees of the polynomials are set equal to the sample size.The empirical beta copula is a genuine copula and has been shown to outperform the classical empirical copula and the empirical checkerboard copula in terms of bias and variance.For bivariate cases, we also include the empirical Bernstein copula with degrees as suggested in Janssen et al. [15] into the comparison.By setting degrees m 1 = m 2 = m and minimizing the asymptotic pointwise mean squared error with respect to m, Janssen et al. [15] suggested the choice of m in the bivariate case as where and C u j and C u j u j are the first-order and second-order partial derivatives of C, respectively, with respect to u j , j = 1, 2.
Note that even we use the integer part m 0 (u 1 , u 2 ) in practice, it is not necessarily a divisor of n, meaning the empirical Bernstein copula with m = m 0 (u 1 , u 2 ) is not guaranteed to be a genuine copula.We consider the same copula models as in Section 3.3.The choice of degrees m 0 (u 1 , u 2 ) suggested by Janssen et al. [15] is not defined for the independent copula as C u j u j = 0 and it is restricted to bivariate cases, so we only take the empirical Bernstein copula into consideration for the bivariate FGM and Gaussian copulas.All results in Table 3 are based on N = 100 MC replications each of sample sizes n = 25, 50, 100.We compare the performance of the ECBC with flexible degrees (referred to as flexible ECBC), the empirical beta copula (referred to as Beta), and the empirical Bernstein copula with m = m 0 (u 1 , u 2 ) (referred to as Bernstein) using the same performance measures as in Section 3.3.
Table 3 indicates that the ECBC with flexible degrees outperforms the empirical beta copula in terms of variance and mean square error in all cases.Compared to the empirical Bernstein copula with m = m 0 (u 1 , u 2 ), the ECBC with flexible degrees has a smaller bias but the ordering with respect to mean square error is not clear between these two copula estimators.For small samples, the empirical beta copula seems to have the largest variance while the bias of the empirical Bernstein copula with m = m 0 (u 1 , u 2 ) is shown to be the largest even though it uses optimal "true" degree given in (14).
Table 3: Comparison of the ECBC with flexible degrees (referred to as flexible ECBC), the empirical beta copula(referred to as Beta) and the empirical Bernstein copula with m = m 0 (u 1 , u 2 ) (referred to as Bernstein) using three performance measures computed by Monte Carlo simulation based on N = 100 replications for sample size n = 25, 50, 100.

Application to Portfolio Risk Management
Copulas have been widely used in portfolio optimization and risk measurement as they are powerful tools to model the dependence among different assets in a portfolio.The proposed ECBC is capable of estimating multivariate copula and it is straightforward to sample from the estimated copula, so it can be applied to find optimal weights and estimate risk measures for a portfolio with a variety of assets.
We now illustrate the use of ECBC for portfolio risk allocation using real data consisting of a d asset values.Value at risk (VaR) and conditional value at risk (CVaR) (also called expected shortfall (ES)) are common measures of risk in the field of risk management (see, e.g., Jorion [20] and Uryasev [39]).Assume that X is the return of a portfolio or asset (daily log-return of a portfolio of stocks or individual stocks, with positive indicating profit and negative values representing loss) with distribution function F X (•) = Pr[X ≤ x].The VaR of X at the level of α ∈ (0, 1) is defined as while the CVaR (ES) of X is defined as Notice that, if we consider the corresponding loss of the same portfolio represented by Y = −X, then we have ). Mean-CVaR portfolio optimization is a popular portfolio optimization technique introduced by Rockafellar et al. [32].The advantage of Mean-CVaR portfolio optimization is that it calculates VaR and minimizes CVaR simultaneously where the optimization can be formulated as a linear programming problem.
Let x ∈ R d denote a realized return value of d assets in a portfolio, and v ∈ S d = {v ∈ R d : v j ≥ 0, ∀ j, d j=1 v j = 1}, denote the portfolio weights to be determined within the d-dimensional simplex S d .The key to the approach in Rockafellar et al. [32] is the auxiliary function for CVaR taking the form of where l(v, x) = −v T x is a linear loss function and F(x) is the joint distribution function of daily (random) return vector X which we will estimate using our proposed ECBC based empirical Bayes method.It has been shown in Theorem 1 of Rockafellar et al. [32] that for any weights v, H α (v, γ) is convex as a function of γ and is equal to CVaR α (v) at the minimum point.Moreover, VaR α (v) would be the left endpoint of arg min γ H α (v, γ).Moreover, Minimizing CVaR α (v) with respect to v is equivalent to minimizing H α (v, γ) with respect to (v, γ) (e.g., see Theorem 2 of Rockafellar et al. [32] for details).To numerically approximate the integral in (15), it is often good enough to generate M samples from F(•) or its estimate, which can be done by using the proposed empirical Bayes method based on the ECBC.However, a relatively less answered question in finance is how large should we choose M for accurate estimation as the integral in (15) depends on sampling the tail part of F(•) or its estimate.The empirical estimate of H α (v, γ) based on generating x k iid ∼ F or F can be written as where (•) + = max(•, 0).
Proposition 1.In order to achieve an accuracy of > 0 for the MC approximation, it is sufficient to generate M MC samples such that where Σ = Var F (X) and λ max is the largest eigenvalue of Σ.
Proof: By the Law of the Iterated Logarithm (see, e.g., Balsubramani [1]), the deviation of MC approximation from the mean is almost surely bounded by Let Σ = Var F (X) and λ max be the largest eigenvalue of Σ, then we have Notice that Σ and hence λ max can be easily estimated from the observed return values without any modeling assumption as long as n > d , however sparse methods are necessary for large size portfolios when n ≤ d.Next it can be shown that minimizing ( 16) is equivalent to minimizing Thus, along with the linear constraints on the weights v, it can be formulated as a linear programming problem and can be solved using standard convex optimization methods.Conveniently, R function BDportfolio optim within the package PortfolioOptim can be used for this purpose.Following the algorithm in Semenov and Smagulov [36], simulated return values can be obtained using estimated ECBC for portfolio optimization.The complete algorithm is summarized below: Step 1. Transform assets' historical data X t j to pseudo-observations U t j and estimate copula using our proposed method.
Our copula estimator is useful to find optimal weights and estimate risk measures as sampling from the estimated copula is straightforward.Considering the Bernstein copula density given in ( 4) and ( 5), we can obtain samples

, d}
As an example with moderately large dimension, we investigate the time series of daily closing stock prices of 10 top Nasdaq companies: AMZN, FB, GOOGL, AAPL, MSFT, INTC, CSCO, NFLX, CMCSA, and ADBE for the time period from January 1, 2018, to December 31, 2019.This data set consists of 502 observations and can be obtained using R package quantmod.
Suppose we want to find an optimal portfolio of stocks above that minimizes the expected shortfall of the portfolio.First, we convert the price series P t j to log-returns X t j X t j = ln P t j P (t−1) j , t = 1, . . ., T, j ∈ {1, . . ., d}, resulting in T = 501 log-return values for d = 10 assets.Then we follow Steps 1-3 as above to obtain the optimal portfolio weights and the corresponding VaR and CVaR.Similar to Semenov and Smagulov [36], we set the minimum weight to be limited by v j ≥ 0.01, j ∈ {1, . . ., d} to avoid corner portfolio cases.
In Step 1, the posterior mode estimator of the degrees of the proposed ECBC are (209,206,208,206,208,209,211,210,209,208).The largest eigenvalue of the covariance matrix is λ max ≈ 2 × 10 −3 , so we are able to find the value of M that is sufficient for a given accuracy from the relationship in (17).
We set M = 10000 (adequate for an accuracy ≈ 9 × 10 −4 ) and repeat Step 2-3 N = 100 times to quantify estimation uncertainty.For each replicate we conduct portfolio optimization at the level of α ∈ {0.10, 0.05, 0.01} as popularly used.As a result, we are able to obtain the distribution of optimal weights (Fig. 5) and risk measures (6) using simulated data from the estimated copula.
From the boxplots of optimal weights in Fig. 5 we can see that CMCSA has a much higher weight than the other stocks in the mean-CVaR optimal portfolio across different levels.Also, by applying mean-CVaR portfolio optimization to the historical log-return data X t j , we can get estimates of optimal weights and risk measures as well.In Fig. 6, the dashed lines indicate empirical estimates of risk measures using historical data.
We can see from the plots in Fig. 6 that the estimated risk measures from two different methods seem to be fairly close.However, we are able to quantify the uncertainty for all the estimates by repeatedly sampling from the estimated copula.Semenov and Smagulov [36] conduct a similar stability study to report the means and SDs of VaR and CVaR, but they used predetermined weights based on historical data and didn't report the distribution of optimal weights obtained from simulated data.Our copula estimator has shown good performance for relatively small samples and operationally we can generate as many samples as we want from the estimated copula, thus the copula-based method would be more reliable when there is not sufficient historical data.Besides, compared to the empirical estimates, it is possible to estimate VaR and CVaR for much smaller values of levels using the copula-based method.

Concluding Remarks
In this paper, we propose the empirical checkerboard Bernstein copula, which is a nonparametric multivariate copula estimator.It can be considered as an advancement of the empirical Bernstein copula since it is a valid copula with any polynomial degrees for any sample size.For automatic data-dependent dimension-varying degree selections, we further developed an empirical Bayesian method that has been shown to be practically useful.While the proposed copula estimator is shown to be large-sample consistent, it also has a good finite-sample performance.Moreover, it has a beneficial effect on measuring the strength of dependence for large dimensions because the estimates derived from the proposed copula are always within the proper range.
As sampling from the estimated copula is quite straightforward, it is applicable to portfolio optimization and risk measurement where estimation is often done with simulations generated from copulas.We investigate the number of simulations that are good enough to achieve any given accuracy, which has been apparently out of reach in the literature.Furthermore, we are able to provide uncertainty quantification for all the estimates in portfolio risk management.
Under the hierarchical structure of the proposed empirical Bayes model, MCMC methods have been shown to work reasonably fast for relatively large dimensions (d ≤ 50) with a moderate sample size (n = 100).To speed up the MCMC methods for very large sample sizes, it would be of interest to explore some scalable MCMC methods such as divide-and-conquer approaches and subsampling approaches (see, e.g., Quiroz et al. [29], Robert et al. [31], etc.).The codes (written using R software) to implement the procedure is available upon request from the first author and could be made available on a GitHub page following the publication of this paper.

Proof:
We denote the ECBC as B m (C # n ) for simplicity.Also, the empirical Bernstein copula and the Bernstein copula are denoted as B m (C n ) and B m (C), respectively.Let ||g|| = sup u∈[0,1] d g(u) denote the supremum norm of a function g(•) defined on d-dimensional square [0, 1] d .Using the familiar triangle inequality we have

Fig. 1 :
Fig. 1: Estimation of Frank copulas using the ECBC with empirical Bayesian method for choosing proper degrees when sample size n = 100.

Fig. 2 :
Fig.2: Estimation of various copulas using the ECBC with empirical Bayesian method for choosing proper degrees when sample size n = 100.

Fig. 3 :
Fig. 3: Boxplots of the two estimators based on N = 100 replicates and a horizontal line of true multivariate Spearman's rho ρ d for each of four trivariate copulas.

Fig. 4 :
Fig. 4: Choice of dimension-varying degrees obtained by applying the proposed empirical Bayesian method based on N = 100 replications when sample size n = 100.

Fig. 5 :
Fig.5: Distribution of optimal weights obtained from simulated data using the estimated copula at the level of 0.10, 0.05, and 0.01, for a portfolio of d = 10 stocks.

Fig. 6 :
Fig. 6: Distribution of VaR and CVaR obtained from simulated data using the estimated copula at the level of 0.10, 0.05, and 0.01 for a portfolio of d = 10 stocks.Dashed lines indicate empirical estimates using historical data.

Table 1 :
Comparison of two estimates of multivariate Spearman's rho based on N = 100 replications of size n = 100 generated from the independent copula and Clayton copulas with different parameter values when dimension d = 3.

Table 2 :
Comparison of copula estimators with flexible degrees vs. equal degrees using three performance measures computed by Monte Carlo simulation based on N = 100 replications when sample size n = 100.