Modern extreme value statistics
for Utopian extremes

Jordan Richards1∗, Noura Alotaibi1, Daniela Cisneros1, Yan Gong1,3,
Matheus B. Guerrero1,2, Paolo Redondo1, Xuanjie Shao1

11footnotetext:  Statistics Program, Computer, Electrical and Mathematical
Sciences and Engineering (CEMSE) Division, King Abdullah University of Science
and Technology (KAUST), Thuwal, Saudi Arabia.22footnotetext:  Department of
Mathematics, California State University Fullerton, California,
USA.33footnotetext:  Department of Mathematics, Imperial College London,
London, UK. $\;{}^{*}$Corresponding author<EMAIL_ADDRESS>

Abstract

Capturing the extremal behaviour of data often requires bespoke marginal and
dependence models which are grounded in rigorous asymptotic theory, and hence
provide reliable extrapolation into the upper tails of the data-generating
distribution. We present a modern toolbox of four methodological frameworks,
motivated by modern extreme value theory, that can be used to accurately
estimate extreme exceedance probabilities or the corresponding level in either
a univariate or multivariate setting. Our frameworks were used to facilitate
the winning contribution of Team Yalla to the data competition organised for
the 13${}^{\text{th}}$ International Conference on Extreme Value Analysis
(EVA2023). This competition comprised seven teams competing across four
separate sub-challenges, with each requiring the modelling of data simulated
from known, yet highly complex, statistical distributions, and extrapolation
far beyond the range of the available samples in order to predict
probabilities of extreme events. Data were constructed to be representative of
real environmental data, sampled from the fantasy country of “Utopia”.

Keywords: additive models; conditional extremes; neural Bayes estimation; non-
parametric probability estimators; non-stationary extremal dependence;
quantile regression

##  1 Introduction

Extreme value analysis (EVA) concerns inference for events with rare
frequency, large magnitude, and high impact across a number of sectors. In the
financial sector, EVA helps to assess the risk of market crashes and extreme
losses, as demonstrated by Gong and Huser, (2022). In environmental sciences,
EVA is used to understand the characteristics of extreme weather events and
natural disasters, as explored by Shao et al., (2022) and Cisneros et al.,
2023a . In biostatistics, EVA helps to model outbreaks of rare diseases and
epidemics, as demonstrated by Zhang et al., (2023). In neuroscience, EVA
reveals hidden patterns in brain connectivity, with applications shown in
Guerrero et al., (2023) and Redondo et al., (2023). These diverse applications
highlight the versatility and importance of EVA in analysing and managing
extreme events across various scientific fields.

Mitigation of risk associated with extreme events first requires risk
quantification, and this typically proceeds through estimation of exceedance
probabilities or their corresponding exceedance level, i.e., a quantile or
return level. Events of interest are often so extreme that this level exceeds
the range of sampled data, and so EVA leverages asymptotically-justified
distributional models to extrapolate far into the tails of the distribution.
Classical univariate and multivariate extreme value models have been well-
studied in, for example, textbooks by Coles, (2001), Beirlant et al., (2006),
Haan and Ferreira, (2006), and Reiss et al., (2007), and reviews by Segers,
(2012) and Cooley et al., (2019). We exploit these models to provide a toolbox
of four frameworks that can be used to predict extreme exceedance
probabilities or levels for univariate and multivariate data, in both the
presence and absence of covariates. These frameworks combine classical EVA
methods with modern modelling techniques, including additive models (Chavez-
Demoulin and Davison,, 2005; Youngman,, 2019), deep learning-based inference
(Sainsbury-Dale et al., 2023b, ; Richards et al., 2023b, ), non-stationary
conditional extremal dependence models (Heffernan and Tawn,, 2004; Winter et
al.,, 2016), and non-parametric probability estimators (Krupskii and Joe,,
2019).

We are motivated by the challenge organised for the 13${}^{\text{th}}$
International Conference on Extreme Value Analysis (EVA2023), with full
details provided in the editorial by Rohrbeck et al., (2023). In this paper,
we detail the four modelling frameworks used by “Team Yalla” to win the
aforementioned challenge, which comprised four sub-challenges that each
required prediction of exceedance probabilities or quantiles for data
simulated from highly complex statistical models. Whilst the data considered
in this work are simulated, the data were constructed to be representative of
real observations of an environmental process (sampled from the fantasy
country of “Utopia”), and so exhibit realistic characteristics, such as
sparsity, non-stationarity, and missingness. Moreover, as the true values of
the challenge predictands are known, our predictions can be easily validated.
Hence, we expect our proposed frameworks to perform well in practice, with
real data.

The novelty of this work is threefold: i) we propose an amortised neural Bayes
estimator for univariate quantiles; ii) we generalise the non-parametric
multivariate exceedance probability estimator of Krupskii and Joe, (2019); and
iii) we illustrate the efficacy of the extreme value regression models
proposed by Winter et al., (2016) and Youngman, (2019) when applied to
simulated data. The remainder of the paper is organised as follows: Section 2
details the data provided by Rohrbeck et al., (2023) and the corresponding
prediction challenge. Section 3 describes the methodology we adopt to address
the challenge, with Sections 3.1–3.2 and 3.3–3.5 focusing on univariate and
multivariate extremes, respectively. In Section 4, we apply our proposed
methodology to the data. Section 5 provides a concise conclusion and suggests
avenues for further work.

##  2 Data

We now describe the data available for inference, the prediction challenge
associated with each dataset, and the methodology used for assessing
predictions. For full details on the data, its construction, and the
subsequent prediction challenge, see Rohrbeck et al., (2023). We also note
that Rohrbeck et al., (2023) detail the true values for all predictands.

Sections 2.1 and 2.2 detail the available univariate and multivariate data,
respectively. Throughout, we use the notation $Y$ and
$\mathbf{Y}:=(Y_{1},\dots,Y_{d})^{\prime}$ to denote a random response
variable and random $d$-vector of response variables, respectively. Covariates
are denoted by $\mathbf{X}:=(X_{1},\dots,X_{l})^{\prime}\in\mathbb{R}^{l}$ for
$l\in\mathbb{N}$. Where appropriate, we use the subscript
$t\in\\{1,\dots,n\\}$ to denote temporal replicates of the response (i.e.,
$Y_{t}$ or $\mathbf{Y}_{t}:=(Y_{1,t},\dots,Y_{d,t})^{\prime}$) or covariates
(i.e., $\mathbf{X}_{t}:=(X_{1,t},\dots,X_{l,t})^{\prime}$) and lowercase
notation to denote observations. The values of $l$, $d$, and sample size ${n}$
differ throughout the paper.

###  2.1 Univariate setting

In the univariate setting, $d=1$ and $n=21000$. Covariate vector $\mathbf{X}$
comprises $l=8$ variables: wind direction, wind speed, atmosphere, season, and
four unnamed variables (denoted $V_{1},\dots,V_{4})$. All variables are real-
valued except season, which is Boolean. Atmosphere and season exhibit a
cyclical pattern, and are known to repeat values every $21000$ and $300$
observations, respectively. The distribution of the covariates is known to be
stationary with respect to time $t$, and $\\{Y_{t}:t=1,\dots,n\\}$ are known
to be independent conditional on the covariates. Note that $11.7\%$ of
observations $\\{\mathbf{x}_{t}:t=1,\dots,n\\}$ have at least one missing
value amongst all covariates but atmosphere and season; this value is known to
be missing at random.

Sub-challenges C1 and C2 pertain to these univariate data. Sub-challenge C1
requires the prediction of a $50\%$ confidence interval for the $q-$quantile
of $Y\mid(\mathbf{X}=\mathbf{x}^{*}_{j})$ for 100 test covariate sets,
$\\{\mathbf{x}^{*}_{j}:j=1,\dots,100\\},$ and for $q=0.9999$ corresponding to
an extreme quantile. Predictions are assessed based on the actual coverage of
the $100$ produced confidence intervals, that is, how close they are to having
a $50\%$ coverage of the true quantile values.

Sub-challenge C2 requires estimation of the marginal $q$-quantile of $Y$ for
$q=1-(6\times 10^{4})^{-1}$, i.e., estimating $\theta$ such that
$\Pr\\{Y>\theta\\}=1-q$. However, inference for $\theta$ should seek to
minimise the conservative asymmetric loss function

$L(\theta,\hat{\theta})=\begin{cases}0.9(0.99\theta-\hat{\theta}),\quad&\text{if
}0.99\theta>\hat{\theta},\\\ 0,\quad&\text{if }|\theta-\hat{\theta}|\leq
0.01\theta,\\\ 0.1(\hat{\theta}-1.01\theta),\quad&\text{if
}1.01\theta<\hat{\theta},\end{cases}$ (1)

where $\hat{\theta}$ denotes estimates of $\theta$ and $|\cdot|$ denotes the
absolute value; this loss function is illustrated in Figure 2. The loss
function in Eq. (1) provides a larger penalty for under-estimates of $\theta$,
relative to over-estimates, to encourage conservative estimates of the
quantile. Assessment of $\hat{\theta}$ is determined by smaller values of
$L(\theta,\hat{\theta})$, unless $L(\theta,\hat{\theta})=0$; in which case,
$\hat{\theta}$ with smaller values of $|\theta-\hat{\theta}|$ are preferred.

The methodologies for sub-challenges C1 and C2 are described in Section 3.1
and 3.2, respectively, with results for both sub-challenges provided in
Section 4.1.

###  2.2 Multivariate setting

Sub-challenges C3 and C4 concern two multivariate datasets. For C3 and C4, we
have $d=3$ and $d=50$, respectively, and $n=21000$ and $n=10000$,
respectively. For both datasets, the response vector $\mathbf{Y}$ is known to
have all standard Gumbel marginals. The three response variables for C3,
$(Y_{1},Y_{2},Y_{3})$, are accompanied by the atmosphere and season covariates
described in Section 2.1; hence, $l=2$. No covariates accompany the 50
response variables that compose the data for C4 (i.e., $l=0$).

Both sub-challenges require estimation of exceedance probabilities with low
order of magnitude. Sub-challenge C3 requires estimation of

$\displaystyle p_{1}$ $\displaystyle:=\Pr(Y_{1}>6,Y_{2}>6,Y_{3}>6),$ (2)
$\displaystyle p_{2}$ $\displaystyle:=\Pr(Y_{1}>7,Y_{2}>7,Y_{3}<m),$

where $m=-\log{(\log{2})}$ is the median of the standard Gumbel distribution.
Estimates $(\hat{p}_{1},\hat{p_{2}})$ of $(p_{1},p_{2})$ are assessed using a
probability-based scoring rule (see Smith,, 1999), defined by
$\sum\nolimits^{2}_{j=1}|\;p_{j}\log(p_{j}/\hat{p}_{j})+(1-p_{j})\log[(1-p_{j})/(1-\hat{p}_{j})]\;|$,
with smaller values desired.

Sub-challenge C4 requires estimation of two exceedance probabilities of the
form

$\displaystyle p_{3}$ $\displaystyle:=\Pr(Y_{1}>s_{1},\dots,Y_{50}>s_{1}),$
$\displaystyle p_{4}$
$\displaystyle:=\Pr(Y_{1}>s_{1},\dots,Y_{25}>s_{1},Y_{26}>s_{2},\dots,Y_{50}>s_{2}),$
(3)

where $s_{j},j=1,2,$ is the $(1-\phi_{j})$-quantile of the standard Gumbel
distribution with $\phi_{1}=1/300$ and $\phi_{2}=12\phi_{1}$. Hence, $p_{3}$
corresponds to an exceedance probability with all components of $\mathbf{Y}$
being equally extreme, i.e., concurrently exceeding the same marginal
quantile, whilst, for $p_{4}$, the first 25 components of $\mathbf{Y}$ exceed
a higher quantile, i.e., are more extreme, than the latter 25 components.
Assessment of estimates $(\hat{p}_{3},\hat{p}_{4})$ is identical as for
$(\hat{p}_{1},\hat{p}_{2})$. The methodologies for sub-challenges C3 and C4
are described in Sections 3.3 and 3.4–3.5, respectively, with results for both
sub-challenges provided in Section 4.2.

##  3 Methodology

In this section, we present the methodological details of our approaches. To
be specific, Sections 3.1 and 3.2 detail methodology for estimating univariate
extreme quantiles, whilst Sections 3.3-3.5 concern methodology for modelling
extremal dependence. Section 3.1 describes peaks-over-threshold modelling
using the generalised Pareto distribution and generalised additve models.
Section 3.2 describes a likelihood-free neural Bayes estimator for point
estimation of extreme quantiles. Section 3.3 describes a non-stationary model
for conditional extremal dependence, whilst Section 3.4 details a bivariate
extremal dependence measure. Section 3.5 concludes with details of a non-
parametric estimator for tail probabilities.

###  3.1 Peaks-over-threshold models

The upper tails of the distribution of a random variable can be modelled in a
peaks-over-threshold framework using the generalised Pareto distribution
(GPD), see, for example, Pickands, (1975), Davison and Smith, (1990), and
Coles, (2001). For a random variable $Y$, we assume that there exists some
high threshold $u$ such that the distribution of exceedances $(Y-u)\mid(Y>u)$
can be characterised by the GPD, denoted GPD$(\sigma_{u},\xi)$,
$\sigma_{u}>0,\xi\in\mathbb{R}$, with distribution function

$\displaystyle H(y)=\begin{cases}1-(1+\xi y/\sigma_{u})^{-1/\xi},\quad&\xi\neq
0,\\\ 1-\exp(-y/\sigma_{u}),\quad&\xi=0,\end{cases}$ (4)

where $y\geq 0$ for $\xi\geq 0$ and $0\leq y\leq-\sigma_{u}/\xi$ for $\xi<0$.

For regression, we model the conditional distribution
$(Y-u(\mathbf{x}))\mid(Y>u(\mathbf{x}),\mathbf{X}=\mathbf{x})$ as
$\mbox{GPD}(\sigma_{u}(\mathbf{x}),\xi(\mathbf{x}))$, where the exceedance
threshold and GPD scale and shape parameters are functions of covariates. We
follow Chavez-Demoulin and Davison, (2005); Youngman, (2019); Richards et al.,
2023a  and use a generalised additive model (GAM) representation for the
distribution, with $\sigma_{u}(\mathbf{x})$ and $\xi(\mathbf{x})$ modelled via
a basis of splines. The threshold $u(\mathbf{x})$ is taken to be the
intermediate $\lambda$-quantile of $Y\mid(\mathbf{X}=\mathbf{x})$ for
$\lambda<0.9999$ and this is modelled using additive quantile regression
(Fasiolo et al.,, 2021).

###  3.2 Neural point estimation

Neural point estimators, that is, neural networks that are trained to map
input data to parameter point estimates, have shown recent success as a
likelihood-free inference approach for statistical models. Although they have
been typically used for inference with classical spatial processes (see, e.g.,
Cremanns and Roos,, 2017; Zammit-Mangion and Wikle,, 2020; Gerber and Nychka,,
2021) and spatial extremal processes (see, e.g., Lenzi et al.,, 2023; Lenzi
and Rue,, 2023; Sainsbury-Dale et al., 2023a, ; Sainsbury-Dale et al., 2023b,
; Richards et al., 2023b, ), they can be exploited in a univariate setting.
For example, Rai et al., (2023) use a neural point estimator to make inference
with the univariate generalised extreme value distribution (see Coles,, 2001).
We construct a neural point estimator to perform extreme single quantile
estimation for a random variable $Y$. In particular, we follow Sainsbury-Dale
et al., 2023b  and construct a neural Bayes estimator (NBE), which is a neural
point estimator that approximately minimises the Bayes risk.

Define a set of univariate probability distributions $\mathcal{P}$ on a sample
space, taken to be $\mathbb{R}$, which are parameterised by a parameter
${\theta}\in\mathbb{R}$ such that
$\mathcal{P}\equiv\\{P_{{\theta}}:{\theta}\in\Theta\\}$, where $\Theta$ is the
parameter space; then $\mathcal{P}$ defines a parametric statistical model
(see McCullagh,, 2002). Denote
$\mathbf{Y}\equiv({Y}_{1},\dots,{Y}_{m})^{\prime}$ as $m$ mutually independent
realisations from $P_{{\theta}}\in\mathcal{P}$. A point estimator
$\hat{{\theta}}(\cdot)$ for model $\mathcal{P}$ is any mapping from
$\mathbb{R}^{m}$ to $\Theta$, and the output of such an estimator, for a given
${\theta}$ and $\mathbf{Y}$, can be assessed using a non-negative loss
function $L({\theta},\hat{{\theta}}(\mathbf{Y}))$. The risk of this point
estimator, evaluated at ${\theta}$, $R({\theta},\hat{{\theta}}(\cdot)),$ is
the loss $L(\cdot,\cdot)$ averaged over all possible realisations of
$\mathbf{Y}$, that is,

$R({\theta},\hat{{\theta}}(\cdot))\equiv\int_{\mathbb{R}^{m}}L({\theta},\hat{{\theta}}(\mathbf{y}))f(\mathbf{y}\mid{\theta})\mathop{}\\!\mathrm{d}\mathbf{y},$
(5)

where $f(\mathbf{y}\mid{\theta})$ is the density function of the data. We
define the Bayes risk $r_{\pi}(\cdot)$ as the weighted average of Eq. (5) over
all ${\theta}\in\Theta$, with respect to some prior measure $\pi(\cdot)$, as

$r_{\pi}(\hat{{\theta}}(\cdot))\equiv\int_{\Theta}R({\theta},\hat{{\theta}}(\cdot))\mathop{}\\!\mathrm{d}\pi({\theta}).$
(6)

If the estimator $\hat{{\theta}}(\cdot)$ minimises Eq. (6), we term it a Bayes
estimator with respect to the loss $L(\cdot,\cdot)$ and prior measure
$\pi(\cdot)$. Such estimators have attractive theoretical properties under
mild regularity conditions, for example, asymptotic normality, efficiency, and
consistency under the squared error loss (see, e.g., Lehmann and Casella,,
1998). Note that in the context of the prediction task detailed in Section 2,
the parameter ${\theta}$ is taken to be the $q$-quantile of the distribution
$P_{\theta}$ and the loss function is $L(\theta,\hat{\theta})$ in Eq. (1)
(illustrated in Figure 2). Details on the construction of the prior measure
$\pi(\cdot)$ and the statistical model $\mathcal{P}$ will follow.

A neural Bayes estimator (NBE) is a neural network designed to approximately
minimise Eq. (6). Such an estimator must be permutation-invariant, that is,
$\hat{{\theta}}(\mathbf{Y})=\hat{{\theta}}(\mathbf{Y}^{\star})$ where
$\mathbf{Y}^{*}$ denotes any permutation of the independent replicates in
$\mathbf{Y}$ (Sainsbury-Dale et al., 2023b, ). Sainsbury-Dale et al., 2023b
construct NBEs by leveraging the DeepSets neural network architecture (Zaheer
et al.,, 2017). Consider functions $\bm{\psi}:\mathbb{R}\mapsto\mathbb{R}^{q}$
and ${\phi}:\mathbb{R}^{q}\mapsto\mathbb{R}$, and a permutation-invariant set
function $\mathbf{a}:(\mathbb{R}^{q})^{m}\mapsto\mathbb{R}^{q}$, where the
$j$-th component of $\mathbf{a}$, $a_{j}(\cdot),$ returns the element-wise
average over its input set for $j=1,\dots,q$. We represent ${\phi}(\cdot)$ and
$\bm{\psi}(\cdot)$ as neural networks, and collect in
$\bm{\gamma}\equiv(\bm{\gamma}_{\phi}^{\prime},\bm{\gamma}_{\bm{\psi}}^{\prime})^{\prime}$
their estimable “weights” and “biases”. Our NBE is of the form

$\hat{{\theta}}(\mathbf{Y};\bm{\gamma})={\phi}(\mathbf{T}(\mathbf{Y};\bm{\gamma}_{{\phi}});\bm{\gamma}_{{\phi}}),\quad\text{
with
}\quad\mathbf{T}(\mathbf{Y};\bm{\gamma}_{\bm{\psi}})=\mathbf{a}(\\{\bm{\psi}(Y_{t};\bm{\gamma}_{\bm{\psi}}):t=1,\dots,m\\}).$
(7)

We use a densely-connected neural network to model both ${\phi}(\cdot)$ and
$\bm{\psi}(\cdot)$.

The NBE is built by obtaining neural network weights $\bm{\gamma}^{*}$ that
minimise the Bayes risk in the estimator space spanned by
$\hat{{\theta}}(\cdot,\bm{\gamma})$. As we cannot directly evaluate Eq. (6),
it is approximated using Monte Carlo methods. For a set of $K$ parameters in
the set $\vartheta$ drawn from the prior $\pi(\cdot)$, we simulate, for each
${\theta}\in\vartheta$, $J$ sets of $m$ mutually independent realisations from
$P_{{\theta}}$, and collect these in $\mathcal{Y}^{(\bm{\theta})}$. The Bayes
risk in Eq. (6) is then approximated by

$\hat{r}_{\pi}(\hat{{\theta}}(\cdot;\bm{\gamma}))=\frac{1}{K}\sum_{{\theta}\in\vartheta}\frac{1}{J}\sum_{\mathbf{y}\in\mathcal{Y}^{({\theta})}}L({\theta},\hat{{\theta}}(\mathbf{y};\bm{\gamma}))\approx{r}_{\pi}(\hat{{\theta}}(\cdot;\bm{\gamma})).$
(8)

We obtain estimates $\bm{\gamma}^{*}={\rm
argmin}_{\bm{\gamma}}\hat{r}_{\pi}(\hat{{\theta}}(\cdot;\bm{\gamma}))$ using
the package NeuralEstimators (Sainsbury-Dale et al., 2023b, ) in Julia
(Bezanson et al.,, 2017). To reduce computational expense while still
obtaining unbiased estimates of Eq. (6), we follow Richards et al., 2023b  by
setting $J$ in Eq. (8) to $J=1$.

To simulate prior $q$-quantiles $\theta\in\vartheta$ and the associated
training data $\mathbf{y}\in\mathcal{Y}^{(\bm{\theta})}$, we require a class
of reasonable parametric models for the upper tails of our data. To this end,
we exploit the model described in Section 3.1. Let
$(Y-u_{t})\mid(Y>u_{t},\mathbf{X}=\mathbf{x}_{t})\sim\mbox{GPD}(\sigma_{t},\xi_{t})$
for $t=1,\dots,m,$ and let $Y\mid(Y\leq u_{t},\mathbf{X}=\mathbf{x}_{t})\sim
F_{t}^{*}(\cdot)$. To construct a sample $\mathbf{y}$, we first draw $m$
parameter values, $\\{(u_{t},\sigma_{t},\xi_{t}):t=1,\dots,m\\},$ from a
suitable prior distribution, and then marginalise out these covariates to get
realisations of the unconditional $Y$ (see, e.g., Rohrbeck et al.,, 2018). As
there is no closed-form expression for the $q$-quantile $\theta$ of $Y$, we
compute it using Monte Carlo methods.

We begin by sampling $M>m$ values (uniformly, with replacement) from
$\\{(u_{t},\sigma_{t},\xi_{t}):t=1,\dots,m\\}$; for ease of exposition, we
denote this new sample $\\{(u_{j},\sigma_{j},\xi_{j}):j=1,\dots,M\\}$.
Conditional on the $M$ parameter values, we construct a sample of $M$
realisations of $Y$, denoted $\mathbf{y}^{*}:=(y^{*}_{1},\dots,y^{*}_{M})$. To
construct the sample $\mathbf{y}^{*}$, we take, for each $j=1,\dots,M$,
$y^{*}_{j}$ to be a draw from $u_{j}+\mbox{GPD}(\sigma_{j},\xi_{j})$ with
probability $(1-\lambda)$ and a draw from $F_{j}^{*}(\cdot)$, otherwise;
recall that the threshold $u_{j}$ is taken to be the $\lambda$-quantile of
$Y\mid(\mathbf{X}=\mathbf{x}_{j})$. As values of $Y<u_{j}$ are unlikely to
impact the high $q$-quantiles we seek to infer, we take $F_{j}^{*}(\cdot)$ to
be the empirical distribution of all observations
$\\{\tilde{y}_{t}:\tilde{y}_{t}\leq u_{t}\\}$. In our application,
$q=1-(6\times 10^{4})^{-1}$; provided that $M$ is sufficiently large, the
$q$-quantile $\theta$ of the associated distribution of $Y$ (unconditional on
$\mathbf{X}$) can be estimated empirically from $\mathbf{y}^{*}$. Finally, we
sub-sample $m$ replicates $\mathbf{y}:=(y_{1},\dots,y_{m})$ from
$\mathbf{y}^{*}$ and pair this with the quantile $\theta$; this pair provides
a single input in the NBE training data.

###  3.3 Conditional extremal models

Proposed by Heffernan and Tawn, (2004) and later generalised by Heffernan and
Resnick, (2007), the conditional extremes framework models the behaviour of a
random vector, conditional on one of its components being extreme. We adopt
this model as its inference is typically less computationally demanding than
that of other models for multivariate extremal dependence (Huser and
Wadsworth,, 2022). Additionally, it is capable of capturing both asymptotic
dependence and asymptotic independence in a parsimonious manner. To
accommodate non-stationarity with respect to covariates in the extremal
dependence structure of our data, we utilise the extension of the Heffernan
and Tawn, (2004) model proposed by Winter et al., (2016). This extension
represents the dependence parameters as a linear function (subject to some
non-linear link transformation) of the covariates.

Let the vector $\mathbf{Y}_{t}\equiv(Y_{1,t},Y_{2,t},Y_{3,t})^{\prime}$ for
$t\in\\{1,\dots,n\\}$ have standard Laplace margins (Keef et al.,, 2013). It
is noteworthy that our original data are known to possess standard Gumbel
marginals; therefore, we transform these to have standard Laplace margins.
Then, denote by $\mathbf{Y}_{-i,t}$ as the vector $\mathbf{Y}_{t}$ with its
$i$-th component removed. Note that all vector operations hereafter are taken
component-wise. Winter et al., (2016) assume that there exist vectors of
coefficients $\bm{\alpha}_{-i,t}:=\\{\alpha_{j\mid i,t}:j\in(1,2,3)\setminus
i\\}\in[-1,1]^{2}$ and $\bm{\beta}_{-i,t}:=\\{\beta_{j\mid
i,t}:j\in(1,2,3)\setminus i\\}\in[0,1)^{2}$ such that, for
$\mathbf{z}\in\mathbb{R}^{2}$ and $y>0$,

$\Pr\left\\{\frac{\mathbf{Y}_{-i,t}-\bm{\alpha}_{-i,t}Y_{i,t}}{Y_{i,t}^{\bm{\beta}_{-i,t}}}\leq\mathbf{z},Y_{i,t}-u>y\;\;\bigg{|}\;\;Y_{i,t}>u\right\\}\rightarrow
G_{-i,t}(\mathbf{z})\exp(-y),$ (9)

as $u\rightarrow\infty$ with $G_{-i,t}(\cdot)$ a non-degenerate bivariate
distribution function. The values of the dependence parameters
$\bm{\alpha}_{-i,t}$ and $\bm{\beta}_{-i,t}$ determine the strength and class
of extremal dependence exhibited between $Y_{i,t}$ and the corresponding
component of $\mathbf{Y}_{-i,t}$; for details, see Heffernan and Tawn, (2004).
We allow these parameters to vary with covariates $\mathbf{x}_{t}$ by letting

$\tanh^{-1}(\bm{\alpha}_{-i,t}):=\bm{\alpha}^{(0)}_{-i}+\bm{\alpha}^{(1)}_{-i}\mathbf{x}_{t},\quad\text{logit}(\bm{\beta}_{-i,t}):=\bm{\beta}^{(0)}_{-i}+\bm{\beta}^{(1)}_{-i}\mathbf{x}_{t},$
(10)

with coefficients
$\bm{\alpha}^{(0)}_{-i}$,$\;\bm{\beta}^{(0)}_{-i}\in\mathbb{R}^{2}$ and
$\bm{\alpha}^{(1)}_{-i}$,$\;\bm{\beta}^{(1)}_{-i}\in\mathbb{R}^{2\times l},$
where $l$ is the number of covariates. Note that the $\tanh(\cdot)$ and
$\text{logistic}(\cdot)$ link functions are used to ensure that the parameter
values are constrained to their correct ranges.

Modelling follows under the assumption that the limit in Eq. (9) holds in
equality for all $Y_{i,t}>u$ for some sufficiently high threshold $u>0$. In
this case, rearranging Eq. (9) provides the model

$\mathbf{Y}_{-i,t}=\bm{\alpha}_{-i,t}Y_{i,t}+Y_{i,t}^{\bm{\beta}_{-i,t}}\mathbf{Z}_{i,t}\;\big{|}\;(Y_{i,t}>u),$
(11)

where the residual random vector $\mathbf{Z}_{i,t}:=\\{Z_{j\mid
i,t}:j\in(1,2,3)\setminus i\\}\sim G_{-i,t}$ is independent of $Y_{i,t}$. For
inference, we make the working assumption that $G_{-i,t}(\cdot)$ does not
depend on time $t$ and follows a bivariate standard Gaussian copula with
correlation $\rho_{i}\in(-1,1)$ and delta-Laplace margins (Shooter et al.,,
2021). We hereafter drop the subscript $t$ from the notation. A random
variable that follows the delta-Laplace distribution with location, scale, and
shape parameters $\mu\in\mathbb{R},\;\sigma>0,$ and $\delta>0$, respectively,
has density function
$f(z)=\delta(2k\sigma\Gamma(\delta^{-1}))^{-1}\exp\\{-(|z-\mu|/(k\sigma))^{\delta}\\}$
for $k^{2}=\Gamma(\delta^{-1})/\Gamma(3\delta^{-1}),$ where $\Gamma(\cdot)$
denotes the standard gamma function. Note that when $\delta=1$ and $\delta=2$,
we have the Laplace and Gaussian densities, respectively.

Inference proceeds via maximum likelihood estimation. The model is fitted
separately for each conditioning variable $i=1,2,3$. For each $i$, we have
eight parameters in
$\bm{\alpha}^{(0)}_{-i}$,$\;\bm{\alpha}^{(1)}_{-i}$,$\;\bm{\beta}^{(0)}_{-i}$,$\;\bm{\beta}^{(1)}_{-i},$
as well as the seven parameters that characterise $G_{-i}$, that is, the
correlation $\rho_{i}$ and the three marginal parameters for each component of
$\mathbf{Z}_{i,t}$. After estimation of the parameters, we no longer require
the working Gaussian copula assumption for $\mathbf{Z}_{i,t}$. We instead use
observations $\tilde{\mathbf{y}}_{i,t}$ to derive the empirical residual
vector

$\tilde{\mathbf{z}}_{i,t}:=(\tilde{\mathbf{y}}_{i,t}-\hat{\bm{\alpha}}_{-i,t})/{\tilde{\mathbf{y}}_{i,t}}^{\hat{\bm{\beta}}_{-i,t}},$

where $\hat{\bm{\alpha}}_{-i,t}$ and $\hat{\bm{\beta}}_{-i,t}$ denote
estimates of ${\bm{\alpha}_{-i,t}}$ and ${\bm{\beta}_{-i,t}}$, respectively.
Assuming independence across time, we use the empirical residuals to provide
an empirical estimate $\tilde{G}_{-i}(\cdot)$ of $G_{-i}(\cdot)$.

We estimate the necessary exceedance probabilities, $p_{1}$ and $p_{2}$ in Eq.
(2), via the following Monte-Carlo procedure. We first note that, for fixed
$t\in\\{1,\dots,n\\}$, we require realisations of $\mathbf{Y}_{t}$, i.e.,
unconditional on an exceedance. We follow Richards et al., (2022); Richards et
al., 2023c  and obtain these realisations by drawing a realisation from

$\mathbf{Y}_{t}\;\bigg{|}\left(\max_{i=1,2,3}Y_{i,t}>u\right),$ (12)

with probability

$\Pr\left\\{\max_{i=1,2,3}Y_{i,t}>u\right\\},$ (13)

and, otherwise, drawing a realisation of
$\mathbf{Y}_{t}\mid(\max_{i=1,2,3}Y_{i,t}<u)$. As realisations of the latter
are unlikely to significantly impact estimates of $p_{1}$ and $p_{2}$, we draw
them empirically (Richards et al.,, 2022). In practice, we also replace
probability in Eq. (13) with an empirical estimate. In addition, as both
$\mathbf{Y}_{t}\mid(\max_{i=1,2,3}Y_{i,t}<u)$ and Eq. (13) depend on $t$, we
estimate them empirically by assuming stationarity over time.

To simulate from Eq. (12), we must first draw realisations of the conditional
exceedance model, $\mathbf{Y}_{-i,t}\mid(Y_{i,t}>u)$ in Eq. (11), using
Algorithm 1.

Algorithm 1 Simulating from Eq. (11).

For time $t$ and conditioning index $i\in\\{1,2,3\\}$:

  1. 1.

Simulate $E\sim\mbox{Exp}(1)$ and set $y_{i,t}=u+E$.

  2. 2.

Draw a residual vector $\mathbf{z}\sim\tilde{G}_{-i}(\cdot)$.

  3. 3.

Set
$\mathbf{y}_{-i,t}=\hat{\bm{\alpha}}_{-i,t}y_{i,t}+y_{i,t}^{\hat{\bm{\beta}}_{-i,t}}\mathbf{z}$.

We can then combine realisations from the three separate conditional
exceedance models, i.e., $\mathbf{Y}_{-i,t}\mid(Y_{i,t}>u)$ for each
$i=1,2,3$, into a single realisation of Eq. (12) using the importance sampling
(Wadsworth and Tawn,, 2022). Whilst this can be achieved for a single fixed
value of $t$, we note that $p_{1}$ and $p_{2}$ in Eq. (2) do not depend on the
time $t$. Hence, we treat $t$ as random during simulation and average over all
times $t\in\\{1,\dots,n\\}$ to produce approximate realisations of
$\mathbf{Y}$ with $t$ marginalised out. The full simulation algorithm is
detailed in Algorithm 2. Note that we back-transform from standard Laplace to
standard Gumbel margins after simulation, as the original data has the latter.

Algorithm 2 Simulating $\\{(Y_{1,t},Y_{2,t},Y_{3,t}):t\in\\{1,\dots,n\\}\\}$

  1. 1.

For $j=1,\dots,N^{\prime}$ with $N^{\prime}>N$:

    1. (a)

Draw a conditioning index $i_{j}\in\\{1,2,3\\}$ with equal probability.

    2. (b)

Draw a time $t_{j}\in\\{1,\dots,n\\}$ with equal probability.

    3. (c)

Simulate $\mathbf{y}^{(j)}$ using Algorithm 1 with time $t_{j}$ and
conditioning index $i_{j}$.

  2. 2.

Assign each simulated vector
$\mathbf{y}^{(j)}:=(y_{1}^{(j)},y_{2}^{(j)},y_{3}^{(j)})^{\prime}$ an
importance weight of

$\left\\{\sum_{i=1}^{3}\mathbbm{1}\\{{y}^{(j)}_{i}>u\\}\right\\}^{-1},$

for $j=1,\dots,N^{\prime}$, and sub-sample $N$ realisations from the
collection with probabilities proportional to these weights.

  3. 3.

With probability
$n^{-1}\sum^{n}_{t=1}\mathbbm{1}\\{\max_{i=1,2,3}\tilde{y}_{i,t}<u\\}$, re-
sample (with equal probability) from

$\left\\{\tilde{\mathbf{y}}_{t}:t\in\\{1,\dots,n\\},\max\limits_{i=1,2,3}\tilde{y}_{i,t}<u\right\\}.$

  4. 4.

Back-transform sample onto standard Gumbel margins.

###  3.4 Extremal dependence measure

A number of pairwise measures have been proposed to quantify the strength of
extremal dependence between random variables $(Y_{1},Y_{2})$, see, for
example, Heffernan, (2000). We adopt the extremal dependence measure (EDM)
proposed by Resnick, (2004) and Larsson and Resnick, (2012), which has been
extended to a multivariate setting by Cooley and Thibaud, (2019). Let
$(Y_{1},Y_{2})\in[0,\infty)^{2}$ be a regularly-varying random vector with
index $\alpha>0$; for full details on regular variation, see Resnick, (2007).
For some symmetric norm $\|\cdot\|$ on $[0,\infty)^{2}$, define the
transformation
$(R,\bm{\Omega}):=(\|Y_{1},Y_{2}\|,(Y_{1},Y_{2})/\|Y_{1},Y_{2}\|)$. Then,
there exists a sequence $b_{n}\rightarrow\infty$ and constant $c>0$ such that

$n\Pr\\{(b_{n}^{-1}R,\bm{\Omega})\in\cdot)\xrightarrow{v}c\nu_{\alpha}\times
H,$

where $\xrightarrow{v}$ denotes vague convergence, $\nu_{\alpha}$ is a measure
on $(0,\infty]$ such that $\nu_{\alpha}((y,\infty])=y^{-\alpha}$ and the
angular measure $H$ is defined on the unit sphere
$\mathcal{S}_{+}=\\{\mathbf{y}\in[0,\infty)^{2}\setminus\\{\mathbf{0}\\}:\|\mathbf{y}\|=1\\}$.
The extremal dependence measure (Larsson and Resnick,, 2012) between $Y_{1}$
and $Y_{2}$ is

$\mbox{EDM}(Y_{1},Y_{2})=\int_{\mathcal{S}_{+}}\omega_{1}\omega_{2}\mathop{}\\!\mathrm{d}H(\bm{\omega})=\lim_{y\rightarrow\infty}\mathbb{E}\left(\Omega_{1}\Omega_{2}\mid
R>y\right).$

Note that $\mbox{EDM}(Y_{1},Y_{2})\in[0,1]$ quantifies the strength of extreme
dependence between $Y_{1}$ and $Y_{2}$. If $\mbox{EDM}(Y_{1},Y_{2})=0$, then
$H$ concentrates all mass along the axes and hence $Y_{1}$ and $Y_{2}$ are
asymptotically independent (Resnick,, 2007). Conversely,
$\mbox{EDM}(Y_{1},Y_{2})$ is maximal if and only if $(Y_{1},Y_{2})$ has full
asymptotic dependence, or $H$ places all mass on the diagonal. Larsson and
Resnick, (2012) propose the empirical estimator

$\widehat{\operatorname{EDM}}\left(Y_{1},Y_{2}\right)=\frac{1}{N_{n}}\sum_{t=1}^{n}\frac{y_{1,t}}{r_{t}}\frac{y_{2,t}}{r_{t}}\mathbbm{1}{\left\\{r_{t}>u\right\\}},$
(14)

where $\bm{y}_{t}:=\left(y_{1,t},y_{2,t}\right)$, $t=1,\ldots,n$, is an i.i.d
sample from $\left(Y_{1},Y_{2}\right)$, $r_{t}=\left\|\bm{y}_{i}\right\|$, and
$N_{n}=\sum_{t=1}^{n}\mathbbm{1}{\left\\{r_{t}>u\right\\}}$ is the number of
exceedances of $\\{r_{t}:t=1,\dots,n\\}$ above some high threshold $u>0$. In
Section 4.2, we use pairwise EDM estimates to investigate the extremal
dependence structure of a high-dimensional random vector and then decompose it
into subvectors of strongly tail dependent variables.

###  3.5 Non-parametric tail probability estimation

We seek an estimator for exceedance probabilities in high-dimensional
multivariate random vectors $\mathbf{Y}\in\mathbb{R}^{d}$, as described in
Section 2.2. While scalable parametric models for high-dimensional
multivariate extremes do exist, as exemplified by Lederer and Oesting, (2023)
and the review by Engelke and Ivanovs, (2021), they often impose restrictive
assumptions about extremal dependence in $\mathbf{Y}$, particularly as the
dimension $d$ grows large. Examples include regular or hidden regular
variation (Resnick,, 2002). As our goal is point estimation of exceedance
probabilities, rather than a full characterisation of the joint upper tail of
$\mathbf{Y}$, we choose instead to use a non-parametric estimator for
exceedance probabilities that makes few assumptions on the joint upper tails.
In particular, we adopt the non-parametric multivariate tail probability
estimator proposed by Krupskii and Joe, (2019).

Denote by $F_{i}(\cdot)$ the distribution function111In our application,
$F_{i}(\cdot)$ is known to be the standard Gumbel distribution function for
all $i=1,\dots,d$; see Section 2. of $Y_{i}$ and let $U_{i}=F_{i}(Y_{i})$ for
$i=1,\dots,d$. Then, for $\phi\in(0,1)$, let

$p(\phi):=\Pr\\{U_{1}>1-\phi,\ldots,U_{d}>1-\phi\\}=\Pr\\{U_{\text{max}}<\phi\\},$
(15)

where $U_{\text{max}}=\max(U_{1},\dots,U_{d})$. Denote by
$C:[0,1]^{d}\mapsto\mathbb{R}$ the copula associated with $\mathbf{Y}$ and by
$\bar{C}$ the corresponding survival copula, such that
$p(\phi)=\bar{C}(1-\phi,\dots,1-\phi)$. Krupskii and Joe, (2019) construct
estimators of $p(\phi)$ for small $\phi>0$ by making assumptions about the
joint tail decay of $C$. In particular, they assume that $C$ (and $\bar{C}$)
has continuous partial derivatives and, as $\phi\downarrow 0$, that

$\displaystyle p(\phi)$
$\displaystyle=\Pr\\{U_{\text{max}}<\phi\\}=\lambda_{1}\phi+\lambda_{2}\phi^{1/\eta}\ell(\phi)+o\\{\phi^{1/\eta}\ell(\phi)\\},$
$\displaystyle p^{\prime}(\phi)$
$\displaystyle=\lambda_{1}+\lambda_{2}\phi^{1/\eta-1}\ell(\phi)/\eta+o\\{\phi^{1/\eta-1}\ell(\phi)\\},$
(16)

for $\lambda_{1},\lambda_{2}\geq 0,$ and $\eta<1$, and where $\ell(\cdot)$ is
a slowly-varying function at zero such that, for all $s>0$,
$\ell(s\phi)/\ell(\phi)\rightarrow 0$ as $\phi\downarrow 0$; these regularity
assumptions hold for a number of popular parametric copulas (see Krupskii and
Joe,, 2019). We note that, when $d=2$, $\eta$ corresponds to the coefficient
of tail dependence proposed by Ledford and Tawn, (1996); the extension to
$d$-dimensions is described by Eastoe and Tawn, (2012).

Under the model in Eq. (3.5), the parameters $\lambda_{1},\lambda_{2},$ and
$\eta$ can be estimated by defining a small threshold $\phi_{*}>\phi$, and
considering $U_{\text{max}}\mid U_{\text{max}}<\phi_{*}$. Then, as
$\phi\downarrow 0$,

$\Pr\\{U_{\text{max}}<\phi\mid
U_{\text{max}}<\phi_{*}\\}=p(\phi)/p(\phi^{*})\sim
k_{1}\phi+k_{2}\phi^{1/\eta},$ (17)

where $k_{1}=\lambda_{1}/(\lambda_{1}\phi_{*}+\lambda\phi_{*}^{1/\eta})$ and
$k_{2}=(1-k_{1}\phi_{*})/\phi_{*}^{1/\eta}$. By fixing $\phi_{*}$ and assuming
equality between the left and right hand sides of Eq. (17), parameter
estimates $\hat{k}_{1}$ and $\hat{\eta}_{1}$ can be computed using maximum
likelihood methods. Then, our estimate of $p(\phi)$ is

$\widehat{p(\phi)}:=\\{\hat{k}_{1}\phi+(1-\hat{k}_{1}\phi_{*})\phi^{1/\hat{\eta}}\\}\widehat{p(\phi_{*})},$
(18)

where $\widehat{p(\phi_{*})}$ is the empirical estimate of $p(\phi_{*})$.

The estimator in Eq. (18) can be directly used to estimate $p_{3}$ in Eq.
(2.2) by setting $\phi=\phi_{1}=1/300$. However, it cannot be used to estimate
$p_{4}$, as the formulation of the joint survival probability in Eq. (15) does
not account for different levels of marginal tail decay for each component of
$\mathbf{Y}$. In fact, we require an adaptation of $p(\phi)$ such that each
component $U_{i},i=1,\dots,d,$ exceeds a scaled value of the form
$1-c_{i}\phi$ for $0<c_{i}<1/\phi$. Setting $c_{i}=1$ for $i=1,\dots,25$ and
$c_{i}=12$, otherwise, yields $p_{3}$ in Eq. (2.2). We construct such an
estimator by considering, for $\mathbf{c}:=(c_{1},\dots,c_{d})^{\prime}$, the
probability

$p(\phi,\mathbf{c}):=\Pr\\{U_{1}>1-c_{1}\phi,\ldots,U_{d}>1-c_{d}\phi\\}=\Pr\\{U_{\text{max}}(\mathbf{c})<\phi\\},$
(19)

for weighted maxima
$U_{\text{max}}(\mathbf{c}):=\max\\{(1-U_{1})/c_{1},\dots,(1-U_{d})/c_{d}\\}$;
note that when $\mathbf{c}=(1,\dots,1)$, we have equivalence between Eq. (15)
and Eq. (19). An estimator of the form in Eq. (19) was alluded to by Krupskii
and Joe, (2019); however, they did not provide theoretical results for its
asymptotic behaviour. We choose to estimate $p(\phi,\mathbf{c})$ by assuming
that, as $\phi\downarrow 0$, that $\Pr\\{U_{\text{max}}(\mathbf{c})<\phi\\}$
has the same parametric form as $p(\phi)$ in Eq. (3.5). Inference for
$p(\phi,\mathbf{c})$ then follows in a similar manner as for $p(\phi)$, only
with samples of
$U_{\text{max}}(\mathbf{c})\mid(U_{\text{max}}(\mathbf{c})>\phi_{*})$ (rather
than $U_{\text{max}}\mid(U_{\text{max}}>\phi_{*})$) used for inference.

##  4 Results

###  4.1 Univariate quantile estimation

We now describe estimation of univariate conditional quantiles described in
Section 2.1. To estimate the conditional $q-$quantile of
$Y\mid(\mathbf{X}=\mathbf{x})$ for $q=0.9999$, we fit the GPD-GAM model
described in Section 3.1. The structure of $u(\mathbf{x})$,
$\sigma_{u}(\mathbf{x})$, and $\xi(\mathbf{x})$ are optimised, for a fixed
value of $\lambda$, by minimising the model’s BIC for different additive
combinations of linear and smooth functions of the covariates. For all
additive models that we fit, the linear and smooth terms are centered at zero,
and, during model fitting, their contribution to the model is set to zero when
the corresponding covariate value is missing.

We optimise $\lambda$ by using the threshold selection scheme proposed by
Varty et al., (2021) and extended by Murphy et al., (2023). We first follow
Heffernan and Tawn, (2001) and use the quantile and GPD GAM model estimates to
transform all data onto standard exponential margins. Then, we define a grid
of $n$ equally spaced probabilities, $\\{0<\lambda_{1}<...<\lambda_{n}<1$},
with $\lambda_{n}=1-(\lambda_{2}-\lambda_{1})$. Let $F_{E}(\cdot)$ be the
standard exponential distribution function and $\lambda_{*}\in(0,1)$ be a pre-
specified cut-off threshold. For all $i=1,\dots,n$, we define a vector of
weights $\mathbf{w}=(w_{1},\dots,w_{n})^{\prime}$ with components
$w_{i}=F_{E}^{-1}(\lambda_{i})/\sum^{n}_{j=1}F_{E}^{-1}(\lambda_{j})$ for
$\lambda_{i}>\lambda_{*}$ and zero, otherwise. We then define the tail-
weighted standardised mean absolute deviance (twsMAD) by
$(1/n)\sum^{n}_{i=1}w_{i}\lvert\tilde{\theta}(\lambda_{j})-F_{E}^{-1}(\lambda_{i})\rvert,$
where $\tilde{\theta}(\lambda_{j})$ denotes the empirical
$\lambda_{j}$-quantile of the standardised data. We set $\lambda_{*}=0.99$ and
evaluate the twsMAD for a range of candidate $\lambda$ values that subcede
$\lambda_{*}$. We then choose the optimal value which minimises this metric.

Using the BIC and twsMAD, our final model uses $\lambda=0.972$ and, for the
covariates, we include: season and wind speed as linear terms in
$u(\mathbf{x})$, $\sigma_{u}(\mathbf{x}),$ and $\xi(\mathbf{x})$; $V_{2}$,
$V_{3}$, $V_{4}$, and wind direction as smooth terms in $u(\mathbf{x})$ and
$\sigma_{u}(\mathbf{x})$. The median and the $50\%$ confidence intervals for
the test conditional $q$-quantiles (see Section 2.1) are estimated using the
empirical quantiles over $2500$ non-parametric bootstrap samples, and
presented in Figure 1. The true values of the quantiles are provided by
Rohrbeck et al., (2023). We observe good predictions of the conditional
quantiles, with our estimates close to the true values. Our framework attains
a $38\%$ coverage rate, which is a slight underestimation of the normal
$50\%$.

Figure 1: Predicted conditional quantiles against true values for the 100 test
covariate sets for sub-challenge C1. Vertical error bars denote the $50\%$
confidence interval, with the black points denoting the median over bootstrap
samples.

To estimate the required extreme quantile of $Y$, we use the NBE framework
detailed in Section 3.2 with $m=n=21000$ independent replicates. Neural
networks are trained using the Adaptive Moment Estimation (Adam) algorithm
(Kingma and Ba,, 2014) with a mini-batch size of 32 and, in order to improve
numerical stability, input data are standardised prior to training by
subtracting and dividing by the mean and standard deviation (evaluated over
the entire training data set), respectively. We use
$\lvert\vartheta_{\mbox{train}}\rvert:=K=350,000$ and
$\lvert\vartheta_{\mbox{val}}\rvert:=K/5$ parameter values for training and
validation, with $M=4\times 10^{7}$ replicates used to estimate the
theoretical quantiles (see Section 3.2). When producing training data, we use
an empirical prior distribution for
$\\{(u_{t},\sigma_{t},\xi_{t}):t=1,\dots,m\\}$. To construct this prior, we
estimate the GPD-GAM model, as described above but with $\lambda=0.6$, across
750 bootstrap samples222A smaller value of $\lambda$ is used here to increase
the variety of the input data used to train the NBE. A larger value of
$\lambda$ would lead to more observations having to be drawn from the
empirical distribution to construct the training data.. Draws from the
empirical prior on $(u_{t},\sigma_{t},\xi_{t})$ then correspond to a sample
from all bootstrap parameter estimates. To increase the variety of models (and
hence, quantiles) used to train the estimator, we permute the values of
$u_{t}$, $\sigma_{t},$ and $\xi_{t}$ across
$\\{(u_{t},\sigma_{t},\xi_{t}):t=1,\dots,m\\}$ when constructing the training
and validation parameter sets.

The estimator is trained using early stopping (see, e.g., Prechelt,, 2002);
the validation loss is recorded at each epoch during training, which halts if
the validation loss has not decreased for 10 epochs. We choose the optimal
architecture for the neural networks, $\phi$ and $\bm{\psi}$ in Eq. (7), by
minimising the risk for a test set of $1000$ parameter values (not used in
training); this is provided in Table 1.

Table 1: Optimal DNN architecture used in Section 4.1. All layers used rectified linear unit (ReLU) activation functions, except the final layer, which used the identity function.  Function | input dimension | output dimension | parameters
---|---|---|---
$\bm{\psi}$ | [1] | [48] | $49$
[48] | [48] | 2352
$\mathbf{a}$ | [48] | [48] | 0
$\phi$ | [48] | [48] | 2352
[48] | [1] | 49
total trainable parameters: | 4802

NBEs are trained on GPUs randomly selected from within KAUST’s Ibex cluster,
see https://www.hpc.kaust.edu.sa/ibex/gpu_nodes for details (last accessed
13/07/2023).

To illustrate the efficacy of our estimator, Figure 2 presents extreme
quantile estimates for $1000$ test data sets. We observe that estimates are
generally aligned to the left of the diagonal, i.e., the majority of estimates
in Figure 2 are overestimates. This is a consequence of the asymmetric loss
function (also presented in Figure 2), which favours conservative estimates of
the quantile. We estimate the $q$-quantile for $q=1-(6\times 10^{4})^{-1}$ to
be 201.25. An estimated 95% confidence interval of (200.79, 201.73) was
derived using a non-parametric bootstrap; due to the amortised nature of our
neural estimator, non-parametric bootstrap-based uncertainty estimation is
extremely fast, as the estimator does not need to be retrained for every new
sample.

Figure 2: Left: Extreme quantile estimates for 1000 test data sets. Dark blue
points denote those for which $L(\theta,\hat{\theta})=0$. Right: Loss
function, $L(\theta_{0},\hat{\theta})$ in Eq. (1), for $\theta_{0}:=200$.

###  4.2 Extremal dependence estimation

For sub-challenge C3, we estimate $p_{1}$ and $p_{2}$ in Eq. (2) using the
conditional extremes model described in Section 3.3. We perform model
selection to identify the best covariates to include in the linear model for
the dependence parameters, $\bm{\alpha}_{-i,t}$ and $\bm{\beta}_{-i,t}$ in Eq.
(10), by minimising the AIC over the three separate model fits, i.e.,
$\mathbf{Y}_{-i,t}\mid(Y_{i,t}>u),i=1,2,3$, with $u$ set to the $90\%$
standard Laplace quantile. The best fitting model did not include any
transformation (e.g, $\log$ or square) of the covariates, and uses a
homogeneous $\bm{\beta}_{-i,t}$ function, i.e., with
$\bm{\beta}_{-i,t}=\bm{\beta}_{-i}$ for all $t=1,\dots,n,$ and for $i=1,2,3$.

We use a visual diagnostic to optimise the exceedance threshold $u$. We fit
the three conditional models for a grid of $u$ values, use Algorithm 2 to
simulate $N:=10^{7}$ realisations from the fitted model for
$\mathbf{Y}=(Y_{1},Y_{2},Y_{3}),$ i.e., with time $t$ marginalised out, and
then compute realisations of $R:=Y_{1}+Y_{2}+Y_{3}$. The upper tails of the
aggregate variable $R$ are heavily influenced by the extremal dependence
structure in the underlying $(Y_{1},Y_{2},Y_{3})$ (Richards and Tawn,, 2022).
Hence, we can exploit the empirical distribution of the aggregate as a
diagnostic measure for the quality of the dependence model fit. We produce a
Q-Q plot of $R$ and choose the optimal $u$ to be the value which illustrates
the best alignment for extreme quantiles (see, also, Richards et al.,, 2022;
Richards et al., 2023c, ). This is the $95\%$ standard Laplace quantile, and
we illustrate the corresponding Q-Q plot in Figure 3. Table 2 provides the
corresponding parameters estimates
${\hat{\bm{\alpha}}^{(0)}_{-i}},{\hat{\bm{\alpha}}^{(1)}_{-i}},$ and
$\hat{\bm{\beta}}_{-i}$, for $i=1,2,3$.

Table 2: Sub-challenge C3: Parameter estimates ${\hat{\bm{\alpha}}^{(0)}_{-i}},{\hat{\bm{\alpha}}^{(1)}_{-i}},$ and $\hat{\bm{\beta}}_{-i}$, for the conditional extremal dependence model, $\mathbf{Y}_{-i,t}\mid Y_{i,t}>u$, $i=1,2,3$. Here, $\hat{\bm{\alpha}}^{(1,j)}_{-i}$ denotes the $j^{\text{th}}$ column of $\hat{\bm{\alpha}}^{(1)}_{-i}$ with $j=1$ and $j=2$ corresponding to the linear coefficient vector for Season and Atmosphere, respectively. Values in the $j^{\text{th}}$ row give the component of the parameter vector for $\mathbf{Y}_{-i,t}$ that corresponds to $Y_{j,t}$, $j=1,2,3$ (and so are missing for $j=i$).  | $i=1$ | $i=2$ | $i=3$
---|---|---|---
| $\hat{\bm{\alpha}}^{(0)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,1)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,2)}_{-i}$ | $\hat{\bm{\beta}}_{-i}$ | $\hat{\bm{\alpha}}^{(0)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,1)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,2)}_{-i}$ | $\hat{\bm{\beta}}_{-i}$ | $\hat{\bm{\alpha}}^{(0)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,1)}_{-i}$ | $\hat{\bm{\alpha}}^{(1,2)}_{-i}$ | $\hat{\bm{\beta}}_{-i}$
$Y_{1,t}$ | - | - | - | - | $-0.08$ | $0.22$ | $0.08$ | $0.81$ | $0.15$ | $-0.33$ | $0.01$ | $0.97$
$Y_{2,t}$ | $0.56$ | $-0.06$ | $0.07$ | $0.70$ | - | - | - | - | $0.01$ | $-0.07$ | $0.01$ | $0.45$
$Y_{3,t}$ | $0.28$ | $-0.02$ | $0.07$ | $0.82$ | $0.24$ | $-0.14$ | $0.14$ | $0.09$ | - | - | - | -

Figure 3: Top-left panel: Q-Q plot for the simulated and empirical sum
$R=Y_{1}+Y_{2}+Y_{3}$. Other panels provide scatter plots of (black)
observations and (blue) realisations from the fitted conditional extremal
dependence models. Red vertical and horizontal lines correspond to the
exceedance threshold $u$, i.e., the $95\%$ marginal quantile.

Example realisations from the fitted model are illustrated in Figure 3. Good
model fit is illustrated, as the realisations mimic the extremal
characteristics of the observed data. Probabilities $p_{1}$ and $p_{2}$ are
estimated empirically from the aforementioned $N=10^{7}$ realisations of
$(Y_{1},Y_{2},Y_{3})$, and found to be $\hat{p}_{1}=3.13\times 10^{-5}$ and
$\hat{p}_{2}=2.29\times 10^{-5}$, respectively.

For sub-challenge C4, we estimate $p_{3}$ and $p_{4}$ in Eq. (2.2) using the
non-parametric tail probability estimators defined in Eq. (18) and Eq. (19),
respectively. However, we first seek to decompose the high-dimensional random
vector $\mathbf{Y}$ into strongly tail-dependent sub-vectors. Our reasoning is
threefold: i) the estimator in Eq. (19) requires strong assumptions about the
joint upper tail decay of the random vector $\mathbf{Y}$ which are unlikely to
hold for high dimension $d$; ii) reliable estimation of the intermediate
exceedance probability, $\widehat{p(\phi_{*})}$ in Eq. (18), becomes difficult
as $d$ grows large; iii) Krupskii and Joe, (2019) observed through simulations
that the estimator in Eq. (18) provides more accurate estimates when the true
value of $\eta$ is close to one, i.e., when $\mathbf{Y}$ exhibits strong upper
tail dependence.

We identify sub-vectors by estimating the EDM between all pairs of components
of $\mathbf{Y}$, with $\alpha=2$, $\|\cdot\|$ set to the $L_{2}$ norm, and
with $u$ in Eq. (14) taken to be the empirical 0.99-quantile. Estimates are
presented in a heatmap in Figure 4, with values grouped using hierarchical
clustering. An appropriate grouping can be found by visual inspection of
Figure 4; the clustered heatmap produces a block-diagonal matrix with five
blocks (with size ranging from 8–13 components), and with elements of the off-
diagonal blocks close to zero. For less well-behaved data, the number of
clusters (and sub-vectors of $\mathbf{Y}$) can instead be determined by
placing a cut-off at an appropriate point on the clustering dendrogram.

Figure 4: Left: Heatmap of pairwise EDM estimates grouped using hierarchical
clustering, alongside the corresponding dendrogram. Right: Threshold stability
plots of estimated log-exceedance probabilities, $\log(\hat{p_{3}})$ and
$\log(\hat{p_{4}}),$ against the intermediate quantile $\phi_{*}$ using in
their estimation. Red horizontal lines denote the optimal choice of
$\phi_{*}$.

We decompose $\mathbf{Y}$ into the five sub-vectors suggested by Figure 4,
which we denote by $\mathcal{Y}_{i}$ for $i=1,\dots,5$. Estimation of $p_{3}$
then follows under the approximation

$p_{3}=\Pr(Y_{1}>s_{1},\dots,Y_{50}>s_{1})\approx\prod_{i=1}^{5}\Pr\left(\bigcap\limits_{Y\in\mathcal{Y}_{i}}Y>s_{1}\right),$
(20)

and similarly for $p_{4}$ (with the exceedance value $s_{1}$ or $s_{2}$ chosen
appropriately). That is, we assume complete independence between the sub-
vectors of $\mathbf{Y}$ to approximate $p_{3}$ and $p_{4}$. Exceedance
probabilities for each sub-vector are estimated using the methodology
described in Section 3.5. We fix the intermediate quantile level $\phi_{*}$ in
Eq. (18) across all sub-vectors and for estimation of both $p_{3}$ and
$p_{4}$. An optimal value $\hat{\phi}_{*}$ is chosen via a threshold stability
plot (see, e.g., Coles,, 2001). In Figure 4, we present estimates of
$\log(p_{3})$ and $\log(p_{4})$ for a sequence of $\phi_{*}$ values in
$[0.1,0.4]$. We choose $\hat{\phi}_{*}$ as the smallest $\phi_{*}$ such that
there is visual stability in Figure 4 for estimates of both $\log(p_{3})$ and
$\log(p_{4})$ when $\phi_{*}>\hat{\phi}_{*}$. In this case, we adopt
$\hat{\phi}_{*}=0.25$ as the optimal value. We use a non-parametric bootstrap
with 200 samples to assess uncertainty in estimates of $p_{3}$ and $p_{4}$.
The resulting bootstrap median333We submitted the bootstrap median estimates
for the data competition. estimates of $\log(p_{3})$ and $\log(p_{4}$) (and
their estimated 95% confidence intervals) are $-59.38$ $(-68.46,-49.82)$ and
$-55.24$ $(-60.12,-49.79)$, respectively.

##  5 Conclusion

We propose four frameworks for the estimation of exceedance probabilities and
levels associated with extreme events. To estimate extreme conditional
quantiles, we adopt an additive model that represents parameters in a peaks-
over-threshold model as splines. To estimate univariate quantiles in an
unconditional setting, we construct a neural Bayes estimator that estimates a
quantile subject to a conservative loss function. Our new approach to quantile
estimation is amortised, likelihood-free, and requires few parametric
assumptions about the underlying distribution. For multivariate data, we
consider two frameworks: i) in the presence of covariates, we adopt a non-
stationary conditional extremal dependence model to capture linear trends in
extremal dependence parameters, and ii) in the absence of covariates, we
propose a non-parametric estimator for multivariate tail probabilities that
can be applied to high-dimensional ($d=50)$ data via a sparse decomposition
using pairwise extremal dependence measures. We validate these modelling
approaches using simulated data from the EVA2023 data prediction challenge
(Rohrbeck et al.,, 2023).

We illustrated the efficacy of additive GPD regression models for estimating
extreme conditional quantiles. Machine learning-based models offer a
potentially more flexible alternative to additive models. These include
approaches using decision trees (Farkas et al., 2021b, ; Farkas et al., 2021a,
), random forests (Gnecco et al.,, 2022), gradient boosting (Velthoen et al.,,
2023), and neural networks (Pasche and Engelke,, 2022; Richards et al., 2023a,
), all of which have been previously utilised to fit GPD regression models. We
could also consider regression models with alternative asymptotically-
justified extreme value response distributions, e.g., mixture models (Carreau
and Bengio,, 2007), point process models (Richards and Huser,, 2022), or the
extended generalised Pareto distribution (Papastathopoulos and Tawn,, 2013;
Naveau et al.,, 2016), which has been applied in a regression setting by
Carrer and Gaetan, (2022), de Carvalho et al., 2022b , and Cisneros et al.,
2023b .

We designed a neural Bayes estimator (NBE) to perform simulation-based
inference for extreme quantiles. To train this estimator, we constructed an
empirical prior measure $\pi(\cdot)$ for the extreme quantile, by simulating
from a family of pre-defined peaks-over-threshold models. The resulting NBE
was optimal with respect to the prior measure $\pi(\cdot)$. Whilst our
estimator performed well in the application, the underlying truth was a peaks-
over-threshold model (Rohrbeck et al.,, 2023), and, hence, the truth
distribution was contained within the class of distributions covered by
$\pi(\cdot)$. We note that our estimator may not perform as well when the
converse holds, and the true distribution of $Y$ is not approximated well by
any family in $\pi(\cdot)$. Further study of the optimality properties of our
proposed NBE is required and, in particular, its efficacy under model
misspecification.

We proposed an extension of the conditional extremal dependence regression
model (Winter et al.,, 2016) and highlighted its efficacy for estimation of
extreme exceedance probabilities. Our conditional model performed well in
practice, but alternative multivariate extremal dependence regression models,
which rely on more classical assumptions about the extremal behaviour of
$\mathbf{Y}$, e.g., regular or hidden regular variation, could have been
tested (see, e.g., Cooley et al.,, 2012; de Carvalho et al., 2022a, ; Murphy-
Barltrop and Wadsworth,, 2022).

## Acknowledgments

All authors were supported by funding from the King Abdullah University of
Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award
No. OSR-CRG2020-4394. Support from the KAUST Supercomputing Laboratory is
gratefully acknowledged. The authors are grateful to Raphaël Huser and Matthew
Sainsbury-Dale for helpful discussions and to Jonathan Tawn, Christian
Rohrbeck, and Emma Simpson of Lancaster University, University of Bath, and
University College London, respectively, for organisation of the data
challenge that motivated this work.

## References

  * Beirlant et al., (2006) Beirlant, J., Goegebeur, Y., Segers, J., and Teugels, J. L. (2006).  Statistics of Extremes: Theory and Applications.  John Wiley & Sons.
  * Bezanson et al., (2017) Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017).  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98.
  * Carreau and Bengio, (2007) Carreau, J. and Bengio, Y. (2007).  A hybrid Pareto model for conditional density estimation of asymmetric fat-tail data.  In Artificial Intelligence and Statistics, pages 51–58. PMLR.
  * Carrer and Gaetan, (2022) Carrer, N. L. and Gaetan, C. (2022).  Distributional regression models for Extended Generalized Pareto distributions.  arXiv preprint arXiv:2209.04660.
  * Chavez-Demoulin and Davison, (2005) Chavez-Demoulin, V. and Davison, A. C. (2005).  Generalized additive modelling of sample extremes.  Journal of the Royal Statistical Society Series C: Applied Statistics, 54(1):207–222.
  * (6) Cisneros, D., Gong, Y., Yadav, R., Hazra, A., and Huser, R. (2023a).  A combined statistical and machine learning approach for spatial prediction of extreme wildfire frequencies and sizes.  Extremes, 26(2):301–330.
  * (7) Cisneros, D., Richards, J., Dahal, A., Lombardo, L., and Huser, R. (2023b).  Deep graphical regression for jointly moderate and extreme Australian wildfires.  arXiv preprint arXiv:2308.14547.
  * Coles, (2001) Coles, S. (2001).  An Introduction to Statistical Modeling of Extreme Values, volume 208.  Springer.
  * Cooley et al., (2012) Cooley, D., Davis, R. A., and Naveau, P. (2012).  Approximating the conditional density given large observed values via a multivariate extremes framework, with application to environmental data.  The Annals of Applied Statistics, 6(4):1406–1429.
  * Cooley et al., (2019) Cooley, D., Hunter, B. D., and Smith, R. L. (2019).  Univariate and multivariate extremes for the environmental sciences.  In Gelfand, A. E., Fuentes, M., Hoeting, J. A., and Smith, R. L., editors, Handbook of Environmental and Ecological Statistics, chapter 8, pages 153–180. Chapman and Hall/CRC.
  * Cooley and Thibaud, (2019) Cooley, D. and Thibaud, E. (2019).  Decompositions of dependence for high-dimensional extremes.  Biometrika, 106(3):587–604.
  * Cremanns and Roos, (2017) Cremanns, K. and Roos, D. (2017).  Deep Gaussian covariance network.  arXiv preprint arXiv:1710.06202.
  * Davison and Smith, (1990) Davison, A. C. and Smith, R. L. (1990).  Models for exceedances over high thresholds.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 52(3):393–425.
  * (14) de Carvalho, M., Kumukova, A., and Dos Reis, G. (2022a).  Regression-type analysis for multivariate extreme values.  Extremes, 25(4):595–622.
  * (15) de Carvalho, M., Pereira, S., Pereira, P., and de Zea Bermudez, P. (2022b).  An extreme value Bayesian lasso for the conditional left and right tails.  Journal of Agricultural, Biological and Environmental Statistics, 27(2):1–18.
  * Eastoe and Tawn, (2012) Eastoe, E. F. and Tawn, J. A. (2012).  Modelling the distribution of the cluster maxima of exceedances of subasymptotic thresholds.  Biometrika, 99(1):43–55.
  * Engelke and Ivanovs, (2021) Engelke, S. and Ivanovs, J. (2021).  Sparse structures for multivariate extremes.  Annual Review of Statistics and Its Application, 8:241–270.
  * (18) Farkas, S., Heranval, A., Lopez, O., and Thomas, M. (2021a).  Generalized Pareto regression trees for extreme events analysis.  arXiv preprint arXiv:2112.10409.
  * (19) Farkas, S., Lopez, O., and Thomas, M. (2021b).  Cyber claim analysis using generalized Pareto regression trees with applications to insurance.  Insurance: Mathematics and Economics, 98:92–105.
  * Fasiolo et al., (2021) Fasiolo, M., Wood, S. N., Zaffran, M., Nedellec, R., and Goude, Y. (2021).  Fast calibrated additive quantile regression.  Journal of the American Statistical Association, 116(535):1402–1412.
  * Gerber and Nychka, (2021) Gerber, F. and Nychka, D. (2021).  Fast covariance parameter estimation of spatial Gaussian process models using neural networks.  Stat, 10(1):e382.
  * Gnecco et al., (2022) Gnecco, N., Terefe, E. M., and Engelke, S. (2022).  Extremal random forests.  arXiv preprint arXiv:2201.12865.
  * Gong and Huser, (2022) Gong, Y. and Huser, R. (2022).  Asymmetric tail dependence modeling, with application to cryptocurrency market data.  The Annals of Applied Statistics, 16(3):1822–1847.
  * Guerrero et al., (2023) Guerrero, M. B., Huser, R., and Ombao, H. (2023).  Conex–connect: Learning patterns in extremal brain connectivity from multichannel eeg data.  The Annals of Applied Statistics, 17(1):178–198.
  * Haan and Ferreira, (2006) Haan, L. and Ferreira, A. (2006).  Extreme Value Theory: An Introduction, volume 3.  Springer.
  * Heffernan, (2000) Heffernan, J. E. (2000).  A directory of coefficients of tail dependence.  Extremes, 3:279–290.
  * Heffernan and Resnick, (2007) Heffernan, J. E. and Resnick, S. I. (2007).  Limit laws for random vectors with an extreme component.  The Annals of Applied Probability, 17(2):537 – 571.
  * Heffernan and Tawn, (2001) Heffernan, J. E. and Tawn, J. A. (2001).  Extreme value analysis of a large designed experiment: a case study in bulk carrier safety.  Extremes, 4:359–378.
  * Heffernan and Tawn, (2004) Heffernan, J. E. and Tawn, J. A. (2004).  A conditional approach for multivariate extreme values (with discussion).  Journal of the Royal Statistical Society Series B: Methodology, 66(3):497–546.
  * Huser and Wadsworth, (2022) Huser, R. and Wadsworth, J. L. (2022).  Advances in statistical modeling of spatial extremes.  Wiley Interdisciplinary Reviews: Computational Statistics, 14(1):e1537.
  * Keef et al., (2013) Keef, C., Papastathopoulos, I., and Tawn, J. A. (2013).  Estimation of the conditional distribution of a multivariate variable given that one of its components is large: Additional constraints for the Heffernan and Tawn model.  Journal of Multivariate Analysis, 115:396–404.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Krupskii and Joe, (2019) Krupskii, P. and Joe, H. (2019).  Nonparametric estimation of multivariate tail probabilities and tail dependence coefficients.  Journal of Multivariate Analysis, 172:147–161.
  * Larsson and Resnick, (2012) Larsson, M. and Resnick, S. I. (2012).  Extremal dependence measure and extremogram: the regularly varying case.  Extremes, 15(2):231–256.
  * Lederer and Oesting, (2023) Lederer, J. and Oesting, M. (2023).  Extremes in high dimensions: methods and scalable algorithms.  arXiv preprint arXiv:2303.04258.
  * Ledford and Tawn, (1996) Ledford, A. W. and Tawn, J. A. (1996).  Statistics for near independence in multivariate extreme values.  Biometrika, 83(1):169–187.
  * Lehmann and Casella, (1998) Lehmann, E. L. and Casella, G. (1998).  Theory of Point Estimation.  Springer, New York, NY, 2nd edition.
  * Lenzi et al., (2023) Lenzi, A., Bessac, J., Rudi, J., and Stein, M. L. (2023).  Neural networks for parameter estimation in intractable models.  Computational Statistics & Data Analysis, 185:107762.
  * Lenzi and Rue, (2023) Lenzi, A. and Rue, H. (2023).  Towards black-box parameter estimation.  arXiv preprint arXiv:2303.15041.
  * McCullagh, (2002) McCullagh, P. (2002).  What is a statistical model?  The Annals of Statistics, 30(5):1225–1310.
  * Murphy et al., (2023) Murphy, C., Tawn, J. A., and Varty, Z. (2023).  Automated threshold selection and associated inference uncertainty for univariate extremes.  arXiv preprint arXiv:2310.17999.
  * Murphy-Barltrop and Wadsworth, (2022) Murphy-Barltrop, C. and Wadsworth, J. (2022).  Modelling non-stationarity in asymptotically independent extremes.  arXiv preprint arXiv:2203.05860.
  * Naveau et al., (2016) Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016).  Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection.  Water Resources Research, 52(4):2753–2769.
  * Papastathopoulos and Tawn, (2013) Papastathopoulos, I. and Tawn, J. A. (2013).  Extended generalised Pareto models for tail estimation.  Journal of Statistical Planning and Inference, 143(1):131–143.
  * Pasche and Engelke, (2022) Pasche, O. C. and Engelke, S. (2022).  Neural networks for extreme quantile regression with an application to forecasting of flood risk.  arXiv preprint arXiv:2208.07590.
  * Pickands, (1975) Pickands, J. (1975).  Statistical inference using extreme order statistics.  The Annals of Statistics, 3(1):119–131.
  * Prechelt, (2002) Prechelt, L. (2002).  Early stopping-but when?  In Neural Networks: Tricks of the trade, pages 55–69. Springer.
  * Rai et al., (2023) Rai, S., Hoffman, A., Lahiri, S., Nychka, D. W., Sain, S. R., and Bandyopadhyay, S. (2023).  Fast parameter estimation of generalized extreme value distribution using neural networks.  arXiv preprint arXiv:2305.04341.
  * Redondo et al., (2023) Redondo, P. V., Huser, R., and Ombao, H. (2023).  Measuring information transfer between nodes in a brain network through spectral transfer entropy.  arXiv preprint arXiv:2303.06384.
  * Reiss et al., (2007) Reiss, R.-D., Thomas, M., and Reiss, R. (2007).  Statistical Analysis of Extreme Values, volume 3.  Springer.
  * Resnick, (2002) Resnick, S. (2002).  Hidden regular variation, second order regular variation and asymptotic independence.  Extremes, 5:303–336.
  * Resnick, (2004) Resnick, S. (2004).  The extremal dependence measure and asymptotic independence.  Stochastic Models, 20(2):205–227.
  * Resnick, (2007) Resnick, S. I. (2007).  Heavy-Tail Phenomena: Probabilistic and Statistical Modeling.  Springer Science & Business Media.
  * Richards and Huser, (2022) Richards, J. and Huser, R. (2022).  Regression modelling of spatiotemporal extreme U.S. wildfires via partially-interpretable neural networks.  arXiv preprint arXiv:2208.07581.
  * (55) Richards, J., Huser, R., Bevacqua, E., and Zscheischler, J. (2023a).  Insights into the drivers and spatio-temporal trends of extreme Mediterranean wildfires with statistical deep-learning.  Artificial Intelligence for the Earth Systems, 2(4):e220095.
  * (56) Richards, J., Sainsbury-Dale, M., Zammit-Mangion, A., and Huser, R. (2023b).  Neural Bayes estimators for censored peaks-over-threshold models.  arXiv preprint arXiv:2306.15642.
  * Richards and Tawn, (2022) Richards, J. and Tawn, J. A. (2022).  On the tail behaviour of aggregated random variables.  Journal of Multivariate Analysis, 192:105065.
  * Richards et al., (2022) Richards, J., Tawn, J. A., and Brown, S. (2022).  Modelling extremes of spatial aggregates of precipitation using conditional methods.  The Annals of Applied Statistics, 16(4):2693–2713.
  * (59) Richards, J., Tawn, J. A., and Brown, S. (2023c).  Joint estimation of extreme spatially aggregated precipitation at different scales through mixture modelling.  Spatial Statistics, 53:100725.
  * Rohrbeck et al., (2018) Rohrbeck, C., Eastoe, E. F., Frigessi, A., and Tawn, J. A. (2018).  Extreme value modelling of water-related insurance claims.  The Annals of Applied Statistics, 12(1):246–282.
  * Rohrbeck et al., (2023) Rohrbeck, C., Simpson, E. S., and Tawn, J. A. (2023).  Editorial: EVA 2023 data challenge.
  * (62) Sainsbury-Dale, M., Richards, J., Zammit-Mangion, A., and Huser, R. (2023a).  Neural bayes estimators for irregular spatial data using graph neural networks.  arXiv preprint arXiv:2310.02600.
  * (63) Sainsbury-Dale, M., Zammit-Mangion, A., and Huser, R. (2023b).  Likelihood-free parameter estimation with neural Bayes estimators.  The American Statistician (in press).
  * Segers, (2012) Segers, J. (2012).  Max-stable models for multivariate extremes.  REVSTAT–Statistical Journal, 10(1):61–82.
  * Shao et al., (2022) Shao, X., Hazra, A., Richards, J., and Huser, R. (2022).  Flexible modeling of nonstationary extremal dependence using spatially-fused LASSO and ridge penalties.  arXiv preprint arXiv:2210.05792.
  * Shooter et al., (2021) Shooter, R., Tawn, J., Ross, E., and Jonathan, P. (2021).  Basin-wide spatial conditional extremes for severe ocean storms.  Extremes, 24:241–265.
  * Smith, (1999) Smith, R. L. (1999).  Bayesian and frequentist approaches to parametric predictive inference (with discussion).  In Bernardo, J., Berger, J., Dawid, A., and Smith, A., editors, Bayesian Statistics, volume 6, pages 589–612. Oxford University Press, London.
  * Varty et al., (2021) Varty, Z., Tawn, J. A., Atkinson, P. M., and Bierman, S. (2021).  Inference for extreme earthquake magnitudes accounting for a time-varying measurement process.  arXiv preprint arXiv:2102.00884.
  * Velthoen et al., (2023) Velthoen, J., Dombry, C., Cai, J.-J., and Engelke, S. (2023).  Gradient boosting for extreme quantile regression.  Extremes, pages 1–29.
  * Wadsworth and Tawn, (2022) Wadsworth, J. L. and Tawn, J. (2022).  Higher-dimensional spatial extremes via single-site conditioning.  Spatial Statistics, 51:100677.
  * Winter et al., (2016) Winter, H. C., Tawn, J. A., and Brown, S. J. (2016).  Modelling the effect of the El Niño-Southern Oscillation on extreme spatial temperature events over Australia.  The Annals of Applied Statistics, 10(4):2075 – 2101.
  * Youngman, (2019) Youngman, B. D. (2019).  Generalized additive models for exceedances of high thresholds with an application to return level estimation for US wind gusts.  Journal of the American Statistical Association, 114(528):1865–1879.
  * Zaheer et al., (2017) Zaheer, M., Kottur, S., Ravanbakhsh, S., Poczos, B., Salakhutdinov, R. R., and Smola, A. J. (2017).  Deep sets.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Zammit-Mangion and Wikle, (2020) Zammit-Mangion, A. and Wikle, C. K. (2020).  Deep integro-difference equation models for spatio-temporal forecasting.  Spatial Statistics, 37:100408.
  * Zhang et al., (2023) Zhang, Z., Arellano‐Valle, R. B., Genton, M. G., and Huser, R. (2023).  Tractable bayes of skew‐elliptical link models for correlated binary data.  Biometrics, 79(3):1788–1800.