# Deep Extreme Value Copulas for Estimation and Sampling

Ali Hasan1  Equal contribution. Khalil Elkhalil,2 11footnotemark: 1 João M.
Pereira,3 Sina Farsiu,12 Jose H. Blanchet,4 Vahid Tarokh2
1Department of Biomedical Engineering, Duke University
2Department of Electrical and Computer Engineering, Duke University
3Department of Electrical and Computer Engineering, University of Texas at
Austin
4Department of Management Science and Engineering, Stanford University

###### Abstract

We propose a new method for modeling the distribution function of high
dimensional extreme value distributions. The Pickands dependence function
models the relationship between the covariates in the tails, and we learn this
function using a neural network that is designed to satisfy its required
properties. Moreover, we present new methods for recovering the spectral
representation of extreme distributions and propose a generative model for
sampling from extreme copulas. Numerical examples are provided demonstrating
the efficacy and promise of our proposed methods.

##  1 Introduction

Modeling the occurrence of extreme events is an important task in many
disciplines, such as medicine, environmental science, engineering, and
finance. For example, understanding the probability of a patient having an
adverse reaction to medication or the distribution of economic shocks is
critical to mitigating the associated effects of these events. However, these
events are rare in occurrence and often difficult to characterize with
traditional statistical tools. This has been the primary focus of extreme
value theory (EVT), which describes how to extrapolate the occurrence of rare
events outside the range of available data [1]. In the one-dimensional case,
EVT provides remarkably simple models for the distribution of the maximum of
an infinite number of independent and identically distributed (i.i.d) random
variables. The latter result is due to the fundamental result by Fisher-
Tippet-Gnedenko theorem [2].

The problem in high dimensions becomes more difficult due to the challenge of
modeling complex interactions between different extreme variables. These
interactions are often captured by extreme value copulas [3]. Unfortunately,
these are hard to identify and learn due to the lack of data of simultaneous
occurrences of rare events. Moreover, unlike the one-dimensional case,
modeling multivariate extreme value (MEV) distributions does not easily permit
analytical forms of the underlying density, leading to difficulties in
performing inference tasks using conventional methods such as maximum
likelihood estimation (MLE). Due to the pioneering work of Pickands [4], MEV
distributions can be entirely determined by tail dependence functions [3].
Therefore, learning MEV distributions amounts to learning and estimating the
underlying tail dependence function.

Motivated by these challenges, we present an algorithm for learning
multivariate extreme value distributions. More precisely, we propose to model
the Pickands dependence function, which is the restriction of the tail
dependence function on the unit simplex, using flexible parametric models. By
exploiting that Pickands dependence functions turn out to be convex [2], we
consider the use of input convex neural networks (ICNNs), which have the
capabilities to approximate any convex function up to desirable accuracy [5,
6]. This allows us to enforce the convexity property in Pickands dependence
functions while still permitting a flexible model. For training, we exploit a
property discussed in [4], stating that a suitable statistic of the covariates
can be transformed to be distributed as an exponential random variable with
the rate parameter given by the Pickands dependence function. This allows
efficient training of the ICNN by optimizing the likelihood function of the
exponential distribution using an ICNN to parameterize the rate. We provide
numerical results to demonstrate the potential of the proposed Pickands based
ICNN for estimating survival probabilities in the context of both synthetic
and real data as compared to state of the art methods.

On another front, we provide a method for recovering the spectral
representation [7, 8] of MEV distributions which is a key ingredient to create
a generative model that makes it possible to directly sample from the
underlying MEV distribution. This allows us to use existing samplers such as
[9, 10] to sample from the MEV distribution. Furthermore, recovering the
spectral representation from the Pickands dependence function also gives
additional insight for possible clustering of the extremes, as described in
[8].

For completeness, we provide the common literature dealing with the described
challenges in MEV distributions.

Related Work. A number of techniques have been developed to estimate extreme
value copulas from data (see [3] for a self-contained review). Arguably, the
most relevant to the present work is that by Pickands [4] where the notion of
the multivariate extreme value distribution and a non parametric estimator
were first proposed. By appropriately transforming and scaling the marginals,
we can simplify the problem to only considering exponential random variables
with parameters given by the Pickands function (please refer to Section 3 for
more details on such transformations). This allows to use estimators such as
MLE on the transformed variables to learn the Pickands dependence function.

Many follow-up works describe various modifications of the method proposed by
[4]. Additional approaches such as the work in [11, 12] describe similar, but
alternative takes on estimating the dependence function. These methods,
however, do not guarantee the estimate to be a convex function, a necessary
condition of the Pickands dependence function. In [13], the authors consider a
projection of a nonparametric estimator to a convex function represented as a
Bernstein polynomial. However, the number of parameters required significantly
increases with both the amount of data and the dimensionality, making it
difficult for higher dimensional problems or problems with many data points.
Our proposed method can be used with any of the previously mentioned
nonparameteric estimators while also guaranteeing convexity. In practice, we
observed that the estimator in [4] is the most effective. For additional
details, please refer to the review on extreme value copulas in [14].

Main contributions. First, we describe an approach to parameterize Pickands
dependence functions that scales effectively to high dimensions. Second, we
provide numerical results to demonstrate the effectiveness of the proposed
approach for both synthetic and real data. Finally, we propose an algorithm
that permits to recover the spectral representation of Pickands dependence
functions which in turns allows to directly sample from the MEV distribution
without resorting to the underlying density.

Outline. In Section 2, we provide the necessary background on extreme value
theory and its extension to the high dimensional setting. In Section 3, we
describe how to model and train Pickands dependence functions with ICNNs.
Then, in Section 4, we provide details on how to recover the underlying
spectral representations of MEV distributions and propose an algorithm for
sampling. We provide numerical results in Section 5. Finally, we provide some
concluding remarks in Section 6.

##  2 Extreme value copulas and Pickands dependence function

###  2.1 Extreme value theory

Extreme value theory (EVT) is a fundamental tool that characterizes the
behavior of maxima of $n$ independent and identically distributed (i.i.d)
random variables $X_{1},\cdots,X_{n}$ with continuous distribution function
$F$. More precisely, let $M_{n}=\max_{1\leq i\leq n}X_{i}$, we are interested
in finding sequences of real numbers $a_{n}>0$ and $b_{n}$ such that the limit
$\mathbb{P}\left[\left(M_{n}-b_{n}\right)/a_{n}\leq x\right]\to H(x)$ as
$n\to\infty$ is non-degenerate. We then say that $F$ is in the maximum domain
of attraction of $H$ or equivalently $F\in\text{MDA}(H)$. This limit is fully
identified by the generalized extreme value (GEV) distribution given by:

$\begin{split}H_{\xi}(x)=\left\\{\begin{matrix}\exp\left(-(1+\xi
x)^{-1/\xi}\right),&\text{if }\xi\neq 0\\\ \exp{(-e^{-x})},&\text{if
}\xi=0\end{matrix}\right.\end{split}$ (1)

where $1+\xi x>0$ and $\xi$ is the shape parameter indicating the thickness of
the tail. The following theorem [2] due to Fisher, Tippet and Gnedenko states
the fundamental result of EVT.

######  Theorem 1.

If $F\in\text{MDA}(H)$, and if the limit $H$ exists then it belongs to the
class of GEV distributions, i.e. $H=H_{\xi}$ for some real number $\xi$.

###  2.2 Extreme value copulas

After characterizing the limit behavior of maxima of i.i.d random variables,
the natural next step is to characterize the joint behavior of component-wise
maxima of i.i.d random vectors. In more concrete terms, let
$\left(X_{i}^{(1)},\cdots,X_{i}^{(d)}\right)\in\mathbb{R}^{d}$ for
$i\in\\{1,\cdots,n\\}$ be a sample of i.i.d random vectors with common
continuous probability distribution $F$, marginals $F_{1},\cdots,F_{d}$ and
copula $C_{F}$. Recall that $C_{F}:[0,1]^{d}\to[0,1]$ satisfies:

$\displaystyle C_{F}(u_{1},\cdots,u_{d})$
$\displaystyle=\mathbb{P}\left[F_{1}(X^{(1)})\leq
u_{1},\cdots,F_{d}(X^{(d)})\leq u_{d}\right].$

Let the vector of component-wise maxima given by:

$\displaystyle M_{n}=\left(M_{n}^{(1)},\cdots,M_{n}^{(d)}\right),$

where $M_{n}^{(k)}=\max_{1\leq i\leq n}X_{i}^{(k)}$, $k\in\\{1,\cdots,d\\}$.
Similarly, each component-wise maxima $M_{n}^{(k)}$ is normalized with
sequences of real numbers $a_{n}^{(k)}>0$ and $b_{n}^{(k)}$ such that the
corresponding limiting marginal is non-degenerate. Let $C_{n}$ be the copula
of $\bar{M}_{n}$ given by:

$\displaystyle\bar{M}_{n}=\left(\frac{M_{n}^{(1)}-b_{n}^{(1)}}{a_{n}^{(1)}},\cdots,\frac{M_{n}^{(d)}-b_{n}^{(d)}}{a_{n}^{(d)}}\right).$

Then the following relation holds.

$\displaystyle
C_{n}(u_{1},\cdots,u_{d})=C_{F}(u_{1}^{1/n},\cdots,u_{d}^{1/n})^{n},$

for all $(u_{1},\cdots,u_{d})\in[0,1]^{d}$. Similar to the one-dimensional
case, we are interested in finding the limiting copula $C$ of $C_{n}$ as
$n\to\infty$. The limiting copula is called an extreme-value copula and we
similarly say that $C_{F}\in\text{MDA}(C)$. Before stating the main theorem
that characterizes extreme value copulas, it is necessary to provide some
technical definitions related to tail dependence functions.

######  Definition 1 (Tail dependence function).

A function $\ell:[0,\infty)^{d}\to[0,\infty)$ is a tail dependence function if
for all $\left(x_{1},\cdots,x_{d}\right)\in[0,\infty)^{d}$, the following
conditions are satisfied:

  * •

(i) $\ell$ is convex and homogeneous of order $1$, i.e.
$\ell(cx_{1},\cdots,cx_{d})=c.\ell(x_{1},\cdots,x_{d})$, for all $c>0$.

  * •

(ii) $\max_{1\leq k\leq
d}x_{k}\leq\ell(x_{1},\cdots,x_{d})\leq\sum_{k=1}^{d}x_{k}$.

With the above definition of tail dependence functions at hands, we are in a
position to state the following theorem [3] which provides the mathematical
foundation of extreme value copulas.

######  Theorem 2.

If $C$ is a $d-$variate extreme value copula then there exists a tail
dependence function $\ell:[0,\infty)^{d}\to[0,\infty)$ such that:

$\vspace{-.2cm}C\left(u_{1},\cdots,u_{d}\right)=e^{-\ell\left(-\log
u_{1},\cdots,-\log u_{d}\right)},$ (2)

where $\left(u_{1},\cdots,u_{d}\right)\in(0,1]^{d}$. Using the homogeneity
property of $\ell$, the extreme value copula $C$ can be rewritten as:

$\begin{split}&C\left(u_{1},\cdots,u_{d}\right)\\\
&=e^{\left(\sum_{k=1}^{d}\log u_{k}\right)A\left(\frac{\log
u_{1}}{\sum_{k=1}^{d}\log u_{k}},\cdots,\frac{\log u_{d}}{\sum_{k=1}^{d}\log
u_{k}}\right)},\end{split}$ (3)

where $A$ is known as the Pickands dependence function, which can be thought
of as the restriction of $\ell$ to the unit simplex
$\Delta_{d-1}=\\{\mathbf{w}=\left(w_{1},\cdots,w_{d}\right)\in[0,\infty)^{d}:\sum_{k=1}^{d}w_{k}=1\\}$.
The Pickands function $A$ is known to be convex and satisfies:

$\max_{1\leq k\leq d}w_{k}\leq A(w_{1},\cdots,w_{d})\leq 1\vspace{-0.2cm}$ (4)

for all $\mathbf{w}=\left(w_{1},\cdots,w_{d}\right)\in\Delta_{d-1}$. These
properties characterize Pickands dependence functions for $d=2$ but not
necessarily for general $d$ [3].

An equivalent formulation of the Pickands dependence function which proves to
be essential in the design of sampling algorithms (to be discussed in Section
4) is the integral representation given below.

######  Definition 2 (Integral form of the Pickands dependence function [15,
8]).

For any Pickands dependence function $A(\cdot,\cdots,\cdot)$, there exists a
Borel measure $\mu$ on $\Delta_{d-1}$ satisfying
$\int_{\Delta_{d-1}}s_{k}d\mu(\mathbf{s})=1$ for all $k\in\\{1,...,d\\}$ and
$\mathbf{s}=(s_{1},\cdots,s_{d})\in\Delta_{d-1}$ such that

$A(w_{1},...,w_{d})=\int_{\Delta_{d-1}}\max(w_{1}s_{1},...,w_{d}s_{d})d\mu(\mathbf{s}).$
(5)

Next, we provide details on how to model and train Pickands dependence
functions using input convex neural networks.

##  3 Modeling and training Pickands dependence functions with input convex
neural networks

###  3.1 Input convex neural networks

Since the Pickands dependence function $A$ is convex over $\Delta_{d-1}$, we
are interested in finding a flexible parameterization of $A$ that also
enforces convexity over the input domain. Input convex neural networks (ICNNs)
provide the convexity guarantee while also having the capability of
approximating convex functions in the $\sup$ norm [5, 6]. Using ICNNs, we
mitigate issues faced by previous estimators [11, 12, 13] in enforcing
convexity or projecting a nonparametric estimator onto a convex basis. Besides
convexity, we enforce (4) both by designing the architecture carefully and
adding an additional loss term. We provide further details on enforcing (4) in
the supplemental material. Followingly, we provide full details on how to
train the Pickands based ICNN from data.

###  3.2 Training the Pickands dependence function

After modeling the Pickands dependence function with ICNNs, the remaining
challenging task consists of training the model Pickands function
$A(\mathbf{w};\theta)$ to fit the data. One possible approach is to obtain the
underlying probability density function (PDF) and then train the ICNN with
maximum likelihood estimation (MLE). The drawback of such method is the need
to differentiate the $d-$variate CDF in order to obtain the corresponding PDF.
This may be computationally complex particularly in high dimensions. Let
$F_{k}$ denote the marginal CDF 111Fitting the marginals can be done using MLE
[1] or the $L-$moments method [16]. of the $k$th normalized component wise
maxima $\bar{M}_{n}^{(k)}=\frac{M_{n}^{(k)}-b_{n}^{(k)}}{a_{n}^{(k)}}$,
$k\in\\{1,\cdots,d\\}$ and let
$\mathbf{w}=\left(w_{1},\cdots,w_{d}\right)\in\Delta_{d-1}$. We perform the
following transformations on $\bar{M}_{n}^{(k)}$ (Please see Section 3 of
[14]):

$\displaystyle\widetilde{M}_{n}^{(k)}$ $\displaystyle=-\log
F_{k}(\bar{M}_{n}^{(k)}),\>\forall k\in\\{1,\cdots,d\\}.$ (6) $\displaystyle
Z_{w}$ $\displaystyle=\min_{1\leq k\leq d}\widetilde{M}_{n}^{(k)}/w_{k}.$ (7)

Then, the following relation holds:

$\displaystyle\mathbb{P}\left[Z_{w}>z\right]$
$\displaystyle=\mathbb{P}\left[\widetilde{M}_{n}^{(1)}>w_{1}z,\cdots,\widetilde{M}_{n}^{(d)}>w_{d}z\right]$
$\displaystyle=\mathbb{P}\left[F_{k}(\bar{M}_{n}^{(1)})<e^{-zw_{1}},\cdots,F_{k}(\bar{M}_{n}^{(d)})<e^{-zw_{d}}\right]$
$\displaystyle=e^{\left(-\sum_{k=1}^{d}zw_{k}\right)A\left(\frac{-zw_{1}}{-\sum_{k=1}^{d}zw_{k}},\cdots,\frac{-zw_{d}}{-\sum_{k=1}^{d}zw_{k}}\right)}$
$\displaystyle=e^{-zA(\mathbf{w})}.$

This means that the random variables $Z_{w}$ are exponentially distributed
with parameter given by the Pickands dependence function $A(\mathbf{w})$.

In light of above, we can learn the Pickands dependence function
$A(\mathbf{w})$ by fitting the model $A(\mathbf{w};\theta)$ to samples $Z_{w}$
using MLE. This can be done by training the ICNN modeling
$A(\mathbf{w};\theta)$ with stochastic gradient descent (SGD) to match the
data points $Z_{w}$ as follows:

$\displaystyle
A^{\star}(\mathbf{w})=\arg\min_{\theta}\mathbb{E}_{Z_{w}}\mathcal{L}(Z_{w};\theta),$
(8)

where $\mathcal{L}(Z_{w};\theta)=A(\mathbf{w};\theta)Z_{w}-\log
A(\mathbf{w};\theta)$. It should be noted that maximizing this likelihood with
the ICNN guarantees a convex estimate similar to the nonparametric estimator
defined in [4]. Alternative losses could additionally be considered, such as
reformulating the loss with respect to the estimators defined in [12, 11].
However, we have empirically found the MLE approach described in (8) gives
superior performance for most applications. Moreover, it follows naturally
from the original formulation of Pickands. Details of the training procedure
are summarized in Algorithm 1.

Algorithm 1 Training the Pickands dependence function with SGD (Pickands-ICNN)

1: Input:
$\left\\{\left(X_{i}^{(1)},\cdots,X_{i}^{(d)}\right)\right\\}_{i=1}^{N}$,
$N=B\times n$ samples of i.i.d random vectors where $B$ is the number of
blocks of data and $n$ is the size of each block.

2: Take block maxima of size $n$ and obtain component-wise maxima:
$\left\\{\left(M_{n,b}^{(1)},\cdots,M_{n,b}^{(d)}\right)\right\\}_{b=1}^{B}$
where $M_{n,b}^{(k)}=\max_{(b-1)n+1\leq i\leq bn}X_{i}^{(k)}$,
$(k,b)\in\\{1,\cdots,d\\}\times\\{1,\cdots,B\\}$.

3: Fit a GEV to each component-wise maxima $\\{M_{n,b}^{(k)}\\}_{b=1}^{B}$ and
obtain $\\{\bar{M}_{n,b}^{(k)}\\}_{b=1}^{B}$ then estimate marginals $F_{k}$
for each $k\in\\{1,\cdots,d\\}$.

4: Initialize the parameters $\theta$ of the ICNN Repeat:

5: Randomly sample a minibatch of training data
$\\{\bar{M}_{n,b}^{(k)}\\}_{b\in\text{batch}}$ and uniformly sample
$\mathbf{w}\in\Delta_{d-1}$.

6: Transform samples according to Equations (6) and (7) and obtain samples
$\\{Z_{w,b}\\}_{b\in\text{batch}}$.

7: Compute gradient
$\nabla_{\theta}\sum_{b\in\text{batch}}\mathcal{L}\left(Z_{w,b};\theta\right)$.

8: Update $\theta$ with Adam [17] Until convergence Output:
$A(\mathbf{w};\theta_{*})$.

The choice of hyperparameters such as the block size, number of blocks,
learning rate for Adam and batch size are provided in section 5 with settings
depending on the application.

###  3.3 Survival probability estimation

Following the training of the Pickands dependence function
$A(\mathbf{w};\theta_{*})$, our goal is to estimate survival probabilities.
More precisely, let
$\left(\gamma_{1},\cdots,\gamma_{d}\right)\in\mathbb{R}^{d}$ be a
$d-$dimensional vector of thresholds, we are interested in calculating the
following survival probability:

$\begin{split}&\mathbb{P}\left[M_{n}^{(1)}>\gamma_{1},\cdots,M_{n}^{(d)}>\gamma_{d}\right]=\mathbb{P}\left[\bar{M}_{n}^{(1)}>\bar{\gamma}_{1},\cdots,\bar{M}_{n}^{(d)}>\bar{\gamma}_{d}\right],\end{split}$
(9)

where $\bar{\gamma}_{k}=\frac{\gamma_{k}-b_{n}^{(k)}}{a_{n}^{(k)}}$,
$k\in\\{1,\cdots,d\\}$. Note that calculating the survival probability in (9)
using the model extreme value copula is not straightforward. To calculate
this, we use the following change-of-variable trick, which we prove in the
supplementary material.

######  Proposition 1.

Let $G_{k}(x):=F_{k}^{-1}(1-F_{k}(x))$ for $k~{}\in~{}\\{1,\dots,d\\}$, then
the random variables $G_{k}(\bar{M}_{n}^{(k)})$ and $\bar{M}_{n}^{(k)}$ have
the same marginal CDF $F_{k}$, for $k\in\\{1,\dots,d\\}$, and

$\displaystyle\mathbb{P}\left[\bar{M}_{n}^{(1)}>\bar{\gamma}_{1},\cdots,\bar{M}_{n}^{(d)}>\bar{\gamma}_{d}\right]=\mathbb{P}\left[G_{1}(\bar{M}_{n}^{(1)})<G_{1}(\bar{\gamma}_{1}),\cdots,G_{d}(\bar{M}_{n}^{(d)})<G_{d}(\bar{\gamma}_{d})\right].$
(10)

This proposition implies that the transformed variables
$G_{k}(\bar{M}_{n}^{(k)})$ are samples from extreme value distributions,
therefore we train a model Pickands dependence function on these transformed
variables, and finally evaluate the corresponding extreme value copula on
$\left(F_{1}(G_{1}(\bar{\gamma}_{1})),\cdots,F_{d}(G_{d}(\bar{\gamma}_{d}))\right)$.
Details on how to estimate the survival probability in (9) are summarized in
Algorithm 2.

Algorithm 2 Estimating survival probabilities with the Pickands dependence
function

1: Input: $\\{\bar{M}_{n,b}^{(k)}\\}_{b=1}^{B}$, thresholds:
$\left(\gamma_{1},\cdots,\gamma_{d}\right)$.

2: Train a model $A(\mathbf{w};\theta)$ with the transformed variables
$\\{(G_{1}(\bar{M}_{n,b}^{(1)},\dots,G_{d}(\bar{M}_{n,b}^{(d)}))\\}_{b=1}^{B}$
using Algorithm 1 and obtain $A(\mathbf{w};\theta_{*})$.

3: Evaluate the Pickands copula:  $\displaystyle
C\left(1-F_{1}(\bar{\gamma}_{1}),\cdots,1-F_{d}(\bar{\gamma}_{d})\right),$
where $C$ is calculated as in Equation (3) with $A=A(\mathbf{w};\theta_{*})$.

##  4 Sampling from multivariate extreme value distributions

While learning MEV distributions from data is important for estimating
survival probabilities and CDFs, it is also useful to simulate possible
scenarios by sampling from the estimated MEV distribution. We propose a method
for generating samples from an MEV distribution with an arbitrary Pickands
dependence function. To the best of our knowledge there are no general
sampling algorithms for general extreme value copula in arbitrary dimensions.
We exploit the exact sampling algorithms for the infinite dimensional analogue
of MEV distributions known as max-stable processes to sample from general
Pickands copulas [9, 10]. We additionally propose a heuristic for sampling
based on the relationship between max-stable processes and Pickands dependence
functions which we outline in Algorithm 4. To give a brief outline of how we
proceed, we begin by describing max-stable processes and the spectral
representation of stationary max-stable processes. We then describe the
relationship between Pickands functions and stationary max-stable processes
and use this relationship to describe a sampling algorithm. This ultimately
leads us to recast the MEV sampling in terms of previous work on sampling from
max-stable processes.

###  4.1 Pickands dependence and max-stable processes

####  4.1.1 Definition of max-stable processes

We initially review the definition of a max-stable process as defined in [7].
Let $X_{1}(t),X_{2}(t),...,X_{n}(t)$ be independent copies of a continuous
stochastic process on a compact set $\mathcal{T}$. If there exists normalizing
functions $a_{n}(t)>0$ and $b_{n}(t)\in\mathbb{R}$ for all $t\in\mathcal{T}$
such that the limit

$\frac{\max_{1\leq i\leq n}X_{i}(t)-b_{n}(t)}{a_{n}(t)}\to
M(t)\>\>\text{as}\>\>n\to\infty$

is non degenerate then the limiting process $M(t)$ is a max-stable process. A
key property of max-stable processes is that any univariate marginal is
distributed as a GEV distribution and any subset of marginals is distributed
as a MEV distribution. This motivates the relationship to MEV distributions
parameterized by Pickands dependence functions. We then consider this analysis
with respect to stationary max-stable processes and the ensuing spectral
decomposition of stationary max-stable processes.

####  4.1.2 Spectral decomposition of stationary max-stable processes

Stationary max-stable processes can be intuitively interpreted as i.i.d
samples from infinite dimensional extreme value distributions. A stationary
max-stable process can be decomposed by the spectral representation defined in
[7]. Suppose that $M(t)$ has unit Fréchet margins and is stationary. Then,
$M(t)$ can be written as:

$M(t)=\max_{i\geq 1}\xi_{i}Y_{i}^{+}(t),\quad t\in\mathcal{T}.$ (11)

$\\{Y_{i}(t)\\}_{i\geq 1}$ are i.i.d copies of a continuous stochastic process
$Y$ defined on $\mathcal{T}$ such that $\mathbb{E}[Y^{+}(t)]=1$ with
$Y^{+}(t)=\max\\{0,Y(t)\\}$ and $\xi_{i}$ is the $i$th realization of an
independent Poisson point process on $[0,\infty)$ with intensity
$\xi^{-2}d\xi$.

####  4.1.3 Margins of max-stable processes as MEV distributions

From the spectral decomposition of the max-stable process, we can recover a
new representation for the MEV distribution. Motivated by (5) and using the
results in [15, 8], the Pickands dependence function is related to the max
stable process through:

$A(\mathbf{w})=\mathbb{E}[\max_{1\leq k\leq
d}w_{k}Y^{+}(t_{k})],\>\>\mathbf{w}\in\Delta_{d-1}.$

Intuitively, this suggests that a $d-$dimensional set of samples of the max-
stable process are related to the Pickands dependence function. When
restricted to the discrete case of MEV distributions, we can rewrite
$\\{Y^{+}(t_{k})\\}_{k=1}^{d}$ as a $d-$dimensional vector
$\mathbf{y}=(y_{1},\cdots,y_{d})\in\mathbb{R}^{d}_{+}$.

$A(\mathbf{w})=\mathbb{E}_{\mathbf{y}}[\max_{0\leq k\leq
d}w_{k}y_{k}],\>\mathbf{w}\in\Delta_{d-1}.$ (12)

Next, we describe how to generate $\mathbf{y}$ from a given Pickands
dependence function. After learning the generative mechanism for $\mathbf{y}$,
we can use exact sampling algorithms to sample from the MEV distribution.

###  4.2 Learning the underlying generator for Pickands dependence

When given a Pickands dependence function $A$, we propose to recover the
underlying spectral representation by optimizing a generative model for
$\mathbf{y}$ in (12). We model $\mathbf{y}$ as the output of a function
$G(\cdot;\phi)\in\mathbb{R}^{d}_{+}$ with parameters $\phi$, i.e.
$\mathbf{y}=G(\mathbf{z};\phi)$ which maps input samples $\mathbf{z}\sim
p_{z}$ to $\mathbf{y}$, where $p_{z}$ is a distribution that is easy to sample
from (such as a multivariate Gaussian distribution). We then learn the
parameters $\phi$ of the generator by solving the following optimization
problem:

$\displaystyle\min_{\phi}\mathbb{E}_{\mathbf{w}\sim\Delta_{d-1}}$
$\displaystyle\left\|A(\mathbf{w})-\mathbb{E}_{\mathbf{y}}[\max_{0\leq k\leq
d}w_{k}y_{k}]\right\|_{2}^{2}+\left\|\mathbb{E}_{\mathbf{y}}[\mathbf{y}]-\mathbf{1}_{d}\right\|_{2}^{2},$
(13)

where $\mathbf{y}=(y_{1},\cdots,y_{d})=G(\mathbf{z};\phi)$,
$\mathbf{y}\in\mathbb{R}^{d}_{+}$ and $\mathbf{z}\in\mathbb{R}^{k}\sim p_{z}$.
The last expectation is required to enforce the condition of
$\mathbb{E}[Y(\cdot)]=1$ from [15]. The expectations with respect to
$\mathbf{y}$ in (13) are approximated using sample mean with samples from the
generator. Note that optimizing the objective in (13) only requires samples of
$\mathbf{y}$ rather than samples from the MEV distribution, thereby bypassing
additional complexities required by sampling from the full MEV distribution
while training. The full details of how to train the generator are outlined in
Algorithm 3.

Algorithm 3 Training a generator for a Pickands copula

1: Input: $A(\mathbf{w})$, $p_{z}$, tolerance parameter $\epsilon$

2: Initialize parameters $\phi$ of generator $G(\cdot;\phi)$

3: Sample $\\{\mathbf{w}^{(j)}\\}_{j=1}^{N_{\text{simplex}}}$ samples
uniformly over $\Delta_{d-1}$.

4: while
$\sum_{j=1}^{N_{\text{simplex}}}\mathcal{L}(\mathbf{w}^{(j)};\phi)>\epsilon$
do

5: Sample $\\{\mathbf{w}^{(j)}\\}_{j=1}^{N_{\text{simplex}}}$ samples
uniformly over $\Delta_{d-1}$.

6: Sample $\\{\mathbf{y}^{(i)}\\}_{i=1}^{N_{\text{gen}}}$ where
$\mathbf{y}^{(i)}=G(\mathbf{z}^{(i)};\phi),\mathbf{z}^{(i)}\sim p_{z}$ for
$1\leq i\leq N_{\text{gen}}$.

7: Define $\eta(\mathbf{w},\mathbf{y})=\max\\{\mathbf{w}\odot\mathbf{y}\\}$
with $\odot$ denotes the point-wise multiplication.

8: Compute gradient w.r.t $\phi$ of
$\sum_{j=1}^{N_{\text{simplex}}}\mathcal{L}(\mathbf{w}^{(j)};\phi)$ where:
$\displaystyle\mathcal{L}(\mathbf{w}^{(j)};\phi)$
$\displaystyle=\left\|A(\mathbf{w}^{(j)})-\frac{1}{N_{\text{gen}}}\sum_{i=1}^{N_{\text{gen}}}\eta(\mathbf{w}^{(j)},\mathbf{y}^{(i)})\right\|_{2}^{2}+\|\frac{1}{N_{\text{gen}}}\sum_{i=1}^{N_{\text{gen}}}\mathbf{y}^{(i)}-1\|_{2}^{2}$

9: Update $\phi$ using Adam [17].

10: end while

11: Output: $G(.;\phi_{*})$.

Finally, to generate samples from the full MEV distribution, we can consider
methods that only require the generator, such as the ones described in [10,
9]. For simplicity, we provide a heuristic inspired by the method in [9] to
provide a straightforward procedure of sampling which we describe in Algorithm
4. This approach involves repeated sampling from the generator and unit
Fréchet distributions and taking the maximum over the samples. This can be
seen as a simulation of (11) since unit Frechét is equivalent to
$1/\xi,\xi\sim\text{Exp}(1)$ which describes a Poisson process with rate
$\xi^{-2}d\xi$.

Algorithm 4 Heuristic for sampling from a Pickands copula

1: Input: $A(\mathbf{w})$

2: Optimize a generator $G(\cdot;\phi)$ using Algorithm 3.

3: for $i\in\\{1,\ldots,N_{\text{gen}}\\}$ do

4: Generate $\mathbf{y}^{(i)}$ where
$\mathbf{y}^{(i)}=G(\mathbf{z}^{(i)};\phi_{*}),\mathbf{z}^{(i)}\sim p_{z}$.

5: Sample $\xi^{(i)}$ from a unit Fréchet distribution.

6: end for

7: Compute the component-wise maxima as: $M=\max_{1\leq i\leq
N_{\text{gen}}}\\{\xi^{(i)}\odot\mathbf{y}^{(i)}\\}$.

8: Output: $M$.

##  5 Results

In this section, we provide numerical results that compare the performance of
the proposed ICNN-based model for the Pickands dependence function with some
well-known estimators from the literature [11, 12, 4]. We start by evaluating
the performance for estimating the Pickands function for some bi-variate
symmetric and asymmetric logistic models. We then consider higher dimensional
results for the same models. Finally, we consider some benchmark datasets such
as Ozone data discussed in [18], wind gust data, and historical prices of
commodities. For evaluating the efficacy of sampling from a MEV, we consider
experiments for different dimensions and dependence with data generated from
the symmetric and asymmetric logistic models. For comparison, we consider the
following estimators: Pickands [4], CFG [12], and the estimator described in
[11] which we refer to as BDV, for which we give brief description in the
supplementary material. Throughout this section, we refer to the proposed
ICNN-based estimator as Pickands-ICNN. Details regarding properties of ICNNs
can be found in [5], and additional details on architectural choices for all
experiments are given in the supplementary material.

###  5.1 Synthetic data

(a) $A_{\text{SL}}$ MSE ($d=2$)

(b) $A_{\text{ASL}}$ MSE ($d=2$)

(c) $A_{\text{SL}}$ MSE ($\alpha=0.5$)

(d) $A_{\text{ASL}}$ MSE ($\alpha=0.5$)

Figure 1: (1(a), 1(b)) MSE of survival probabilities for $d=2$ with $100$
samples for $A_{\text{SL}}$ (1(a)) and $A_{\text{ASL}}$ (1(b)). Thresholds are
above the $75$th percentile. (1(c), 1(d)) Comparison of
$||\hat{A}(w)-A(w)||_{2}^{2}$ for different estimators $\hat{A}$ for different
$\alpha=0.5$ and for dimensions $d=\\{16,64,256,512,1024\\}$ for
$A_{\text{SL}}$ (1(c)) and $A_{\text{ASL}}$ (1(d)). Results are averages of 5
runs.

We consider some canonical families of extreme value distributions known as
symmetric logistic family where the underlying Pickands function is given by
[3]:

$\displaystyle
A_{\text{SL}}(\mathbf{w})=\left(\sum_{k=1}^{d}w_{k}^{1/\alpha}\right)^{\alpha},\>\mathbf{w}\in\Delta_{d-1},$
(14)

where $\alpha\in(0,1]$ is a parameter modeling the degree of dependence
between variables ranging from complete dependence ($\alpha=0$) to complete
independence ($\alpha=1$). We also consider a generalization of this copula,
the asymmetric logistic family:

$\displaystyle
A_{\text{ASL}}(\mathbf{w})=\sum_{b\in\mathcal{P}_{d}}\bigg{(}\sum_{i\in
b}(\lambda_{i,b}w_{i})^{1/\alpha_{b}}\bigg{)}^{\alpha_{b}},$ (15)

where $\mathcal{P}_{d}$ is the power set of $\\{1,...,d\\}$. The dependence
parameters $\alpha_{b}\in(0,1]$ are defined for $b\in\mathcal{P}_{d}$ except
for all singleton sets. The asymmetry parameters $\lambda_{i,b}\in[0,1]$ are
randomly generated to satisfy: $\sum_{b\in\\{c\in\mathcal{P}_{d}\>:\>i\in
c\\}}\lambda_{i,b}=1$. Exact sampling of distributions of this type are
described in [19]. Note that for both the symmetric and asymmetric copulas,
the marginals are distributed according to the standard Fréchet distribution.
In this case, we skip Steps 2 and 3 of Algorithm 1 and consider empirical
marginals given by the ranks of different samples.

#### Estimation of survival probabilities

We start by comparing the MSE of survival probabilities for $d=2$ where the
true Pickands dependence function is given by the symmetric or asymmetric
model described respectively in (14) and (15) for different degrees of
dependence $\alpha$. We compute the exact values of the survival probability
and consider survival probabilities associated with margins above the 75th
percentile. As shown in Figures 1(a) and 1(b), the proposed Pickands-ICNN
estimator achieves the lowest MSE performance for almost all degrees of
dependence $\alpha$ for both the the symmetric and asymmetric logistic models.
This is due in part to the flexibility of the proposed Pickands-ICNN and the
enforcing of convexity through the ICNN.

#### Discrepancy in estimated Pickands dependence functions

We additionally showcase the ability of the proposed method to model high
dimensional extreme value distributions. To do this, we train the Pickands-
ICNN with data for dimensions $\\{16,64,256,512,1024\\}$ and compare the
estimated $\hat{A}$ to the true $A$. Specifically, we compute the mean squared
error between the true Pickands dependence and the estimated one for 1000
uniformly sampled points in $\Delta_{d-1}$. The results are illustrated for
$\alpha=0.5$ in Figures 1(d) and 1(c). As shown by Figures 1(d) and 1(c), the
proposed Pickands ICNN-based estimator achieves the lowest MSE performance for
most values of the dimension $d$. The results suggest the proposed method is
effective in recovering the Pickands dependence for high dimensional data.

###  5.2 Real data

In order to determine the efficacy of the proposed method on real data, we
consider data on extreme ozone levels $(d=4)$, wind gusts $(d=10)$, and
commodity prices $(d=3)$. Details for each dataset are in the supplemental
materials. The main challenge associated to real data is the lack of a ground
truth Pickands function for comparison purposes. We therefore resort to
computing the accuracy of the model estimate w.r.t the empirical estimate over
a series of thresholds. This is quantified as:

$\frac{1}{|Q|}\sum_{\mathbf{\gamma}\in
Q}\left[\frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}_{\\{M_{n,b}\geq\mathbf{\gamma}\\}}-P_{\theta}(M_{n}\geq\mathbf{\gamma})\right]^{2},$
(16)

where $M_{n,b}=\left(M_{n,b}^{(1)},\cdots,M_{n,b}^{(d)}\right)$ is the
$d-$dimensional vector of point-wise maxima and $P_{\theta}$ is the estimated
survival probability computed using Algorithm 2 and $Q$ is a set of thresholds
to consider. In general, we are interested in thresholds corresponding to
extremes in the tail; however, due to the lack of data in that region, we
estimate accuracy in areas where more data observations are available.

| Pickands | CFG | BDV | Proposed
---|---|---|---|---
Ozone | $3.16(\pm 4.27)\text{\footnotesize$\times 10$}^{-3}$ | $1.59(\pm 3.21)\text{\footnotesize$\times 10$}^{-1}$ | $2.98(\pm 4.05)\text{\footnotesize$\times 10$}^{-3}$ | $\bf 1.54(\pm 2.04)\text{\footnotesize$\times 10$}^{-3}$
Wind | $2.57(\pm 4.53)\text{\footnotesize$\times 10$}^{-3}$ | $7.31(\pm 21.8)\text{\footnotesize$\times 10$}^{-2}$ | $8.02(\pm 14.0)\text{\footnotesize$\times 10$}^{-4}$ | $\bf 7.44(\pm 12.9)\text{\footnotesize$\times 10$}^{-4}$
Commodities | $3.12(\pm 4.68)\text{\footnotesize$\times 10$}^{-3}$ | $3.12(\pm 4.64)\text{\footnotesize$\times 10$}^{-3}$ | $\bf 2.97(\pm 4.26)\text{\footnotesize$\times 10$}^{-3}$ | $3.03(\pm 3.92)\text{\footnotesize$\times 10$}^{-3}$
Table 1: MSE performance on the test set for different estimators and
datasets. All data are for the 75th percentile and above.

(a) CFG for generated vs real samples (SL)

(b) CFG for generated vs real samples (ASL)

Figure 2: MSE (solid) and relative error (dashed) of CFG estimate for 1000
samples and 1000 simplex points for $d=225$ at various $\alpha\in(0,1)$ for
$A_{\text{SL}}$ (2(a)) and $A_{\text{ASL}}$ (2(b)) for learned sampled data
(red) and exact sampled (blue).

###  5.3 Sampling from the copula

To determine the efficacy of sampling from an arbitrary Pickands copula, we
consider two synthetic examples using the previously described analytical MEV
distributions. In this experiment, we train the generator $G(.;\phi)$ in (13)
based on the given Pickands dependence function. To determine the efficacy of
the sampling algorithm, we estimate the Pickands dependence function using the
nonparameteric CFG estimator with empirical marginals and compare to the CFG
estimate of samples generated through the exact simulation method. We consider
the CFG estimator due to its ubiquity in the literature and its well regarded
status as a standard estimator for the Pickands dependence function. The
results for generating 225 dimensional samples are shown in figures 2(a) and
2(b). The figures suggest that our sampling algorithm performs comparatively
well for both distributions considered, with the worst performance occurring
in the nearly independent cases ($\alpha=1$).

(a) SL CFG MSE

(b) ASL CFG MSE

Figure 3: Comparison of MSE (solid) and relative error (dashed) for 1000
samples and 1000 simplex points of CFG estimate for Pickands function of
learned sampled data (red) and exact sampled (blue). Figures are for
$d=\\{2,16,256,784,1024\\}$ dimensional samples with $\alpha=0.5$.

Furthermore, the error of the CFG estimate for the proposed sampling method
(red) and the exact sampling (blue) follow very similar trends in errors,
suggesting the proposed sampling method is effectively representing the true
distribution. The experiments are then repeated where we vary the dimension of
the samples and compare with real data samples in Figures 3(a) and 3(b). The
dependence parameter $\alpha$ is fixed to be $\alpha=0.5$ for both
experiments. The results suggest the proposed sampling method is comparable to
the exact sampling method, with certain dimensions performing better than
others. This is likely due to the need for more careful choices of
hyperparameters for different dimensions, but we leave that analysis as future
work.

##  6 Concluding remarks

In this paper, we introduced a new method for modeling Pickands dependence
functions for arbitrary dimensions using input convex neural networks. This
allows for better representations of the underlying MEV distribution due to
the representational power of ICNNs. We additionally present a method for
recovering the spectral representation of MEV distributions, which makes it
possible to create a generative model that samples from the underlying MEV
distribution. Numerical results have been provided to empirically demonstrate
the effectiveness of the proposed methods in their respective tasks. However,
it is worth mentioning that these methods can still be further be improved by
considering the following limitations.

### Limitations of Pickands-ICNN

ICNNs can guarantee convexity in approximating Pickands dependence functions
for high dimensions while maintaining a reasonable number of parameters.
However, it remains a challenge to optimize the parameters using standard deep
learning techniques. Specifically, choosing appropriate hyperparameters is a
difficult and opaque task that requires additional care. This is an attribute
where nonparametric methods are beneficial, at the cost of being unable to
guarantee convexity. Additional progress on understanding properties of
training deep neural networks should improve the representational capabilities
of the ICNN given its theoretical potential to approximate convex functions to
arbitrary precision.

### Future Work

The proposed methods have possible applications in a variety of modeling
situations. For instance, understanding the spectral measure is important for
finding groups of variables that are extreme simultaneously [8]. Therefore, an
important but non-trivial extension of the present work lies in capturing
clusters of variables and understanding the joint interaction between them.
Finally, applications of high dimensional extremes are important in
understanding robustness properties of neural networks [20]. The proposed work
provides foundation for such extensions.

## Acknowledgements

Material in this paper is based upon work supported by the Air Force Office of
Scientific Research under award number FA9550-20-1-0397. AH was supported by
NSF Graduate Research Fellowship.

## References

  * [1] Paul Embrechts, Thomas Mikosch, and Claudia Klüppelberg.  Modelling Extremal Events: For Insurance and Finance.  Springer-Verlag, Berlin, Heidelberg, 1997.
  * [2] Laurens de Haan and Ana Ferreira.  Extreme Value Theory: An Introduction (Springer Series in Operations Research and Financial Engineering).  Springer, 2010.
  * [3] Gordon Gudendorf and Johan Segers.  Extreme-value copulas.  In Copula Theory and Its Applications. Springer Berlin Heidelberg, 2010.
  * [4] J Pickands.  Multivariate extreme value distributions, bull. int. statist.  1981\.
  * [5] Brandon Amos, Lei Xu, and J Zico Kolter.  Input convex neural networks.  In International Conference on Machine Learning, pages 146–155. PMLR, 2017.
  * [6] Yize Chen, Yuanyuan Shi, and Baosen Zhang.  Optimal control via neural networks: A convex approach.  arXiv preprint arXiv:1805.11835, 2018.
  * [7] Laurens De Haan et al.  A spectral representation for max-stable processes.  The annals of probability, 12(4):1194–1204, 1984.
  * [8] Johan Segers.  Max-stable models for multivariate extremes.  arXiv preprint arXiv:1204.0332, 2012.
  * [9] Clément Dombry, Sebastian Engelke, and Marco Oesting.  Exact simulation of max-stable processes.  Biometrika, 103(2):303–317, 2016.
  * [10] Zhipeng Liu, Jose H Blanchet, AB Dieker, and Thomas Mikosch.  On optimal exact simulation of max-stable and related random fields.  arXiv preprint arXiv:1609.06001, 2016.
  * [11] Axel Bücher, Holger Dette, Stanislav Volgushev, et al.  New estimators of the pickands dependence function and a test for extreme-value dependence.  The Annals of Statistics, 39(4):1963–2006, 2011.
  * [12] Philippe Capéraà, A-L Fougères, and Christian Genest.  A nonparametric estimation procedure for bivariate extreme value copulas.  Biometrika, 84(3):567–577, 1997.
  * [13] Giulia Marcon, SA Padoan, Philippe Naveau, Pietro Muliere, and Johan Segers.  Multivariate nonparametric estimation of the pickands dependence function using bernstein polynomials.  Journal of Statistical Planning and Inference, 183:1–17, 2017.
  * [14] Gordon Gudendorf and Johan Segers.  Nonparametric estimation of multivariate extreme-value copulas.  Journal of Statistical Planning and Inference, 142(12):3073 – 3085, 2012.
  * [15] Johan Segers.  Nonparametric inference for max-stable dependence.  Statistical Science, 27(2):193–196, 2012.
  * [16] J. R. M. Hosking.  L-moments: Analysis and estimation of distributions using linear combinations of order statistics.  Journal of the Royal Statistical Society. Series B (Methodological), 52(1):105–124, 1990.
  * [17] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In http://arxiv.org/abs/1412.6980, 2014.
  * [18] Javier Rojo Jiménez, Enrique Villa-Diharce, and Miguel Flores.  Nonparametric estimation of the dependence function in bivariate extreme value distributions.  Journal of Multivariate Analysis, 76(2):159 – 191, 2001.
  * [19] Alec Stephenson.  Simulating multivariate extreme value distributions of logistic type.  Extremes, 6(1):49–59, 2003.
  * [20] Tsui-Wei Weng, Huan Zhang, Pin-Yu Chen, Jinfeng Yi, Dong Su, Yupeng Gao, Cho-Jui Hsieh, and Luca Daniel.  Evaluating the robustness of neural networks: An extreme value theory approach.  arXiv preprint arXiv:1801.10578, 2018.
  * [21] Gordon Gudendorf and Johan Segers.  Nonparametric estimation of an extreme-value copula in arbitrary dimensions.  Journal of multivariate analysis, 102(1):37–47, 2011.

##  Appendix A Proof of Proposition 1

With the change-of-variable $G_{k}(x):=F_{k}^{-1}(1-F_{k}(x))$ for
$k\in\\{1,\cdots,d\\}$, it easily follows that the random variables
$G_{k}(\bar{M}_{n}^{(k)})\sim F_{k}$. Moreover, the survival probability can
be written as:

$\displaystyle\mathbb{P}\left[\bar{M}_{n}^{(1)}>\bar{\gamma}_{1},\cdots,\bar{M}_{n}^{(d)}>\bar{\gamma}_{d}\right]$
$\displaystyle=\mathbb{P}\left[1-F_{1}(\bar{M}_{n}^{(1)})<1-F_{1}(\bar{\gamma}_{1}),\cdots,1-F_{d}(\bar{M}_{n}^{(d)})<1-F_{d}(\bar{\gamma}_{d})\right]$
$\displaystyle=\mathbb{P}\left[G_{1}(\bar{M}_{n}^{(1)})<G_{1}(\bar{\gamma}_{1}),\cdots,G_{d}(\bar{M}_{n}^{(d)})<G_{d}(\bar{\gamma}_{d})\right]$
$\displaystyle=C\left(1-F_{1}(\bar{\gamma}_{1}),\cdots,1-F_{d}(\bar{\gamma}_{d})\right),$

where $C$ is the copula of
$\left(G_{1}(\bar{M}_{n}^{(1)}),\cdots,G_{d}(\bar{M}_{n}^{(d)})\right)$.

##  Appendix B Pickands, CFG and BDV estimators

#### Pickands Estimator

The Pickands estimator [4] is built following the transformations (6) and (7)
in the paper. The estimator is obtained by exactly maximizing the likelihood
(Equation (8) in the paper) resulting in the following non-parametric
estimate:

$\displaystyle\widehat{A}_{\text{Pickands}}(\mathbf{w})=\left(\frac{1}{B}\sum_{i=1}^{B}Z_{w,i}\right)^{-1}.$
(17)

#### CFG Estimator

The CFG estimator [12] is constructed following the observation:

$\displaystyle\mathbb{E}\log Z_{w}=-\log A(\mathbf{w})-\gamma,$

where $\gamma=-\int_{0}^{\infty}\log xe^{-x}dx$ denotes the Euler’s contant.
The CFG estimator is thus given by:

$\displaystyle\widehat{A}_{\text{CFG}}(\mathbf{w})=\exp\left[-\gamma-\frac{1}{B}\sum_{i=1}^{B}\log
Z_{w,i}\right].$ (18)

In our main submission we use a similar estimator, with the correction term
presented in [21]:

$\widehat{A}_{\text{CFG,C}}(\mathbf{w})=\exp\left(\log\widehat{A}_{\text{CFG}}(\mathbf{w})-\sum_{k=1}^{d}w_{k}\log\left(\widehat{A}_{\text{CFG}}(\mathbf{e}_{k})\right)\right),$
(19)

where $\mathbf{e}_{k}$ is the $k$-th canonical basis vector.

#### BDV Estimator

We propose an $d$-dimensional extension to the bivariate estimator described
in [11]. We begin by defining the minimum distance estimator between the true
CDF, $C(\mathbf{u})$ and the one estimated by the Pickands function
$A(\mathbf{w})$.

$\displaystyle\int_{[0,1]^{d}}\left[\log C(\mathbf{u})-\sum_{k=1}^{d}\log
u_{k}A\left(\frac{\log(\mathbf{u})}{\sum_{k}\log
u_{k}}\right)\right]^{2}\,d\mathbf{u}$ (20)
$\displaystyle=\int_{\Delta_{d-1}}\int_{0}^{1}(\log
C(y^{w_{1}},\dots,y^{w_{d}})-\log(y)A(\mathbf{w}))^{2}(-\log(y))^{d-1}\,dy\,d\mathbf{w}.$
(21)

We have

$\displaystyle\hat{C}(y^{w_{1}},\dots,y^{w_{d}})$
$\displaystyle=\frac{1}{B}\sum_{i=1}^{B}\mathbf{1}(F_{1}(\bar{M}_{n,i}^{(1)})\leq
y^{w_{1}},\cdots,F_{d}(\bar{M}_{n,i}^{(d)})\leq y^{w_{d}})$ (22)
$\displaystyle=\frac{1}{B}\sum_{i=1}^{B}\mathbf{1}(F_{1}(\bar{M}_{n,i}^{(1)})^{\frac{1}{w_{1}}}\leq
y,\cdots,F_{d}(\bar{M}_{n,i}^{(d)})^{\frac{1}{w_{d}}}\leq y)$ (23)
$\displaystyle=\frac{1}{B}\sum_{i=1}^{B}\mathbf{1}(\max_{1\leq k\leq
d}F_{k}(\bar{M}_{n,i}^{(k)})^{\frac{1}{w_{k}}}\leq y)$ (24)
$\displaystyle=\frac{1}{B}\sum_{i=1}^{B}\mathbf{1}\left(\Gamma_{w,i}\leq
y\right),$ (25)

where $\Gamma_{w,i}=\exp(-Z_{w,i})$. Now, if we reorder these so that
$\Gamma_{w,1}\leq\cdots\leq\Gamma_{w,B}$, we have that

$\hat{C}(y^{w_{1}},\dots,y^{w_{d}})=\left\\{\begin{aligned}
0&\;\>\text{if}\;\>y<\Gamma_{w,1},\\\
\frac{i}{B}&\;\>\text{if}\;\>\Gamma_{w,i}\leq
y<\Gamma_{w,i+1},\,i\in\\{1,\dots,B-1\\},\\\
1&\;\>\text{if}\;\>\Gamma_{w,B}\leq y.\end{aligned}\right.$ (26)

Because $\log\hat{C}(\cdots)$ is not defined if $y<\Gamma_{w,1}$, the
following modified estimator is considered in [11].

$\tilde{C}(y^{w_{1}},\dots,y^{w_{d}}):=\max\left\\{C(y^{w_{1}},\dots,y^{w_{d}}),B^{-\gamma}\right\\},$
(27)

where $\gamma$ is any positive real greater or equal than $\frac{1}{2}$. For
convenience, we choose $\gamma=1$ so that:

$\tilde{C}(y^{w_{1}},\dots,y^{w_{d}})=\left\\{\begin{aligned}
\frac{1}{B}&\;\>\text{if}\;\>y<\Gamma_{w,2},\\\
\frac{i}{B}&\;\>\text{if}\;\>\Gamma_{w,i}\leq
y<\Gamma_{w,i+1},\,i\in\\{2,\dots,B-1\\},\\\
1&\;\>\text{if}\;\>\Gamma_{w,B}\leq y.\end{aligned}\right.$ (28)

Finally, as in [11], for any positive weight function
$h:(0,1)\to\mathbb{R}_{0}^{+}$, let $h^{*}(y):=h(y)(\log y)^{2}$,

$B_{h}:=\int_{0}^{1}h^{*}(y)\,dy\quad\text{and}\quad
g(x):=-B_{h}^{-1}\int_{0}^{x}\frac{h^{*}(y)}{\log y}\,dy.$ (29)

Then, letting $\Gamma_{w,0}=0$, $\Gamma_{w,B+1}=1$, we define the BDV
estimator $\widehat{A}_{\text{BDV,}h}$ as follows

$\displaystyle\widehat{A}_{\text{BDV,}h}(\mathbf{w})$
$\displaystyle=B_{h}^{-1}\int_{0}^{1}\frac{\log\tilde{C}(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}h^{*}(y)\,dy$ (30)
$\displaystyle=B_{h}^{-1}\sum_{i=0}^{B}\int_{\Gamma_{w,i}}^{\Gamma_{w,i+1}}\frac{\log\tilde{C}(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}h^{*}(y)\,dy$ (31)
$\displaystyle=-\log\frac{1}{n}g(\Gamma_{w,2})-\sum_{i=2}^{n}\log\frac{i}{n}\left(g(\Gamma_{w,i+1})-g(\Gamma_{w,i})\right)$
(32)
$\displaystyle=-\sum_{i=2}^{n}\log\frac{i-1}{n}g(\Gamma_{w,i})+\sum_{i=2}^{n}\log\frac{i}{n}g(\Gamma_{w,i})$
(33)
$\displaystyle=\sum_{i=2}^{n}\log\left(1+\frac{1}{i-1}\right)g(\Gamma_{w,i})$
(34)

In our main submission, we use a slightly modified estimator, which proved to
have superior performance in our experiments. Recall that, if $A$ is a
Pickand’s dependence function, we have $\max(\mathbf{w})\leq A(\mathbf{w})\leq
1$, which implies the true copula verifies:

$\max(\mathbf{w})\leq\frac{\log C(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}=A(\mathbf{w})\leq 1.$ (35)

Accordingly, we let

$\operatorname{clamp}_{a,b}(x):=\left\\{\begin{array}[]{rl}a&\text{if }x\leq
a,\\\ x&\text{if }a<x<b,\\\ b&\text{if }x\geq b,\end{array}\right.$ (36)

and define

$\breve{C}(y^{w_{1}},\dots,y^{w_{d}})=\exp\left(\operatorname{clamp}_{\log
y,\max(\mathbf{w})\log y}\log\hat{C}(y^{w_{1}},\dots,y^{w_{d}})\right),$ (37)

and

$\displaystyle\widehat{A}_{\text{BDV,MM,}h}(\mathbf{w})$
$\displaystyle=B_{h}^{-1}\int_{0}^{1}\frac{\log\breve{C}(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}h^{*}(y)\,dy$ (38)
$\displaystyle=B_{h}^{-1}\sum_{i=0}^{B}\int_{\Gamma_{w,i}}^{\Gamma_{w,i+1}}\frac{\log\breve{C}(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}h^{*}(y)\,dy$ (39)

Letting
$\Gamma_{w,i}^{(\ell)}=\operatorname{clamp}_{\Gamma_{w,i},\Gamma_{w,i+1}}\left(\left(\frac{i}{n}\right)^{\frac{1}{\max(\mathbf{w})}}\right)$,
$\Gamma_{w,i}^{(u)}=\operatorname{clamp}_{\Gamma_{w,i},\Gamma_{w,i+1}}\left(\frac{i}{n}\right)$
for $i\in\\{0,\dots,B\\}$ and $\eta(x)=B_{h}^{-1}\int_{0}^{x}h^{*}(y)\,dy$, we
have

$\displaystyle\widehat{A}_{\text{BDV,MM,}h}(\mathbf{w})$
$\displaystyle=B_{h}^{-1}\sum_{i=0}^{B}\int_{\Gamma_{w,i}}^{\Gamma_{w,i+1}}\frac{\log\breve{C}(y^{w_{1}},\dots,y^{w_{d}})}{\log
y}h^{*}(y)\,dy$ (40)
$\displaystyle=B_{h}^{-1}\sum_{i=0}^{B}\int_{\Gamma_{w,i}}^{\Gamma_{w,i}^{(\ell)}}\max(\mathbf{w})h^{*}(y)\,dy+\int_{\Gamma_{w,i}^{(\ell)}}^{\Gamma_{w,i}^{(u)}}\log\frac{i}{n}\frac{h^{*}(y)}{\log
y}\,dy+\int_{\Gamma_{w,i}^{(u)}}^{\Gamma_{w,i+1}}h^{*}(y)\,dy$ (41)
$\displaystyle=\sum_{i=0}^{B}\max(\mathbf{w})\left(\eta(\Gamma_{w,i}^{(\ell)})-\eta(\Gamma_{w,i})\right)-\log\frac{i}{n}\left(g(\Gamma_{w,i}^{(u)})-g(\Gamma_{w,i}^{(\ell)})\right)+\eta(\Gamma_{w,i+1})-\eta(\Gamma_{w,i}^{(u)})$
(42)

In our main submission, we use $\widehat{A}_{\text{BDV,MM,}h}$ with
$h(y)=\frac{1}{\log(y)}$.

##  Appendix C Data Description

#### Synthetic data

For the synthetic data experiments we consider samples of 100 points from each
respective distribution. We use the full dataset for the batch size during
training. We additionally sample 1000 points from the simplex for each data
point during training.

#### Ozone data

We consider ozone levels measured at 4 different stations in Sequoia National
Park from data that can be downloaded from the National Park Service official
website 222https://ard-request.air-resource.com/data.aspx. The 4 stations are
located at Ash Mountain, Lower Kaweah, Grant Grove and Lookout Point. We train
the different models on daily maxima of ozone levels at the 4 different
stations for the period from January 1984 to December 1996. To encounter the
effect of seasonality, we do not train over the whole period, but we train
different models on a single month (training month, e.g. April of each year)
and compute accuracy on the consecutive month (validation month e.g. May of
the same year). The accuracy is averaged with the specific validation month of
each year over the whole period. For the experiments, we consider the
following pair of (training/validation) months: (April/May), (May/June),
(June/July), and (July/August).

#### California wind data

Figure 4: Locations (yellow circles) of weather stations sampled for the wind
speed experiments. Figure generated via Google Maps.

We are interested in modeling the extremal relationship of wind gusts between
different locations in California. We consider 10 locations in California
illustrated in 4. We obtained the data from the Remote Automated Weather
Station (RAWS) archive available at at the online
repository333https://raws.dri.edu/index.html. The RAWS data are collected from
various time intervals from December 1989 to December 2020. We consider only
the time points that occur in the intersection of all the data collected and
where all values are valid (i.e. not NaNs or missing). Similarly to the ozone
data, and in an effort to reduce seasonality, we consider the daily max wind
gust for the different locations for a single month over all the years the
data were collected. To evaluate the proposed method, we train and test on
data from consecutive months and repeat for multiple sets of months in our
dataset. (training/validation months): (June/July), (July/August),
(August/September), and (September/October).

#### Commodities data

We consider the extreme dependency between different commodities such as
Copper, Nickel and Zinc. We collect data of daily prices of the different
commodities from January 2010 to December 2020 as published in
444https://www.investing.com/commodities/. For training, we consider monthly
maxima of daily prices log returns over a period of two consecutive years. We
validate the performance by evaluating accuracy over next three years. For
example we consider the following pairs of [(training years),(validation
years)]: [(2010, 2011), (2012, 2013, 2014)], [(2012, 2013), (2014, 2015,
2016)], [(2014, 2015), (2016, 2017, 2018)], [(2016, 2017), (2018, 2019,
2020)].

##  Appendix D Further details on experiments

#### Architecture Details

For learning the Pickands dependence function, in all real data experiments we
used 24 width and 3 depth ICNNs. For synthetic data, we used 16 width and 4
depth ICNNs. All activations used were LeakyReLU activation functions. We
initialize the Pickands dependence function to be the output of the ICNN, i.e.
$A(\mathbf{w};\theta)=\text{ICNN}(\mathbf{w};\theta)$ and then force it to be
1 on the corners. To accomplish that, let $\mathbf{e}_{k}$ be the $k$th
$d-$dimensional canonical basis vector. We then set $\mathbf{b}$ to be the
vector with $k$th component $b_{k}=\text{ICNN}(\mathbf{e}_{k};\theta)$ for
$k\in\\{1,\cdots,d\\}$ and set
$A(\mathbf{w};\theta)=\text{ICNN}(\mathbf{w};\theta)-\langle\mathbf{w},\mathbf{b}\rangle+1$.
Therefore, when $\mathbf{w}=\mathbf{e}_{k}$,
$A(\mathbf{w};\theta)=\text{ICNN}(\mathbf{e}_{k};\theta)-b_{k}+1=1$, $\forall
k\in\\{1,\cdots,d\\}$. We finally apply a ReLU activation at the end to ensure
all outputs are positive.

For the sampling experiments, we vary the size of $p_{z}$ depending on the
dimensionality of the data. The latent sizes were $\\{2,16,64,64,64,256\\}$
for $d=\\{2,16,225,256,784,1024\\}$ respectively. The generator is a basic
multi layer perceptron (MLP) with ReLU activations. For all experiments, we
use a width 256 and depth 2 MLP for the generator. The final output is ensured
to be positive through a final ReLU operation.

#### Hyperparameter Tuning

For learning the Pickands dependence experiments, we used the Adam [17]
optimizer for optimizing all parameters with learning rate $1\times 10^{-4}$
with exponential decay of $0.9998$ per epoch. Each model was trained for 1000
epochs.

For the sampling experiments, the generator was trained using Adam with
learning rate $1\times 10^{-4}$, $\beta_{1}=0.5$ and $\beta_{2}=0.99$ with
exponential decay on the learning rate of $0.99998$. Models for the generator
were trained for 3000 epochs

##  Appendix E Sampling examples

(a) $\alpha=0$

(b) $\alpha=1$

(c) $\alpha=0.5$

Figure 5: Comparison of samples for $A=A_{SL}$. Blue points are from the
proposed method and orange points are true samples.

In Figure 5 we show examples of samples from the symmetric logistic
distribution for different values of the dependence parameter $\alpha$.
Samples from our heuristic and learned generator (blue) are compared to the
exact samples (orange) provided by the method in [19].