Grid-Uniform Copulas and Rectangle Exchanges: Bayesian Model and Inference for a Rich Class of Copula Functions

Copula-based models provide a great deal of flexibility in modelling multivariate distributions, allowing for the specifications of models for the marginal distributions separately from the dependence structure (copula) that links them to form a joint distribution. Choosing a class of copula models is not a trivial task and its misspecification can lead to wrong conclusions. We introduce a novel class of grid-uniform copula functions, which is dense in the space of all continuous copula functions in a Hellinger sense. We propose a Bayesian model based on this class and develop an automatic Markov chain Monte Carlo algorithm for exploring the corresponding posterior distribution. The methodology is illustrated by means of simulated data and compared to the main existing approach.


Introduction
One of the primary interests of statistical analysis of multivariate data is the study of how random variables relate to each other.Among the various ways to express the relationship between random variables, one of the more flexible ones is the use of a marginals-copula representation, which provides a way to separate multivariate distributions into their single variate marginals and a function which represents their association structure, the copula function (Nelsen, 2007;Smith, 2011;Joe, 2014).For any d-variate distribution H : R d → [0, 1] with d > 1 and marginals F 1 , F 2 , . . ., F d , a copula is a function C : R d → [0, 1] such that where x = (x 1 , . . ., x d ) ∈ R d .Sklar's theorem (Sklar, 1959;Faugeras, 2013;Nelsen, 2007) is a classical result which states that this function C exists for any multivariate distribution H, and that if H is continuous, C is unique.Copula functions are themselves multivariate probability distributions supported on the unit hyper-cube and are such that the single variate marginals are uniform (Nelsen, 2007).
Modeling phenomena with copulas involves specifying models for marginals F 1 , F 2 , . . ., F d , and for C(x) separately.Modeling F i (x i ) is a matter of modeling single variable probability distributions, which is a well understood topic and that has been studied with great depth in both frequentist and Bayesian contexts (see, e.g., Müller et al., 2015).Our primary interest here is the modeling of the copula function C(x).There is a large body of literature studying parametric copula models arising from standard multivariate distributions such as the Gaussian and multivariate-t (Nelsen, 2007;Pitt et al., 2006;Joe, 2014;Choroś et al., 2010;Smith, 2011).Some popular copula models are members of the Archimedean family, which are single-parameter copulas satisfying for a function φ(• | θ) with a single parameter, known as the generator function of the Archimedean copula Nelsen (2007), where φ −1 is its inverse function.Archimedean copulas have been used for both frequentist and Bayesian analyses in the literature (McNeil and Nešlehová, 2009;Genest et al., 2011;Kaewsompong et al., 2020).
Choosing a class of copula models is not a trivial task and constraining the inference to parametric copula models can lead to wrong conclusions, because it reduces our ability to represent relationships between random variables.Motivated by these facts, different flexible approaches have been discussed in the literature.Non-parametric approaches for discrete data have been attempted in numerous ways, as described in Genest and Nešlehová (2007) and Yang et al. (2019).
In the context of continuous data, classical semiparametric approaches have been discussed by Genest et al. (1995), while classical nonparametric approaches can be traced back to Deheuvels (1979), and can be found in Mukhopadhyay and Parzen (2020).The classical methods commonly rely on the use on partial-or pseudo-likelihood and do not allow for a proper quantification of the uncertainties associated to the lack of knowledge of the marginal distributions.Furthermore, these approaches cannot be employed for modelling the association structure of latent variables in the context of hierarchical models.
Flexible model-based approaches for copula functions can be found in Guillotte and Perron (2012), Wu et al. (2013a), Wu et al. (2013b), and Ning and Shephard (2018).Guillotte and Perron (2012) proposed an interesting semi-parametric Bayesian approach for bivariate copulas based on a finite-dimensional approximation.This approximation is structured from a partition of the unit interval based of intervals of the same length [(i − 1)/m, i/m], where m is the number of intervals and i ∈ {1, . . ., m}.Their proposal is constructed using this partition and indicator functions for the corresponding intervals.The density of the copula is constructed via a mixture of the cross products of the indicator functions, resulting in a locally uniform distribution over the unit square, parameterized by a doubly stochastic matrix.Taking advantage of the properties of doubly stochastic matrices, the authors proceed to develop both a conjugate Jeffreys prior and a Markov chain Monte Carlo (MCMC) algorithm to sample from the posterior distribution.The approach is flexible but is difficult to generalize to not equally-spaced partitions and higher dimensions because of its reliance on the properties of doubly stochastic matrices.Wu et al. (2013a) proposed a model based on mixtures of Gaussian copulas.A Gaussian copula function is given by where Φ R is the CDF of a d-variate normal distribution with mean zero, variance one and covariance matrix R, arising from the corresponding correlation function, and Φ is the CDF of a standard normal distribution function.The authors claim that mixtures of bivariate Gaussian copulas can approximate any continuous bivariate copula function.Unfortunately, the Gaussian copula kernel is not rich enough to form a dense class.For instance, the density function of a bivariate Gaussian copula has the property that c R (u 1 , u 2 ) = c R (u 2 , u 1 ).Therefore, the density of a mixture of arbitrarily many Gaussian copulas also has this feature and, thus, cannot approximate an asymmetrical copula function, such as the asymmetrical t copula described by Church (2012).Wu et al. (2013b) extended the idea of a mixture of copula functions to the class of skew-normal copulas.However, there is no evidence that this class of copula functions is dense in the space of all copula functions.Finally, Ning and Shephard (2018) employed Dirichlet-based Polya trees models to propose a fully non-parametric Bayesian approach to modeling copula functions in any number of dimensions and a method for conjugate posterior simulation from the resulting posterior.This attractive result is, however, significantly marred by the flaw that the simulations from the posterior distribution are not themselves copula functions.
In this paper, we introduce a novel class of grid-uniform copula functions, which is dense in the space of all continuous copula functions.We propose a Bayesian model based on this class and develop an automatic MCMC algorithm for exploring the corresponding posterior distribution, allowing for the flexible modelling of continuous joint distributions.The paper is organized as follows.In Section 2 we introduce the class of grid-uniform copula functions and state its main properties, including its ability to approximate any given continuous copula function.In Section 3 we propose a Bayesian model based on the class of grid-uniform copulas and describe the MCMC algorithm.In Section 4, we illustrate the behavior of the model by means of analyses of simulated data.A final discussion concludes the article.

Definition
We begin by defining the class of copulas on which we develop our proposal.Let ρ be an orthogonal grid of [0, 1] d .Specifically, let ρ i be an ordered collection of points in [0, 1], and set Let ν ρ be the collection of sets formed by ρ, which are indexed by their upper right (or higher dimensional equivalent) corner.Now, let F be a probability measure defined on an appropriate space and B a measurable set such that F (B) > 0. We denote by F | B to the restriction of F to B defined by ρ-uniform copula if it is ρ-uniform and its one-dimensional marginal distributions are uniform.
A grid-uniform copula can be completely described by specifying the grid ρ and the probabilities for every B ∈ ν ρ .Hence, for each grid ρ, the space of grid-uniform copulas over this grid are a compact finite dimensional domain.For a grid ρ and a distribution C, we will use C ρ to denote the grid-uniform version of C, which assigns to each set B ∈ ν ρ the probability assigned by C to that set, i.e., C ρ (B) = C(B) for every B ∈ ν ρ .It is easy to see that if C is a continuous copula, then C ρ is also a grid-uniform copula and that the CDF of C ρ and C coincide for every y ∈ ρ.

Richness of grid-uniform copulas
We now prove that the class of grid-uniform copulas is sufficiently rich to approximate any arbitrary continuous copula function.
Theorem 1.Let C be an arbitrary copula which is absolutely continuous with respect to Lebesgue measure.Then for every > 0, there exists a grid-uniform copula D such that the Hellinger distance between D and C is smaller than , H(D, C) < .
PROOF: First, we will prove the theorem for the case when C admits a continuous density, denoted by c.For each grid ρ, let c ρ be the density of the the grid-uniform copula C ρ .Notice now that for every set in ν ρ , there is a point y such that c(y) = c ρ (y).The reason for this is that over the set, c and c ρ are two continuous functions with the same integral.Notice also that since c is continuous, it is bounded above by some bound b.Also, since c is continuous, it is uniformly continuous, i.e., for all points x, y and > 0, there is δ > 0 such that if ||x − y|| < δ, then |c(x) − c(y)| ≤ 2 /b.
We can now pick a grid ρ such that for all points x, y within the same set of ρ , ||x − y|| < δ.
By the above observation, it follows that |c It follows that the squared Hellinger distance between C and C ρ , is given by For the case when C does not admit a continuous density, we proceed with a somewhat similar idea.Let c be a density of C. Let E be the set of discontinuities of c.For every grid ρ, let ν ρ 1 be the collection of sets which have a discontinuity and ν ρ 2 be the collection of sets which do not have a discontinuity.Now, there is a grid ρ 1 such that the measure that C assigns to ν ρ 1 1 is less than /3, for any > 0. Furthermore, for every set in ν ρ 1 2 , C is uniformly continuous and has an upper bound b, so by the same argument used in the case when c was continuous, we can find a new grid ρ 2 such that |c ρ 2 (x) − c(x)| < 2 /9b for all x contained in one of the sets of ν ρ 2 .Now let ρ be any grid which is a refinement of ρ 1 and ρ 2 .It follows that the squared Hellinger distance between C and C ρ , is given by which completes the proof of the theorem.
We may note that none of the steps in the proof actually make use of the fact that C has uniform marginals, and it can easily be extended to prove that grid-uniform distributions can approximate any continuous distribution over a rectangular support.

Measures of association
In this section we provide the expression for two important measures of dependence for griduniform copula models.Specifically, we provide the expression for Kendall's tau and Spearman's rho, which are considered the best alternatives to the linear correlation coefficient as a measure of dependence for non-elliptical distributions, for which the linear correlation coefficient is inappropriate and often misleading.
Let C be a ρ-uniform copula function.Let C (i,j) be the bivariate marginal copula of C for the variables in the coordinates i and j.Let a (i,k) , k = 0, . . ., m i , be the kth element in ρ i and b (j,l) , l = 0, . . ., m j , be the lth element in ρ j .It is straightforward to show that C (i,j) is a ρ (i,j) -uniform copula function, where ρ (i,j) = (ρ i , ρ j ) and where S (a(i,k),b(j,l)) is the collection of sets in ν ρ such that ith coordinate of the index of the set is a (i,k) and the jth coordinate of the index of the set is b (j,l) .Spearman's rho, β, and Kendall's tau, τ , for the variables in the coordinates i and j is given by respectively, where .

Rectangle exchanges
We now introduce a class of transformations on grid-uniform copulas referred to as rectangle exchanges, which have the following important properties: (i) a rectangle exchange on a ρ-uniform copula produces another ρ-uniform copula, and (ii) given a grid ρ, and two ρ-uniform copulas C and D, there is a finite sequence of rectangle exchanges which can transform C into D.
Definition 2. Let ρ be a grid on [0, 1] d and C be a grid ρ-uniform copula function.The function C * is the result of a rectangle exchange of C, if C * is constructed using the following steps: (1) Set C * = C, and pick i and j in the set {1, . . ., d}, such that i < j and the cardinality of ρ i and ρ j is greater than or equal to 2. Also, for all k ∈ {1, . . ., d} \ {i, j} pick point (3) Set where l, m ∈ {1, 2}.
(4) Pick some in the interval The rectangle exchange operation is illustrated in two-dimensions in Figure 1.
[Figure 1 about here.] Rectangle exchanges are a closed operation on grid-uniform copulas since they conserve the uniformity of all marginals.
Lemma 1.Let G be ρ-uniform copula and G the resulting rectangle exchange of G, Then G is also a ρ-uniform copula.
PROOF: We will call the two coordinates of the two dimensional plane on which the exchange was performed, coordinate i and coordinate j, respectively.For considering marginals of coordinates other than i and j, we observe that all slices contain either no changed set, or all of the changed sets.For slices without any changed set, we note that the probability has not changed, and for slices with all of the changed sets, the magnitude of the change is − + − = 0, so the marginals along these coordinates remain uniform.
For coordinate i and j, we can assume, without loss of generality, that we are only worried about coordinate i, and that i is the first coordinate in the two dimensional slice on which the rectangle exchange was performed.Then slices in direction i come in two varieties: slices other than the ones along a 1 and a 2 contain no changed sets, so the probability of the slice remains unchanged.
The slice along a 1 contains the sets ν ρ p ( a 1 ,b 1 ) and ν ρ p ( a 1 ,b 2 ) so the probability of the slice changes by − + = 0, and slices along a 2 are similar.Hence, rectangle exchanges conserve uniformity for all marginals.
Another interesting property of grid-uniform copula functions is that starting from an uniform distribution on [0, 1] d , it is possible to reach any given grid-uniform copula, say C 0 , by doing certain type of operations.Specifically for a given grid-uniform copula C, we will refer as a grid division on C, to the addition of a division along any of the coordinates, such that the sets that are not divided retain their probabilities, and those that are divided distribute their probability in proportion to their volume.In other words, the resulting copula function arising from a grid division of C is identical to C, but is mapped onto a more refined grid.
Lemma 2. Let U be a uniform distribution in [0, 1] d , and let C be an arbitrary grid-uniform copula function.Then, there is a finite sequence of grid divisions and rectangle exchanges which will transform U into C.
PROOF: We will prove this for the two-dimensional case.For higher dimensions, the proof is identical, but the notation is more complex.We proceed by induction on the size of the grid.For C on grids up to 2 × 2 the result is obvious, in particular, on any grid of size 1 × d, the only grid-uniform copula is a completely uniform distribution.Now we assume that it is possible to transform U into Q for all grid-uniform copulas Q on grids of size m × k.Assume that C is a grid-uniform copula on a grid of size m × (k + 1) and we will prove that it is possible to transform Q into C.It is also necessary to consider the case of splitting the other coordinate, that is, the case where C is (m + 1) × k.However, the proof is identical.
Consider the distribution C such that C (A) = C(A) for A ⊆ i j<k ν ij and where the last two sets of each row have had the probability distributed uniformly between them.We note that C is a grid-uniform copula on an m × k grid.We will now show a sequence of steps, starting at C which will lead to C. Begin by splitting the last column along the last division in the grid of C, thus making the grids equal.Now we have a distribution which has the same grid as C and which is equal except for 2m sets.We will update C to be this new distribution.We note that for ).Hence, if the only transformations performed are rectangle exchanges, and if all of those that include ν ρ i,k also include ν ρ i,k+1 then this sum will remain constant.It can, therefore, be concluded that if only rectangle exchanges of this kind are used, it is sufficient to adjust the probability of the kth column.Now we will label the number of sets in the kth that are different in C and C as q.Observe that the sum of the column is the same in both C and C because they are copulas.Therefore, if C = C at least one set must have a probability which is higher than the probability of C and one set must have a lower probability.Pick the index of any of the sets with higher probability and call it α, and any of the lower sets and call it β.Now perform a rectangle exchange with This will yield a new grid-uniform copula with at most q − 1 probabilities in the kth column different from those in C.
The proof is completed by iterating this procedure as necessary.
These results allow us to prove that via rectangle exchanges is possible to generate the full space of grid-uniform copula functions.
Theorem 2. The C 1 and C 2 be two ρ-uniform copulas.There is a finite sequence of rectangle exchanges to transform C 1 into C 2 .
PROOF: Lemma 1 proves that if a distribution can be created by performing a rectangle exchange then it is a ρ-uniform copula.Lemma 2 shows us that it is possible to reach an arbitrary ρ-uniform copula, C 0 , from a Uniform distribution by means of rectangle exchanges and grid-divisions; however the grid divisions are not actually necessary.To see this, note that the uniform distribution U is already a ρ-uniform copula.We can follow the procedure from the previous lemma if we consider U in it's ρ-uniform representation.The difference is that some rectangle exchanges in the procedure described previously were performed between sets which are the union of several elements of ν ρ .At the time of the exchange, these unions of sets can be seen as a single superset with probability distributed uniformly.We can achieve the result of the larger exchange by performing exchanges with the component sets, adjusting the exchanged probability ( ) in proportion to the volume of the set, the details of exactly how to do this are explained next.
Let A l , A k , A m and A n be the sets involved in a larger rectangle exchange and let A i = A i,j sets which result from adding grid divisions.For any set A i , we will refer to its probability before the exchange as P i and its probability after the exchange as Q i .We proceed to describe how this same probability distribution of the rectangle exchange on A i can be achieved by performing smaller rectangle exchanges on the sets in A i,j .
We will consider the case where a single additional grid division is performed.In this case exactly two sets, A i are split.With no loss of generality, we can assume that A l and A k were split.
Note also that P l,1 = ξP l and P k,1 = ξP k where ξ is the volume of A l,1 divided by the volume of A l (and matches with A k and A k,1 ).Hence Q l,1 = P l,1 + ξζ and Q k,1 = P k, 1 − ξζ.We therefore perform a rectangle exchange with A m , A n , A l,1 , and A k,1 exchanging the probability ξζ and then another with A m , A m , A l,2 , and A k,2 exchanging the probability (1 − ξ)ζ.Thus, we have shown how a single rectangle exchange can be emulated by performing a rectangle exchange over a grid with an additional split.For refinements of the grid which add more than one additional split, add the splits one at a time.
The remaining issue is to prove that if C 0 and C 1 are both ρ-uniform copulas, then it is possible to perform rectangle exchanges to get from C 0 to C 1 .One way to do this is to reverse the steps to arrive at C 0 from a uniform copula, and then go from the uniform copula to C 1 , which completes the proof of the theorem.
There is something surprising going on here.Intuition would lead us to believe that rectangle exchanges would work for 2 dimensions, whereas higher dimensions would require parallelepiped exchanges.However, the surprising fact is that the above theorem is true regardless of the number of dimensions.In essence, the apparently very complex problem of exploring the space of grid-uniform copulas is solved in any number of dimensions by repeated transformations of the sort illustrated in Figure 1.
3 Bayesian modeling and inference using grid-uniform copulas Our ultimate objective is to use grid-uniform copulas to perform Bayesian statistical inference.We state the components of the Bayesian model in this section.Assume that we observe an independent and identically distributed (i.i.d.) sample of size n from a d-variate continuous distribution ∼ H, where H(y) = C(F 1 (y 1 ), F 2 (y 2 ), . . ., F d (y d )), with F 1 , . . ., F n being the marginal distributions of H, and C is the corresponding copula function.We model C as a grid-uniform copula function.Under the grid-uniform copula model, the log-likelihood function is given by: where

Grid-uniform prior models
Let ρ be a grid on [0, 1] d .Let C 0 be an arbitrary reference copula function and α > 0. We propose prior models for grid-uniform copula functions of the form where D be a suitable distance for probability distributions, and C ρ is the space of ρ-uniform copulas.Many choices for D could be considered.One choice that provides a simple interpretation of the hyper-parameters is the squared-L 2 distance.Let c 0 and c be densities for C 0 and C, respectively.Let B 1 , . . ., B p be the sets included in ν ρ .Under the squared-L 2 distance the grid-uniform prior model is given by where c l = B l c(x)dx λ(B l ) , with λ(A) being the Lebesgue measure of the set A. The prior takes the form of a truncated |ν ρ |-variate normal random distribution, centered at the ρ-uniform version of C 0 , C 0,ρ , and precision matrix given by α × I |ν ρ | .
C 0,ρ plays the role a centering parameter under the proposed prior and corresponds to the prior mode.On the other hand, α plays the role of a precision parameter, since as α → +∞, the prior variance var(π(C|ρ, α, C 0 )) → 0. The impact of α scales with the size of the sets in the grid, meaning that for fine grids, α may have to be very large.To facilitate the prior elicitation process we consider the parameterization α = α |ν ρ | .For purposes of calculation, we note that using C 0,ρ instead of C 0 as a reference copula produces the exact same prior, and calculating the prior using C 0,ρ is much simpler, since its value is constant throughout each cell in the grid.
Assigning D to the squared-L 2 norm provides nicely interpretable parameters.However, it does not allow for the incorporation of prior information on the degree of smoothness of the copula function.For continuous copulas, it is often reasonable to expect that the value of a copula density at a point is similar to the value at nearby points.Let V ρ be a set containing the elements of ν ρ in a given order.Let W be a symmetric matrix in which each entry W i,j encodes information about the spatial relationship of the sets V ρ i and V ρ j .Finally, let D W be a diagonal matrix with Borrowing ideas from models commonly used in spatial statistics and the nature of the grid-uniform model, we propose to take where γ > 0. Under this distance, the grid-uniform prior is given by where the order induced by V ρ .Notice that the proposed prior corresponds to a truncated Gaussian conditional autoregressive (CAR) model, which allows for spatial correlation of nearby values, as specified by a smoothing parameter γ > 0. It is worth noting that D is only a distance under certain conditions on W and γ (see, e.g., Banerjee et al., 2014).
A popular option is to set W such that W i,j = 1 d ij , where d ij is the distance between the centroids of V ρ i and V ρ j (Schmidt and Nobre, 2014).Another option is to set W as an adjacency matrix, where W i,j = 1 if the sets V ρ i and V ρ j are grid-neighbors (in the usual intuitive sense), and W i,j = 0 otherwise.For practical purposes, the sparseness of the adjacency matrix produces a prior which is faster to compute, and is a reasonable choice under most circumstances.When W is the adjacency matrix, then (D W ) i,i = |N V ρ i |, where N B is the collection of grid-neighbors of the set B and the distance reduces to the following expression Expressions for D as described above are not distances for all values of γ.To force D to be a distance, we can pick γ ∈ ( 1 λ 1 , 1 λn ), where λ 1 and λ n are the smallest and largest eigenvalues of W .As discussed by (see, e.g., Banerjee et al., 2014, section 6.4.3.3), the spatial correlation is low unless γ is close to 1.Because of this, a popular alternative is to consider γ = 1, which is known as the Intrinsic CAR (ICAR) model.Under the ICAR, D is given by In addition, when W is the adjacency matrix, the distance reduces to In general, when γ = 1 the latter expression is not a distance, but only to a pseudo-metric, since adding a constant to either c or c 0 does not change the value of D(C, C 0 ).However, since C and C 0 are both restricted to the space of ρ-uniform copulas, D does define a distance on the corresponding domain.Finally, note that all of the prior models share the algebraic structure of a truncated Gaussian distribution centered at C ρ 0 .In fact, all but the ICAR are exactly truncated Gaussian distributions.Therefore, the interpretation of C 0 and α remains intact.

On the choice of hyper-parameters
The prior depends on the choice of the grid ρ.The grid plays an equivalent role to the knots in the context of nonparametric regression based on splines.Rather than attempting to optimize the choice of a grid of reduced size and "well" located divisions, here we follow the approach proposed by Eilers and Marx (1996) in the context of penalized spline regression.Specifically, we consider an equally spaced and fine grid, along a penalization induced by an ICAR model, given in expression (3).The precise spacing of the grid can be chosen in relation to the available computational resources.We have found that on mid-range modern hardware (a 3.9 GHz processor) a good posterior estimate for a 50 × 50 grid (2401 free sets in the grid) can be computed in about 24 hours, whereas for a 10 × 10 grid (81 free sets) a posterior estimate can be computed in about 2 minutes.This is due not only to the greater computational cost of calculating the prior and likelihood functions, but also since a larger number of MCMC iterations are required.
The parameters α and C 0 have clear interpretations, and when prior information is available, it can be used to inform their choice.For situations when such information is not readily available, we propose suitable defaults.In this setting, selecting a single default C 0 around which to center the prior is difficult because, once specified, a single centering distribution may affect inference unduly.For instance, the use of the independent copula is highly informative because the lack of dependence is itself an extreme form of association structure.Rather than selecting a single centering copula function, one option is to consider a mixture of grid-uniform copula models by allowing the parameters of the centering copula function to be random.One possible choice is to use the Gaussian copula family given by parametrized by the correlation matrix R.
Choosing a prior for R is delicate since π(C|ρ, α, C 0 ) is known only up to a proportionality term, and this term depends on C 0 (and hence on R).We can write out the full prior for C as where N (R) is a normalizing constant.By default, we consider a conditional prior for R, such that . This choice is computationally convenient for simulation, as described in Section 3.3.3.It is difficult to find a closed form expression for this prior, but it is possible to characterize its behavior by means of simulation.Figure 3 illustrate the form of the prior in the bivariate case by considering 10 × 10 and 20 × 20 grids.We observe that this prior does depend slightly on ρ and α, but the overall distribution is not greatly affected by the changing to a grid that has four times as many cells.

about here.]
A similar procedure could be done to observe the prior for R in a higher dimensional case, but the results of this simulation of a prior over a high dimensional correlation matrix are quite difficult to interpret.Section 3.3.3describes the algorithm to perform this simulation in the general case, but since the prior is difficult to understand, our implementation uses a fixed independent C 0 in dimensions higher than 2.
It is also possible to use a similar structure to allow α to be random, but the proportionality term for π(C|ρ, α, C 0 ) also depends on α.A joint prior on α, R with no closed form expression would be nearly impossible to interpret.

An automatic MCMC algorithm
We propose a Metropolis-within Gibbs algorithm for exploring the posterior distribution of the copula functions and the parameters associated with the marginal distributions.The precise law of this selection does not matter so long as it is independent of C (b) ρ and every pair of coordinates has a positive probability of being selected.In practice, we will select them uniformly.b) Pick a 1 and a 2 from ρ d 1 and pick b 1 and b 2 from ρ d 2 .Also, for all k ∈ {1, . . ., d} \ {d 1 , d 2 }, pick x k ∈ ρ k , and set p (a l ,bm) = (x 1 , . . ., x i−1 , a l , , x i+1 , . . ., x j−1 , b m , x i+1 , . . ., x d ), where l, m ∈ {1, 2}.The precise law of these selections does not matter so long as it is independent of C (b) ρ and every rectangle along the selected coordinates has positive probability of being selected.In practice, we will select them uniformly.c) Pick uniformly in the interval and ρ ) to the candidate generating distribution induced by random rectangle exchange described by steps a) -d).An interesting property of this candidate generating distribution is that it is symmetric, which simplifies the computation of the acceptance probability.This is explained by the uniform selection of in the valid interval.The MH acceptance probability of the candidate C ρ is given by ρ , C 0 )) can also be further simplified to: In practice, we have found that this MH behaves well, with acceptance rates around 23%.

Updating the marginal distributions
There is no single technique which will work efficiently for the updating of the parameters of the marginal distributions, and tuning the posterior sampler may be difficult.However, there are some algorithms which are effective for a broad scope of distributions, and that can do reasonably good posterior exploration without having to worry about tuning.A good starting point is the t-walk (Christen et al., 2010), which is a general purpose sampler for parametric continuous distributions.
The t-walk is a MH algorithm which adapts to the scale of the target distribution and can sample well from most finite dimensional continuous distributions without tuning.

Updating the centering copula hyper-parameter
When working with a hierarchical prior that establishes a prior distribution for C 0 which depends on R, updating R can be done by adding a kernel to the MCMC chain.To update C 0 , we use a variation of the metropolized hit-and-run algorithm (Chen and Dey, 1998), which makes proposals that are always valid correlation matrices.
In our specific case, we allow δ to be a pre-specified tuning parameter.We consider values between 0.3 and 1.0, as discussed by Chen and Dey (1998).To update R, we propose a move from + H, by picking H as follows: (1) Let ξ (i) be the least eigenvalue of R (i) (2) Pick a sequence of i.i.d.standard normal variables z 1,2 , z 1,3 , . . .z d−1,d . ( . Also set h i,i = 0 for all i, and for i > j set h i,j = h j,i .Set The acceptance probability is given by where Our prior is selected so that log(N (R)) − log(N (R)) cancels and we are not hampered by our inability to calculate it.
Note that for a two dimensional copula, R depends only the single dimensional correlation coefficient, r and the hit and run algorithm reduces to a standard Random Walk Metropolis kernel (Robert and Casella, 2013).

Illustrations
We illustrate the behavior of the proposed model by means of the analysis of simulated data.
Functions implementing the MCMC algorithms employed in these analyses were written in Julia and are available upon request to the authors.

Estimation of parametric and non-standard copula functions
To illustrate that the proposal model does not overfit the data when a parametric copula model holds and that is able to capture deviations from the standard parametric models with finite sample sizes, we consider bivariate models with Gaussian (0, 1) marginals, under the following copula functions: • Model 1: A Clayton copula with parameter θ = 3, given by • Model 2: A Gaussian copula with correlation 0.5.
Figures 4, 5, and 6 display the true models under consideration.For each model, we simulate a single data set of size N = 500, 1000, 5000 and 10000.For each simulated dataset we fit our proposed model by considering a 50 × 50 grid, the hierarchically centered prior, with the ICAR correlation structure described in Section 3.1, and α = 400.In these analyses we assume the marginals distributions to be known.We create a Markov chain of (conservative) of size 2,000,000 using the automatic algorithm described in Section 3.3.We consider a burn-in period of 20,000 and a thinning of 1,000.Figures 4, 5, and 6 show the posterior mean under the different models and sample sizes.
[Figure 4 about here.] [Figure 5 about here.] [Figure 6 about here.] Figure 7 displays the posterior mean and 95% credibility intervals for the Hellinger distance to the true model under the different models and sample sizes.The results show that adequate estimates for complex true models can be obtained, even for reduced sample sizes, and that when the copula model is simple, the proposed model does not overfit the data.The results also show that the posterior mean gets closer to the true model when the sample size increases and that the posterior distribution concentrates around the true model as the sample size increases.

Comparison with existing approaches
We compare our proposal with the flat prior, which is one of those proposed by Guillotte and Perron (2012).There are similarities between our model and their flat proposal.The model proposed by Guillotte and Perron (2012) is a specific case of a grid-uniform copula, restricted to two dimensions and with grids that are necessarily evenly spaced.When these conditions hold, the models differs in way the prior probability mass is assigned.The work of Guillotte and Perron (2012) focused mainly on reference priors, whereas our prior is designed to share information on neighboring sets.
We compare our proposal with the flat prior proposed by Guillotte and Perron (2012) here.The flat prior of Guillotte and Perron (2012) can be thought as a limiting case of our model when α −→ 0.
We compare the models under settings considered by Guillotte and Perron (2012) in the evaluation of their proposal.We consider a two dimensional problem, a 6 × 6 grid, and a Gaussian, Gumbel and Clayton copula.We set the parameters of the different copula models such that they imply a similar association structure for the two variables.In particular, we set them such that they have the same Kendall's τ , and consider τ = 0.05, 0.35, 0.50, and 0.64.We consider three sample sizes in the comparison, N = 30, N = 100, N = 400, and N = 800.
We performed a Monte Carlo study, considering 100 replicates for each model and sample size.
For each data set we fit our proposal with the hierarchically centered prior, the ICAR correlation structure described in Section 3.1, and α = 40.The performance of the models was evaluated by computing the mean integrated squared error between the posterior mean and the true data generating copula model.The results are presented in Table 1.
[Table 1 about here.] The results illustrate that the proposed model outperform the flat prior across the board.As expected, the biggest differences between models are observed at small sample sizes; the larger the sample size, the smaller the difference between models regardless of the association structure.
Furthermore, for a given sample size, our model tends to produce better results than the flat prior as the level of association increases.

Concluding remarks
Flexible inference of copula functions had mainly relied on partial likelihood or pseudo-likelihood methods.This approach is useful in some cases.However, they do not allow for a proper quantification of the uncertainties associated to the lack of knowledge of the marginal distributions and cannot be employed for modelling the association structure of latent variables in the context of hierarchical models.We have proposed a novel and rich family of copula functions that can overcome these problems, the class of grid-uniform copula functions.We prove that this class is dense in the space of all continuous copula functions in a Hellinger sense.
We proposed a hierarchically centered prior distribution based on the proposed family, borrowing ideas from spatial statistics.We have described a class of transformations on grid-uniform copulas which is closed in the space of grid-uniform copula functions and that is able to span the complete space of grid-uniform copula functions in finite number of steps, starting from any point in the space.This family of transformations, referred to as rectangle exchanges, is employed to develop an automatic MCMC algorithm for exploring the corresponding posterior distribution.We have illustrated the behavior of the proposal and compared it with the approach proposed by Guillotte and Perron (2012).By considering similar simulation settings to the ones considered by Guillotte and Perron (2012), we show that our proposal outperforms their flat model when the posterior mean is the point estimator and mean integrated squared error is considered as a model comparison criteria.
The proposed prior model can be extended in different ways.The current proposal depends on a user-specified grid ρ.The size of the grid and the location of the points may have an important influence in the resulting model.For instance, equally spaced grids can lead to over fitting of the copula function in areas where few data points are "observed".On the other hand, they can lead to under fitting in areas where more data points are "observed".The study of strategies for the estimation of the optimal size and location of the grid is the subject of ongoing research.
The proposed model suffers from the curse of dimensionality.For a sample of size 1000, for a grip-uniform prior with a 10 × 10 grid it takes only a few seconds to generate a Markov chain of length 20,000 using the automatic MCMC algorithm and an i5 processor.We have also been able to use the proposed model in dimensions up to six.However, the implementation of the models in high dimensions and with fine grids would result in an explosion of parameters that need to be updated, which makes the implementation of this approach practically impossible.The study of marginal versions of the model, where the copula probabilities are integrated out of the model is also subject of ongoing research.
Finally, the extension of the model to handle mixed, discrete and continuous, variables and to copula regression problems is also subject of ongoing research.
B) is the Lebesgue measure of the set B, I A (B) is the indicator function that takes the value 1 if B ∈ A, and 0 otherwise.
Figure 2 illustrates the role of the parameters of the prior model.The figure displays the mean and credible interval of the prior distribution of the copula density at different points of the sample space.[Figure 2 about here.] 3.3.1 Updating C using a random rectangle exchange as a proposal Rectangle exchanges provide us with a way to explore the space of ρ-uniform copulas.This movement can be used to generate proposals in the context of an Metropolis-Hastings (MH) algorithm.Let C (b) ρ be the grid uniform copula which corresponds to the current state of the chain.Given, C (b) ρ , we propose the candidate C ρ using the following random rectangle exchange: a) Set C ρ = C (b) ρ and pick d 1 and d 2 randomly from the set {1, . . ., d}, and such that d 1 < d 2 .

Figure 1 :
Figure 1: Rectangle exchange -Illustration of a rectangle exchange for a ρ-uniform bivariate copula function C, where ρ 1 and ρ 2 have 7 equally-spaced points.Panel (a) illustrates step (1), where i and j are picked from {1, . . ., d}, and x 3 is picked from ρ 3 .Panel (b) illustrates step (2), where i, j ∈ {1, . . ., d} and b 1 , b 2 ∈ ρ j are picked and where the rectangles to be exchanged are selected.Panel (c) illustrates step (3), showing the assigned sets and the mass assigned to them by C. Panel (d) illustrates step (5), where the mass of the new copula function C * is computed.

Figure 2 :
Figure 2: Grid-uniform prior model -Prior mean (circle) and 95% equal-tail credible interval (vertical line) of the copula density evaluated at 100 equidistant points of the bivariate sample space.The results are shown for different values of α = α |ν ρ | .In all cases, the centering copula model, C 0 , is a bivariate Gaussian copula function with correlation equals to 0.5.Panel (a), (b), (c), and (d), provides the results for α = 2, 20, 200, and 2000, respectively.In each panel, the prior mode of the copula density is represented by triangles.

Figure 6 :
Figure 6: Model 3 (Copula of a Gaussian mixture model).Posterior mean of the bivariate density function.Panel (a) -(d) show the results for N = 500, 1000, 5000, and 10000, respectively.Panel (e) displays the true model.

Figure 7 :
Figure 7: Posterior mean (point) and 95% credibility interval (vertical bar) for the Hellinger distance to the true joint distribution, for different sample sizes.Panel (a) displays the results for Model 1 (Clayton copula).Panel (b) displays the results for Model 2 (Gaussian copula).Finally, panel (c) displays the results for Model 3 (a copula of a mixture of Gaussian distributions).

Table 1 :
Guillotte and Perron (2012)or (×10 3 ) for the posterior mean of the copula function under our default default hierarchical prior (proposal) and under the flat prior proposed byGuillotte and Perron (2012).The results are presented for for different true copula model, sample size (N ), and value of Kendall's τ .