<figure><img src="image_1.png" /></figure>

HAL
open science

# ExceedGAN: Simulation above extreme thresholds using Generative Adversarial Networks

Michaël Allouche, Stéphane Girard, Emmanuel Gobet

<figure><img src="image_2.png" /></figure>

## To cite this version:

Michaël Allouche, Stéphane Girard, Emmanuel Gobet. ExceedGAN: Simulation above extreme thresholds using Generative Adversarial Networks. Extremes, In press. hal-05044516v2

HAL Id: hal-05044516

https://inria.hal.science/hal-05044516v2

Submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ExceedGAN: Simulation above extreme thresholds using Generative Adversarial Networks

Michaël Allouche¹†, Stéphane Girard²*†, Emmanuel Gobet³†

1,3Kaiko - Quantitative Data, 33 rue du Louvre, Paris, 75002, France.

2* Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, Grenoble, 38000, France.

³LPSM, Sorbonne Université, Case courrier 158, 4 Place Jussieu, Paris Cedex, 75252, France.

*Corresponding author(s). E-mail(s): stephane.girard@inria.fr;

Contributing authors: michael.allouche@kaiko.com;

emmanuel.gobet@sorbonne-universite.fr;

†These authors contributed equally to this work.

## Abstract

This paper devises a novel neural-inspired approach for simulating multivariate extremes. Specifically, we propose a GAN-based generative model for sampling multivariate data exceeding large thresholds, giving rise to what we refer to as the ExceedGAN algorithm. Our approach is based on approximating marginal log-quantile functions using feedforward neural networks with eLU activation functions specifically introduced for bias correction. An error bound is provided on the margins, assuming a *J*th order condition from extreme value theory. The numerical experiments illustrate that ExceedGAN outperforms competitors, both on synthetic and real-world data sets.

**Keywords:** Extremes, Generative AI, GAN (Generative Adversarial Network), Neural networks, Bias correction

**MSC Classification:** 62G32, 68T07, 62G09
# 1 Introduction

## Context and statement of the problem

This paper aims at investigating the use of generative methods based on neural networks to simulate extreme events and to assess the associated risks. Risk management is indeed an activity of the utmost importance to ensure the long-term sustainability and resilience of systems of any nature, and in this context, assessing the risk associated with extreme events is a major concern; see applications in insurance [1], finance [2], engineering [3] and health [4] for instance. These extreme events are rare by nature, leading to a few observations only and consequently, making the risk assessment difficult in practice. In this work, we focus on the distribution of the excess of a D-dimensional random vector $X$ given that any marginal of $X$ exceeds a large threshold $u$, see Figure 1 for an illustration in dimension $D = 2$. This is usually referred to as the exceedance distribution. In the one-dimensional case $D = 1$, it is well known that, for a large enough threshold $u$, the random variable $X$ given that $X > u$ approximately follows a Generalized Pareto Distribution [5, Theorem 4.1] including position, scale and shape parameters to be estimated, for instance using a numerical maximization of the log-likelihood. In a general multi-dimensional setting, the description of the conditional distribution $X | X > u$ is much more difficult and simple approximations are no longer available (here and in the sequel, the inequality is meant componentwise). See de Haan and Ferreira [6, Chapter 8] for an introduction to the statistical modeling of multivariate extremes. Our purpose is to design a conditional generative model for such a situation, with a single generative model available for any extreme threshold $u$ represented as a continuous conditioning variable. We specialize our study to the challenging setting of heavy-tailed $X$.

Generally speaking, to compute tail statistics on the distribution of $X | X > u$, two families of methods are usually considered. First, one can rely on extreme value theory and use extrapolation procedures to analytically compute extreme risk metrics that are related to the tail. We refer for instance to Weissman estimator dedicated to extreme quantiles [7] and its adaptations to other risk measures such as expectiles [8] or extremiles [9] in the unidimensional framework. An example of extension to the $D > 1$ framework can be found in [10]. Alternatively, we focus on a second line of works which consists in generating additional extreme samples to compute the tail statistics of interest. One possibility is to use a stochastic model and select a sampling scheme within the various existing ones (importance sampling [11, Chapter 4], MCMC with splitting [12], or interacting particles system [13]). However, these methods suffer from a high computational complexity (especially in high dimension) and from the problem of tuning critical hyper-parameters. In addition, these methods assume that a stochastic model is already calibrated, which can be a major issue in case of scarce observations. Another possibility is to use methods from generative artificial intelligence. They allow to get extra samples in an efficient way [14] by leveraging the available samples and without recourse to a pre-specified stochastic model. They are based on the approximate equality in distribution:

$$X \stackrel{d}{\sim} G(Z), \qquad (1)$$
<figure><img src="image_4.png" /><figcaption>**Fig. 1:** Scatterplot of $n = 1000$ simulated observations in dimension $D = 2$ from a bivariate Gumbel copula (dependence parameter $\mu = 2$) and Burr margins parametrized by $(\gamma_1 = 0.5, \rho_1 = -1)$ and $(\gamma_2 = 0.5, \rho_2 = -2)$; see Section 4.1 for details. Two extreme regions $Q(\delta_n)$ are represented by green $(\delta_n = (0.1, 0.1)^\top)$ and red $(\delta_n = (0.05, 0.05)^\top)$ dashed lines; see (21) for their definitions.</figcaption></figure>

for a given (easy-to-sample) latent random variable $Z$ and where $G$ is a (possibly quite
complicate) generating function to learn from observations. This is the path we are
taking in this work. The most popular neural network (NN) generative models are
the Variational Auto-Encoder (VAE) [15] and the Generative Adversarial Network
(GAN) [16], on which we shall focus, see Section 2.3 below and [17, Chapter 17]
for introductions to GANs. More recent models such as Normalizing flows [18] and
Diffusion models [19, 20] have gained some popularity. However, the authors of [21]
have demonstrated both theoretically and from a practical point of view that GANs
cannot simulate data in the distribution tails without a proper adaptation to the
extreme value framework: It means that (1) holds in the bulk of the distribution but
hardly in the tails, see Section 2.4 hereafter for details. To this end, they introduced
the EV-GAN method which is designed to deal with heavy-tailed distributions for
$X$, based on first-order and second-order characterizations of the Fréchet maximum
domain of attraction [6, Section 2.3]. Nonetheless, observe that designing a generative
model for $X$ is not equivalent to designing a generative model for $X \mid X > u$: of
course, one could sample $X$ and then, owing to the acceptance-rejection algorithm, the
samples above $u$ would stand for samples of $X \mid X > u$. This procedure is yet far from
being efficient, especially for extreme events, since only a fraction $\mathbb{P}\{X > u\} \ll 1$ of
the samples are kept.
## Main contributions

Our key contributions are threefold:

* We design a generative model for the excess $X | X > u$ over a given $u$ - called fixed-level ExceedGAN. It is based on a feedforward neural network with eLU activation functions dedicated to bias correction under a $J$th order condition from extreme value theory.

* We establish an accuracy result (Theorem 3.2) on the log-quantile of the excess over the threshold $u$. It is established that it can be uniformly approximated by a neural network with $J(J - 1)$ neurons. The approximation rate is also provided.

* We extend the previous generative model to any $u$ with a single learning procedure referred to as level-varying ExceedGAN.

In our experiments, we show that the new ExceedGAN model outperforms the state-of-the-art model EV-GAN. In contrast with EV-GAN which uses ReLU activation functions, the new ExceedGAN makes use of eLU activation functions. This idea is inspired by the work of [22] which combines higher-order extreme value conditions and eLU NNs to estimate extreme quantiles with a significantly reduced bias.

## Outline

In Section 2, we briefly review the main concepts of generative modeling with neural networks, in preparation for our main contributions stated and proved in Section 3. Section 4 gathers the numerical experiments both on simulated and real cryptocurrency data sets.

# 2 Generative modeling

This section reviews the main ingredients in generative modeling that will be used to develop our method. For a more complete overview, see [23].

## 2.1 Theoretical framework

Given data $\{X_1, \dots, X_n\}$ assumed to be independently sampled from an unknown distribution $p_X$ on a measurable space $(\mathcal{X}, \mathcal{B}_\mathcal{X})$, the goal of generative modeling is to find a function $G: \mathcal{Z} \to \mathcal{X}$, called generator, and a probability distribution $p_Z$, called latent distribution and defined on some set $\mathcal{Z}$, such that

$$G(\mathbf{Z}) \stackrel{d}{=} \mathbf{X}, \quad \mathbf{Z} \sim p_{\mathbf{Z}}. \qquad (2)$$

The existence of $G$ and $p_Z$ is provided by Kuratowski's Theorem [24, Proposition 7.15] also called the measurable isomorphism Theorem [25, p. 7]: If $(\mathcal{Z}, p_Z)$ and $(\mathcal{X}, p_X)$ are two Polish probability spaces without atom, then there exists a (non-unique) measurable bijection $G: \mathcal{Z} \to \mathcal{X}$ such that $G_{\#}p_Z = p_X$ and $G_{\#}^{-1}p_X = p_Z$, where, for all $E \subset \mathcal{X}$, $G_{\#}p_Z(E) := p_Z(\{z \in \mathcal{Z}: G(z) \in E\}) = p_Z(G^{-1}(E))$ stands for the push-forward measure of $p_Z$ by $G$. From an optimal transport perspective [26], the
push-forward operator $G_{\#}$ transfers the mass of the probability measure $p_Z$ on $\mathcal{Z}$ into $\mathcal{X}$ by applying the map $G$.

In our work, $\mathbf{X}$ and $\mathbf{Z}$ take values respectively in the Euclidean spaces $\mathbb{R}^D$ and $\mathbb{R}^q$, where $q$ is a hyper-parameter of the generative method. We focus on a parametric family of generators $\{G_{\theta}\}_{\theta \in \Theta}$ where $\Theta \subset \mathbb{R}^p$ and $\{p_{\theta}\}_{\theta \in \Theta}$ denotes the set of associated parametric distributions such that $G_{\theta}(\mathbf{Z}) \sim p_{\theta}$. The problem reduces to finding the “best” parameter $\theta^*$ such that $p_{\theta^*}$ and $p_{\mathbf{X}}$ are as close as possible, or equivalently

$$ \mathbf{X} \stackrel{\text{d}}{\approx} G_{\theta^*}(\mathbf{Z}). \qquad (3) $$

Given $\{\mathbf{X}_1, \dots, \mathbf{X}_n\}$, the generative modeling framework mainly consists of combining two ingredients:

* The parametrization $G_{\theta}$ and the latent distribution $p_Z$ to use as inputs;

* The distance or the similarity criterion between the probability distributions $p_{\theta}$ and $p_{\mathbf{X}}$ (or its empirical version based on the data set), as well as the optimization process that defines the optimal $\theta$.

Kuratowski's theorem theoretically provides a large flexibility in the choice of the latent distribution $p_Z$. In practice, easy-to-simulate distributions are considered such as the Gaussian and the uniform distributions. As a consequence, the modeling effort is usually put on the generator $G_{\theta}: \mathbb{R}^q \to \mathbb{R}^D$, which is made of transformations of neural networks.

## 2.2 Neural networks

A neural network (NN) is a non-linear function built with a fixed number of neurons, each one representing a function, and distributed across several hidden layers. Neurons are scaled and translated in the network by parameters called, respectively, weights and biases. Among many different existing neural network architectures, we consider the classical one-hidden layer feedforward neural network $NN_{\theta}: \mathbb{R}^q \to \mathbb{R}$ composed of $K$ neurons

$$ NN_{\theta}(\mathbf{z}) = b^{(2)} + \sum_{k=1}^{K} w_{k}^{(2)} \sigma(\langle \mathbf{w}_{k}^{(1)}, \mathbf{z} \rangle + b_{k}^{(1)}), \qquad (4) $$

with parameters

$$ \theta := \{b^{(2)}\} \cup \{\mathbf{w}_k^{(1)}, w_k^{(2)}, b_k^{(1)}\}_{k=1}^K \in \Theta := \mathbb{R} \times (\mathbb{R}^q \times \mathbb{R} \times \mathbb{R})^K, $$

and where $\langle \cdot, \cdot \rangle$ is a scalar product on $\mathbb{R}^q$, and $\sigma: \mathbb{R} \to \mathbb{R}$ is a non-linear “activation” function. The total number of parameters is $p = K(q + 2) + 1$. A usual activation function is the Rectified Linear Unit (ReLU) given by $\sigma^R(x) = \max(x, 0)$, but in our extreme value setting, the eLU activation function

$$ \sigma^{\mathrm{E}}(x) = x\mathrm{I}(x > 0) + (e^x - 1)\mathrm{I}(x \le 0) \qquad (5) $$
plays a crucial role, as firstly noticed in [22]. The generator can be easily extended to D-dimensional outputs by replication of (4), that is by considering the D-dimensional vector $(NN_{\theta}^{(1)},...,NN_{\theta}^{(D)})$ with various uses of parameters $\theta$ for each coordinate. Note that, in the case of conditional generative modeling, the input dimension $q$ is increased by the number of conditioning variables: For level-varying ExceedGAN, the exceedance threshold $u$ is taken as an input, so that the input dimension of the NN is $q+D$.

The ability of $NN_{\theta}$, defined in (4), to approximate a given function $\varphi$ as $K \to \infty$, has been studied under various assumptions on $\varphi$ and with respect to different norms. Nowadays, the state-of-the-art result is the following [27]:

**Theorem 2.1** (Universal approximation theorem). Suppose $\varphi$ is a continuous function on a compact space $Z \subset \mathbb{R}^q$ and $\sigma$ is not a polynomial. Then, for any $\varepsilon > 0$, there exists a one-hidden layer neural network $NN_{\theta}$ (for some $K$ depending on $\varepsilon$) defined in (4) such that

$$ \sup_{z \in Z} |\varphi(z) - NN_{\theta}(z)| < \varepsilon. $$

The key assumption of the universal approximation theorem is that the function $\varphi$ to approximate is continuous on a compact set and thus bounded. We shall see in Section 2.4 the consequences for simulating heavy-tailed random variables in a generative modeling framework.

## 2.3 Generative Adversarial Networks (GANs)

A GAN aims at approximating the unknown generator $G: \mathbb{R}^q \to \mathbb{R}^D$ through a parametric family of non-linear functions $\{G_{\theta}: \mathbb{R}^q \to \mathbb{R}^D\}_{\theta \in \Theta}$ and to estimate the optimal parameter $\theta^*$ from a data set $\{X_1,...,X_n\}$ of independent samples from the unknown distribution $p_X$. The estimation is performed by optimizing an objective function which can be interpreted as an adversarial game between a generator and a discriminator also chosen in a parametric family of non-linear functions $\{D_{\phi}: \mathbb{R}^D \to [0,1]\}_{\phi \in \Phi}$, where $\Phi$ stands for the set in which the discriminator parameters take values. In other words, $D_{\phi}(X)$ represents the probability that an observation $x \in \mathbb{R}^D$ is drawn from $X \sim p_X$. Both the generator and the discriminator are non-linear functions (usually built from neural networks) with opposing objectives: The former tries to generate realistic data, while the latter tries to discriminate between synthetic generated data and real observations. In [16], this optimization problem is defined as:

$$ \arg\min_{\theta \in \Theta} \max_{\phi \in \Phi} (E_{p_X}\{\log D_{\phi}(X)\} + E_{p_Z}\{\log [1 - D_{\phi}\{G_{\theta}(Z)\}]\}), \quad (6) $$

where $Z$ denotes the latent variable, see (2). The goal of the discriminator is to maximize the criterion with respect to $\phi$ which is achieved when both $D_{\phi}(X) \simeq 1$ (realistic outputs) and $D_{\phi}\{G_{\theta}(Z)\} \simeq 0$ (fake outputs). In contrast, the generator aims at minimizing the above criterion with respect to $\theta$ in such a way that $G_{\theta}(Z) \stackrel{d}{=} X$ leading to $D_{\phi}\{G_{\theta}(Z)\} \simeq 1$ (increasing the realism of fake outputs). The analysis
of the statistical fluctuations of the estimators obtained by considering the empirical counterpart of the above optimization problem is conducted in [28]. While the GAN's optimization criteria provides theoretical guarantees to learn the appropriate data distribution [16, Theorem 1], finding the optimal generator and discriminator is known to be very challenging.

## 2.4 Generative modeling for extremes: Challenges and solutions

To exemplify the difficulty of designing a generative model for heavy tail distributions, let us first consider a real-valued random variable $X$ with cumulative distribution function $F_X$. In this one-dimensional setting $D = q = 1$, the problem (2) benefits from an explicit solution based on the quantile function $F_X^{-1}(z) := \inf\{x \in \mathbb{R} : F_X(x) \ge z\} =: q_X(z)$. The inversion method [29] shows that one can set $G := F_X^{-1}$ and $Z \sim U([0, 1])$. When $p_X$ has an infinite upper endpoint, $F_X^{-1}(z) \to +\infty$ as $z \to 1$ so that $F_X^{-1}$ does not fulfill the assumptions of Theorem 2.1: a neural network (4) cannot uniformly approximate $G$. Moreover, since $Z$ is a bounded random variable, when the activation function $\sigma$ is continuous, $NN_\theta(Z)$ is also a bounded random variable. Similarly, if $Z$ is a Gaussian variable and $\sigma$ is sub-polynomial (like a ReLU activation function), then $NN_\theta(Z)$ is dominated by a polynomial of Gaussian random variables and as such, it is light-tailed [30]. See [23] for illustrations on simulated and real data. To overcome these problems, an adapted GAN procedure based on distribution shifting is proposed in [31], to reduce the lack of training data in the tails. In addition to this heuristic, three main directions have been investigated to adapt GANs to heavy tails: A preprocessing of the data to get rid of the tail heaviness [32, 33], the use of heavy-tailed latent variables [34, 35], and the introduction of new parametrizations [21] to adapt the optimization problem (6) to the heavy-tail situation.

Our approach belongs to this third line of works. Similarly to [21], we propose in the next section a new parametrization of GANs adapted to the heavy-tail framework, but including a higher order approximation of the distribution tail by eLU neural networks. Before that, let us mention recent works dedicated to generative modeling for extremes using alternative architectures to GANs. On the VAE side, it is shown in [36, Corollary 7] that a VAE built with piecewise linear activation functions and Gaussian distributions cannot reproduce heavy-tailed margins. The authors thus consider a univariate heavy-tailed distribution to sample a radius, and, conditionally on the latter, a direction is sampled from a multivariate normal distribution. The product of the two yields the desired multivariate regularly-varying vector. In [37], a VAE is designed to sample spatial extremes with flexible and non-stationary dependence structures. Normalizing flows have been used instead in [38] to model the marginals by transforming Gaussian tails into generalized Pareto distribution (GPD) tails while [39, 40] take advantage of geometric approaches to model multivariate extremes, capturing both the marginal distributions and the dependence structure. See also [41] for a comparison study of different neural network generative models dedicated to multivariate extremes through the modeling of angular variables.
# 3 Bias-reduced simulation of extremes with GANs

We begin in Section 3.1 by describing the approximation result underpinning the simulation of exceedances in the univariate case with a GAN using eLU activation functions. Then, it is shown in Section 3.2 how this principle is extended to the multivariate case of interest.

## 3.1 Univariate neural network approximation

Let us start with the one-dimensional framework of Section 2.4. Given an independent sample {$X_1, \dots, X_n$} from $F_X$, we focus on the simulation using GANs of exceedances from $X$ defined by

$$Y \stackrel{d}{=} X \text{ given } X > u. \qquad (7)$$

We consider an extreme threshold $u$ and to emphasize the extreme behavior, we allow $u \to +\infty$ in a way it can depend on the size $n$ of the data set:

$$u_n = F_X^{-1}(1 - \delta_n) \qquad (8)$$

where $\delta_n \in (0, 1)$ is a probability sequence such that $\delta_n \to 0$ as $n \to \infty$. Hence, now we focus on $Y = Y(\delta_n)$ whose simulation can rely on the inversion method [29], as shown in the next lemma.

**Lemma 3.1.** Assume that the distribution of $X$ is continuous. Then, for all $\delta_n \in (0, 1)$,

$$Y(\delta_n) \stackrel{d}{=} q_{Y(\delta_n)}(1 - Z) = F_X^{-1}(1 - \delta_n Z),$$

with $Z \sim \mathcal{U}([0, 1])$.

*Proof of Lemma 3.1.* For any $y > F_X^{-1}(1 - \delta_n)$,

$$\begin{aligned} \bar{F}_{Y(\delta_n)}(y) &:= \mathbb{P}\{Y(\delta_n) > y\} = \mathbb{P}\{X > y \mid X > F_X^{-1}(1 - \delta_n)\} \\ &= \frac{\mathbb{P}\{X > y\}}{\mathbb{P}\{X > F_X^{-1}(1 - \delta_n)\}} = \frac{\bar{F}_X(y)}{\delta_n}, \end{aligned}$$

while $\bar{F}_{Y(\delta_n)}(y) = 1$ if $y \le F_X^{-1}(1 - \delta_n)$. Similarly, for any $y > F_X^{-1}(1 - \delta_n)$,

$$\bar{F}_{F_X^{-1}(1-\delta_n Z)}(y) := \mathbb{P}\{F_X^{-1}(1-\delta_n Z) > y\} = \mathbb{P}\{Z < \frac{\bar{F}_X(y)}{\delta_n}\} = \frac{\bar{F}_X(y)}{\delta_n},$$

and $\bar{F}_{F_X^{-1}(1-\delta_n Z)}(y) = 1$ provided $y \le F_X^{-1}(1 - \delta_n)$. The conclusion follows. $\square$

In practice, one may then define a GAN with generator based on $\hat{F}_{Y(\delta_n)}(\cdot)$, an estimator of $F_{Y(\delta_n)}^{-1}(\cdot)$, and a latent variable $Z \sim \mathcal{U}([0, 1])$ by setting $\tilde{Y}(\delta_n) = \hat{F}_X^{-1}(1 - \delta_n Z)$. Yet, when $n\delta_n \to c \ge 0$ as $n \to \infty$, the sample maximum $\max(X_1, \dots, X_n)$ has an asymptotic probability $\exp(-zc)$ to be smaller than $F_X^{-1}(1 - \delta_n z)$, for any $z \in (0, 1)$. Consequently, for small values of $z$, $F_X^{-1}(1 - \delta_n z)$ is an extreme quantile, likely to be larger than the sample maximum, whose estimation requires dedicated extrapolation techniques.
Assuming that $F_X$ is in the maximum domain of attraction of Fréchet, it is known from [6, Theorem 1.2.1] and [6, Proposition B.1.9.9] that the tail quantile function

$$U_X(t) := F_X^{-1}(1 - 1/t),$$

defined for all $t > 1$, is regularly-varying with index $\gamma > 0$ (this property is denoted by $U_X \in RV_\gamma$ in the sequel), i.e.

$$U_X(t) = t^\gamma L(t),$$

where $\gamma$ is the so-called tail-index and $L \in RV_0$ is a slowly-varying function at infinity, i.e. $L$ is positive and, for all $y > 0$,

$$\lim_{t \to \infty} \frac{L(yt)}{L(t)} = 1. \qquad (9)$$

The index $\gamma$ tunes the tail heaviness of $F_X$: the larger the index, the heavier the right tail. Examples include the (generalized) Pareto, Burr, Fisher, Inverse gamma, and Student $t$ distributions. The idea underpinning the estimation is to take advantage of $U_X \in RV_\gamma$ to establish a link between the extreme quantiles of interest

$$q_{Y(\delta_n)}(1-z) = F_X^{-1}(1 - \delta_n z) = U_X(1/(\delta_n z)),$$

and the threshold $u_n = F_X^{-1}(1-\delta_n) = U_X(1/\delta_n)$. To this end, introduce the log-spacing function defined as

$$f(x_1, x_2) := \log U_X(\exp(x_1 + x_2)) - \log U_X(\exp(x_2)) = \gamma x_1 + \varphi(x_1, x_2), \quad (10)$$

with $x_1, x_2 \ge 0$ and

$$\varphi(x_1, x_2) := \log L(\exp(x_1 + x_2)) - \log L(\exp(x_2)). \qquad (11)$$

Considering $x_1 = \log(1/z)$ and $x_2 = \log(1/\delta_n)$, it immediately follows that the quantity of interest in Lemma 3.1 can be rewritten as

$$q_{Y(\delta_n)}(1 - \cdot) = F_X^{-1}(1 - \delta_n \cdot) = F_X^{-1}(1 - \delta_n) \exp(f(\log(1/\cdot), \log(1/\delta_n))). \quad (12)$$

In order to build high order approximations of $\varphi$ defined in (11), consider a $J$-th order condition, introduced in [42] for all $J \ge 2$, on the tail quantile function. Assume there exist $\gamma > 0$ and, for all $j \in \{2, \dots, J\}$, $\rho_j \le 0$ as well as positive or negative functions $A_j$ with $A_j(t) \to 0$ as $t \to \infty$, $|A_j| \in RV_{\rho_j}$, such that

$$\log U_X(ty) - \log U_X(t) = \gamma \log y + \sum_{j=2}^{J} \prod_{\ell=2}^{j} A_{\ell}(t) R_j(y) + o\left(\prod_{j=2}^{J} A_j(t)\right), \quad (13)$$
as $t \to \infty$ for all $y > 0$, where:

$$R_j(y) = \int_1^y y^{\rho_2-1} \int_1^{y_2} y^{\rho_3-1} \dots \int_1^{y_j-1} y^{\rho_j-1} dy_j \dots dy_3 dy_2, \quad (14)$$

and suppose that

$$A_j(t) = c_j t^{\rho_j}, \quad (15)$$

where $c_j \neq 0$ and $\rho_j < 0$ for $j \in \{2, \dots, J\}$. Clearly, when $J = 2$, one recovers the second-order condition [6, Section 2.3]. Moreover, $J = 3$ and $J = 4$ yield back respectively the third-order [43, 44] and fourth-order conditions [45]. In the following, we show that the key parameter to quantify the approximation convergence rate is

$$\bar{\rho}_J := \rho_2 + \dots + \rho_J. \quad (16)$$

A NN approximation of $\varphi$ is introduced using $J(J-1)$ eLU activation functions (5):

$$\varphi^{NN_J}(x_1, x_2; \theta) := \sum_{i=1}^{J(J-1)/2} w_i^{(1)} \left\{ \sigma^E \left( w_i^{(2)} x_1 + w_i^{(3)} x_2 \right) - \sigma^E \left( w_i^{(4)} x_2 \right) \right\}, \quad (17)$$

$x_1, x_2 \ge 0$, for some $\theta = \{(w_i^{(1)}, w_i^{(2)}, w_i^{(3)}, w_i^{(4)}), i \in 1, \dots, J(J-1)/2\} \in \Theta := (\mathbb{R} \times \mathbb{R}_-^3)^{J(J-1)/2}$. The NN approximation of the log-spacing function (10) is then defined as

$$f^{NN_J}(x_1, x_2; \phi) := w^{(0)} x_1 + \varphi^{NN_J}(x_1, x_2; \theta), \quad (18)$$

where $\phi := (w^{(0)}, \theta) \in \Phi := \mathbb{R}_+ \times \Theta$. Combining (12) with (18), the NN approximation of the extreme quantile is defined as

$$q_{Y(\delta_n)}^{NN_J}(1-z; \phi) := F_X^{-1}(1-\delta_n) \exp \left\{ f^{NN_J}(\log(1/z), \log(1/\delta_n); \phi) \right\}. \quad (19)$$

The approximating NN includes $d = J(J-1)$ neurons and $p \le CJ^2$ parameters, where $C > 0$ is a constant independent of $J$. The next result states that a proper choice of parameter $\phi^*$ yields an accurate approximation of $q_{Y(\delta_n)}(1-z)$, uniformly in $z$, with a controlled error.

**Theorem 3.2.** Assume the $J$th order condition (13) holds together with (15) for some $J \ge 2$. There exists a parameter $\phi^* \in \Phi$, and a threshold $t_0 \in (0, 1)$ such that the one hidden-layer NN approximation (19) with $\phi = \phi^*$ of the extreme quantile of $Y(\delta_n)$ defined in (7) verifies

$$\sup_{z \in (0,1]} \left| \log q_{Y(\delta_n)}(1-z) - \log q_{Y(\delta_n)}^{NN_J}(1-z; \phi^*) \right| \le |\bar{\rho}_J \bar{c}_J| \delta_n^{-\bar{\rho}_J},$$
for all $0 < \delta_n \le t_0$, and where $\bar{c}_J = c_2 \times \cdots \times c_j$.

*Proof of Theorem 3.2.* We start by recalling Allouche et al. [22, Lemma 6], with a slight adaptation that follows from an easy inspection of its proof.

**Lemma 3.3.** Assume the Jth order condition (13) holds together with (15) for some $J \ge 2$. Let $\bar{\rho}_J = \rho_2 + \cdots + \rho_J$ and $\bar{c}_J = c_2 \times \cdots \times c_J$. There exists a parameter $\theta^* \in \Theta$ such that, for all $\varepsilon > 0$, there exists $x_\varepsilon > 0$ such that

$$|\varphi(x_1, x_2) - \varphi^{NN_J}(x_1, x_2; \theta^*)| \le \varepsilon \exp(x_1(\bar{\rho}_J + \varepsilon)) |\bar{c}_J| \exp(x_2\bar{\rho}_J),$$

for all $x_1 \ge 0$ and $x_2 \ge x_\varepsilon$.

Choose arbitrarily an $\varepsilon \in (0, -\bar{\rho}_J)$ and define $x_\varepsilon$ as in the above lemma. Choose $w^{(0),*} = \gamma$ and set $\phi^* = (w^{(0),*}, \theta^*)$. We have, in view of the triangle inequality, for all $z \in (0, 1]$,

$$\begin{aligned}
& \left| \log F_X^{-1}(1 - \delta_n z) - \log q_{Y(\delta_n)}^{NN_J}(1 - z; \phi^*) \right| \\
& = \left| f(\log(1/z), \log(1/\delta_n)) - f^{NN_J}(\log(1/z), \log(1/\delta_n); \phi^*) \right| \\
& \le \left| \gamma - w^{(0),*} \right| \log(1/z) + \left| \varphi(\log(1/z), \log(1/\delta_n)) - \varphi^{NN_J}(\log(1/z), \log(1/\delta_n); \theta^*) \right| \\
& \le \varepsilon z^{-(\bar{\rho}_J + \varepsilon)} |\bar{c}_J| \delta_n^{-\bar{\rho}_J},
\end{aligned}$$

for any $\delta_n \le \exp(-x_\varepsilon)$. Remarking that $z < 1$ gives the announced error bound. $\square$

The log-quantile of the excess over the threshold $u_n = F_X^{-1}(1 - \delta_n)$ can thus be uniformly approximated at the $\delta_n^{-\bar{\rho}_J}$ rate by a neural network with $J(J-1)$ neurons.

## 3.2 Multivariate neural network simulation

In view of (19), in the one-dimensional framework ($D = q = 1$), we define the NN extreme quantile estimator as

$$\hat{q}_{Y(\delta_n)}^{NN_J}(1 - z; \hat{\phi}) := X_{n-k+1,n} \exp \left\{ f^{NN_J}(\log(1/z), \log(1/\delta_n); \hat{\phi}) \right\}, \quad (20)$$

where $\hat{\phi}$ is the estimator of $\phi^*$ computed through a GAN-type optimization problem, see below. Here, $k = \lfloor n\delta_n \rfloor$ so that the order statistic $X_{n-k+1,n}$ is the empirical counterpart of the intermediate quantile $u_n = F_X^{-1}(1 - \delta_n)$. Combining with Lemma 3.1, one readily obtains a generator of exceedances in the set $[u_n, \infty)$. In the multidimensional setting $q > 1$ and $D > 1$, the multivariate threshold $\mathbf{u}_n$ is parametrized by component-wise probability levels $\delta_n^{(m)}$,

$$u_n^{(m)} = F_{X^{(m)}}^{-1}(1 - \delta_n^{(m)}), \quad m \in \{1, \dots, D\},$$

and the goal is thus to generate exceedances in the upper orthant

$$Q(\delta_n) = \{\mathbf{x} \in \mathbb{R}^D : x^{(1)} > F_{X^{(1)}}^{-1}(1 - \delta_n^{(1)}), \dots, x^{(D)} > F_{X^{(D)}}^{-1}(1 - \delta_n^{(D)})\}. \quad (21)$$
The strategy of both approximation and estimation steps consists in preserving the same parametric form for each marginal component, and in mixing the latent components to generate dependence between the D coordinates: the mth coordinate of the ExceedGAN generator is written as, with $\mathbf{z} = (z^{(1)}, \dots, z^{(q)})$,

$$ G_{\hat{\phi}, \delta_n^{(m)}}^{\text{EX}, (m)}(\mathbf{z}) = X_{n-k^{(m)}+1, n}^{(m)} \exp \left( \sigma^\mathbb{R} \left[ \sum_{j=1}^q \hat{w}_{j,m}^{(0)} \log(1/z^{(j)}) \right. \right. \\ \left. + \sum_{i=1}^{J^{(m)}(J^{(m)}-1)/2} \hat{w}_{i,m}^{(1)} \left\{ \sigma^\mathbb{E} \left( \sum_{j=1}^q \hat{w}_{i,j,m}^{(2)} \log(1/z^{(j)}) + \sum_{m'=1}^D \hat{w}_{i,m'}^{(3)} \log(1/\delta_n^{(m')}) \right) \right. \right. \\ \left. - \sigma^\mathbb{E} \left( \sum_{m'=1}^D \hat{w}_{i,m'}^{(4)} \log(1/\delta_n^{(m')}) \right) \right\} \right] \right), \quad (22) $$

with $k^{(m)} = \lfloor n\delta_n^{(m)} \rfloor$ and where $J^{(m)}$ denotes the order of condition (13) assumed on the mth margin, $m \in \{1, \dots, D\}$. In the fixed-level ExceedGAN method, the estimated parameter $\hat{\phi}$ is computed by plugging the generator (22) into the GAN optimization problem (6), while in the level-varying ExceedGAN method, the optimization framework associated with a conditional GAN [46] is adopted:

$$ \arg\min_{\theta \in \Theta} \max_{\phi \in \Phi} \left( E_{p_U} \{ E_{p_{Y(u_n)}} \{ \log D_{\phi}^{\text{EX}}(\mathbf{Y}, \delta_n) \} \} \right. \\ \left. + E_{p_U} \{ E_{p_Z} \{ \log [1 - D_{\phi}^{\text{EX}} \{ G_{\theta}^{\text{EX}}(\mathbf{Z}, \delta_n), \delta_n \} ] \} \} \right). \quad (23) $$

Here, $p_U$ is the uniform distribution on the hypercube $(a, 1)^D$ with $a \in (0, 1)$ to be chosen by the user while $p_{Y(u_n)}$ denotes the distribution of $\mathbf{Y}(u_n)$. See Figure 2 for an illustration of the architecture of the ExceedGAN generator and Algorithm 1 for details on the training process associated with the level-varying ExceedGAN approach.

# 4 Numerical experiments

The experiments on simulated and real data are reported in Section 4.1 and Section 4.2 respectively. In both cases, the Adam optimizer [47] is implemented for model training and parameters optimization, this is nowadays a standard method for stochastic gradient descent. It is used with default parameters $\beta_1 = 0.9$ and $\beta_2 = 0.999$ for the exponential decay rates. All experiments have been conducted on an Apple M2 Pro with 16Go Ram memory. The code was implemented in Python 3.8.2 using the library PyTorch 1.7.1 for the GANs' training.

## 4.1 Validation on synthetic data

We consider an experimental design similar to the one of [21, Section 4] where the explored hyperparameter ranges are reported in [21, Table 5], except that the focus is on simulations of upper tail values in dimension $D = 2$, i.e. in the upper quadrant. Denoting by $C$ the copula of $\mathbf{X}$, one has $\mathbb{P}(\mathbf{X} \in Q(\delta_n)) = \hat{C}(\delta_n^{(1)}, \delta_n^{(2)})$, where $\hat{C}$ is the
<figure><img src="image_5.png" /><figcaption>Fig. 2: Generator of the ExceedGAN with $q = 3$ and $D = 2$.</figcaption></figure>

survival copula associated with C, see [48, Section 2.6]. Using the bivariate Gumbel copula:

$$C(u, v) = \exp \left[ - \left\{ (\log 1/u)^\mu + (\log 1/v)^\mu \right\}^{1/\mu} \right], (u, v) \in (0, 1]^2, \mu > 0,$$

we simulate N exceedances in $Q(\delta_n)$ with $\delta_n = (0.1, 0.1)^\top$ and $\delta_n = (0.05, 0.05)^\top$ using a basic acceptance-rejection method. The same values of the dependence parameter $\mu \in \{1.1, 2, 10\}$ are used leading to a Kendall's tau $\tau \in \{0.1, 0.5, 0.9\}$. The two margins are chosen to be Burr distributed, with cumulative distribution function

$$F_m(x) = 1 - (1 + x^{-\rho_2^{(m)}/\gamma})^{1/\rho_2^{(m)}}, \quad x \ge 0, m \in \{1, 2\},$$

common tail-index $\gamma \in \{0.3, 0.5, 0.9\}$ and second-order parameters $(\rho_2^{(1)}, \rho_2^{(2)}) \in \{(-1, -2), (-1, -3), (-2, -3)\}$. Summarizing, one obtains $2 \times 3 \times 3 \times 3 = 54$ simulation settings.

Performance comparisons are made between classical GAN, EV-GAN [21] and the proposed ExceedGAN through the Mean Squared Logarithmic Error (MSLE)

$$\text{MSLE}(\delta_n) = \frac{1}{ND\{\log(2)\}^2} \sum_{m=1}^{D} \sum_{i=1}^{N} \left\{ \log Y_{i,N}^{(m)}(\delta_n^{(m)}) - \log \tilde{Y}_{i,N}^{(m)}(\delta_n^{(m)}) \right\}^2 \quad (24)$$

between the N original exceedances $\{\mathbf{Y}_1(\delta_n), \dots, \mathbf{Y}_N(\delta_n)\}$ and the generated ones $\{\tilde{\mathbf{Y}}_1(\delta_n), \dots, \tilde{\mathbf{Y}}_N(\delta_n)\}$. Note that the MSLE is the usual mean squared error criterion computed on a logarithm scale to deal with heavy-tailed random variables, see also
**Theorem 3.2.** The log(2) factor is introduced so that a 100% relative error on the log-marginals (i.e. $\tilde{Y}_{i,N}^{(m)}(\delta_n^{(m)}) = 2Y_{i,N}^{(m)}(\delta_n^{(m)})$ or $\tilde{Y}_{i,N}^{(m)}(\delta_n^{(m)}) = Y_{i,N}^{(m)}(\delta_n^{(m)})/2$ for all $i \in \{1, \dots, N\}$ and $m \in \{1, \dots, D\}$) yields MSLE($\delta_n$) = 1.

Two experimental frameworks are considered. First, we aim at comparing the performance of the three generative models GAN, EV-GAN and fixed-level ExceedGAN to learn the excess distribution in two extreme regions $Q(\delta_n)$ with $\delta_n = (0.1, 0.1)^\top$ and $\delta_n = (0.05, 0.05)^\top$ using respectively $N = 1000$ and $N = 250$ training points. The models are trained during 10000 iterations and every 10 iterations the MSLE is computed on the training set (thus giving 1000 values). The parameters of each model associated with the best results (i.e. the lowest MSLE) among the 1000 checkpoints are selected. Note that one could choose another performance criteria (such as the Wasserstein distance for instance) which would lead to other selected models. For both regions, the performance of a selected model is assessed by computing the MSLE between the generated samples in the appropriate region through acceptance-rejection sampling (with a common input noise for all competitors) and a testing set of 10000 points. Note that, here, fixed-level ExceedGAN is trained thanks to the same procedure as GAN and EV-GAN with fixed anchor point $u_n$. The objective is to highlight the performance of the eLU parametrization for simulating extreme values.

Second, we focus on the ability of level-varying ExceedGAN described in Algorithm 1 to learn the excess distribution for all $\delta_n \in [a, 1]^2$ with $N = 100000$ training points and $a = 1/2$. Note that the training time is in average 10% longer than the three other models, mainly because of Step 2 in Algorithm 1. All models take approximately one hour to run. The model selection procedure is similar to the previous one, except that we compute on the training set an average of MSLEs on the two regions $Q(\delta_n)$ with $\delta_n = (0.1, 0.1)^\top$ and $\delta_n = (0.05, 0.05)^\top$. The model performance is assessed on the same testing set as in study 1) in order to compare the results between the fixed-level and level-varying versions of ExceedGAN.

Table 1 shows that fixed-level ExceedGAN outperforms both classical GAN and EV-GAN in 44 out of the 54 considered configurations. EV-GAN provides the best result in 9 out of the 10 remaining cases. Fixed-level ExceedGAN is particularly effective for heavy tails ($\gamma \ge 0.5$), yielding the smallest MSLE in 32 out of the such situations. As expected, level-varying ExceedGAN is less suitable than the fixed-level version with a total of 19 best results out of 54, while fixed-level ExceedGAN and EV-GAN provide the best result respectively in the 31 and 4 remaining cases. Note that the level-varying model seems particularly efficient in the extreme region associated with $\delta_n = (0.05, 0.05)^\top$ reaching 13 best results out of its total 19 top scores. In contrast, the classical GAN is not adapted to the simulation of extreme values: It never provides the best result, with a very disappointing behavior for $\gamma = 0.9$. Focusing on the considered parameters ($\gamma, \mu, \rho$), the MLSE criterion associated with ExceedGAN is increasing when the tail gets heavier (large $\gamma$) requiring a more distant extrapolation, decreasing (except for the level-varying version) when the dependence increases (large $\mu$) since the simulated points are less spread out in the space, and mostly increasing with the bias parameters $\rho_1$ and $\rho_2$.
Algorithm 1: Level-varying ExceedGAN training

**Input:** $m_U$, batch size of the conditional variable, $m_X$, batch size of the data for each conditional value $a$, left support of the conditional distribution
**Output:** $(\hat{\theta}, \hat{\phi})$, trained parameters
1 **for** number of iterations **do**
&nbsp;&nbsp;&nbsp;&nbsp;sample a minibatch of $m_U$ levels $\{\delta_k^{(1)} \sim \mathcal{U}(a, 1), \dots, \delta_k^{(D)} \sim \mathcal{U}(a, 1)\}_{k=1}^{m_U}$
&nbsp;&nbsp;&nbsp;&nbsp;and store the associated thresholds $\{u_k^{(1)}, \dots, u_k^{(D)}\}_{k=1}^{m_U}$
&nbsp;&nbsp;&nbsp;&nbsp;**for** $k = 1 : m_U$ **do**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**do**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sample a minibatch of $m_X$ data $\{\mathbf{x}_{i,k} \in Q(\delta_k)\}_{i=1}^{m_X}$ with replacement
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update the discriminator by computing the gradient
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\nabla_{\phi} \frac{1}{m_X m_U} \sum_{k=1}^{m_U} \sum_{i=1}^{m_X} [\log(D_{\phi}^{\text{EX}}(\mathbf{x}_{i,k}, \delta_k)) + \log(1 - D_{\phi}^{\text{EX}}(G_{\theta}^{\text{EX}}(\mathbf{z}_{i,k}, \delta_k), \delta_k))]$,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;with $\mathbf{z}_{i,k} \sim p_{\mathbf{z}}$.
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update the generator by computing the gradient
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\nabla_{\theta} \frac{1}{m_X m_U} \sum_{k=1}^{m_U} \sum_{i=1}^{m_X} [\log(1 - D_{\theta}^{\text{EX}}(G_{\theta}^{\text{EX}}(\mathbf{z}_{i,k}, \delta_k), \delta_k))]$,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;with $\mathbf{z}_{i,k} \sim p_{\mathbf{z}}$.
2 **end for**

Figures 3 and 4 illustrate the solid performance of both ExceedGAN versions for simulating in heavy distribution tails. While EV-GAN manages to fit well the distribution in the intermediate upper quadrant associated with $\delta_n = (0.1, 0.1)^\top$, it underestimates the tails in the extreme one where $\delta_n = (0.05, 0.05)^\top$. Additionally, these graphs highlight the challenges raised in Section 2.4 for a classical GAN to learn heavy distribution tails.

As a conclusion, this numerical study illustrates the dominance of the ExceedGAN in a selection of heavy-tailed cases, with a benefit to use the fixed-level version when enough observations are available in the tails; but also with the level-varying version as a competitive alternative that performs well in all regions.

## 4.2 Illustration on cryptocurrency data

The considered data are the daily log-returns of a portfolio composed by five cryptocurrency pairs (BTC/USD, ETH/USD, SOL/USD, XRP/USD, DOGE/USD) coming from the Kaiko Benchmark Reference Rates<fnref n="1" /> during more than four years between 2021-03-24 and 2025-10-26. Filtering out on the joint negative log-returns and taking their absolute values, one obtains a dataset of $n = 432$ observations. The heavy-tail behavior of the data is illustrated in Figure 5 with estimated tail-indices $\hat{\gamma}_{\text{btc}} \approx \hat{\gamma}_{\text{eth}} \approx$

<footnote n="1">¹https://www.kaiko.com/indices/benchmarks</footnote>
<table><caption>**Table 1:** Comparison between the best GAN (first column), EV-GAN (second column), fixed-level ExceedGAN (third column) and level-varying ExceedGAN (fourth column) results on 10 000 simulated test data for the 27 combinations of parameters using the MSLE criterion in the regions associated with $\delta_n = (0.1, 0.1)^T$ and $\delta_n = (0.05, 0.05)^T$. MSLE greater than 1 are not reported. Values are scaled by $10^2$ and the best results are emphasized in bold.</caption><thead><tr><th rowspan="2">γ</th><th rowspan="2">(ρ<sub>1</sub>, ρ<sub>2</sub>)</th><th colspan="4">μ = 1.1</th><th colspan="4">μ = 2</th><th colspan="4">μ = 10</th></tr><tr><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th></tr></thead><tbody><tr><td rowspan="3">0.3</td><td>(-1, -2)</td><td>0.46</td><td>0.52</td><td><strong>0.30</strong></td><td>0.38</td><td>0.23</td><td>0.26</td><td><strong>0.11</strong></td><td>0.46</td><td>0.34</td><td>0.97</td><td>0.18</td><td><strong>0.16</strong></td></tr><tr><td>(-1, -3)</td><td>0.98</td><td>0.36</td><td><strong>0.35</strong></td><td>3.28</td><td>0.43</td><td>0.24</td><td><strong>0.13</strong></td><td>0.31</td><td>0.44</td><td>0.84</td><td>0.22</td><td><strong>0.20</strong></td></tr><tr><td>(-2, -3)</td><td>0.53</td><td>0.35</td><td>0.39</td><td><strong>0.24</strong></td><td>0.14</td><td>0.24</td><td><strong>0.11</strong></td><td>0.76</td><td>0.29</td><td>1.29</td><td><strong>0.05</strong></td><td>0.12</td></tr><tr><td rowspan="3">0.5</td><td>(-1, -2)</td><td>4.59</td><td>1.49</td><td>1.07</td><td><strong>1.03</strong></td><td>3.18</td><td>1.11</td><td><strong>0.54</strong></td><td>6.72</td><td>2.64</td><td>1.33</td><td><strong>0.39</strong></td><td>1.01</td></tr><tr><td>(-1, -3)</td><td>3.93</td><td>1.21</td><td>1.19</td><td><strong>1.16</strong></td><td>3.98</td><td>1.41</td><td>1.30</td><td><strong>1.19</strong></td><td>2.38</td><td>1.50</td><td><strong>0.22</strong></td><td>2.60</td></tr><tr><td>(-2, -3)</td><td>6.03</td><td><strong>1.09</strong></td><td>1.67</td><td>2.45</td><td>1.22</td><td>3.05</td><td><strong>0.69</strong></td><td>2.53</td><td>1.93</td><td>1.30</td><td><strong>0.51</strong></td><td>2.09</td></tr><tr><td rowspan="3">0.9</td><td>(-1, -2)</td><td>61.92</td><td><strong>6.13</strong></td><td>6.68</td><td>13.54</td><td>35.82</td><td>3.77</td><td><strong>2.11</strong></td><td>15.72</td><td>28.88</td><td>4.83</td><td><strong>1.32</strong></td><td>86.05</td></tr><tr><td>(-1, -3)</td><td>38.65</td><td><strong>4.50</strong></td><td>4.85</td><td>25.65</td><td>44.17</td><td>3.82</td><td><strong>2.59</strong></td><td>13.13</td><td>24.56</td><td>5.69</td><td><strong>1.20</strong></td><td>84.16</td></tr><tr><td>(-2, -3)</td><td>41.19</td><td>5.49</td><td><strong>5.05</strong></td><td>21.92</td><td>45.03</td><td>4.11</td><td><strong>1.86</strong></td><td>15.15</td><td>36.71</td><td>3.73</td><td><strong>1.56</strong></td><td>83.37</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">γ</th><th rowspan="2">(ρ<sub>1</sub>, ρ<sub>2</sub>)</th><th colspan="4">μ = 1.1</th><th colspan="4">μ = 2</th><th colspan="4">μ = 10</th></tr><tr><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th><th>GAN</th><th>EV-GAN</th><th>ExceedGAN</th></tr></thead><tbody><tr><td rowspan="3">0.3</td><td>(-1, -2)</td><td>2.98</td><td>1.38</td><td>1.14</td><td><strong>0.36</strong></td><td>2.10</td><td>0.86</td><td>1.91</td><td><strong>0.36</strong></td><td>1.23</td><td>1.90</td><td>1.22</td><td><strong>0.62</strong></td></tr><tr><td>(-1, -3)</td><td>2.75</td><td>2.76</td><td><strong>0.56</strong></td><td>2.60</td><td>2.35</td><td>1.28</td><td>1.50</td><td><strong>0.89</strong></td><td>1.35</td><td>1.82</td><td>1.27</td><td><strong>0.34</strong></td></tr><tr><td>(-2, -3)</td><td>2.34</td><td>1.12</td><td>2.10</td><td><strong>0.17</strong></td><td>1.88</td><td>1.39</td><td>1.47</td><td><strong>0.46</strong></td><td>0.75</td><td>1.88</td><td>1.19</td><td><strong>0.43</strong></td></tr><tr><td rowspan="3">0.5</td><td>(-1, -2)</td><td>18.58</td><td>6.19</td><td>2.21</td><td><strong>1.18</strong></td><td>13.88</td><td>5.22</td><td><strong>2.21</strong></td><td>6.64</td><td>8.87</td><td>3.76</td><td>2.77</td><td><strong>1.60</strong></td></tr><tr><td>(-1, -3)</td><td>15.96</td><td>4.71</td><td>5.55</td><td><strong>1.57</strong></td><td>15.08</td><td>3.20</td><td><strong>2.17</strong></td><td>2.40</td><td>9.06</td><td>4.37</td><td><strong>2.64</strong></td><td>3.16</td></tr><tr><td>(-2, -3)</td><td>20.91</td><td>7.16</td><td>5.59</td><td><strong>1.11</strong></td><td>14.23</td><td>4.08</td><td><strong>2.16</strong></td><td>3.93</td><td>6.24</td><td>5.04</td><td><strong>2.18</strong></td><td>5.28</td></tr><tr><td rowspan="3">0.9</td><td>(-1, -2)</td><td>—</td><td>16.38</td><td>13.70</td><td><strong>12.97</strong></td><td>—</td><td>9.27</td><td><strong>7.81</strong></td><td>15.72</td><td>90.41</td><td>15.66</td><td><strong>7.10</strong></td><td>88.78</td></tr><tr><td>(-1, -3)</td><td>—</td><td>22.34</td><td>17.54</td><td><strong>14.21</strong></td><td>—</td><td>12.06</td><td><strong>6.82</strong></td><td>12.91</td><td>—</td><td>15.62</td><td><strong>7.37</strong></td><td>86.84</td></tr><tr><td>(-2, -3)</td><td>—</td><td>20.21</td><td><strong>14.31</strong></td><td>18.13</td><td>—</td><td>11.90</td><td><strong>9.96</strong></td><td>20.70</td><td>—</td><td>34.13</td><td><strong>7.46</strong></td><td>82.66</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 3:** Example of results on simulated bivariate data from Burr margins with $\gamma \in \{0.3, 0.5, 0.9\}$ (from left to right), $(\rho_1, \rho_2) = (-1, -3)$ and Gumbel copula with $\mu = 2$. Comparison between the test set generated from the real distribution (black crosses) and the best model according to the MSLE criterion in the upper quadrant $Q(\delta_n)$ with $\delta_n = (0.1, 0.1)^\top$: original GAN (first row, blue), EV-GAN (second row, orange), fixed-level ExceedGAN (third row, green) and level-varying ExceedGAN (fourth row, red). The axes have been transformed to the log-scale.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4:** Example of results on simulated bivariate data from Burr margins with $\gamma \in \{0.3, 0.5, 0.9\}$ (from left to right), $(\rho_1, \rho_2) = (-1, -3)$ and Gumbel copula with $\mu = 2$. Comparison between the test set generated from the real distribution (black crosses) and the best model according to the MSLE criterion in the upper quadrant $Q(\delta_n)$ with $\delta_n = (0.05, 0.05)^\top$: original GAN (first row, blue), EV-GAN (second row, orange), fixed-level ExceedGAN (third row, green) and level-varying ExceedGAN (fourth row, red). The axes have been transformed to the log-scale.</figcaption></figure>
$\hat{\gamma}_{\text{xrp}} \approx 0.34$, $\hat{\gamma}_{\text{sol}} \approx 0.40$, $\hat{\gamma}_{\text{doge}} \approx 0.29$ which is coherent with previous statistical analysis on cryptocurrency data [49]. Consider the testing set as the $N = 32$ exceedances located in the upper orthant $Q(\boldsymbol{\delta}_n)$ at the level $\boldsymbol{\delta}_n = (0.8, 0.8, 0.8, 0.8, 0.8)^{\top}$, and consider as the training set the lower orthant (i.e. the set where all the margins are lower than the thresholds) containing 268 observations. The challenge is to sample observations in an extreme region where no training point belongs to.

We aim at comparing the performances of the level-varying ExceedGAN and a multivariate Gumbel copula with Pareto margins where the latter are fitted by maximum likelihood estimation and where the dependence structure is parametrized by the estimated dependence coefficient $\hat{\mu} = 1.43$ (leading to $\hat{\tau} = 0.30$) on the training set. Once both models are trained, we sample from each of them 10 000 replications of $N$ points in $Q(\boldsymbol{\delta}_n)$ and estimate the MSLE (24) with respect to the testing set. While the ExceedGAN will naturally sample in $Q(\boldsymbol{\delta}_n)$, a acceptance-rejection method is used for the copula model.

The MSLE results of all replications are illustrated in Figure 6 with respect to the number of samples required to obtain $N$ points in the extreme region $Q(\boldsymbol{\delta}_n)$. It appears that although the Copula model has a slightly smaller error than the ExceedGAN in average (MSLE of 8.5% vs 12%), it requires to sample in average 50 times more points. Therefore the ExceedGAN model represents a better cost-effective model to sample in extreme regions.

Finally, note that this experiment is conducted under the assumption that the portfolio return distribution is stationary over the 2021-2025 period. We refer to [50] for a method to capture the temporal dynamics in the tail dependence structure.

# 5 Conclusion

We have introduced a new generative model dedicated to extremes which is able to learn the distribution of the excess. A bias reduction technique is derived from extreme value theory to overcome the difficulties inherent in extrapolation using generative models. Using an appropriate eLU basis functions, we recover the uniform convergence rate on the margins derived in [22] for neural network extreme quantile approximation and we illustrate the solid performance of the proposed model against another state of the art competitor on numerous numerical experiments.

It would be interesting to numerically compare ExceedGAN with the method recently introduced in [51]. It relies on a combination of univariate extreme value analysis to model the marginal distribution tails with a GAN modeling of the angular distribution.

Extensions to other risk measure estimation can be considered using such a data augmentation technique with light numerical resources. Such a generative model indeed allows one to sample in extreme regions and then estimate risk measures by their empirical counterpart without recourse to extrapolation techniques. Additionally, it would be interesting to extend this work to address the extreme dependence structure as proposed for instance in [36, 52].
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 5:** Illustration on real data. Log quantile-quantile plot $\log((n + 1)/i)$ $\mapsto \log X_{n-i+1,n}^{(m)}$, $i \in \{1, \dots, \lfloor n/10 \rfloor\}$, for each marginal $m \in \{1, \dots, 5\}$ associated respectively with (a): BTC/USD, (b): ETH/USD, (c): XRP/USD, (d): SOL/USD, (e): DOGE/USD. The regression lines are superimposed.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 6:** Scatter plot on 10,000 replications of MSLEs with respect to the number of samples required to obtain *N* points in the testing set *Q*(δ*n*) with δ*n* = (0.8, 0.8, 0.8, 0.8, 0.8)⊤ from both a level-varying ExceedGAN (blue) and the Gumbel copula (orange). The red triangles represent the empirical mean for both models.</figcaption></figure>
References

[1] Asmussen, S., Albrecher, H.: Ruin Probabilities, 2nd edn. Advanced Series on Statistical Science & Applied Probability, 14. World Scientific Publishing, Hackensack, NJ (2010)

[2] Embrechts, P., Klüppelberg, C., Mikosch, T.: Modelling Extremal Events: for Insurance and Finance vol. 33. Springer, Berlin (2013)

[3] Prandini, M., Watkins, O.J.: Probabilistic aircraft conflict detection. HYBRIDGE WP3: Reachability analysis for probabilistic hybrid systems (2005)

[4] Cirillo, P., Taleb, N.N.: Tail risk of contagious diseases. Nature Physics 16(6), 606–613 (2020)

[5] Coles, S.: An Introduction to Statistical Modeling of Extreme Values. Springer, London (2001)

[6] de Haan, L., Ferreira, A.: Extreme Value Theory. Springer Series in Operations Research and Financial Engineering. Springer, New York (2006)

[7] Weissman, I.: Estimation of parameters and large quantiles based on the $k$ largest observations. Journal of the American Statistical Association 73(364), 812–815 (1978)

[8] Girard, S., Stupfler, G., Usseglio-Carleve, A.: Extreme conditional expectile estimation in heavy-tailed heteroscedastic regression models. The Annals of Statistics 49(6), 3358–3382 (2021)

[9] Daouia, A., Gijbels, I., Stupfler, G.: Extremiles: A new perspective on asymmetric least squares. Journal of the American Statistical Association 114(527), 1366–1381 (2019)

[10] He, Y., Einmahl, J.H.J.: Estimation of extreme depth-based quantile regions. Journal of the Royal Statistical Society Series B 79(2), 449–461 (2016)

[11] Bucklew, J.A.: Introduction to Rare Event Simulation. Springer Series in Statistics. Springer, New York (2004)

[12] Gobet, E., Liu, G.: Rare event simulation using reversible shaking transformations. SIAM Journal on Scientific Computing 37(5), 2295–2316 (2015)

[13] Del Moral, P., Garnier, J.: Genealogical particle analysis of rare events. Annals of Applied Probability 15(4), 2496–2534 (2005)

[14] Foster, D.: Generative Deep Learning: Teaching Machines to Paint, Write, Compose, and Play. O'Reilly Media, Beijing-Boston-Farnham-Sebastopol-Tokyo (2019)
[15] Kingma, D.P., Welling, M.: Auto-encoding variational Bayes. In: 2nd International Conference on Learning Representation (Banff) (2014)

[16] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. In: Advances in Neural Information Processing Systems (Montreal), vol. 27, pp. 2672–2680 (2014)

[17] Bishop, C.M., Bishop, H.: Deep Learning: Foundations and Concepts. Springer, New York (2023)

[18] Rezende, D., Mohamed, S.: Variational inference with normalizing flows. In: Proceedings of the 32nd International Conference on Machine Learning (Lille), vol. 37, pp. 1530–1538 (2015). PMLR

[19] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. In: Advances in Neural Information Processing Systems, vol. 33, pp. 6840–6851 (2020)

[20] Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., Ganguli, S.: Deep unsupervised learning using nonequilibrium thermodynamics. In: Proceedings of the 32nd International Conference on Machine Learning (Lille), vol. 37, pp. 2256–2265 (2015). PMLR

[21] Allouche, M., Girard, S., Gobet, E.: EV-GAN: Simulation of extreme events with ReLU neural networks. Journal of Machine Learning Research 23(150), 1–39 (2022)

[22] Allouche, M., Girard, S., Gobet, E.: Estimation of extreme quantiles from heavy-tailed distributions with neural networks. Statistics and Computing 34(12) (2024)

[23] Allouche, M., Girard, S., Gobet, E.: On the simulation of extreme events with neural networks. In: de Carvalho, M., Huser, R., Naveau, P., Reich, B.J. (eds.) Handbook on Statistics of Extremes. Chapman & Hall/CRC, Boca Raton, FL (2026)

[24] Bertsekas, D.P., Shreve, S.E.: Stochastic Optimal Control. Mathematics in Science and Engineering, vol. 139. Academic Press, Inc., New York-London (1978)

[25] Villani, C.: Optimal Transport: Old and New. Springer, Berlin (2009)

[26] Peyré, G., Cuturi, M.: Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning 11(5-6), 355–607 (2019)

[27] Pinkus, A.: Approximation theory of the MLP model in neural networks. In: Acta Numerica vol. 8, pp. 143–195. Cambridge University Press, Cambridge (1999)

[28] Biau, G., Cadre, B., Sangnier, M., Tanielian, U.: Some theoretical properties of
GANs. The Annals of Statistics 48(3), 1539-1566 (2020)

[29] Eckhardt, R.: Stam Ulam, John Von Neumann and the Monte-Carlo method. Los Alamos Science **Special Issue**, 131–143 (1987)

[30] Chamakh, L., Gobet, E., Szabo, Z.: Orlicz random Fourier feature. Journal of Machine Learning Research **21**, 1–37 (2020)

[31] Bhatia, S., Jain, A., Hooi, B.: Exgan: Adversarial generation of extreme samples. In: Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, pp. 6750–6758 (2021)

[32] Boulaguiem, Y., Zscheischler, J., Vignotto, E., Wiel, K., Engelke, S.: Modeling and simulating spatial extremes by combining extreme value theory with generative adversarial networks. Environmental Data Science **1**, 5 (2022)

[33] Wiese, M., Knobloch, R., Korn, R., Kretschmer, P.: Quant GANs: deep generation of financial time series. Quantitative Finance **20**(9), 1419–1440 (2020)

[34] Feder, R.M., Berger, P., Stein, G.: Nonlinear 3D cosmic web simulation with heavy-tailed generative adversarial networks. Physical Review D. **102**(10), 103504–18 (2020)

[35] Huster, T., Cohen, J.E.J., Lin, Z., Chan, K., Kamhoua, C., Leslie, N., Chiang, C.Y.J., Sekar, V.: Pareto GAN: Extending the representational power of GANs to heavy-tailed distributions. In: International Conference on Machine Learning, pp. 4523–4532 (2021). PMLR

[36] Lafon, N., Naveau, P., Fablet, R.: A VAE approach to sample multivariate extremes. arXiv preprint arXiv:2306.10987 (2023)

[37] Zhang, L., Ma, X., Wikle, C.K., Huser, R.: Flexible and efficient spatial extremes emulation via variational autoencoders. arXiv preprint arXiv:2307.08079 (2024)

[38] Hickling, T., Prangle, D.: Flexible tails for normalising flows, with application to the modelling of financial return data. In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 374–386 (2023)

[39] Hu, C., Castro-Camilo, D.: GPDFlow: Generative multivariate threshold exceedance modeling via normalizing flows. arXiv preprint arXiv:2503.11822 (2025)

[40] De Monte, L., Huser, R., Papastathopoulos, I., Richards, J.: Generative modelling of multivariate geometric extremes using normalising flows. arXiv preprint arXiv: 2505.02957 (2025)

[41] Wessel, J.B., Murphy-Barltrop, C.J., Simpson, E.S.: A comparison of generative deep learning methods for multivariate angular simulation. arXiv preprint arXiv:
2504.21505 (2025)

[42] Wang, X.Q., Cheng, S.H.: General regular variation of the n-th order and 2nd order Edgeworth expansions of the extreme value distribution. II. Acta Mathematica Sinica (Engl. Ser.) **22**(1), 27–40 (2006)

[43] Fraga Alves, I., Haan, L., Lin, T.: Third order extended regular variation. Publications de l'Institut Mathématique **80**(94), 109–120 (2006)

[44] Neves, C.: From extended regular variation to regular variation with application in extreme value statistics. Journal of Mathematical Analysis and Applications **355**(1), 216–230 (2009)

[45] Goegebeur, Y., Wet, T.: Estimation of the third-order parameter in extreme value statistics. Test **21**(2), 330–354 (2012)

[46] Mirza, M., Osindero, S.: Conditional generative adversarial nets. arXiv preprint arXiv:1411.1784 (2014)

[47] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. In: International Conference on Learning Representation (San Diego) (2015)

[48] Nelsen, R.: An Introduction to Copulas, 2nd edn. Springer Series in Statistics. Springer, New York (2006)

[49] Allouche, M., Echenim, M., Gobet, E., Maurice, A.-C.: Statistical error bounds for weighted mean and median, with application to robust aggregation of cryptocurrency data. Mathematical Finance **35**(4), 760–778 (2025)

[50] de Carvalho, M., Palacios Ramirez, K.V.: Semiparametric Bayesian modelling of nonstationary joint extremes: How do big tech's extreme losses behave? Journal of the Royal Statistical Society Series C **74**(2), 447–465 (2025)

[51] Lhaut, S., Rootzén, H., Segers, J.: Wasserstein-Aitchison GAN for angular measures of multivariate extremes. arXiv preprint arXiv:2504.21438 (2025)

[52] Girard, S., Gobet, E., Pachebat, J.: HTGAN: Heavy-tail GAN for multivariate dependent extremes via latent-dimensional control. International Journal of Computer Mathematics (2025). To appear.

**Acknowledgments.** This work is partially supported by the French National Research Agency (ANR) under the grant ANR-23-CE40-0009 and the grant ANR "BLOCKFI: Blockchain and Decentralized Finance" - 2024 - CES38. S. Girard and E. Gobet acknowledge the support of the Chair Stress Test, Risk Management and Financial Steering, led by the French Ecole Polytechnique and its Foundation and sponsored by BNP Paribas.