# Stochastic generator for rainfall with a Hawkes process marked by an extended generalized Pareto and a vine copula

Antoine Chapon, Taha B.M.J. Ouarda, Nathalie Bertrand

► To cite this version:

Antoine Chapon, Taha B.M.J. Ouarda, Nathalie Bertrand. Stochastic generator for rainfall with a Hawkes process marked by an extended generalized Pareto and a vine copula. Environmental Modelling and Software, 2025, 191, pp.106490. 10.1016/j.envsoft.2025.106490. hal-05074179

HAL Id: hal-05074179

https://asnr.hal.science/hal-05074179v1

Submitted on 19 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stochastic generator for rainfall with a Hawkes process marked by an extended generalized Pareto and a vine copula

Antoine Chapon<sup>a,b,∗</sup>, Taha B.M.J. Ouarda<sup>a</sup>, Nathalie Bertrand<sup>b</sup>

## ARTICLE INFO

Dataset link: https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=overview

**Keywords:**

Rainfall

Stochastic generator

Extreme values

Extended generalized Pareto

Vine copula

## ABSTRACT

A stochastic generator for rainfall is built from a Hawkes process, which is modeling the occurrence and serial correlation of non-zero rainfall values. Hawkes processes are suited to model intermittent signals, which is the case of rainfall at a fine enough observation frequency. This Hawkes process has a two-scale intensity function accounting for two orders of clustering in rainfall time series. The rainfall amount of each non-zero value is modeled by an extended generalized Pareto (EGP) distribution with the whole range of rainfall as support, from low to extreme values. New parametric EGP forms adapted to high frequency rainfall time series are defined. The Hawkes process only models the serial correlation of occurrences but not that of the amounts. A conditional version of the EGP is hence developed by adding a copula, modeling the temporal dependence of rainfall amounts. A subsettable canonical vine copula models this dependency for multiple time lags, while accounting for the intermittency of non-zero rainfall values. An application to a 40 yr time series of hourly rainfall in France is presented. Simulations from the model reproduce adequately the marginal distribution of rainfall, the temporal clustering of events, and the autocorrelation. The simulations are also able to reproduce the intensity-duration-frequency relation of the IDF extreme value model, showing that this stochastic generator is suitable for risk assessment of duration-dependent extremes.

## 1. Introduction

Rainfall modeling is important for flood risk assessment. Rainfall extremes are duration-dependent because the risk depends on both the amount and the duration of an event. Extreme value distributions and stochastic weather generators (WGs) are two types of statistical models suitable to model extreme rainfall. In recent years, WGs applied to rainfall have integrated extreme value distributions to properly account for its upper tail (Beneyto et al., 2023), which brings the two types of models closer. The opposite is also true, as expanding the classic extreme value model of the generalized Pareto (GP) can end up transforming it into a WG.

The common extreme value framework for a duration-dependent variable is the intensity-duration-frequency (IDF) model (Langousis and Veneziano, 2007; Ouarda et al., 2018), which considers the amount of rainfall (called intensity in this model) as a stochastic process conditional on both frequency and duration. The IDF model is an extension of the generalized extreme value (GEV) distribution, which considers the maxima over long blocks of time (i.e., block maxima) as a stochastic process. The IDF links these block maxima to both frequency, as in the GEV, and duration. Therefore, the duration of events is not treated

as a stochastic variable in the IDF model. For non duration-dependent variables, the alternative to the GEV model is the GP distribution, which models independent extreme values above a high threshold. The advantage of the GP over the GEV is that the model can be inferred on more observations because the extremes are not defined as the single maxima in each long temporal window, but rather exceedances of a threshold, which improves the accuracy of the inferences (Coles, 2001). However the two drawbacks of the GP is the difficulty to select an appropriate threshold, and the requirement that exceedances are independent, which often involves keeping only the maxima of each cluster of exceedances (Lang et al., 1999; Pan et al., 2023).

The threshold selection for the GP can be simplified by adding flexibility to the distribution with an extended generalized Pareto (EGP) (Papastathopoulos and Tawn, 2013). Selecting an appropriate threshold for an EGP becomes less critical than for the simpler GP because the added flexibility can adapt to different thresholds. The threshold can now be set at a sub-asymptotic level where both extremes and non-extreme values exceed it. For a variable with positive support, such as rainfall, the threshold can even be lowered down to zero, if the EGP is flexible enough, while retaining a valid extreme value

<footnote><sup>a</sup> Institut national de la recherche scientifique, Centre Eau Terre Environnement, 490 rue de la Couronne, Québec, G1K 9A9, Canada</footnote>

<footnote><sup>b</sup> Autorité de sûreté nucléaire et de radioprotection, 31 avenue de la Division Leclerc, Fontenay-aux-Roses, 92260, France</footnote>

<footnote>* Corresponding author at: Institut national de la recherche scientifique, Centre Eau Terre Environnement, 490 rue de la Couronne, Québec, G1K 9A9, Canada. E-mail address: antoine.chapon@inrs.ca (A. Chapon).</footnote>
distribution in the upper tail (Naveau et al., 2016). Therefore, the EGP corrects one of the two drawbacks of the GP.

The second drawback of the GP requiring independent observations can be corrected by modeling the serial correlation of threshold exceedances. Li et al. (2021) developed such model by considering the threshold exceedances as points in a Hawkes process. This class of point process models the autocorrelation of points by conditioning the intensity on previous points (Hawkes, 2018). In other words, the Hawkes process models the temporal clustering of threshold exceedances, so every exceedance can now be kept and not only a subset of independent ones. A mark can be associated with each point, which is a value distributed by an EGP in the model of Li et al. (2021).

Therefore, for a variable with positive support such as rainfall, a Hawkes process marked by an EGP corrects the two drawbacks of the GP and allows the model to be fitted on every observation, while remaining an extreme value model. Since this marked Hawkes process models the whole time series of observations, it is also able to simulate new time series, and as such is a stochastic WG.

Risk assessment from extreme value distributions is commonly expressed as the quantile of a given probability, or the probability of a given quantile, with the corresponding uncertainties. However the risk is not caused directly by the rainfall. It is rather caused by its transformation into runoff, and ultimately by the water level reached at a vulnerable location. As such, a sensible way to express this risk is with respect to a design variable, which can be the maximum water level reached at a vulnerable location without failure. A WG can simulate many time series, which can then be the input of a rainfall-runoff model, and a probabilistic risk assessment can be obtained by the ratio of simulations leading to a failure. Simulations from the fitted model can extrapolate the upper extremes beyond the range of observed values through the EGP.

This work presents a WG built from a Hawkes process marked by an EGP with the whole range of rainfall as support. Novelties of this model are the definition of new parametric EGPs for high frequency rainfall, a two-scale intensity Hawkes process intensity function to account for two orders of temporal clustering, and the addition of a subsettable canonical vine copula to the EGP to model the serial correlation of marks for several time lags. This WG is designed as a threshold-based equivalent to the block maxima-based IDF model, as the GP is to the GEV. For this model to be a potential alternative to the IDF, it must be able to reproduce the relation between amount, duration and frequency of extreme events modeled by the IDF. To evaluate if this is the case, the IDF model is used as a diagnostic by applying it to the observations and a simulated time series of the same length.

Note that it is common to refer to the amount of rainfall for a given time interval as the rainfall intensity, for example in the IDF model, but to avoid any confusion the term “intensity” will be used to refer to point processes. Note also that the temporal point process discussed here is different from the point process representation of an extreme value distribution (Smith, 2003), since the amount of rainfall of a given timestep is not a dimension of the point process, but is instead modeled by a mark.

Several WGs for rainfall or other hydrological variables have been developed in recent years. Evin et al. (2018) built a multisite WG for daily precipitation with an EGP for the marginal distribution and a Markov chain for the occurrence of non-zero values. Papalexiou et al. (2023) also developed a multisite parametric WG for daily precipitation, with a combination of multivariate Gaussian process for the spatial dependence and distributions on [0, 1) having probability mass at 0 to model the intermittency of precipitation. Ahn (2020) developed a WG for precipitation and temperature by separating short term and long term variability. In their work the short term variability is modeled by a parametric distribution with an EGP and a spatial copula, while the long term variability is modeled by wavelet decomposition. Lee et al. (2020) used a deep learning model with recurrent layers in time to simulate monthly streamflow. Some of these models estimate

the parameters separately for each month to account for the annual variability (Evin et al., 2018; Papalexiou et al., 2023).

The motivation behind the present work is to develop a new type of WG for hourly rainfall by combining EGP, Hawkes process and vine copula in a novel framework. Combining these three statistical methods permits to build a WG with an extreme value distribution that also accounts for the duration and intermittency of rainfall events. The Hawkes process modeling the intermittency could be removed from the model to build a WG for a non-intermittent variable, for example temperature. Non-stationarity of the rainfall distribution is left out of the scope of this work. The scope is also restricted to the WG and does not cover the rainfall-runoff model needed to obtain a flood risk assessment by simulation. These restrictions of the scope result that the proposed model is not operational, and would require additional developments to this end.

The model presented in this work is not specific to a particular location and should a priori be suitable for any location where the majority of precipitations are rainfall. The model is demonstrated with a 40-year time series of hourly rainfall in France, but longer or shorter time series could be used, granted that it is long enough for a proper estimation of the shape parameter of the EGP.

The presentation of the method is split into four parts. Section 2 presents new parametric forms for the EGP which are adapted to high frequency rainfall. Section 3 presents a Hawkes process with a novel intensity function for rainfall time series. Section 4 presents the addition of a vine copula to the EGP-marked Hawkes process to account for the serial correlation of rainfall amount. Section 5 presents the inference method via likelihood and how to simulate from the model for risk assessment. An application to an hourly time series is presented in Section 7. Section 8 concludes.

## 2. Distribution of hourly rainfall

### 2.1. Extended generalized pareto

Extreme values above a sufficiently high threshold follow the GP distribution. The cdf of the GP is

$$G(m|\sigma, \xi) = \begin{cases} 1 - \left(1 + \frac{\xi m}{\sigma}\right)^{-1/\xi} & \text{for } \xi \neq 0, \\ 1 - \exp\left(-\frac{m}{\sigma}\right) & \text{for } \xi = 0, \end{cases} \quad (1)$$

where $\sigma > 0$ is the scale parameter and $\xi$ is the shape parameter.

The EGP is a class of distributions where the probability integral transform of a distribution modifies and adds flexibility to the GP (Papathathopoulos and Tawn, 2013). This added flexibility permits to lower the threshold below the domain of the GP, to have a distribution for both extreme and non-extreme values.

The cdf of the EGP is

$$K(m|\sigma, \xi, \kappa) = F_V\{G(m|\sigma, \xi)|\kappa\}, \quad (2)$$

where $F_V$ is the distribution used as probability integral transform and $\kappa$ are the parameters of this transform. The corresponding pdf is

$$k(m) = f_V\{G(m)\} g(m). \quad (3)$$

For it to be a valid distribution, $F_V$ must have [0, 1] as support and meet the condition

$$\lim_{v \to 1} f_V(v) = a, \quad (4)$$

where $a > 0$ is finite (Gamet and Jalbert, 2022), so that the tail index $\xi$ of the nested GP is unaffected by the transformation.

The most popular EGP uses the power function $F_V(v) = v^\kappa$ with $\kappa > 0$ as transform, thereafter named the EGP-power. Naveau et al. (2016) applied the EGP-power to model the whole range of hourly rainfall. The EGP-power is suited for rainfall, as it increases the density at the lower tail compared to the GP when $0 < \kappa < 1$. However, the one
parameter power transform can lack flexibility for rainfall time series with high observation frequency. Gamet and Jalbert (2022) proposed to build EGPs with truncated distributions as $F_V$ transforms, and defined one suitable for rainfall.

Gamet and Jalbert (2022) also added the condition

$$ \lim_{v \to 1} f'_{V}(v) = 0, \qquad (5) $$

on the derivative of the transform's density, which ensures that the GP is asymptotically unmodified by $F_V$.

## 2.2. EGPs defined by truncated distributions

A truncated distribution $T(x)$ with support $[l, u]$, is defined from a distribution $F(x)$ by

$$ T(x) = \frac{F[\max\{\min(x, u), l\}] - F(l)}{F(u) - F(l)}, \qquad (6) $$

where $l$ and $u$ are the lower and upper boundaries of the truncation, respectively (Nadarajah and Kotz, 2006). The corresponding density is given by

$$ t(x) = \begin{cases} \frac{f(x)}{F(u) - F(l)}, & \text{if } l \le x \le u, \\ 0, & \text{otherwise.} \end{cases} \qquad (7) $$

The EGP-beta of Gamet and Jalbert (2022) is defined with a truncated beta distribution. The pdf of its probability integral transform is given by

$$ f_V(v|\kappa, l) = \frac{f_B\{(u-l)v+l\}(u-l)}{B(u)-B(l)}, \qquad (8) $$

where $f_B$ and $B$ are the pdf and cdf of the beta distribution, respectively. The beta distribution has two parameters $\kappa = (\kappa_1, \kappa_2)$. Gamet and Jalbert (2022) parametrized the EGP-beta_1 with $\kappa_1 = \kappa_2$ and fixed truncation boundaries $l = 1/32$ and $u = 1/2$, resulting in a single parameter for the probability integral transform. The lower truncation value of $l = 1/32$ was empirically defined by Gamet and Jalbert (2022). When $\kappa_1 = \kappa_2$ the derivative $f'_B(x)$ is always 0 at $x = 1/2$, so the condition of Eq. (5) is met when $u = 1/2$.

A more flexible EGP-beta_3 with three parameters can be defined by having $l$ as a parameter, with $0 \le l < u$, and allowing $\kappa_1 \ne \kappa_2$. The upper truncation is then obtained by

$$ u = \frac{\kappa_1 \kappa_2 - 1}{\kappa_2 - 2}, \qquad (9) $$

with $f'_B(u) = 0$.

Following the same principle as the EGP-beta_3, an EGP is defined with a truncated generalized beta of the first kind (McDonald, 1984). The pdf of the generalized beta of the first kind is given by

$$ f_{gb}(x|a, b, c) = c B(a, b)^{-1} x^{a-1} (1-x^c)^{b-1}, \qquad (10) $$

where $a > 0$, $b > 0$, and $c > 0$. The corresponding cdf is given by

$$ F_{gb} = B(x^c|a, b). \qquad (11) $$

The probability integral transform of the EGP-genbeta is defined with $(a, b, c) = \kappa$. Its pdf and cdf are given by

$$ f_V(v|\kappa, l) = \frac{f_{gb}\{(u-l)v+l\}(u-l)}{F_{gb}(u)-F_{gb}(l)} \qquad (12) $$

and

$$ F_V(v|\kappa, l) = \frac{F_{gb}(v) - F_{gb}(l)}{F_{gb}(u) - F_{gb}(l)}, \qquad (13) $$

respectively, where the upper truncation bound $u$ is given by

$$ u = \left( \frac{\kappa_1 \kappa_3 - 1}{\kappa_1 \kappa_3 + \kappa_2 \kappa_3 - \kappa_3 - 1} \right)^{1/\kappa_3}, \qquad (14) $$

with $f'_{gb}(u) = 0$. The EGP-genbeta has four parameters $(\kappa_1, \kappa_2, \kappa_3, l)$.

## 3. Rainfall occurrence and intermittency

In a rainfall time series, the non-zero timesteps can be seen as occurrences, and thus can be modeled by a temporal point process. Since rainfall is recorded at a discrete time interval, this point process is discrete as well. Most rainfall events will be observed as several consecutive non-zero values, and similarly most dry periods will span several timesteps. Therefore, from the perspective of the discrete point process, both occurrences and non-occurrences must be modeled as dependent and clustered in time.

A point process is characterized by its intensity function, which gives the expected number of points in an interval. In a discrete setting, this interval is the time series timestep and only one point can occur at each timestep, so the intensity gives the probability of occurrence. The simplest point process is the Poisson process which has a constant intensity, resulting in independent occurrences which would be inappropriate for rainfall. A point process suitable for rainfall must have an intensity function able to account for the dependence of occurrences.

This temporal point process only models the occurrence or non-occurrence of rainfall in a binary way. A mark, here modeled by an EGP, accounts for the amount of rainfall of each occurrence. In a general way, the mark of a point process carries extra dimensions that are not directly included in the intensity function. In this work, the intensity of the point process has only a temporal dimension, and the rainfall amount is a second dimension whose stochastic behavior is not a dimension of the point process intensity. However the value of past marks can have an impact on the intensity of the point process. As an example, high amounts of rainfall could increase the short-term probability of occurrence.

### 3.1. Self-exciting Hawkes process for time series

Hawkes processes are a family of point processes able to account for temporal dependence with a variable intensity conditional on past occurrences (Hawkes, 2018). The varying intensity $\lambda(t)$ of a discrete Hawkes process is given by

$$ \lambda(t) = \mu + \sum_{i \mid t < t} f(t - i), \qquad (15) $$

where $\mu$ is a base level intensity and the function $f(\cdot)$ controls the influence of past occurrences at times $i \in [1, t-1]$ on the intensity at time $t$. A common function for $f(\cdot)$ is the exponential kernel, which exponentially decreases the influence of past occurrences that are further in time from $t$, asymptotically reaching 0. The intensity function is then given by

$$ \lambda(t) = \mu + \sum_{i \mid t < t} \gamma_1 \exp \{\gamma_2 - \gamma_2(t - i)\}, \qquad (16) $$

where $\gamma_1 > 0$ controls how much past occurrences influence the intensity and $\gamma_2 > 0$ controls the decay rate of this influence. The exponential kernel results in a self-exciting Hawkes process, in which occurrences increase the short-term probability of future occurrences. In discrete time, this self-exciting point process models clusters of consecutive occurrences.

New occurrences are possible without the influence of past occurrences (or an almost null influence in the case of the exponential kernel) thanks to the base intensity $\mu$. This component of the intensity which does not depend on the past is typically modeled as a constant but could also be itself varying, for example as a function of time (Hawkes, 2018).

The rainfall amount of each occurrence is a mark modeled by an EGP. The influence of past marks on the intensity is added via a two-parameter impact function

$$ M(i) = \delta_1 m_i^{\delta_2}, \qquad (17) $$

where $m_i$ is the mark value at time $i$, $\delta_1 > 0$ and $\delta_2 > 0$.
<figure><img src="image_3.png" /><figcaption>Fig. 1. 2000 h of rainfall in France (49° N 0° E). Marks are here defined as the exceedances of the u = 0.5 mm.h⁻¹ threshold.</figcaption></figure>

If the linear influence of past occurrences directly yields the intensity, as in Eqs. (15) and (16), the Hawkes process is said to be linear. Instead in a non-linear Hawkes process, the linear form gives a function of the intensity (Hawkes, 2018). In this work, the generalized logistic type II distribution is used as link function to obtain a non-linear Hawkes process.

Combining the impact function and the link function, the intensity of the marked non-linear Hawkes process is given by

$$ \lambda(t) = F_{\omega} \left[ \mu + \sum_{i|i<t} \exp \{\gamma - \gamma(t-i)\} M(i) \right]. \quad (18) $$

$F_\omega$ is the cdf of the generalized logistic type II, given by

$$ F_{\omega}(x) = 1 - \frac{\exp\{-\omega_3(x - \omega_1)/\omega_2\}}{[1 + \exp\{-(x - \omega_1)/\omega_2\}]^{\omega_3}}, \quad (19) $$

where $\omega_1$ is a location parameter, $\omega_2 > 0$ is a scale parameter and $\omega_3 > 0$ is a shape parameter. The exponential kernel in Eq. (18) was simplified compared to Eq. (16) by removing the $\gamma_1$ parameter, since it becomes redundant with the extra flexibility of $F_\omega$.

## 3.2. Two-scale intensity function

Fig. 1 presents 2000 hours of rainfall in France (ERA5 total precipitation at 49° N 0° E). Some occurrences are isolated, but most are clustered in rainfall events spanning up to 10 h in this example. These multi-hour rainfall events are not randomly distributed in time but are also regrouped in second-order clusters separated by dryer periods, as can be clearly seen during the month of April in Fig. 1. The first-order clusters in the time series are an artifact due to how rainfall is discretely observed, but the second-order clusters are caused by the low atmospheric pressure events spanning a few days, during which several rainfall events can occur. Similarly the long dry periods separating the second-order clusters are caused by high atmospheric pressure events. A self-exciting Hawkes process with a single kernel, such as the one of Eq. (18), can only model the first-order of clustering in the time series. The arrival of these clusters would be controlled by the base intensity $\mu$, resulting in them being distributed according to a Poisson process. Instead the arrivals of first-order clusters must also be clustered in time, to model the high and low pressure events.

To model this second-order clustering the base intensity $\mu$ is replaced by a second kernel which affects the probability of arrival of first-order clusters. This two-scale intensity function is given by

$$ \lambda(t) = \mu_{\beta} \exp \left[ \beta \left\{ \mu_{\alpha} - \mu_{\alpha}^{*}(t) \right\} \right] + \sum_{i|i<t} f(t-i), \quad (20) $$

where $\mu_\beta > 0$, $\beta > 0$, and

$$ \mu_{\alpha}^{*}(t) = \mu_{\alpha} \left[ 1 - \sum_{i|i<t} \exp \{\alpha - \alpha(t-i)\} \right], \quad (21) $$

where $\mu_\alpha > 0$ and $\alpha > 0$. The component $\mu_\alpha^*(t)$ is inhibited by occurrences, but overall the second kernel also acts a self-exciting part of the intensity.

Combined with the impact and the link functions, the intensity function is given by

$$ \lambda(t) = F_{\omega} \left[ \mu_{\beta} \exp \left[ \beta \left\{ \mu_{\alpha} - \mu_{\alpha}^{*}(t) \right\} \right] + \sum_{i|i<t} \exp \{\gamma - \gamma(t-i)\} M(i) \right], \quad (22) $$

with ten parameters.

## 4. Serial correlation of rainfall amounts

The Hawkes process described in the previous section models the serial correlation of occurrences but does not model the marks serial correlation. Reproducing the serial correlation of marks is necessary since rainfall extremes are duration-dependent. This is achieved by adding a copula to the EGP to obtain a distribution of the mark at time t conditional on past marks. Copulas are d-dimensional distributions with support on [0, 1]d modeling the dependence structure between these dimensions, to the exclusion of their marginal distributions. To model the marks serial correlation, the dimensions of the copula are the EGP marginal probabilities at times $t-0, ..., t-d+1$.

### 4.1. Conditional EGP with copulas

In the context of the EGP, adding a copula can be seen as applying a second probability integral transform to the GP, on top of the first transform of the EGP. In the case of a pair-copula (i.e., two-dimensional copula), this is done through its conditional distribution, also called the h-function (Aas et al., 2009), given by

$$ h(b|a) = F(b|a) = \frac{\partial C(a, b)}{\partial a}, \quad (23) $$

where C is the cdf of the copula. Note that despite the lowercase notation, the h-function is a cdf, here giving the probability of b given a, both having uniform margins.

The conditional cdf of the EGP is given by

$$ \begin{aligned} W(m_t|m_{t-1}) &= h[F_V\{G(m_t)\}|F_V\{G(m_{t-1})\}] \\ &= h[K(m_t)|K(m_{t-1})], \end{aligned} \quad (24) $$

where $m_t$ and $m_{t-1}$ are the marks at time t and $t-1$, respectively. The corresponding pdf is given by

$$ \begin{aligned} w(m_t|m_{t-1}) &= c[F_V\{G(m_t)\}, F_V\{G(m_{t-1})\}] f_V\{G(m_t)\} g(m_t) \\ &= c[K(m_t), K(m_{t-1})] k(m_t), \end{aligned} \quad (25) $$

where c is the pdf of the pair-copula.

Since rainfall time series are intermittent, not every mark at time t has a non-zero mark at $t-1$. To evaluate the likelihood of a Hawkes process marked with a copula, the pdf $w(m_t|m_{t-1})$ is used for marks with a previous non-zero mark, and $k(m_t)$ otherwise, the latter being the unconditional EGP. Likewise simulation is done with the quantile functions $W^{-1}$ (with the inverse h-function) and $K^{-1}$.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of a four dimension canonical vine copula modeling the serial dependence up to lag three. The vine copula is organized in the trees T₁ to T₃ from top to bottom. The edges between the boxes represent the pair-copulas (for example the edge 23|1 is the pair-copula between dimensions 2 and 3, conditional on dimension 1). The unconditional pair-copulas are in T₁, all linked to the first dimension of the mark at time t. The matrix representation of the vine copula is displayed on the right.</figcaption></figure>

## 4.2. Canonical vine copula for intermittent time series

The conditional EGP can also be constructed with a copula in more than two dimensions to directly model the serial correlation for longer lags. The structure of dependency between these dimensions could be too restricted with a single parametric copula, so a vine copula is used instead. Vine copulas are copulas in three or more dimensions built from an assemblage of pair-copulas (Aas et al., 2009). Any combination of parametric pair-copulas can be used for the $d(d-1)/2$ pair-copulas of a $d$-dimensional vine copula, which offers great flexibility. The structure of a vine copula is organized in so called trees, with unconditional pair-copulas in the first tree and conditional ones on subsequent trees (see Fig. 2). For vine copulas in four or more dimensions, different structures are possible depending on which dimensions are linked by unconditional or conditional pair-copulas.

The canonical vine copula structure has the particularity of having one dimension linked to all the other dimensions with unconditional pair-copulas (Fig. 2). By using the canonical structure to model marks serial correlation, the mark at time $t$ can have an unconditioned pair-copula with marks at times $t-1, ..., t-d+1$. The canonical vine copula is thus similar to a Markov process of order $d-1$.

Another property of canonical vine copulas that makes them suitable for time series is subsetability. Since the rainfall is intermittent, it is possible that not all the timesteps $t-1, ..., t-d+1$ have a non-zero mark. In that case the full $d$-dimensional canonical vine copula cannot be used. However, as long as the first dimension has a non-zero mark (which will always be the case when the copula is used), the copula can be subsetted to a smaller canonical vine copula. Fig. 3 gives the example of a subset of the canonical vine copula of Fig. 2 when one dimension is removed, in that case if there is a zero mark at time $t-2$. The remaining pair-copulas still constitute a valid canonical vine copula. These subsets are used both to compute the likelihood and to simulate. The vine copula reduces to a pair-copula if only one lag timestep has a non-zero mark, and is removed if all lags have zero marks.

Let $w_d(m_t)$ be the pdf of an EGP conditional on a canonical vine copula subsetted according to the non-zero $m_{t-1}, ..., m_{t-d+1}$ marks. In other words, $w_d(m_t)$ can range from $k(m_t)$ if every $m_{t-1}, ..., m_{t-d+1}$ marks are zero, up to an EGP conditional on the full $d$-dimensional canonical

vine copula if every $d-1$ previous marks are non-zero, with every valid subsets of the canonical vine copula in between. In the case of the full vine copula, the conditional EGP pdf is given by

$$w_d(m_t|m_{t-1}, \dots, m_{t-d+1}) = c_d\{K(m_t), \dots, K(m_{t-d+1})\} k(m_t), \quad (26)$$

where $c_d$ is the pdf of the $d$-dimensional vine copula. The corresponding cdf is

$$W_d(m_t|m_{t-1}, \dots, m_{t-d+1}) = C_{1|2,\dots,d}\{K(m_t)|K(m_{t-1}), \dots, K(m_{t-d+1})\}, \quad (27)$$

where $C_{1|2,\dots,d}$ is the conditional distribution of the first dimension of the $d$-dimensional canonical vine copula given the values of the other dimensions.

A convenient representation of vine copulas is a square matrix indicating the set of conditional dimensions for each pair-copula (Morales Naples, 2009). The matrix representation of a 4-dimension canonical vine copula is

$$\begin{bmatrix} 1 & 1 & 1 & 1 \\ 2 & 2 & 2 & 2 \\ 3 & 3 & 3 & 4 \end{bmatrix}.$$

This matrix is read by column, with the number on the diagonal and another number in the same column indicating a pair-copula linking these two dimensions. Each other number above them in the same column indicates the conditioning dimensions of this pair-copula. For instance, the numbers 2 and 4 in bold in the matrix indicate a pair-copula, which is conditioned on the first dimension denoted in italic, but not on the third dimension (this is the 24|1 pair-copula in Fig. 2). When the triangular matrix is in the upper corner, the subsets of a canonical vine copula are obtained by removing the rows and columns of the removed dimensions (which can be any dimensions except for the first one). Fig. 3 shows the matrix representation of the subsetted 4-dimensional canonical vine copula on the right, with the third row and column emptied. By removing these empty rows and columns, the triangular matrix of a smaller vine copula is obtained. The pair-copulas in the subset keep the same parametric form and parameter values as in the full vine copula.

The parametric pair-copula considered in the vine copula must cover different types of dependency to ensure that the vine copula is
<figure><img src="image_5.png" /><figcaption>Fig. 3. Subset of the canonical vine copula of Fig. 2 when the third dimension is removed.</figcaption></figure>

flexible enough. This includes parametric pair-copulas with asymptotic tail dependence or asymptotic tail independence, and pair-copulas with negative dependence. The parametric pair-copulas considered in this application are the Gaussian, Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7 and BB8 pair-copulas, which all have closed form h-functions (Joe, 2015). Rotations by 90, 180 and 270° of the pair-copulas are also considered to further increase the possibilities.

# 5. Model inference

## 5.1. Fitting by maximum likelihood

The model is fitted by maximum likelihood estimate. This could be done in one step in principle, but since the model has many parameters it becomes unpractical beyond a certain data size. For the application on 40 years of hourly data, the fitting was carried out in two steps. First the EGP parameters were estimated, then the Hawkes process and copula parameters. This two-steps fitting between marginal distributions and dependence structure is common in copula modeling (Größer and Okhrin, 2022). Here both the copula and the Hawkes process are part of the dependence structure.

In the case of rainfall, the model can in principle model the whole range with an EGP threshold at 0, provided that the EGP is sufficiently flexible. However rainfall time series are discretized due to instrument precision, with values rounded to 0.1 mm.h⁻¹ for example, and sometimes censored below a certain threshold. This issue was already noted by Naveau et al. (2016) in the case of the EGP, but also causes issues for the copula, as the discretization can artificially increase the lower tail dependence. To circumvent this, a low threshold can be used instead of having it to 0. In that case, the marks are the threshold exceedances. This threshold must remain low enough to make no significant difference in the context of flood risk assessment.

Let $u \ge 0$ be the threshold and $x$ a time series of $T$ timesteps. The marks $m$ are defined as

$$m(t) = \max(x_t - u, 0), \qquad (28)$$

for $t \in 1, \dots, T$.

The full model log-likelihood is given by

$$\log(\mathcal{L}) = \sum_{t|m_t>0} \log\{\lambda(t)\} + \log\{w_d(m_t)\} + \sum_{t|m_t=0} \log\{1 - \lambda(t)\}. \qquad (29)$$

If the two-step optimization is used instead of Eq. (29), the log-likelihood is given by $\log(\mathcal{L}) = \log(\mathcal{L}_1) + \log(\mathcal{L}_2)$, with

$$\log(\mathcal{L}_1) = \sum_{t|m_t>0} \log\{k(m_t)\}, \qquad (30)$$

to estimate $\hat{\kappa}$ the parameters of $K$ (including the lower truncation bound $l$ if the EGP transform $f_V$ in $K$ is Eq. (8) or (12)). The second term is then fitted conditionally on the estimates of the first term, with

$$\log(\mathcal{L}_2|\hat{\kappa}) = \sum_{t|m_t>0} \log\{\lambda(t)\} + \log[c_d\{K(m_t|\hat{\kappa}), \dots, K(m_{t-d+1}|\hat{\kappa})\}k(m_t|\hat{\kappa})] \\ + \sum_{t|m_t=0} \log\{1 - \lambda(t)\}. \qquad (31)$$

where $c_d$ is subsetted according to the $m_{t-1}, \dots, m_{t-d+1}$ non-zero previous marks.

Since the intensity of the Hawkes process at time $t$ depends on all past observations and that their order in time matters, the evaluation of the likelihood cannot be paralleled in principle. However with long time series at high frequency, evaluating the likelihood on a few separate long blocks of the time series should not make any significant difference in the inference, since the length of these blocks remains much larger than the duration for which the second-order kernel in Eq. (21) is significantly above 0. The small error induced by the parallelism could be reduced by having the block slightly overlapping, so that the intensity and copula at the beginning of each block could still be computed on some previous observations, but without computing twice the likelihood of the overlapping observations. This correction has not been implemented because it would make an insignificant difference in this application, given the length of the time series.

## 5.2. Simulation

Risk assessment from the model is obtained through simulation with the algorithm 1, where $W_d^{-1}$ is the quantile function of the conditional EGP, with the copula subsetted according to the non-zero previously simulated values in $m_{t-1}, \dots, m_{t-d+1}$. The first $d-1$ values could either be set to 0, or the copula can be subsetted to accommodate for the beginning of the simulation.

In the case of a two-dimensional copula, $W^{-1}$ uses an inverse h-function. Instead if the EGP is combined with a canonical vine copula in $W_d^{-1}$, a random value must be generated from the conditional cdf of the
**Data:** Length *T* of the time series to simulate.

**Result:** A time series of marks *m* of length *T*.

for *t* ← 1 to *T* do

p = λ(t)
a ~ U(0, 1)
if p > a then
  b ~ U(0, 1)
  m_t = W_d⁻¹(b|m_{t-1}, ..., m_{t-d+1})
else
  m_t = 0
end

end

**Algorithm 1:** Simulation from the discrete Hawkes process marked by a conditional EGP.

vine copula $C_{1|2,...,d}$. The standard simulation algorithm for a canonical vine copula is recursive, with random values being generated from the first to the last dimension, conditionally on the values generated for the previous dimensions (Joe, 2015). Therefore, it would only be possible to adapt this algorithm to obtain the conditional distribution of the last dimension $C_{d|1,...,d-1}$, but the values simulated for the first dimension would be uniformly distributed. Instead, a value is generated by rejection sampling in the first dimension of the pdf of the canonical vine copula, with the values for the other dimensions given by $K(m_{t-1}), ..., K(m_{t-d+1})$.

Unlike the likelihood evaluation, the algorithm 1 cannot be paralleled. However the risk assessment from the model would be obtained from several simulated time series, so different simulations could be ran in parallel.

# 6. Application

## 6.1. Dataset

An application of the model is presented with a time series of 40 years of hourly rainfall. This time series is the total precipitation of the ERA5 reanalysis (Hersbach et al., 2018), from 1980 to 2019, at the location 49° N 0° E. This location is in mainland France, with an oceanic climate in the Köppen classification. Rainfall in this climate type is frequent, including storms which are responsible for heavy rainfall events. Snowfall is rare in oceanic climates, and does not significantly contribute to flood risk as it would in other climates, so not accounting for this different type of precipitation in the model is not an issue for this application. The applicability of the model is therefore restricted to areas where the only type of precipitation contributing to flood risk is rainfall.

ERA5 total precipitations have numerical artifacts in very low values (Fig. 4), which is a different issue from the discretization of station data mentioned in Section 5.1, but could similarly impair the inferences. Therefore, a low threshold of $u = 0.5 \text{ mm.h}^{-1}$ is used, which avoids the eventual spurious correlation in the lower tail of the copula and reduces computational requirement for both parameter estimation and simulation.

Plotting the 40 years of precipitation along the day of the year reveals strong seasonal variability (Fig. 5), as it is expected for rainfall in most regions. This variability is also skewed, with the average values increasing faster during the first half of the year than they decrease during the second half (red curve on Fig. 5). In this case, upper values are higher during summer than winter. The seasonal effect is considered for the scale parameter $\sigma$, with a skewed sine wave. This season-dependent parameter is given by

$$ \sigma(t) = \begin{cases} \sigma_0 + \eta \psi^{-1} \tanh \left\{ \frac{\psi \sin(s_t - \phi)}{1 - \psi \cos(s_t - \phi)} \right\} & \text{if } \psi \neq 0, \\ \sigma_0 + \eta \sin(s_t - \phi) & \text{if } \psi = 0, \end{cases} \quad (32) $$

where $s_t$ in the time of the year mapped to [0, 2π], $\sigma_0 > 0$, $\eta > 0$ is the amplitude of the seasonal effect, $\psi$ controls the skewness of the sine wave, and $\phi$ is the phase.

## 6.2. EGP evaluation

The EGP-beta₃ with a seasonally varying scale parameter $\sigma$ is strongly favored over a time-invariant model by the Akaike information criterion (AIC) and Bayesian information criterion (BIC) (not shown). Fig. 6 presents the histogram of the marks (i.e., exceedances above $u$) along the EGP-beta₃ density for the minimum and maximum values of $\sigma(t)$. The difference between the two curves shows how a stationary model would underestimate the highest rainfall marks, compared to this seasonal model.

The shape parameter $\xi$ (i.e., tail index) of an extreme value distribution is the most critical one, as it controls the heaviness of the upper tail. This parameter is generally positive for rainfall, which corresponds to a tail heavier than the exponentially tailed special case $\xi = 0$. Serinaldi and Kilsby (2014) found that the value of $\xi$ for the GP should asymptotically be in the interval (0.061, 0.097) with infinite sample size. In this work, the EGP-beta₃ is fitted on almost 17 000 exceedances above the $u = 0.5 \text{ mm.h}^{-1}$ (4.8% of the 40 years of hourly time series), so the sample size is quite large from the perspective of the GP, even though only a fraction of these observations are extreme. The estimate of the tail index is $\hat{\xi} = 0.076$, which is in the range given by Serinaldi and Kilsby (2014). The estimates of the other parameters are indicated in Table 2. If a constant $\sigma$ is used instead of the seasonal $\sigma(t)$ of Eq. (32), the estimate becomes $\hat{\xi} = 0.148$, which could be an overestimation caused by $\xi$ compensating for the lack of flexibility in $\sigma$. For $u = 0.5 \text{ mm.h}^{-1}$, these results regarding the tail index are similar with the three other EGPs presented in Section 2: the EGP-power, EGP-beta₁ and EGP-genbeta.

The four EGPs are fitted with a lower threshold of $u = 0.05 \text{ mm.h}^{-1}$ to test if the extra flexibility of the EGP forms with more parameters is useful. In each case, the EGP is fitted both with a constant and seasonal parameter $\sigma$. Table 1 presents the estimated tail index $\hat{\xi}$ of each fitting. With this lower threshold, the lack of flexibility of the one-parameter EGP-power leads to a largely overestimated $\hat{\xi}$ of about 0.44, regardless of a constant or seasonally varying scale parameter $\sigma$. The EGP-beta₁ also lacks flexibility, with $\hat{\xi} = 0.178$ for both forms of $\sigma$. This EGP form is the only one with worse values of the information criteria AIC and BIC for the seasonal $\sigma(t)$. The fitting quality improves significantly with the extra flexibility provided by the three- and four-parameter EGP-beta₃ and EGP-genbeta, respectively, with $\hat{\xi}$ estimates of about 0.09 falling in the (0.061, 0.097) interval of Serinaldi and Kilsby (2014). These two EPGs also have lower values of AIC and BIC (note that Table 1 only reports the number of parameters of the probability transform of the EGP $F_v$, but the information criteria are computed according to the number of parameters of $K$, including the extra parameters when the seasonal $\sigma(t)$ is used). For this time series and a $u = 0.05 \text{ mm.h}^{-1}$ threshold, the EGP-genbeta does not significantly improve the fitting over the EGP-beta₃, but it can be assumed that the more flexible EGP-genbeta could be useful for rainfall time series at higher frequency than hourly. Interestingly, the results between the constant $\sigma$ and seasonal $\sigma(t)$ EGPs are similar with this lower threshold, which indicates that the more complex seasonal form does not affect the estimated $\hat{\xi}$. The purpose of this comparison of the EGPs with a lower threshold is to demonstrate the potential usefulness of the more complex ones, but the EGP-beta₃ with $u = 0.5 \text{ mm.h}^{-1}$ discussed in the previous paragraph is used for the Hawkes process in the rest of the present work.

## 6.3. Hawkes process evaluation

Whether using the full likelihood of the Hawkes process in Eq. (29) or the two-step likelihood with Eqs. (30) and (31), the parameters of the unconditional EGP $K$ can be estimated separately with Eq. (31)
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Histogram (a) of the ERA5 hourly total precipitation time series. Empirical probabilities (b, c) of the precipitation plotted for 1 h time lag. The red line indicates the $u = 0.5$ mm.h$^{-1}$ threshold. The subplot c shows the region above the threshold in b. There are also values around $10^{-12}$ mm.h$^{-1}$, that are not displayed. Only 3000 values are displayed in b and c for readability.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** 40 years of precipitation along the day of the year. The values smoothed along the day of the year (red curve) highlight the seasonal variability, which is also visible in the upper extremes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Histogram of the marks and density of the EGP-$\beta_3$ with a seasonally varying $\sigma(t)$. The EGP density is displayed for the minimum and maximum values of $\sigma$ (which explains the apparent discrepancy with the histogram). The marks are standardized so are not in mm.h$^{-1}$, and the upper values are not displayed for readability.</figcaption></figure>

<caption>**Table 1**</caption>

Inferences for four different EGP distributions with a lower threshold of $u = 0.05$ mm.h$^{-1}$.

<table><thead><tr><th>F<sub>v</sub> form</th><th>F<sub>v</sub> parameters</th><th>seasonal σ</th><th>ê<sub>c</sub></th><th>AIC</th><th>BIC</th></tr></thead><tbody><tr><td rowspan="2">power</td><td rowspan="2">1</td><td>no</td><td>0.438</td><td>143 072</td><td>143 104</td></tr><tr><td>yes</td><td>0.448</td><td>142 945</td><td>143 004</td></tr><tr><td rowspan="2">EGP-$\beta_1$</td><td rowspan="2">1</td><td>no</td><td>0.178</td><td>141 655</td><td>141 686</td></tr><tr><td>yes</td><td>0.178</td><td>141 663</td><td>141 722</td></tr><tr><td rowspan="2">EGP-$\beta_3$</td><td rowspan="2">3</td><td>no</td><td>0.096</td><td>141 514</td><td>141 564</td></tr><tr><td>yes</td><td>0.109</td><td>141 418</td><td>141 495</td></tr><tr><td rowspan="2">genbeta</td><td rowspan="2">4</td><td>no</td><td>0.091</td><td>141 514</td><td>141 573</td></tr><tr><td>yes</td><td>0.094</td><td>141 416</td><td>141 502</td></tr></tbody></table>

for an exploratory analysis guiding the selection of the pair-copulas forms in the canonical vine copula. Software such as the R package *rvinecopulib* (Nagler and Vatter, 2023) can automatically select the best pair-copula for a dataset among different parametric copulas, including their rotations. The probabilities of the non-zero marks are computed with the estimates of the EGP with $K(m|\hat{k})$, then the best parametric form, including rotations, can be selected for each $d(d-1)/2$ pair-copulas of the canonical vine copula. As this automatic selection is only possible on the full $d$-dimensional vine copula in *rvinecopulib*, it is only done on the observations at time $t$ having non-zero marks up to $t-d+1$. This restriction reduces the amount of data usable as $d$ increases. This exploratory analysis is performed for $d = 4$ to consider the lags up to $t-3$ in the vine copula. The unrotated BB8 pair-copula is selected for each pair-copula, among the nine forms and their rotations considered (see Section 4.2 for the list of parametric pair-copulas tested). The BB8 pair-copula is asymmetric between its lower and upper tail, with both tails being asymptotically independent (Joe, 2015).

After this exploratory analysis, the parameters of the Hawkes process and a 4-dimensional canonical vine copula for $W_d$ are jointly estimated with the likelihood in Eq. (31). In this application $d = 4$ because modeling up to lag $t-3$ gave better results than smaller vine copulas, while allowing stable parameter estimation. A principled size selection for the vine copula was not studied. Table 2 presents the estimates for the 30 parameters of the model, regrouped into its subcomponents.

Fig. 7 presents the second-order $\alpha$ and first-order $\gamma$ kernels of the two-scale Hawkes intensity function (Eqs. (20) and (21)). These kernels behave as expected, since the second-order one stays significantly above 0 for a much longer duration than the first-order one. The parametrization of (18) with $\gamma$ appearing twice forces a value of 1 at lag 1, which simplifies the model fitting. A similar parametrization is used for the second-order kernel $\alpha$ (Eq. (21)). The memory kernels are infinite, but in practice they drop to almost 0 at some point. To reduce the computational requirements, the memory kernels were censored for values below $10^{-4}$. Forcing them to have a value of 1 at lag 1 also ensures that this censoring stays consistent for different parameter
<table><caption>**Table 2**
Estimates of the 30 parameters of the model.</caption><thead><tr><th colspan="6">GP with seasonal σ(t) (Eqs. (1) and (32))</th></tr><tr><th>σ<sub>0</sub></th><th>η</th><th>ψ</th><th>φ</th><th>ξ</th></tr></thead><tbody><tr><td>2.073</td><td>0.575</td><td>0.382</td><td>2.112</td><td>0.076</td></tr><tr><th colspan="6">EGP-beta<sub>3</sub> transform F<sub>V</sub> (Eq. (8))</th></tr><tr><th>κ<sub>1</sub></th><th>κ<sub>2</sub></th><th>l</th><th></th><th></th></tr><tr><td>0.679</td><td>0.359</td><td>0.093</td><td></td><td></td></tr><tr><th colspan="6">Hawkes intensity (Eqs. (20) and (21))</th></tr><tr><th>μ<sub>α</sub></th><th>α</th><th>μ<sub>β</sub></th><th>β</th><th>γ</th></tr><tr><td>0.051</td><td>0.043</td><td>6.549</td><td>0.051</td><td>2.426</td></tr><tr><th colspan="6">impact and link functions (Eqs. (17) and (19))</th></tr><tr><th>δ<sub>1</sub></th><th>δ<sub>2</sub></th><th>ω<sub>1</sub></th><th>ω<sub>2</sub></th><th>ω<sub>3</sub></th></tr><tr><td>1.983</td><td>0.273</td><td>6.524</td><td>0.005</td><td>0.002</td></tr><tr><th colspan="6">canonical vine copula (all BB8 pair-copulas)</th></tr><tr><th>θ<sub>12</sub></th><th>θ<sub>13</sub></th><th>θ<sub>14</sub></th><th>θ<sub>23|1</sub></th><th>θ<sub>24|1</sub></th><th>θ<sub>34|12</sub></th></tr><tr><td>[2.587]<br>[0.832]</td><td>[6.277]<br>[0.238]</td><td>[5.246]<br>[0.163]</td><td>[4.204]<br>[0.697]</td><td>[6.700]<br>[0.306]</td><td>[2.683]<br>[0.684]</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Second-order α and first-order γ kernels of the two-scale Hawkes intensity. Note that these kernels are infinite, here only displayed for lags up to 100 and 10, respectively.</figcaption></figure>

values (the censoring at a given threshold would be relatively higher or lower if the kernel had a value lower or higher, respectively, than 1 at lag 1).

Fig. 8 presents 2 000 hours of simulated values (slightly less than three months). The top panel of Fig. 8 is the second-order clustering component $\mu_a^*$ of the intensity (Eq. (21)). The middle panel is the full intensity $\lambda$ controlling the probability of occurrence. The bottom panel is the marks (exceedances of the $u = 0.5 \text{ mm.h}^{-1}$ threshold). Occurrences inhibit the second-order component $\mu_a^*$, which in turn increases the probability of first-order clusters. In other words, the inhibition of $\mu_a^*$ permits the succession of first-order clusters, consequently regrouping them in second-order cluster. These two orders of clustering are visible during the first half of the simulations. The second-order component stays at its maximum of $\mu_\alpha = 0.051$ for long periods in the second half of the simulation, during which the probability of occurrence without past influence reaches its lowest value with $\mu_\beta = 6.524$ (recall that this value is not a direct probability since the intensity in Eq. (22) is nonlinear). These periods of high $\mu_a^*$ correspond to dry spells in the marks (which are here only relatively dry since the threshold $u$ is above 0). Note that a particularly dry period in the simulations is displayed here to illustrate the two-scale intensity.

Fig. 9 presents the density of the three unconditional pair-copulas of the first tree $T_1$ of the 4-dimensional canonical vine copula (see Fig. 2). The correlation is higher between times $t$ and $t-1$ modeled by the pair-copula $c_{12}$, than the subsequent lags $t-2$ and $t-3$ modeled by $c_{13}$ and $c_{14}$, respectively.

Fig. 10 presents the autocorrelation functions of the observations and 40 years of simulation. The autocorrelation is slightly underestimated at lag 1, then has a slight systematic overestimation from lag 4 and higher.

The IDF model is fitted on block maxima of both observations and simulations (40 years of hourly time series in both cases), to assess whether the Hawkes process reproduces the relation between amount of rainfall, duration and probability. Here the IDF model is only used as a diagnostic for the WG, and not as risk assessment. The parametrization of the IDF model is first selected based on the best fit on the observations. As yearly block maxima are used, and for simplicity, the seasonal cycle is omitted. The parametrization of the IDF model is

$$G(z, d|\bar{\mu}, \sigma_0, \xi, \theta, \eta) = \exp \left[ - \left\{ 1 + \xi \left( \frac{z}{\sigma_0(d + \theta)^{-\eta}} - \bar{\mu} \right) \right\}^{-1/\xi} \right], \quad (33)$$

where $z$ are the block maxima for the duration $d$, $\bar{\mu} > 0$, $\sigma_0 > 0$, $\xi \in \mathbb{R}$, $\xi \neq 0$, $\theta \ge 0$ and $\eta \in (0, 1]$. See Ulrich et al. (2020) for details on this IDF parametrization, which is implemented in the R package *IDF*. The same IDF parametrization is then fitted on the simulations from the Hawkes process. Since the vine copula in $W_d$ (Eq. (27)) is the primary component of the Hawkes process accounting for the duration of events, a simpler Hawkes process without copula was fitted and simulated from as comparison, with the marks modeled only by an unconditional EGP *K* (Eq. (3)). Note that $d$ in $W_d$ refers to the vine copula dimension, while it refers to durations in the IDF model. The three IDF models are fitted for durations $d \in 1, \dots, 120$ hours. Table 3 presents their estimated parameters. The tail index is estimated at $\xi = -0.05$ for simulations from the model with vine copula, which is close to the estimate of $-0.07$ for the observations. However the tail index is lower at $-0.16$ for the simulations without vine copula. When IDF curves are plotted on a double log scale, the parameters $\theta$ and $\eta$ (Eq. (33)) control the curvature of the IDF curves for lower durations and the slope for longer durations, respectively (Ulrich et al., 2020). For both IDF models from the simulations, $\hat{\theta} = 0$, whereas it is 1.08 for the observations. Likewise, both simulation models have an estimated $\hat{\eta}$ of about 0.6, while it is 0.8 for the observations.

Fig. 11 presents the IDF curves for the yearly nonexceedance probabilities {0.9, 0.99, 0.999}. The curves from the simulations of the model without vine copula are systematically lower than the observations' curves, whereas the curves of the simulations from the full model with a vine copula are at the same level as the observations' curves, which is explained by the tail index $\xi$ (Table 3). This indicates that the canonical vine copula in the conditional mark distribution $W_d$ significantly improves the Hawkes process. If the autocorrelation of marks was not accounted for and only the autocorrelation of occurrences was modeled (i.e., in a Hawkes process with a mark distribution not conditional on past marks), the risk assessment obtained from the simulations would be underestimated. The observations' IDF curves display a curvature for the lower durations (here, under 10 h), whereas there is none for both simulations, which is explained by the difference in the estimates of $\theta$. Likewise observations' and simulations' curves have a different slope for higher durations, which is explained by the different estimates of $\eta$, and results in the probabilities being overestimated for longer durations in the simulations from the full Hawkes process. This overestimation at longer durations is probably linked to the overestimation in the autocorrelation seen in Fig. 10. These discrepancies between observations and simulations from the full model indicate that some patterns are not reproduced, and the current model is not reliable for longer durations.

## 7. Discussion and perspectives

The two major differences between the Hawkes process of the present work and the one of Li et al. (2021) are the addition of the vine copula to obtain a conditional distribution of the mark modeling their autocorrelation, and the use of a two-scale intensity function.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Second-order clustering component $\mu_a^*$ (top panel) and full intensity $\lambda$ (middle panel) of Eq. (22), and marks (bottom panel) for 2000 hours of simulation.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Density of the unconditional pair-copulas in the first tree $T_1$ (top row of Fig. 2), modeling the dependence between the mark at $t$ and previous marks at $t-1$ to $t-3$ from left to right, with the pair-copulas $c_{12}$ to $c_{14}$, respectively. The three are unrotated BB8 pair-copulas.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Autocorrelation functions of the observations (exceedances of the threshold $u = 0.5 \text{ mm.h}^{-1}$) and simulations (both 40 years of hourly rainfall).</figcaption></figure>

Their Hawkes process was applied to daily river flow, with the sign of the river flow reversed, so that exceedances of the threshold (i.e., non-zero marks) corresponded to low flow values and hence to droughts (as opposed to this work where zero marks correspond to dry periods). In the model of Li et al. (2021), the scale parameter $\sigma$ of the EGP for the

<caption>**Table 3**</caption>

Estimates of the IDF models for the observations, the simulations without vine copula and the simulations from the full Hawkes process with a 4-dimensional canonical vine copula. Note that $\theta$ is not a pair-copula parameter here, but a parameter of Eq. (33).

<table><thead><tr><th>Model</th><th>μ̃</th><th>σ<sub>0</sub></th><th>ξ</th><th>θ</th><th>η</th></tr></thead><tbody><tr><td>obs</td><td>3.87</td><td>2.67</td><td>-0.07</td><td>1.08</td><td>0.80</td></tr><tr><td>sim<sub>no cop</sub></td><td>3.82</td><td>1.88</td><td>-0.16</td><td>0.00</td><td>0.63</td></tr><tr><td>sim<sub>4d cvc</sub></td><td>5.12</td><td>1.00</td><td>-0.05</td><td>0.00</td><td>0.64</td></tr></tbody></table>

marks is linearly dependent on the self-exciting part of the intensity function. This varying $\sigma$ and the autocorrelation of occurrences will result in some autocorrelation of the whole time series in simulation, but not enough since the autocorrelation of exceedances (i.e., marks) is not explicitly modeled. In the dataset of Li et al. (2021), clusters of events almost only happen during the summer months, with only a few isolated non-zero values during winter. They accounted for this by having a partially deterministic intensity function, with the self-exciting component removed for winter months, during which the non-zero marks become Poisson distributed. It appears logical to have the yearly cycle (i.e., seasonal variability) as a deterministic component of the model, as it is the case for the scale parameter $\sigma(t)$ (Eq. (32)) in the proposed model, but the resulting EGP and Hawkes intensity remained fully stochastic in the present work. In the model of Li et al. (2021), clusters of events can only happen from May to November, with this cutoff not being estimated through the likelihood but a priori selected. The two-scale intensity (Eqs. (20) and (21)) would not be appropriate for the river flow of Li et al. (2021), because the second-order of
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** IDF curves for the observations (obs, solid lines), simulations from the Hawkes process without vine copula ($sim_{no\ cop}$, short dash) and simulations from the full model with a 4-dimensional canonical vine copula ($sim_{4d\ cvc}$, long dash). The models are adjusted for durations of 1 to 120 h. Both axes are in log.</figcaption></figure>

clustering in the dataset is the yearly cycle, which must be considered deterministic. However, the Hawkes process can remain fully stochastic if the yearly cycle is integrated in the model as a covariate of the intensity.

Only two clustering orders were considered in the proposed model, with the first-order being an artifact of the discrete time series and the second one being linked to the short-term meteorological variability. Higher orders of clustering are most probably present in the data. Precipitations in France are affected by the atmospheric pressure patterns of the North Atlantic Oscillation (NAO) (Hurrell et al., 2003), on a long-term time scale and large spatial scale. Massei et al. (2017) showed that taking into account large-scale and low frequency (i.e., lower than the yearly cycle) atmospheric pressure patterns improves the modeling of local-scale precipitation. For many regions, the El Niño-Southern Oscillation (ENSO) would be another climate pattern to consider for precipitation modeling (Sun et al., 2015; Ouarda et al., 2021). Such low frequency phenomenon could potentially affect the probability of rainfall, and thus could not be incorporated in the Hawkes process intensity. They could also not affect occurrences but affect the amount of rainfall via the marks distribution, or affect both occurrences and amounts. A WG using climate indices as covariate would require modeling them, then to either simulate or forecast the indices and subsequently conditionally simulate the rainfall. Lee et al. (2023) compared different forecasting methods for the ENSO and Pacific Decadal Oscillation indices.

The EGP-beta₃ and EGP-genbeta are defined to respect the condition of Eq. (5) proposed by Gamet and Jalbert (2022), but a more flexible version of these two EGPs could be defined by having the upper truncation boundary *u* as a parameter to estimate. This modification would result in the EGP only respecting the condition of Eq. (4). In the case of rainfall, the extra flexibility is required for the lower tail, and even more so as the threshold is lowered close to 0, or as the observation frequency increases, so for the application this extra flexibility for the upper tail was not required.

Despite the 4-dimensional canonical vine copula used in the mark distribution $W_d$, the autocorrelation function of the simulations did not fully reproduce the autocorrelation of the observations (Fig. 10). The choice of vine copula dimension and its parametric pair-copulas significantly affected the autocorrelation function of the simulations, so the lack of total adequacy of the vine copula must be the main reason

for this discrepancy. Although the IDF diagnostics (Fig. 11) show that the vine copula improves the model, discrepancies between simulations and observations are also apparent here. A characteristic of rainfall data is the peak of density at the lower tail (Fig. 6), which could disproportionately affect the vine copula. This issue can be compounded by the artifact that can be present in the lower values of rainfall time series, either caused by the discretization of the observations or by numerical artifacts in reanalysis data (as in the ERA5 total precipitation). This could be corrected by a biased likelihood for the copula, giving more weight to the upper rainfall values.

A limitation of the current model is the a priori selection of the pair-copulas in the canonical vine copula, only justified by an exploratory analysis. Instead, the parametric pair-copulas could be selected jointly with the rest of the model by a reversible jump Markov chain Monte Carlo (RJMCMC) algorithm (Green, 1995; El Adlouni and Ouarda, 2009). This Bayesian inference algorithm swaps the parametric pair-copulas in the vine copula to give the posterior distributions of the parametric form and the parameters of each pair-copula (Min and Czado, 2010; Gruber and Czado, 2018; Chapon et al., 2023). Rotations of the pair-copulas could also be considered by a RJMCMC algorithm. Another a priori choice was the dimension of the canonical vine copula, which could also be included in a RJMCMC algorithm. Such algorithm was considered to be outside the scope of the present work and was not incorporated in the proposed model.

The model was fitted in Stan through gradient descent optimization rather than Bayesian sampling to keep the computation time low. However, the current implementation can be inferred by both methods without modification of the code. The only modification needed would be to select parameter priors allowing stable inferences by the Hamiltonian Monte Carlo algorithm used in Stan. This would allow taking into account uncertainties of both the parameters and the process by simulating from the model different time series with different samples from the parameters' posteriors.

The subsets of the canonical vine copula were applied to an intermittent time series, but they could also be useful to handle missing values, for either intermittent or continuous time series. This would allow the use of all the available data to compute the likelihood, without a prior imputation step.

The first-order kernel $\gamma$ drops from 1 at $t-1$ to 0.088 at $t-2$, then to almost 0 for further lags (Fig. 7). The second-order kernel $\alpha$ remains
significantly above 0 for a much longer time. Therefore, a simplification of the model could be to have a first order Markov process instead of an infinite kernel for the first-order of clustering, in a hybrid of Hawkes and Markov processes. Taking only *t* − 1 into account in this Markov process could allow a more complex form than a typical kernel used in a Hawkes process (exponential kernels in this work).

As mentioned in Section 3.2, the first-order of clustering in rainfall time series, of several subsequent non-zero values, does not correspond to clusters in the real world. These clusters are an artifact due to the temporal discretization of rainfall observations. If the duration of rainfall events is omitted, each event can be considered as a point in time, which can be modeled as a point process. These events are not independent, so they must be modeled by a Hawkes process rather than a Poisson process. The duration of events would be accounted for in the mark, along with the amount, and a single kernel continuous Hawkes process would be adapted. In other words, the clustering scale which has physical meaning as a real cluster in the two-scale discrete Hawkes process is the second-order one. Therefore, a single kernel discrete Hawkes process for a rainfall time series only models the artificial clustering and misses the real world clustering (granted that no other part of the model accounts for it).

8. Conclusion

The Hawkes process presented in this article is both an extreme value model, extending the framework of the GP, and a stochastic WG able to simulate time series. The distinction between these two types of models is only a matter of perspective, as inferences from an extreme value model can be obtained through simulation, and the marginal distribution of a WG needs the extreme value theory to account for the heavy upper tail of rainfall. The application of the model to an hourly time series of 40 years shows that the Hawkes process reproduces the pattern between the amount and the duration of extreme events.

The three parts of the model are an EGP distribution for the non-zero values, a Hawkes process for the intermittency of rainfall and a copula for the temporal dependence. A new EGP parametrization for hourly rainfall is proposed. A two-scale intensity function for the Hawkes process models the two scales of clustering of non-zero values in hourly rainfall time series. The temporal dependence is modeled by a canonical vine copula, which can be subsetted for timesteps without rain.

Floods and droughts are the two types of extreme events at the two tails of the rainfall distribution. Even thought the focus was put on the upper extremes throughout the article, the Hawkes process simultaneously models both tails of the distribution. The low threshold of $u = 0.5 \text{ mm.h}^{-1}$ in the application was only to circumvent numerical artifacts in the ERA5 dataset, and most stations hourly rainfall time series would be censored below a value higher than $u = 0.5 \text{ mm.h}^{-1}$ so the model is effectively applicable with a threshold of $u = 0$. Ward et al. (2020) advocate for considering droughts and floods together for risk reduction strategies, because of their many interactions. Jointly modeling them is the first step in this regard. The ability of a WG to handle both types of extremes is an advantage over the IDF model.

Code availability

The model is implemented in Stan (Carpenter et al., 2017), with the code available at https://github.com/antoinechapon/swg_param.

CRediT authorship contribution statement

**Antoine Chapon:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Formal analysis, Conceptualization. **Taha B.M.J. Ouarda:** Writing – review & editing, Supervision, Project administration, Funding acquisition. **Nathalie Bertrand:** Writing – review & editing, Supervision, Funding acquisition.

Acknowledgments

The authors thank the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Research Chairs Program, and the French Nuclear Safety and Radiation Protection Authority (ASNR) for funding this research.

Declaration of competing interest

The authors certify that they have NO affiliations with or involve-
ment in any organization or entity with any financial interest (such
as honoraria; educational grants; participation in speakers' bureaus;
membership, employment, consultancies, stock ownership, or other
equity interest; and expert testimony or patent-licensing arrangements),
or non-financial interest (such as personal or professional relationships,
affiliations, knowledge or beliefs) in the subject matter or materials
discussed in this manuscript.

Data availability

The ERA5 hourly total precipitation dataset is available at https://c
ds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=ov
review.

References

Aas, K., Czado, C., Frigessi, A., Bakken, H., 2009. Pair-copula constructions of multiple dependence. Insurance Math. Econom. (ISSN: 01676687) 44 (2), 182-198. http://dx.doi.org/10.1016/j.insmatheco.2007.02.001, URL https://linkinghub.elsevier.com/retrieve/pii/S0167668707000194.

Ahn, K.-H., 2020. Coupled annual and daily multivariate and multisite stochastic weather generator to preserve low- and high-frequency variability to assess climate vulnerability. J. Hydrol. (ISSN: 00221694) 581, 124443. http://dx.doi.org/10.1016/j.jhydrol.2019.124443, URL https://linkinghub.elsevier.com/retrieve/pii/S0022169419311783.

Beneyto, C., Aranda, J.A., Francés, F., 2023. Exploring the uncertainty of weather generators' extreme estimates in different practical available information scenarios. Hydrol. Sci. J. 68 (9), 1203-1212. http://dx.doi.org/10.1080/02626667.2023.2208754, ISSN: 0262-6667, 2150-3435 URL https://www.tandfonline.com/doi/full/10.1080/02626667.2023.2208754.

Carpenter, B., Gelman, A., Hoffman, M.D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., Riddell, A., 2017. Stan : A probabilistic programming language. J. Stat. Softw. (ISSN: 1548-7660) 76 (1), http://dx.doi.org/10.18637/jss.v076.i01, URL http://www.jstatsoft.org/v76/i01/.

Chapon, A., Ouarda, T.B., Hamdi, Y., 2023. Imputation of missing values in environ-
mental time series by D-vine copulas. Weather. Clim. Extrem. (ISSN: 22120947)
100591. http://dx.doi.org/10.1016/j.wace.2023.100591, URL https://linkinghub.
elsevier.com/retrieve/pii/S2212094723000440.

Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer, London, ISBN: 978-1-84996-874-4.

El Adlouni, S., Ouarda, T.B.M.J., 2009. Joint Bayesian model selection and parameter estimation of the generalized extreme value model with covariates using birth-death Markov chain Monte Carlo. Water Resour. Res. 45 (6), http://dx.doi.org/10.1029/2007WR006427.

Evin, G., Favre, A.-C., Hingray, B., 2018. Stochastic generation of multi-site daily precipitation focusing on extreme events. Hydrol. Earth Syst. Sci. (ISSN: 1607-7938) 22 (1), 655-672. http://dx.doi.org/10.5194/hess-22-655-2018, URL https://hess.copernicus.org/articles/22/655/2018/.

Gamet, P., Jalbert, J., 2022. A flexible extended generalized Pareto distribution for tail estimation. Environmetrics 33 (6), http://dx.doi.org/10.1002/env.2744, ISSN: 1180-4009, 1099-095X URL https://onlinelibrary.wiley.com/doi/10.1002/env.2744.

Green, P., 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 84 (4), 711-732. http://dx.doi.org/10.1093/BIOMET/82.4.711.

Größer, J., Okhrin, O., 2022. Copulae: An overview and recent developments. WIREs Comput. Stat. 14 (3), http://dx.doi.org/10.1002/wics.1557, ISSN: 1939-5108, 1939-0068 URL https://onlinelibrary.wiley.com/doi/10.1002/wics.1557.

Gruber, L.F., Czado, C., 2018. Bayesian model selection of regular vine copulas. Bayesian Anal. (ISSN: 1936-0975) 13 (4), http://dx.doi.org/10.1214/17-BA1089, URL https://projecteuclid.org/journals/bayesian-analysis/volume-13/issue-4/Bayesi an-Model-Selection-of-Regular-Vine-Copulas/10.1214/17-BA1089.full.
Hawkes, A.G., 2018. Hawkes processes and their applications to finance: a review. Quant. Finance 18 (2), 193-198. http://dx.doi.org/10.1080/14697688.2017.1403131, ISSN: 1469-7688, 1469-7696 URL https://www.tandfonline.com/doi/full/10.1080/14697688.2017.1403131.

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J.-N., 2018. ERA5 hourly data on single levels from 1959 to present. Copernic. Clim. Chang. Serv. (C3S) Clim. Data Store (CDS) http://dx.doi.org/10.24381/eds.adbb2d47.

Hurrell, J.W., Kushnir, Y., Ottersen, G., Visbeck, M., 2003. An overview of the North Atlantic Oscillation. In: The North Atlantic Oscillation: Climatic Significance and Environmental Impact. American Geophysical Union (AGU), ISBN: 978-1-118-66903-7, pp. 1-35. http://dx.doi.org/10.1029/134GM01, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/134GM01, _eprint https://onlinelibrary.wiley.com/doi/pdf/10.1029/134GM01.

Joe, H., 2015. Dependence modeling with copulas. In: Monographs on statistics and applied probability, (No. 134), CRC Press, Taylor & Francis Group, Boca Raton, ISBN: 978-1-4665-8322-1.

Lang, M., Ouarda, T., Bobée, B., 1999. Towards operational guidelines for over-threshold modeling. J. Hydrol. (ISSN: 00221694) 225 (3-4), 103-117. http://dx.doi.org/10.1016/S0022-1694(99)00167-5, URL https://linkinghub.elsevier.com/retrieve/pii/S0022169499001675.

Langousis, A., Veneziano, D., 2007. Intensity-duration-frequency curves from scaling representations of rainfall. Water Resour. Res. 43 (2), 2006WR005245. http://dx.doi.org/10.1029/2006WR005245, ISSN: 0043-1397, 1944-7973 URL https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2006WR005245.

Lee, T., Ouarda, T.B.M.J., Seidou, O., 2023. Characterizing and forecasting climate indices using time series models. Theor. Appl. Climatol. 152 (1-2), 455-471. http://dx.doi.org/10.1007/s00704-023-04434-z, ISSN: 0177-798X, 1434-4483 URL https://link.springer.com/10.1007/s00704-023-04434-z.

Lee, T., Shin, J.-Y., Kim, J.-S., Singh, V.P., 2020. Stochastic simulation on reproducing long-term memory of hydroclimatological variables using deep learning model. J. Hydrol. (ISSN: 00221694) 582, 124540. http://dx.doi.org/10.1016/j.jhydrol.2019.124540, URL https://linkinghub.elsevier.com/retrieve/pii/S0022169419312752.

Li, X., Genest, C., Jalbert, J., 2021. A self-exciting marked point process model for drought analysis. Environmetrics (ISSN: 1099-095X) 32 (8), e2697. http://dx.doi.org/10.1002/env.2697, URL http://onlinelibrary.wiley.com/doi/abs/10.1002/env.2697.

Massei, N., Dieppois, B., Hannah, D.M., Lavers, D.A., Fossa, M., Laignel, B., Debret, M., 2017. Multi-time-scale hydroclimate dynamics of a regional watershed and links to large-scale atmospheric circulation: Application to the Seine river catchment, France. J. Hydrol. 546, 262-275. http://dx.doi.org/10.1016/j.jhydrol.2017.01.008.

McDonald, J.B., 1984. Some generalized functions for the size distribution of income. Econometrica (ISSN: 00129682) 52 (3), 647. http://dx.doi.org/10.2307/1913469, URL https://www.jstor.org/stable/1913469?origin=crossref.

Min, A., Czado, C., 2010. Bayesian inference for multivariate copulas using pair-copula constructions. J. Financ. Econ. 8 (4), 511-546. http://dx.doi.org/10.1093/jfinec/nbp031.

Morales Napoles, O., 2009. Bayesian Belief Nets and Vines in Aviation Safety and Other Applications (Ph.D. thesis). Technische Universität Delft.

Nadarajah, S., Kotz, S., 2006. R programs for computing truncated distributions. J. Stat. Softw. (ISSN: 1548-7660) 16, http://dx.doi.org/10.18637/jss.v016.c02, URL http://www.jstatsoft.org/v16/c02/.

Nagler, T., Vatter, T., 2023. rvinecopulib: High performance algorithms for vine copula modeling. URL https://CRAN.R-project.org/package=rvinecopulib.

Naveau, P., Huser, R., Ribereau, P., Hannart, A., 2016. Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection. Water Resour. Res. (ISSN: 1944-7973) 52 (4), 2753-2769. http://dx.doi.org/10.1002/2015WR018552, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/2015WR018552_eprint: http://dx.doi.org/10.1002/2015WR018552.

Ouarda, T.B.M.J., Charron, C., Mahdi, S., Yousef, L.A., 2021. Climate teleconnections: interannual variability, and evolution of the rainfall regime in a tropical Caribbean island: case study of Barbados. Theor. Appl. Climatol. 145 (1-2), 619-638. http://dx.doi.org/10.1007/s00704-021-03653-6, ISSN: 0177-798X, 1434-4483 URL https://link.springer.com/10.1007/s00704-021-03653-6.

Ouarda, T.B.M.J., Yousef, L.A., Charron, C., 2018. Non-stationary intensity-duration-frequency curves integrating information concerning teleconnections and climate change. Int. J. Climatol. 39 (4), 2306-2323. http://dx.doi.org/10.1002/joc.5953.

Pan, X., Yildirim, G., Rahman, A., Haddad, K., Ouarda, T.B.M.J., 2023. Peaks-over-threshold-based regional flood frequency analysis using regularised linear models. Water (ISSN: 2073-4441) 15 (21), 3808. http://dx.doi.org/10.3390/w15213808, URL https://www.mdpi.com/2073-4441/15/21/3808.

Papalexiou, S.M., Sinalaldi, F., Clark, M.P., 2023. Large-domain multisite precipitation generation: Operational blueprint and demonstration for 1,000 sites. Water Resour. Res. 59 (3), e2022WR034094. http://dx.doi.org/10.1029/2022WR034094, ISSN: 0043-1397, 1944-7973 URL https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022WR034094.

Papastathopoulos, I., Tawn, J.A., 2013. Extended generalised Pareto models for tail estimation. J. Statist. Plann. Inference (ISSN: 0378-3758) 143 (1), 131-143. http://dx.doi.org/10.1016/j.jspi.2012.07.001, URL https://www.sciencedirect.com/science/article/pii/S0378375812002388.

Sinalaldi, F., Kilsby, C.G., 2014. Rainfall extremes: Toward reconciliation after the battle of distributions. Water Resour. Res. 50 (1), 336-352. http://dx.doi.org/10.1002/2013WR014211, ISSN: 0043-1397, 1944-7973 URL https://onlinelibrary.wiley.com/doi/10.1002/2013WR014211.

Smith, R., 2003. Statistics of extremes, with applications in environment, insurance, and finance. In: Extreme Values in Finance, Telecommunications, and the Environment, first ed. ISBN: 978-0-429-20451-7, p. 78, URL https://wwwsemanticscholar.org/paper/Statistics-of-Extremes%2C-with-Applications-in-and-Smith/9fe6a18337ec5ee2a520e73795e2db1ac46f343f.

Sun, X., Renard, B., Thyer, M., Westra, S., Lang, M., 2015. A global analysis of the asymmetric effect of ENSO on extreme precipitation. J. Hydrol. (ISSN: 00221694) 530, 51-65. http://dx.doi.org/10.1016/j.jhydrol.2015.09.016, URL https://linkinghub.elsevier.com/retrieve/pii/S002216941500699X.

Ulrich, J., Jurado, O.E., Peter, M., Scheibel, M., Rust, H.W., 2020. Estimating IDF curves consistently over durations with spatial covariates. Water 12 (11), http://dx.doi.org/10.3390/w12113119.

Ward, P.J., De Ruiter, M.C., Mård, J., Schröter, K., Van Loon, A., Veldkamp, T., Von Uexkull, N., Wanders, N., AghaKouchak, A., Arnbjerg-Nielsen, K., Capewell, L., Carmen Llasat, M., Day, R., Dewals, B., Di Baldassarre, G., Huning, L.S., Kreibich, H., Mazzoleni, M., Savelli, E., Teutschbein, C., Van Den Berg, H., Van Der Heijden, A., Vincken, J.M., Waterloo, M.J., Wens, M., 2020. The need to integrate flood and drought disaster risk reduction strategies. Water Secur. (ISSN: 24683124) 11, 100070. http://dx.doi.org/10.1016/j.wasec.2020.100070, URL https://linkinghub.elsevier.com/retrieve/pii/S2468312420300109.