Implicit Generative Copulas
Year: 2022

Abstract:
Copulas are a powerful tool for modeling multivariate distributions as they allow to separately estimate the univariate marginal distributions and the joint dependency structure. However, known parametric copulas offer limited flexibility especially in high dimensions, while commonly used non-parametric methods suffer from the curse of dimensionality. A popular remedy is to construct a tree-based hierarchy of conditional bivariate copulas. In this paper, we propose a flexible, yet conceptually simple alternative based on implicit generative neural networks. The key challenge is to ensure marginal uniformity of the estimated copula distribution. We achieve this by learning a multivariate latent distribution with unspecified marginals but the desired dependency structure. By applying the probability integral transform, we can then obtain samples from the high-dimensional copula distribution without relying on parametric assumptions or the need to find a suitable tree structure. Experiments on synthetic and real data from finance, physics, and image generation demonstrate the performance of this approach.

text: Copulas are a powerful tool for modeling multivariate distributions as they allow to separately estimate the univariate marginal distributions and the joint dependency structure. However, known parametric copulas offer limited flexibility especially in high dimensions, while commonly used non-parametric methods suffer from the curse of dimensionality. A popular remedy is to construct a tree-based hierarchy of conditional bivariate copulas. In this paper, we propose a flexible, yet conceptually simple alternative based on implicit generative neural networks. The key challenge is to ensure marginal uniformity of the estimated copula distribution. We achieve this by learning a multivariate latent distribution with unspecified marginals but the desired dependency structure. By applying the probability integral transform, we can then obtain samples from the high-dimensional copula distribution without relying on parametric assumptions or the need to find a suitable tree structure. Experiments on synthetic and real data from finance, physics, and image generation demonstrate the performance of this approach.

heading: 1 Introduction
text: Many approaches are available to reliably model univariate distributions of real-valued variables. One can use various parametric distributions or employ flexible, non-parametric methods, e.g., by using empirical distribution functions, kernel density estimation (KDE), or quantile regression [22]. Modeling distributions in higher dimensions is a much harder task. Parametric approaches using, e.g., Gaussian distributions, are common but inflexible. Many non-parametric approaches become infeasible in higher dimensions due to the curse of dimensionality [18].
Recent successes in modeling multivariate distributions have been achieved with deep neural networkbased likelihood-free implicit generative models [28], such as Generative Adversarial Networks [13] (GANs) and Generative Moment Matching Networks (GMMNs) [10; 24].
Copulas are a tool to decouple the modeling of the univariate marginal distributions from modeling the high-dimensional joint dependency structure [19]. This allows to obtain high-quality marginal models with the mentioned univariate techniques. Moreover, the (conditional) marginals can be tailored to the problem at hand and many fields have developed sophisticated domain-specific univariate models for this task, e.g., for probabilistic forecasts in finance [4], weather [36], or energy [17]. These models can then be combined with a suitable copula structure to enable the simulation of multivariate quantities, e.g., [34; 29; 41]. In the machine learning literature copulas have been used to increase the flexibility of Bayesian networks [11], to model multi-agent coordination in reinforcement learning [43], or for image generation via the Vine Copula Autoencoder [40]. Another popular application of copulas is synthetic tabular data generation [33; 27]. Note that in many of these applications, the main goal is to sample from a multivariate distribution and the copula is used as a building block for a generative model.
35th Conference on Neural Information Processing Systems (NeurIPS 2021), virtual.
High-dimensional copula distributions are commonly modeled via parametric structures such as the Gaussian or the student-t copula [8]. [25] propose to enhance the subclass of Archimedean copulas by learning the generator functions via deep neural networks. However, the family of Archimedean copulas is of limited use in higher dimension as it makes restrictive symmetry assumptions. A more flexible, semi-parametric state-of-the art approach is to use vine copulas which are built from trees of bivariate pair-copulas [1]. These pair-copulas can again be parametric, such as the Gumbel or Clayton copula, or non-parametric, e.g., by using bivariate kernel density estimation [32].
Directly applying implicit generative modeling to the task of estimating copula distributions is not straight-forward since copula distributions are required to have uniform marginals. Such an approach has been proposed by [23; 16], but without ensuring the uniformity property, at least not for finite sample sizes when the model is not expected to exactly fit the true copula distribution. Ensuring the marginal uniformity of the learned copula distribution is important since deviations from uniformity will result in unwanted alterations of the marginal distributions in the data space.
In this paper we show how to design and train implicit generative copula (IGC) models to match the dependency structure of given data. A key challenge is to ensure the marginal uniformity of the estimated distribution. We achieve this through first learning a latent distribution with unspecified marginals but the same dependency structure as the training data. We then obtain the desired copula model by applying the probability integral transform component-wise to the latent distribution. The IGC model and the data distribution are matched in copula space using the energy distance [39]. During training, the probability integral transform is approximated through a differentiable softrank layer which allows to use gradient-based methods.

heading: Our contributions
text: • We propose the first universal, non-parametric model for estimating high-dimensional copula distributions with guaranteed uniformity of the marginal distributions.
• We show how flexible likelihood-free implicit generative models based on deep neural networks can be trained for this task through the use of a differentiable softrank layer.
• Compared with the state-of-the-art semi-parametric vine copula approach, we demonstrate similar or improved performance for different tasks.
We present the proposed copula model in Section 2. Model training is described in Section 3. Section 4 shows experimental results for synthetic and real data from finance, physics, and image generation. We conclude in Section 5.

heading: 2 Proposed IGC Model
text: We introduce our IGC model following the schematic shown in Figure 1. Figure 2 provides an exemplary application for a bivariate two-component Gaussian mixture data distribution.
Copula Basics The vector-valued continuous random variable X ∈ RD represents the data source of our task. We denote its distribution by P and the cumulative distribution function (cdf) of X by FX(x). Moreover, let FXd(xd) be the cdf of the marginal distribution of Xd, the d-th component of X, for d = 1, . . . , D.
Each sample vector x ∈ RD can be mapped to a vector u ∈ [0, 1]D by defining as ud the value of xd with respect to the corresponding marginal cdf of Xd, i.e., ud = FXd(xd) for all d = 1, . . . , D. This operation is called the probability integral transform (PIT). In the copula literature, values obtained by the PIT are called pseudo-observations [19]. The random variable U follows a so-called copula distribution PC by construction, i.e., the distribution is defined on the unit space [0, 1]D and its marginals are uniform. The joint cdf of U is typically called the copula, which we denote by FU(u) here. Sklar’s theorem states that such such a copula function exists for all random variables X and it holds that FX(X) = FU(FX1(X1), . . . , FXd(XD)), i.e., any multivariate distribution can be expressed in terms of its marginals and its copula [38].
IGC Model In this work, we aim at defining a flexible, non-parametric model for copula distributions QCθ in high dimensions. Such a model can either be used to (approximately) represent the
copula distribution PC of the unit space vectors U or the distribution P of the original data vectors X, in which case the U vectors have to be transformed component-wise with the inverse cdfs of the components Xd, i.e., xd = F−1Xd (ud).
A flexible model class for learning high-dimensional probability distributions is the class of implicit generative models [28]. Here, latent random variables with a simple and known distribution are transformed via a parameterizable mapping, typically a deep neural network, to a complex distribution that approximates the distribution of the training data. A straight-forward application of this framework to model copula distributions, as done in [16] and [23], is difficult. The output of the generative model can be guaranteed to lie in [0, 1]D by applying an appropriate output layer, e.g., using a sigmoid function. However, guaranteeing the desired uniformity of the marginal distributions is not straight-forward. We have first experimented with additional training losses that penalize deviations from marginal uniformity. However, a simpler approach without any hyper-parameters for weighting additional cost terms and with a guarantee of marginal uniformity is the following two-step procedure.
In the first step, we model the distribution of the vector-valued latent random variable Y ∈ RD. To this end, we start with random variables Z ∈ RK from a simple base distribution that we choose as zero mean, unit variance Gaussian here. The tuneable map gθ : RK → RD with parameters θ is then used to transform the noise samples z to the latent samples y as
y = gθ(z). (1)
We denote the resulting probability distribution of Y by Qθ and the marginal cdfs by FYd(yd).
In a second step, we transform the latent variables Y into unit space vectors V. To this end, we define component-wise for d = 1, . . . , D
vd = FYd(yd). (2)
This transformation step is analogous to going form X to U. It guarantees both that the sample values vd from the model lie in the unit interval [0, 1] and that they are uniformly distributed, regardless of the distribution Qθ. We refer to the distribution of V as the model copula distribution and denote it by QCθ .
Interpretation Note that the distributions Qθ and P might be very different. More precisely, the mapping (1) can result in arbitrary marginal distributions but an optimal model of the distribution P would have
PC = QCθ . (3)
Given a sufficiently rich function class for gθ any distribution of Y can be modeled with small error [26]. This statement naturally extends to universal approximation properties for the dependency structure of Y, i.e., the unit space random vectors V. The approximation of the data copula distribution PC with our model copula distribution QCθ will thus be arbitrarily close if the generator function gθ is flexible enough and if the optimal θ is selected.

heading: 3 Training Procedure
text: Problem statement We aim at empirically estimating an IGC model QCθ from data, i.e., we want to determine the optimal parameter vector θ given a set of N training samples x1, . . . ,xN from X. The underlying target is to optimally match the copula model QCθ to the (typically unknown) copula distribution PC from which the samples were generated.
Parameter estimation We first transform the given data samples into unit space vectors in [0, 1]D, before starting the actual training procedure. Since we typically do not know the exact marginal cdfs FXd of the components Xd, we estimate them from the given data. To this end, we use the empirical cdfs separately for each component and define for n = 1, . . . , N and d = 1, . . . , D,
und = 1
N N∑ i=1 1[xid ≤ xnd ]. (4)
Here, 1[A] is the indicator function of event A, i.e., it is 1 iff A is true. Note that this is an unbiased and consistent estimator for the true cdf [42].
We then aim at matching the copula distribution QCθ to the unit space vector samples u 1, . . . ,uN . We base the inference of our likelihood-free model on the energy distance ED2 [39] between distributions PC and QCθ ,
ED2(PC , QCθ ) = 2E∥U−V∥ − E∥U−U′∥ − E∥V −V′∥. (5)
Here, U,U′ and V,V′ are independent copies of a random vector with distribution PC and QCθ , respectively. It holds that ED2(PC , QCθ ) = 0 if and only if P C = QCθ [39].
In practice, we use a sample-based approximation of the energy distance. We draw M samples v1(θ), . . . ,vM (θ) from QCθ and minimize the loss function
LED2(θ) = 1
NM N∑ n=1 M∑ m=1 ∥un − vm(θ)∥2 − 1 2M(M − 1) M∑ m=1 M∑ m′=1 ∥vm(θ)− vm ′ (θ)∥2. (6)
Since PC is independent of θ, the second term of (5) does not need to be included. Note that the energy distance is an instance of the more general maximum mean discrepancy (MMD) measure [15; 37]. We also experimented with MMD losses based on the Gaussian kernel but found the energy distance to be more robust as it does not require to tune the kernel bandwidth.
To generate samples from our copula model QCθ , we first draw samples y 1, . . . ,yM from the implicit generative model Qθ . To obtain vectors in unit space, we then again require the cdfs FYd(yd) of the marginal distributions of the components Yd of Y, d = 1, . . . , D. However, these functions are not known during model training and will likely change after each gradient step. Hence, we again use the empirical distribution functions to define for m = 1, . . . ,M and d = 1, . . . , D,
vmd (θ) = 1
M M∑ j=1 1[yjd(θ) ≤ y m d (θ)]. (7)
This provides us with an unbiased estimate of the current cdf during training. However, the indicator operation in (7) is not differentiable and hence will not allow the gradients to flow through this
operation. Therefore, during training, we replace (7) with a softrank layer based on a scaled sigmoid function
vmd (θ) = 1
M
( 0.5 + M∑ j=1
1
1 + exp(α(ymd (θ)− y j d(θ)))
) , (8)
where α is a scaling constant [35]. For sufficiently large α, (8) provides a close approximation to the empirical marginal cdfs at the current training step. A larger M will result in a finer approximation of the true marginal cdfs but comes with the increased cost for computing the ranks which requires O(M2) operations for all samples.
Sampling from the trained model Once we have completed training and thus have fixed θ, we have also fixed the component-wise cdfs FYd(yd). Now, the operation to estimate FYd(yd) does not need to be differentiable anymore and we can chose any available method to estimate the univariate marginal distributions based on samples from Qθ. We choose to simply draw a very large set of samples from Qθ and then use the empirical marginal cdfs, i.e., for each sampled value ytd we store the corresponding value utd according to (4), t = 1, . . . , T . Given a new realization of Yd, we obtain its approximate cdf value F̂Yd(yd) by interpolation of the stored values. Other more compact approximations like univariate kernel density estimation would also be feasible.
In sum, we obtain a sample vi from QCθ at test time by first sampling a realization z i from N(0, I), transform it via yi = gθ(zi), and then apply the component-wise PIT approximation to obtain vi = [F̂Y1(y i 1), . . . , F̂Yd(y i D)]. Given estimates for the component-wise marginal cdfs of the data F̂X1 , . . . , F̂Xd we can derive a sample s i in data space via si = [F̂−1X1 (v i 1), . . . , F̂ −1 Xd (viD)].

heading: 4 Experiments
text: In the following we empirically demonstrate the capabilities of IGC models on a series of experiments on synthetic and real data with increasing complexity. Additional experiments can be found Appendix B.
Implementation For all experiments except the image generation task, we use a fully connected neural network with two layers, 100 units per layer, ReLU activation functions, and train for 500 epochs. For the image generation experiment we use a three layer, fully connected neural network with 200 neurons in each layer, and train for 100 epochs. In all cases we train with a batch size of Nbatch = 100 and generate M = 200 samples from the model per batch. The number of noise distributions is set as K = 3D and T = 106. We use the Adam optimizer [20] with default parameters.
All experiments besides the training of the autoencoder models were carried out on a desktop PC with a Intel Core i7-7700 3.60Ghz CPU and 8GB RAM. For the training of the autoencoders we used Google Colab [14]. Training times for all experiments are in the range of a few minutes except for the image generation task. Details are provided in Appendix B. We use tensorflow with the Keras API for the neural networks [2; 7]. For copula modeling, we use pyvinecopulib [31] in Python and the R package kdecopula [30]. Our code is available from https://github.com/TimCJanke/igc.
Evaluation We evaluate the models by comparing the distance between the learned and the true copula. More specifically, we use the integrated squared error (ISE) in unit space, i.e.,
ISE = ∫ [0,1]D (FU(w)− FV(w))2dw ≈ 1 N N∑ n=1 (FU(u n)− FV(un))2, (9)
where FU is the joint cdf of U, the true copula, and FV the joint cdf of the learned model. Since analytical integration is not possible, we approximate the integral with an empirical sum and use the data vectors u1, ...,uN for this purpose.
The copula function FU is not known for real data. Instead, we use the empirical cdf of the unit space data vectors in (9), i.e., for w ∈ [0, 1]D we set
F̂U(w) = 1
N N∑ n=1 D∏ d=1 1[und ≤ wd]. (10)
Due to the curse of dimensionality this approximation of the joint cdf will have a coarse structure in high dimensions. This is why we only use the values at the data vectors u1, ...,uN in (9).
For parametric copula models the distribution function FV is analytic. For the non-parametric models, including IGC, we again resort to approximation. We draw 105 samples from the copula distribution model and use their empirical cdf in unit space in (9). We found that this procedure provided reliable and stable estimates for the underlying copula in dimensions D ≤ 5.

heading: 4.1 Synthetic data
text: Learning bivariate parametric copulas We first show that IGC models are able to emulate bivariate parametric copulas. We compare our IGC approach to two non-parametric copula density estimation techniques, kernel density estimation with beta kernels (BETA) [6] and the transformation local likelihood estimator (TLL) [12]. These are state-of-the art non-parametric estimators for paircopulas [30]. The compact support of the beta distribution on [0, 1] is convenient for copula density estimation. The idea of the TLL approach is to first transform the data using the inverse normal CDF such that the data is supported on R2. Then the density is estimated via local regression using linear (TLL1) or cubic polynomials (TLL2) on a fixed grid of points. Finally the estimated density is transformed back to unit space. Additionally we report the performance of a parametric model. In this approach a copula family is selected from a set of copulas based on the Bayesian information criterion (BIC). This set comprises the Independence, Gaussian, Student-t, Clayton, Gumbel, Frank, Joe, BB1, and BB7 copula. Parameters are then estimated via maximum likelihood.
We run experiments for data sets generated from a Student-t, a Gumbel, and a Clayton copula as well as the copula resulting from a two-component Gaussian mixture distribution as shown in Figure 2. For each one of these, we sampled a random parameter set and then generated 1000 training samples from the resulting model. We repeated each experiment 25 times. Appendix B contains the details on the considered parameter ranges and the exact sampling procedure.
Figure 2 shows the data and the model with intermediate steps for one instance of the Gaussian mixture case. The fitted copula matches the data distribution in unit space very well. Figure 3 presents aggregated ISE results for the different setups. The IGC model shows comparable performance scores to TLL1 and TLL2 for the Student-t, Clayton, and Gumbel copulas, and better performance than the BETA approach. For the more complex Gaussian mixture test case, IGC shows superior performance and a much lower variance than all baseline methods. Notably, the parametric approach with BIC-based model selection results in large variance in accuracy for the Clayton and Gumbel data, most likely because the wrong copula family is selected in some cases.
Learning multivariate vine copulas We now turn to the problem of estimating copulas in dimensions D > 2. To this end, we conduct a similar experiment as before using 5-dimensional vine copulas as data generating distribution. For each of 25 repetitions, we first sample a random tree structure using pyvinecopulib’s RVineStructure.Simulate() method. Then, we randomly assign parametric pair-copulas with random parameters to each edge in the tree. The considered families of bivariate copulas are Independence, Gaussian, Student-t, Clayton, Gumbel, Frank, Joe, BB1, and BB7. Further details are available in Appendix B. We generate 5000 samples from the resulting vine copula
model and use these to train our IGC model as well as the benchmark models. As benchmarks, we use a vine copula model with TLL2 pair-copulas only and a vine copula model which can select all parametric copula families named above as well as the TLL2. Additionally, we estimate a parametric Gaussian copula. The parameters are estimated via maximum likelihood and the copula families are selected using the BIC. The selection of the vine structure is based on the Dissmann algorithm [9]. We evaluate all models with the ISE at 10000 data points sampled from the true model.
The results of the simulation study are presented in Figure 4. The IGC model has the lowest mean ISE (0.0697) followed by the TLL2-Vine model (0.2093), while the latter has a lower median ISE (0.0339 vs. 0.045). This is the result of some large outliers of the ISE for the vine models that do not occur for the IGC model. These are most likely caused by a poorly selected vine structure. Interestingly, the results specifically imply that the data is on average better modeled by the IGC model than by the vine approach although this model class entails the true data generating model. However, recovering the true model does not seem to be an easy task.
4.2 Real data
Exchange rates Copulas are a popular tool in financial risk management as they can be used to estimate the distribution of the multivariate returns over different assets [34] under the assumption of a stationary copula. We consider a data set of size N = 5844 that contains 15 years of daily exchange rates between the US-Dollar and the Canadian Dollar, the Euro, the British Pound, the Swiss Franc, and the Japanese Yen. The data was obtained from the R package qrm_data. We preprocess the data to filter out the effects of temporal dependencies. To this end, we fit an AR(1)-GARCH(1,1) [4] process with Student-t innovations to the time series of the daily returns. We then obtain the standardized residuals from the AR-GARCH models and transform these observations to the unit space using the empirical cdfs. Figure 5a shows a kernel density estimate for two selected dimensions of the resulting data set in unit space, namely the US-Dollar/Euro and US-Dollar/Pound exchange rates. While being non-trivial and multi-modal, the distribution is strongly concentrated along the diagonal.
We then estimate a Gaussian copula, a vine copula, a vine copula with only TLL2 pair-copulas, a GMMN as proposed by [16], and our IGC model for this data. The GMMN model has the same loss function, architecture, and hyper parameters as the IGC model but uses a sigmoid activation in the final layer. To ensure a test data set of appropriate size, we use a 5-fold cross-validation scheme where 20% of the data is used for training and 80% for testing. For the IGC and the GMMN model we additionally use five different random initializations per fold for the neural network weights, i.e., we report means and standard deviations from 25 values for the IGC and the GMMN model and 5 values for the other methods. The results are presented in Table 1. The IGC model clearly outperforms the Gaussian baseline and also the GMMN model. However, the vine copula models show lower average ISEs. This might be due to the symmetric nature of the data which can be approximated well by Student-t and TLL pair-copulas.

heading: MAGIC Gamma Telescopes data In
text: order to test our approach on a more complex dependency structure, we consider the MAGIC (Major Atmospheric Gamma-ray Imaging Cherenkov) Telescopes data set available from the UCI repository (https://archive.ics.uci.edu/ml/ datasets/MAGIC+Gamma+Telescope). This data was also used for benchmarking non-parametric copula estimation techniques in [32]. We only consider the observations classified as gamma and the 5 variables fLength, Width, fConc, fM3Long,
fM3Trans. The size of the data set is N = 12332. We use the empirical cdfs to transform the
observations to the unit space. Figure 5b showcases the dependency structure of the data for two selected data dimensions. The observed structures are highly asymmetric. We again fit the data with a Gaussian copula, a vine copula with all available pair-copulas as above, a vine copula with TLL2 pair-copulas, and a GMMN model with a sigmoid output layer [16] as benchmarks and use the same 5-fold CV strategy as before, i.e., in each fold we use 20% as training data, 80% as test data, and report the average ISE. The results are presented in Table 1. The IGC model achieves the lowest average ISE. Again the GMMN model shows a substantially larger mean ISE with a much larger standard deviation. As only TLL2 pair-copulas were selected by the BIC criterion, the scores for both vine models are the same. The Gaussian copula is clearly not suited for such complex data as the average ISE is 5 times larger than for the other approaches.
Copula Autoencoders [40] introduced the Vine Copula Autoencoder, a generative model which uses vine copulas for ex-post density estimation of the latent space of a trained autoencoder. In a first step, an autoencoder is trained on a data set to learn a low dimensional representation of the data. After training, the encoder part of the network is used to map the training data to the autoencoder’s latent space representation. Next, the uni-variate marginal distributions are estimated, e.g., by using the empirical cdfs or kernel density estimation, and the compressed data is mapped to unit space. After fitting a copula model to the observations in unit space, one can sample from the fitted copula model, apply the inverse marginal cdfs, and map the simulated data back to the image space using the decoder network of the autoencoder in order to generate new data samples.
In the following, we present results for the FashionMNIST [44] data set. We train a convolutional autoencoder on the entire training data of 60000 samples with a latent space of dimension 25. Details on the architecture are found in Appendix B. After training the autoencoder, we estimate the empirical cdfs of the compressed data using the training data. See Figures 5c and 5d for exemplary visualizations of the resulting bivariate data densities in unit space. Notably, these distributions are more complex than the ones from the previous experiments. We fit this data with a Gaussian copula, a vine copula with TLL2 pair-copulas, a GMMN with sigmoid output layer
[16], and an IGC model. We also test an independence copula, i.e., we assume independence over the latent space. Additionally we report results for a standard variational autoencoder (VAE) [21] with the same architecture.
Exemplary samples from all models and the test set are presented in Figure 6. Sampling with the independence copula, i.e., assuming no dependency structure, leads to images with many artifacts. The Gaussian copula produces better images, but still produces some artifacts. Samples from the vine copula or the IGC model are comparable in quality. They show smaller details and few implausible artifacts. The images generated by the VAE are blurry and show much less variation in pixel intensity than the test data.
In 25 dimensions it is not feasible anymore to compute and compare empirical cdfs as required for ISE scoring. We thus resort to the MMD [15] for the numerical evaluation since it is commonly used as a simple and robust measure for evaluating image generation models [46]. We generate 10000 images from each model and compare those to the 10000 test set images by computing the MMD with a Gaussian kernel. We use the bandwidths 1.5, 350, and 8 for the latent unit space, the latent data space, and the image space, respectively. The bandwidths were selected based on the median heuristic proposed in [15].
The results are found in Table 2. We provide p-values for these results using the test proposed by [5] in Appendix B. The IGC and GMMN models perform better than all other models for the latent unit space distribution as well as the image space. Interestingly the GMMN model shows a relatively large MMD value for the latent data space. This could be caused by deviations from marginal uniformity of the learned copula distribution which alter the marginal distributions in the data space.

heading: 5 Conclusion
text: We have proposed the first fully non-parametric model framework for copulas in higher dimensions that guarantees uniformity of the marginal distributions. The proposed IGC approach, which is based on an implicit generative step and a differentiable ranking transformation, is structurally simple. Yet, if the generator class is sufficiently complex, any copula dependency structure can be modeled. The model can be well implemented with standard deep learning frameworks. For various data sets we have shown a modeling performance on par or above other state-of-the-art approaches, especially, the vine copula approach.
IGC models should be further investigated in various ways. First, it would be straight-forward to condition the model on external factors, either by including such values into the input of the generator network or by reparameterization of the noise distributions. This would allow to describe contextdependent changes of the dependency structure. Second, many more applications for copulas should be examined, e.g., synthetic tabular data generation [45]. Third, the ranking operation could probably be sped up from O(M2) to O(M logM) using ideas from [3]. Finally, the use of adversarial training schemes could also be investigated.

heading: Acknowledgments and Disclosure of Funding
text: This work has been performed in the context of the LOEWE center emergenCITY.

Review summaries:
The paper proposes an implicit method (IGC) of estimating a copula. The suggested approach is elegant and from the presented results the method seems to work on par with existing non-parametric vine copula approaches. Previous works have made an attempt to use GANs for implicit copula estimation with difficulties, however, in IGC the authors suggest techniques to address the requirements for uniform margins with a simple two-step procedure leveraging the probability integral transform. The general idea and contributions are clearly presented and I find this work valuable due to its simplicity and I believe its efficiency can be further exploited in other copula-related/inspired applications.

The authors propose to model general copula models by training a neural network that samples from a latent space and maps samples to the high dimensional unit cube. Fitting the model is done by transforming the generated samples to have uniform margins by empirically computing the ranks. Then, the model is optimized by minimizing the MMD distance w.r.t true uniform data samples.