# Neural Networks for Extreme Quantile Regression
with an Application to Forecasting of Flood Risk

Olivier C. Pasche
Research Center for Statistics,
University of Geneva, Switzerland,
<EMAIL_ADDRESS>
&Sebastian Engelke
Research Center for Statistics,
University of Geneva, Switzerland,
<EMAIL_ADDRESS>

###### Abstract

Risk assessment for extreme events requires accurate estimation of high
quantiles that go beyond the range of historical observations. When the risk
depends on the values of observed predictors, regression techniques are used
to interpolate in the predictor space. We propose the EQRN model that combines
tools from neural networks and extreme value theory into a method capable of
extrapolation in the presence of complex predictor dependence. Neural networks
can naturally incorporate additional structure in the data. We develop a
recurrent version of EQRN that is able to capture complex sequential
dependence in time series. We apply this method to forecasting of flood risk
in the Swiss Aare catchment. It exploits information from multiple covariates
in space and time to provide one-day-ahead predictions of return levels and
exceedances probabilities. This output complements the static return level
from a traditional extreme value analysis and the predictions are able to
adapt to distributional shifts as experienced in a changing climate. Our model
can help authorities to manage flooding more effectively and to minimize their
disastrous impacts through early warning systems.

_Keywords_ Extreme value theory $\cdot$ Generalized Pareto distribution
$\cdot$ Machine learning $\cdot$ Prediction $\cdot$ Recurrent neural network

##  1 Introduction

Risk assessment is concerned with the analysis of rare events, which have
small occurrence probabilities but carry the potential of serious impacts on
our health, the environment or the economy. Examples for such extreme events
are floods in hydrology, crises in the financial system or heatwaves in a
changing climate. In theses applications, the quantity of interest is
typically a univariate response variable $Y$ representing the random risk
factor. The goal is to estimate a quantile $Q(\tau)=F^{-1}_{Y}(\tau)$ at level
$\tau\in[0,1]$, where we denote by $F^{-1}_{Y}$ the generalized inverse of the
distribution function of $Y$.

Since for risk quantification the level $\tau$ is usually very close to 1 so
that the quantile $Q(\tau)$ goes beyond the range of the data, the classical
approach is to model the tail of the distribution of $Y$ using extrapolation
results from extreme value theory. Two main approaches exist. When $Y$
represents, say, a daily quantity, then the generalized Pareto distribution
(GPD) can be used to approximate the tail above a high threshold $u$ by
(Balkema and de Haan, 1974; Pickands III, 1975)

$\displaystyle\mathbb{P}(Y>y)\approx\mathbb{P}(Y>u)\left(1+\xi\frac{y-u}{\sigma(u)}\right)_{+}^{-1/\xi},\quad
y\geq u,$ (1)

where $\xi\in\mathbb{R}$ and $\sigma(u)>0$ are the shape and scale parameters,
and the second factor on the right-hand side is the GPD approximation to
$\mathbb{P}(Y>y\mid Y>u)$. On the other hand, if $Y$ represents an annual
maxima, then the generalized extreme value distribution provides a good fit
(Fisher and Tippett, 1928). In hydrology and climate science, risk is often
assessed as the $T$-year return level $Q^{T}$, that is, the size of an event
that is exceeded in average once every $T$ years. If $Y$ represent a quantity
with $n_{Y}$ independent recordings per year (e.g., $n_{Y}=365$ for daily data
and $n_{Y}=1$ for annual maxima), then the $T$-year return level is the
quantile $Q^{T}=Q(1-1/(n_{Y}T))$.

Figure 1 shows the river catchment of a gauging station on the Aare river in
Bern, Switzerland. It is part of the Aare–Rhein basin where flooding is a
major economic and safety concern (Andres et al., 2021). The Swiss Federal
Office for the Environment (FOEN) provides recordings of daily average
discharges throughout the country. For risk assessment, they report the
100-year return levels using the GEV method for annual maxima and the GPD for
large daily discharges. The horizontal lines in Figure 2 show estimates of
this return level at the Bernese station based on data from the years
1930–1958. Both methods give very similar results.

Figure 1: Topographic map of water catchment of the gauging station in
Bern–Schönau (62) on the Aare in Switzerland. Another gauging station upstream
in Gsteig (43) on the Lütschine river and six meteorological stations with
precipitation measurements (triangles) are also shown.

The disadvantage of such an unconditional approach is twofold. First, the
return level is static and unable to reflect changes in the size of extreme
floods over time, which can occur for instance due to climate change. For the
Bernese station, for instance, a structural break has been observed in the
nineties without a clearly defined cause111see flood report of the FOEN at
https://www.hydrodaten.admin.ch/en/2135.html, making classical extreme value
modeling challenging. Second, while the return level $Q^{T}$ is relevant for
construction of long-term flood infrastructure, it can not be used to assess
the risk of flooding at a given day. Such forecasting of extreme events is
crucial for early warning systems. Indeed, the probability of exceeding on a
particular day a given high threshold, say the (constant) 100-year return
level, depends on many covariates $\mathbf{X}$ such as the river flows
upstream and precipitation in the catchments during the preceding days and
weeks.

In this paper we therefore advocate a conditional version of return levels
defined as the conditional quantile of $Y$ given a vector of observed
covariates $\mathbf{X}=\mathbf{x}$, that is,

$Q_{\mathbf{x}}(\tau)=F^{-1}_{Y\mid\mathbf{X}=\mathbf{x}}(\tau).$ (2)

The interpretation of such a conditional $T$-year return level
$Q^{T}_{\mathbf{x}}=Q_{\mathbf{x}}(1-1/(n_{Y}T))$ is different from the
unconditional return level $Q^{T}$. Since $Q^{T}_{\mathbf{x}}$ depends on the
exact configuration of the covariates $\mathbf{x}$, one can see a conditional
$T$-year return level as the size of an event that is exceeded in average once
every $T$ years, if the covariate vector $\mathbf{X}$ of all observations of
$Y$ had exactly the same value $\mathbf{x}$. A more precise interpretation is
to see $Q^{T}_{\mathbf{x}}$ as the value that is exceeded in the next time
step with probability $1/(n_{Y}T)$, and we refer to it as the conditional
quantile with return period $T$ years; for a comprehensive discussion of
return levels and quantiles see Bücher and Zhou (2018).

The top panel of Figure 2 shows one-day-ahead forecasts of such conditional
100-year quantiles $Q^{100}_{\mathbf{x}}$ for the Bernese river data based the
method of this paper fitted to the years 1930–1958. We see that, as opposed to
the unconditional return level $Q^{100}$, the conditional quantile changes
from day to day depending on past precipitation and river flows. In fact, on
August 21, the day before the first exceedance of the unconditional $Q^{100}$
during the serious flood of August 2005 in Bern, the size of the conditional
100-year event predicted for the next day was much higher than on other days,
and the forecasted conditional probability of such an exceedance (bottom panel
of Figure 2) was 920 times larger than the static 100-year probability. Both
outputs can be used as triggers for early warnings and additional flood
management measures. In particular, the days when an exceedance above
$Q^{100}$ is likely can be effectively pinpointed thanks to their temporal
sparsity in the probability forecast.

Figure 2: Top: Daily average discharge (points) at the Bern–Schönau station
(62) and one-day-ahead EQRN forecasts of conditional 100-year quantiles (solid
line) during the 2005 flood. Horizontal lines show unconditional $Q^{100}$
based on GEV (dashed) and GPD (dotted). Bottom: one-day-ahead EQRN forecast of
conditional probability of exceeding the GEV estimated $Q^{100}$ as a ratio to
the unconditional probability. The vertical line indicates August 22, the day
of the first exceedance.

Forecasting extreme events is notoriously difficult due to the low occurrence
probabilities involved. This is particularly true for hydrological models,
which are typically used by national agencies, and which do not use explicit
tail extrapolation. In the aftermath of the 2005 flood, the FOEN published a
detailed analysis of the internal forecasting procedures during this event
(Bezzola and Hegg, 2007). The report showed that too late warnings lead to
more severe consequences since forecasters did not trust the predicted
precipitation amounts during this extreme scenario and underestimated the
flood risk. This gives another motivation why our statistical model and its
output in Figure 2 could have been a helpful tool during this flood.

Conditional quantiles (2) are studied in the field of quantile regression
where many flexible methods exist (Meinshausen, 2006; Cannon, 2011; Zhang et
al., 2019; Athey et al., 2019). While they work well for quantile levels
within the data range, they break down for extreme values of $\tau$ close to
1. Such extreme quantile regression relies on extrapolation results as in (1),
where extreme value parameters such as the scale $\sigma(\mathbf{x})$ and
shape $\xi(\mathbf{x})$ may be modeled as functions of the covariates through
linear models (Wang et al., 2012; Li and Wang, 2019), generalized additive
models (Chavez-Demoulin and Davison, 2005; Youngman, 2019) or kernel methods
(Daouia et al., 2011; Gardes and Stupfler, 2019; Velthoen et al., 2019). In
order to overcome limitations of additive and kernel based methods in higher
dimensions, more recently, flexible tree-based methods have been combined with
the GPD extrapolation for predicting extreme conditional quantiles on complex
data (Velthoen et al., 2021; Gnecco et al., 2022). Tree-based methods have the
advantage of requiring little tuning for good prediction performance. However,
they can not incorporate additional structure of the data as encountered in
time series or spatial applications.

The goal of our work is to combine ideas from extreme value theory and machine
learning to propose an extreme quantile regression model that has the ability
to extrapolate in direction of the response $Y$ and to model complex covariate
dependencies in the predictors $\mathbf{X}$. We propose an extreme quantile
regression network (EQRN) that models covariate-dependent GPD parameters
$\sigma(\mathbf{x})$ and $\xi(\mathbf{x})$ as outputs of a neural network.
Conditional quantile estimates at the desired extreme level are then readily
derived from the estimated conditional tail distribution. Neural networks are
known for their ability to model complex dependencies and to approximate any
measurable function arbitrarily well (Hornik, 1991). The second advantage is
versatility. The deep learning literature is rich on network architectures,
activation functions and regularization methods. In particular, convolutions
produce shift-invariant models for covariates with spatial dependencies, and
recurrent architectures provide models for sequentially dependent observations
such as time series. In our application of flood forecasting with time
dependent data, recurrent neural networks (Werbos, 1988; Elman, 1990) are of
particular interest.

The main output of our EQRN model for sequentially dependent data is the one-
day-ahead risk forecast, either in terms of conditional quantiles or
exceedances probabilities. Thanks to the GPD approximation, the method can
extrapolate beyond the range of the data as illustrated in Figure 2. The
strength of our recurrent model lies in the ability to exploit information
from multiple covariates and capture complex time dependence. It is therefore
an effective early warning tool even for unprecedented, record-shattering
events, like the 2005 floods. This is of particular importance in a non-
stationary system where climate change makes extreme events increasingly
likely (Fischer et al., 2021).

The paper is organized as follows. In Section 2 we provide background on
quantile regression, extreme value theory and neural networks. We propose our
EQRN model for both independent and sequentially dependent data in Section 3.
Section 4 contains a simulation study to assess the performance of our
approach in comparison to existing methods. In Section 5 we describe the Swiss
river data, apply our methodology to forecasting of flood risk and discuss the
implications of the results. Section 6 concludes with a brief discussion.

##  2 Background

###  2.1 Quantile Regression

In the classical quantile regression setup, we observe an independent and
identically distributed sample
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ of the random
vector $(\mathbf{X},Y)$, where $Y$ is the real-valued response variable and
$\mathbf{X}$ is a vector of $p$ covariates (or predictors). One aims at
predicting the conditional quantile $Q_{\mathbf{x}}(\tau)$ defined in (2) of
$Y$ given $\mathbf{X}=\mathbf{x}$, for some predictor value
$\mathbf{x}\in\mathbb{R}^{p}$ and probability level $\tau\in(0,1)$ of
interest.

Analogously to regression that minimizes the mean squared error, quantile
regression minimizes the quantile loss,

$Q_{\mathbf{x}}(\tau)=\operatorname*{arg\,min}_{q}\mathbb{E}[\rho_{\tau}(Y-q)\mid\mathbf{X}=\mathbf{x}],$
(3)

where $\rho_{\tau}(t):=t(\tau-\mathds{1}_{\\{t<0\\}})$ is the quantile check
function (Koenker and Bassett, 1978). Many parametric and non-parametric
quantile regression models exist, including linear models (Chernozhukov,
2005), random forests (Athey et al., 2019) and neural networks (Cannon, 2011;
Zhang et al., 2019). They yield a conditional quantile estimate by minimizing
the empirical quantile loss over the training sample $\mathcal{D}$, that is,

$\hat{Q}_{\mathbf{x}}(\tau)=\operatorname*{arg\,min}_{q_{\tau}\in\mathcal{M}}\dfrac{1}{n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-q_{\tau}(\mathbf{x}_{i})),$
(4)

where $\mathcal{M}$ is the set of possible quantile functions
$q_{\tau}(\cdot)$ characterized by the model.

Classical methods for quantile regression that rely on the quantile loss (3)
perform well for “moderate” probability levels $\tau$. To define what that
means exactly, we typically let $\tau_{n}$ depend on the sample size $n$. The
expected number of exceedances of $y_{i}$ over the respective conditional
quantile $Q_{\mathbf{x}_{i}}(\tau_{n})$, $i=1,\dots,n$, is given by
$n(1-\tau_{n})$. With moderately extreme, or intermediate, we refer to a
sequence $\tau_{n}\to 1$ with $n(1-\tau_{n})\to\infty$, meaning that the
quantile goes to the upper endpoint of the distribution but there are more and
more exceedances with growing sample size $n$. On the other hand, we call a
quantile level $\tau_{n}\to 1$ extreme if $n(1-\tau_{n})\to c\in[0,\infty)$,
that is, there are finitely many, or possibly zero, exceedances over
$Q_{\mathbf{x}_{i}}(\tau_{n})$ in the sample. In this situation, classical
quantile regression methods do not perform well due to the scarcity of
observations in the tail of the response.

The left panel of Figure 3 illustrates this issue for a sample
$y_{1},\dots,y_{n}$ with $n=1000$ and no covariates. The dashed line shows for
different quantile levels $\tau_{n}$ the empirical quantile estimates,
obtained by solving the respective quantile loss function without covariates.
It can be seen that as soon as the number of exceedances $n(1-\tau_{n})<1$,
that is, $\tau_{n}>99.9\%$, there is a significant bias compared to the true
quantiles (solid line). The reason is that the empirical estimates can not
predict higher than the largest observation. When covariates are present, this
issue persists since quantile regression relies on solving the empirical
quantile loss (4).

In the sequel we omit the dependence on $n$ and write $\tau$ instead of
$\tau_{n}$. Intermediate quantile levels will be denoted by $\tau_{0}$.

###  2.2 Generalized Pareto Distribution

In order to predict well on extreme quantiles, a method should rely on
asymptotic results from extreme value theory for accurate extrapolation beyond
the sample. In particular, we rely on the generalized Pareto distribution
(GPD). In the presence of covariates, it arises as an approximation of the
tail of the distribution of $Y\mid\mathbf{X}=\mathbf{x}$. More precisely, we
use a conditional version of (1),

$\displaystyle\mathbb{P}(Y>y\mid\mathbf{X}=\mathbf{x})\approx(1-\tau_{0})\left(1+\xi(\mathbf{x})\frac{y-u(\mathbf{x})}{\sigma(\mathbf{x})}\right)_{+}^{-1/\xi(\mathbf{x})},\quad
y>u(\mathbf{x}),$ (5)

where the threshold $u(\mathbf{x})$ is chosen as an intermediate quantile
$Q_{\mathbf{x}}(\tau_{0})$ at level $\tau_{0}\in(0,1)$ close to 1, and the
shape $\xi(\mathbf{x})\in\mathbb{R}$ and scale $\sigma(\mathbf{x})>0$ depend
on the covariates; here we omit the dependence of $\sigma(\mathbf{x})$ on the
intermediate level $\tau_{0}$ in the notation. This approximation holds under
weak conditions on the tail of $Y\mid\mathbf{X}=\mathbf{x}$; see (Balkema and
de Haan, 1974; Pickands III, 1975) for the precise statement. The shape
parameter $\xi(\mathbf{x})$ is important since it encodes the tail heaviness
of the response: if it is positive, the response has a heavy-tailed
distribution such as Pareto or Student-$t$; if it is zero, the response is
light-tailed such as a Gaussian or exponential; if it is negative then the
response has a finite upper endpoint.

In order to predict an extreme quantile at level $\tau>\tau_{0}$ from
approximation (5), we can invert this expression to find

$Q_{\mathbf{x}}(\tau):=Q_{\mathbf{x}}(\tau_{0})+\frac{\sigma(\mathbf{x})}{\xi(\mathbf{x})}\left[\left(\dfrac{1-\tau_{0}}{1-\tau}\right)^{\xi(\mathbf{x})}-1\right].$
(6)

This shows that an estimate $\hat{Q}_{\mathbf{x}}(\tau)$ of an extreme
quantile requires estimates of the intermediate quantile
$\hat{Q}_{\mathbf{x}}(\tau_{0})$ and of the conditional GPD parameters
$\hat{\xi}(\mathbf{x})$ and $\hat{\sigma}(\mathbf{x})$ as functions of the
predictor vector. For the intermediate quantile function we can use any of the
existing methods for quantile regression since they work well for this
moderate quantile level as discussed above. Estimation of the GPD parameters
can be done by specifying a parametric or non-parametric model. We will use
neural networks for this purpose, which are introduced in the next section.

The green line in Figure 3 shows estimates $\hat{Q}(\tau)$ for different
quantile levels $\tau$ using the approximation (6) without covariate
dependence, with empirical intermediate quantile at $\tau_{0}=90\%$ and GPD
parameters estimated with maximum likelihood. It can be seen that the
extrapolation solves the bias issue of empirical methods.

$\mathbf{x}$$\mathbf{x}^{(1)}$$\mathbf{x}^{(2)}$$\mathbf{x}^{(L)}$$g_{\mathcal{W}}(\mathbf{x})$$y$$\ell$

Figure 3: Left: true $\tau$-quantiles (solid line) compared to empirical
estimates (dashed line) and GPD based estimates (dotted line) for moderate to
extreme probability levels (log-scale) for sample size 1,000. Estimates are
averages over 100 trials. Right: multi-layer perceptron flowchart from input
$\mathbf{x}$ to output $g_{\mathcal{W}}(\mathbf{x})$, with loss function
$\ell$ and corresponding response $y$.

###  2.3 Neural Networks and Conditional Density Estimation

The literature on neural networks is vast and existing methods are being
improved constantly. We concentrate in this section on well-established
techniques that are most relevant for our purpose of modeling extreme
quantiles.

A multi-layer perceptron (MLP) or fully-connected feed-forward neural network
model is a parametric family of non-linear functions
$g_{\mathcal{W}}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{q}$ that map a
$p$-dimensional input $\mathbf{x}$ to a $q$-dimensional output by

$\mathbf{x}\mapsto\mathbf{x}^{(L+1)},\;\text{ with
}\;\mathbf{x}^{(l)}=\sigma^{l}\left(W^{l}\mathbf{x}^{(l-1)}+b^{l}\right)\;\forall
l=1,\ldots,L+1,$ (7)

where $\mathbf{x}^{(0)}=\mathbf{x}$. The number of hidden layers
$L\in\mathbb{N}$, the hidden layer dimensions
$h_{1},\ldots,h_{L}\in\mathbb{N}$ and the choice of activation functions
$\sigma^{l}:\mathbb{R}^{h_{l}}\rightarrow\mathbb{R}^{h_{l}},l=1,\ldots,L+1$
(with $h_{0}=p$ and $h_{L+1}=q$) are hyperparameters that need be chosen for
instance by cross-validation. The set of trainable parameters to be inferred
from data cotains all weights and bias terms of the network, that is,
$\mathcal{W}=\left\\{(W^{l},b^{l});l=1,\ldots,L+1\right\\}$, with
$W^{l}\in\mathbb{R}^{h_{l}\times h_{l-1}}$ and $b^{l}\in\mathbb{R}^{h_{l}}$.
Figure 3 shows a schematic illustration of the transformations inside the MLP.

In the general setting, $p$ is the number of features or covariates considered
in the model and $q$ depends on the task at hand. In order to train a model, a
loss function $\ell:\mathbb{R}\times\mathbb{R}\to[0,\infty)$ is required that
maps a tuple $(y,g_{\mathcal{W}}(\mathbf{x}))$ of response and prediction to a
positive number quantifying their discrepancy. Common tasks include mean
regression with $q=1$ and squared error loss, quantile regression with $q=1$
and quantile loss, and classification with $q$ equal to the number of possible
classes and cross-entropy as loss.

For conditional density estimation we suppose that $Y$ follows a distribution
with parametric probability density $f_{Y}(\cdot;{\theta})$ and parameter
${\theta}={\theta}(\mathbf{x})$ depending on the vector
$\mathbf{X}=\mathbf{x}$. Conditional density estimation networks are neural
networks that aim at outputting conditional estimates
$g_{\mathcal{W}}(\mathbf{x})={\theta}(\mathbf{x})$ based on realizations of
$\mathbf{X}=\mathbf{x}$ as input (e.g., Cannon, 2012). In this setting, $p$ is
the dimension of $\mathbf{X}$ and $q$ is the dimension of ${\theta}$. The loss
function is the deviance or negative log-likelihood loss
$\ell(y,{{\theta}}(\mathbf{x}))=-\log f_{Y}(y;{{\theta}}(\mathbf{x}))$.

To train a neural network on the training data set
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ we find the
optimal parameter values minimizing the average empirical loss, that is,

$\hat{\mathcal{W}}\in\operatorname*{arg\,min}_{\mathcal{W}}\frac{1}{n}\sum_{i=1}^{n}\ell(y_{i},g_{\mathcal{W}}(\mathbf{x}_{i})).$
(8)

This is generally achieved via backpropagation using mini-batches and
optimization algorithms such as the well-performing gradient descent variants
(Kingma and Ba, 2014; Tieleman and Hinton, 2012; Duchi et al., 2011). Since
neural networks are typically overparameterized, overfitting has to be
prevented with regularization methods such as $L_{2}$ weight penalties for
narrow networks and dropout (Srivastava et al., 2014) for deeper structures.
As the optimization problem (8) is often non-convex, local-minima convergence
is an issue. Restricting training to a subset of $\mathcal{D}$ and keeping the
rest to track the validation loss at the end of each epoch helps to avoid
local minima by learning rate decay. Restarting training with different
initializations and keeping the best fit in terms of validation loss often
leads to lower minima. Early stopping based on the validation loss is another
measure against overfitting. The final validation loss is used for model
selection and the choice of optimal hyperparameters; for more details on
fitting of neural networks see Goodfellow et al. (2016).

When observations are dependent in space or time, generalizations of the MLP
exist to account for these particular structures. Convolutional neural
networks use the neighborhood information in images or spatial observations to
render neural networks more parsimonious and effective (LeCun et al., 2015).
We concentrate here on methods for sequential dependence that arises typically
in time series $\left\\{(\mathbf{X}_{t},Y_{t})\right\\}_{t=1}^{T}$. For this
type of data, recurrent architectures of the network allow to capture
dependence between observations. A simple recurrent neural network (RNN) layer
(Werbos, 1988; Elman, 1990) takes as input a vector $\mathbf{x}_{t}$ and
outputs the hidden recurrent state vector

$\mathbf{h}_{t}=\tanh(W_{\rm xh}\mathbf{x}_{t}+W_{\rm
hh}\mathbf{h}_{t-1}+\mathbf{b}_{\rm h}),$

depending both on $\mathbf{x}_{t}$ and the hidden state $\mathbf{h}_{t-1}$
recursively resulting from the previous inputs $\mathbf{x}_{t-1}$ and
$\mathbf{h}_{t-2}$ in the sequence. Here and in the sequel, the bias vectors
$b_{\cdot}$ and weight matrices $W_{\cdot\cdot}$, indexed by the input and
output variables, are the trainable parameters. This model has then been
improved by the addition of a gating cell state $\mathbf{c}_{t}$, to avoid
vanishing gradient issues, as well as a forget gate $\mathbf{f}_{t}$, an input
gate $\mathbf{i}_{t}$ and and output gate $\mathbf{o}_{t}$, to control both
short- and long-term dependencies in the sequence. This yields the long short-
term memory (LSTM) layer (Hochreiter and Schmidhuber, 1997; Gers et al., 2000,
2003; Jozefowicz et al., 2015)

$\begin{array}[]{ll}\mathbf{i}_{t}=\sigma(W_{\rm xi}\mathbf{x}_{t}+W_{\rm
hi}\mathbf{h}_{t-1}+\mathbf{b}_{\rm i}),&\mathbf{f}_{t}=\sigma(W_{\rm
xf}\mathbf{x}_{t}+W_{\rm hf}\mathbf{h}_{t-1}+\mathbf{b}_{\rm f}),\\\
\mathbf{g}_{t}=\tanh(W_{\rm xg}\mathbf{x}_{t}+W_{\rm
hg}\mathbf{h}_{t-1}+\mathbf{b}_{\rm g}),&\mathbf{o}_{t}=\sigma(W_{\rm
xo}\mathbf{x}_{t}+W_{\rm ho}\mathbf{h}_{t-1}+\mathbf{b}_{\rm o}),\\\
\mathbf{c}_{t}=\mathbf{f}_{t}\odot\mathbf{c}_{t-1}+\mathbf{i}_{t}\odot\mathbf{g}_{t},&\mathbf{h}_{t}=\mathbf{o}_{t}\odot\tanh(\mathbf{c}_{t}),\end{array}$
(9)

where $\sigma(z)=1/(1+\exp(-z))$ is the sigmoid activation and $\odot$ is the
Hadamard (or componentwise) product; see Figure 4 for an illustration. The
common dimension of the vectors defined in (9) is a hyperparameter of the
layer. The input of this layer
$\tilde{\mathbf{x}}_{t}:=(\mathbf{x}_{t-s},\ldots,\mathbf{x}_{t-1})$ can
include predictors from the past to model longer dependencies, where
$s\in\mathbb{N}$ determines the time horizon.

The LSTM model has been simplified by Cho et al. (2014) into the gated
recurrent unit (GRU) layer, which has become a popular alternative. A multi
layer recurrent network is obtained by considering the $\mathbf{h}_{t}$ as a
sequence of inputs for the following recurrent layers; see Figure S.1 in
Supplementary Material S.1. Usually a fully connected layer as in (7) is used
to map the hidden state of the final recurrent layer to the network output
$\tilde{g}_{\mathcal{W}}(\tilde{\mathbf{x}}_{t})$.

$\mathbf{x}_{t-1}$$\mathbf{x}_{t-2}$$\mathbf{x}_{t-3}$$\mathbf{x}_{t-s}$LSTMLSTMLSTMLSTM$\mathbf{h}_{t-1}$$\mathbf{h}_{t-2}$$\mathbf{h}_{t-3}$$\mathbf{c}_{t-1}$$\mathbf{c}_{t-2}$$\mathbf{c}_{t-3}$$\tilde{g}_{\mathcal{W}}(\tilde{\mathbf{x}}_{t})$$y_{t}$$\ell$
Figure 4: Single-layer LSTM network flowchart from input
$\tilde{\mathbf{x}}_{t}:=(\mathbf{x}_{t-s},\ldots,\mathbf{x}_{t-1})$ to output
$\tilde{g}_{\mathcal{W}}(\tilde{\mathbf{x}}_{t})$, with loss evaluation. The
LSTM cells represent the transformation in (9).

##  3 Extreme Quantile Regression Neural Networks

In this section we propose a new methodology that combines the extrapolation
power of the GPD model with the high-dimensional predictor space capabilities
and flexibility of neural networks to obtain accurate estimates for quantile
functions $Q_{\mathbf{x}}(\tau)$ at extreme levels $\tau$. Let
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ be the
training data set. Estimation of conditional extreme quantiles
$\hat{Q}_{\mathbf{x}}(\tau)$ using (6) requires estimators for the
intermediate quantile function $\hat{Q}_{\mathbf{x}}(\tau_{0})$ with
$\tau_{0}<\tau$ and the GPD parameter $\sigma(\mathbf{x})$ and
$\xi(\mathbf{x})$. It is customary to proceed in two steps.

First, we model the intermediate quantile at level $\tau_{0}$ using classical
quantile regression methods. We then define the conditional exceedances

$z_{i}:=y_{i}-\hat{Q}_{\mathbf{x}_{i}}(\tau_{0}),\quad
i\in\mathcal{I}:=\\{i=1,\ldots,n\,:\;y_{i}>\hat{Q}_{\mathbf{x}_{i}}(\tau_{0})\\}.$

The intermediate probability $\tau_{0}$ should be chosen low enough to allow
for stable estimation of $Q_{\mathbf{x}}(\tau_{0})$ with classical empirical
methods, but high enough for the approximation in (5) to be accurate, so that
the exceedances $z_{i}$ are approximate samples of a GPD.

In the second step we estimate the GPD parameters $\sigma(\mathbf{x})$ and
$\xi(\mathbf{x})$ based on the set of exceedances $z_{i}$, $i\in\mathcal{I}$.
Modelling these parameters directly in the extrapolation formula (6) may lead
to strong dependence between the estimates and numerical instabilities. We
therefore rely on an orthogonal reparametrization that has diagonal Fisher
information matrix. As for the standard asymptotic GPD likelihood properties,
the latter is well-defined for the GPD model when $\xi(\mathbf{x})>-0.5$, and
the reparametrization

$(\sigma(\mathbf{x}),\xi(\mathbf{x}))\mapsto(\nu(\mathbf{x}),\xi(\mathbf{x})),\qquad\nu(\mathbf{x}):=\sigma(\mathbf{x})(\xi(\mathbf{x})+1),$

yields the desired orthogonality (Cox and Reid, 1987; Chavez-Demoulin and
Davison, 2005). In our experiments this reparametrization significantly
improves stability and convergence in every considered setting.

In this section we propose a flexible neural network model for the
orthogonalized GPD parameters $\nu(\mathbf{x};\mathcal{W})$ and
$\xi(\mathbf{x};\mathcal{W})$, where $\mathcal{W}$ denotes the collection of
all model parameters. This can be seen as conditional density estimation with
output dimension $q=2$ where the parametric family is the GPD model with
parameters $\theta=(\nu,\xi)$ depending on the covariate
$\mathbf{X}=\mathbf{x}$. In general, an estimate of the model parameters
$\hat{\mathcal{W}}$ is thus obtained as a minimizer of the GPD deviance loss
over the training exceedances

$\hat{\mathcal{W}}\in\operatorname*{arg\,min}_{\mathcal{W}}\sum_{i\in\mathcal{I}}\ell_{\rm
OGPD}\\{z_{i};\hat{\nu}(\mathbf{x}_{i};\mathcal{W}),\hat{\xi}(\mathbf{x}_{i};\mathcal{W})\\},$
(10)

where the deviance or negative log-likelihood of the GPD in terms of the
orthogonal reparametrization is

$\ell_{\rm
OGPD}(z;\nu,\xi)=\left(1+\dfrac{1}{\xi}\right)\log\left\\{1+\xi\dfrac{(\xi+1)z}{\nu}\right\\}+\log(\nu)-\log(\xi+1).$
(11)

In the next two sub-sections we discuss the details of the model for
independent observations and for time series data with sequential dependence,
respectively.

###  3.1 Independent Observations

We first consider the case where the training data
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ is a set of
independent, identically distributed observations of $(\mathbf{X},Y)$. The
goal in this case is the estimation of the conditional quantile
$Q_{\mathbf{x}}(\tau)$ for a predictor value $\mathbf{X}=\mathbf{x}$ at an
extreme level $\tau>0$. For the first step of estimating the intermediate
quantile function, in principle any classical quantile regression method can
be used. In order to avoid overfitting and obtain unbiased generalization
error estimates from the training set $\mathcal{D}$, the predicted
$\hat{Q}_{\mathbf{x}_{i}}(\tau_{0})$, $i=1,\dots,n$, should be constructed out
of training sample. This is achievable in two ways. Using bagging methods such
as generalized random forests (Athey et al., 2019) is a convenient choice
since they allow for out-of-bag predictions where only a single fit on
$\mathcal{D}$ is required. For other methods such as quantile regression
neural networks (Cannon, 2011), out of training sample predictions can be
obtained in a foldwise manner similar to cross-validation. In the sequel, we
assume that the intermediate quantile, and thus the exceedances, are given.

In the second step, we propose to model the orthogonalized GPD parameters
$\nu(\mathbf{x})$ and $\xi(\mathbf{x})$ by a fully-connected feed-forward
neural network with parameter vector $\mathcal{W}$ and deviance loss function
as in (10); see Section 2.3 for details. Choices of the network structure such
as the number of neurons, the number of layers and activation functions, are
hyperparameters, denoted by $\Theta$, to be selected. We provide sensible
default values in our implementation, but one can also choose them in a data-
driven way based on a validation set.

The only restrictions are on the output activation functions, since
${\nu}(\mathbf{x})$ should be strictly positive. We find the exponential
function or the ${\rm SELU}$ activation (Klambauer et al., 2017) shifted above
zero to be good choices. Regarding the output activation for ${\xi}$, no
strict restrictions apply and the identity would be the natural choice.
However, standard likelihood regularity properties are not satisfied for the
GPD model when $\xi\leq-0.5$, which very rarely occurs in practice. We observe
that smoothly restricting the shape estimates, for example between $-0.5$ and
$0.7$ with the activation $x\mapsto 0.6\tanh{(x)}+0.1$, helps improving
training stability. This avoids aberrant $\xi$ estimates in early stages of
the training and still covers almost all practical cases. In many situations
it is reasonable to assume that only the scale $\nu(\mathbf{x})$ varies
locally but $\xi(\mathbf{x})\equiv\xi$ is constant (e.g., Kinsvater et al.,
2016). This can be achieved by restricting the network so that the shape
output only depends on a bias term.

Algorithm 1 summarizes our extreme quantile regression network (EQRN) for
independent observations, which takes the intermediate quantiles and training
data as input and outputs the extreme quantile at a desired test predictor
value $\mathbf{x}\in\mathbb{R}^{p}$ and level $\tau>\tau_{0}$. Optionally, the
conditional GPD parameters can also be obtained.

Algorithm 1 EQRN for independent observations

The tuning parameters $\Theta$ for the conditional GPD density estimation
network $g_{\mathcal{W}}$ and the intermediate quantile model
$\hat{Q}_{\cdot}(\tau_{0})$ capable of out of sample prediction are pre-
specified. The training data
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ and test
covariates $\mathbf{x}$ are observed. Let $\tau\in(\tau_{0},1)$ be the desired
probability level.

1:procedure EQRN-Fit($\mathcal{D}$, $\Theta$, $\hat{Q}_{\cdot}(\tau_{0})$)

2:
$\mathcal{I}\leftarrow\\{i=1,\ldots,n\,:\;y_{i}>\hat{Q}_{\mathbf{x}_{i}}(\tau_{0})\\}$

3: $z_{i}\leftarrow y_{i}-\hat{Q}_{\mathbf{x}_{i}}(\tau_{0})\;\;\forall
i\in\mathcal{I}$

4: $\mathcal{T},\mathcal{V}\leftarrow$ RandomValidationSplit($\mathcal{I}$)
$\triangleright$ If no validation:
$\mathcal{T}=\mathcal{I},\mathcal{V}=\emptyset$

5: $\hat{\mathcal{W}}\leftarrow$ InitializeNetworkWeights($\Theta$)

6: for $e=1$ to maximum number of epochs $E$ do

7: for all $\mathcal{B}\in$ GetMiniBatches($\mathcal{T}$) do

8: $\\{(\hat{\nu}_{i},\hat{\xi}_{i})\\}_{i\in\mathcal{B}}\leftarrow$
$g_{\hat{\mathcal{W}}}(\mathbf{x}_{\mathcal{B}},\hat{Q}_{\mathbf{x}_{\mathcal{B}}}(\tau_{0}))$

9: $\ell\leftarrow\sum_{i\in\mathcal{B}}\ell_{\rm
OGPD}(z_{i},\hat{\nu}_{i},\hat{\xi}_{i})/\left\lvert\mathcal{B}\right\rvert$

10: $\hat{\mathcal{W}}\leftarrow$ BackPropUpdate($\ell$, $\hat{\mathcal{W}}$,
$\mathbf{x}_{\mathcal{B}}$, $\hat{Q}_{\mathbf{x}_{\mathcal{B}}}(\tau_{0})$,
$\Theta$)

11: stop if $\mathcal{V}\neq\emptyset$ and
LossNotImproving($\hat{\mathcal{W}}$, $\mathbf{x}_{\mathcal{V}}$,
$\hat{Q}_{\mathbf{x}_{\mathcal{V}}}(\tau_{0})$, $z_{\mathcal{V}}$)

12: output $\hat{\mathcal{W}}$

13:procedure EQRN-Predict($\mathbf{x}$, $\tau$, $\hat{\mathcal{W}}$,
$\hat{Q}_{\cdot}(\tau_{0})$)

14: $\\{\hat{\nu}(\mathbf{x}),\hat{\xi}(\mathbf{x})\\}\leftarrow$
$g_{\hat{\mathcal{W}}}(\mathbf{x},\hat{Q}_{\mathbf{x}}(\tau_{0}))$

15:
$\hat{\sigma}(\mathbf{x})\leftarrow\hat{\nu}(\mathbf{x})/\\{\hat{\xi}(\mathbf{x})+1\\}$

16: compute $\hat{Q}_{\mathbf{x}}(\tau)$ w.r.t. $\hat{\sigma}(\mathbf{x})$,
$\hat{\xi}(\mathbf{x})$, $\hat{Q}_{\mathbf{x}}(\tau_{0})$, $\tau$ and
$\tau_{0}$ using equation (6)

17: output $\hat{Q}_{\mathbf{x}}(\tau)$, and optionally
$\\{\hat{\sigma}(\mathbf{x}),\hat{\xi}(\mathbf{x})\\}$

The conditional GPD estimation in the second step relies on the exceedances to
ensure that only information of the tail is used for extrapolation. There may
however be residual information in the moderately extreme observations that
should not be discarded. We propose to use the intermediate quantiles
$\hat{Q}_{\mathbf{x}_{i}}(\tau_{0})$, as an additional feature in the
conditional density estimation. This feature engineering seems to consistently
and significantly improve the accuracy of the final prediction
$\hat{Q}_{\mathbf{x}}(\tau)$ on test data in our simulations. This idea is to
some degree related to stacked learning (Breiman, 1996; Wolpert, 1992) and is
achieved by considering the $(p+1)$-dimensional vector
$(\mathbf{x}_{i},\hat{Q}_{\mathbf{x}_{i}}(\tau_{0}))$, $i=1,\dots,n,$ instead
of $\mathbf{x}_{i}$ as predictors for the network input. To avoid overfitting
it is again important that the intermediate quantile estimates are constructed
out of training sample. This new feature also improves other extreme quantile
regression methods such as the GBEX model (Velthoen et al., 2021), as observed
in our simulations in Section 4.

Classical backpropagation and optimization procedures are performed to find
$\hat{\mathcal{W}}$. To avoid overfitting and local minima convergence, the
validation loss can be tracked as discussed in Section 2.3. The
hyperparameters $\Theta$ of this network include choices of optimization
algorithms and regularization, for instance. More details on the function
calls in the algorithm can be found in Supplementary Material S.2.

Since the true quantile $Q_{\mathbf{x}}(\tau)$ is unknown in real-world data,
we can not assess the performance of $\hat{Q}_{\mathbf{x}}(\tau)$ with metrics
such as mean squared or absolute errors. As illustrated in Section 2.1, the
quantile loss is also unreliable due to the data scarcity at extreme
quantiles. We therefore choose the final validation loss based on the GPD
deviance to compare different choices of hyperparameters, as it is the most
reliable surrogate metric.

###  3.2 Sequential Dependence

In many applications the observations are not independent but display
sequential dependence such as in time series. In this case, we denote the
training data by
$\mathcal{D}=\left\\{(\mathbf{x}_{t},y_{t})\right\\}_{t=1}^{T}$, which are
observed sequentially from a time series
$\left\\{(\mathbf{X}_{t},Y_{t})\right\\}_{t=1}^{T}$. The goal here is
different than in the case of independent observations. Indeed, we would like
to predict as well as possible high quantiles of the response $Y_{u}$ at some
time point $u$ one step in the future based on all past information
$\tilde{\mathbf{X}}_{u}:=\left\\{(\mathbf{X}_{t},Y_{t})\right\\}_{t<u}$.
Therefore the target is

$Q_{\tilde{\mathbf{x}}_{u}}(\tau):=F^{-1}_{Y_{u}\mid\tilde{\mathbf{X}}_{u}=\tilde{\mathbf{x}}_{u}}(\tau),$
(12)

where $\tilde{\mathbf{x}}_{u}:=\left\\{(\mathbf{x}_{t},y_{t})\right\\}_{t<u}$
are observations that are not necessarily part of the training set. In this
section we propose a recurrent neural network to solve this task. Several
principles are the same as in the case of independent observations, such as
the choice of output activation functions. We thus focus on the differences to
the independent case.

While in principle it is still possbile to use any classical quantile
regression method to model the intermediate conditional quantiles at level
$\tau_{0}$, we recommend using quantile regression neural networks (Cannon,
2011; Zhang et al., 2019) with a recurrent architecture. These models are
specifically designed for sequential dependence and can easily adapt to
varying sequence lengths of the input features $\tilde{\mathbf{X}}_{t}$. In
our experiments, recurrent quantile regression neural networks consistently
outperformed generalized random forests in the presence of sequential
dependence. Although a varying sequence length of input features is possible,
we restrict this length to a fixed horizon $s\ll T$, for computational
efficiency. Thus, for any time point $t$ we define its past by
$\tilde{\mathbf{x}}_{t}=\\{(\mathbf{x}_{j},y_{j})\\}_{j=t-s}^{t-1}$. For
simplicity, we denote our augmented training set by
$\tilde{\mathcal{D}}=\\{(\tilde{\mathbf{x}}_{t},y_{t})\\}_{t=s+1}^{T}$.

Algorithm 2 EQRN for sequential observations

The tuning parameters $\Theta$ for the recurrent conditional GPD density
estimation network $\tilde{g}_{\mathcal{W}}$, the intermediate quantile model
$\hat{Q}_{\cdot}(\tau_{0})$ capable of out of sample prediction and horizon
$s$ are pre-specified. The training data
$\tilde{\mathcal{D}}=\left\\{(\tilde{\mathbf{x}}_{t},y_{t})\right\\}_{t=s+1}^{T}$
and test covariates $\tilde{\mathbf{x}}_{u}$ are observed. Let
$\tau\in(\tau_{0},1)$ be the desired probability level.

1:procedure EQRN-Fit($\tilde{\mathcal{D}}$, $\tau_{0}$, $s$, $\Theta$,
$\hat{Q}_{\cdot}(\tau_{0})$)

2: $z_{t}\leftarrow
y_{t}-\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau_{0})\;\;\forall
t\in\mathcal{I}:=\\{t=s+1,\ldots,T\,:\;y_{t}>\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau_{0})\\}$

3: $\mathcal{T},\mathcal{V}\leftarrow$
SequentialValidationSplit($\mathcal{I}$)

4: $\hat{\mathcal{W}}\leftarrow$ InitializeRecurrentNetWeights($\Theta$)

5: for $e=1$ to maximum number of epochs $E$ do

6: for all $\mathcal{B}\in$ GetMiniBatches($\mathcal{T}$) do

7: $\\{(\hat{\nu}_{t},\hat{\xi}_{t})\\}_{t\in\mathcal{B}}\leftarrow$
$\tilde{g}_{\hat{\mathcal{W}}}(\tilde{\mathbf{x}}_{\mathcal{B}},\hat{Q}_{\tilde{\mathbf{x}}_{\mathcal{B}}}(\tau_{0}))$

8: $\ell\leftarrow\sum_{t\in\mathcal{B}}\ell_{\rm
OGPD}(z_{t},\hat{\nu}_{t},\hat{\xi}_{t})/\left\lvert\mathcal{B}\right\rvert$

9: $\hat{\mathcal{W}}\leftarrow$ BackPropUpdate($\ell$, $\hat{\mathcal{W}}$,
$\tilde{\mathbf{x}}_{\mathcal{B}}$,
$\hat{Q}_{\tilde{\mathbf{x}}_{\mathcal{B}}}(\tau_{0})$, $\Theta$)

10: stop if $\mathcal{V}\neq\emptyset$ and
LossNotImproving($\hat{\mathcal{W}}$, $\tilde{\mathbf{x}}_{\mathcal{V}}$,
$\hat{Q}_{\tilde{\mathbf{x}}_{\mathcal{V}}}(\tau_{0})$, $z_{\mathcal{V}}$)

11: output $\hat{\mathcal{W}}$

12:procedure EQRN-Predict($\tilde{\mathbf{x}}_{u}$, $\tau$,
$\hat{\mathcal{W}}$, $\hat{Q}_{\cdot}(\tau_{0})$)

13:
$\\{\hat{\nu}(\tilde{\mathbf{x}}_{u}),\hat{\xi}(\tilde{\mathbf{x}}_{u})\\}\leftarrow$
$\tilde{g}_{\hat{\mathcal{W}}}(\tilde{\mathbf{x}}_{u},\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau_{0}))$

14:
$\hat{\sigma}(\tilde{\mathbf{x}}_{u})\leftarrow\hat{\nu}(\tilde{\mathbf{x}}_{u})/\\{\hat{\xi}(\tilde{\mathbf{x}}_{u})+1\\}$

15: compute $\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau)$ w.r.t.
$\hat{\sigma}(\tilde{\mathbf{x}}_{u})$, $\hat{\xi}(\tilde{\mathbf{x}}_{u})$,
$\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau_{0})$, $\tau$ and $\tau_{0}$ using
equation (6)

16: output $\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau)$, and optionally
$\\{\hat{\sigma}(\tilde{\mathbf{x}}_{u}),\hat{\xi}(\tilde{\mathbf{x}}_{u})\\}$

Algorithm 2 summarizes our EQRN for sequential data. For the estimation of the
conditional GPD parameters, the main difference compared to the independent
model is the use of a recurrent architecture to capture the sequential nature
of the data. If validation splitting is used during training, the split should
preserve the sequential structure instead of being performed randomly. For the
use of the intermediate quantile as a feature, two approaches seem relevant.
The first one is to only use $\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau_{0})$ as a
separate additional input to $\tilde{\mathbf{x}}_{t}$ in the network. The
second approach is to also use past intermediate information by considering
$\\{(\mathbf{x}_{j},y_{j},\hat{Q}_{\tilde{\mathbf{x}}_{j}}(\tau_{0}))\\}_{j=t-s}^{t-1}$
instead of $\tilde{\mathbf{x}}_{t}$ as input features to model the GPD
parameters $\nu(\tilde{\mathbf{x}}_{t})$ and $\xi(\tilde{\mathbf{x}}_{t})$. We
prefer the second approach as it can pass more information to the tail model.
More details on the function calls in the algorithm can be found in
Supplementary Material S.2.

The training data consists of time points $\\{1,\dots,T\\}$, while the test
data uses information from time points $\\{u-s,\dots,u-1\\}$ to predict at
time $u$. These two intervals are typically disjoint when the model was fitted
in the past and is applied for prediction in the present. The prediction model
can of course be used to predict on the training data when $u\leq T$, but such
predictions might be overly precise since $y_{u}$ was used in the training
procedure.

##  4 Simulation Study

###  4.1 Setup

In this section we assess the accuracy of our EQRN model in predicting extreme
conditional quantiles on simulated data, and compare it to existing state-of-
the-art methods. The aim is to study a simplified version of the application,
which motivates the simulation setup and modeling choices. We focus on the
case of sequentially dependent data; a simulation study for independent data
can be found in Supplementary Material S.3.

The main competitors from the extreme value literature in terms of flexibility
are the generalized additive models (EGAM) (Youngman, 2019) and gradient
boosting for extreme quantile regression (GBEX) (Velthoen et al., 2021), which
both use conditional GPD modeling. For sequentially dependent data, we also
consider the extreme quantile autoregression (EXQAR) (Li and Wang, 2019) as,
although assuming linear quantile dependence, the model is designed for time
series. As benchmark, we consider an unconditional GPD model that ignores the
covariate dependence altogether, and a semi-conditional GPD model that uses
the covariate only for the intermediate quantile. We include results for the
generalized random forests for quantile regression (GRF) (Athey et al., 2019),
which does not use extrapolation for high quantiles. The training data
$\mathcal{D}=\left\\{(\mathbf{x}_{t},y_{t})\right\\}_{t=1}^{T}$, with
$T=\numprint{7000}$, are sequentially generated from the time series

$\begin{cases}Y_{t}=\sigma_{t}\left\lvert\varepsilon_{t}^{Y}\right\rvert,\quad
X_{t}=0.4\cdot
X_{t-1}+\left\lvert\varepsilon_{t}^{X}\right\rvert,\quad\varepsilon_{t}^{Y},\varepsilon_{t}^{X}\sim\mathcal{N}(0,1),\\\
\begin{split}\sigma_{t}^{2}=1&+0.1\cdot\\{2Y_{t-1}^{2}+Y_{t-2}^{2}+Y_{t-3}^{2}+Y_{t-4}^{2}+Y_{t-5}^{2}\\}+\\\
&+0.1\cdot\\{3X_{t-1}^{2}+2X_{t-2}^{2}+X_{t-3}^{2}+X_{t-4}^{2}+X_{t-5}^{2}\\}.\end{split}\end{cases}$
(13)

Figure S.5 in Supplementary Material S.4 shows part of the simulated data. To
have a fair comparison, all methods use the same covariate vectors
$\tilde{\mathbf{x}}_{t}=\\{(\mathbf{x}_{j},y_{j})\\}_{j=t-s}^{t-1}$ with
$s=10$. For the methods that use covariate dependent intermediate quantiles,
we use the same estimates $\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau_{0})$ with
$\tau_{0}=80\%$ from a recurrent quantile regression neural (QRN) network; for
a sensitivity analysis of the choice of $\tau_{0}$ see Supplementary Material
S.4. The best QRN structure and hyperparameters are chosen based on validation
quantile loss. For the methods that use covariate dependent GPD parameters, we
also use $\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau_{0})$ as additional covariate.
Although designed for univariate time series, we adapted EXQAR to accept
several covariate sequences.

For the EQRN model, $\numprint{2000}$ observations are kept for validation
tracking and thus only the remaining $\numprint{5000}$ are effectively used
for weight training. The best choices for EQRN hyperparameters are made based
on validation loss by performing a grid search over a set of possible values
and network architectures. All other models are fitted on the whole training
data set as they do not use validation loss tracking. The best set of
hyperparameters for GBEX (tree depths, learning rate and number of trees) are
chosen using cross-validation, and the ground truth for whether the shape is
constant is given to EGAM. For EXQAR, we use $\delta_{2n}=n^{-0.9}$ as
recommended by the authors, but set $\delta_{1n}=1-\tau_{0}$ thus increasing
the number of quantile pseudo-observations used for inference and allowing
better comparison with the other methods.

The predictions of all models are evaluated by their root mean squared error
(RMSE) compared to the true conditional quantiles on a newly generated test
data set that follows the same distribution as the training data in (13).

###  4.2 Results

For estimation of the conditional GPD parameters with our recurrent EQRN
model, we consider both LSTM and GRU structures with one to three recurrent
layers and hidden dimensions between $32$ and $256$. As the networks are not
too deep, $L_{2}$ penalty was chosen over dropout for regularization during
training, with possible penalty
$\lambda\in\\{0,10^{-6},10^{-5},10^{-4},10^{-3}\\}$. Both constant and
covariate dependent shape parameter outputs are considered. The model with
minimum validation loss is a single LSTM layer with hidden dimension $128$,
followed by the usual fully connected output layer, with constant shape and
$\lambda=10^{-4}$. As comparison, we also retain predictions for the best
unpenalized network with $\lambda=0$ and fixed shape, which has two LSTM
layers of hidden size $128$.

The left panel of Figure 5 shows the RMSE of best penalized and unpenalized
EQRN models compared to the competitors as a function of the quantile level
$\tau$. We observe that for the lowest level $\tau=\tau_{0}=80\%$, all
structured GPD models have the same performance since they use the same
intermediate quantiles. GRF and the unconditional model have already higher
errors since they are not able to capture the sequential dependence at the
intermediate level sufficiently. For growing quantile levels $\tau$, the
errors of the covariate dependent GPD models start to diverge. This is due to
the differences in modeling flexibility in terms of the GPD parameters of each
method. We observe a similar behaviour for EXQAR, as its linear quantile
dependence is not flexible enough. Our EQRN method based on recurrent neural
networks seems to be best at modeling sequential tail dependence.

Figure 5: Left: root mean squared error between predicted and true conditional
quantiles at different probability levels $\tau$ (log-scale), for the selected
EQRN models and competitors. Center-right: true versus predicted quantiles at
probability level $\tau=99.95\%$ for the best unpenalized (middle) and
penalized (right) EQRN models (dots), compared to the semi-conditional
estimates (crosses).

Figure 5 also shows the predicted quantiles
$\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau)$ on the test data compared to the true
$Q_{\tilde{\mathbf{x}}_{u}}(\tau)$ for a fixed $\tau=99.95\%$ for the best
penalized and unpenalized EQRN models. In general, both models seem to perform
well in predicting the high conditional quantiles. The weight penalty seems to
mainly affect the larger quantile predictions. Compared to the unpenalized
model, we observe that the reduction in the variance of the predictions comes
at the cost of a bias for larger quantile values. This bias-variance trade-off
is typical with penalization. The poor performance of the semi-conditional
estimates highlights the added value of covariate dependence in the GPD
parameters.

Additional results such as the impact of the intermediate level $\tau_{0}$ on
the EQRN accuracy as well as the corresponding quantile R squared coefficient
and bias-variance decomposition of the RMSEs presented in Figure 5 are
discussed in Supplementary Material S.4.

##  5 Application

###  5.1 Motivation

Flood risk is a major natural hazard in Europe, which causes huge economic
damage and endangers human lives. There is longstanding interest in
statistical methods of extreme value theory for hydrology (e.g., Katz et al.,
2002; Keef et al., 2009; Asadi et al., 2015; Engelke and Hitz, 2020), and
national agencies commonly use them to assess the long-term risk of flooding
in cities, at power plants and other key locations. Return levels with long
return periods can be estimated using the GEV distribution for annual maxima
or the GPD model for daily threshold exceedances. The output then guides
effective long-term flood management measures.

An example for an important location in Switzerland is the gauging station in
Bern on the Aare river, which is shown within its water catchment in Figure 1.
The Swiss Federal Office for the Environment (FOEN) monitors the Aare and we
use daily average discharges (in $\text{m}^{3}\text{s}^{-1}$) in Bern and
another upstream station together with recordings of daily precipitation (in
mm) at six locations in the Bern catchment; see Figure 1 for details. All time
series are available in the period from 1930–2014 and can be obtained from the
FOEN222https://www.hydrodaten.admin.ch/ (for discharges) and
MeteoSwiss333https://gate.meteoswiss.ch/idaweb (for precipitation). Figure 6
shows an excerpt for the two river stations and one precipitation gauge.

To illustrate possible drawbacks of a classical extreme value analysis, the
left panel of Figure 7 shows the annual maxima of river discharges at the
Bernese station on the Aare. The dashed line is the estimated 100-year return
level based on the GEV approximation using the training period from 1930–1958.
One can see that starting from the year 1999 there are several exceedances
over this return level, somewhat contradicting the fact that it should only be
exceeded in average once in 100 years. The solid line is the same return level
based on data from 1930–$y$, where $y\in\\{1959,\dots,2014\\}$ denotes end of
the training period. While the predictions are fairly stable until 1999, an
extreme value analysis performed after that year would yield much higher
values for the 100-year return level. Conversely, historical estimates before
such a break-point would severely underestimate the flood risk. In general,
distributional shifts can be due to climate change, changes in the river
system or other structural breaks in factors influencing discharge at this
location. For the Bernese station, the FOEN indeed reports a significant
break-point in extreme discharges in the nineties, but acknowledges that a
clear cause can not be identified444see flood report of the FOEN at
https://www.hydrodaten.admin.ch/en/2135.html. One factor may be a multi-
decadal variability of flood occurrence as described in Schmocker-Fackel and
Naef (2010).

Figure 6: Daily average discharge observations at Bern–Schönau (62) and at the
upstream station at Gsteig (42), and daily precipitation at the closest
meteorological station to Bern (BEP), over five years; see Figure 1 for
geographical locations of the gauging stations. Figure 7: Left: annual maxima
of daily average discharges (points) at Bern–Schönau (62) together with the
unconditional 100-year return level based on GEV fitted on the training data
of 1930–1958 (dashed line), and the evolution the same return level (solid
line) using data from 1930–$y$, where $y\in\\{1959,\dots,2014\\}$ denotes end
of the period. Middle: evolution of the validation loss (solid line) of the
selected EQRN network for the river discharge data as a function of the
training epoch; the dashed line shows the validation loss of the semi-
conditional model. Right: number of observations exceeding the EQRN
predictions on the test set (dotted line) compared to the expected number of
exceedances (solid line) for different probability levels (log-scale).

The aim of our methodology is complementary to classical extreme value
analysis and addresses this issue with static return levels. We apply our EQRN
model to estimate one-day-ahead extreme quantiles of the river discharge
conditionally on previous observations of discharge and precipitation in the
catchment. This allows the forecasting of flood risk even in non-stationary
systems such as a changing climate. The strength of our approach lies in the
ability to exploit information from multiple covariates and capture the
complex time dependence. Even in situations where the causes for structural
changes are unknown, our method implicitly accounts for them through their
effects on the covariates. The output of the model can help practitioners and
authorities to manage flooding more effectively and help minimizing their
disastrous impacts by early warning systems.

To illustrate our methodology and show its effectiveness in comparison to
classical forecasting approaches, we consider in more detail the flood in
August 2005 in Switzerland. At the Bernese gauging station, it was the largest
event since beginning of the recordings, and it caused severe economic damage
across large parts of the country and the loss of several lives.

###  5.2 Model Specification

The whole data set consists of 31,046 daily observations
$(\mathbf{x}_{t},y_{t})$ between 1930–2014. The response $y_{t}$ is the daily
average discharge at the Bernese gauging station on the Aare, and the
covariates $\mathbf{x}_{t}\in\mathbb{R}^{p}$, $p=7$, consist of discharge at
another upstream station and daily precipitation measurements from six
locations in the same catchment; see Figures 1 and 6. The discharges show
significant seasonality, both in trend and variance, with largest extremes
only appearing in the summer. We do not reduce artificially the non-
stationarity via classical approaches from times series analysis (e.g.,
Cleveland et al., 1990) as we believe the seasonality and other trends are
captured through the covariates.

We split the data into training and test sets. The first $T=10,349$
observations in the period between 1930–1958 are used to train the models,
whereof the first three-quarters serve the parameter estimation, and the
remaining quarter is a validation set to determine hyperparameters (sets
$\mathcal{T}$ and $\mathcal{V}$ in Algorithm 2, respectively). The test set
contains 20,697 observations from 1958–2014, which is used for neither fitting
nor selection of parameters, but only to evaluate the model performance on an
independent time period. We choose this rather small proportion of training
data to study the ability of the model to adapt to possible non-stationarity
over time without refitting. In particular the model weights are not updated
with any information from data after 1958, even for forecasts in the study of
the 2005 flood of interest. A large test set is also required to evaluate
extreme properties of the data distribution. As augmented covariates at time
$t$ for the recurrent neural network models we use the $s=10$ preceding days
and set $\tilde{\mathbf{x}}_{t}=\\{(\mathbf{x}_{j},y_{j})\\}_{j=t-s}^{t-1}$.
The augmented training set is then
$\tilde{\mathcal{D}}=\\{(\tilde{\mathbf{x}}_{t},y_{t})\\}_{t=s+1}^{T}$. Only
observations of the $p+1=8$ variables during the preceding ten days are used
to predict one-day-ahead for a new test time point.

As discussed in Section 3, we perform two steps for estimation of the
conditional tail model. First, we fit an intermediate quantile regression
model to estimate $Q_{\tilde{\mathbf{x}}_{t}}(\tau_{0})$. As in the
simulations with sequential dependence, we choose a recurrent QRN for this
purpose and set $\tau_{0}=0.8$. For the second step, we include the
intermediate quantile estimates $\hat{Q}_{\tilde{\mathbf{x}}_{j}}(\tau_{0})$
during the same time horizon $j=t-s,\dots,t-1$ as additional covariates, and,
slightly abusing notation, we denote $\tilde{\mathbf{x}}_{t}$ as the new
covariate vector; see Section 3.2. A recurrent EQRN is then fitted to the
exceedances for estimation of the conditional GPD parameters
$\sigma(\tilde{\mathbf{x}}_{t})$ and $\xi(\tilde{\mathbf{x}}_{t})$; see
Algorithm 2. Similarly to the study in Section 4.2, a grid-search is performed
on the training data to select the best hyperparameters and architectures for
both recurrent neural network models based on validation losses.

The final model chosen to regress the intermediate quantiles is a QRN with two
LSTM layers of dimension $256$, followed by the usual fully connected layer,
and $L_{2}$ weight penalty with parameter $\lambda=10^{-6}$. The chosen EQRN
model has two LSTM layers of dimension $16$, followed by a fully connected
layer, and $L_{2}$ weight penalty with parameter $\lambda=10^{-6}$.

###  5.3 Results

The middle panel of Figure 7 shows the validation loss (solid line) of the
selected EQRN model as a function of the training epoch. It can be seen that
already after a few epochs, the method has a lower loss than the simple, semi-
conditional model with constant GPD parameter $\sigma$ and $\xi$ (dashed
line). This shows that the GPD distribution varies with the predictor values
$\tilde{\mathbf{x}}_{t}$ and that a flexible model is beneficial.

The main output of the EQRN model are the extreme quantile estimates
$\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau)$ for a level $\tau$ and a time point
$u$ of interest, conditionally on the past covariates
$\tilde{\mathbf{x}}_{u}$. These one-day-ahead risk forecasts are shown as a
function of time on the test set in the top panel of Figure 2. We observe that
the model is able to extrapolate beyond the range of the data, since the event
shown in the plot is unprecedented and the predictions still anticipate the
first exceedance of the unconditional 100-year return level $Q^{100}$.

An unconditional $\tau$-quantile is defined as the value that is exceeded by a
proportion of $1-\tau$ of the data. An analogous property holds for
conditional quantiles in data with sequential dependence, which yields a
natural model assessment tool. On population level, if
$(\tilde{\mathbf{X}}_{t},Y_{t})_{t=1}^{T}$ is the random time series with
augmented covariate vectors, then the expected number of exceedances over the
true conditional $\tau$-quantiles $Q_{\tilde{\mathbf{X}}_{t}}(\tau)$ is

$\mathbb{E}\sum_{t=1}^{T}1\\{Y_{t}>Q_{\tilde{\mathbf{X}}_{t}}(\tau)\\}=(1-\tau)T.$

Consequently, plugging in the data $(\tilde{\mathbf{x}}_{t},y_{t})_{t=1}^{T}$
and estimates $\hat{Q}_{\tilde{\mathbf{x}}_{t}}(\tau)$ from a quantile
regression method, the equation should approximately hold if the model well
calibrated. Such a model assessment plot is shown in the right-hand panel of
Figure 7 for our EQRN fit as a function of the quantile level $\tau$. We
observe that the model is fairly well calibrated, with a slight bias towards
more exceedances than expected.

An additional output are the corresponding GPD parameters
$\sigma(\tilde{\mathbf{x}}_{u})$ and $\xi(\tilde{\mathbf{x}}_{u})$, which
together with the intermediate quantile
$\hat{Q}_{\tilde{\mathbf{x}}_{u}}(\tau_{0})$ specify the whole tail of the
distribution of $Y_{u}\mid\tilde{\mathbf{X}}_{u}=\tilde{\mathbf{x}}_{u}$
according to (5). For a given threshold level of interest $Q$ we can plot the
flood risk for the next day as the one-day-ahead forecast of the exceedance
probability over $Q$, that is, an estimate of the function
$u\mapsto\mathbb{P}(Y_{u}>Q\mid\tilde{\mathbf{X}}_{u}=\tilde{\mathbf{x}}_{u})$.
The bottom panel of Figure 2 shows the EQRN based estimate of this function on
the test set as a ratio to the unconditional $\mathbb{P}(Y_{u}>Q)$, where the
threshold $Q$ is chosen as the static 100-year return level $Q^{100}$ based on
the GEV distribution fitted on the training set; this threshold is relevant
since it is often used to determine the height of dams for flood management.
It results in a daily measure of how likely the exceedance on the next day is
compared to what was expected unconditionally. Times with large predicted
probability ratios are apparently times of imminent danger that can be used as
triggers for early warning systems or additional flood management measures.

As an example, one may issue a warning when the forecasted conditional
probability of exceeding $Q^{100}$ is, say, a hundred times larger than the
baseline unconditional probability $\mathbb{P}(Y_{u}>Q^{100})$. In the test
data there are four time clusters that typically last several days where
$Q^{100}$ is exceeded; see Figure 8 for one of these events. Applying this
early warning system, in all four of these cases a timely warning would have
been issued on the days preceding the first exceedance of the cluster. Such a
decision rule would lead to an average of only $1.3$ warnings for clusters of
exceedances per year on the test set, and is therefore not overly
conservative.

We consider the period of the 2005 flood in Bern in more detail. Figure 8
shows the two discharge time series and precipitation at the closest
meteorological station before and after the event. On the evening of August
21, 2005, the day preceding the first exceedance of the 100-year GEV return
level (horizontal dashed line), the prediction of our EQRN already indicated a
sudden increase in the probability of this exceedance; see bottom panel of
Figure 2. Equivalently, the blue point in Figure 8 shows the increased value
of the conditional 100-year quantile predicted by the model. The diamonds mark
the observations of the previous ten days that were used for this prediction.

In this case, the high precipitation values on August 21 and the preceding
days seem to have driven this prediction, possibly together with high values
of the rivers. It is interesting to note that a similar situation on August 2
has not resulted in a “flood warning” since the forecasted exceedance
probability and return level are not exceptionally high – as a matter of fact,
there was no exceedance on the next day. This means that the predictions are
driven by a complex combination of the risk factors $\tilde{\mathbf{X}}_{t}$
in space and time that are well captured by the recurrent architecture of the
EQRN.

Figure 8: Discharge at Bern–Schönau (62) and Gsteig (42) and precipitation at
BEP during the period of the 2005 flood, where diamonds indicate covariates
used for prediction of the 100-year conditional quantile (triangle) on August
22, 2005; other precipitation stations are not shown. The top panel also shows
unconditional 100-year return level (dashed line) fitted on the training data
and predictions on other days (solid line).

At the time, the FOEN used an adapted version of the hydrological model HBV
(Lindström et al., 1997) for forecasting river discharges based on several
inputs such as precipitation forecasts. However, the forecasts prior to this
event underestimated the flood risk and resulted in too late warnings.
Physical models for discharges and precipitation do not use explicit
extrapolation in the extreme tails and often have poor performance in the
largest predictions. In fact, during the 2005 flood, a main reason for the
late warnings was that forecasters did not trust the predicted precipitation
amounts during this extreme scenario. In the aftermath of this flood the FOEN
therefore published a detailed analysis of the internal forecasting procedures
(Bezzola and Hegg, 2007). The exact forecasts from that time are not
available, but the above discussion of the results shows that our statistical
methodology is a competitive alternative to physical models for the
forecasting of flood risk. We also note that our model uses much less
information as it relies only on observed discharge and precipitation and does
not require forecasts of atmospheric variables.

In Section S.5 of the Supplementary Material we compare and discuss the
forecasts of our EQRN method with those of some of the competitors. Figures
S.8-S.11 show that all methods seem to capture at least some of the temporal
structure based on the past covariates. Except for the GBEX method, the
forecasts of the competitors suffer however from a low sensitivity to changes
in the conditional tail or from a too erratic behavior as a function of time.
Overall, the recurrent structure of our EQRN method therefore seems to be the
best suited model for this kind of sequentially dependent data. Since in real-
world applications the true quantiles are unknown, direct computation of the
prediction error is difficult and model assessment plots as in the right-hand
panel of Figure 7 are crucial. This highlights the importance of simulation
studies to evaluate and compare quantitatively the accuracy of different
methods. In particular in situations with temporal dependence, our EQRN method
clearly outperforms the competitors (e.g., Figure 5). This is another
indicator to trust the EQRN forecasts in applications.

##  6 Conclusion

Our EQRN model combines extrapolation results from extreme value theory with
the prediction power of neural networks. It provides a flexible and versatile
method for extreme quantile regression that is capable of prediction beyond
the range of the data in the presence of a large number of covariates. The
main focus in this paper was the case of sequential dependence to develop a
tool for risk forecasting in time series that can be used for effective early
warning systems in flood management. Our model already performs well in
issuing sparse warnings for the days with increased risk of flooding as
illustrated in case study of the Aare catchment in Switzerland. A further
improvement could be attained by using additional covariates as input of the
model. This could include observations of variables that are typically used in
hydrological models such soil moisture, or forecasts of atmospheric variables
such as precipitation and temperature.

Many other applications seem pertinent. Even in the case of independent data,
our simulation study in Supplementary Material S.3 shows that neural networks
outperform tree-based methods such as ERF and GBEX if the quantile function is
more complex. Applications are financial risk assessment in insurance
companies or banks. For spatial data or images, convolutional neural networks
(LeCun et al., 2015) are known to perform extremely well in capturing
neighborhood structures. Our EQRN method can therefore be applied to quantify
the risk of climate extremes where the predictor space contains spatio-
temporal observations of meteorological variables (e.g., Boulaguiem et al.,
2022).

The price for the high flexibility of machine learning methods, which focus on
prediction accuracy, is a limited statistical interpretability. There is
active research on the construction of prediction intervals for black-box
methods, for instance through conformal inference (Lei et al., 2018; Romano et
al., 2019). How such techniques can be adapted to assess uncertainty for
extreme quantile regression is an interesting future research question.

#### Acknowledgements

The authors would like to thank Daniel Viviroli for his helpful comments.

#### Funding

Both authors were supported by Swiss National Science Foundation Eccellenza
Grant 186858.

## Supplementary material

#### Supplementary Results

The appended Supplementary Material contains additional information on
Algorithms 1 and 2, the simulation study on independent data, additional
results for the simulation study on dependent data, an analysis of the EQRN
sensitivity to the intermediate probability level and competitor approaches to
the application.

#### Source Code

The R implementation to reproduce the methodology, experiments and results
presented in this paper and more are available on
https://github.com/opasche/EQRN_Results. An easy-to-use open-source “EQRN” R
package is also under development, and will be made available on
https://github.com/opasche/EQRN.

## References

  * Andres et al. (2021) Andres, N., Steeb, N., Badoux, A., and Hegg, C., editors.  _Extremhochwasser an der Aare. Hauptbericht Projekt EXAR. Methodik und Resultate. [Extreme flooding of the Aare. Main Report on the EXAR Project. Methodology and results.]_ , volume 104 of _WSL Berichte_ , 2021\.
  * Asadi et al. (2015) Asadi, P., Davison, A. C., and Engelke, S.  Extremes on river networks.  _Ann. Appl. Stat._ , 9(4):2023–2050, 2015.
  * Athey et al. (2019) Athey, S., Tibshirani, J., and Wager, S.  Generalized random forests.  _Ann. Stat._ , 47(2):1148–1178, 2019.
  * Balkema and de Haan (1974) Balkema, A. A. and de Haan, L.  Residual Life Time at Great Age.  _Ann. Probab._ , 2(5):792 – 804, 1974.
  * Bücher and Zhou (2018) Bücher, A. and Zhou, C.  A horse racing between the block maxima method and the peak-over-threshold approach, 2018.  ArXiv:1807.00282.
  * Bezzola and Hegg (2007) Bezzola, G. R. and Hegg, C.  Ereignisanalyse Hochwasser 2005, Teil 1 — Prozesse, Schäden und erste Einordnung [Event analysis of the 2005 flood, Part 1 — Processes, damage and initial classification].  Technical report, Federal Office for the Environment FOEN, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, 2007.  Umwelt-Wissen Nr. 0707. 215 S.
  * Boulaguiem et al. (2022) Boulaguiem, Y., Zscheischler, J., Vignotto, E., van der Wiel, K., and Engelke, S.  Modeling and simulating spatial extremes by combining extreme value theory with generative adversarial networks.  _Environ. Data Sci._ , 1:e5, 2022.
  * Breiman (1996) Breiman, L.  Stacked Regressions.  _Mach. Learn._ , 24:49–64, 1996.
  * Cannon (2011) Cannon, A. J.  Quantile regression neural networks: Implementation in R and application to precipitation downscaling.  _Comput. Geosci._ , 37(9):1277–1284, 2011.
  * Cannon (2012) Cannon, A. J.  Neural networks for probabilistic environmental prediction: Conditional Density Estimation Network Creation and Evaluation (CaDENCE) in R.  _Comput. Geosci._ , 41:126–135, 2012.
  * Chavez-Demoulin and Davison (2005) Chavez-Demoulin, V. and Davison, A. C.  Generalized additive modelling of sample extremes.  _J. R. Stat. Soc. C_ , 54(1):207–222, 2005\.
  * Chernozhukov (2005) Chernozhukov, V.  Extremal quantile regression.  _Ann. Stat._ , 33(2):806 – 839, 2005.
  * Cho et al. (2014) Cho, K., van Merrienboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y.  Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation.  _CoRR_ , 2014.  ArXiv:1406.1078.
  * Cleveland et al. (1990) Cleveland, R. B., Cleveland, W. S., McRae, J. E., and Terpenning, I.  STL: A Seasonal-Trend Decomposition Procedure Based on Loess.  _J. Off. Stat._ , 6(1):3–73, 1990.
  * Cox and Reid (1987) Cox, D. R. and Reid, N.  Parameter Orthogonality and Approximate Conditional Inference (with discussion).  _J. R. Stat. Soc. B_ , 49:1–39, 1987.
  * Daouia et al. (2011) Daouia, A., Gardes, L., Girard, S., and Lekina, A.  Kernel estimators of extreme level curves.  _TEST_ , 20(2):311–333, 2011.
  * Duchi et al. (2011) Duchi, J., Hazan, E., and Singer, Y.  Adaptive Subgradient Methods for Online Learning and Stochastic Optimization.  _J. Mach. Learn. Res._ , 12:2121–2159, 2011.
  * Elman (1990) Elman, J. L.  Finding Structure in Time.  _Cogn. Sci._ , 14(2):179–211, 1990.
  * Engelke and Hitz (2020) Engelke, S. and Hitz, A.  Graphical models for extremes (with discussion).  _J. R. Stat. Soc. B_ , 82:871–932, 2020.
  * Fischer et al. (2021) Fischer, E., Sippel, S., and Knutti, R.  Increasing probability of record-shattering climate extremes.  _Nat. Clim. Chang._ , 11:689–695, 2021.
  * Fisher and Tippett (1928) Fisher, R. A. and Tippett, L. H. C.  Limiting forms of the frequency distribution of the largest or smallest member of a sample.  In _Math. Proc. Camb. Philos. Soc._ , volume 24, pages 180–190. Cambridge University Press, 1928.
  * Gardes and Stupfler (2019) Gardes, L. and Stupfler, G.  An integrated functional Weissman estimator for conditional extreme quantiles.  _REVSTAT_ , 17(1):109–144, 2019.
  * Gers et al. (2000) Gers, F. A., Schmidhuber, J., and Cummins, F.  Learning to Forget: Continual Prediction with LSTM.  _Neural Comput._ , 12(10):2451–2471, 2000.
  * Gers et al. (2003) Gers, F. A., Schraudolph, N. N., and Schmidhuber, J.  Learning precise timing with LSTM recurrent networks.  _J. Mach. Learn. Res._ , 3(Aug):115–143, 2003\.
  * Gnecco et al. (2022) Gnecco, N., Terefe, E. M., and Engelke, S.  Extremal Random Forests, 2022.  ArXiv:2201.12865.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., and Courville, A.  _Deep Learning_.  MIT Press, 2016.
  * Halton (1964) Halton, J. H.  Algorithm 247: Radical-inverse quasi-random point sequence.  _Commun. ACM_ , 7(12):701–702, dec 1964.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long Short-Term Memory.  _Neural Comput._ , 9(8):1735–1780, 1997.
  * Hornik (1991) Hornik, K.  Approximation capabilities of multilayer feedforward networks.  _Neural Netw._ , 4(2):251–257, 1991.
  * Jozefowicz et al. (2015) Jozefowicz, R., Zaremba, W., and Sutskever, I.  An Empirical Exploration of Recurrent Network Architectures.  In _Proc. 32nd Int. Conf. Mach. Learn._ , volume 37, pages 2342–2350, 2015.
  * Katz et al. (2002) Katz, R. W., Parlange, M. B., and Naveau, P.  Statistics of extremes in hydrology.  _Adv. Water Resour._ , 25:1287–1304, 2002.
  * Keef et al. (2009) Keef, C., Tawn, J., and Svensson, C.  Spatial risk assessment for extreme river flows.  _J. R. Stat. Soc. C_ , 58:601–618, 2009.
  * Kingma and Ba (2014) Kingma, D. P. and Ba, J.  Adam: A Method for Stochastic Optimization.  _3rd Int. Conf. Learn. Repres._ , 2014.
  * Kinsvater et al. (2016) Kinsvater, P., Fried, R., and Lilienthal, J.  Regional extreme value index estimation and a test of tail homogeneity.  _Environmetrics_ , 27(2):103–115, 2016.
  * Klambauer et al. (2017) Klambauer, G., Unterthiner, T., Mayr, A., and Hochreiter, S.  Self-Normalizing Neural Networks.  In _Proc. 31st Int. Conf. Neural Inf. Process. Syst._ , pages 972–981, 2017.
  * Koenker and Bassett (1978) Koenker, R. and Bassett, G.  Regression Quantiles.  _Econometrica_ , 46(1):33–50, 1978.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., and Hinton, G.  Deep learning.  _Nature_ , 521:436–444, 2015.
  * Lei et al. (2018) Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L.  Distribution-Free Predictive Inference for Regression.  _J. Am. Stat. Assoc._ , 113(523):1094–1111, 2018.
  * Li and Wang (2019) Li, D. and Wang, H. J.  Extreme Quantile Estimation for Autoregressive Models.  _J. Bus. Econ. Stat._ , 37(4):661–670, 2019\.
  * Lindström et al. (1997) Lindström, G., Johansson, B., Persson, M., Gardelin, M., and Bergström, S.  Development and test of the distributed HBV-96 hydrological model.  _J. Hydrol._ , 201:272–288, 1997.
  * Meinshausen (2006) Meinshausen, N.  Quantile Regression Forests.  _J. Mach. Learn. Res._ , 7(35):983–999, 2006\.
  * Pickands III (1975) Pickands III, J.  Statistical Inference Using Extreme Order Statistics.  _Ann. Stat._ , 3(1):119 – 131, 1975.
  * Romano et al. (2019) Romano, Y., Patterson, E., and Candes, E.  Conformalized Quantile Regression.  In _Adv. Neural Inf. Process. Syst._ , volume 32, 2019.
  * Schmocker-Fackel and Naef (2010) Schmocker-Fackel, P. and Naef, F.  Changes in flood frequencies in Switzerland since 1500.  _Hydrol. Earth Syst. Sci._ , 14(8):1581–1594, 2010.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.  Dropout: a Simple Way to Prevent Neural Networks from Overfitting.  _J. Mach. Learn. Res._ , 15:1929–1958, 2014.
  * Tieleman and Hinton (2012) Tieleman, T. and Hinton, G.  Lecture 6.5 – RMSProp.  Technical report, COURSERA: Neural Networks for Machine Learning, 2012\.
  * Velthoen et al. (2019) Velthoen, J., Cai, J.-J., Jongbloed, G., and Schmeits, M.  Improving precipitation forecasts using extreme quantile regression.  _Extremes_ , 22(4):599–622, 2019.
  * Velthoen et al. (2021) Velthoen, J., Dombry, C., Cai, J.-J., and Engelke, S.  Gradient boosting for extreme quantile regression, 2021.  ArXiv:2103.00808.
  * Wang et al. (2012) Wang, H. J., Li, D., and He, X.  Estimation of High Conditional Quantiles for Heavy-Tailed Distributions.  _J. Am. Stat. Assoc._ , 107(500):1453–1464, 2012.
  * Werbos (1988) Werbos, P. J.  Generalization of backpropagation with application to a recurrent gas market model.  _Neural Netw._ , 1(4):339–356, 1988.
  * Wolpert (1992) Wolpert, D. H.  Stacked generalization.  _Neural Netw._ , 5(2):241–259, 1992.
  * Youngman (2019) Youngman, B. D.  Generalized Additive Models for Exceedances of High Thresholds With an Application to Return Level Estimation for U.S. Wind Gusts.  _J. Am. Stat. Assoc._ , 114(528):1865–1879, 2019.
  * Zhang et al. (2019) Zhang, W., Quan, H., and Srinivasan, D.  An Improved Quantile Regression Neural Network for Probabilistic Load Forecasting.  _IEEE Trans. Smart Grid_ , 10(4):4425–4434, 2019.

SUPPLEMENTARY MATERIAL TO
“Neural Networks for Extreme Quantile Regression with an Application to
Forecasting of Flood Risk”

##  Appendix S.1 Additional LSTM Illustration

Figure S.1 shows a schematic representation of a multi-layer LSTM network.

##  Appendix S.2 Details on Algorithms 1 and 2

Algorithms 1 and 2 contain some abbreviated function calls. We give some
details here:

  1. RandomValidationSplit($\mathcal{I}$): For independent data, splits the index set $\mathcal{I}$ randomly into training set $\mathcal{T}$ and validation set $\mathcal{V}$ with pre-specified proportions.

  2. SequentialValidationSplit($\mathcal{I}$): For squential data, splits the index set $\mathcal{I}$ sequentially into training set $\mathcal{T}$ and validation set $\mathcal{V}$ with pre-specified proportions, such that all observations in $\mathcal{T}$ are before $\mathcal{V}$ in time.

  3. InitializeNetworkWeights($\Theta$) / InitializeRecurrentNetWeights($\Theta$): Initializes the weights of the GPD (recurrent) neural network randomly. The number of weights is determined by $\Theta$.

  4. GetMiniBatches($\mathcal{T}$): Splits the training set $\mathcal{T}$ into mini-batches for stochastic gradient descent.

  5. BackPropUpdate($\ell$, $\hat{\mathcal{W}}$, $\mathbf{x}_{\mathcal{B}}$, $\hat{Q}_{\mathbf{x}_{\mathcal{B}}}(\tau_{0})$, $\Theta$): Updates the parameter vector by a gradient step computed by backpropagation; may involve regularization methods such as $L_{2}$-penalty or dropout, specified in the hyperparameters $\Theta$.

  6. LossNotImproving($\hat{\mathcal{W}}$, $\mathbf{x}_{\mathcal{V}}$, $\hat{Q}_{\mathbf{x}_{\mathcal{V}}}(\tau_{0})$, $z_{\mathcal{V}}$): If validation loss is tracked, then also a stopping criterion (e.g., for early stopping) is specified and this function returns TRUE if this criterion is attained, indicating that the validation loss is not improving anymore.

$\mathbf{x}_{t-1}$$\mathbf{x}_{t-2}$$\mathbf{x}_{t-3}$$\mathbf{x}_{t-s}$LSTM1LSTM1LSTM1LSTM1$\mathbf{h}^{(1)}_{t-1}$$\mathbf{h}^{(1)}_{t-2}$$\mathbf{h}^{(1)}_{t-3}$$\mathbf{c}^{(1)}_{t-1}$$\mathbf{c}^{(1)}_{t-2}$$\mathbf{c}^{(1)}_{t-3}$LSTMLLSTMLLSTMLLSTML$\mathbf{h}^{(L)}_{t-1}$$\mathbf{h}^{(L)}_{t-2}$$\mathbf{h}^{(L)}_{t-3}$$\mathbf{c}^{(L)}_{t-1}$$\mathbf{c}^{(L)}_{t-2}$$\mathbf{c}^{(L)}_{t-3}$$\tilde{g}_{\mathcal{W}}(\tilde{\mathbf{x}}_{t})$$y_{t}$$\ell$
Figure S.1: Multi-layer LSTM network flowchart from input
$\tilde{\mathbf{x}}_{t}:=(\mathbf{x}_{t-s},\ldots,\mathbf{x}_{t-1})$ to output
$\tilde{g}_{\mathcal{W}}(\tilde{\mathbf{x}}_{t})$, with loss evaluation. The
LSTM cells represent the transformation in (9).

##  Appendix S.3 Simulation Study for Independent Observations

Three data generating models with independent observations are considered. The
training sample
$\mathcal{D}=\left\\{(\mathbf{x}_{i},y_{i})\right\\}_{i=1}^{n}$ is drawn from

$\begin{cases}\mathbf{X}\sim\mathcal{U}\left([-1,1]^{p}\right),\\\
Y\mid\mathbf{X}=\mathbf{x}\sim\sigma(\mathbf{x})\cdot
t_{\alpha(\mathbf{x})},\end{cases}$ (S.1)

with $p=10$,
$\alpha(\mathbf{x})=1/\xi(\mathbf{x}):=7\cdot\\{1+\exp(4x_{1}+1.2)\\}^{-1}+3$,
and three different models for $\sigma(\mathbf{x})$:

  1. Model 1: $\sigma(\mathbf{x}):=1+6\phi(x_{1},x_{2})$, where $\phi$ is the bivariate Gaussian density with correlation $0.9$,

  2. Model 2: $\sigma(\mathbf{x}):=4+3\cos(7\left\|(x_{1},x_{2})^{\top}\right\|_{2}+3)$,

  3. Model 3: $\sigma(\mathbf{x}):=4+3\cos(6\left\|\mathbf{x}\right\|_{2}+3.5)$.

The expression of $\alpha(\mathbf{x})$ and of $\sigma(\mathbf{x})$ in Model 1
are the same as in the simulation study in Velthoen et al. (2021). The choices
for Models 2 and 3 are desgined to study more complex covariate dependencies.
The constants are chosen to have positive scale values and enough variation
for the inference task to be interesting.

As the focus of the experiments is on the extremal part of the model and since
all extreme value competitors use the same intermediate quantile estimates, we
use the true $Q_{\mathbf{x}}(\tau_{0})$ as intermediate quantiles. Two main
types of network architectures are considered for EQRN. The first type is an
MLPs with $\tanh$ activation functions, narrow structures with between one and
four hidden layers and optional $L_{2}$ weight penalty as training
regularization. The second type is self-normalizing networks (Klambauer et
al., 2017) using SELU activation functions, deeper architectures with between
three and eight hidden layers and optional alpha dropout as regularization.
This type of network is designed to maintain unit variance and zero mean
across layers in deeper networks trained for regression tasks, in order to
avoid vanishing and exploding gradient issues. Results suggest that the
additional flexibility of the second type was not necessary for the three
tasks at hand as they never yielded better validation score than the best
$\tanh$ models. Table 1 summarizes the hyperparameters of the chosen EQRN
networks.

To evaluate the accuracy of the best models over the full feature space
$\mathcal{X}=[-1,1]^{p}$, we use the integrated squared error (ISE) between
the prediction $\hat{Q}_{\mathbf{x}}(\tau)$ and the true quantile
$Q_{\mathbf{x}}(\tau)$,

$\int_{\mathcal{X}}\left(\hat{Q}_{\mathbf{x}}(\tau)-Q_{\mathbf{x}}(\tau)\right)^{2}d\mathbf{x}.$
(S.2)

We generate test features using a Halton sequence (Halton, 1964) and compute
the MSE between the corresponding predicted and true response quantiles, to
estimate the $p$-dimensional integral.

Figure S.2 shows the accuracy of EQRN and the competitor methods for an
increasingly large $\tau$ for the three data models. For every model, EQRN
outperforms all competitors with a difference in accuracy increasing with
$\tau$. EGAM seems to suffer from the large dimension of the feature space at
large $\tau$, both when the actual quantile depends on only two or all of the
ten features. The difference in accuracy of EQRN and GBEX compared to the
unconditional and semi-conditional models is particularly significant for
Models 1 and 3, and EQRN generally outperforms GBEX especially at high
probability levels.

Table 1: Hyperparameters of the EQRN networks with best validation loss for the three independent data models. | Hidden activation function | Hidden layer dimensions |  $L_{2}$ penalty
---|---|---|---
Model 1 | $\tanh$ | $(128,128,128)$ | $10^{-5}$
Model 2 | $\tanh$ | $(20,10,10)$ | $10^{-5}$
Model 3 | $\tanh$ | $(10,10,10)$ | $0$
Figure S.2: Root integrated squared error between predicted and true
conditional quantiles at different probability levels $\tau$ (log-scale) for
the selected EQRN model and competitors, for data Models 1–3 (left to right).
The cropped-out RMISE for EGAM at level 0.9999 are around 43, 115 and 150,
respectively.

We also define the quantile R squared of $\hat{Q}_{\mathbf{x}}(\tau)$ over the
sample $\mathcal{D}$ as

$R^{2}_{\tau}:=1-\dfrac{\sum_{i=1}^{n}\left(Q_{\mathbf{x}_{i}}(\tau)-\hat{Q}_{\mathbf{x}_{i}}(\tau)\right)^{2}}{\sum_{i=1}^{n}\left(Q_{\mathbf{x}_{i}}(\tau)-\overline{Q_{\mathcal{D}}(\tau)}\right)^{2}},\quad\text{
with
}\quad\overline{Q_{\mathcal{D}}(\tau)}:=\dfrac{1}{n}\sum_{i=1}^{n}Q_{\mathbf{x}_{i}}(\tau).$
(S.3)

The definition is similar to the classical R squared coefficient of
determination in regression, but the true conditional quantile values are used
as targets instead of the response observations. The $R^{2}_{\tau}$ is
essentially the reversed MSE normalized by the variance of the true
conditional quantile. A value close to unity indicates a very low MSE compared
to the quantile variance and negative values indicate a MSE larger than the
quantile variance.

Figure S.3 shows the quantile R squared, the biases and the residual standard
deviations of the same respective quantile predictions compared to the truth.
The R squared lead to the same conclusions as the RISE. Regarding the bias-
variance decomposition of the RISE, it seems that the variance term is
dominating the square bias. Although EQRN is here not the least biased model
for large $\tau$ values, it is it’s dominating performance in terms of
residual variance that leads to it having the lowest RISE values for the three
models.

Figure S.3: Quantile R squared, bias and residual standard deviation of the
predicted quantiles compared to the truth at different probability levels
(log-scale) for the selected EQRN model and competitors, for data Models 1–3
(top to bottom).

Figure S.4 shows the predicted $\hat{Q}_{\mathbf{x}}(0.9995)$ for EQRN and the
competitors, as a function of the two significant covariates $(x_{1},x_{2})$
for Model 1. At that extreme level, EQRN still seems to capture the true
conditional quantile function quite well, although the predictions show some
residual noise. GBEX shows an elliptical stepwise approximation behaviour. and
seems to underestimate the smaller quantiles and overestimate the largest
quantiles. EGAM and GRF fail drastically at recovering the conditional
quantile function. The semi-conditional estimates are a translation of the
intermediate quantiles, which in particular fail to capture the varying shape
along $X_{1}$.

Figure S.4: Conditional quantile predictions of iEQRN and competitor models at
probability level $\tau=0.9995$, shown as a function of $X_{1}$ and $X_{2}$
for Model 1.

##  Appendix S.4 Simulation Study for Sequentially Dependent Data

The main results of the simulation study on sequentially dependent data is
presented in the main paper. This section discusses additional results. Figure
S.5 shows part of the sequential data simulated from the generating process
described in the main paper.

Figure S.5: First $\numprint{1000}$ observations of the sequential data
simulated from (13).

Figure S.6 shows the quantile R squared (S.3), the biases and the residual
standard deviations of the quantile predictions compared to the truth, for the
two selected EQRN models and competitors. The R squared evolution again shows
EQRN is the model that best captures the covariate sequential dependence in
the tail, as it outperforms all competitors with a difference in accuracy
increasing with $\tau$. In terms of bias, the penalized EQRN here scores
similar values as EXQAR, and both EQRN versions outperform all other methods.
The unpenalized EQRN has the lowest residual variance, closely followed by
both GBEX and the penalized EQRN, although GBEX has a bad accuracy overall,
due to its large bias.

Figure S.6: Quantile R squared, bias and residual standard deviation of the
predicted quantiles compared to the truth at different probability levels
(log-scale) for the selected EQRN model and competitors, for the sequential
data model.

Figure S.7 shows the impact of the intermediate level $\tau_{0}$ on the
accuracy of the best EQRN model. The value $\tau_{0}=0.8$ used in the rest of
the analysis leads to the best RMSE. This relatively low value shows that
$\tau_{0}$ can be chosen much lower than for the classical unconditional GPD
model. An intuitive explanation for this fact is the following. In a situation
without covariates, the choice of $\tau_{0}$ is a trade-off between
approximation bias (which favours larger thresholds) and variance (which
favours lower thresholds); see Figure 3 in the main document. For covariate-
dependent data, the distribution of the exceedances varies and more data is
needed to accurately capture this function of the covariates. The variance
becomes more important than the approximation bias and therefore lower
thresholds are preferable. Moreover, a possible bias can be “corrected” by the
high flexibility of the model.

The final accuracy seems in fact to not be too sensitive to the choice for
$\tau_{0}$ compared to the network’s grid-searched hyperparameters, discussed
in the main analysis, as the differences in RMSE for $\tau_{0}$ values close
to the optimum are relatively small.

Figure S.7: Boxplot of the absolute residuals of the quantile predictions from
the selected EQRN model (blue) and their RMSE (red diamond) at probability
levels $0.995$ (left), $0.999$ (middle) and $0.9995$ (right) for different
choices of intermediate probability level $\tau_{0}$, for the sequential data
model.

##  Appendix S.5 Application: Competitor Results

The main results from our application to forecasting of flood risk in
Switzerland using our proposed EQRN methodology are presented in the main
paper. This section discusses and compares additional results using the
competitor methods, adapted to provide the same type of forecast as the EQRN
approach, with a focus on the 2005 flood event (see Figure 2 in the main
paper).

We first observe that the predictions have roughly the same behaviour, which
shows that all methods capture at least some of the temporal structure based
on the past covariates. The semi-conditional method (Figure S.8) is clearly
not flexible enough since the predictions only show a weak sensitivity to the
changes in covariates. The reason is that, while the intermediate quantile is
covariate-dependent, the GPD parameters are constant over time. We conclude
that a covariate-dependent model for the tail is required in this application.
Figure S.9 shows predictions from the EXQAR model (Li and Wang, 2019). This
model is more sensitive to changes in the covariates, but the regression
function looks fairly erratic, with sudden spikes at some time points. Those
spikes are here due to unusually small shape estimates in combination with a
large scale estimate. This might be caused by an instability in the estimation
of the local moments used in the model. EGAM (Figure S.10) seems to severely
overestimate the river flow at the main event. The best competing model seems
to be the GBEX (Figure S.11), as it yields a smooth prediction curve with
pronounced spike at the main event. Comparing this with our EQRN method
(Figure 2 in the main paper), it reacts slightly later and its risk forecast
(probability ratio) on the day before the first exceedance is significantly
lower than for EQRN

Figure S.8: Top: Daily average discharge (points) at the Bern–Schönau station
(62) and one-day-ahead semi-conditional forecasts of conditional 100-year
quantiles (solid line) during the 2005 flood. Horizontal lines show
unconditional $Q^{100}$ based on GEV (dashed) and GPD (dotted). Bottom: one-
day-ahead forecast of conditional probability of exceeding the GEV estimated
$Q^{100}$ as a ratio to the unconditional probability, using the semi-
conditional parameter forecast. The vertical line indicates August 22, the day
of the first exceedance Figure S.9: Top: Daily average discharge (points) at
the Bern–Schönau station (62) and one-day-ahead EXQAR forecasts of conditional
100-year quantiles (solid line) during the 2005 flood. Horizontal lines show
unconditional $Q^{100}$ based on GEV (dashed) and GPD (dotted). Bottom: one-
day-ahead forecast of conditional probability of exceeding the GEV estimated
$Q^{100}$ as a ratio to the unconditional probability, using the EXQAR
forecast. The vertical line indicates August 22, the day of the first
exceedance Figure S.10: Top: Daily average discharge (points) at the
Bern–Schönau station (62) and one-day-ahead EGAM forecasts of conditional
100-year quantiles (solid line) during the 2005 flood. Horizontal lines show
unconditional $Q^{100}$ based on GEV (dashed) and GPD (dotted). Bottom: one-
day-ahead forecast of conditional probability of exceeding the GEV estimated
$Q^{100}$ as a ratio to the unconditional probability, using the EGAM
parameter forecast. The vertical line indicates August 22, the day of the
first exceedance. Figure S.11: Top: Daily average discharge (points) at the
Bern–Schönau station (62) and one-day-ahead GBEX forecasts of conditional
100-year quantiles (solid line) during the 2005 flood. Horizontal lines show
unconditional $Q^{100}$ based on GEV (dashed) and GPD (dotted). Bottom: one-
day-ahead forecast of conditional probability of exceeding the GEV estimated
$Q^{100}$ as a ratio to the unconditional probability, using the GBEX
parameter forecast. The vertical line indicates August 22, the day of the
first exceedance.

This qualitative way of checking the model is important since in real-world
applications, the true extreme quantiles are unknown. Therefore only some
quantitative model checks can be performed (e.g., right-hand panel of Figure 7
in the main paper). This highlights the importance of the simulation studies,
which allow us to evaluate in several settings which methods are more
accurate. In particular, in situations with temporal dependence, our EQRN
method outperforms the competitors (e.g., Figure 5 in the main paper).