# Modelling and simulating spatial extremes by combining extreme value theory
with generative adversarial networks

Younes Boulaguiem1111Corresponding author<EMAIL_ADDRESS>Jakob
Zscheischler2,3,4 Edoardo Vignotto1 Karin van der Wiel5 Sebastian Engelke1

(1Geneva School of Economics and Management, University of Geneva, Boulevard
du Pont d’Arve 40, Geneva, 1205, Switzerland
2Department of Computational Hydrosystems, Helmholtz Centre for Environmental
Research, Permoserstr. 15, 04318 Leipzig, Germany
3Climate and Environmental Physics, University of Bern, Sidlerstrasse 5, 3012
Bern, Switzerland
4Oschger Centre for Climate Change Research, University of Bern, 3012 Bern,
Switzerland
5Royal Netherlands Meteorological Institute, De Bilt, The Netherlands
)

###### Abstract

Modelling dependencies between climate extremes is important for climate risk
assessment, for instance when allocating emergency management funds. In
statistics, multivariate extreme value theory is often used to model spatial
extremes. However, most commonly used approaches require strong assumptions
and are either too simplistic or over-parameterized. From a machine learning
perspective, Generative Adversarial Networks (GANs) are a powerful tool to
model dependencies in high-dimensional spaces. Yet in the standard setting,
GANs do not well represent dependencies in the extremes. Here we combine GANs
with extreme value theory (evtGAN) to model spatial dependencies in summer
maxima of temperature and winter maxima in precipitation over a large part of
western Europe. We use data from a stationary 2000-year climate model
simulation to validate the approach and explore its sensitivity to small
sample sizes. Our results show that evtGAN outperforms classical GANs and
standard statistical approaches to model spatial extremes. Already with about
50 years of data, which corresponds to commonly available climate records, we
obtain reasonably good performance. In general, dependencies between
temperature extremes are better captured than dependencies between
precipitation extremes due to the high spatial coherence in temperature
fields. Our approach can be applied to other climate variables and can be used
to emulate climate models when running very long simulations to determine
dependencies in the extremes is deemed infeasible.

##### Keywords

Extreme value theory; generative adversarial networks; spatial extremes;
climate model simulations

##  1 Introduction

Understanding and modelling climate extremes such as floods, heatwaves and
heavy precipitation, is of paramount importance because they often lead to
severe impacts on our socio-economic system (Field et al.,, 2012). Many
impacts are associated with compounding drivers, for instance multiple hazards
occurring at the same time at the same or different locations affecting the
same system (Leonard et al.,, 2014; Zscheischler et al.,, 2018). Ignoring
potential dependencies between multiple hazards can lead to severe
misspecification of the associated risk (Zscheischler and Seneviratne,, 2017;
Hillier and Dixon,, 2020). However, estimating dependencies between extremes
is a challenging task, requiring large amounts of data and/or suitable
approaches to model the phenomena of interest in a computationally feasible
time.

A particularly challenging class of events are spatially compounding events
(Zscheischler et al.,, 2020). Spatially compounding events occur when multiple
connected locations are affected by the same or different hazards within a
limited time window, thereby causing an impact. The compounding is established
via a system capable of spatial integration, which accumulates hazard impacts
in spatially distant locations. For instance, the global food systems are
vulnerable to multiple co-occurring droughts and heatwaves in key crop-
producing regions of the world (Anderson et al.,, 2019; Mehrabi and
Ramankutty,, 2019). Similarly, spatially extensive floods (Jongman et al.,,
2014) or sequences of cyclones (Raymond et al.,, 2020) can deplete emergency
response funds. Climate extremes can be correlated over very large distances
due to teleconnections (Boers et al.,, 2019) but modelling such dependencies
is challenging.

One approach to tackle the challenge of spatially correlated extremes is to
create a large amount of data by running very long simulations with state-of-
the-art climate models, which have the physical spatial dependencies built
into them. For instance, over the recent years for a number of climate models
large ensembles have been created (Deser et al.,, 2020). However, these
simulations typically have rather coarse resolution, are usually not
stationary in time and are very expensive to run.

Extreme value theory provides mathematically justified models for the tail
region of a multivariate distribution $(X_{1},\dots,X_{d})$, $d\geq 2$. This
enables the extrapolation beyond the range of the data and the accurate
estimation of the small probabilities of rare (multivariate) events.
Statistical methods building on this theory are popular in a broad range of
domains such as meteorology (Le et al.,, 2018), climate science (Naveau et
al.,, 2020) and finance (Poon et al.,, 2004). Applications have so far been
limited to low-dimensional settings for several reasons. On the one hand, even
for moderately large dimensions $d$, the fitting and simulation of parametric
models is computationally intensive because it requires computing complex
likelihoods (Dombry et al.,, 2017). On the other hand, the extremal dependence
structures in applications are difficult to model and the existing approaches
are either simplistic or over-parameterized. In spatial applications, for
instance, Euclidean distance is used to parametrize the isotropic dependence
of stationary max-stable random fields (Blanchet and Davison,, 2011). Most
real world datasets that cover larger spatial domains do, however, feature
non-stationarities in space that cannot be captured by the stationarity
assumed in current geostatistical models; see Engelke and Ivanovs, (2021) for
a review on recent extreme value methods in higher dimensions.

Methods from machine learning are well-suited for complex and non stationary
data sets. Their loss functions are, however, typically designed with the
purpose to predict well in the bulk of the distribution. It is therefore a
difficult problem to construct approaches with an accurate performance outside
the range of the training data. In prediction problems, one possibility is to
adapt the loss function to make the algorithm more sensitive to extreme
values, as for instance done in quantile regression forests (Meinshausen,,
2006) or extreme gradient boosting (Velthoen et al.,, 2021) for the prediction
of high conditional quantiles. Jalalzai et al., (2018) discuss this problem
for classification in extreme regions and propose a new notion of asymptotic
risk that helps to define classifiers with good generalization capacity beyond
the observations’ range in the predictor space.

Rather than predicting an unknown response, in this work we are interested in
a generative model that is able to learn a high-dimensional distribution
$(X_{1},\dots,X_{d})$ both in the bulk and in the extremes. We concentrate on
Generative Adversarial Networks (GANs), which are known to be competitive
models for multivariate density estimation. While classical applications of
GANs are often in the field of image analysis (Zhang et al.,, 2017; Karras et
al.,, 2017; Choi et al.,, 2018), they have been used much more broadly in
domains such as finance (Efimov et al.,, 2020), fraud detection (Zheng et
al.,, 2019), speech recognition (Sahu et al.,, 2019) and medicine (Schlegl et
al.,, 2017).

There are two different aspects that characterize the extremal properties of
multivariate data: the univariate tail behaviour of each margin and the
extremal dependence between the largest observations. For GANs, it has
recently been shown (Wiese et al.,, 2019; Huster et al.,, 2021) that the
marginal distributions of the generated samples are either bounded or light-
tailed if the input noise is uniform or Gaussian, respectively. Huster et al.,
(2021) propose to use heavy-tailed input to overcome this issue. Bhatia et
al., (2020) propose the ExGAN algorithm that uses conditional GANs to perform
importance sampling of extreme scenarios. The main difference to our approach
is that ExGAN simulates single extreme events, while our model will rely on
block maxima and therefore does not suffer from serial correlation and
seasonality issues. Moreover, ExGAN does not model the marginal distributions
with extreme value theory and, since the input noise of ExGAN is Gaussian,
this may be problematic especially for heavy-tailed data sets in view of the
results of Wiese et al., (2019) and Huster et al., (2021) mentioned above.

In this paper, we concentrate on modeling the extremal properties of both the
marginal distributions and the dependence structure in a realistic way even
for high-dimensional and spatially non-stationary data. Our model, called
evtGAN, combines the asymptotic theory of extremes with the flexibility of
GANs to overcome the limitations of classical statistical approaches to low
dimensional and stationary data. By using extrapolation for the marginals from
extreme value theory, our method can be applied to data sets with arbitrary
combinations of light-tailed and heavy-tailed distributions. We use a
stationary 2000-year climate simulation to test and validate our approach.

In sum, in this paper we demonstrate how our method combines the best of
extreme value theory and generative adversarial networks in order to
efficiently and accurately learn complex multidimensional distributions with
non-trivial extremal dependence structures. Furthermore, one can easily sample
from these learned distributions, which facilitates the study of complex data
such as spatially distributed climate extremes with spatially highly
heterogeneous dependence structure. The remainder of the article is structured
as follows: in Section 2 we briefly lay out the theoretical background of
extreme value theory and generative adversarial networks, and present the main
algorithm of our methodology. Section 3 contains a description of the data and
the model’s architecture and hyper-parameters used for its training, and a
presentation of the obtained results, which we then discuss in Section 4.
Finally, concluding remarks are laid out in Section 5.

##  2 Methodology

In the climate community there is a great need for methods able to efficiently
provide an empirical description of the climate system, including extreme
events, starting from as few ensemble runs as possible (Castruccio et al.,,
2019; Deser et al.,, 2020). Emulation techniques that tackle this challenge
have been proposed in the last few years, usually focusing on the spatial
correlation properties of the phenomena of interest (McKinnon et al.,, 2017;
McKinnon and Deser,, 2018; Link et al.,, 2019; Brunner et al.,, 2021). In this
work, our aim is to build an emulator specifically designed for extreme events
that is able to reproduce the spatial tail dependencies of the data and
extrapolate outside the range of the training samples. Our method does not
require to discard any simulations and is therefore highly efficient. In this
section we first recall some background on extreme value theory and generative
adversarial networks, and then describe our algorithm called evtGAN. Similar
to a copula approach (e.g. Brunner et al.,, 2021), it relies on the idea to
disentangle the modeling of the marginal distributions and the extremal
dependence structure in order to use extreme value theory for the former and
generative adversarial networks for the latter.

###  2.1 Extreme value theory

Extreme value theory provides mathematically justified methods to analyse the
tail of a $d$-dimensional random vector $\mathbf{X}=(X_{1},\dots,X_{d})$. This
branch of statistical theory has been widely applied in climate science, for
example to study droughts (Burke et al.,, 2010), floods (Asadi et al.,, 2018),
heatwaves (Tanarhte et al.,, 2015) and extreme rainfalls (Le et al.,, 2018).

For multivariate data, there are two different aspects that determine the
quality of an extrapolation method: the univariate tail behaviour of each
margin and the extremal dependence between the largest observations. We
explain how to measure and model those based on componentwise maxima. Let
$\mathbf{X}_{i}=(X_{i1},\dots,X_{id})$, $i=1,\dots,k$, be a sample of
$\mathbf{X}$ and let $\mathbf{M}_{k}=(M_{k1},\dots,M_{kd})$, where
$M_{kj}=\max(X_{1j},\dots,X_{kj})$ is the maximum over the $j$th margin.

####  2.1.1 Univariate theory

The Fisher–Tippett–Gnedenko theorem (e.g., Coles et al.,, 2001) describes the
limit behaviour of the univariate maximum $M_{kj}$. It states that if there
exist sequences $a_{kj}\in\mathbb{R}$, $b_{kj}>0$ such that

$\displaystyle\lim_{k\to\infty}\mathbb{P}\left(\frac{M_{kj}-a_{kj}}{b_{kj}}\leq
z\right)=G_{j}(z),\qquad n\rightarrow\infty,$ (1)

then if $G_{j}$ is the distribution function of a non-degenerate random
variable $Z_{j}$, it is in the class of generalized extreme value (GEV)
distributions with

$G_{j}(z)=\exp\left[-\left\\{1+\xi_{j}\left(\frac{z-\mu_{j}}{\sigma_{j}}\right)\right\\}_{+}^{-1/\xi_{j}}\right],\quad
z\in\mathbb{R},$

where $x_{+}$ denotes the positive part of a real number $x$, and
$\xi_{j}\in\mathbb{R}$, $\mu_{j}\in\mathbb{R}$ and $\sigma_{j}>0$ are the
shape, location and scale parameters, respectively. The shape parameter is the
most important parameter since it indicates whether the $j$th margin is heavy-
tailed ($\xi_{j}>0$), light-tailed ($\xi=0$) or whether it has a finite uper
end-point ($\xi<0$).

It can be shown that under mild conditions on the distribution of margin
$X_{j}$, appropriate sequences exist for (1) to hold. The above theorem
therefore suggests to fit a generalized extreme value distribution
$\widehat{G}_{j}$ to maxima taken over blocks of the same lengths, which is
common practice when modelling yearly maxima. The shape, location and scale
parameters can then be estimated in different ways, including moment-based
estimators Hosking, (1985), Bayesian methods (Yoon et al.,, 2010) and maximum
likelihood estimation (Hosking et al.,, 1985). We use the latter in this work.
This allows extrapolation in the direction of the $j$th marginal since the
size of a $T$-year event can be approximated by the $(1-1/T)$-quantile of the
distribution $\widehat{G}_{j}$, even if $T$ is larger than the length of the
data record.

####  2.1.2 Multivariate theory

For multivariate data, the correct extrapolation of the tail not only depends
on correct models for the marginal extremes, but also on whether the
dependence between marginally large values is well captured. For two
components $X_{i}$ and $X_{j}$ with limits $Z_{i}$ and $Z_{j}$ of the
corresponding maxima in (1), this dependence is summarized by the extremal
correlation $\chi_{ij}\in[0,1]$ (Schlather and Tawn,, 2003)

$\chi_{ij}=\lim_{q\to 1}\mathbb{P}(F_{i}(X_{i})>q\mid
F_{j}(X_{j})>q)\in[0,1].$

The larger $\chi_{ij}$, the stronger the dependence in the extremes between
the two components, that is, between $Z_{i}$ and $Z_{j}$. If $\chi_{ij}>0$ we
speak of asymptotic dependence, and otherwise of asymptotic independence.
While the extremal correlation is a useful summary statistic, it does not
reflect the complete extremal dependence structure between $X_{i}$ and
$X_{j}$. Under asymptotic dependence, the latter can be characterized by the
so-called spectral distribution (de Haan and Resnick,, 1977). Let
$\tilde{X}_{j}=-1/\log F_{j}(X_{j})$, $j=1,\dots,d$, be the data normalized to
standard Fréchet margins. The spectral distribution describes the extremal
angle of $(X_{i},X_{j})$, given that the radius exceeds a high threshold, that
is,

$H(x)=\lim_{u\to\infty}\mathbb{P}(\tilde{X}_{i}/(\tilde{X}_{i}+\tilde{X}_{j})\leq
x\mid\tilde{X}_{i}+\tilde{X}_{j}>u),\quad x\in[0,1].$

Under strong extremal dependence, the spectral distribution centers around
$1/2$; under weak extremal independence, it has mass close to the boundary
points $0$ and $1$.

In multivariate extreme value theory, a popular way to ensure a correct
extrapolation of tail dependence under asymptotic dependence is by modeling
the joint distribution $\mathbf{Z}=(Z_{1},\dots,Z_{d})$ of the limits in (1)
as a max-stable distribution with multivariate distribution function

$\mathbb{P}(G_{1}(Z_{1})\leq z_{1},\dots,G_{d}(Z_{d})\leq
z_{d})=\exp\\{-\ell(z_{1},\dots,z_{d})\\},\quad z_{1},\dots,z_{d}\in[0,1],$
(2)

where $\ell$ is the so-called stable tail dependence function (e.g., De Haan
and Ferreira,, 2007). This is a copula approach in the sense that the right-
hand side is independent of the original marginal distributions $G_{j}$ and
only describes the dependence structure. In practice, the $G_{j}$ are replaced
by estimates $\hat{G}_{j}$ that are fitted to the data first. The right-hand
side of (2) is also called an extreme value copula (e.g., Gudendorf and
Segers,, 2010).

A natural extension of max-stable distributions to the spatial setting is
given by max-stable processes $\\{Z(t):t\in\mathbb{R}^{d}\\}$ (e.g., Davison
et al.,, 2012), where the maxima $Z_{i}=Z(t_{i})$ are observed at spatial
locations $t_{i}$, $i=1,\dots,d$. This type of models have been widely applied
in many different areas such as flood risk (Asadi et al.,, 2015), heat waves
(Engelke et al.,, 2019) and extreme precipitations (Buishand et al.,, 2008). A
popular parametric model for a max-stable process $Z$ is the Brown–Resnick
process (Kabluchko et al.,, 2009; Engelke et al.,, 2015). This model class is
parameterized by a variogram function $\gamma$ on $\mathbb{R}^{d}$. In this
case, each bivariate distribution is in the so-called Hüsler–Reiss family
(Hüsler and Reiss,, 1989), and the corresponding dependence parameter is the
variogram function evaluated at the distance between the two locations,

$\lambda_{ij}=\gamma(\|t_{i}-t_{j}\|),\qquad i,j=1,\dots,d.$

One can show that the extremal correlation coefficient for this model is given
by $\chi_{ij}=2-2\Phi(\sqrt{\lambda_{ij}}/2)$, where $\Phi$ is the standard
normal distribution function. A popular parametric family is the class of
fractal variograms that can be written as $\gamma_{\alpha,s}(h)=h^{\alpha}/s$,
$\alpha\in(0,2]$, $s>0$. While these methods enjoy many interesting
properties, they are often quite complex to fit (Dombry et al.,, 2017) and
impose strong assumptions on the spatial phenomena of interest, such as
spatial stationarity, isotropic behaviour and asymptotic dependence between
all pairs of locations.

A simple way of fitting a Brown–Resnick process to data is to compute the
empirical versions $\hat{\chi}_{ij}$ of the extremal correlation coefficients
and then numerically find the parameter values (e.g., $\alpha$ and $s$ for the
fractal variogram family) whose implied model extremal coefficients minimise
the squared error compared to the empirical ones.

###  2.2 Generative adversarial network

Generative adversarial networks (Goodfellow et al.,, 2014) are attracting
great interest thanks to their ability of learning complex multivariate
distributions that are possibly supported on lower-dimensional manifolds.
These methods are generative since they allow the generation of new realistic
samples from the learnt distribution that can be very different from the
training samples.

Given observations $\mathbf{U}_{1},\dots,\mathbf{U}_{n}$ from some
$d$-dimensional unknown target distribution $p_{data}$, generative adversarial
networks can be seen as a game opposing two agents, a discriminator $D$ and a
generator $G$. The generator $G$ takes as input a random vector
$\mathbf{Y}=(Y_{1},\dots,Y_{p})$ from a $p$-dimensional latent space and
transforms it to a new sample $\mathbf{U}^{*}=G(\mathbf{Y})$. The components
$Y_{j}$, $j=1,\dots,p$, are independent and typically have a standard uniform
or Gaussian distribution.

The discriminator $D$ has to decide whether a new sample $\mathbf{U}^{*}$
generated by $G$ is fake and following the distribution $p_{G}$, or coming
from the real observations with distribution $p_{data}$. The discriminator
expresses its guesses with a value between 0 and 1 corresponding to its
predicted probability of the sample coming from the real data distribution
$p_{data}$. Both the generator and the discriminator become better during the
game, in which $D$ is trained to maximise the probability of correctly
labelling samples from both sources, and $G$ is trained to minimise $D$’s
performance and thus learns to generate more realistic samples.

Mathematically, the optimization problem is a two-player minimax game with
cross-entropy objective function

$\min\limits_{G}\max\limits_{D}\mathbb{E}_{\mathbf{U}\sim
p_{data}}[\log(D(\mathbf{U}))]+\mathbb{E}_{\mathbf{Y}\sim
p_{\mathbf{Y}}}[\log(1-D(G(\mathbf{Y}))],$ (3)

where $\mathbf{Y}$ is a random vector sampled from a latent space. In
equilibrium, the optimal generator satisfies $p_{G}=p_{data}$ (Goodfellow et
al.,, 2014).

In practice, the discriminator and the generator are modeled through feed-
forward neural networks. While the equilibrium with $p_{G}=p_{data}$ is only
guaranteed thanks to the convex-concave property in the non-parametric
formulation (3), the results remain very good as long as suitable adjustments
to the losses, training algorithm and overall architecture are made to improve
on the stability and convergence of training. For instance, a standard
adjustment to the generator’s loss function was already proposed in the GANs
original paper (Goodfellow et al.,, 2014), and consists on training the
generator to minimise
$-\log\left(D\left(G\left(\mathbf{Y}\right)\right)\right)$ rather than
$\log\left(1-D\left(G\left(\mathbf{Y}\right)\right)\right)$. This new loss
function is called the non-saturating heuristic loss, and was suggested to
mitigate close to null gradients that occur in early training. Indeed, as the
generator is still poorly trained, the discriminator can easily detect its
samples and therefore associates values close to zero to them. This results in
a very slow improvement of the generator during backpropagation as the
gradients associated with values around
$D\left(G\left(\mathbf{Y}\right)\right)=0$ are too small, and thus its loss
$L\left(G\right)=\log\left(1-D\left(G\left(\mathbf{Y}\right)\right)\right)$
saturates. Training the generator to minimise
$-\log\left(D\left(G\left(\mathbf{Y}\right)\right)\right)$ rather than
$\log\left(1-D\left(G\left(\mathbf{Y}\right)\right)\right)$ provides stronger
gradients in early training without compromising on the dynamics of the two
agents.

The architecture of the generator $G$ and the discriminator $D$ can be
specifically designed for the structure of the data under consideration. In
this direction, Radford et al., (2015) introduces the use of the convolutional
neural network as an add-on feature to the GAN yielding a deep convolutional
generative adversarial network model (DCGAN). It is considered the preferred
standard for dealing with image data and more generally any high dimensional
vector of observations describing a complex underlying dependence structure,
showing superior performance in the representation and recognition fields.

###  2.3 evtGAN

Here we give an overview of evtGAN. Let $\mathbf{X}_{1},\dots,\mathbf{X}_{m}$
be independent observations of dimensionality $d$, where
$\mathbf{X}_{i}=\left(X_{i1},\ldots,X_{id}\right)$ and $m=k\cdot n$. Now let
$\mathbf{Z}_{1},\dots,\mathbf{Z}_{n}$ be independent observations obtained as
block maxima taken sequentially over the span of $k$ observations, i.e.,
$\mathbf{Z}_{1}=\max(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$,
$\mathbf{Z}_{2}=\max(\mathbf{X}_{k+1},\dots,\mathbf{X}_{2k})$ and so on (see
also Section 2.1), where the maxima are taken in each component. The evtGAN
algorithm (Algorithm 1) below takes these empirical block maxima as input and
follows a copula approach where marginal distributions and the dependence
structure are treated separately. In particular, this allows us to impose
different amounts of parametric assumptions on the margins and the dependence.

It is known that classical GANs that are trained with bounded or light tailed
noise input distributions in the latent space will also generate bounded or
light-tailed samples, respectively (see Wiese et al.,, 2019; Huster et al.,,
2021). That means that they are not well-suited for extrapolating in a
realistic way outside the range of the data. For this reason, we suppose that
the approximation of the margins by generalized extreme value distributions
$G_{j}$ as in (1) holds, which allows for extrapolation beyond the data range.
On the other hand, we do not make explicit use of the assumption of a
multivariate max-stable distribution as in (2). This has two reasons. First,
while the univariate limit (1) holds under very weak assumptions, the
multivariate max-stability as in (2) requires much stronger assumptions (see
Resnick,, 2008) and may be too restrictive for cases with asymptotic
independence (e.g., Wadsworth and Tawn,, 2012). Second, even if the data
follow a multivariate max-stable distribution (2), the probabilistic structure
is difficult to impose on a flexible machine learning model such as a GAN.

Input: Observations $\mathbf{Z}_{i}=(Z_{i1},\dots,Z_{id})$, $i=1,\dots,n$.

1: For $j=1,\dots,d$, fit a GEV distribution $\widehat{G}_{j}$ to the data
$Z_{1j},\dots,Z_{nj}$ with estimated parameters
$(\hat{\mu}_{j},\hat{\sigma}_{j},\hat{\xi}_{j})$.

2: Normalize all margins empirically to a standard uniform distribution to
obtain pseudo observations
$\mathbf{U}_{i}=(\widehat{F}_{1}(Z_{i1}),\dots,\widehat{F}_{d}(Z_{id})),\quad
i=1,\dots,n,$ where $\widehat{F}_{j}$ is the empirical distribution function
of the $Z_{1j},\dots,Z_{nj}$.

3: Train a DCGAN $G$ on the normalized data
$\mathbf{U}_{1},\dots,\mathbf{U}_{n}$ based on the loss in equation (3).

4: Generate $n^{*}$ new data points
$\mathbf{U}_{1}^{*},\dots,\mathbf{U}_{n}^{*}$ from $G$ with uniform margins,
$i=1,\dots,n^{*}$.

5: Normalize back to the scale of the original observations
$\mathbf{Z}_{i}^{*}=(\widehat{G}_{1}^{-1}(U_{i1}^{*}),\dots,\widehat{G}_{d}^{-1}(U_{id}^{*})),\quad
i=1,\dots,n^{*}.$

Output: Set of new generated observations
$\mathbf{Z}^{*}_{i}=(Z_{i1},\dots,Z^{*}_{id})$, $i=1,\dots,n^{*}$.

Algorithm 1 evtGAN

The margins are first fitted by a parametric GEV distribution (line 1 in
Algorithm 1). They are then normalized to (approximate) pseudo-observations by
applying the empirical distributions functions to each margin (line 2). This
standardization to uniform margins stabilizes the fitting of the GAN.
Alternatively, we could use the fitted GEV distributions for normalization but
this seems to give slightly worse results in practice. The pseudo-observations
contain the extremal dependence structure of the original observations. Since
this dependence structure is very complex in general and our focus is to
reproduce gridded spatial fields, we do not rely on restrictive parametric
assumptions but rather learn it by a highly flexible DCGAN (line 3). From the
fitted model we can efficiently simulate any number of new pseudo-observations
that have the correct extremal dependence (line 4). Finally, we transform the
generated pseudo-observations back to the original scale to have realistic
properties of the new samples also in terms of the margins (line 5).

##  3 Application

We assess the ability of the proposed method to correctly model the spatial
tail dependence structure between climate extremes. Because observational
records are relatively short and contain temporal trends, here we rely on 2000
years of climate model output that is representative of present-day weather
over western Europe. In particular, we apply our approach to summer
temperature and winter precipitation maxima, yielding 2000 and 1999 maxima for
temperature and precipitations, respectively. We then compare the performance
of our algorithm to the Brown–Resnick model, a state-of-the-art statistical
model for spatial extremes. We use the fractal variogram family and the
fitting described in Section 2.1.2. In order to retrieve a non-biased estimate
of the error that each of the method makes, we divide the data into a training
set of 50 observations, and the rest is taken as a test set. Both methods are
fitted using only the training data set and the evaluation is made considering
the information contained in the test set as ground truth. When not stated
otherwise, the results from the Brown–Resnick model are analytical, while the
results from evtGAN are obtained simulating $10^{\prime}000$ data points.

###  3.1 Data

The model experiment, which uses large ensemble simulations with the EC-Earth
global climate model (v2.3, Hazeleger et al.,, 2012) has been widely used in
climate impact studies (Van der Wiel et al., 2019a, ; Van Kempen et al.,,
2021; Tschumi et al.,, 2021; Vogel et al.,, 2021; Van der Wiel et al.,, 2021)
and was originally designed to investigate the influence of natural
variability on climate extremes. A detailed description of these climate
simulations is provided in Van der Wiel et al., 2019b , here we only provide a
short overview. The large ensemble contains 2000 years of daily weather data,
representative of the present-day climate. Present-day was defined by the
observed value of global mean surface temperature (GMST) over the period
2011-2015 (HadCRUT4 data, Morice et al.,, 2012); the five-year time slice in
long transient model simulations (forced by historical and Representative
Concentration Pathway (RCP) 8.5 scenarios) with the same GMST was simulated
repeatedly. To create the large ensemble, at the start of the time slice,
twenty-five ensemble members were branched off from the sixteen transient
runs. Each ensemble member was integrated for the five-year time slice period.
Differences between ensemble members were forced by choosing different seeds
in the atmospheric stochastic perturbations (Buizza et al.,, 1999). This
resulted in a total of $16\times 25\times 5=2000$ years of meteorological
data, at T159 horizontal resolution (approximately 1∘), among which we
selected temperature (Kelvin) and precipitations (meter/day) for this paper.

We choose an area such that it is big enough to be relevant for climate
application while being small enough to ensure a fair comparison of our
approach and the classical approach in statistics to model spatial tail
dependence, the Brown–Resnick process. Our analysis thus focuses on a total of
$18\times 22$ grid points covering most of western Europe. For that area we
compute for each grid point the summer temperature maxima and winter
precipitation maxima.

###  3.2 Architecture and hyper-parameters

In our application, we consider the non-saturating heuristic loss (Goodfellow
et al.,, 2014) for the generator, and the standard empirical cross-entropy
loss for the discriminator. Training is done iteratively, where the
discriminator is allowed to train two times longer than the generator. A
bigger number of discriminator training iterations per generator training
iteration was initially considered but it did not achieve better results for
the considered sample sizes and input dimensions but only resulted in longer
training. We incorporate an exponential moving average scheme (Gidel et al.,,
2018) to the training algorithm as it has demonstrated more stability, a far
better convergence and improved results.

Since we focus on reproducing gridded spatial fields, we make use of the DCGAN
model (Radford et al.,, 2015) as an architectural design to take advantage of
convolution (Fukushima,, 1980), and for a more accurate analysis of the input
data, we apply zero-padding of a single layer to the train set observations
before they are fed to the discriminator (albawi2017understanding), increasing
the dimension of its input to $20\times 24$. We use classical regularisation
tools such as drop-out (Srivastava et al.,, 2014), batch-normalisation (Ioffe
and Szegedy,, 2015) and the Adam optimizer (Kingma and Ba,, 2014) with a
learning rate of $2\times 10^{-4}$ and batch size of $50$ to train the neural
network, and a total of $30,000$ training epochs. The value of the exponential
moving average parameter in the training algorithm is set to $\alpha=0.9$. The
values of the hyper-parameters are the result of an extensive heuristic
tuning, and are as such only valid for the dimension of the input and the
considered sample sizes involved in training the network. That being said, we
are fairly confident that a minimal intervention in tuning the hyper-
parameters/architecture would be required if for instance higher resolution
inputs were to be used, as these models are fairly robust to such changes.
Indeed, models based on the general GANs paradigm are traditionally trained on
much higher dimensional objects such as images, and they learn well using an
architecture that is very close to the one we have adopted. One may argue that
the underlying distributions defining real world images are "simpler", and
therefore easier to learn, which is why we propose the evtGAN model to address
the increased complexity of the underlying distribution of the data at hand:
initially, we started with a coarser grid (lower resolution input), as we
considered France only, then decided to cover a larger part of Europe as the
results were very good. In this change, the model’s architecture needed to be
adjusted only by adding a few layers to the generator. We therefore remain
confident that a higher resolution grid would only entail a rather small
hyper-parameter/architectural adjustment for the results to remain
satisfactory. For a given input resolution, increasing the complexity of the
model (using a deeper one) just increased training time and pushed convergence
further down the line. Figures 1 and 2 illustrate the final architectures of
the discriminator and the generator, as well as the values of some of the
hyper-parameters used for training.

Figure 1: Discriminator’s architecture and values of hyper-parameters used for
training. lrelu stands for the leaky relu function and its argument
corresponds to the value of the slope when $x<0$; drop stands for drop-out
(see Section 3.2) and its argument corresponds to the drop-out probability; F
stands for the kernel’s dimension, and S for the value of the stride; FC
stands for fully connected, and sigmoid for the sigmoid activation function
Figure 2: Generator’s architecture and values of hyper-parameters used for
training. lrelu stands for the leaky relu function and its argument
corresponds to the value of the slope when $x<0$; drop stands for drop-out
(see Section 3.2) and its argument corresponds to the drop-out probability; F
stands for the kernel’s dimension, and S for the value of the stride; BN
stands for batch-normalisation (see Section 3.2); FC stands for fully
connected

###  3.3 Results

A particularity of evtGAN is that it decomposes the full high-dimensional
distribution of the data into its marginal distributions and its dependence
structure and processes them separately. We thus first report estimated key
parameters of the marginal distributions. For temperature, the location
parameter $\mu$ is typically higher over land than over oceans (Fig. 3a).
Furthermore, there is a trend towards lower values of $\mu$ for more northern
regions, illustrating the well-known latitudinal temperature gradient from
south to north. The scale parameter $\sigma$, representing the width of the
distribution, is also higher over land than over the ocean, with a fairly
homogeneous distribution over land (Fig. 3b). The shape parameter $\xi$ is
well below 0 for most regions except some areas in the Atlantic north of Spain
and in the Mediterranean at the north African coast, indicating a bounded tail
(Fig. 3c).

Figure 3: The generalized extreme value distribution parameters estimated for
each grid point for temperature (a-c) and precipitation (d-f) extremes. On the
left (a, d) the mean parameter $\mu$, in the middle (b, e) the scale parameter
$\sigma$ and on the right (c, f) the shape parameter $\xi$

For precipitation, the location parameter $\mu$ is more similar between land
and ocean compared to temperature but is a little bit more heterogeneous in
space (Fig. 3d), illustrating orographic effects on precipitation extremes.
For the scale parameter $\sigma$ there is also no clear land-sea contrast but
areas with relatively much higher values (corresponding to higher interannual
variability in extremes) occur over the Western part of the Iberian peninsula
and along a coastal band at the northern Mediterranean coast (Fig. 3e). The
shape parameter $\xi$ is spatially quite heterogeneous, with values ranging
from -0.15 to 0.25 (Fig. 3f). Overall, the parameter of the extreme value
distributions of precipitation extremes show much higher spatial heterogeneity
than the one for temperature extremes, suggesting that it might be more
difficult to learn a model that represents well all tail dependencies for the
entire region.

Next we look at examples of bivariate scatterplots of simulated temperature
extremes from the different approaches for three cases with varying tail
dependence (Fig. 4). The three rows correspond to three pairs of grid points
with weak, mild and strong tail dependence, respectively. The train sets (Fig.
4a, f, k) illustrate the distributions from which the models are tasked to
learn while the test sets (Fig. 4b, g, l) are used as ground truth. As can be
clearly seen, extremes are much less likely to co-occur between locations that
are characterized by weak tail dependence (Fig. 4b) compared to locations that
are characterized by strong tail dependence (Fig. 4l). Purely based on visual
judgement, it seems that evtGAN is able to characterize the different tail
dependencies relatively well and can simulate samples outside of the range of
the train set (Fig. 4c, h, m) whereas DCGAN only simulates distributions
bounded to the range of the train set (Fig. 4d, i, n) and Brown–Resnick tends
to overestimate tail dependence in cases where tail dependence is weak (Fig.
4e, j).

Figure 4: Bivariate plots of temperature extremes based on train and test sets
and simulated through the different presented methods. Shown are selected
pairs of locations with varying tail dependence. Columns from left to right:
train, test, evtGAN, DCGAN and Brown–Resnick. From top to bottom: weak tail
dependence (a-e), mild tail dependence (f-j), strong tail dependence (k-o).
The colours represent the empirically estimated density: the spectrum goes
from yellow to grey, reflecting a decrease in density

The corresponding figure for precipitation extremes is shown in Fig. 5.
Conclusions are similar as for temperature extremes. The evtGAN simulates
distributions with different tail dependencies and samples that are far
outside the range of the train set (Fig. 5c, h, m). DCGAN simulates
distributions that are bounded to the train set range (Fig. 5d, i, n). On the
other hand, Brown–Resnick overestimates tail dependence in the case of mild
tail dependence (Fig. 5j), and underestimates tail dependence when tail
dependence is strong (Fig. 5o).

Figure 5: Bivariate plots of precipitation extremes based on train and test
sets and simulated through the different presented methods. Shown are selected
pairs of locations with varying tail dependence. Columns from left to right:
train, test, evtGAN, DCGAN and Brown–Resnick. From top to bottom: weak tail
dependence (a-e), mild tail dependence (f-j), strong tail dependence (k-o).
The colours represent the empirically estimated density: the spectrum goes
from yellow to grey, reflecting a decrease in density

A scatterplot of bivariate extremal correlations between 100 randomly selected
locations estimated from the train set, evtGAN and Brown–Resnick,
respectively, against estimates based on the test set (1950 samples) is shown
in Fig. 6. The estimates derived directly from the train sets (Fig. 6a, d) are
the benchmark, and by design better performance is not possible. Clearly,
pairs of locations with stronger tail dependence are much more likely for
temperature (Fig. 6a) than for precipitation (Fig. 6d), confirming the
impression obtained from the spatial homogeneity of the parameters in the
extreme value distributions (Fig. 3). Furthermore, evtGAN seems to better
capture the observed relationship. Brown–Resnick has difficulties in
particular with pairs that have weak or no tail dependence (extremal
correlation equals 0, lower left in the figures), which is consistent with
Fig. 4 and Fig. 5.

Figure 6: Scatter plots of the extremal correlations for temperature extremes
(a-c) and precipitation extremes (d-f) between $100$ randomly selected
locations. The x-axes always show the estimates based on the test set.
Estimates on the y-axes are based on the train set (a, d), on the output of
evtGAN (b, d) and on the Brown–Resnick model

To further explore the behaviour of evtGAN and Brown–Resnick in dealing with
different levels of dependence, Fig. 7 shows a comparison of the estimated
spectral distribution for the same three pairs of chosen grid points
characterized by weak (a, d), mild (b, e) and strong tail dependence (c, f),
respectively. The magenta bars show the estimates based on the 1950 samples in
the test set. The estimates of evtGAN are very close to the ground truth for
all cases, i.e., weak, mild and strong dependence (red lines in Fig. 7),
except for some bias in the case of mild dependence in precipitation extremes
(Fig. 7e). In contrast, the performance of the Brown–Resnick model (black
lines in Fig. 7) is much more variable. It captures relatively well the two
pairs with weak tail dependence (Fig. 7a, d) and does a decent job for strong
dependence in temperature extremes (Fig. 7c) but fails completely for the
remaining three cases (Fig. 7b,e, f). Furthermore, in contrast to
Brown–Resnick, evtGAN is able to represent asymmetric dependence structures,
as is evident from Fig. 7.

Figure 7: Spectral distributions for a threshold of $0.95$ for selected pairs
of locations with varying tail dependence for temperature (a-c) and
precipitation (d-f). (a, d) weak tail dependence, (b, e) mild tail dependence,
(c, f) strong tail dependence. In red the kernel density estimate of the
evtGAN, in blue the Brown–Resnick model and in magenta bars the ground truth

We finally present a sensitivity analysis of the performance of evtGAN for
different realistic sample sizes, namely 30, 50 and 100, while learning for
30,000 epochs (Fig. 8). Overall the error tends to decrease the longer the
models learn and it does not seem to matter whether we evaluate the
performance on the train set or the test set (difference between black and red
lines in Fig. 8 is small). For small sample sizes ($n=30$) the best
performance might be reached for epochs smaller than 30,000 (Fig. 8a, b) and
the $l^{2}$ norm between the extremal coefficients of evtGAN and the train set
could be used as a stopping criterion. An improvement in performance with
increasing sample size is clearly visible (compare Fig. 8a, b with Fig. 8e,
f). Furthermore, the model errors in temperature extremes are smaller than the
ones for precipitation extremes (left versus right column in Fig. 8).

Figure 8: Model performance versus training epochs for different sample sizes
in evtGAN for temperature (a, c, e) and precipitation (b, d, f) extremes. A
number of observations equal to $n_{test}=N-n_{train}$ was sampled from
evtGAN, where $N=2000$ for temperature, and $N=1999$ for precipitation. The
mean $l^{2}$ norms for train (black) and test set (red) are defined as
$C_{te}=||\chi_{evtG}-\chi_{te}||_{2}$,
$C_{tr}=||\chi_{evtG}-\chi_{tr}||_{2}$, where $\chi_{evtG}$, $\chi_{te}$ and
$\chi_{tr}$ denote the vectors of extremal correlations calculated on the
samples from evtGAN, the test set and the train set, respectively

##  4 Discussion

The use of evtGAN combines the best of two worlds: correct extrapolation based
on extreme value theory on the one hand, and flexible dependence modeling
through GANs on the other hand. GANs are an excellent tool to model complex
dependencies in high-dimensional spaces. However, they are typically not
tailored to model extremes in the marginals well. Indeed, for a standard DCGAN
implementation where the marginals are not estimated by GEV distributions but
empirically transformed, Figs. 4d, i, n and Figs. 5d, i, n show that the
generated samples are bounded by the range of the training sample. For an
accurate extrapolation that resembles the marginal distributions of the test
set, extreme value theory is required (panels c, e, h, j, m and o of Figs. 4
and 5).

On the other hand, classical methods of spatial extreme value theory such as
the Brown–Resnick process have accurate extrapolation properties for the
marginal distributions. However, for an application to a spatially
heterogeneous data set on a large domain (Fig. 3) their parametric assumption
may be too restrictive. Indeed, Figs. 6c and f show a clear bias of the
Brown–Resnick model in terms of bivariate extremal correlations, which is
particularly visible for pairs with weak extremal dependence. Another
indication for this bias can be seen in Figs. 5j and o where the distributions
of the Brown–Resnick samples differ strongly from the test set distributions.
The evtGAN does not make prior assumptions on spatial stationarity or isotropy
and therefore it does not exhibit a bias (points in Figs. 6b and e are
centered around the diagonal). This is particularly noteworthy since modeling
complex non-stationarities for extremes is very difficult with parametric
models (Huser and Genton,, 2016; Engelke and Hitz,, 2020). Considering the
fitted bivariate distributions of evtGAN and Brown–Resnick underlines this
point.

The spectral distributions of the Brown–Resnick model are restricted to a
parametric class of functions, which, for instance, are symmetric about $1/2$.
The blue lines in Fig. 7 show that this is too restrictive for our data since
the strength of dependence is not correctly modelled (Fig. 7b) or the
asymmetry is not captured (Fig. 7c,f). The evtGAN (red lines) on the other
hand can model weak and strong dependences, and it even adapts to possible
asymmetries. This is also evident from the scatterplots in Fig. 4 and Fig. 5,
where the shape of the Brown–Resnick samples are restricted to a parametric
sub-class of distributions, the so-called Hüsler–Reiss family (Hüsler and
Reiss,, 1989).

A further restriction of classical spatial max-stable processes is the fact
that they are limited to modeling asymptotic dependence. For precipitation, it
can be seen in Fig. 6d-f that most of the test extremal correlations are close
to zero, indicating asymptotic independence. While the evtGAN is able to
capture this fairly well (Fig. 6e), the Brown–Ressnick model always has
positive extremal correlations, explaining the bias in the bottom left corner
of Fig. 6f. A spatial asymptotically independent model (e.g., Wadsworth and
Tawn,, 2012) would be a possible remedy for this, but such processes would
still suffer from the limitation induced by non-stationarity and asymmetry
described above.

Overall, evtGAN tends to perform better in capturing dependencies between
temperature extremes than precipitation extremes (Fig. 8). This is likely
related to the fact that extremes in temperature are more spatially coherent
(Keggenhoff et al.,, 2014) (Fig. 3).

##  5 Conclusions

Understanding and evaluating the risk associated with extreme events is of
primal importance for society, as recently emphasized in the 6th Assessment
Report of the Intergovernmental Panel on Climate Change (Seneviratne et al.,,
2021). Extreme event analysis and impact assessments are often limited by the
available sample sizes. Furthermore, simulations with complex climate models
are very expensive. Here we combine a machine learning approach with extreme
value theory to model complex spatial dependencies between extreme events in
temperature and precipitation across Europe based on a limited sample size. We
demonstrate that this hybrid approach outperforms the typically used approach
in multivariate extreme value theory and can well represent the marginal
distributions and extremal dependencies across spatially distributed climate
extremes. The approach can be easily adapted to other types of extremes and
used to create large sample sizes that are often required for climate risk
assessments.

##### Funding Statement

This work was funded by the Swiss National Science Foundation (grant nos.
179876, 189908 and 186858) and the Helmholtz Initiative and Networking Fund
(Young Investigator Group COMPOUNDX; grant agreement no. VH-NG-1537).

##### Competing Interests

None

##### Data Availability Statement

The temperature and precipitation maxima alongside the scripts to reproduce
the results in this study are available at
https://doi.org/10.5281/zenodo.5821485.

##### Author Contributions

Conceptualization: Y.B.; E.V.; S.E. Methodology: Y.B; E.V.; S.E. Formal
analysis: Y.B. Software: Y.B.; E.V. Data curation: K.v.d.W. Data
visualisation: Y.B. Supervision: J.Z.; S.E. Writing original draft: J.Z.;
E.V.; S.E. All authors approved the final submitted draft.

## References

  * Anderson et al., (2019) Anderson, W. B., Seager, R., Baethgen, W., Cane, M., and You, L. (2019).  Synchronous crop failures and climate-forced production variability.  Science Advances, 5(7):eaaw1976.
  * Asadi et al., (2015) Asadi, P., Davison, A. C., Engelke, S., et al. (2015).  Extremes on river networks.  The Annals of Applied Statistics, 9(4):2023–2050.
  * Asadi et al., (2018) Asadi, P., Engelke, S., and Davison, A. C. (2018).  Optimal regionalization of extreme value distributions for flood estimation.  Journal of Hydrology, 556:182–193.
  * Bhatia et al., (2020) Bhatia, S., Jain, A., and Hooi, B. (2020).  Exgan: Adversarial generation of extreme samples.  arXiv preprint arXiv:2009.08454.
  * Blanchet and Davison, (2011) Blanchet, J. and Davison, A. C. (2011).  Spatial modeling of extreme snow depth.  Ann. Appl. Stat., 5:1699–1725.
  * Boers et al., (2019) Boers, N., Goswami, B., Rheinwalt, A., Bookhagen, B., Hoskins, B., and Kurths, J. (2019).  Complex networks reveal global pattern of extreme-rainfall teleconnections.  Nature, 566(7744):373–377.
  * Brunner et al., (2021) Brunner, M. I., Gilleland, E., and Wood, A. W. (2021).  Space–time dependence of compound hot–dry events in the United States: assessment using a multi-site multi-variable weather generator.  Earth System Dynamics, 12(2):621–634.
  * Buishand et al., (2008) Buishand, T., De Haan, L., Zhou, C., et al. (2008).  On spatial extremes: with application to a rainfall problem.  The Annals of Applied Statistics, 2(2):624–642.
  * Buizza et al., (1999) Buizza, R., Milleer, M., and Palmer, T. N. (1999).  Stochastic representation of model uncertainties in the ECMWF ensemble prediction system.  Quarterly Journal of the Royal Meteorological Society, 125(560):2887–2908.
  * Burke et al., (2010) Burke, E. J., Perry, R. H., and Brown, S. J. (2010).  An extreme value analysis of uk drought and projections of change in the future.  Journal of Hydrology, 388(1-2):131–143.
  * Castruccio et al., (2019) Castruccio, S., Hu, Z., Sanderson, B., Karspeck, A., and Hammerling, D. (2019).  Reproducing internal variability with few ensemble runs.  Journal of Climate, 32(24):8511–8522.
  * Choi et al., (2018) Choi, Y., Choi, M., Kim, M., Ha, J.-W., Kim, S., and Choo, J. (2018).  Stargan: Unified generative adversarial networks for multi-domain image-to-image translation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 8789–8797.
  * Coles et al., (2001) Coles, S., Bawa, J., Trenner, L., and Dorazio, P. (2001).  An introduction to statistical modeling of extreme values, volume 208.  Springer.
  * Davison et al., (2012) Davison, A. C., Padoan, S. A., and Ribatet, M. (2012).  Statistical modeling of spatial extremes.  Statistical Science, pages 161–186.
  * De Haan and Ferreira, (2007) De Haan, L. and Ferreira, A. (2007).  Extreme value theory: an introduction.  Springer Science & Business Media.
  * de Haan and Resnick, (1977) de Haan, L. and Resnick, S. I. (1977).  Limit theory for multivariate sample extremes.  Z. Wahrscheinlichkeitstheorie und Verw. Gebiete, 40(4):317–337.
  * Deser et al., (2020) Deser, C., Lehner, F., Rodgers, K., Ault, T., Delworth, T., DiNezio, P., Fiore, A., Frankignoul, C., Fyfe, J., Horton, D., et al. (2020).  Insights from Earth system model initial-condition large ensembles and future prospects.  Nature Climate Change, pages 1–10.
  * Dombry et al., (2017) Dombry, C., Engelke, S., and Oesting, M. (2017).  Bayesian inference for multivariate extreme value distributions.  Electron. J. Stat., 11:4813–4844.
  * Efimov et al., (2020) Efimov, D., Xu, D., Kong, L., Nefedov, A., and Anandakrishnan, A. (2020).  Using generative adversarial networks to synthesize artificial financial datasets.  arXiv preprint arXiv:2002.02271.
  * Engelke et al., (2019) Engelke, S., de Fondeville, R., and Oesting, M. (2019).  Extremal behaviour of aggregated data with an application to downscaling.  Biometrika, 106:127–144.
  * Engelke and Hitz, (2020) Engelke, S. and Hitz, A. S. (2020).  Graphical models for extremes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(4):871–932.
  * Engelke and Ivanovs, (2021) Engelke, S. and Ivanovs, J. (2021).  Sparse structures for multivariate extremes.  Annual Review of Statistics and Its Application, 8:241–270.
  * Engelke et al., (2015) Engelke, S., Malinowski, A., Kabluchko, Z., and Schlather, M. (2015).  Estimation of Hüsler–Reiss distributions and Brown–Resnick processes.  J. R. Stat. Soc. Ser. B Stat. Methodol., 77:239–265.
  * Field et al., (2012) Field, C. B., Barros, V., and Stocker, T. F. (2012).  Cambridge University Press, Cambridge.
  * Fukushima, (1980) Fukushima, K. (1980).  Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position.  Biological cybernetics, 36:193–202.
  * Gidel et al., (2018) Gidel, G., Berard, H., Vignoud, G., Vincent, P., and Lacoste-Julien, S. (2018).  A variational inequality perspective on generative adversarial networks.  arXiv preprint arXiv:1802.10551.
  * Goodfellow et al., (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. (2014).  Generative adversarial nets.  In Advances in neural information processing systems, pages 2672–2680.
  * Gudendorf and Segers, (2010) Gudendorf, G. and Segers, J. (2010).  Extreme-value copulas.  In Copula Theory and Its Applications, pages 127–145. Springer.
  * Hazeleger et al., (2012) Hazeleger, W., Wang, X., Severijns, C., Ştefănescu, S., Bintanja, R., Sterl, A., Wyser, K., Semmler, T., Yang, S., Van den Hurk, B., et al. (2012).  EC-Earth v2. 2: description and validation of a new seamless earth system prediction model.  Climate dynamics, 39(11):2611–2629.
  * Hillier and Dixon, (2020) Hillier, J. K. and Dixon, R. S. (2020).  Seasonal impact-based mapping of compound hazards.  Environmental Research Letters, 15(11):114013.
  * Hosking, (1985) Hosking, J. (1985).  Algorithm as 215: Maximum-likelihood estimation of the parameters of the generalized extreme-value distribution.  Journal of the Royal Statistical Society. Series C (Applied Statistics), 34(3):301–310.
  * Hosking et al., (1985) Hosking, J., Wallis, J. R., and Wood, E. F. (1985).  Estimation of the generalized extreme-value distribution by the method of probability-weighted moments.  Technometrics, 27(3):251–261.
  * Huser and Genton, (2016) Huser, R. and Genton, M. G. (2016).  Non-stationary dependence structures for spatial extremes.  J. Agric. Biol. Environ. Stat., 21(3):470–491.
  * Hüsler and Reiss, (1989) Hüsler, J. and Reiss, R.-D. (1989).  Maxima of normal random vectors: between independence and complete dependence.  Statist. Probab. Lett., 7:283–286.
  * Huster et al., (2021) Huster, T., Cohen, J. E. J., Lin, Z., Chan, K., Kamhoua, C., Leslie, N., Chiang, C.-Y. J., and Sekar, V. (2021).  Pareto gan: Extending the representational power of gans to heavy-tailed distributions.
  * Ioffe and Szegedy, (2015) Ioffe, S. and Szegedy, C. (2015).  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  arXiv preprint arXiv:1502.03167.
  * Jalalzai et al., (2018) Jalalzai, H., Clémençon, S., and Sabourin, A. (2018).  On binary classification in extreme regions.  In NeurIPS, pages 3096–3104.
  * Jongman et al., (2014) Jongman, B., Hochrainer-Stigler, S., Feyen, L., Aerts, J. C. J. H., Mechler, R., Botzen, W. J. W., Bouwer, L. M., Pflug, G., Rojas, R., and Ward, P. J. (2014).  Increasing stress on disaster-risk finance due to large floods.  Nature Climate Change, 4(4):264–268.
  * Kabluchko et al., (2009) Kabluchko, Z., Schlather, M., and de Haan, L. (2009).  Stationary max-stable fields associated to negative definite functions.  Ann. Probab., 37:2042–2065.
  * Karras et al., (2017) Karras, T., Aila, T., Laine, S., and Lehtinen, J. (2017).  Progressive growing of gans for improved quality, stability, and variation.  arXiv preprint arXiv:1710.10196.
  * Keggenhoff et al., (2014) Keggenhoff, I., Elizbarashvili, M., Amiri-Farahani, A., and King, L. (2014).  Trends in daily temperature and precipitation extremes over georgia, 1971–2010.  Weather and Climate Extremes, 4:75–85.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Le et al., (2018) Le, P. D., Davison, A. C., Engelke, S., Leonard, M., and Westra, S. (2018).  Dependence properties of spatial rainfall extremes and areal reduction factors.  Journal of hydrology, 565:711–719.
  * Leonard et al., (2014) Leonard, M., Westra, S., Phatak, A., Lambert, M., van den Hurk, B., McInnes, K., Risbey, J., Schuster, S., Jakob, D., and Stafford-Smith, M. (2014).  A compound event framework for understanding extreme impacts.  Wiley Interdisciplinary Reviews: Climate Change, 5(1):113–128.
  * Link et al., (2019) Link, R., Snyder, A., Lynch, C., Hartin, C., Kravitz, B., and Bond-Lamberty, B. (2019).  Fldgen v1. 0: an emulator with internal variability and space–time correlation for earth system models.  Geoscientific Model Development (Online), 12(PNNL-SA-132926).
  * McKinnon and Deser, (2018) McKinnon, K. A. and Deser, C. (2018).  Internal variability and regional climate trends in an observational large ensemble.  Journal of Climate, 31(17):6783–6802.
  * McKinnon et al., (2017) McKinnon, K. A., Poppick, A., Dunn-Sigouin, E., and Deser, C. (2017).  An “observational large ensemble” to compare observed and modeled temperature trend uncertainty due to internal variability.  Journal of Climate, 30(19):7585–7598.
  * Mehrabi and Ramankutty, (2019) Mehrabi, Z. and Ramankutty, N. (2019).  Synchronized failure of global crop production.  Nature Ecology & Evolution, 3(5):780–786.
  * Meinshausen, (2006) Meinshausen, N. (2006).  Quantile regression forests.  Journal of Machine Learning Research, 7(Jun):983–999.
  * Morice et al., (2012) Morice, C. P., Kennedy, J. J., Rayner, N. A., and Jones, P. D. (2012).  Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 data set.  Journal of Geophysical Research: Atmospheres, 117(D8).
  * Naveau et al., (2020) Naveau, P., Hannart, A., and Ribes, A. (2020).  Statistical methods for extreme event attribution in climate science.  Annual Review of Statistics and Its Application, 7:89–110.
  * Poon et al., (2004) Poon, S.-H., Rockinger, M., and Tawn, J. (2004).  Extreme value dependence in financial markets: Diagnostics, models, and financial implications.  The Review of Financial Studies, 17:581–610.
  * Radford et al., (2015) Radford, A., Metz, L., and Chintala, S. (2015).  Unsupervised representation learning with deep convolutional generative adversarial networks.  arXiv preprint arXiv:1511.06434.
  * Raymond et al., (2020) Raymond, C., Horton, R. M., Zscheischler, J., Martius, O., AghaKouchak, A., Balch, J., Bowen, S. G., Camargo, S. J., Hess, J., Kornhuber, K., Oppenheimer, M., Ruane, A. C., Wahl, T., and White, K. (2020).  Understanding and managing connected extreme events.  Nature Climate Change, 10(7):611–621.
  * Resnick, (2008) Resnick, S. I. (2008).  Extreme values, regular variation and point processes.  Springer, New York.
  * Sahu et al., (2019) Sahu, S., Gupta, R., and Espy-Wilson, C. (2019).  Modeling feature representations for affective speech using generative adversarial networks.  arXiv preprint arXiv:1911.00030.
  * Schlather and Tawn, (2003) Schlather, M. and Tawn, J. A. (2003).  A dependence measure for multivariate and spatial extreme values: Properties and inference.  Biometrika, 90:139–156.
  * Schlegl et al., (2017) Schlegl, T., Seeböck, P., Waldstein, S. M., Schmidt-Erfurth, U., and Langs, G. (2017).  Unsupervised anomaly detection with generative adversarial networks to guide marker discovery.  In International conference on information processing in medical imaging, pages 146–157. Springer.
  * Seneviratne et al., (2021) Seneviratne, S. I., Zhang, X., Adnan, M., Badi, W., Dereczynski, C., Di Luca, A., Ghosh, S., Iskandar, I., Kossin, J., Lewis, S., Otto, F., Pinto, I., Satoh, M., Vicente-Serrano, S. M., Wehner, M., and Zhou, B. (2021).  Weather and Climate Extreme Events in a Changing Climate.  In Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., editors, Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, and New York, NY, USA.
  * Srivastava et al., (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014).  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research, 15:1929–1958.
  * Tanarhte et al., (2015) Tanarhte, M., Hadjinicolaou, P., and Lelieveld, J. (2015).  Heat wave characteristics in the eastern mediterranean and middle east using extreme value theory.  Climate Research, 63:99–113.
  * Tschumi et al., (2021) Tschumi, E., Lienert, S., van der Wiel, Joos, F., and Zscheischler, J. (2021).  A climate database with varying drought-heat signatures for climate impact modelling.  Geoscience Data Journal.
  * Van der Wiel et al., (2021) Van der Wiel, K., Lenderink, G., and De Vries, H. (2021).  Physical storylines of future European drought events like 2018 based on ensemble climate modelling.  Weather and Climate Extremes, 33:100350.
  * (64) Van der Wiel, K., Stoop, L. P., Van Zuijlen, B., Blackport, R., Van den Broek, M., and Selten, F. (2019a).  Meteorological conditions leading to extreme low variable renewable energy production and extreme high energy shortfall.  Renewable and Sustainable Energy Reviews, 111:261–275.
  * (65) Van der Wiel, K., Wanders, N., Selten, F., and Bierkens, M. (2019b).  Added value of large ensemble simulations for assessing extreme river discharge in a 2 C warmer world.  Geophysical Research Letters, 46(4):2093–2102.
  * Van Kempen et al., (2021) Van Kempen, G., Van der Wiel, K., and Melsen, L. A. (2021).  The impact of hydrological model structure on the simulation of extreme runoff events.  Natural Hazards and Earth System Sciences, 21(3):961–976.
  * Velthoen et al., (2021) Velthoen, J., Dombry, C., Cai, J.-J., and Engelke, S. (2021).  Gradient boosting for extreme quantile regression.
  * Vogel et al., (2021) Vogel, J., Rivoire, P., Deidda, C., Rahimi, L., Sauter, C. A., Tschumi, E., van der Wiel, K., Zhang, T., and Zscheischler, J. (2021).  Identifying meteorological drivers of extreme impacts: an application to simulated crop yields.  Earth System Dynamics, 12(1):151–172.
  * Wadsworth and Tawn, (2012) Wadsworth, J. L. and Tawn, J. A. (2012).  Dependence modelling for spatial extremes.  Biometrika, 99(2):253–272.
  * Wiese et al., (2019) Wiese, M., Knobloch, R., and Korn, R. (2019).  Copula & marginal flows: Disentangling the marginal from its joint.
  * Yoon et al., (2010) Yoon, S., Cho, W., Heo, J.-H., and Kim, C. E. (2010).  A full bayesian approach to generalized maximum likelihood estimation of generalized extreme value distribution.  Stochastic environmental research and risk assessment, 24(5):761–770.
  * Zhang et al., (2017) Zhang, H., Xu, T., Li, H., Zhang, S., Wang, X., Huang, X., and Metaxas, D. N. (2017).  Stackgan: Text to photo-realistic image synthesis with stacked generative adversarial networks.  In Proceedings of the IEEE international conference on computer vision, pages 5907–5915.
  * Zheng et al., (2019) Zheng, P., Yuan, S., Wu, X., Li, J., and Lu, A. (2019).  One-class adversarial nets for fraud detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1286–1293.
  * Zscheischler et al., (2020) Zscheischler, J., Martius, O., Westra, S., Bevacqua, E., R., C., Horton, R. M., van den Hurk, B., AghaKouchak, A., Jézéquel, A., Mahecha, M. D., Maraun, D., Ramos, A. M., Ridder, N., Thiery, W., and Vignotto, E. (2020).  A typology of compound weather and climate events.  Nature Reviews Earth & Environment, pages 333–347.
  * Zscheischler and Seneviratne, (2017) Zscheischler, J. and Seneviratne, S. I. (2017).  Dependence of drivers affects risks associated with compound events.  Science Advances, 3(6):e1700263.
  * Zscheischler et al., (2018) Zscheischler, J., Westra, S., Van Den Hurk, B. J., Seneviratne, S. I., Ward, P. J., Pitman, A., AghaKouchak, A., Bresch, D. N., Leonard, M., Wahl, T., and Zhang, X. (2018).  Future climate risk from compound events.  Nature Climate Change, 8(6):469–477.