# Assessing Extreme Risk using Stochastic Simulation of Extremes

Nisrine Madhar1, Juliette Legrand2 and Maud Thomas3

###### Abstract

Risk management is particularly concerned with extreme events, but analysing
these events is often hindered by the scarcity of data, especially in a
multivariate context. This data scarcity complicates risk management efforts.
Various tools can assess the risk posed by extreme events, even under
extraordinary circumstances. This paper studies the evaluation of univariate
risk for a given risk factor using metrics that account for its asymptotic
dependence on other risk factors. Data availability is crucial, particularly
for extreme events where it is often limited by the nature of the phenomenon
itself, making estimation challenging. To address this issue, two non-
parametric simulation algorithms based on multivariate extreme theory are
developed. These algorithms aim to extend a sample of extremes jointly and
conditionally for asymptotically dependent variables using stochastic
simulation and multivariate Generalised Pareto Distributions. The approach is
illustrated with numerical analyses of both simulated and real data to assess
the accuracy of extreme risk metric estimations.

1 Université Paris Cité, CNRS, Laboratoire de Probabilités, Statistique et
Modélisation, LPSM, F-75013 Paris, France,
2 Univ Brest, CNRS UMR 6205, Laboratoire de Mathématiques de Bretagne
Atlantique, France
3 Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et
Modélisation, LPSM, 4 place Jussieu, F-75005 Paris, France,

E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Keywords – Multivariate Generalised Pareto distributions, Risk management,
Simulation of multivariate extremes, Tail risk metrics

##  1 Introduction

Risk management is of crucial importance in various sectors. It involves the
identification, assessment, monitoring and mitigation of potential risks,
regardless of the area of application. Risk is a common factor to various
domains, yet its manifestation is specific to each sector. In climatology,
meteorological and marine hazards such as droughts, floods and landslides can
cause significant material damages and endanger the affected areas. The
identification and anticipation of these events is thus of paramount
importance in order to prevent them. In finance, market movements can result
to substantial losses. It is thus essential to anticipate and mitigate these
risks, which can be achieved through in-depth analysis of market trends and
the implementation of tailored risk management strategies. In the sequel, any
variable that could result in a loss or damage is referred to as a risk
factor. For instance, in climatology, a risk factor could be any physical
quantity, such as wave heights (Legrand et al., 2023), wind gusts (Goegebeur
et al., 2024) or precipitation levels (Gründemann et al., 2023). In finance,
risk factors are typically market parameters, such as interest rates or
exchange rates, which may induce potential losses for financial institutions
if they experience unfavourable fluctuations (see e.g. McNeil et al., 2015).
The risk factors can be represented by random variables that quantify the
magnitude of potential losses. It is typically the case that the magnitudes
with the strongest impact result from events that occur with a very low
probability of occurrence. This type of event is often referred to as a tail
risk (McNeil et al., 2015). A primary objective of risk managers is to
quantify, for a given target risk factor, the associated tail risk.

Modern risk management has witnessed a surge in the development of tail risk
measures and methodologies for their estimation. The simplest and most common
measure often considered for evaluating tail risks is the Value at Risk (VaR),
which represents the quantile at a level $\alpha\in(0,1)$ of a given random
variable. Nevertheless, Artzner et al. (1999) have identified significant
shortcomings of the $\mathrm{VaR}$ approach, indicating that tail risk
assessment through $\mathrm{VaR}$ estimation could result in an inaccurate
quantification of the risk. Several limitations of the $\mathrm{VaR}$ have
been identified in the literature. In particular, this risk metrics does not
account for the severity of the target risk factor, which may be a loss in
financial applications (Yamai and Yoshiba, 2005), or the magnitude of a
climatic variable in environmental applications, e.g. the magnitude of daily
precipitation (Gründemann et al., 2023). In the literature, a number of
alternative approaches have been proposed (Bellini and Di Bernardino, 2017;
Daouia et al., 2018; Tasche, 2002). In this paper, we focus on three risk
metrics (introduced in Section 2), referred to as tail risk metrics (TRMs),
which are defined conditionally on the tail distribution of a given random
variable, or target risk factor.

Extreme Value Theory (EVT) provides the theoretical foundations for the
analysis and modelling of tail distributions, making it an ideal framework for
the estimation of TRMs (see e.g. Resnick, 2007). In this paper, we propose two
non-parametric simulation schemes for multivariate extremes by extending the
procedures developed by Legrand et al. (2023) in the bivariate case to higher
dimensions. The first algorithm—joint simulation of multivariate
extremes—allows for the estimation of TRMs and the second one—conditional
simulation of multivariate extremes—enables the inference of quantities
involving some conditional tail distribution.

Several parametric approaches based on EVT have been proposed for the
estimation of TRMs (see e.g. McNeil et al., 2015; Singh et al., 2013).
However, the parametric setting is quite restrictive in the multivariate
framework. Indeed, in contrast to the univariate framework, the limiting
distribution of multivariate extremes is no longer parametric (Rootzén and
Tajvidi, 2006). To address this issue, Rootzén et al. (2018) have introduced
several parametrisations and stochastic representations that enable the
calibration of parametric models whose densities could be expressed through
explicit formulas. Given the diversity of potential models, which may include
nested models, model selection is necessary and can prove to be challenging
and time-consuming. An alternative approach would be to model the univariate
margins using univariate extreme value theory, and then to model the
dependence structure of extremes through the typical tools for dependence
modelling, such as copulas. For extreme modelling, Clayton and Gumbel copula
allow for the modelling of lower and upper tail dependence, respectively and
are characterised by a single parameter (see e.g. Nelsen, 2006). This implies
uniform dependence across the considered risk factors, which may not be a
valid assumption. Another limitation of empirical estimates is that for high
risk scenarios, the number of available tail observations becomes exceedingly
limited.

In this paper, we thus rely on a non-parametric approach to expand the number
of observations above the extreme level of interest based on joint simulation
of multivariate extremes, by generalising the approach of Legrand et al.
(2023). Consequently, in addition to ensuring more reliable estimations, it is
also possible to extrapolate beyond the range of observed data, which is not
feasible when considering only historical data, and to go beyond the bivariate
case, which is typically discussed in the literature.

A by-product of this first procedure concerns the estimation of quantities
describing the tail of some conditional distribution, which could be analogous
to an extremal regression approach. Existing approaches that address the issue
of linear regression in the context of extreme values include quantile
regression (see e.g. Pasche and Engelke, 2022; Gnecco et al., 2022, for recent
developments). However, the performance of these approaches is highly
dependent on the available data. It is crucial to acknowledge that the
inherent limitations of the existing methodologies can be attributed to two
key factors: the dimensionality of the random vector of risk factors and the
number of available observations of this random vector. This latter issue is,
in fact, the very same problem that was encountered for TRMs estimation.
However, the joint simulation algorithm does not apply in this context. Once
more, we build upon the work of Legrand et al. (2023), who introduced a second
algorithm for conditional simulation of bivariate extremes. In this paper, we
extend their approach to higher dimensions. This second non-parametric
algorithm enables the generation of new samples of a given random variable
conditionally on joint event of simultaneous extreme risk factors, thus
facilitating the estimation of quantities involving the conditional
distribution of a target risk factor given the observation of other risk
factors with which the target risk factor has a strong tail dependence.

The outline of this paper is as follows. After a formal presentation of the
TRMs of interest in Section 2, the joint and conditional simulation algorithms
for multivariate extremes are developed in Section 3. Their performances are
evaluated first on synthetic data in Section 4 and then on real data in
Section 5.

#### Notations.

Throughout the paper, we use the following notations for vectors. Symbols in
bold denote $d$-dimensional vectors. For example, a $d$-dimension random
vector is denoted by $\bm{a}=\left(a_{1},\ldots,a_{d}\right)$, and
$\bm{a}_{-j}$ denotes the vector $\bm{a}$ deprived from its $j$-th component
for $j=1,\ldots,d$. Operations and relations are meant component-wise, that
is, for example, $\bm{a}\geq\bm{b}=\left(a_{1}\geq b_{1},\ldots,a_{d}\geq
b_{d}\right)$.

##  2 Tail Related Risk metrics

In this section, we formally define the three TRMs considered in this paper.
Hereinafter, $\bm{X}=(X_{1},\ldots,X_{d})$ denotes a $d$-dimensional vector
with marginal continuous probability density function $f_{j}$, for
$j=1,\dots,d$.

As all TRMs are based on the $\mathrm{VaR}$, it is first necessary to recall
its definition. The $\mathrm{VaR}$ at level $\alpha\in(0,1)$ of $X_{j}$,
denoted $\textrm{VaR}_{\alpha}(X_{j})$, is defined as the
$(1-\alpha)$-quantile of the distribution of $X_{j}$, that is

$\mathrm{VaR}_{\alpha}(X_{j})=\inf\\{x\in\mathbb{R}:\mathbb{P}\left(X_{j}>x\right)\leq\alpha\\}\,$

with $\alpha=\alpha_{n}\xrightarrow[n\to\infty]{}0$. The
$\mathrm{VaR}_{\alpha}$ represents the quantity that will be exceeded with
probability $\alpha$. In other words, it corresponds to an extreme quantile of
$X_{j}$. It should be noted that the term “Value at Risk” is primarily used in
the field of financial risk management (e.g. McNeil et al., 2015); in other
domains, it is often referred to as the “return level” (for more details, see
Coles, 2001).

Several limitations of the $\mathrm{VaR}$ have been identified in the
literature. In particular, this TRMs does not account for the severity of the
random variable $X_{j}$, which may be a loss in financial applications (Yamai
and Yoshiba, 2005), or the magnitude of a climatic variable in environmental
applications, e.g. the magnitude of daily precipitation (Gründemann et al.,
2023). To address this shortcoming, alternative TRMs have been proposed,
including the expected shortfall ($\mathrm{ES}$) (Artzner et al., 1999). The
$\mathrm{ES}$ at level $\alpha\in(0,1)$ is defined as

$\mathrm{ES}_{\alpha}(X_{j})=\mathbb{E}\left[X_{j}\mid
X_{j}>\mathrm{VaR}_{\alpha}(X_{j})\right]=\frac{1}{\alpha}\int_{\mathrm{VaR}_{\alpha}(X_{j})}^{\infty}xf_{j}(x)\mathop{}\\!\mathrm{d}x\,.$
(1)

In words, the $\mathrm{ES}_{\alpha}$ corresponds to the expected loss above
the $\mathrm{VaR}_{\alpha}$ (for more details, see e.g. Artzner et al., 1999;
Rockafellar et al., 2000). The $\mathrm{ES}$ is a useful tool for gaining
insights into the severity of losses above a high quantile (i.e. the
$\mathrm{VaR}$). It is therefore commonly used in financial institutions for
the calculation of the minimum capital requirements for market risk as
specified by financial regulators (see e.g. Basel Committee on Banking
Supervision, 2013). Given its high profile, numerous techniques proposed for
$\mathrm{ES}$ estimation have been proposed (see e.g. Nadarajah et al., 2014,
for a comprehensive review). In a parametric modelling approach, $\mathrm{ES}$
is computed as the empirical mean value of observations above the
$\mathrm{VaR}$ which is estimated as the quantile of an extreme value
distribution (see e.g. Coles, 2001). Subsequently, more advanced estimation
approaches were suggested (see e.g. McNeil and Frey, 2000; Singh et al.,
2013). One limitation of $\mathrm{ES}$ lies in its univariate risk assessment
that ignores the potential asymptotic dependence that $X_{j}$ could exhibit
with other risk factors. To address this issue, we consider two alternative
risk metrics accounting for the asymptotic dependence in their univariate risk
evaluation of a target risk factor $X_{j}$.

We thus introduce a second TRM, namely the marginal expected shortfall (MES)
that has been intensively studied (see e.g. Cai et al., 2015; Acharya et al.,
2017). Cai et al. (2015) define the MES as the following conditional
expectation $\mathbb{E}\left[X_{j}\mid S(\bm{X})\geq v\right]$, where $S$ is a
given statistic of the random vector $\bm{X}$ (e.g. the sum, the minimum, the
maximum of $\bm{X}$) and $v\in\mathbb{R}$ a threshold defining the occurrence
of an extreme event. Building on the aforementioned definition, we propose to
study a novel version by considering a conditioning on a joint event of
simultaneous extreme risk factors deprived from the target risk factor
$X_{j}$. That is, we define the multivariate marginal expected shortfall
($\mathrm{MMES}$) of $X_{j}$ at level $\alpha$ as follows

$\mathrm{MMES}_{\alpha}(X_{j};\bm{X})=\mathbb{E}\left[X_{j}\mid\bm{X}_{-j}\geq\bm{v}^{\alpha}_{-j}\right]=\int_{\mathbb{R}}x_{j}f_{X_{j}\mid\bm{X}_{-j}\geq\bm{v}^{\alpha}_{-j}}(x_{j})\mathop{}\\!\mathrm{d}x_{j},$
(2)

where $\bm{v}^{\alpha}=\mathrm{VaR}_{\alpha}(\bm{X})\in\mathbb{R}^{d}$ and
$f_{X_{j}\mid\bm{X}_{-j}\geq\bm{v}^{\alpha}_{-j}}$ denotes the conditional
density of $X_{j}$ given the event $\bm{X}_{-j}\geq\bm{v}^{\alpha}_{-j}$.
Through this formulation the aim is to capture the behaviour of $X_{j}$ when
similar risk factors reach extreme levels.

In a bivariate setting, authors have investigated the incorporation of a
dependent risk factor in order to improve the quantification of tail risk
associated with a target risk factor (see e.g. Josaphat and Syuhada, 2021;
Goegebeur et al., 2024). In this paper, we adopt the definition of dependent
conditional tail expectation ($\mathrm{DCTE}$), as defined in the bivariate
case by Goegebeur et al. (2024). We extend this definition to the general
multivariate case as follows

$\mathrm{DCTE}_{\alpha}(X_{j};\bm{X})=\mathbb{E}\left[X_{j}|\bm{X}\geq\bm{v}^{\alpha}\right]=\int_{v^{\alpha}_{j}}^{\infty}x_{j}f_{X_{j}\mid\bm{X}\geq\bm{v}^{\alpha}}(x_{j})\mathop{}\\!\mathrm{d}x_{j}.$
(3)

This third TRMs quantifies the risk of a target risk factor when similar risk
factors and the target risk itself reach extreme levels simultaneously. This
metric could be relevant in a wide range of applications, including wind gust
analysis (Goegebeur et al., 2024) and car insurance policies (Syuhada et al.,
2022).

Josaphat and Syuhada (2021) have proposed a parametric estimator of
$\mathrm{DCTE}$ in the bivariate setting based on a model combining a Pareto
distribution and a specific copula. However, the conditioning event considered
by Josaphat and Syuhada (2021) constrains the fluctuations of $X_{j}$ to a
specific interval, which hinders the extrapolation of the estimates beyond the
observed data range. Syuhada et al. (2022) have showed that the non-parametric
estimators were more accurate than parametric ones at estimating DCTE.
Finally, Goegebeur et al. (2024) have introduced an estimator using EVT
arguments through a two-step approach, where the $\mathrm{DCTE}$ is first
estimated at some intermediate level in order to then extrapolate the
estimation at extreme levels. However, this approach is limited to the
bivariate case.

In this paper, we propose two non-parametric simulation schemes for
multivariate extremes. The first scheme is for the joint simulation of
multivariate extreme events, which will allow for the empirical estimation of
our three TRMs of interest (Equations (1), (2) and (3)). The second scheme is
for the conditional simulation of multivariate extreme events, which will be
applied to the inference of quantities involving some conditional tail
distribution.

##  3 Extreme Scenario Simulation

This section outlines the two simulation algorithms developed in this study.
These algorithms are natural extensions of the algorithms presented in
(Legrand et al., 2023), where the analysis is restricted to the bivariate
setting.

The simulation procedure presented in Legrand et al. (2023) is based on the
stochastic representation of a standard multivariate generalised Pareto (MGP)
vector proposed by Rootzén et al. (2018). Building on the stochastic
representation presented in Section 3.1, we propose two simulation algorithms.
A joint simulation algorithm generating simulated samples of a MGP vector,
while a conditional simulation algorithm is presented for simulating samples
of a component of the MGP vector conditionally on the other components. Both
algorithms are illustrated through numerical applications on synthetic data in
Section 4, and on real data in Section 5.

###  3.1 Multivariate Generalised Pareto Vectors

Let $\bm{X}$ be a $d$-dimensional random vector. For a vector of suitably
chosen thresholds $\bm{u}\in\mathbb{R}^{d}$, we define the vector of excesses
above $\bm{u}$ as $\bm{X}-\bm{u}$ given that $\bm{X}\not\leq\bm{u}$, where
$\bm{X}\not\leq\bm{0}$ means that at least one of the components of $\bm{X}$
is positive. If $\bm{X}\not\leq\bm{u}$, $\bm{X}$ is said to be extreme.
Rootzén and Tajvidi (2006) have shown that, under certain conditions, the
conditional excesses $\bm{X}-\bm{u}\mid\bm{X}\not\leq\bm{u}$ converge
asymptotically to a non-degenerate distribution $H$ that belongs to the family
of MGP distributions as $\bm{u}\to\bm{\infty}$. As in the univariate EVT (e.g.
Coles, 2001), for $\bm{u}$ large enough, we may approximate the conditional
excesses, namely $\bm{X}-\bm{u}\mid\bm{X}\not\leq\bm{u}$, by a random vector
$\bm{Z}$ which is distributed as a MGP distribution, $H$.

For such vector $\bm{Z}$, the marginals $H_{j}$, $j=1,\ldots,d$, are in
general not univariate generalised Pareto (GP) distributions. However, their
restrictions to the positive subset are GP distributed (see e.g. Rootzén et
al., 2018). That is, for all $j=1,\ldots,d$,

$H_{j}^{+}(x)=\mathbb{P}\left[Z_{j}>x\mid
Z_{j}>0\right]=\left(1+\gamma_{j}x/\sigma_{j}\right)_{+}^{-1/\gamma_{j}},\mbox{
for }x\geq 0\mbox{ such that }\sigma_{j}+\gamma_{j}x>0,$ (4)

where $a_{+}=\max(a,0)$, $\sigma_{j}>0$ is the j-th marginal scale parameter
and $\gamma_{j}\in\mathbb{R}$ is the j-th marginal shape parameter. If
$\gamma_{j}=0$, we use the classical convention by taking the limit as
$\gamma_{j}\to 0$ in Equation (4).

In the sequel, we denote $\bm{\sigma}=(\sigma_{1},\ldots,\sigma_{d})$ and
$\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{d})$ the vectors of marginal scale and
shape parameters of a MGP vector. When $\bm{\sigma}=\bm{1}$ and
$\bm{\gamma}=\bm{0}$, we say that $\bm{Z}$ is distributed according to a
standard MGP distribution. Conversely, a standard MGP vector can be
transformed into a general MGP vector using the following simple
transformation: $Y$ is a MGP vector with parameters
$(\bm{\sigma},\bm{\gamma})$ if and only if
$Y=\bm{\sigma}(e^{\bm{\gamma}\bm{Z}}-1)/\bm{\gamma}$, with $\bm{Z}$ a standard
MGP vector. In light of this result, we can therefore concentrate solely on
the special case of standard MGP vectors.

Note that any vector $\bm{X}$ can be transformed into a standard MGP vector by
following the two steps outlined below:

  1. 1.

Standardise the data $\bm{X}$ to exponential margins $\bm{X}^{E}$, using the
probability integral transform on each component

$X^{E}_{j}=\log\left(1-F_{j}(X_{j})\right),\quad\mbox{for }j=1,\ldots,d,$

where $F_{j}$ is the cumulative distribution function (c.d.f.) associated with
the $j$-th risk factor for $j=1,\ldots,d$. $F_{j}$ can be either a parametric
c.d.f. that fits the data or the empirical version of c.d.f..

  2. 2.

Compute the multivariate excess vector as defined in Rootzén and Tajvidi
(2006) as follows

$\bm{Z}=\bm{X}^{E}-\bm{u}^{E}\mid\bm{X}^{E}\not\leq\bm{u}^{E}\,,$ (5)

where $\bm{u}^{E}$ is a suitably chosen threshold on the exponential scale.

The choice of the thresholds $u_{j}$, for $j=1,\ldots,d$ can be understood as
a compromise between bias and variance. The smaller the threshold, the less
valid the asymptotic approximation, which leads to bias. On the other hand, a
threshold too high will generate few excesses to fit the model, which leads to
high variance. In practice, threshold selection is a challenging task. The
existing methods for the choice of the threshold $u$ rely on graphical
diagnostics or on computational approaches based on supplementary conditions
(that depend on unknown parameters) on the underlying distribution function
$F$ (see Scarrott and MacDonald, 2012).

Rootzén et al. (2018) have derived different stochastic representations of MGP
distributions for which explicit density formulas can be obtained (see also
Kiriliouk et al., 2019). In particular, they have shown that a standard MGP
vector $\bm{Z}$ can be decomposed as follows

$\bm{Z}=E+\bm{T}-\max\left(\bm{T}\right),$ (6)

where $E$ is a unit exponential variable and $\bm{T}$ a $d$-dimensional random
vector independent of $E$.

In Equation (6), it is necessary to ensure that the components of the original
vector $\bm{X}$ are asymptotically dependent (see e.g. Coles, 2001), simply
meaning that the large values of the components $X_{j}$ occur simultaneously.
In the applications that we have in mind, this hypothesis will be verified
since the risk factors that will be considered are precisely those exhibiting
asymptotic dependence.

From Equation (6), we derive the cornerstone of our two simulation algorithms
by simply considering the following multivariate difference

$\Delta^{j,k}=Z_{j}-Z_{k}=T_{j}-T_{k},\mbox{ for all }j,k=1,\ldots,d.$ (7)

Subsequently, Equation (6) can be rewritten as follows

$Z_{j}=E+\sum_{k=1,k\neq j}^{d}\Delta^{j,k}\prod_{\ell=1,\ell\neq
k}^{d}\mathbf{1}_{\Delta^{\ell,k}<0},\mbox{ for all }j=1,\ldots,d,$ (8)

where $\mathbf{1}_{\cdot}$ denotes the indicator function.

Furthermore, for a fixed $q\in\\{1,\dots,d\\}$, the $q$-th vector of
differences is defined as follows

$\bm{\Delta}^{(q)}=\left(\Delta^{q,1},\dots,\Delta^{q,d}\right)\in\mathbb{R}^{d}.$

With Representation (8) at our disposal, we may proceed to the simulation of a
MGP vector $\bm{Z}$, via the simulation of $\bm{\Delta}^{(q)}$ for a fixed
$q\in\\{1,\dots,d\\}$ using bootstrapping techniques combined with a rejection
algorithm in some specific cases.

###  3.2 Non-parametric joint MGP simulation

In this section, without any loss of generality, we fix $q$ equal to 1. We
then describe the joint simulation procedure, which is a generalisation of the
algorithm originally proposed by Legrand et al. (2023) for the case where
$d=2$.

For a given sample $(\bm{Z}_{i})_{1\leq i\leq n}$ of independent and
identically distributed (i.i.d.) replicates of a MGP distributed vector
$\bm{Z}$, the joint MGP simulation algorithm (see Algorithm 1) performs the
stochastic generation of $m\in\mathbb{N}$ new i.i.d. copies of $\bm{Z}$,
denoted $(\widetilde{\bm{Z}}_{\ell})_{1\leq\ell\leq m}$.

Algorithm 1 is as follows. First, we simulate unit exponential random
variables. Then, independently on this first step, we simulate new
realisations $(\widetilde{\bm{\Delta}}^{(1)}_{\ell})_{1\leq\ell\leq m}$ of the
differences $\bm{\Delta}^{(1)}$ using a non-parametric bootstrap approach,
i.e. by resampling among the observed differences $\bm{\Delta}^{(1)}$.
Finally, we use the stochastic relation in Equation (8) to merge both
simulation steps, in order to generate new realisations of MGP vectors, noting
that $\Delta^{r,s}=\Delta^{1,s}-\Delta^{1,r}$, for all $1\leq r,s\leq d$.

Algorithm 1 Non-parametric joint MGP simulation.

1: Observations $\left(\bm{Z}_{i}\right)_{1\leq i\leq
n}=\left(Z_{i,1},\ldots,Z_{i,d}\right)_{1\leq i\leq n}$ from a standard MGP
distribution;

2: A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$

3: A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ procedure A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ procedure

4: A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ A standard MGP simulated sample
$\left(\widetilde{\bm{Z}}_{m}\right)_{1\leq\ell\leq
m}=\left(\widetilde{Z}_{\ell,1},\ldots,\widetilde{Z}_{\ell,d}\right)_{1\leq\ell\leq
m}$ procedure $\Delta_{i}^{1,k}\leftarrow Z_{i,1}-Z_{i,k}$, for $1\leq i\leq
n$ and $1\leq k\leq d$ $\Delta_{i}^{1,k}\leftarrow Z_{i,1}-Z_{i,k}$, for
$1\leq i\leq n$ and $1\leq k\leq d$

5: $\Delta_{i}^{1,k}\leftarrow Z_{i,1}-Z_{i,k}$, for $1\leq i\leq n$ and
$1\leq k\leq d$ Generate
$E_{1},\ldots,E_{m}\overset{\text{i.i.d.}}{\sim}\mathrm{Exp}(1)$, and
independent of $(\Delta_{i}^{1,k})_{1\leq i\leq n,1\leq k\leq d}$ Generate
$E_{1},\ldots,E_{m}\overset{\text{i.i.d.}}{\sim}\mathrm{Exp}(1)$, and
independent of $(\Delta_{i}^{1,k})_{1\leq i\leq n,1\leq k\leq d}$

6: Generate $E_{1},\ldots,E_{m}\overset{\text{i.i.d.}}{\sim}\mathrm{Exp}(1)$,
and independent of $(\Delta_{i}^{1,k})_{1\leq i\leq n,1\leq k\leq d}$ Generate
a $m$-bootstrap sample
$\widetilde{\bm{\Delta}}^{(1)}=\left(\widetilde{\bm{\Delta}}^{(1)}_{\ell}\right)_{\ell=1,\ldots,m}$
from $\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$ Generate a
$m$-bootstrap sample
$\widetilde{\bm{\Delta}}^{(1)}=\left(\widetilde{\bm{\Delta}}^{(1)}_{\ell}\right)_{\ell=1,\ldots,m}$
from $\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$

7: Generate a $m$-bootstrap sample
$\widetilde{\bm{\Delta}}^{(1)}=\left(\widetilde{\bm{\Delta}}^{(1)}_{\ell}\right)_{\ell=1,\ldots,m}$
from $\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$
$\widetilde{\Delta}_{\ell}^{r,s}\leftarrow\widetilde{\Delta}_{\ell}^{1,s}-\widetilde{\Delta}_{\ell}^{1,r}$,
for $1\leq\ell\leq m$ and all $1\leq r,s\leq d$
$\widetilde{\Delta}_{\ell}^{r,s}\leftarrow\widetilde{\Delta}_{\ell}^{1,s}-\widetilde{\Delta}_{\ell}^{1,r}$,
for $1\leq\ell\leq m$ and all $1\leq r,s\leq d$

8:
$\widetilde{\Delta}_{\ell}^{r,s}\leftarrow\widetilde{\Delta}_{\ell}^{1,s}-\widetilde{\Delta}_{\ell}^{1,r}$,
for $1\leq\ell\leq m$ and all $1\leq r,s\leq d$
$\widetilde{Z}_{\ell,j}\leftarrow E_{\ell}+\sum_{s=1,s\neq
j}^{d}\widetilde{\Delta}^{j,s}_{\ell}\prod_{r=1,r\neq
s}^{d}\mathbf{1}_{\widetilde{\Delta}^{r,s}_{\ell}<0}$ for all $1\leq\ell\leq
m$ and $1\leq j\leq d$ $\widetilde{Z}_{\ell,j}\leftarrow
E_{\ell}+\sum_{s=1,s\neq j}^{d}\widetilde{\Delta}^{j,s}_{\ell}\prod_{r=1,r\neq
s}^{d}\mathbf{1}_{\widetilde{\Delta}^{r,s}_{\ell}<0}$ for all $1\leq\ell\leq
m$ and $1\leq j\leq d$

9: $\widetilde{Z}_{\ell,j}\leftarrow E_{\ell}+\sum_{s=1,s\neq
j}^{d}\widetilde{\Delta}^{j,s}_{\ell}\prod_{r=1,r\neq
s}^{d}\mathbf{1}_{\widetilde{\Delta}^{r,s}_{\ell}<0}$ for all $1\leq\ell\leq
m$ and $1\leq j\leq d$ $\widetilde{Z}_{\ell,j}\leftarrow
E_{\ell}+\sum_{s=1,s\neq j}^{d}\widetilde{\Delta}^{j,s}_{\ell}\prod_{r=1,r\neq
s}^{d}\mathbf{1}_{\widetilde{\Delta}^{r,s}_{\ell}<0}$ for all $1\leq\ell\leq
m$ and $1\leq j\leq d$ end procedure

Proposition 7.1 (see Section 7) guarantees that the samples simulated with
Algorithm 1 are indeed distributed according to a standard MGP distribution,
for all $q\in\\{1,\dots,d\\}$, and in particular for $q=1$. The proof can be
found in Section 7 with the statement of the proposition.

We then provide a numerical illustration of Algorithm 1 in which we consider a
3-dimensional MGP vector $\bm{Z}$, with $\bm{T}$ (see Equation (6))
distributed according to a centred multivariate Gaussian distribution with
correlation coefficients $\rho_{1,2}=0.4,\rho_{1,3}=0.8$ and $\rho_{2,3}=0.1$.

|  |
---|---|---
a) | b) | c)
|  |
d) | e) | f)
Figure 1: Bivariate representations of the original sample of size 2,000 of a
MGP vector $\bm{Z}\in\mathbb{R}^{3}$ with $\bm{T}$ (see Equation (6))
distributed according to a centred multivariate Gaussian distribution with
correlation coefficients $\rho_{1,2}=0.4,\rho_{1,3}=0.8$ and $\rho_{2,3}=0.1$
and the simulated sample $\widetilde{\bm{Z}}$ of size $10,000$. Figures a), b)
and c) display the QQ plots of ($Z_{1}$, $\widetilde{Z}_{1}$), ($Z_{2}$,
$\widetilde{Z}_{2}$) and ($Z_{3}$, $\widetilde{Z}_{3}$), respectively, and
their associated 95% point-wise confidence intervals based on $1,000$
bootstrap replications. Figures d), e) et f) represent scatter plots of the
components ($Z_{1}$, $Z_{2}$), ($Z_{1}$, $Z_{3}$) and ($Z_{2}$, $Z_{3}$) of
the original sample (black crosses). On each scatter plot, the points of the
simulated sample (red circles) are added.

Figure 1 presents a bivariate representation of a simulated sample of size
$10,000$ and the original sample of the 3-dimensional vector $\bm{Z}$ of size
$2,000$. The top row (Figures a)-c)) of Figure 1 depicts the quantile-quantile
(QQ) plots between the original observations and simulated observations for
each component. The QQ plots show a good fit between the marginal
distributions of the simulated samples $\widetilde{Z}_{j}$ and the marginal
distributions of the original sample $Z_{j}$ for $j\in\\{1,2,3\\}$.
Consequently, the QQ plots indicate that the joint simulation algorithm is
capable of generating samples with the same marginal distributions as the ones
of the original sample. The scatter plots in the bottom row of Figure 1 show
the simulated samples of $\widetilde{\bm{Z}}$ in red and the original sample
of $\bm{Z}$ in black. Each column of the second row focuses on a single plane
of the $3$-dimensional space, and hence on a bundle of the component of
$\bm{Z}$ and $\widetilde{\bm{Z}}$. With regard to the dependence structure,
the scatter plots of Figure 1 show that the shape displayed by the simulated
sample seems to be similar to the one of the original sample. It can be
observed that the simulated sample comprises a greater number of observations
in the extreme regions.

###  3.3 Non-parametric conditional MGP simulation

Thanks to the representation of the univariate components $Z_{j}$ of a MGP
vector $\bm{Z}$ by Equation (8), we were able to derive our joint simulation
algorithm. In the same way, but using the definition of $\Delta^{j,k}$ in
Equation (7), we propose a second algorithm, namely the conditional MGP
simulation algorithm (see Algorithm 2). In fact, for practical application, we
are also interested in generating observations of $Z_{j}$, when $\bm{Z}_{-j}$
is observed at some extreme level $\bm{z}_{-j}$. To be able to simulate
$Z_{j}$ conditionally on the rest of the components $\bm{Z}_{-j}$, the
conditional distribution of $\Delta^{q,j}$ given $\bm{Z}_{-j}$ must be
computed for any fixed $q=1,\ldots,d$ with $q\neq j$. Then, realisations of
$\Delta^{q,j}$ are drawn, which can be used to compute the new observations of
$Z_{j}$ according to Equation (8).

The procedure of non-parametric conditional MGP simulation is based on the
conditional distribution of $\Delta^{q,j}$ given $\bm{Z}_{-j}$. Denoting
$z_{\star}\coloneqq\max\bm{z}_{-j}$, the conditional distribution of
$\Delta^{q,j}$ given $\bm{Z}_{-j}$ can then be expressed as follows

  1. 1.

If $z_{\star}>0$ and $z_{\star}=z_{q}$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(0)+I_{2}(0)}\left(\mathbf{1}_{\delta^{q,j}>0}+\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta^{q,j}\leq
0}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right).$

  2. 2.

If $z_{\star}>0$ and $z_{\star}\neq z_{q}$, denoting $\delta_{\star}\coloneqq
z_{q}-z_{\star}$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q}-z_{\star})+I_{2}(z_{q}-z_{\star})}\left(\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta_{\star}>\delta^{q,j}}+\mathrm{e}^{\delta_{\star}}\mathbf{1}_{\delta_{\star}\leq\delta^{q,j}}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right).$

  3. 3.

If $z_{\star}\leq 0$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q})}\mathrm{e}^{\delta^{q,j}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{\delta^{q,j}<z_{q}}.$

where

  * •

$\bm{\delta}^{(q)}=(\delta^{(q,1)},\dots,\delta^{(q,d)})$

  * •

for any $x\in\mathbb{R}$,

$I_{1}(x)\coloneqq\int_{-\infty}^{x}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\quad\text{
and }\quad
I_{2}(x)\coloneqq\mathrm{e}^{x}\int_{x}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta$

where $\bar{\bm{\delta}}^{(q)}$ has the same components as $\bm{\delta}^{(q)}$
except for $\delta^{q,j}$ which is replaced by $\zeta=z_{q}-z_{j}$.

The theoretical result and the proof can be found in Section 7.2 (Proposition
7.4). It can be observed that in Cases 1 and 3, the conditional distribution
of $\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}$ is independent of $\bm{z}_{-j}$.
Hence, in these specific cases, bootstrapping will be used for the conditional
simulation of $\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}$. However, in Case 2
($z_{\star}>0$ and $z_{\star}\neq z_{q}$), the conditional distribution of
$\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}$ does depend on $\bm{z}_{-j}$. In
that case, bootstrapping needs to be combined with the rejection sampling
method. However, it is not common practice to use the rejection algorithm to
sample a univariate variable $\Delta^{q,j}$ through observations of random
vector $\bm{\Delta}^{(q)}$. Therefore, Propositions 7.5 and 7.6 are devoted to
the rejection algorithm and its application in the conditional MGP simulation
procedure. Finally, Algorithm 2 outlines the procedure in question.

Algorithm 2 Non-parametric conditional MGP simulation

1: Observations $\left(\bm{\Delta}^{(1)}_{i}\right)_{1\leq i\leq n}$;
$j\in\\{2,\dots,d\\}$ the index of risk factor of interest; a realisation
$\bm{z}_{-j}$ of $\bm{Z}_{-j}$

2: A simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$

3: A simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure

4: A simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then A
simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then

5: A simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then A
simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then A
simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then if
$z_{\star}=z_{1}$ then A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then A
simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then if
$z_{\star}=z_{1}$ then

6: A simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then A
simulated sample $\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$
conditionally on $\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ A simulated sample
$\left(\widetilde{Z}_{\ell,j}\right)_{1\leq\ell\leq m}$ conditionally on
$\bm{Z}_{-j}=\bm{z}_{-j}$ procedure if $\exists k\neq j\mid z_{k}>0$ then if
$z_{\star}=z_{1}$ then Define $\bm{\Delta}^{(1)}_{|z_{\star}=z_{1}}$, the
subset of $\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$ such that
$\max\bm{z}_{-j}=z_{1}$ Define $\bm{\Delta}^{(1)}_{|z_{\star}=z_{1}}$, the
subset of $\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$ such that
$\max\bm{z}_{-j}=z_{1}$

7: Define $\bm{\Delta}^{(1)}_{|z_{\star}=z_{1}}$, the subset of
$\left(\bm{\Delta}_{i}^{(1)}\right)_{1\leq i\leq n}$ such that
$\max\bm{z}_{-j}=z_{1}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$

8: Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else

9: Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ else Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else for $1\leq\ell\leq m$ do Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else for $1\leq\ell\leq m$ do

10: Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
Bootstrap $m$ realisations $\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq
m}$, from $\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of
$\bm{Z}_{-j}$ Bootstrap $m$ realisations
$\left(\Delta_{\ell}^{1,j}\right)_{1\leq\ell\leq m}$, from
$\left(\Delta^{1,j}_{|z_{\star}=z_{1}}\right)$ independently of $\bm{Z}_{-j}$
else for $1\leq\ell\leq m$ do Sample a realisation $\Delta^{1,j}_{m}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$
Sample a realisation $\Delta^{1,j}_{m}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$

11: Sample a realisation $\Delta^{1,j}_{m}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$
Compute $\Delta_{\star}:=z_{1}-z_{\star}$ Compute
$\Delta_{\star}:=z_{1}-z_{\star}$

12: Compute $\Delta_{\star}:=z_{1}-z_{\star}$ Generate a random number
$u\in[0,1]$ Generate a random number $u\in[0,1]$

13: Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$
Generate a random number $u\in[0,1]$ while
$u>\left(e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta_{\star}>\Delta^{1,j}_{\ell}}+e^{\Delta_{\star}}\mathbf{1}_{\Delta_{\star}\leq\Delta^{1,j}_{\ell}}\right)$
do Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$
while
$u>\left(e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta_{\star}>\Delta^{1,j}_{\ell}}+e^{\Delta_{\star}}\mathbf{1}_{\Delta_{\star}\leq\Delta^{1,j}_{\ell}}\right)$
do

14: Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$
while
$u>\left(e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta_{\star}>\Delta^{1,j}_{\ell}}+e^{\Delta_{\star}}\mathbf{1}_{\Delta_{\star}\leq\Delta^{1,j}_{\ell}}\right)$
do Repeat last two steps Repeat last two steps

15: Repeat last two steps Repeat last two steps end while Repeat last two
steps end while

16: Repeat last two steps end while Repeat last two steps end while end for
Repeat last two steps end while end for

17: Repeat last two steps end while end for Repeat last two steps end while
end for end if Repeat last two steps end while end for end if

18: Repeat last two steps end while end for end if Repeat last two steps end
while end for end if Repeat last two steps end while end for end if Repeat
last two steps end while end for end if else Repeat last two steps end while
end for end if Repeat last two steps end while end for end if Repeat last two
steps end while end for end if else

19: Repeat last two steps end while end for end if Repeat last two steps end
while end for end if Repeat last two steps end while end for end if else
Repeat last two steps end while end for end if Repeat last two steps end while
end for end if Repeat last two steps end while end for end if else Repeat last
two steps end while end for end if Repeat last two steps end while end for end
if Repeat last two steps end while end for end if else for $1\leq\ell\leq m$
do Repeat last two steps end while end for end if Repeat last two steps end
while end for end if Repeat last two steps end while end for end if else
Repeat last two steps end while end for end if Repeat last two steps end while
end for end if Repeat last two steps end while end for end if else for
$1\leq\ell\leq m$ do

20: Repeat last two steps end while end for end if Repeat last two steps end
while end for end if Repeat last two steps end while end for end if else
Repeat last two steps end while end for end if Repeat last two steps end while
end for end if Repeat last two steps end while end for end if else for
$1\leq\ell\leq m$ do Sample a realisation $\Delta^{1,j}_{\ell}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$
Sample a realisation $\Delta^{1,j}_{\ell}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$

21: Sample a realisation $\Delta^{1,j}_{\ell}$ from
$\left(\Delta^{1,j}_{i}\right)_{1\leq i\leq n}$ independently of $\bm{Z}_{-j}$
Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$

22: Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$
Generate a random number $u\in[0,1]$ while
$u>e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta^{1,j}_{\ell}<z_{q}}$ do Generate
a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$ while
$u>e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta^{1,j}_{\ell}<z_{q}}$ do

23: Generate a random number $u\in[0,1]$ Generate a random number $u\in[0,1]$
while $u>e^{\Delta^{1,j}_{\ell}}\mathbf{1}_{\Delta^{1,j}_{\ell}<z_{q}}$ do
Repeat last two steps Repeat last two steps

24: Repeat last two steps Repeat last two steps end while Repeat last two
steps end while

25: Repeat last two steps end while Repeat last two steps end while end for
Repeat last two steps end while end for

26: Repeat last two steps end while end for Repeat last two steps end while
end for end if Repeat last two steps end while end for end if

27: Repeat last two steps end while end for end if
$\widetilde{Z}_{\ell,j}\leftarrow z_{1}-\Delta_{\ell}^{1,j}$ for all
$1\leq\ell\leq m$ $\widetilde{Z}_{\ell,j}\leftarrow z_{1}-\Delta_{\ell}^{1,j}$
for all $1\leq\ell\leq m$

28: $\widetilde{Z}_{\ell,j}\leftarrow z_{1}-\Delta_{\ell}^{1,j}$ for all
$1\leq\ell\leq m$ $\widetilde{Z}_{\ell,j}\leftarrow z_{1}-\Delta_{\ell}^{1,j}$
for all $1\leq\ell\leq m$ end procedure

We illustrate the performances of the conditional MGP simulation Algorithm 2
through the same parametric example than in Section 3.3, that is a
3-dimensional MGP vector $\bm{Z}$, with $\bm{T}$ (see Equation (6))
distributed according to a centred multivariate Gaussian distribution, but
with different correlation coefficients $\rho_{1,2}=0.6,\rho_{1,3}=0.8$ and
$\rho_{2,3}=0.5$. We consider the conditional distribution of
$Z_{2}\mid\bm{Z}_{-2}=\bm{z}_{-2}$. The numerical experiment consists in
generating realisations of the $2$-th component of a 3-dimensional standard
MGP vector given observations $\bm{z}_{-2}=(z_{1},z_{3})$ for the original
sample $\bm{Z}$ of size 2,000 and a sample simulated using Algorithm 2 of size
10,000. We illustrate each case (Cases 1, 2 and 3) depending on
$\bm{z}_{\star}$. For each simulation case, the empirical density of simulated
samples of $Z_{2}\mid\bm{Z}_{-2}=\bm{z}_{-2}$ is compared to the theoretical
density. Figure 2 shows that regardless the value of $\bm{z}_{-2}$, the
histogram shape and the kernel estimates of the density of the simulated
samples of $Z_{2}\mid\bm{Z}_{-2}=\bm{z}_{-2}$ match systematically the
theoretical density curve.

|  |
---|---|---
$\bm{z}_{-2}=(0.54,0.31)$ | $\bm{z}_{-2}=(1.14,1.03)$ | $\bm{z}_{-2}=(0.21,0.11)$
|  a) Case 1: $z_{\star}>0$ and $z_{\star}=z_{1}$ |
|  |
$\bm{z}_{-2}=(0.17,0.22)$ | $\bm{z}_{-2}=(0.24,0.79)$ | $\bm{z}_{-2}=(1.37,1.40)$
|  b) Case 2: $z_{\star}>0$ and $z_{\star}\neq z_{1}$ |
|  |
$\bm{z}_{-2}=(-0.42,-0.35)$ | $\bm{z}_{-2}=(-0.86,-1.00)$ | $\bm{z}_{-2}=(-0.67,-0.49)$
|  c) Case 3: $z_{\star}\leq 0$ |
Figure 2: Histograms of simulated samples $Z_{2}\mid\bm{Z}_{-2}=\bm{z}_{-2}$
of size $m=10,000$. Kernel density estimates (curves in black) and and
theoretical densities of the original sample (curves in red). Figures a)
correspond to Case 1 with $z_{\star}>0$ and $z_{\star}=z_{1}$, Figures b) to
Case 2 with $z_{\star}>0$ and $z_{\star}\neq z_{1}$ and Figures c) to Case 3
with $z_{\star}<0$.

##  4 Illustrations on simulated data

As motivated in the introduction, the main goal of the MGP-based simulation
approaches is to address the problem of scarcity of available observations in
extreme regions. These regions will take different forms depending on the
application of interest.

In the following, we illustrate the two simulation algorithms derived in
Section 3 on a numerical setting inspired by the field of finance and more
specifically financial returns (details are given below). This simulation
framework is presented in Section 4.1. The joint MGP simulation algorithm
(Algorithm 1) is then applied in Section 4.2 in order to estimate the three
TRMs presented in Section 2 from the simulated data. Finally, in Section 4.3,
the conditional MGP simulation algorithm (Algorithm 2) is applied to estimate
financial shocks, defined as a specific conditional expectation.

###  4.1 Simulation framework

We consider a 3-dimensional random vector $\bm{X}=(X_{1},X_{2},X_{3})$ with
marginal Student’s $t$-distribution with a relatively low degree of freedom
$\nu$. This choice of distribution for the marginals is chosen so as to mimic
financial returns, which are typically heavy tailed. Indeed, the tail of the
Student’s $t$-distribution is governed by $1/\nu$. So, the smaller $\nu$ the
heavier the tail.

To fully characterise the joint distribution of $\bm{X}$, we need to define
the dependence structure in addition to the marginal distribution. Sklar’s
theorem (Sklar, 1959) states that given the continuous joint distribution
$\bm{F}$ of a random vector $\bm{X}$ with marginal cumulative distribution
functions $F_{1},F_{2},F_{3}$, there exists a unique copula $\mathcal{C}$ such
that

$\displaystyle\bm{F}(x_{1},x_{2},x_{3})=\mathcal{C}\left(F_{1}(x_{1}),F_{2}(x_{2}),F_{3}(x_{3})\right).$

Copulas are powerful tools for modelling the dependence structure of a random
vector separately from the marginal distributions.

Let us also recall that an underlying assumption of our simulation framework
is that the components of $\bm{X}$ are asymptotically dependent. To ensure
that this hypothesis is satisfied, we consider the Gumbel copula (Nelsen,
2006) to obtain dependent extremes in the upper tail. In the $3$-dimensional
case, the Gumbel copula is defined as follows

$\displaystyle\mathcal{C}(\bm{y}):=\exp\left(-\left(\sum_{i=1}^{3}\left[-\log(y_{i})\right]^{\theta}\right)^{1/\theta}\right),$
(9)

where $\theta\geq 1$ is the copula parameter. The larger $\theta$, the
stronger the asymptotic dependence structure between the components of
$\bm{X}$. Consequently, the pair $\left(\bm{\nu},\theta\right)$, where
$\bm{\nu}=(\nu_{1},\nu_{2},\nu_{3})$ represents the vector of degrees of
freedom associated with each component of $\bm{X}$ and $\theta$ parameterise
the structure of dependence through Equation (9), fully characterises the
joint distribution of $\bm{X}$.

The numerical experiments, unless otherwise specified, are performed on
simulated data sets $\mathcal{D}\in\mathbb{R}^{1500\times 3}$, with
$\nu_{1}=2$, $\nu_{2}=3$, $\nu_{3}=2.5$. From this parametric framework, we
can derive the theoretical values of the TRMs and infer quantities describing
the tail of some conditional distribution, that will be used as benchmarks in
the following sections.

As previously mentioned, Algorithms 1 and 2 have to run on standard MGP
vectors. To this end, the vectors $\bm{X}$ are transformed to standard MGP
vectors $\bm{Z}$, following the procedure presented in Section 3.1.
Conversely, once the algorithms have been run on $\bm{Z}$, the resulting
simulated sample outputs of the algorithms $\widetilde{\bm{Z}}$ are projected
back into the original space through

$\widetilde{X}_{j}=F_{j}^{-1}\left(1-e^{\widetilde{Z}_{j}+u^{E}_{j}}\right),\mbox{
for }j\in\\{1,\ldots,d\\}$

where $F^{-1}_{j}$ is the quantile function associated with the $j$-th risk
factor for $j=1,\ldots,d$. $F^{-1}_{j}$ can be the quantile function of any
parametric distribution that has been fitted to the data, in our simulation
framework, that is the Student’s $t$-distribution quantile distribution. Thus,
any quantity to be estimated on the $\bm{X}$-scale can be estimated on the
$\widetilde{\bm{X}}$-samples.

For comparison purpose, the estimation procedure is performed on the following
three different samples.

  * •

Original sample $\mathcal{D}$ generated using the parametric framework above,
denoted Orig;

  * •

Simulated sample $\mathcal{\widetilde{D}}$ obtained using the joint simulation
algorithm (Algorithm 1), denoted Simu;

  * •

Extended sample $\mathcal{D}\cup\mathcal{\widetilde{D}}$ corresponding to the
concatenation of the original and simulated samples, denoted Ext.

###  4.2 Tail-related risk metrics estimation through the joint simulation of
multivariate extremes

In this section, we investigate the performance of our joint simulation
procedure (Algorithm 1) in estimating TRMs, as defined in Equations (1), (2)
and (3), with respect to the several parameters. Particular attention is paid
to the strength of the asymptotic dependence, which can be controlled through
the copula parameter $\theta$ and the level $\alpha$ (see Supplementary
material, Section 1).

In the following, TRMs estimates are compared to their respective theoretical
benchmark values (see Supplementary Material, Section 2). For the sake of
brevity, we only considered the estimation procedure for the component
$X_{1}$. However, it is important to note that it could also be performed for
the other components, $X_{2}$ and $X_{3}$.

In order to run the joint simulation algorithm (Algorithm 1), it is necessary
to specify the size of the simulated samples $m$. The simulation procedure was
run with $m=10,000$ which yielded stable results while exhibiting a reasonable
computational cost. Further results regarding the performance of the
simulation algorithm with respect to the size $m$ can be found in the
Supplementary Material (Section 3).

The estimation of the $\mathrm{VaR}$ is of paramount importance in estimating
the TRM, as all TRMs are defined from the $\mathrm{VaR}$ (see Section 2).
However, this is not the focus of our current investigation. The performance
of our approach is evaluated using the theoretical $\mathrm{VaR}$, that is to
say, the $\mathrm{VaR}$ computed as a Student’s $t$-quantile with the
corresponding degrees of freedom. Other techniques for estimating the
$\mathrm{VaR}$ have been investigated, including an empirical approach and the
estimation of the $\mathrm{VaR}$ as a univariate extreme value distribution
quantile. The results of these analyses can be found in the Supplementary
Material (Section 4). In particular, it can be seen that when the theoretical
$\mathrm{VaR}$ is unknown, which is often the case in real-data applications,
the estimation of the $\mathrm{VaR}$ using a univariate extreme value approach
appears to be a good candidate.

Given that the $\mathrm{VaR}$ estimation technique and $m$ are fixed, we can
now proceed to investigate the performance of our approach with respect to
$\alpha$ and $\theta$. Moreover, since the estimation errors are sensitive to
the sample $\mathcal{D}$, it is necessary to ascertain that the approach
performance is not specific to the original sample $\mathcal{D}$. To this end,
50 distinct original samples $\mathcal{D}$ were simulated (of sample size $n$
and with the same joint distribution as outlined in Section 4.1). Then, for
each original sample $\mathcal{D}$, a set of 50 simulated samples
$\widetilde{\mathcal{D}}$ of size $m=10,000$ was generated following the joint
simulation algorithm (Algorithm 1). The TRMs are estimated on each of 50
original, simulated and extended samples, along with the relative errors on
these TRMs approximations.

Figure 3 and Table 1 summarise the results for TRMs estimations at different
levels $\alpha\in\\{0.9975,0.999,0.9997\\}$ and for different copula
parameters $\theta\in\\{1.3,2.6,7.3\\}$. Boxplots of the distributions of the
relative errors are presented, showing that the relative errors for the TRMs
estimations on the simulated samples Simu and on the extended samples Ext are
smaller and present less dispersion than those made on the original samples
Orig. It is also important to ascertain that the quality of TRMs do not
deteriorate when $\alpha$ gets closer to 1. Figure 3 shows TRMs estimations at
level $\alpha\in\\{0.9975,0.999,0.9997\\}$, one may observe that the
estimation on Simu and Ext samples induce fairly low approximation errors
regardless of the level $\alpha$. We notice that the joint simulation of
multivariate extremes using Algorithm 1 improves not only the estimation of
$\mathrm{MMES}$ and $\mathrm{DCTE}$ (defined with the joint distribution of
$\bm{X}$) but also the estimation of $\mathrm{ES}$ even though it only evolves
the marginal distribution of $X_{j}$.

|  |
---|---|---
a) $\mathrm{ES}_{0.9975}$ |  b) $\mathrm{ES}_{0.999}$ |  c) $\mathrm{ES}_{0.9997}$
|  |
d) $\mathrm{MMES}_{0.9975}$ |  e) $\mathrm{MMES}_{0.999}$ |  f) $\mathrm{MMES}_{0.9997}$
|  |
g) $\mathrm{DCTE}_{0.9975}$ |  h) $\mathrm{DCTE}_{0.999}$ |  i) $\mathrm{DCTE}_{0.9997}$
Figure 3: Distribution of relative approximation errors of the estimations of TRMs on 50 original samples (grey oblique lines), 50 simulated samples (red dots) and 50 extended samples (yellow grid) for the $\mathrm{ES}_{\alpha}$ (Figures a), b), c)), the $\mathrm{MMES}_{\alpha}$ (Figures d), e), f)) and the $\mathrm{DCTE}_{\alpha}$ (Figures g), h), i)) at varying level $\alpha\in\\{0.9975,0.999,0.9997\\}$ with respect to copula parameter $\theta\in\\{1.3,2.6,7.3\\}$. On each graph, the dashed blue horizontal line corresponds to the associated theoretical value. | $\alpha$ | $\mathrm{ES}_{\alpha}$ | $\mathrm{MMES}_{\alpha}$ | $\mathrm{DCTE}_{\alpha}$
---|---|---|---
Orig | Simu | Orig | Simu | Orig | Simu
0.9975 | 3.7 (1.6) | 87.4 (5.2) | 0.82 (0.9) | 33.9 (2.2) | 1.3 (1.3) | 42.9 (3.3)
0.999 | 1.7 (1.4) | 35.0 (2.3) | 0.4 (0.6) | 13.5 (1.1) | 0.5 (0.7) | 17.2 (1.4)
0.9997 | 0.6 (0.7) | 10.4 (0.8) | 0.2 (0.49) | 4.0 (0.4) | 0.2 (0.5) | 5.1 (0.5)
a) $\theta=1.3$
| $\alpha$ | $\mathrm{ES}_{\alpha}$ | $\mathrm{MMES}_{\alpha}$ | $\mathrm{DCTE}_{\alpha}$
---|---|---|---
Orig | Simu | Orig | Simu | Orig | Simu
0.9975 | 3.9 (2.4) | 114.3 (6.6) | 2.3 (1.6) | 74.1 (3.9) | 2.7 (1.7) | 83.2 (4.6)
0.999 | 1.6 (1.2) | 46.0 (2.8) | 0.9 (0.9) | 29.8 (1.7) | 1.0 (0.9) | 33.5 (2.1)
0.9997 | 0.4 (0.6) | 13.6 (0.7) | 0.2 (0.5) | 8.8 (0.53) | 0.3 (0.5) | 10.0 (0.6)
b) $\theta=2.6$
| $\alpha$ | $\mathrm{ES}_{\alpha}$ | $\mathrm{MMES}_{\alpha}$ | $\mathrm{DCTE}_{\alpha}$
---|---|---|---
Orig | Simu | Orig | Simu | Orig | Simu
0.9975 | 4.1 (2.1) | 144.8 (10.0) | 3.4 (1.9) | 126.0 (8.5) | 3.5 (2.0) | 131.0 (8.9)
0.999 | 1.6 (1.4) | 57.7 (4.1) | 1.5 (1.3) | 50.2 (3.38) | 1.5 (1.4) | 52.3 (3.6)
0.9997 | 0.4 (0.6) | 17.3 (10.0) | 0.3 (0.5) | 15.1 (1.1) | 0.3 (0.5) | 15.7 (1.2)
c) $\theta=7.3$
Table 1: Mean values and standard deviations in brackets of the number of
observations on which the TRMs are estimated on the original samples (Orig)
and on the simulated samples (Simu) for different values of
$\alpha\in\\{0.9975,0.999,0.9997\\}$ and different copula parameters a)
$\theta=1.3$, b) $\theta=2.6$ and c) $\theta=1.3$.

The accuracy of the TRMs estimates is attributable to the enlargement of the
sample size on which TRMs are being estimated compared to the original sample
size. Table 1 shows the mean value of the sample size on which the three TRMs
are estimated when considering the original samples and simulated samples for
different levels $\alpha\in\\{0.9975,0.999,0.9997\\}$ and for different values
of copula parameters $\theta\in\\{1.3,2.6,7.3\\}$. Overall, it can be noted
that for the considered levels of $\alpha$, the TRMs estimation on the
original samples is conducted with limited data, that is, here with at most 5
observations. In contrast, our simulation approach allows for a larger sample
size to estimate the various TRMs, extending the original observations to 100
observations in some cases. As expected, the number of points above the
$\mathrm{VaR}$ increases when the asymptotic dependence is stronger.

###  4.3 Conditional expectation estimation through the conditional simulation
of multivariate extremes

In Section 4.2, we proposed a numerical illustration of our first Algorithm 1,
the joint simulation algorithm, through the empirical estimation of TRMs. As
previously stated in the introduction, the second algorithm 2 can be used to
infer quantities describing the tail of some conditional distribution. Since
the TRMs under investigation can be expressed as conditional expectations, one
might be inclined to consider the conditional simulation algorithm for their
estimation. However, accurate estimation of TRMs with the conditional
simulation algorithm (Algorithm 2) is not possible. Indeed, for a MGP vector
$\bm{Z}$, the aforementioned algorithm generates samples of
$Z_{j}\mid\bm{Z}_{-j}=\bm{z}_{-j}$ for $j=1,\ldots,d$. However, when
considering the estimation of $\mathrm{MMES}$ for instance, we are attempting
to estimate the expectation of $Z_{j}\mid\bm{Z}_{-j}\geq\bm{z}_{-j}$ (see
Equation (2)).

Nevertheless, the conditional simulation algorithm could be employed in other
contexts. For a vector of risk factors $\bm{X}$, we consider here the
inference of the expectation of risk factor $X_{j}$ given that the other risk
factors $\bm{X}_{-j}$ are observed, that is $\bm{X}_{-j}=\bm{x}_{-j}$, defined
as follows

$\mu(\bm{x}_{-j})=\mathbb{E}\left[X_{j}\mid\bm{X}_{-j}=\bm{x}_{-j}\right].$
(10)

Algorithm 2, the conditional simulation algorithm, is illustrated on this
specific application using the simulated data generated with the parametric
framework specified in Section 4.1. In this simulation framework, the
conditional distribution of $X_{j}$ given $\bm{X}_{-j}=\bm{x}_{-j}$ can be
derived explicitly, thereby enabling the derivation of $\mu(\bm{x}_{-j})$ as
defined in Equation (10). Hence, a benchmark value is available, and can be
used to measure the accuracy of the empirical estimations using simulated
samples generated with Algorithm 2.

For each original sample $\mathcal{D}$, we consider the estimation of the
conditional expectation denoted by $\mu(\bm{x}_{-2})$ (see Equation (10)).
This estimation is based on the conditioning observations $\bm{x}_{-2}$, the
original sample $\mathcal{D}$ and the generated conditional simulated samples
of $X_{2}\mid\bm{X}_{-2}=\bm{x}_{-2}$. Once the vector $\bm{X}$ has been
transformed into a MGP vector $\bm{Z}$, the conditioning observations are
selected such that each conditional simulation case (see Section 3.3),
$z_{\star}>0$, $z_{\star}=z_{1}$ (Case 1), $z_{\star}>0$, $z_{\star}\neq
z_{1}$ (Case 2) and $z_{\star}\leq 0$ (Case 3), is represented. For each
$\bm{z}_{-2}$, 10 simulated samples $Z_{2}\mid\bm{Z}_{-2}=\bm{z}_{-2}$ of size
$m=10,000$ are generated, and $\mu(\bm{x}_{-2})$ is estimated as the empirical
mean of each simulated sample. The deviation of the mean over these estimates
from the reference value are computed for each conditioning observation
$\bm{x}_{-2}$. Then, Figure 4 depicts the distribution of the mean absolute
error across 50 original samples $\mathcal{D}$ for each conditional simulation
case, with the conditioning observations being increasingly extreme. In other
words, the selected conditioning observations are in the bulk of the
distribution for Figure 4 a), in the tail of the distribution for Figure 4 b),
and even further in the tail in Figure 4 c). For comparison purposes, a linear
regression model with $X_{2}$ as the response variable and $\bm{X}_{-2}$ as
the covariates is fitted. It is evident that the conditional simulation errors
are of smaller magnitude than those induced by linear regression, which result
in substantial errors. It should be noted that while errors exhibit large
amplitude as the conditioning observation becomes extreme with linear
regression, the situation is reversed when the conditional simulation is
considered. The lowest errors are observed for the most extreme conditioning
observations, illustrated in Figure 4 c).

|
---|---
a) $\alpha=0.975$ |  b) $\alpha=0.99$
c) $\alpha\rightarrow 1$ (i.e. the maximum)
Figure 4: Distribution of relative errors of the deviations from the reference
value of $\mathbb{E}\left[X_{2}\mid\bm{X}_{-2}=\bm{x}_{-2}\right]$ with
respect to each conditional simulation case, using conditional simulation
(dark grey boxplot) and linear regression (light orange boxplot) for
increasingly extreme conditioning observations, that is $\bm{x}_{-2}$ chosen
as three different $\alpha$-quantile of all-observations.

##  5 Illustrations on real financial data

In the previous section, the two non-parametric simulation algorithms 1 and 2
were illustrated on simulated data. These illustrations have demonstrated the
advantages of these procedures for the estimation of TRMs and the inference of
quantities involving the tail of some conditional distribution. This section
presents a more practical illustration based on real financial data.

These data consist of weekly negative returns on stocks from three UK banks:
HSBC, Lloyds and RBS. The corresponding time series are denoted $X_{1}$,
$X_{2}$ and $X_{3}$ respectively. The sample size of each time series is
$n=470$, with temporal resolution spanning from 29/10/2007 to 17/10/2016. The
data used in this study were initially extracted from Yahoo Finance, and
considered in (Kiriliouk et al., 2019).

As in Section 4, both algorithms are applied to this real data set
$\bm{X}=(X_{1},X_{2},X_{3})$ in two stages. Firstly, Algorithm 1 is employed
in order to estimate TRMs (Section 5.1). Secondly, Algorithm 2 is used to
estimate the conditional expectation (Equation (10), Section 5.2).

Prior to this, and as done with the simulated data, the first step is to
standardise the data to common margins (see Section 3.1). In our case, we have
chosen to fit a Student’s $t$-distribution on each risk factor, which is often
used with heavy-tailed financial data (see e.g. McNeil et al., 2015). However,
other choices could have been considered. Goodness-of-fit plots can be found
in the Supplementary Material (see Section 5) along with the estimated
parameters of the fitted Student’s $t$-distributions.

###  5.1 Tail related risk metrics estimations

A set of $100$ new samples of weekly negative returns of size $m=10,000$ were
generated using our joint simulation algorithm (Algorithm 1) in order to
perform empirical estimation of the three TRMs of interest. The results are
displayed in Table 2 along with the estimations on the original sample Orig,
i.e. without increasing the amount of available data. The performance
assessment conducted under the parametric framework has shown that TRMs
estimates on the simulated samples Simu were more precise than estimations on
the extended sample Ext, hence, we exclusively discuss results from the
simulated samples Simu in this section. Table 2 shows that the joint
simulation approach allows for the estimation of the TRMs that would be
unfeasible if only the original sample was considered. This is to be expected,
given that the objective of the simulation approach is to increase the number
of available observations in extreme regions in order to enable the
performance of empirical estimations in a reliable manner.

Stock | $\mathrm{ES}_{0.9975}$ | $\mathrm{MMES}_{0.9975}$ | $\mathrm{DCTE}_{0.9975}$
---|---|---|---
Orig | Simu | Orig | Simu | Orig | Simu
HSBC | 0.24 (-) | 0.28 (0.015) | NA | 0.29 (0.022) | 0.11 (-) | 0.32 (0.026)
LL | 0.55 (-) | 0.83 (0.078) | NA | 0.91 (0.122) | NA | 1.04 (0.150)
RBS | 0.62 (-) | 0.56 (0.030) | NA | 0.62 (0.053) | NA | 0.66 (0.058)
Table 2: Empirical mean estimates, and standard deviation in brackets, of the
TRMs for the three negative returns of interest, using the original data Orig
and using the simulation algorithm with $100$ generated samples of size
$m=10,000$ each Simu. For each TRM, a threshold level $\alpha=0.9975$ is
considered. NA meaning that there are no observations above the theoretical
$\mathrm{VaR}$, so that the computation of the TRM was not possible.

###  5.2 Conditional expectation estimations

We now turn to the second application of our simulation framework, namely the
estimation of the conditional expectation
$\mu(\bm{x}_{-j})=\mathbb{E}\left[X_{j}\mid\bm{X}_{-j}=\bm{x}_{-j}\right]$,
for $j\in\\{1,2,3\\}$. For illustrative purposes, the conditioning
observations $(\bm{x}_{-j})_{j\in\\{1,2,3\\}}$ are selected in such a way that
they cover and illustrate the three possible cases of the conditional
simulation algorithm (see Section 3.3). However, in practice, if this approach
is used, for instance, for shock inference, these values would typically be
provided by regulators.

In order to simulate samples of $X_{j}\mid\bm{X}_{-j}=\bm{x}_{-j}$, Algorithm
2 considers the differences $\Delta^{j,q}$ with $q\neq j$ (see Equation (8)).
The algorithm was presented with $q=1$ and thus $j\neq 1$. But for $j=1$, any
$q\neq 1$ can be selected. In our case, to simulate
$X_{1}\mid\bm{X}_{-1}=\bm{x}_{-1}$, $q$ was chosen to be equal to 3. To
illustrate the three conditional simulation cases identified in Section 3.3,
we consider values of the components of the conditioning vector $\bm{x}_{-j}$
as quantiles at different extreme levels of the marginal distribution of the
$X_{j}$. For each conditional simulation case and each stock, we consider only
a single conditioning observation vector, and we generate 100 samples of size
$m=10,000$ using the conditioning simulation Algorithm 2. Then
$\mu(\bm{x}_{-j})$ is estimated as the empirical mean of each simulated sample
with Algorithm 2. The mean and the standard deviation over the 100 estimates
are displayed in Table 3. As in the numerical illustration in Section 4.3, we
compare the results obtained through our approach with those obtained through
linear regression. The results are depicted in Table 3.

Table 3 shows that while the shocks estimates for the HSBC negative returns
using conditional simulation are comparable to those obtained with linear
regression, the discrepancy between these two shock estimated is more
pronounced for Llyods and RBS negative returns. These two negative returns
exhibit a heavier tail than HSBC. However, this type of distribution is
precisely the one where linear regression is not a suitable approach, and
other methodologies should be considered. The results obtained under the
simulation framework of Section 4.1, where the theoretical value was used as a
reference, indicated that the approximation errors were relatively low, in
contrast to those produced by linear regression. Consequently, for real data,
it might be preferable to use our shock estimates over those provided by
linear regression.

|  |  Case 1 ($z_{\star}>0,z_{\star}=z_{q}$) |  Case 2 ($z_{\star}>0,z_{\star}\neq z_{q}$) |  Case 3 ($z_{\star}\leq 0$)
---|---|---|---
HSBC |  0.1526 ($10^{-3}$) |  0.1094 ($10^{-4}$) |  0.0436 ($10^{-4}$)
LLyods |  0.3031 ($10^{-3}$) |  0.0880 ($10^{-3}$) |  0.0793 ($10^{-3}$)
RBS |  0.2421 ($10^{-3}$) |  0.3273 ($10^{-3}$) |  0.0951 ($10^{-4}$)
a) Conditional Simulation
|  |  Case 1 ($z_{\star}>0,z_{\star}=z_{q}$) |  Case 2 ($z_{\star}>0,z_{\star}\neq z_{q}$) |  Case 3 ($z_{\star}\leq 0$)
---|---|---|---
HSBC | 0.1481 | 0.1051 | 0.0187
LLyods | 0.2042 | 0.1417 | 0.0295
RBS | 0.1475 | 0.3931 | 0.0478
b) Linear Regression
Table 3: Empirical mean estimates, and standard deviation in brackets when
available, of the estimated shocks
$\mu(\bm{x}_{-j})=\mathbb{E}\left[X_{j}\mid\bm{X}_{-j}=\bm{x}_{-j}\right]$,
$j\in\\{1,2,3\\}$, with $X_{1}$ (resp. $X_{2}$ and $X_{3}$) corresponding to
negative returns from HSBC (resp. Lloyds and RBS). Estimations are performed
using our conditional simulation approach (Algorithm 2) and b) linear
regression. For both inference approaches are illustrated through the three
cases identified for the conditional simulation approach.

##  6 Conclusion and Perspectives

One of the primary concerns when attempting to estimate risks at high levels
is the scarcity of available data. In fact, data sparseness makes any
estimation prone to instabilities and inconsistencies, thereby limiting the
effectiveness of this estimate as a trustworthy risk metric. Another
consequence of considering solely historical data is the fact that the
information contained in historical observations is incomplete. Essentially,
this approach represents only past crisis information, implying that upcoming
crises cannot exceed the severity of those that have already been experienced.

The issue was addressed in this paper by the development of two non-parametric
simulation approaches of multivariate extremes. Both algorithms are driven by
multivariate extreme value theory. The first algorithm, referred to as the
joint simulation algorithm, generates MGP vectors. The second algorithm,
referred to as the conditional simulation algorithm, generates samples of a
component of MGP vectors conditionally on the others components being extreme.
The aforementioned algorithms have enabled the expansion of the available data
set, and the use of the newly generated observations has resulted in the
accurate estimation of tail risk metrics.

The performance of each of the algorithms has been demonstrated through two
potential applications on both simulated and real data. The first is the
estimation of TRMs for the joint simulation algorithm, and the second one is
the inference of quantities involving some conditional tail distribution for
the conditional simulation algorithm.

The aim of this study was to enhance the quantity of observations in extreme
regions. Following the conclusion of the parametric numerical experiment, it
was established that even empirical methods applied to data simulated with our
algorithms can yield accurate risk measurement estimates. Consequently, future
works could focus on the combination of our simulation approaches with more
sophisticated estimation techniques than empirical ones (e.g. Padoan et al.,
2023; Davison et al., 2023).

#### Codes

All the codes and data are publicly available at
https://github.com/MadharNisrine/MultivariateExtremeSimulator.git.

##  7 Theoretical support

###  7.1 Theoretical support for Algorithm 1

The following Proposition 7.1 guarantees that the simulated samples through
Algorithm 1 are actually distributed according to a standard MGP distribution,
for all $q\in\\{1,\dots,d\\}$, and in particular for $q=1$.

######  Proposition 7.1.

For $j\in\\{1,\ldots,d\\}$, let $F$ be the distribution function of
$\bm{\Delta}^{(j)}_{1},\ldots,\bm{\Delta}^{(j)}_{n}$, and $\widetilde{F}_{nm}$
the empirical distribution function of
$\widetilde{\bm{\Delta}}^{(j)}_{1},\ldots,\widetilde{\bm{\Delta}}^{(j)}_{m}$.
If $\widetilde{F}_{nm}$ converges in distribution to $F$, as $n$ and $m$ tend
to infinity, then $\left(\widetilde{\bm{Z}}_{\ell}\right)_{1\leq\ell\leq m}$
converges in distribution to a standard MGP vector distributed as the sample
$\left(\bm{Z}_{i}\right)_{1\leq i\leq n}$.

The two following technical lemmas are useful to prove Propostion 7.1. Their
proofs are straightforward.

######  Lemma 7.2.

$\displaystyle\min_{j}\left(z_{j}-\mathcal{D}_{j}\right)=\min_{j}\left(z_{1}+\Delta^{q,1},z_{2}+\Delta^{q,2},\ldots,z_{d}+\Delta^{q,d}\right)-\min_{j}\left(\Delta^{q,1},\Delta^{q,2},\ldots,\Delta^{q,d}\right),$

where $\mathcal{D}_{j}=\sum_{k=1,k\neq j}^{d}\Delta^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\Delta^{l,k}<0}$.

######  Lemma 7.3.

$\displaystyle\min_{j}\left(z_{1}+\Delta^{q,1},z_{2}+\Delta^{q,2},\ldots,z_{d}+\Delta^{q,d}\right)-\min_{k}\left(\Delta^{q,1},\Delta^{q,2},\ldots,\Delta^{q,d}\right)$
$\displaystyle=-\max_{j}\left(T_{1}-z_{1},T_{2}-z_{2},\ldots,T_{d}-z_{d}\right)+\max_{j}\left(T_{1},T_{2},\ldots,T_{d}\right)$

###### Proof.

Proposition 7.1 and Algorithm 1 are written with $q=1$, but since the proof
does not change with the value of $q$, we prove the lemma with a generic value
$q=1,\ldots,d$. Let $\widetilde{\bm{Z}}=(Z_{1},\ldots,Z_{d})$ be a
d-dimensional vector resulting from Algorithm 1. We wish to prove that
$\widetilde{\bm{Z}}$ converges in distribution towards a standard MGP
distribution. Then, $\widetilde{\bm{Z}}$ satisfies Equation (8)

$\widetilde{Z}_{j}=E+\sum_{k=1,k\neq
j}^{d}\widetilde{\Delta}^{j,k}\prod_{\ell=1,\ell\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{\ell,k}<0},\mbox{ for all
}j=1,\ldots,d,$

where $E$ is a standard exponential distribution and
$\widetilde{\Delta}^{j,k}=\widetilde{Z}_{j}-\widetilde{Z}_{k}$, for all
$j,k=1,\ldots,d$. Recall also the notation : for a fixed
$q\in\\{1,\dots,d\\}$, the $q$-th vector of differences is defined as follows

$\widetilde{\bm{\Delta}}^{(q)}=\left(\widetilde{\Delta}^{q,1},\dots,\widetilde{\Delta}^{q,d}\right)\in\mathbb{R}^{d}.$

Then, the multivariate distribution of $\widetilde{\bm{Z}}$ is given by

$\displaystyle\mathbb{P}\left(\widetilde{\bm{Z}}\leq\bm{z}\right)$
$\displaystyle=\mathbb{P}\left(E\leq z_{1}-\sum_{k=1,k\neq
1}^{d}\widetilde{\Delta}^{1,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{l,k}\leq 0},\ldots,E\leq
z_{d}-\sum_{k=1,k\neq d}^{d}\widetilde{\Delta}^{d,j}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{l,k}\leq 0}\right)$
$\displaystyle=\mathbb{P}\left(E\leq\min_{j}\left\\{z_{j}-\sum_{k=1,k\neq
j}^{d}\widetilde{\Delta}^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{l,k}\leq 0}\right\\}\right)$
$\displaystyle=\int\mathbb{P}\left(\left.E\leq\min_{j}\left\\{z_{j}-\sum_{k=1,k\neq
j}^{d}\widetilde{\Delta}^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{l,k}\leq
0}\right\\}\right|\widetilde{\bm{\Delta}}^{(q)}=\widetilde{\bm{\delta}}^{(q)}\right)f_{\widetilde{\bm{\Delta}}^{(q)}}(\widetilde{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\widetilde{\bm{\delta}}^{(q)}$
$\displaystyle=\int\mathbb{P}\left(E\leq\min_{j}\left\\{z_{j}-\sum_{k=1,k\neq
j}^{d}\widetilde{\delta}^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\delta}^{l,k}\leq
0}\right\\}\right)f_{\widetilde{\bm{\Delta}}^{(q)}}(\widetilde{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\widetilde{\bm{\delta}}^{(q)}$
$\displaystyle=\int\left[1-\min\left(1,\exp\left[-\min_{j}\left\\{z_{j}-\sum_{k=1,k\neq
j}^{d}\widetilde{\delta}^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\delta}^{l,k}\leq
0}\right\\}\right]\right)\right]f_{\widetilde{\bm{\Delta}}^{(q)}}(\widetilde{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\widetilde{\bm{\delta}}^{(q)}$
$\displaystyle=1-\mathbb{E}_{\widetilde{\bm{\Delta}}^{(q)}}\left[\min\left(1,\exp\left[-\min_{j}\left\\{z_{j}-\sum_{k=1,k\neq
j}^{d}\widetilde{\Delta}^{j,k}\prod_{l=1,l\neq
k}^{d}\mathbf{1}_{\widetilde{\Delta}^{l,k}\leq
0}\right\\}\right]\right)\right]$
$\displaystyle=1-\mathbb{E}_{\widetilde{\bm{\Delta}}^{(q)}}\left[\min\left(1,e^{-\min_{j}\left(z_{1}+\widetilde{\Delta}^{q,1},z_{2}+\widetilde{\Delta}^{q,2},\ldots,z_{d}+\widetilde{\Delta}^{q,d}\right)+\min_{j}\left(\widetilde{\Delta}^{q,1},\widetilde{\Delta}^{q,2},\ldots,\widetilde{\Delta}^{q,d}\right)}\right)\right]$

Last equation arises from Lemma 7.2. Then, assuming that the Portemanteau
theorem holds in dimension $d>2$,

$\displaystyle\mathbb{P}\left(\left.\widetilde{\bm{Z}}<\bm{z}\right\lvert\bm{\Delta}_{1}^{(q)},\ldots,\bm{\Delta}_{n}^{(q)}\right)$
$\displaystyle\rightarrow
1-\mathbb{E}_{\bm{\Delta}^{(q)}}\left[\min\left(1,e^{-\min_{j}\left(z_{1}+\Delta^{q,1},z_{2}+\Delta^{q,2},\ldots,z_{d}+\Delta^{q,d}\right)+\min_{j}\left(\Delta^{q,1},\Delta^{q,2},\ldots,\Delta^{q,d}\right)}\right)\right]$

Now, noting that the term in the exponent could be written alternatively as
described in Lemma 7.3 injecting this result into the exponent term gives

$\displaystyle\mathbb{P}\left(\left.\widetilde{\bm{Z}}<\bm{z}\right\lvert\bm{\Delta}_{1}^{(q)},\ldots,\bm{\Delta}_{n}^{(q)}\right)$
$\displaystyle\rightarrow
1-\mathbb{E}_{\bm{\Delta}^{(q)}}\left[\min\left(1,e^{\max_{j}\left(T_{1}-z_{1},T_{2}-z_{2},\ldots,T_{d}-z_{d}\right)-\max_{j}\left(T_{1},T_{2},\ldots,T_{d}\right)}\right)\right]$
$\displaystyle=1-\mathbb{E}_{\bm{\Delta}^{(q)}}\left[\min\left(1,e^{\max\left(\bm{T}-\bm{z}\right)-\max\left(\bm{T}\right)}\right)\right],$

hence, one recovers the cumulative distribution function of a multivariate GP
distributions as defined in (Rootzén et al., 2018). ∎

###  7.2 Theoretical support for Algorithm 2

Proposition 7.4 derives the density of the conditional distribution of
$\Delta^{q,j}$ given $\bm{Z}_{-j}=\bm{z}_{-j}$ needed for Algorithm 2.

######  Proposition 7.4.

For $q\in\\{1,\ldots,d\\}$ with $q\neq j$, let $z_{\star}=\max\bm{z}_{-j}$.
The conditional distribution of $\Delta^{q,j}$ given $\bm{Z}_{-j}=\bm{z}_{-j}$
is given by

  1. 1.

If $z_{\star}>0$ and $z_{\star}=z_{q}$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(0)+I_{2}(0)}\left(\mathbf{1}_{\delta^{q,j}>0}+\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta^{q,j}\leq
0}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right).$

  2. 2.

If $z_{\star}>0$ and $z_{\star}\neq z_{q}$, denoting $\delta_{\star}\coloneqq
z_{q}-z_{\star}$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q}-z_{\star})+I_{2}(z_{q}-z_{\star})}\left(\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta_{\star}>\delta^{q,j}}+\mathrm{e}^{\delta_{\star}}\mathbf{1}_{\delta_{\star}\leq\delta^{q,j}}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right).$

  3. 3.

If $z_{\star}\leq 0$,

$f_{\Delta^{q,j}\mid\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q})}\mathrm{e}^{\delta^{q,j}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{\delta^{q,j}<z_{q}}.$

where

  * •

$\bm{\delta}^{(q)}=(\delta^{(q,1)},\dots,\delta^{(q,d)})$

  * •

for any $x\in\mathbb{R}$,

$I_{1}(x)\coloneqq\int_{-\infty}^{x}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\quad\text{
and }\quad
I_{2}(x)\coloneqq\mathrm{e}^{x}\int_{x}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta,$

where $\bar{\bm{\delta}}^{(q)}$ has the same components as $\bm{\delta}^{(q)}$
except for $\delta^{q,j}$ which is replaced by $\zeta=z_{q}-z_{j}$.

###### Proof.

The conditional distribution of $Z_{j}$ given $\bm{Z}_{-j}=\bm{z}_{-j}$ is
given by

$f_{Z_{j}|\bm{Z}_{-j}=\bm{z}_{-j}}(z_{j})=f_{\bm{Z}}(\bm{z})/f_{\bm{Z}_{-j}}(\bm{z}_{-j}),$

where the joint distribution function of $\bm{Z}$ is given by

$\displaystyle
f_{\bm{Z}}(\bm{z})=\mathrm{e}^{-\max\\{\bm{z}\\}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{\bm{z}\nleqslant
0}$

for a fixed $q\in\\{1,\dots,d\\}$ such that $q\neq j$.

We distinguish between two cases depending on the sign of
$z_{\star}:=\max\bm{z}_{-j}$.

#### If $z_{\star}>0$.

First, the marginal distribution of $\bm{Z}_{-j}$ is obtained by

$\displaystyle f_{\bm{Z}_{-j}}(\bm{z}_{-j})$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{e}^{-\max\\{\bm{z}\\}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathop{}\\!\mathrm{d}z_{j}$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{e}^{-\max\\{z_{1}-z_{q}+\zeta,z_{2}-z_{q}+\zeta,\ldots,z_{d}-z_{1}+\zeta\\}-z_{q}+\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta,$

where $\bar{\bm{\delta}}^{(q)}$ has the same components as $\bm{\delta}^{(q)}$
except for $\delta^{q,j}$ which is replaced by $\zeta=z_{q}-z_{j}$. The last
integral is split in two integrals of $f_{\bm{Z}}(\bm{z})$ on $\mathbb{R}_{-}$
and $\mathbb{R}_{+}$, respectively. We now compute each integral separately.

For the integral $\mathbb{R}_{+}$, let us note that when
$\zeta=z_{q}-z_{j}>0$, then $z_{\star}>z_{j}$, thus
${\max\\{z_{1}-z_{q}+\zeta,z_{2}-z_{q}+\zeta,\ldots,z_{d}-z_{q}+\zeta\\}+z_{q}-\zeta=\max(\bm{z})=z_{\star}}$,
which yields

$\int_{0}^{\infty}\mathrm{e}^{-\max\\{z_{1}-z_{q}+\zeta,z_{2}-z_{q}+\zeta,\ldots,z_{d}-z_{q}+\zeta\\}-z_{q}+\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta=\mathrm{e}^{-z_{\star}}\int_{0}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta.$

For the integral $\mathbb{R}_{-}$, we split the integral on the intervals
$(-\infty,z_{q}-z_{\star}]$ and $[z_{q}-z_{\star},0]$. On the interval
$(-\infty,z_{q}-z_{\star}]$, $z_{j}\geq z_{\star}$, so that
$\max(\bm{z})=z_{j}$ and on the interval $[z_{q}-z_{\star},0]$, $z_{j}\leq
z_{\star}$, so that $\max(\bm{z})=z_{\star}$. Thus, we obtain

$\displaystyle\int_{-\infty}^{0}\mathrm{e}^{-\max\\{z_{1}-z_{q}+\zeta,z_{2}-z_{q}+\zeta,\ldots,z_{d}-z_{q}+\zeta\\}-z_{q}+\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta$
$\displaystyle=$
$\displaystyle\left(\mathrm{e}^{-z_{q}}\int_{-\infty}^{z_{q}-z_{\star}}e^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta+e^{-z_{\star}}\int_{z_{q}-z_{\star}}^{0}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\right)\mathbf{1}_{z_{\star}\neq
z_{q}}$
$\displaystyle+\left(e^{-z_{q}}\int_{-\infty}^{0}e^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\right)\mathbf{1}_{z_{\star}=z_{q}}\,.$

Hence, the marginal distribution of $\bm{Z}_{-j}$ is equal to

  1. 1.

If $z_{\star}>0$ and $z_{\star}=z_{q}$

$f_{\bm{Z}_{-j}}(\bm{z}_{-j})=\mathrm{e}^{-z_{q}}\left(\int_{-\infty}^{0}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta+\int_{0}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\right)=\mathrm{e}^{-z_{q}}\left(I_{1}(0)+I_{2}(0)\right)$

  2. 2.

If $z_{\star}>0$ and $z_{\star}\neq z_{q}$

$\displaystyle f_{\bm{Z}_{-j}}(\bm{z}_{-j})$
$\displaystyle=\mathrm{e}^{-z_{q}}\left(\int_{-\infty}^{z_{q}-z_{\star}}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta+\mathrm{e}^{z_{q}-z_{\star}}\int_{z_{q}-z_{\star}}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\right)$
$\displaystyle=\mathrm{e}^{-z_{q}}\left(I_{1}(z_{q}-z_{\star})+I_{2}(z_{q}-z_{\star})\right)$

where we have denoted, for any $x\in\mathbb{R}$

$I_{1}(x)=\int_{-\infty}^{x}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\quad\text{
and }\quad
I_{2}(x)=\mathrm{e}^{x}\int_{x}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta$

The conditional distribution of $Z_{j}$ given $\bm{Z}_{-j}=\bm{z}_{-j}$ is
then obtained as

  1. 1.

If $z_{\star}>0$ and $z_{\star}=z_{q}$

$f_{Z_{j}|\bm{Z}_{-j}=\bm{z}_{-j}}(z_{j})=\frac{1}{I_{1}(0)+I_{2}(0)}\left(\mathbf{1}_{z_{q}>z_{j}}+\mathrm{e}^{z_{q}-z_{j}}\mathbf{1}_{z_{q}\leq
z_{j}}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)$

  2. 2.

If $z_{\star}>0$ and $z_{\star}\neq z_{q}$

$f_{Z_{j}|\bm{Z}_{-j}=\bm{z}_{-j}}(z_{j})=\frac{1}{I_{1}(z_{q}-z_{\star})+I_{2}(z_{q}-z_{\star})}\left(\mathrm{e}^{(z_{q}-z_{j})}\mathbf{1}_{z_{j}>z_{\star}}+\mathrm{e}^{(z_{q}-z_{\star})}\mathbf{1}_{z_{j}\leq
z_{\star}}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)$

Yielding,

  1. 1.

If $z_{\star}>0$ and $z_{\star}=z_{q}$

$f_{\Delta^{q,j}|\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(0)+I_{2}(0)}\left(\mathbf{1}_{\delta^{q,j}>0}+\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta^{q,j}\leq
0}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)$

  2. 2.

If $z_{\star}>0$ and $z_{\star}\neq z_{q}$

$f_{\Delta^{q,j}|\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q}-z_{\star})+I_{2}(z_{q}-z_{\star})}\left(\mathrm{e}^{\delta^{q,j}}\mathbf{1}_{\delta_{\star}>\delta^{q,j}}+\mathrm{e}^{\delta_{\star}}\mathbf{1}_{\delta_{\star}\leq\delta^{q,j}}\right)f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)$

where $\delta_{\star}=z_{q}-z_{\star}$.

Note that while the conditional distribution does not depend on $\bm{z}_{-j}$
when $z_{\star}=z_{q}$, it does depend on $\bm{z}_{-j}$ when $z_{\star}\neq
z_{q}$.

#### If $z_{\star}\leq 0$ (Case 3).

The joint distribution is given by

$f_{\bm{Z}}(\bm{z})=e^{-\max\\{\bm{z}\\}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{z_{j}>0}.$

Thus, the marginal distribution of $\bm{Z}_{-j}$ is obtained via

$f_{\bm{Z}_{-j}}(\bm{z}_{-j})=\mathrm{e}^{-z_{q}}\int_{-\infty}^{z_{q}}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}({\bar{\bm{\delta}}^{(q)}})\mathop{}\\!\mathrm{d}\zeta,$

the conditional distribution is then given by

$f_{Z_{j}|\bm{Z}_{-j}=\bm{z}_{-j}}(z_{j})=\frac{1}{I_{1}(z_{q})}\mathrm{e}^{z_{q}-z_{j}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{z_{j}>0}\,.$

Finally,

$f_{\Delta^{q,j}|\bm{Z}_{-j}=\bm{z}_{-j}}(\delta^{q,j})=\frac{1}{I_{1}(z_{q})}\mathrm{e}^{\delta^{q,j}}f_{\bm{\Delta}^{(q)}}\left(\bm{\delta}^{(q)}\right)\mathbf{1}_{\delta_{j}<z_{q}}\,.$

∎

Proposition 7.5 guarantees that the use of the rejection sample in a
multivariate context.

######  Proposition 7.5.

Let $V$ be a random variable with density $f$ and let $\bm{Z}$ be a
$d$-dimensional random vector of density $g$ from which simulation of samples
is known. Assuming that there exits a function $R$ in $\mathbb{R}^{d-1}$ such
that $\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}<\infty$ and
$z_{j}$-independent. Then, for all $v\in\mathbb{R}$,

$\mathbb{P}\left(V<v\right)=\mathbb{P}\left(Z_{j}<v\lvert
U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\right)\,,$

where $U$ is a uniform random variable on $[0,1]$.

###### Proof.

Let $v\in\mathbb{R}$.

$\displaystyle\mathbb{P}\left(Z_{j}<v\mid
U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\right)$
$\displaystyle=\frac{\mathbb{P}\left(Z_{j}<v,U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\right)}{\mathbb{P}\left(U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\right)}$
$\displaystyle=\frac{\int\mathbb{P}\left(Z_{j}<v,U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\mid\bm{Z}=\bm{z}\right)g(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}}{\int\mathbb{P}\left(U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\mid\bm{Z}=\bm{z}\right)g(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}}$
$\displaystyle=\frac{\int\mathbf{1}_{z_{j}<v}\frac{f(z_{j})}{R(\bm{z}_{-j})g(\bm{z})}g(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}}{\int\frac{f(z_{j})}{R(\bm{z}_{-j})g(\bm{z})}g(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}}$
$\displaystyle=\frac{\int_{-\infty}^{v}\mathbf{1}_{z_{j}<v}f(z_{j})\left(\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}\right)\mathop{}\\!\mathrm{d}z_{j}}{\int
f(z_{j})\left(\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}\right)\mathop{}\\!\mathrm{d}z_{j}}\,.$

Assuming that
$\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}<\infty$ and
$z_{j}$-independent, it follows that

$\mathbb{P}\left(Z_{j}<v\lvert
U<\frac{f(Z_{j})}{R(\bm{Z}_{-j})g(\bm{Z})}\right)=\frac{\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}\int_{-\infty}^{v}f(z_{j})\mathop{}\\!\mathrm{d}z_{j}}{\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}\int
f(z_{j})\mathop{}\\!\mathrm{d}z_{j}}=\int_{-\infty}^{v}f(z_{j})\mathop{}\\!\mathrm{d}z_{j}=\mathbb{P}\left(V\leq
v\right)\,.$

∎

Now, we need to specify the rejection constant $R(\bm{z}_{-j})$ is our context
and verify that
${\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}<\infty}$ and that
it is $z_{j}$-independent. $R(\bm{z}_{-j})$ is defined as follows

$R(\bm{z}_{-j})=R(\bm{\delta}_{-j})=\sup_{\delta^{q,j}}\frac{f_{\Delta^{q,j}|\bm{Z}_{-j}=\bm{z}_{-j}}\left(\delta^{q,j}\right)}{f_{\bm{\Delta}^{(q)}}(\bm{\delta}^{(q)})}$

where $\delta^{q,j}=z_{q}-z_{j}$. Thanks to Proposition 7.4, we obtained that

$R(\bm{z}_{-j})=\begin{cases}I_{1}(0)+I_{2}(0)&\text{ if $z_{\star}>0$ and
$z_{\star}=z_{q}$}\\\ I_{1}(z_{1}+z_{\star})+I_{2}(z_{1}+z_{\star})&\text{ if
$z_{\star}>0$ and $z_{\star}\neq z_{q}$}\\\ I_{1}(z_{q})&\text{ if
$z_{\star}\leq 0$ and $z_{j}>0$}\end{cases}$ (11)

######  Proposition 7.6.

$R(\bm{z}_{-j})$ chosen as in Equation (11) verifies
${\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}<\infty}$ and that
it is $z_{j}$-independent.

###### Proof.

Consider the first case where $\bm{z}_{\star}=z_{q}>0$. Then,

$\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}=\int\left(\int_{-\infty}^{0}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta+\int_{0}^{\infty}f_{\bm{\Delta}^{(q)}}\left(\bar{\bm{\delta}}^{(q)}\right)\mathop{}\\!\mathrm{d}\zeta\right)\mathop{}\\!\mathrm{d}\bm{z}_{-j}\,.$

Knowing that $f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})$ is a probability
density function,

$\displaystyle\int_{0}^{\infty}f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\zeta<\int_{\mathbb{R}}f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\zeta<+\infty\,.$
(12)

Similarly ,

$\displaystyle\int_{-\infty}^{0}\mathrm{e}^{\zeta}f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\zeta$
$\displaystyle=\int_{0}^{+\infty}\mathrm{e}^{-\xi}f_{\bm{\Delta}^{(q)}}(\bar{\delta}^{q,1},\ldots,-\xi,\ldots,\bar{\delta}^{q,d})\mathop{}\\!\mathrm{d}\xi$
$\displaystyle<\int_{0}^{+\infty}f_{\bm{\Delta}^{(q)}}(\bar{\delta}^{q,1},\ldots,-\xi,\ldots,\bar{\delta}^{q,d})\mathop{}\\!\mathrm{d}\xi$
$\displaystyle<\int_{\mathbb{R}}f_{\bm{\Delta}^{(q)}}(\bar{\delta}^{q,1},\ldots,-\xi,\ldots,\bar{\delta}^{q,d})\mathop{}\\!\mathrm{d}\xi<+\infty$
(13)

Inequalities (12) and (7.2) show that $R(\bm{z}_{-j})$ is well defined and
finite. Now, we show that
${\int(R(\bm{z}_{-j}))^{-1}\mathop{}\\!\mathrm{d}\bm{z}_{-j}<\infty}$, which
can be shown through the application of this same argument iteratively on each
component of $\bm{z}_{-j}$, which gives

$\displaystyle\int_{\mathbb{R}}\ldots\int_{\mathbb{R}}\int_{0}^{+\infty}f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}\zeta\mathop{}\\!\mathrm{d}z_{d}\ldots\mathop{}\\!\mathrm{d}z_{1}<\int_{\mathbb{R}}\cdots\int_{\mathbb{R}}\int_{\mathbb{R}}f_{\bm{\Delta}^{(q)}}(\bar{\bm{\delta}}^{(q)})\mathop{}\\!\mathrm{d}u\mathop{}\\!\mathrm{d}z_{d}\ldots\mathop{}\\!\mathrm{d}z_{1}=1.$

These results held true for the remaining cases. ∎

## References

  * Acharya et al. (2017) Acharya, V. V., L. H. Pedersen, T. Philippon, and M. Richardson (2017).  Measuring systemic risk.  The review of financial studies 30(1), 2–47.
  * Artzner et al. (1999) Artzner, P., F. Delbaen, J.-M. Eber, and D. Heath (1999).  Coherent measures of risk.  Mathematical finance 9(3), 203–228.
  * Basel Committee on Banking Supervision (2013) Basel Committee on Banking Supervision (2013).  Fundamental review of The Trading Book: A revised market risk framework.  Basel Committee on Banking Supervision.
  * Bellini and Di Bernardino (2017) Bellini, F. and E. Di Bernardino (2017).  Risk management with expectiles.  The European Journal of Finance 23(6), 487–506.
  * Cai et al. (2015) Cai, J.-J., J. H. Einmahl, L. Haan, and C. Zhou (2015).  Estimation of the marginal expected shortfall: the mean when a related variable is extreme.  Journal of the Royal Statistical Society Series B: Statistical Methodology 77(2), 417–442.
  * Coles (2001) Coles, S. (2001).  An introduction to statistical modeling of extreme values, Volume 208.  Springer.
  * Daouia et al. (2018) Daouia, A., S. Girard, and G. Stupfler (2018).  Estimation of tail risk based on extreme expectiles.  Journal of the Royal Statistical Society Series B: Statistical Methodology 80(2), 263–292.
  * Davison et al. (2023) Davison, A. C., S. A. Padoan, and G. Stupfler (2023).  Tail risk inference via expectiles in heavy-tailed time series.  Journal of Business & Economic Statistics 41(3), 876–889.
  * Gnecco et al. (2022) Gnecco, N., E. M. Terefe, and S. Engelke (2022).  Extremal random forests.  arXiv preprint arXiv:2201.12865.
  * Goegebeur et al. (2024) Goegebeur, Y., A. Guillou, and J. Qin (2024).  Dependent conditional tail expectation for extreme levels.  Stochastic Processes and their Applications, 104330.
  * Gründemann et al. (2023) Gründemann, G. J., E. Zorzetto, H. E. Beck, M. Schleiss, N. Van de Giesen, M. Marani, and R. J. van der Ent (2023).  Extreme precipitation return levels for multiple durations on a global scale.  Journal of Hydrology 621, 129558.
  * Josaphat and Syuhada (2021) Josaphat, B. P. and K. Syuhada (2021).  Dependent conditional value-at-risk for aggregate risk models.  Heliyon 7(7), e07492.
  * Kiriliouk et al. (2019) Kiriliouk, A., H. Rootzén, J. Segers, and J. L. Wadsworth (2019).  Peaks over thresholds modeling with multivariate generalized pareto distributions.  Technometrics 61(1), 123–135.
  * Legrand et al. (2023) Legrand, J., P. Ailliot, P. Naveau, and N. Raillard (2023).  Joint stochastic simulation of extreme coastal and offshore significant wave heights.  The Annals of Applied Statistics 17(4), 3363 – 3383.
  * McNeil and Frey (2000) McNeil, A. J. and R. Frey (2000).  Estimation of tail-related risk measures for heteroscedastic financial time series: an extreme value approach.  Journal of empirical finance 7(3-4), 271–300.
  * McNeil et al. (2015) McNeil, A. J., R. Frey, and P. Embrechts (2015).  Quantitative risk management: concepts, techniques and tools-revised edition.  Princeton university press.
  * Nadarajah et al. (2014) Nadarajah, S., B. Zhang, and S. Chan (2014).  Estimation methods for expected shortfall.  Quantitative Finance 14(2), 271–291.
  * Nelsen (2006) Nelsen, R. B. (2006).  An introduction to copulas.  Springer.
  * Padoan et al. (2023) Padoan, S. A., S. Rizzelli, and M. Schiavone (2023).  Marginal expected shortfall inference under multivariate regular variation.  arXiv preprint arXiv:2304.07578.
  * Pasche and Engelke (2022) Pasche, O. C. and S. Engelke (2022).  Neural networks for extreme quantile regression with an application to forecasting of flood risk.  arXiv preprint arXiv:2208.07590.
  * Resnick (2007) Resnick, S. I. (2007).  Heavy-tail phenomena: probabilistic and statistical modeling.  Springer Science & Business Media.
  * Rockafellar et al. (2000) Rockafellar, R. T., S. Uryasev, et al. (2000).  Optimization of conditional value-at-risk.  Journal of risk, 21–42.
  * Rootzén et al. (2018) Rootzén, H., J. Segers, and J. L. Wadsworth (2018).  Multivariate generalized pareto distributions: Parametrizations, representations, and properties.  Journal of Multivariate Analysis 165, 117–131.
  * Rootzén and Tajvidi (2006) Rootzén, H. and N. Tajvidi (2006).  Multivariate generalized pareto distributions.  Bernoulli 12(5), 917–930.
  * Scarrott and MacDonald (2012) Scarrott, C. and A. MacDonald (2012).  A review of extreme value threshold estimation and uncertainty quantification.  REVSTAT-Statistical journal 10(1), 33–60.
  * Singh et al. (2013) Singh, A. K., D. E. Allen, and P. J. Robert (2013).  Extreme market risk and extreme value theory.  Mathematics and computers in simulation 94, 310–328.
  * Sklar (1959) Sklar, M. (1959).  Fonctions de répartition à n dimensions et leurs marges.  In Annales de l’ISUP, Volume 8, pp. 229–231.
  * Syuhada et al. (2022) Syuhada, K., O. Neswan, and B. P. Josaphat (2022).  Estimating copula-based extension of tail value-at-risk and its application in insurance claim.  Risks 10(6), 113.
  * Tasche (2002) Tasche, D. (2002).  Expected shortfall and beyond.  Journal of Banking & Finance 26(7), 1519–1533.
  * Yamai and Yoshiba (2005) Yamai, Y. and T. Yoshiba (2005).  Value-at-risk versus expected shortfall: A practical perspective.  Journal of Banking & Finance 29(4), 997–1015.