# Generative Archimedean Copulas

Yuting Ng Department of Electrical and Computer Engineering
Duke University  Ali Hasan  Department of Biomedical Engineering
Duke University  Khalil Elkhalil  Department of Electrical and Computer
Engineering
Duke University  Vahid Tarokh  Department of Electrical and Computer
Engineering
Duke University

###### Abstract

We propose a new generative modeling technique for learning multidimensional
cumulative distribution functions (CDFs) in the form of copulas. Specifically,
we consider certain classes of copulas known as Archimedean and hierarchical
Archimedean copulas, popular for their parsimonious representation and ability
to model different tail dependencies. We consider their representation as
mixture models with Laplace transforms of latent random variables from
generative neural networks. This alternative representation allows for
computational efficiencies and easy sampling, especially in high dimensions.
We describe multiple methods for optimizing the network parameters. Finally,
we present empirical results that demonstrate the efficacy of our proposed
method in learning multidimensional CDFs and its computational efficiency
compared to existing methods.

##  1 Introduction

Copulas are a special class of cumulative distribution functions (CDFs) that
model the dependencies between multiple random variables in isolation from
their marginals [Nelsen, 2010, Joe, 2014]. Copulas have found applications in
many areas, including hydrology [Genest and Favre, 2007] and finance
[Cherubini et al., 2004]. In finance, for example, more expressive modeling
using copulas of the joint distribution of two stocks can result in more pairs
trading opportunities [Stander et al., 2013, Liew and Wu, 2013].

In machine learning, copulas have been used to create new distributions,
increasing the flexibility of modeling multivariate dependencies [Wilson and
Ghahramani, 2010, Elidan, 2010, Huang and Frey, 2011, Tagasovska et al., 2019,
Wiese et al., 2019, Kamthe et al., 2021, Chilinski and Silva, 2020]. The
utility of copulas can be attributed to their powerful representation
capabilities, ease of use and intuitive decomposition into marginals and a
dependence function. However, many challenges related to parameterization and
estimation are still unsolved.

A particularly useful class of copulas are known as Archimedean copulas, which
endow a specific structure for representing the dependence function in terms
of a one-dimensional _generator_ function. Most work involving Archimedean
copulas consider different parameterizations for this generator function.
Parameterizations of the generator function have generally been limited to
simple forms, since complicated generator functions lead to difficulties in
computing the copula density, a necessary component for maximum likelihood
estimation. Ling et al. [2020] proposed parameterizing the generator function
as a neural network, but ran into computational difficulties for dimensions
greater than 5. We consider an alternative construction based on a mixture
representation with latent random variables, first proposed by Marshall and
Olkin [1988], wherein we parameterize a latent distribution, whose Laplace
transform acts as the generator function, with a generative neural network.
Depending on the application, this latent variable is sometimes known as a
_resilience_ or _frailty_ parameter. Using this construction, we can scale
computations to higher dimensions and bypass numerical issues involved with
automatic differentiation.

Employing the Laplace transform to a learned latent model also provides
important benefits beyond computational efficiency and numerical stability.
When sampling from the copula using established approaches [Marshall and
Olkin, 1988, McNeil, 2008, Hering et al., 2010], knowledge of the latent
distribution is necessary. Parameterizing the latent distribution with a
generative neural network allows for efficient sampling after training.

Archimedean copulas can also be extended to the so-called hierarchical (or
nested) Archimedean copulas, where multiple generators are used in conjunction
to increase expressiveness of the model [Joe, 1997]. This architecture
mitigates a central deficiency of vanilla Archimedean copulas — the assumed
symmetry in the dependence structure. We use a construction based on Lévy
subordinators, i.e. non-decreasing Lévy processes such as the compound Poisson
process, first proposed by Hering et al. [2010], and parameterize the Lévy
subordinators using generative neural networks. We also use Laplace transforms
as in the vanilla Archimedean copula to obtain the generator functions and
subsequently recover a richer class of copulas.

### Related Work

Our part of work on Archimedean copulas is related to [Ling et al., 2020],
where a neural network is proposed to represent the generator function of an
Archimedean copula. We propose instead a generative neural network to
represent the latent random variable, whose Laplace transform gives the
generator function of the Archimedean copula. We then approximate the Laplace
transform with the empirical Laplace transform using samples from the
generative neural network. We note that there exist prior work that replace
the Laplace transform with the empirical Laplace transform, such as those on
the estimation of compound Poisson processes and distribution goodness-of-fit
tests. These can be found in [Csörgő and Teugels, 1990, Henze et al., 2012],
but do not consider/employ neural networks.

Existing semiparametric methods for Archimedean copulas are mainly
concentrated on two dimensional cases, and their efficacy in higher dimensions
remains unclear [Hernández-Lobato and Suárez, 2011, Hoyos-Argüelles and Nieto-
Barajas, 2020]. Other work on the mixture representation with a latent random
variable is limited to cases of known distributions that can be sampled and
for which the Laplace transform can be calculated, since it is often
challenging to find and sample from a distribution corresponding to arbitrary
Laplace transforms [McNeil, 2008, Hofert, 2008].

Our part of work on hierachical Archimedean copulas is inspired by [Hering et
al., 2010] who recognized that sufficient nesting conditions of hierarchical
Archimedean copulas may be satisfied using Lévy subordinators. We then let the
increments associated with the Lévy measure of the Lévy subordinator be the
output of a generative neural network, and compute its integral in the Laplace
exponent as an expectation with samples from the generative neural network.
Related work parameterizing the Lévy measure with a neural network can be
found in [Xu and Darve, 2020], but the integral is approximated as a Riemann
sum, and it does not relate to hierarchical Archimedean copulas.

Other related works combine one-parameter families of Archimedean copulas,
usually in a homogeneous manner, where all components are from the same
family. It is challenging to combine Archimedean copulas from different
families due to the nesting conditions. For example, the Clayton and Gumbel
copulas are not compatible for nesting [McNeil, 2008]. Thus, related works on
heterogeneous Archimedean copulas have resulted in limited combinations of
different families [McNeil, 2008, Hofert, 2008, Savu and Trede, 2010, Okhrin
et al., 2013, Górecki et al., 2017].

### Main Contributions

First, we propose to use a generative neural network to represent the latent
random variable, whose Laplace transform provides the generator function of an
Archimedean copula. This allows approximation of the Laplace transform with
its empirical version through samples from the generative neural network.
Computing higher-order derivatives using the properties of the empirical
Laplace transform additionally allows scalability to higher-dimensional data.
Second, we extend this concept to modeling hierarchical Archimedean copulas
with Lévy subordinators. We represent the Lévy measure of a Lévy subordinator
with a generative neural network and compute its Laplace exponent using
samples from the generative neural network. We then propose three methods for
training: maximum likelihood with the copula density, goodness-of-fit with the
Cramér-von Mises statistic, and adversarial training by minimizing a
divergence between true samples from data and fake samples from the copula.
Finally, we adapt existing Marshall-Olkin type efficient sampling algorithms
to our parameterization with generative neural networks. The source code for
this paper may be found at https://github.com/yutingng/gen-AC.

### Outline

Section 2 provides the mathematical background on copulas, Archimedean copulas
and hierarchical Archimedean copulas. Section 3 discusses modeling, sampling
and training generative Archimedean copulas. Section 4 extends the
construction to hierarchical Archimedean copulas. Section 5 shows our
experiment results on learning known Archimedean and hierarchical Archimedean
copulas that have different tail dependencies. We also compare its flexibility
in fitting real-world data to commonly-used one-parameter families. In
addition, we show its computational efficiency and sampling in higher-
dimensions. Finally, we conclude the paper in Section 6.

##  2 Background

We begin by describing the necessary background on copulas. A _copula_ is a
multivariate cumulative distribution function (CDF) where all univariate
margins are uniform, i.e. it is the CDF of a vector of dependent uniform
random variables. Multidimensional dependence modeling with copulas is based
on a theorem due to Sklar [1959] which gives a general representation of a
multivariate CDF as a composition of its univariate margins and a copula.

######  Theorem 1 (Sklar’s theorem).

For a $d-$variate cumulative distribution function $F$, with $j$th univariate
margin $F_{j}$, and $j^{th}$ quantile function $F_{j}^{-1}$, the copula
associated with $F$ is a cumulative distribution function
$C:[0,1]^{d}\to[0,1]$ with $\text{U}(0,1)$ margins satisfying:

$\displaystyle F(\mathbf{x})$
$\displaystyle=C(F_{1}(x_{1}),\cdots,F_{d}(x_{d})),\;\mathbf{x}\in\mathbb{R}^{d},$
(1) $\displaystyle C(\mathbf{u})$
$\displaystyle=F(F_{1}^{-1}(u_{1}),\cdots,F_{d}^{-1}(u_{d})),\;\mathbf{u}\in[0,1]^{d}.$
(2)

In addition, if $F$ is continuous, then $C$ is unique.

Moreover, due to Sklar’s theorem, every CDF endows such a decomposition. Thus,
copulas allow characterization of the multivariate dependence between the
random variables $X_{1},\cdots,X_{d}$ separately from their univariate margins
$F_{1},\cdots,F_{d}$ [Nelsen, 2010, Joe, 2014].

###  2.1 Archimedean Copulas

An important class of copulas are the Archimedean copulas, due to their ease
of construction and ability to represent different tail dependencies. An
Archimedean copula is defined as:

$C(\mathbf{u})=\varphi\left(\varphi^{-1}(u_{1})+\cdots+\varphi^{-1}(u_{d})\right),$
(3)

with density:

$\displaystyle c(\mathbf{u})$
$\displaystyle=\frac{\partial^{d}C(u_{1},\cdots,u_{d})}{\partial
u_{1}\cdots\partial u_{d}}$ (4)
$\displaystyle=\frac{\varphi^{(d)}(\varphi^{-1}(u_{1})+\cdots+\varphi^{-1}(u_{d}))}{\prod_{i=1}^{d}\varphi^{\prime}(\varphi^{-1}(u_{i}))}.$

For the above expression to be a valid copula for all $d$, the one-dimensional
function $\varphi:[0,\infty)\to[0,1]$, known as the _generator_ of the
Archimedean copula must satisfy:

  * •

$\varphi(0)=1,\;\varphi(\infty)=0$,

  * •

$\varphi$ is _completely monotone_ ,
i.e. $(-1)^{k}\varphi^{(k)}\geq 0$ for all $k\in\\{0,1,2,\cdots\\}$.

The criteria that $\varphi$ is completely monotone, i.e. its derivatives
change signs, guarantees positiveness of the copula density [Kimberling,
1974]. The class of completely monotone $\varphi$ coincides with the class of
Laplace-Stieltjes transforms (henceforth simply Laplace transforms) of a
positive random variable [Bernstein, 1929, Widder, 1941].

######  Theorem 2 (Bernstein [1929] and Widder [1941]).

$\varphi$ is completely monotone and $\varphi(0)=1$ if and only if $\varphi$
is the Laplace transform of a positive random variable,

$\varphi(x)=\int_{0}^{\infty}e^{-xs}dF_{M}(s),$ (5)

where $M>0$ is a positive random variable with Laplace transform $\varphi$.

Conversely, a probabilistic construction of the Archimedean copula as a
mixture model, with the variables being conditionally independent given a
positive latent random variable, leads to the Laplace transform representation
for $\varphi$. For a given $d$, $\varphi$ may come from a broader class of
functions than Laplace transforms [McNeil and Nešlehová, 2009]. However, if
$\varphi$ is not a Laplace transform, the simple mixture representation fails
[Marshall and Olkin, 1988]. In the mixture representation, the latent
variable, depending on its application, is known as a _resilience_ or
_frailty_ parameter [Marshall and Olkin, 1988, Joe, 1997]. Common Archimedean
copulas such as the Ali-Mikhail-Haq, Clayton, Frank, Gumbel and Joe copulas
can be respectively derived from geometric, gamma, logarithmic, stable, and
Sibuya latent distributions. The mixture representation also leads to
efficient sampling algorithms [Marshall and Olkin, 1988, McNeil, 2008].

We restate the probabilistic construction and sampling algorithm in the
supplementary material.

###  2.2 Hierarchical Archimedean Copulas

While Archimedean copulas have been widely employed, the functional symmetry
of the Archimedean copula implies exchangeability of the underlying dependence
structure, which is sometimes not realistic. Hierarchical (or nested)
Archimedean copulas are popular for overcoming this drawback [Joe, 1997].

Figure 1: Hierarchical Archimedean copula with $J$ nested, possibly
hierarchical, Archimedean copulas.

In this case, the copula can be written as:

$C(\mathbf{u})=C_{\varphi_{0}}(C_{\varphi_{1}}(\mathbf{u}_{1}),\cdots,C_{\varphi_{J}}(\mathbf{u}_{J})),$
(6)

where $C_{\varphi_{j}}$ are nested, possibly hierarchical, Archimedean copulas
with generators $\varphi_{j}$, inputs
$\mathbf{u}_{j}=[u_{j,1},\cdots,u_{j,d_{j}}]$, $j\in\\{1,\cdots,J\\}$, and
$\mathbf{u}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{J}]$, $d_{1}+\cdots+d_{J}=d$.

For the above expression to be a valid copula, additional sufficient nesting
conditions, derived from the nested mixture representation, first given in
[Joe, 1997] and restated in [McNeil, 2008] for nesting to arbitrary depth are:

  * •

$\varphi_{j}$ for all $j\in\\{0,1,\cdots,J\\}$ are completely monotone,

  * •

$(\varphi_{0}^{-1}\circ\varphi_{j})^{\prime}$ for $j\in\\{1,\cdots,J\\}$ are
completely monotone.

The criteria that $(\varphi_{0}^{-1}\circ\varphi_{j})^{\prime}$ are completely
monotone come from the composition of an _outer generator_ $\varphi_{0}$ and
an _inner generator_ $\varphi_{j}$ to produce a completely monotone Laplace
transform _nested generator_ of the form
$e^{-M\varphi_{0}^{-1}\circ\varphi_{j}}$, where $M$ is distributed with
Laplace transform $\varphi_{0}$ [Joe, 1997, McNeil, 2008]. This criteria is
addressed in [Hering et al., 2010] using Lévy subordinators, i.e. non-
decreasing Lévy processes such as the compound Poisson process, by recognizing
that the Laplace transform of Lévy subordinators at a given ‘time’ $t\geq 0$
have the form $e^{-t\psi_{j}}$, where the Laplace exponent $\psi_{j}$ has
completely monotone derivative. Conversely, a probabilistic construction by
combining Lévy subordinators evaluated at common ‘time’ $t=M$, leads to a
well-defined hierarchical Archimedean copula with an efficient sampling
algorithm [Hering et al., 2010]. We restate the probabilistic construction
from [Hering et al., 2010] in the supplementary material.

Thus for a given _outer generator_ $\varphi_{0}$, a _compatible inner
generator_ $\varphi_{j}$ can be modeled as a composition of the outer
generator and the Laplace exponent $\psi_{j}$ of a Lévy subordinator:

$\varphi_{j}(x)=(\varphi_{0}\circ\psi_{j})(x),$ (7)

where the Laplace exponent $\psi_{j}:[0,\infty)\to[0,\infty)$ of a Lévy
subordinator has a convenient representation with drift $\mu_{j}\geq 0$ and
Lévy measure $\nu_{j}$ on $(0,\infty)$ due to the Lévy-Khintchine theorem
[Sato, 1999]:

$\psi_{j}(x)=\mu_{j}x+\int_{0}^{\infty}(1-e^{-xs})\nu_{j}(ds).$ (8)

A popular Lévy subordinator is the compound Poisson process with drift
$\mu_{j}\geq 0$, jump intensity $\beta_{j}>0$ and jump size distribution
determined by its Laplace transform $\varphi_{M_{j}}$. In this case, the
Laplace exponent has the following expression:

$\displaystyle\psi_{j}(x)$
$\displaystyle=\mu_{j}x+\beta_{j}(1-\varphi_{M_{j}}(x))$ (9)
$\displaystyle=\mu_{j}x+\beta_{j}(1-\int_{0}^{\infty}e^{-xs}dF_{M_{j}}(s)),$

where $M_{j}>0$ is a positive random variable with Laplace transform
$\varphi_{M_{j}}$ characterizing the jump sizes of the compound Poisson
process.

In addition, we choose $\mu_{j}>0$ to satisfy the condition
$\varphi_{j}(\infty)=(\varphi_{0}\circ\psi_{j})(\infty)=0$ such that
$\varphi_{j}$ is a valid generator of an Archimedean copula.

##  3 Generative Archimedean Copulas

Motivated by the probabilistic construction of the Archimedean copula, we
propose to learn the distribution of the positive latent variable by
approximating its Laplace transform using samples from a generative neural
network.

###  3.1 Modeling the Latent Variable with a Generative Neural Network

We let $M$ be the output of a generative neural network such that samples
$M\sim F_{M}$ are computed as $M=G(\epsilon;\theta)$, where $G(.;\theta)$
represents the generative neural network with parameters $\theta$ and
$\epsilon$ is a source of randomness. Unlike the modeling of monotone
functions with neural networks [Chilinski and Silva, 2020], there is no
restriction on the weights and intermediate activations of $G(.;\theta)$. In
this preliminary work, the network architecture is a multilayer perceptron. To
guarantee that $M$ is a positive random variable, we use $\exp(.)$ as the
output activation.

We then approximate the Laplace transform with its empirical version using $L$
samples of $M$ from $G(.;\theta)$ as:

$\varphi(x)=\int_{0}^{\infty}e^{-xs}dF_{M}(s)=\mathbb{E}_{M}[e^{-Mx}]\approx\frac{1}{L}\sum_{l=1}^{L}e^{-M_{l}x}.$
(10)

Derivatives of the Laplace transform are similarly approximated with their
empirical version as:

$\varphi^{(k)}(x)=\mathbb{E}_{M}[(-M)^{k}e^{-Mx}]\approx\frac{1}{L}\sum_{l=1}^{L}(-M_{l})^{k}e^{-M_{l}x}.$
(11)

Subsequently, we replace instances of $\varphi$ and $\varphi^{(k)}$ in the
copula distribution, density and sampling algorithm with their sample
approximations computed as in (10) and (11).

###  3.2 Generating Samples from the Archimedean Copula

We modify existing Marshall-Olkin type sampling algorithms [Marshall and
Olkin, 1988, McNeil, 2008] to our parameterization with generative neural
networks, as detailed in Algorithm 1 and Figure 2, on the next page.

This sampling method is efficient as it only requires sampling unit
exponential random variables $E_{j}\sim\text{Exp}(1),\;j\in\\{1,\cdots,d\\}$
and a latent random variable $M=G(\epsilon;\theta)$. In addition, unlike the
conditional sampling method, this sampling method does not require
differentiation of the copula distribution to get the conditional distribution
and does not require inversion of the conditional distribution.

Algorithm 1 Sampling Generative Archimedean Copulas

Input: $G(.;\theta)$,

1: Sample $M$ as $M=G(\epsilon;\theta).$

2: Sample i.i.d. $E_{j}\sim\text{Exp}(1),j\in\\{1,\cdots,d\\}$.

3: Approximate $\varphi$ with samples $\\{M_{l}\\}_{l=1}^{L}$, where
$M_{l}=G(\epsilon_{l};\theta)$, as in (10).

4: Compute $\mathbf{U}$ where $U_{j}=\varphi(E_{j}/M),j\in\\{1,\cdots,d\\}$.

Output: $\mathbf{U}$.

Figure 2: Sampling Archimedean copulas using generative neural networks and
Marshall-Olkin type sampling.

###  3.3 Training Methods

An important consideration when modeling CDFs is the optimization procedure
for fitting the model to data. We describe multiple methods for fitting the
model to data with various performance and efficiency trade-offs.

####  3.3.1 Training with Maximum Likelihood

We consider training through maximum likelihood by minimizing the negative log
likelihood with backpropagation gradient descent on the model parameters,
similar to the proposal in [Ling et al., 2020]. However, since the copula
models the CDF, differentiation is required to obtain the copula density.
Unlike [Ling et al., 2020] that computes the copula density from the copula
distribution using automatic differentiation, we compute the copula density
from its analytical expression in (4) using the properties of the Laplace
transform for computing higher-order derivatives in (11). For increasing
dimensions, computing higher-order derivatives using the Laplace transform
representation instead of automatic differentiation leads to a significant
speed up in computation.

For the computation of $\varphi^{-1}$ and its derivative with respect to model
parameters, we borrow the method in [Ling et al., 2020]. The inverse is
computed using Newton’s root-finding method. The derivatives are computed from
the derivatives of $\varphi$ then supplemented to backpropagation.

####  3.3.2 Training with Goodness-of-Fit

To circumvent computing the copula density, the model may also be fitted to
data via minimum distance criterions used in goodness-of-fit tests [Genest et
al., 2009]. Though not statistically efficient compared to maximum likelihood
estimation, minimum distance estimation is significantly less computationally
intensive.

We consider the Cramér-von Mises statistic [Cramér, 1928] to measure a
discrepancy between the model copula $C_{\theta}$ and the empirical copula
$C_{N}$:

$S_{N}=\frac{1}{N}\sum_{i=1}^{N}(C_{\theta}(\mathbf{u}_{i})-C_{N}(\mathbf{u}_{i}))^{2},$
(12)

where $\mathbf{u}_{i}$ is an observation of the margins, $N$ is the number of
observations and $C_{N}$ is the empirical copula given by:

$C_{N}(\mathbf{u})=\frac{1}{N}\sum_{i=1}^{N}\mathbbm{1}\\{u_{i,1}\leq
u_{1},\cdots,u_{i,d}\leq u_{d}\\}.$ (13)

####  3.3.3 Adversarial Training with Samples

An alternative way to train the model is by minimizing a divergence between
true samples from data and fake samples from the model copula, similar to
generative adversarial networks (GANs) [Goodfellow et al., 2014]. In this
case, we solve the minimax problem in GANs where the generating network must
satisfy an Archimedean copula. This is another method that allows training
without computing the copula density.

We create a discriminative neural network $D(.;\phi)$ with parameters $\phi$
and $\text{sigmoid}(.)$ output activation to distinguish between true samples
from data and fake samples from the copula. We then minimize the Jensen-
Shannon loss between true and fake samples as in [Goodfellow et al., 2014]:

$\min_{\theta}\max_{\phi}\mathbf{E}_{U\sim\text{data}}[\log(D(U;\phi))]+\mathbf{E}_{\tilde{U}\sim
C}[\log(1-D(\tilde{U};\phi))],$ (14)

where $\tilde{U}\sim C$ is generated via the sampling method described in
Algorithm 1 using the latent random variable represented as the output of the
generative neural network $G(.;\theta)$ with parameters $\theta$ as discussed
in Section 3.1.

##  4 Generative Hierarchical Archimedean Copulas

In the following, we extend the application of generative neural networks to
hierarchical Archimedean copulas. We present our results for two levels of
hierarchy, but our construction extends to nesting with more levels.

###  4.1 Modeling the Laplace Exponent with a Generative Neural Network

For a given outer generator $\varphi_{0}$, the inner generator
$\varphi_{j},j\in\\{1,\cdots,J\\}$ is obtained as the composition
$\varphi_{j}=\varphi_{0}\circ\psi_{j}$, as in (7), where $\psi_{j}$ is the
Laplace exponent of a compound Poisson process with Lévy-Khintchine
representation, as in (9). We let the drift $\mu_{j}>0$ and the jump intensity
$\beta_{j}>0$ be trainable parameters with $\exp(.)$ output activation. We let
the jump size $M_{j}>0$ be the output of a generative neural network
$G(.;\theta_{j})$ with parameters $\theta_{j}$ and $\exp(.)$ output
activation. We then compute the Laplace transform $\varphi_{M_{j}}$ and its
derivatives $\varphi_{M_{j}}^{(k)}$ using samples from $G(.;\theta_{j})$ as in
(10) and (11).

###  4.2 Generating Samples from the Hierarchical Archimedean Copula

We modify the Marshall-Olkin type algorithm given in [Hering et al., 2010] to
work with our parameterization using generative neural networks. We first
describe sampling of a compound Poisson process in Algorithm 2. We then
describe sampling of a generative hierarchical Archimedean copula in Algorithm
3. A sample from the hierarchical Archimedean copula is obtained by combining
compound Poisson processes evaluated at a common ‘time’ $t=M$, where $M$ is
the random variable with distribution given by the Laplace transform outer
generator $\varphi_{0}$.

Algorithm 2 Sampling compound Poisson process with jump sizes parameterized by
generative neural network

Input: $t,\mu_{j},\beta_{j},G(.;\theta_{j})$,

1: Sample $N_{j}(t)\sim\text{Pois}(\beta_{j}t)$, i.e. the number of jumps by
time $t$ of a Poisson random variable with rate $\beta_{j}$.

2: Sample $N_{j}(t)$ samples of $M_{j}$ from $G(.;\theta_{j})$.

3: Compute $\Lambda_{j}(t)=\mu_{j}t+\sum_{i=1}^{N_{j}(t)}M_{j,i}$.

Output: $\Lambda_{j}(t)$.

Algorithm 3 Sampling Generative Hierarchical Archimedean Copulas

Input: $G(.;\theta_{0}),\\{\mu_{j},\beta_{j},G(.;\theta_{j})\\}_{j=1}^{J}$,

1: Sample $t=M$ from $G(.;\theta_{0})$.

2: Approximate $\varphi_{0}$ with samples from $G(.;\theta_{0})$, as in (10).

3: for $j\in\\{1,\cdots,J\\}$ do

4: Sample $\Lambda_{j}(M)$, the compound Poisson process at $t=M$, following
Algorithm 2.

5: Approximate $\psi_{j}$ with samples from $G(.;\theta_{j})$, as in (10).

6: end for

7: Sample i.i.d. $E_{j,i}\sim\text{Exp}(1)$, $j\in\\{1,\cdots,J\\}$,
$i\in\\{1,\cdots,d_{j}\\}$.

8: Compute $\mathbf{U}=(U_{1,1},\cdots,U_{J,d_{J}})$ as
$U_{j,i}=(\varphi_{0}\circ\psi_{j})(E_{j,i}/\Lambda_{j}(M)),j\in\\{1,\cdots,J\\},i\in\\{1,\cdots,d_{J}\\}$.

Output: $\mathbf{U}$.

###  4.3 Training with Goodness-of-Fit and Maximum Likelihood

We first fit the outer generator $\varphi_{0}$, fix it, then fit the inner
generators $\varphi_{j}=\varphi_{0}\circ\psi_{j}$. Fixing the outer generator
then optimizing the inner generator provides additional numerical stability
during training. In our experiments, the outer generator was trained using
minimium distance estimation with the empirical copula based Cramér-von Mises
statistic in (12) and empirical copulas on $C_{\varphi_{j}}$. An alternative
method may be to train the outer generator using a composite likelihood with
bivariate margins since bivariate margins are Archimedean with generator given
by the outer generator. The inner generators were trained using maximum
likelihood estimation with copula densities $c_{\varphi_{j}}$ in (4).

##  5 Experiments

###  5.1 Generative Archimedean Copula

####  5.1.1 Learning Bivariate Copulas with Different Tail Dependencies and
Fitting Real-World Data

Following the experiment setup in [Ling et al., 2020], we consider the
Clayton, Frank, and Joe copulas, chosen for their different tail dependencies,
and the following real-world data sets: Boston housing, Intel-Microsoft (INTC-
MSFT) stocks and Google-Facebook (GOOG-FB) stocks. We applied the three
training methods discussed earlier: maximum likelihood, goodness-of-fit and
adversarial training. All training methods were implemented in PyTorch and
converged within 10k epochs. Experiment details are given in the supplementary
material.

The negative log-likelihoods from learning known copulas are reported in Table
2. We use the following shorthands ‘GT’, ‘ACNet’, ‘MLE’, ‘CvM’, ’GAN’ to
respectively denote ground truth, ACNet [Ling et al., 2020], and generative
Archimedean copulas trained with maximum likelihood, goodness-of-fit and
adversarial training. The negative log-likelihoods from fitting real-world
data are reported in Table 2, where the log-likelihood of the best-fit single
parameter copula (chosen from Clayton, Frank, Joe and Gumbel, as in [Ling et
al., 2020]), with shorthand ‘BF’ is reported in place of the ground truth. The
proposed generative Archimedean copulas achieved comparable performance to
ACNet in terms of log-likelihood scores. In addition, out of the three
methods, training with maximum likelihood achieved the best results; however,
its increased computation cost, due to computing derivatives and inverses,
motivates the use of the proposed alternative losses.

Table 1: Negative log-likelihoods of learning known copulas

| Benchmark | Generative AC
---|---|---
Dataset | GT | ACNet | MLE | CvM | GAN
Clayton | -0.94 | -0.92 | -0.89 | -0.86 | -0.89
Frank | -0.90 | -0.88 | -0.89 | -0.86 | -0.89
Joe | -0.51 | -0.49 | -0.48 | -0.35 | -0.47

Table 2: Negative log-likelihoods of fitting real-world data

| Benchmark | Generative AC
---|---|---
Dataset | BF | ACNet | MLE | CvM | GAN
Boston | -0.30 | -0.27 | -0.29 | -0.30 | -0.28
INTC-MSFT | -0.19 | -0.20 | -0.16 | -0.15 | -0.17
GOOG-FB | -0.93 | -0.96 | -0.95 | -0.92 | -0.94

(a)

(b)

Figure 3: Samples from ground truth and learned copulas fitted with maximum
likelihood, goodness-of-fit and adversarial training. In (a), the copulas are
Clayton, Frank and Joe. In (b), the datasets are Boston housing, Intel-
Microsoft stocks and Google-Facebook stocks.

Samples from the learned copulas are compared to the ground truth in Figure 3.
We additionally note the differences in sampling time between our method and
the conditional sampling method used in ACNet [Ling et al., 2020]. The time to
generate 3000 samples using our method was on average $3.8\times 10^{-2}$
seconds. In comparison, the conditional sampling method via automatic
differentiation of the copula distribution followed by inversion of the
conditional distribution, takes on average $1.98\times 10^{+2}$ seconds, the
difference on the order of 3 magnitudes.

####  5.1.2 Learning Latent Distributions

The generative neural network was able to learn the latent Gamma distributions
whose Laplace transforms give the generator functions of Clayton copulas. We
show the learned latent distributions for Clayton copulas with parameters 1,
3, 5, 8 in Figure 4.

Figure 4: Gamma latent distributions of Clayton copulas with parameters 1, 3,
5, 8, learned in solid blue; ground truth in dashed black.

####  5.1.3 Learning Higher-Dimensional Copulas

While ACNet faces numerical issues for dimensions $d\geq 5$ due to repeated
automatic differentiation when computing the copula density [Ling et al.,
2020], the Laplace transform representation allows efficient computation of
higher-order derivatives without automatic differentiation.

In addition to the bivariate copulas in Section 5.1.1, we fitted Clayton,
Frank and Joe copulas for 10 and 20 dimensions. The negative log-likelihoods
are given in Table 3. When compared to the ground truth negative log-
likelihoods for 10-dimensional and 20-dimensional datasets, the learned
negative log-likelihoods were off by 2%. During our experiments, we could not
obtain a reasonably trained ACNet for high dimensions due to the computational
complexity.

Table 3: Negative log-likelihoods of learning higher-dimensional copulas | Ground Truth | Generative AC
---|---|---
Dataset | 10-dim | 20-dim | 10-dim | 20-dim
Clayton | -10.6 | -23.2 | -10.4 | -22.8
Frank | -10.4 | -23.1 | -10.4 | -23.1
Joe | -5.4 | -12.2 | -5.3 | -12.0

Moreover, while the CPU runtimes of ACNet for computing the copula density
increases exponentially with dimensions, the CPU runtimes of computing the
copula density using the Laplace transform representation increases linearly
with dimensions, as shown in Figure 5.

Figure 5: CPU runtimes for computing the likelihoods of 3000 samples from
generative Archimedean copula in solid blue; ACNet [Ling et al., 2020] in
dashed black.

###  5.2 Hierarchical Archimedean Copula

(a)

(b)

Figure 6: Samples, displayed as mirrors on the diagonal, from generative
hierarchical Archimedean copulas above in blue and from ground truth below in
black. Each plot is a bivariate margin $(U_{i},U_{j})$. In (a), a homogeneous
nested Clayton copula. In (b), a heterogeneous hierarchical Archimedean copula
with a Clayton outer generator combined with inner generators ‘12’ and ‘19’,
numbering following [Nelsen, 2010, Górecki et al., 2017]

.

We demonstrate that our model can represent more complex dependence
structures, beyond the exchangeability implied by the functional symmetry of
Archimedean copulas, and learn hierarchical Archimedean copulas.

We experiment with fitting a four-variate hierarchical Archimedean copula
$C_{\varphi_{0}}(C_{\varphi_{1}}(u_{1},u_{2}),C_{\varphi_{2}}(u_{3},u_{4}))$.
The ground truth was generated using the state-of-the-art HACopula Toolbox
[Górecki et al., 2017]. Samples from the learned copulas are compared to the
ground truth in Figure 6(b). In (a),
$C_{\varphi_{0}},C_{\varphi_{1}},C_{\varphi_{2}}$ are Clayton copulas with
parameters 1, 3, and 8. We let the outer generator be a generative Archimedean
copula. In (b), $C_{\varphi_{0}},C_{\varphi_{1}},C_{\varphi_{2}}$ are Clayton,
‘12’ and ‘19’ with parameters 0.5, 3, and 1. Since our model is compatible
with outer generators of other forms, we let the outer generator be a one-
parameter Clayton copula instead of a generative Archimedean copula.

##  6 Conclusions

We modeled Archimedean and hierarchical Archimedean copulas with generative
neural networks based on their probabilistic constructions as mixture and
nested mixture models with latent random variables. We gave efficient sampling
algorithms for sampling from the generative Archimedean and hierarchical
Archimedean copulas. We also described three methods for fitting the model to
data: maximum likelihood with the copula density, goodness-of-fit with the
empirical copula-based Cramér von-Mises statistic and adversarial training by
minimizing a divergence between true samples from data and fake samples from
the copula. Empirically, the generative Archimedean copula was able to learn
known copulas with different tail dependencies and fit real-world data. We
also showed an extension to higher-dimensional data using hierarchical
Archimedean copulas. Future work includes an end-to-end application such as
pairs trading and architecture selection for the generative neural network.

#### Acknowledgements

This work was supported in part by the Air Force Office of Scientific Research
under award number FA9550-20-1-0397.

## References

  * Bernstein [1929] Serge Bernstein.  Sur les fonctions absolument monotones.  _Acta Mathematica_ , 52:1 – 66, 1929.  10.1007/BF02592679.  URL https://doi.org/10.1007/BF02592679.
  * Cherubini et al. [2004] Umberto Cherubini, Elisa Luciano, and Walter Vecchiato.  _Copula Methods in Finance_.  The Wiley Finance Series. Wiley, 2004.  ISBN 978-0-470-86344-2.
  * Chilinski and Silva [2020] Pawel Chilinski and Ricardo Silva.  Neural likelihoods via cumulative distribution functions.  In Jonas Peters and David Sontag, editors, _Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI)_ , volume 124 of _Proceedings of Machine Learning Research_ , pages 420–429. PMLR, 03–06 Aug 2020.  URL http://proceedings.mlr.press/v124/chilinski20a.html.
  * Cramér [1928] Harald Cramér.  On the composition of elementary errors.  _Scandinavian Actuarial Journal_ , 1928(1):13–74, 1928.  10.1080/03461238.1928.10416862.  URL https://doi.org/10.1080/03461238.1928.10416862.
  * Csörgő and Teugels [1990] Sándor Csörgő and Jef L. Teugels.  Empirical laplace transform and approximation of compound distributions.  _Journal of Applied Probability_ , 27(1):88–101, 1990.  ISSN 00219002.  URL http://www.jstor.org/stable/3214597.
  * Elidan [2010] Gal Elidan.  Copula bayesian networks.  In J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, editors, _Advances in Neural Information Processing Systems_ , volume 23, pages 559–567. Curran Associates, Inc., 2010.  URL https://proceedings.neurips.cc/paper/2010/file/2a79ea27c279e471f4d180b08d62b00a-Paper.pdf.
  * Genest and Favre [2007] Christian Genest and Anne-Catherine Favre.  Everything you always wanted to know about copula modeling but were afraid to ask.  _Journal of Hydrologic Engineering_ , 12(4):347–368, 2007.  10.1061/(ASCE)1084-0699(2007)12:4(347).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%291084-0699%282007%2912%3A4%28347%29.
  * Genest et al. [2009] Christian Genest, Bruno Remillard, and David Beaudoin.  Goodness-of-fit tests for copulas: A review and a power study.  _Insurance: Mathematics and Economics_ , 44(2):199–213, 2009.  URL https://EconPapers.repec.org/RePEc:eee:insuma:v:44:y:2009:i:2:p:199-213.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 27, pages 2672–2680. Curran Associates, Inc., 2014.  URL https://proceedings.neurips.cc/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf.
  * Górecki et al. [2017] Jan Górecki, Marius Hofert, and Martin Holeňa.  On structure, family and parameter estimation of hierarchical archimedean copulas.  _Journal of Statistical Computation and Simulation_ , 87(17):3261–3324, 2017.  10.1080/00949655.2017.1365148.  URL https://doi.org/10.1080/00949655.2017.1365148.
  * Henze et al. [2012] Norbert Henze, Simos G. Meintanis, and Bruno Ebner.  Goodness-of-fit tests for the gamma distribution based on the empirical laplace transform.  _Communications in Statistics - Theory and Methods_ , 41(9):1543–1556, 2012.  10.1080/03610926.2010.542851.  URL https://doi.org/10.1080/03610926.2010.542851.
  * Hering et al. [2010] Christian Hering, Marius Hofert, Jan-Frederik Mai, and Matthias Scherer.  Constructing hierarchical archimedean copulas with lévy subordinators.  _Journal of Multivariate Analysis_ , 101(6):1428–1433, 2010.  ISSN 0047-259X.  https://doi.org/10.1016/j.jmva.2009.10.005.  URL https://www.sciencedirect.com/science/article/pii/S0047259X09001961.
  * Hernández-Lobato and Suárez [2011] José Miguel Hernández-Lobato and Alberto Suárez.  Semiparametric bivariate archimedean copulas.  _Computational Statistics & Data Analysis_, 55(6):2038–2058, 2011.
  * Hofert [2008] Marius Hofert.  Sampling archimedean copulas.  _Computational Statistics & Data Analysis_, 52(12):5163–5174, 2008.  ISSN 0167-9473.  https://doi.org/10.1016/j.csda.2008.05.019.  URL https://www.sciencedirect.com/science/article/pii/S0167947308002910.
  * Hoyos-Argüelles and Nieto-Barajas [2020] Ricardo Hoyos-Argüelles and Luis Nieto-Barajas.  A bayesian semiparametric archimedean copula.  _Journal of Statistical Planning and Inference_ , 206:298–311, 2020.  ISSN 0378-3758.  https://doi.org/10.1016/j.jspi.2019.09.015.  URL https://www.sciencedirect.com/science/article/pii/S0378375819301004.
  * Huang and Frey [2011] Jim C. Huang and Brendan J. Frey.  Cumulative distribution networks and the derivative-sum-product algorithm: Models and inference for cumulative distribution functions on graphs.  _Journal of Machine Learning Research_ , 12(10):301–348, 2011.
  * Joe [1997] Harry Joe.  _Multivariate models and dependence concepts_.  Chapman and Hall/CRC, May 1997.  ISBN 9780367803896.  10.1201/9780367803896.  URL https://www.taylorfrancis.com/books/9781466581432.
  * Joe [2014] Harry Joe.  _Dependence modeling with copulas_.  Chapman and Hall/CRC, June 2014.  ISBN 9781466583238.  10.1201/b17116.  URL https://www.taylorfrancis.com/books/9781466583238.
  * Kamthe et al. [2021] Sanket Kamthe, Samuel Assefa, and Marc Deisenroth.  Copula flows for synthetic data generation, 2021.  URL https://arxiv.org/abs/2101.00598.
  * Kimberling [1974] Clark H. Kimberling.  A probabilistic interpretation of complete monotonicity.  _aequationes mathematicae_ , 10(2):152–164, 1974\.  10.1007/BF01832852.  URL https://doi.org/10.1007/BF01832852.
  * Liew and Wu [2013] Rong Qi Liew and Yuan Wu.  Pairs trading: A copula approach.  _Journal of Derivatives & Hedge Funds_, 19(1):12–30, 2013.  10.1057/jdhf.2013.1.  URL https://doi.org/10.1057/jdhf.2013.1.
  * Ling et al. [2020] Chun Kai Ling, Fei Fang, and J. Zico Kolter.  Deep archimedean copulas.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 1535–1545. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/10eb6500bd1e4a3704818012a1593cc3-Paper.pdf.
  * Marshall and Olkin [1988] Albert W. Marshall and Ingram Olkin.  Families of multivariate distributions.  _Journal of the American Statistical Association_ , 83(403):834–841, 1988.  10.1080/01621459.1988.10478671.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1988.10478671.
  * McNeil [2008] Alexander J. McNeil.  Sampling nested archimedean copulas.  _Journal of Statistical Computation and Simulation_ , 78(6):567–581, 2008.  10.1080/00949650701255834.  URL https://doi.org/10.1080/00949650701255834.
  * McNeil and Nešlehová [2009] Alexander J. McNeil and Johanna Nešlehová.  Multivariate Archimedean copulas, d-monotone functions and $l$1-norm symmetric distributions.  _The Annals of Statistics_ , 37(5B):3059 – 3097, 2009.  10.1214/07-AOS556.  URL https://doi.org/10.1214/07-AOS556.
  * Nelsen [2010] Roger B. Nelsen.  _An introduction to Copulas_.  Springer series in statistics. Springer New York, New York, NY, 2. ed. 2006. corr. 2. pr. softcover version of original hardcover edition 2006 edition, 2010.  ISBN 9780387286785 9781441921093.  OCLC: 700190717.
  * Okhrin et al. [2013] Ostap Okhrin, Yarema Okhrin, and Wolfgang Schmid.  On the structure and estimation of hierarchical archimedean copulas.  _Journal of Econometrics_ , 173(2):189–204, 2013\.  URL https://EconPapers.repec.org/RePEc:eee:econom:v:173:y:2013:i:2:p:189-204.
  * Sato [1999] Ken-iti Sato.  _Lévy processes and infinitely divisible distributions_.  Number 68 in Cambridge studies in advanced mathematics. Cambridge University Press, 1999.  ISBN 9780521553025.
  * Savu and Trede [2010] Cornelia Savu and Mark Trede.  Hierarchies of archimedean copulas.  _Quantitative Finance_ , 10(3):295–304, 2010\.  10.1080/14697680902821733.  URL https://doi.org/10.1080/14697680902821733.
  * Sklar [1959] Abe Sklar.  Fonctions de répartition à n dimensions et leurs marges.  8:229–231, 1959.
  * Stander et al. [2013] Yolanda Stander, Daniel Marais, and Ilse Botha.  Trading strategies with copulas.  _Journal of Economic and Financial Sciences_ , 6(1):83–107, 2013.  10.10520/EJC135921.
  * Tagasovska et al. [2019] Natasa Tagasovska, Damien Ackerer, and Thibault Vatter.  Copulas as high-dimensional generative models: Vine copula autoencoders.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/15e122e839dfdaa7ce969536f94aecf6-Paper.pdf.
  * Widder [1941] David Vernon Widder.  _Laplace Transform (PMS-6)_.  Princeton University Press, 1941.  doi:10.1515/9781400876457.  URL https://doi.org/10.1515/9781400876457.
  * Wiese et al. [2019] Magnus Wiese, Robert Knobloch, and Ralf Korn.  Copula & marginal flows: Disentangling the marginal from its joint, 2019\.  URL https://arxiv.org/abs/1907.03361.
  * Wilson and Ghahramani [2010] Andrew G Wilson and Zoubin Ghahramani.  Copula processes.  In J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, editors, _Advances in Neural Information Processing Systems_ , volume 23, pages 2460–2468. Curran Associates, Inc., 2010.  URL https://proceedings.neurips.cc/paper/2010/file/fc8001f834f6a5f0561080d134d53d29-Paper.pdf.
  * Xu and Darve [2020] Kailai Xu and Eric Darve.  Calibrating multivariate Lévy processes with neural networks.  In Jianfeng Lu and Rachel Ward, editors, _Proceedings of The First Mathematical and Scientific Machine Learning Conference_ , volume 107 of _Proceedings of Machine Learning Research_ , pages 207–220, Princeton University, Princeton, NJ, USA, 20–24 Jul 2020. PMLR.  URL http://proceedings.mlr.press/v107/xu20a.html.

##  7 Supplementary Material

###  7.1 Probabilistic Construction of Archimedean and Hierarchical
Archimedean Copulas

Copulas can be derived from cumulative distribution functions (CDFs) via
Sklar’s theorem, i.e. specify a joint CDF $F$, compute univariate CDFs
$F_{1},\cdots,F_{d}$ from the joint CDF, then obtain the copula as
$C(\mathbf{u})=F(F_{1}^{-1}(u_{1}),\cdots,F_{d}^{-1}(u_{d})),\mathbf{u}\in[0,1]^{d}$.
Sklar’s theorem also applies to survival functions, i.e. for joint survival
function
$\bar{F}(\mathbf{x})=P(X_{1}>x_{1},\cdots,X_{d}>x_{d}),\mathbf{x}\in\mathbb{R}^{d}$,
with univariate survival functions $\bar{F}_{1},\cdots,\bar{F}_{d}$ where
$\bar{F}_{j}=P(X_{j}>x_{j})$, the copula which couples $\bar{F}$ to
$\bar{F}_{1},\cdots,\bar{F}_{d}$ is called the _survival copula_ and is given
as the copula $C$ for which
$\bar{F}(x)=C(\bar{F}_{1}(x_{1}),\cdots,\bar{F}_{d}(x_{d}))$.

####  7.1.1 Archimedean Copulas

We restate the probabilistic construction found in [Joe, 2014] Chapter 3.2,
following [Marshall and Olkin, 1988]:

Let $G_{1},\cdots,G_{d}$ be univariate CDFs. Let $Q\sim F_{Q}$ be a positive
random variable with Laplace transform $\varphi_{Q}$, let $X_{1},\cdots,X_{d}$
be dependent random variables that are conditionally independent given $Q=q$
such that $[X_{j}|Q=q]\sim G_{j}^{q},\;q>0$.

The joint CDF is:

$F(x_{1},\cdots,x_{d})=\int_{0}^{\infty}G_{1}^{q}(x_{1})\cdots
G_{d}^{q}(x_{d})dF_{Q}(q)=\varphi_{Q}(-\log G_{1}(x_{1})-\cdots-\log
G_{d}(x_{d})),$ (15)

with univariate CDFs obtained from the joint CDF as:

$F_{j}(x_{j})=\int_{0}^{\infty}G_{j}^{q}(x_{j})dF_{Q}(q)=\varphi_{Q}(-\log
G_{j}(x_{j})),\;j\in\\{1,\cdots,d\\},$ (16)

and inverse:

$F_{j}^{-1}(u_{j})=G_{j}^{-1}(\exp\\{-\varphi_{Q}^{-1}(u_{j})\\}),\;u_{j}\in(0,1),\;j\in\\{1,\cdots,d\\},$
(17)

such that the copula via Sklar’s theorem is:

$C(\mathbf{u})=F(F_{1}^{-1}(u_{1}),\cdots,F_{d}^{-1}(u_{d}))=\varphi_{Q}(\varphi_{Q}^{-1}(u_{1})+\cdots+\varphi_{Q}^{-1}(u_{d})),\;\mathbf{u}\in[0,1]^{d}.$
(18)

The multivariate extension of bivariate Archimedean copulas was introduced in
[Kimberling, 1974] with the condition that the above expression is a valid
copula for any $d$ whenever $\varphi$, known as the _generator_ of the
Archimedean copula, is _completely monotone_ , i.e. the Laplace transform of a
positive random variable [Bernstein, 1929, Widder, 1941]. The mixture
representation with Laplace transform generators and an efficient algorithm
for sampling from the mixture representation was subsequently given in
[Marshall and Olkin, 1988].

We restate the sampling algorithm found in [McNeil, 2008], following [Marshall
and Olkin, 1988]:

Consider $(U_{1},\cdots,U_{d})=(\varphi(E_{1}/M),\cdots,\varphi(E_{d}/M))$,
where $(E_{1},\cdots,E_{d})\sim i.i.d.\;\text{Exp}(1)$ are independent and
identically distributed unit exponentials and $M\sim F_{M}$ is a positive
random variable with Laplace transform $\varphi_{M}$.

$\begin{split}P(U_{1}\leq u_{1},\cdots,U_{d}\leq
u_{d})&=\int_{0}^{\infty}P(U_{1}\leq u_{1},\cdots,U_{d}\leq
u_{d}|M=s)dF_{M}(s)\\\
&=\int_{0}^{\infty}e^{-s(\varphi^{-1}(u_{1})+\cdots+\varphi^{-1}(u_{d})}dF_{M}(s)\\\
&=\varphi_{M}(\varphi_{M}^{-1}(u_{1})+\cdots+\varphi_{M}^{-1}(u_{d})).\end{split}$
(19)

Thus an algorithm for sampling $\mathbf{U}\sim C$ is to sample $M\sim F_{M}$
with Laplace transform $\varphi_{M}$, sample $(E_{1},\cdots,E_{d})\sim
i.i.d.\;\text{Exp}(1)$, then compute
$(U_{1},\cdots,U_{d})=(\varphi_{M}(E_{1}/M),\cdots,\varphi_{M}(E_{d}/M))$.

####  7.1.2 Hierarchical Archimedean Copulas

A simple nested mixture representation involving Laplace transform generators
was introduced in [Joe, 1997]. Conditions for the nested copula to be a valid
copula, called _sufficient nesting conditions_ , was derived based on the
composition of an _outer generator_ $\varphi_{0}$ and an _inner generator_
$\varphi_{1}$ to get a completely monotone Laplace transform _nested
generator_ $e^{-\nu_{0}\varphi_{0}^{-1}\circ\varphi_{1}}$, where $\nu_{0}$ is
a positive random variable with Laplace transform $\varphi_{0}$, such that
$\varphi_{0},\varphi_{1}$ and $(\varphi_{0}^{-1}\circ\varphi_{1})^{\prime}$
are completely monotone.

We illustrate with a simple three-dimensional example found in [McNeil, 2008],
following [Joe, 1997]. Consider the hierarchical Archimedean copula:

$C(u_{1},u_{2},u_{3})=\varphi_{0}(\varphi_{0}^{-1}(u_{1})+\varphi_{0}^{-1}\circ\varphi_{1}(\varphi_{1}^{-1}(u_{2})+\varphi_{1}^{-1}(u_{3}))),$
(20)

where $\varphi_{0},\varphi_{1}$ are Laplace transform generators of
Archimedean copulas. We would like to express the above as a mixture of
conditionally independent CDFs. Let $G_{0}$ be a distribution with Laplace
transform $\varphi_{0}$:

$\displaystyle C(u_{1},u_{2},u_{3})$
$\displaystyle=\varphi_{0}(\varphi_{0}^{-1}(u_{1})+\varphi_{0}^{-1}\circ\varphi_{1}(\varphi_{1}^{-1}(u_{2})+\varphi_{1}^{-1}(u_{3})))$
(21)
$\displaystyle=\int_{0}^{\infty}e^{-\nu_{0}\varphi_{0}^{-1}(u_{1})}e^{-\nu_{0}\varphi_{0}^{-1}\circ\varphi_{1}(\varphi_{1}^{-1}(u_{2})+\varphi_{1}^{-1}(u_{3}))}dG_{0}(\nu_{0})$
$\displaystyle=\int_{0}^{\infty}F^{\nu_{0}}_{0}(u_{1})C_{01}(F_{0}^{\nu_{0}}(u_{2}),F_{0}^{\nu_{0}}(u_{3});\nu_{0})dG_{0}(\nu_{0}),$

where $F_{0}(\cdot):=e^{-\varphi_{0}^{-1}(\cdot)}$ and $F^{\nu_{0}}_{0}$ is a
valid CDF for any $\nu_{0}>0$. In addition, $C_{01}(\cdot;\nu)$ is an
Archimedean copula with Laplace transform generator
$\varphi_{01}(\cdot;\nu_{0})=e^{-{\nu_{0}}\varphi_{0}^{-1}\circ\varphi_{1}(\cdot)}$
and generator inverse
$\varphi_{01}^{-1}(\cdot;\nu_{0})=\varphi_{1}^{-1}\circ\varphi_{0}(-\log(\cdot)/\nu_{0})$,
such that $C_{01}(\cdot;\nu_{0})$ taking marginals $F_{0}^{\nu_{0}}(u_{2})$
and $F_{0}^{\nu_{0}}(u_{3})$ as inputs gives:

$\displaystyle C_{01}(F_{0}^{\nu_{0}}(u_{2}),F_{0}^{\nu_{0}}(u_{3});\nu_{0})$
$\displaystyle=\varphi_{01}(\varphi_{01}^{-1}(F_{0}^{\nu_{0}}(u_{2});\nu_{0})+\varphi_{01}^{-1}(F_{0}^{\nu_{0}}(u_{3});\nu_{0});\nu_{0})$
(22)
$\displaystyle=e^{-{\nu_{0}}\varphi_{0}^{-1}\circ\varphi_{1}(\varphi_{1}^{-1}\circ\varphi_{0}(-\log(e^{-{\nu_{0}}\varphi_{0}^{-1}(u_{2})})/{\nu_{0}})+\varphi_{1}^{-1}\circ\varphi_{0}(-\log(e^{-{\nu_{0}}\varphi_{0}^{-1}(u_{3})})/{\nu_{0}}))}$
$\displaystyle=e^{-{\nu_{0}}\varphi_{0}^{-1}\circ\varphi_{1}(\varphi_{1}^{-1}(u_{2})+\varphi_{1}^{-1}(u_{3}))}.$

The completely monotone property of the Laplace transform generator
$\varphi_{01}(\cdot;\nu)=e^{-\nu\varphi_{0}^{-1}\circ\varphi_{1}(\cdot)}$ then
implies $(\varphi_{0}^{-1}\circ\varphi_{1})^{\prime}$ is completely monotone.
In addition, letting $G_{01}(\cdot;\nu_{0})$ be a distribution with Laplace
transform $\varphi_{01}(\cdot;\nu_{0})$, we express the hierarchical
Archimedean copula as a nested mixture of conditionally independent CDFs:

$\displaystyle C(u_{1},u_{2},u_{3})$
$\displaystyle=\int_{0}^{\infty}F^{\nu_{0}}_{0}(u_{1})C_{01}(F_{0}^{\nu_{0}}(u_{2}),F_{0}^{\nu_{0}}(u_{3});\nu_{0})dG_{0}(\nu_{0})$
(23)
$\displaystyle=\int_{0}^{\infty}F^{\nu_{0}}_{0}(u_{1})\int_{0}^{\infty}e^{-{\nu_{01}}\varphi_{01}^{-1}(F_{0}^{\nu_{0}}(u_{2}))}e^{-{\nu_{01}}\varphi_{01}^{-1}(F_{0}^{\nu_{0}}(u_{3}))}dG_{01}({\nu_{01}};\nu_{0})dG_{0}(\nu_{0})$
$\displaystyle=\int_{0}^{\infty}F^{\nu_{0}}_{0}(u_{1})\int_{0}^{\infty}e^{-{\nu_{01}}\varphi_{1}^{-1}(u_{2})}e^{-{\nu_{01}}\varphi_{1}^{-1}(u_{3})}dG_{01}({\nu_{01}};\nu_{0})dG_{0}(\nu_{0})$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}F^{\nu_{0}}_{0}(u_{1})F^{\nu_{01}}_{1}(u_{2})F^{\nu_{01}}_{1}(u_{3})dG_{01}({\nu_{01}};\nu_{0})dG_{0}(\nu_{0}),$

where $F_{1}(\cdot):=e^{-\varphi_{1}^{-1}(\cdot)}$ and $F^{\nu_{01}}_{1}$ is a
valid CDF for any $\nu_{01}>0$.

This construction and condition were restated for nesting to arbitrary depth
in [McNeil, 2008].

Based on the mixture representation, McNeil [2008] also provided algorithms
for sampling nested Clayton and nested Gumbel copulas. It was also showed that
Clayton and Gumbel copulas are unfortunately not compatible for nesting. The
challenge was to find combinations of known distributions with
$\varphi_{0},\varphi_{1}$ and $e^{-\nu_{0}\varphi_{0}^{-1}\circ\varphi_{1}}$
as their Laplace transforms. Sampling using McNeil [2008]’s algorithm for
nested Ali-Mikhail-Haq, nested Frank, nested Joe, more parametric families and
numerical inversion of Laplace transform was by [Hofert, 2008]. It was
subsequently recognized in [Hering et al., 2010] that the sufficient nesting
condition for $(\varphi_{0}^{-1}\circ\varphi_{1})^{\prime}$ to be completely
monotone can be satisfied by letting $\varphi_{1}=\varphi_{0}\circ\psi_{1}$,
where $\psi_{1}$ is the Laplace exponent, with completely monotone derivative,
in the Laplace transform of Lévy subordinators.

We restate the probabilistic construction with Lévy subordinators from [Hering
et al., 2010]:

For each ‘time’ $t\geq 0$, the Laplace transform of a Lévy subordinator
$\Lambda_{t}$, i.e. a non-decreasing Lévy processes such as the compound
Poisson process, is given as:

$\mathbf{E}[e^{-x\Lambda_{t}}]=e^{-t\psi(x)},$ (24)

where $\psi(x)$ is the Laplace exponent.

Consider
$(\frac{E_{1,1}}{\Lambda_{1}(M)},\cdots,\frac{E_{1,d_{1}}}{\Lambda_{1}(M)},\cdots,\cdots,\cdots,\frac{E_{J,1}}{\Lambda_{J}(M)},\cdots,\frac{E_{J,d_{J}}}{\Lambda_{J}(M)})$,
where $E_{j,i}\sim i.i.d.\;\text{Exp}(1)$, $\Lambda_{j}$ are Lévy
subordinators with Laplace exponents $\psi_{j}$, and $\Lambda_{j}$ are
evaluated at a common ‘time’ $t=M$, where M is a positive random variable with
Laplace transform $\varphi_{0}$. The hierarchical Archimedean copula is then
constructed using the survival analog of Sklar’s theorem.

The joint survival function is:

$\begin{split}P(\frac{E_{j,i}}{\Lambda_{j}(M)}>x_{j,i}\text{\;for
all}j,i)&=\mathbf{E}[\prod_{j=1}^{J}\prod_{i=1}^{d_{J}}e^{-\Lambda_{j}(M)x_{j,i}}]\\\
&=\mathbf{E}[\prod_{j=1}^{J}e^{-\Lambda_{j}(M)\sum_{i=1}^{d_{J}}x_{j,i}}]\\\
&=\mathbf{E}[\prod_{j=1}^{J}e^{-M\psi_{j}(\sum_{i=1}^{d_{j}}x_{j,i})}]\\\
&=\mathbf{E}[e^{-M\sum_{j=1}^{J}\psi_{j}(\sum_{i=1}^{d_{j}}x_{j,i})}]\\\
&=\varphi_{0}(\sum_{j=1}^{J}\varphi_{0}^{-1}\circ(\varphi_{0}\circ\psi_{j})(\sum_{i=1}^{d_{j}}x_{j,i})),\end{split}$
(25)

and each component $\frac{E_{j,i}}{\Lambda_{j}(M)}$ has survival function:

$P(\frac{E_{j,i}}{\Lambda_{j}(M)}>x)=\mathbf{E}[e^{-x\Lambda_{j}(M)}]=\mathbf{E}[e^{-M\psi_{j}(x)}]=(\varphi_{0}\circ\psi_{j})(x).$
(26)

Using the survival analog of Sklar’s theorem, given the above univariate
survival functions, the hierarchical Archimedean copula $C$ with outer
generator $\varphi_{0}$ and inner generators
$\varphi_{j}=\varphi_{0}\circ\psi_{j}$, we recover the above joint survival
function.

###  7.2 Experiment Details

Following the experiment setup in [Ling et al., 2020], the commonly-used
copulas were the Clayton, Frank and Joe copulas, each governed by a single
parameter and chosen to be 5, 15, and 3 respectively. Each dataset had 2000
train and 1000 test points. The real-world data were the Boston housing,
Intel-Microsoft (INTC-MSFT) stocks and Google-Facebook (GOOG-FB) stocks. Each
dataset was divided into train and test points in a 3:1 ratio, then rank-
normalized to get approximately uniform margins.

Similar to the experiment parameters in [Ling et al., 2020], the tolerance for
Newton’s root-finding method was 1e-10. The generative neural network was a
multilayer perceptron of comparable size, 2 hidden layers, each of width 10.
We used $\text{U}(0,1)$ as the input source of randomness, default weight
initialization, LeakyReLU intermediate activations and $\exp(.)$ output
activation. For training with maximum likelihood, we used the same
optimization parameters: stochastic gradient descent (SGD) with learning rate
1e-5 and momentum 0.9 on sum of log-likelihoods. For training with goodness-
of-fit, we used SGD with learning rate 1e-3 and momentum 0.9. For adversarial
training, we used Adam with learning rate 1e-4, momentum 0.9 and betas (0.5,
0.999). The discriminative neural network had a single hidden layer of width
20, default weight initialization, LeakyReLU intermediate activations and
$\text{sigmoid}(.)$ output activation. All training methods used the same
batch size of 200 and converged within 10k epochs. We reported the results at
10k epoch. Experiments were conducted using PyTorch, on a 2.7 GHz Intel Core
i7 with 16 GB of RAM.

To reduce computation complexity during training, we used a smaller number
$L=100$ samples from the generative neural network to approximate the Laplace
transforms. To increase inference accuracy for evaluation, we used a larger
number $L=1000$ samples from the generative neural network to approximate the
Laplace transforms.

####  7.2.1 Enforcing Structural Properties

Compared to vanilla GAN [Goodfellow et al., 2014], our generating network must
satisfy an Archimedean copula. We show this via the training progression for
learning a Clayton copula in Figure 7.

Figure 7: Training progression at epochs 0, 500, 5000 and 10000, for learning
a Clayton copula. Samples from our copula, shown on top, must satisfy an
Archimedean copula, while that from a vanilla GAN, shown below, may not.


* Massoni et al. [2002] J Massoni, R Saurel, B Nkonga, and R Abgrall.  Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and eulerian methods for interface problems between compressible fluids with heat transfer.  _International journal of heat and mass transfer_ , 45(6):1287–1307, 2002.  ISSN 0017-9310.
  * Michael and Nikiforakis [2016] L. Michael and N. Nikiforakis.  A hybrid formulation for the numerical simulation of condensed phase explosives.  _Journal of Computational Physics_ , 316:193–217, 2016\.  ISSN 0021-9991.
  * Michael and Nikiforakis [2018] L. Michael and N. Nikiforakis.  A multi-physics methodology for the simulation of reactive flow and elastoplastic structural response.  _Journal of Computational Physics_ , 367(C):1–27, 2018.  ISSN 0021-9991.
  * Michael et al. [2019] Louisa Michael, Stephen T. Millmore, and Nikolaos Nikiforakis.  A multi-physics methodology for four-states of matter.  2019\.
  * Miller and Colella [2002] G.H. Miller and P. Colella.  A conservative three-dimensional Eulerian method for coupled solid–fluid shock capturing.  _Journal of Computational Physics_ , 183(1):26 – 82, 2002.  ISSN 0021-9991.  doi: https://doi.org/10.1006/jcph.2002.7158.  URL http://www.sciencedirect.com/science/article/pii/S0021999102971585.
  * Millett [2015] Jeremy C. F Millett.  Modifications of the response of materials to shock loading by age hardening.  _Metallurgical and Materials Transactions A_ , 46(10):4506–4517, 2015.  ISSN 1073-5623.
  * Moler and Loan [2003] Cleve Moler and Charles Van Loan.  Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later.  _SIAM review_ , 45(1):3–49, 2003.  ISSN 0036-1445.
  * Murakami [2012] S. Murakami.  _Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture_.  Springer, frist edition. edition, 2012.  ISBN 978-94-007-2665-9.
  * Ndanou et al. [2015] S. Ndanou, N. Favrie, and S. Gavrilyuk.  Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation.  _Journal of Computational Physics_ , 295:523–555, 2015\.  ISSN 0021-9991.
  * Ottosen [2005] Niels Saabye Ottosen.  _The mechanics of constitutive modeling / Niels Saabye Ottosen, Matti Ristinmaa._ Elsevier, Amsterdam ; London, 1st ed. edition, 2005.  ISBN 008044606X.
  * Petitpas et al. [2009] Fabien Petitpas, Jacques Massoni, Richard Saurel, Emmanuel Lapebie, and Laurent Munier.  Diffuse interface model for high speed cavitating underwater systems.  _International Journal of Multiphase Flow_ , 35(8):747–759, 2009.  ISSN 0301-9322.
  * Pirondi and Bonora [2003] A Pirondi and N Bonora.  Modeling ductile damage under fully reversed cycling.  _Computational Materials Science_ , 26:129–141, 2003.  ISSN 0927-0256.
  * Poirier [2000] J.-P. Poirier.  _Introduction to the physics of the Earth’s interior_.  Cambridge University Press, second edition edition, 2000.
  * Rezaee-Hajidehi et al. [2020] M Rezaee-Hajidehi, K Tuma, and S Stupkiewicz.  A note on pade approximants of tensor logarithm with application to hencky-type hyperelasticity.  _Computational mechanics_ , 2020.  ISSN 0178-7675.
  * Sambasivan and Udaykumar [2009] Shiv Kumar Sambasivan and H. S Udaykumar.  Ghost fluid method for strong shock interactions part 1: Fluid-fluid interfaces.  _AIAA Journal_ , 47(12):2907–2922, 2009.  ISSN 0001-1452.
  * Saurel and Abgrall [1999] R. Saurel and R. Abgrall.  A simple method for compressible multifluid flow.  _SIAM Journal on Scientific Computing_ , 21(3), 1999.  ISSN 1064-8275.
  * Schoch et al. [2013] Stefan Schoch, Kevin Nordin-Bates, and Nikolaos Nikiforakis.  An Eulerian algorithm for coupled simulations of elastoplastic-solids and condensed-phase explosives.  _Journal of Computational Physics_ , 252:163–194, 2013\.  ISSN 0021-9991.
  * Tavelli et al. [2020] Maurizio Tavelli, Simone Chiocchetti, Evgeniy Romenski, Alice-Agnes Gabriel, and Michael Dumbser.  Space-time adaptive ader discontinuous galerkin schemes for nonlinear hyperelasticity with material failure.  _Journal of computational physics_ , 422:109758, 2020.  ISSN 0021-9991.
  * Thomas et al. [2016] S. Thomas, L. Veeser, W. Turley, and R. Hixson.  Comparisons of CTH simulations with measured wave profiles for simple flyer plate experiments.  _Journal of Dynamic Behavior of Materials_ , 2(3):365–371, 2016.  ISSN 2199-7446.
  * Vitali et al. [2012] E. Vitali, I. N. Lomov, T. H. Antoun, and D. H. Fujino.  An extended Eulerian method for contacts in Godunov formulations.  _International Journal for Numerical Methods in Engineering_ , 92:1139–1156, 2012.
  * Vitali and Benson [2012] Efrem Vitali and David Benson.  Modeling localized failure with arbitrary Lagrangian Eulerian methods.  _Computational Mechanics_ , 49(2):197–212, 2012\.  ISSN 0178-7675.
  * Vogler et al. [2003] Tracy J Vogler, Tom F Thornhill, William D Reinhart, Lalit C Chhabildas, Dennis E Grady, Leonard T Wilson, Omar A Hurricane, and Anne Sunwoo.  Fragmentation of materials in expanding tube experiments.  _International Journal of Impact Engineering_ , 29(1):735–746, 2003.  ISSN 0734-743X.
  * Wallis et al. [2021] Tim Wallis, Philip T. Barton, and Nikolaos Nikiforakis.  A diffuse interface model of reactive-fluids and solid-dynamics.  _Computers & Structures_, 254:106578, 2021.  ISSN 0045-7949.  doi: https://doi.org/10.1016/j.compstruc.2021.106578.  URL https://www.sciencedirect.com/science/article/pii/S0045794921001000.
  * Wang et al. [2006] C. W. Wang, T. G. Liu, and B. C. Khoo.  A real ghost fluid method for the simulation of multimedium compressible flow.  _SIAM Journal on Scientific Computing_ , 28(1):278–302, 2006.  ISSN 1064-8275.
  * Xiao et al. [2005] F. Xiao, Y. Honma, and T. Kono.  A simple algebraic interface capturing scheme using hyperbolic tangent function.  _International Journal for Numerical Methods in Fluids_ , 48(9):1023–1040, 7 2005.  ISSN 1097-0363.  doi: 10.1002/fld.975.  URL https://doi.org/10.1002/fld.975.
  * Youngs [1984] David Youngs.  An interface tracking method for a 3d Eulerian hydrodynamics code.  Technical report, AWRE, 1984.
  * Zhang et al. [2019] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale.  AMReX: a framework for block-structured adaptive mesh refinement.  _Journal of Open Source Software_ , 4(37):1370, May 2019.  doi: 10.21105/joss.01370.  URL https://doi.org/10.21105/joss.01370.