Dealing with seasonal variability and inter-site dependence in regional flood frequency analysis

TL;DR

This paper introduces two novel regional estimators for high river flow quantiles that incorporate seasonal variability and spatial dependence, improving flood risk assessment accuracy.

Contribution

It proposes parametric and semi-parametric estimators that handle seasonal and spatial features, with proven asymptotic normality and no reliance on complex dependence models.

Findings

Estimators effectively incorporate seasonal variability.

Asymptotic normality enables uncertainty quantification.

Application to German river data demonstrates practical utility.

Abstract

This paper considers the regional estimation of high quantiles of annual maximal river flow distributions $F$, an important problem from flood frequency analysis. Even though this particular problem has been addressed by many papers, less attention has been payed to incorporating seasonal variability and spatial dependence into the methods. We are going to discuss two regional estimators of high quantiles of local distributions $F$ that are able to deal with these important features, namely, a parametric approach based on so-called two-component extreme value distributions and a semi-parametric approach based on regional estimation of a tail index. The asymptotic normality of the estimators is derived for both procedures, which for instance enables us to account for estimation uncertainty without the need of parametric dependence models or bootstrap procedures. A comprehensive simulation study is conducted and our main findings are illustrated on river flow series from the Mulde basin in Germany, where people have suffered several times from severe floods over the last 100 years.

Equations

$${\mathcal{H}}_{0,\text{IF}}:F_j^{-1}=s_j\cdot G_{\theta }^{-1}\forall j=1,\dots ,d,$$

$$G_{\theta }(x)=\exp \left[-{\left(1+\xi \frac{x-\mu }{\sigma }\right)}^{-1/\xi }\right]\text{ for }1+\xi \frac{x-\mu }{\sigma }>0,$$

$$F(x)=\mathbb{P}(\max \{W,S\}\le x)=\mathbb{P}(W\le x,S\le x)=F_w(x)\cdot F_s(x)\approx G_{{\theta }_w}(x)\cdot G_{{\theta }_s}(x)$$

$$\displaystyle\left.\begin{array}[]{r}X_{1,1},\ X_{2,1},\ X_{3,1},\ X_{4,1},\ X_{5,1},\ \ldots,\ X_{n,1}~{}\sim F_{1}\\ X_{a_{2}+1,2},\ X_{a_{2}+2,2},\ X_{a_{2}+3,2},\ \ldots,\ X_{n,2}~{}\sim F_{2}\\ \ddots\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\vdots\\ X_{a_{d}+1,d},\ X_{a_{d}+2,d},\ \ldots,\ X_{n,d}~{}\sim F_{d}\end{array}\right\}~{}\begin{array}[]{c}\text{observations}\\ \text{from}~{}d~{}\text{sites}\end{array}$$

$$F(x)=\mathbb{P}(\max \{W,S\}\le x)=\mathbb{P}(W\le x,S\le x)=\mathbb{P}(W\le x)\cdot \mathbb{P}(S\le x)$$

$${\theta }^{\ast }=({\mu }^{\ast },{\sigma }^{\ast },{\xi }^{\ast })=\underset{\mu ,\xi \in \mathbb{R},\sigma >0}{\mathrm{argmin}}{d}_{KL}(F,G_{\theta })=(2.554,1.235,0.305)$$

$$q_p=F^{-1}(p)\text{ by }{\hat{q}}_p={\hat{q}}_{p,n}={\left(G_{{\hat{\theta }}_w}\cdot G_{{\hat{\theta }}_s}\right)}^{-1}(p),$$

$$\displaystyle\Sigma_{w}=\lim_{n\rightarrow\infty}{\rm Var}\left[\sqrt{n}\big{(}\hat{\theta}_{w}-\theta_{w}\big{)}\right]~{}\text{ and }~{}\Sigma_{s}=\lim_{n\rightarrow\infty}{\rm Var}\left[\sqrt{n}\big{(}\hat{\theta}_{s}-\theta_{s}\big{)}\right].$$

$$\displaystyle\sigma^{2}_{p}=\lim_{n\rightarrow\infty}{\rm Var}\left[\sqrt{n}\big{(}\hat{q}_{p}-q_{p}\big{)}\right]=\frac{G^{2}_{\theta_{s}}(q_{p})\cdot J_{\theta_{w}}(q_{p})\Sigma_{w}J_{\theta_{w}}(q_{p})^{\prime}+G^{2}_{\theta_{w}}(q_{p})\cdot J_{\theta_{s}}(q_{p})\Sigma_{s}J_{\theta_{s}}(q_{p})^{\prime}}{\left[g_{\theta_{w}}(q_{p})\cdot G_{\theta_{s}}(q_{p})+G_{\theta_{w}}(q_{p})\cdot g_{\theta_{s}}(q_{p})\right]^{2}},$$

$$I_{1-\alpha }=\left[{\hat{q}}_p-z_{1-\alpha /2}\cdot \frac{{\hat{\sigma }}_p}{\sqrt{n}},{\hat{q}}_p+z_{1-\alpha /2}\cdot \frac{{\hat{\sigma }}_p}{\sqrt{n}}\right]$$

$$\bar{F}(x)=1-F(x)=x^{-1/\gamma }\cdot L(x),x>0,$$

$$\underset{t\to \infty }{\lim }\frac{U(tx)}{U(t)}=x^{\gamma }\text{ for all }x>0.$$

$$F_j(x)=\mathbb{P}(X_{i,j}\le x)=1-x^{-1/{\gamma }_j}{L}_j(x),x>0,j=1,\dots ,d,$$

$${\mathcal{H}}_{0,\text{evi}}:{\gamma }_1=\dots ={\gamma }_d=\gamma \text{ for some }\gamma >0.$$

$${\hat{F}}_1^{-1}(p)={\hat{u}}_1\cdot {\left(\frac{k_1}{n_1(1-p)}\right)}^{\hat{\gamma }(\mathbf{w})}\text{ with }\hat{\gamma }(\mathbf{w})=\sum _{j=1}^d{w}_j{\hat{\gamma }}_j,$$

$${\hat{F}}_1(x)=1-\frac{k_1}{n_1}\cdot {\left(\frac{x}{{\hat{u}}_1}\right)}^{-1/\hat{\gamma }(\mathbf{w})}.$$

$$\mathrm{Var}\left[\hat{\gamma }(\mathbf{w})\right]\approx \frac{{\gamma }^2}{k_1}{\mathbf{w}}^{\prime }\Sigma \mathbf{w}$$

$${\mathbf{w}}_{\text{opt}}=\underset{\mathbf{w}\in \mathbf{W}}{\mathrm{argmin}}\frac{{\gamma }^2}{k_1}{\mathbf{w}}^{\prime }\Sigma \mathbf{w}={\left(1^{\prime }{\Sigma }^{-1}1\right)}^{-1}\cdot {\Sigma }^{-1}1,$$

$$C_{\mathbf{\theta },\mathbf{c}}=C_{{\theta }_1}({\mathbf{u}}^{\mathbf{c}})\cdot C_{{\theta }_2}({\mathbf{u}}^{1-\mathbf{c}}),\mathbf{\theta }=({\theta }_1,{\theta }_2)\in [1,\infty {)}^2,\mathbf{c}\in [0,1{]}^d,$$

$$F_j(x)=F_j^{(b)}(x)={\left[2\cdot T_{1/\xi }\left(\left\{1+\xi \frac{x-{\mu }_j}{{\sigma }_j}\right\}\cdot T_{1/\xi }^{-1}\left(1-\frac{1}{2b}\right)\right)-1\right]}^b,$$

$${\mu }_j+\frac{{\sigma }_j}{\xi }\left(\frac{\max \{|Z_1|,\dots ,|Z_b|\}}{a_b}-1\right)\text{with }{a}_b=T_{1/\xi }^{-1}\left(1-\frac{1}{2b}\right),$$

$${\beta }_k={\beta }_k(\theta )={\int }_{\mathbb{R}}x\cdot G_{\theta }^k(x)d{G}_{\theta }(x)={\int }_{(0,1)}{G}_{\theta }^{-1}(u)\cdot u^{k}du.$$

$$\displaystyle\left\{\begin{array}[]{cll}\frac{3^{\xi}-1}{2^{\xi}-1}&=\frac{3\beta_{2}-\beta_{0}}{2\beta_{1}-\beta_{0}}&=\frac{\lambda_{3}}{2\lambda_{2}}+\frac{3}{2}\\ \sigma&=\frac{(2\beta_{1}-\beta_{0})\xi}{\Gamma(1-\xi)(2^{\xi}-1)}&=\frac{\lambda_{2}\xi}{\Gamma(1-\xi)(2^{\xi}-1)}\\ \mu&=\beta_{0}+\frac{\sigma}{\xi}\left[1-\Gamma(1-\xi)\right]&=\lambda_{1}+\frac{\sigma}{\xi}\left[1-\Gamma(1-\xi)\right]\end{array}\right.,$$

$$\xi =-7.859\cdot h(\mathbf{\beta })-2.9554\cdot h(\mathbf{\beta }{)}^2\text{ with }h(\mathbf{\beta })=\frac{2{\beta }_1-{\beta }_0}{3{\beta }_2-{\beta }_0}-\frac{\log 2}{\log 3}=\frac{2}{3+{\lambda }_3/{\lambda }_2}-\frac{\log 2}{\log 3}$$

$$\displaystyle\left\{\begin{array}[]{cll}\frac{5\cdot 4^{\xi}-12\cdot 3^{\xi}+9\cdot 2^{\xi}-2}{3^{\xi}-2^{\xi+1}+1}&=\frac{2\cdot(18\beta_{2}-9\beta_{1}+\beta_{0}-10\beta_{3})}{4\beta_{1}-\beta_{0}-3\beta_{2}}&=\frac{9\lambda_{3}^{(0,1)}}{4\lambda_{2}^{(0,1)}}\\ \sigma&=\frac{4\beta_{1}-\beta_{0}-3\beta_{2}}{\Gamma(-\xi)\cdot(3^{\xi}-2^{\xi+1}+1)}&=\frac{2/3\cdot\lambda_{2}^{(0,1)}}{\Gamma(-\xi)\cdot(3^{\xi}-2^{\xi+1}+1)}\\ \mu&=2(\beta_{0}-\beta_{1})+\frac{\sigma}{\xi}-\frac{\sigma\cdot\Gamma(-\xi)}{(2^{\xi}-2)^{-1}}&=\lambda_{1}^{(0,1)}+\frac{\sigma}{\xi}-\frac{\sigma\cdot\Gamma(-\xi)}{(2^{\xi}-2)^{-1}}\\ \end{array}\right.$$

$$\xi =-8.567394\cdot h(\mathbf{\beta })+0.675969\cdot h(\mathbf{\beta }{)}^2\text{ with }h(\mathbf{\beta })=\frac{6{\beta }_1-3/2{\beta }_0-9/2{\beta }_2}{15{\beta }_2-5/3{\beta }_0-40/3{\beta }_3}-\frac{2\log 2-\log 3}{3\log 3-2\log 4}$$

$${\hat{\beta }}_{k,j,r_j,n}=\frac{1}{n_j}\sum _{i=1}^{n_j}{X}_{a_j+i}\cdot F_{j,a_j+1:n}^k(X_{a_j+i}),k=0,\dots ,K,j=1,\dots ,d$$

$$\mathbf{\beta }={\left({\beta }_0(F_1),\dots ,{\beta }_K(F_1),{\beta }_0(F_2),\dots ,{\beta }_K(F_2),\dots \dots ,{\beta }_0(F_d),\dots ,{\beta }_K(F_d)\right)}^{\prime }\in {\mathbb{R}}^{d\cdot K}.$$

$$\mathbb{E}\left[X_j{F}_j^k(X_j)X_{\ell }{F}_{\ell }^m(X_{\ell })\right]<\infty \text{ for all }1\le j,\ell \le d\text{ and }0\le k,m<K$$

$$\sqrt{n}\left({\hat{\mathbf{\beta }}}_{\mathbf{r},n}-\mathbf{\beta }\right)\stackrel{D}{⟶}\mathcal{N}\left(0,{\Sigma }_{\mathbf{r}}\right),$$

Full text

Dealing with seasonal variability and inter-site dependence in regional flood frequency analysis

Paul Kinsvater 111Corresponding author. Technische Universität Dortmund, Fakultät Statistik. E-mail: kinsvater@statistik.tu-dortmund.de. ,  Friederike Deiters   and Roland Fried

Abstract

This paper considers the regional estimation of high quantiles of annual maximal river flow distributions $F$, an important problem from flood frequency analysis. Even though this particular problem has been addressed by many papers, less attention has been payed to incorporating seasonal variability and spatial dependence into the methods. We are going to discuss two regional estimators of high quantiles of local distributions $F$ that are able to deal with these important features, namely, a parametric approach based on so-called two-component extreme value distributions and a semi-parametric approach based on regional estimation of a tail index. The asymptotic normality of the estimators is derived for both procedures, which for instance enables us to account for estimation uncertainty without the need of parametric dependence models or bootstrap procedures. A comprehensive simulation study is conducted and our main findings are illustrated on river flow series from the Mulde basin in Germany, where people have suffered several times from severe floods over the last 100 years.

1 Introduction

Flood frequency analysis (FFA) deals with the estimation of river flow distributions. A flow value is the amount of water (in cubic meter per second, m3/s) passing a measurement station. The ultimate objective is to determine the design of future flood protection systems, for example, the height of a dam for some predefined non-failure probability. Common rules found in many official guidelines (e.g., DWA, 2012) focus on the distribution $F(x)=\mathbb{P}(X\leq x)$ of annual maximal flows $X$ at the site of interest and determine the height of a dam according to a high quantile $q=F^{-1}(p)$, for a given $p\in(0,1)$. Often experts have to deal with rather high probabilities $p\geq 0.99$, depending on the safety-relevance of the local environment.

In practice the distribution $F$ can be estimated from annual maximal flows $X_{1},\ldots,X_{n}$ of the past $n$ years. A serious problem in FFA is that only small sample lengths $n$ are available at each site. For a potential reduction of the estimation uncertainty, so-called regional methods are applied. These procedures combine observations from a set of sites $j\in\{1,2,\ldots,d\}$ sharing similar site characteristics related to the flood magnitude. For instance, the Index Flood method proposed in the seminal work by Dalrymple (1960) is based on the hypothesis that the quantile functions are all identical up to an unknown local scale $s_{j}=s(F_{j})>0$,

[TABLE]

with distribution $G_{\theta}$ known up to a finite-dimensional parameter $\theta$. Thanks to the extreme value theorem (Fisher and Tippett, 1928), a recommended choice for the parametric model $G_{\theta}$ is a generalized extreme value (GEV) distribution function

[TABLE]

with parameters $\theta=(\mu,\sigma,\xi)^{\prime}\in\mathbb{R}\times\mathbb{R}_{+}\times\mathbb{R}$ called location, scale and shape. These distributions arise as the only possible non-degenerate limit laws of standardized block maxima $X^{(b)}=\max\{Z_{1},\ldots,Z_{b}\}$ over independent and identically distributed variables $Z_{i}$ as the block size $b$ tends to infinity.

However, it is also important to account for inter-site dependence. Besides the fact that stronger positive dependence between stations reduces the accuracy of regional estimators (Stedinger, 1983), we need to estimate the correlation between local statistics to be able to consistently estimate the precision of regional approaches (Martins and Stedinger, 2002; Cunderlik and Burn, 2006). Contrary to the previous references, we do not apply Monte Carlo or Bootstrap procedures. Instead we use direct estimates of correlation based on recent theory developed by Lilienthal et al. (2016).

Treating an annual maximal flow $X=\max\{Z_{\text{Jan}},Z_{\text{Feb}},\ldots,Z_{\text{Dec}}\}$ as a maximum over twelve independent and identically distributed monthly maximal flows $Z_{\text{month}}$ motivates the assumption that $F$ equals $G_{\theta}$ approximately for some unknown $\theta=(\mu,\sigma,\xi)^{\prime}$ due to the finite block length of $b=12$. We even might assume that $F=G_{\theta}$ holds, which allows to estimate $F$ by parametric methods. A theoretical justification of this simplification is studied by Ferreira and de Haan (2015) for estimation based on $L$-moments and by Bücher and Segers (2015) for maximum likelihood estimation, where the asymptotic normality of the methods is derived for block maxima instead of exact GEV distributions (with block size tending to infinity).

The former simplification of identically distributed monthly maximal flows $Z_{\text{Jan}},\ldots,Z_{\text{Dec}}$ is not realistic. River flows, similar to temperature and rainfall, are subject to seasonal variability. Flows in the winter/spring season are fed by large masses of melting snow, while floods in the summer/autumn period usually are caused by short but heavy rainfalls. For empirical evidence of the seasonal variability we refer to Figure 1. The boxplots illustrate the distribution of monthly maximal flows for each specific month at station Lichtenwalde. The plot is computed from observations during the years 1910–2011. A seasonal pattern in location (median) and dispersion (interquartile range) of the boxes becomes apparent. This confirms that different physical mechanisms (melting snow, heavy rainfall) cause different river flow distributions. Summing up, it is not clear whether the assumption $F\approx G_{\theta}$ is justified, at least approximately, for annual maximal flow distributions.

An approach that is particularly designed to account for the different flood generating mechanisms proceeds as follows: Suppose that our observations include seasonal maxima $W$ and $S$ from two disjoint seasons, e.g. hydrological winter and summer, such that $X=\max\{W,S\}$. Both, $W$ and $S$, represent maximal values that are computed from more homogeneous periods of, say, six monthly maximal observations with distribution functions $\mathbb{P}(W\leq x)=F_{w}(x)$ and $\mathbb{P}(S\leq x)=F_{s}(x)$. We thus may assume that $F_{w}\approx G_{\theta_{w}}$ and $F_{s}\approx G_{\theta_{s}}$ for some parameter vectors $\theta_{w}$ and $\theta_{s}$. Furthermore, assuming independence between $W$ and $S$, we obtain

[TABLE]

for the distribution function $F$ of annual maximal flows. Again, we may assume $F(x)=G_{\theta_{w}}(x)\cdot G_{\theta_{s}}(x)$ for practical reasons. A distribution function $F=G_{\theta_{w}}\cdot G_{\theta_{s}}$ is called two-component GEV (Rossi et al., 1984; Gabriele and Arnell, 1991; Strupczewski et al., 2012; Rulfová et al., 2016).

In the mentioned references it is stated that such a seasonal approach potentially leads to more plausible estimates of high quantiles of $F$, where plausibility may be interpreted as a smaller number of observations being classified as outliers (Rossi et al., 1984). On the other hand, seasonal GEV estimation means that maxima over blocks with a smaller size are considered, say, $b=6$ instead of $b=12$. Is it still plausible to assume that $F_{w}\approx G_{\theta_{w}}$ and $F_{s}\approx G_{\theta_{s}}$ are GEV distributions, even though the block size $b$ is small and still there might be some minor variability left within each season? If not, what kind of model should we choose instead?

As a compromise between parametric and purely non-parametric approaches we also present a regional estimator of high quantiles under semi-parametric model assumptions based on estimation of a regional extreme value index (Kinsvater et al., 2016). For that purpose we are taking things a step backward: We avoid the GEV assumption by ignoring the fact that our observations $X$ are maximal values and simply treat them as realizations from any heavy-tailed distribution $F$. Instead of fitting the whole distribution $F$, we only focus on the estimation of an approximately parametric right tail. Similar to the two-component GEV model, this semi-parametric procedure leads us to more plausible estimates of high quantiles, since solely the largest observations from $F$ are taken effectively into account.

Our main results can be summarized as follows:

Theory: We discuss two new approaches for joint estimation of high quantiles of annual maximal flow distributions from a homogeneous group of $d$ stations: An estimator based on the two-component GEV model (Rulfová et al., 2016) and one based on semi-parametric assumptions. We derive asymptotic theory for both approaches under a very general framework (inter-site dependence, site-varying sample lengths) in order to be able to evaluate the estimation uncertainty of the procedures. This, for instance, leads to asymptotically optimal weights in the computation of regional estimates. 2. 2.

Simulation: Accounting for seasonality improves the estimation efficiency of high quantile estimators. If seasonal maximal flows, say, $W$ and $S$, are available, we recommend the parametric approach based on two-component distributions. Otherwise, if only annual data $X$ is available, the semi-parametric estimator still is rather competitive, provided the dimension $d$ of the homogeneous group is sufficiently large $(d\geq 5)$. 3. 3.

Application: We consider river flows from a homogeneous group of $d=8$ stations from the Mulde river basin, where winter and summer flows have considerably different tails. Our two-component GEV and semi-parametric procedures are able to deal with the seasonal variability, leading to more plausible estimates of high quantiles.

The severe flood in the summer of 2002 plays an outstanding role in our data set. The corresponding empirical return period/level combination is the only point in Figure 6 that is far away from the estimated curves. A classical GEV approach suggests that an event of at least the strength of 2002 occurs every 400 years on average. On the other hand, our two-component GEV approach suggests a return period of around 200 years, implying much higher risk of severe floods.

The rest of the paper is organized as follows. We begin with a short description of the data in Section 2. The next two sections introduce the models and related statistical methods, starting with two-component GEV distributions in Section 3 and followed by semi-parametric inference in Section 4. Section 5 reports a simulation study and in Section 6 we apply the methods to maximal flows from several river stations located at the Mulde basin in Germany. Technical details and proofs are deferred to Appendix A and B.

2 Data

Regional frequency analysis (RFA) combines observations from many sites within a region, even though the primary interest is still the estimation of quantiles at one target site.

In our application, see Section 6, we will have a look at a dataset which contains monthly maximal flows (in m3/s) from $d=8$ stations located in the Mulde river basin (Saxony, Germany). The sample lengths of these sites range from 83 to 102 hydrological years (in Germany from November to October). The two seasons arise by splitting the hydrological year into the hydrological winter from November to April and summer from May to October; so each season contains 6 months. These two seasons are defined according to official guidelines (DWA, 2012). The winter and summer maximal flows can be regarded as independent for any specific site, but simultaneous observations at sites from the same area are dependent. Within a site our data is complete, so we do not have to deal with missing values here. The structure of our data is described as follows:

Suppose that the random vector $\mathbf{X}=(X_{1},\ldots,X_{d})^{\prime}$ represents annual maximal river flows at the $d$ stations, for any generic year. We denote the local distribution of site $j$ by $F_{j}(x)=\mathbb{P}(X_{j}\leq x)$. Unlike in many other regional approaches, we will not consider the components as independent. Let $\mathbf{X}_{i}=(X_{i,1},\ldots,X_{i,d})^{\prime}$, with time index $i=1,\ldots,n$, be independent copies of $\mathbf{X}$, where $\{1,\ldots,n\}$ covers the observation period. In regional settings it is unlikely that the observation period is the same for all $d$ sites. Often the beginning of recordings is different for different sites because of the different ages of the measurement stations. If we arrange these $d$ sites according to their local sample lengths $n_{j}$, $j=1,\ldots,d$, such that $n=n_{1}\geq n_{2}\geq\ldots\geq n_{d}$, we obtain the observation scheme

[TABLE]

where $a_{j}\in\{0,\ldots,n-1\}$ with $0=a_{1}\leq a_{2}\leq\ldots\leq a_{d}$ and the local sample lengths $n_{j}=n-a_{j}$ for $j=1,\ldots,d$. Each row $X_{i,j}$, $i=a_{j}+1,\ldots,n,$ contains solely observations of one site $j\in\{1,\ldots,d\}$. It is important to consider scheme (10) in order to be able to take into account the dependence between local estimates computed from the scheme.

In fact, each of the observations from the previous scheme can be considered as a maximum $X_{i,j}=\max\{W_{i,j},S_{i,j}\}$ of corresponding winter and summer maximal flows, which is also relevant for the data set analyzed in Section 6. That means, we have a scheme as in (10) for winter and summer maximal flows and combine them to the scheme of the annual maximal flows in (10).

3 Parametric estimation under seasonal variability

In this section we assume that the available data consists of winter and summer maximal flows $W$ and $S$, with $W$ and $S$ being stochastically independent. Our focus is on estimation of high quantiles $q_{p}=F^{-1}(p)$ of the distribution $F(x)=\mathbb{P}(X\leq x)$ of an annual maximum $X=\max\{W,S\}$. By independence of $W$ and $S$, we obtain a factorization

[TABLE]

and, because both $W$ and $S$ are maxima, we will also assume that these two variables are GEV distributed. We set $\mathbb{P}(W\leq x)=G_{\theta_{w}}(x)$ and $\mathbb{P}(S\leq x)=G_{\theta_{s}}(x)$ for some unknown seasonal parameters $\theta_{w},\theta_{s}$ and with $G_{\theta}$ defined in (2). A distribution function $F=G_{\theta_{w}}\cdot G_{\theta_{s}}$ is called two-component GEV.

3.1 Some remarks on two-component GEVs

The parameters $\theta=(\mu,\sigma,\xi)$ of the GEV model $G_{\theta}$ defined in (2) represent location, scale and shape, respectively. These distributions are unimodal and the shape $\xi\in\mathbb{R}$ determines the right tail behavior of $G_{\theta}$. The right tail is bounded, of exponential or of polynomial order for $\xi<0$, $\xi=0$ or $\xi>0$, respectively. The heaviness of the right tail, which may be considered as the risk of extraordinary large realizations from $G_{\theta}$, increases with increasing $\xi$. In extreme value theory this shape is also called extreme value index (EVI), denoted by $\gamma=\xi$. More generally, an EVI $\gamma$ can be determined for any distribution $F$ in the maximum domain of attraction. For instance, the normal distribution with its exponential tail behavior has an EVI of $\gamma=0$. Students $t_{\nu}$ and Fishers $F_{m,k}$ have a polynomial right tail with EVI $\gamma=1/\nu>0$ and $\gamma=2/k>0$, respectively. A more comprehensive list of parametric models and corresponding EVIs can be found in Beirlant et al. (2006, Chap. 2).

In case of positive shapes $\xi_{w},\xi_{s}>0$, which is likely to hold in our data application, the EVI of $F=G_{\theta_{w}}\cdot G_{\theta_{s}}$ is given by $\gamma=\max\{\xi_{w},\xi_{s}\}$, i.e., the right tail behavior of $F$ is determined by that of the heaviest component (Kinsvater, 2016, Cor. 1.11).

To fix ideas, suppose that our annual maximum $X$ follows a two-component GEV $F=G_{\theta_{w}}\cdot G_{\theta_{s}}$ with seasonal parameters $\theta_{w}=(\mu_{w},\sigma_{w},\xi_{w})=(2,1,0.2)$ and $\theta_{s}=(\mu_{s},\sigma_{s},\xi_{s})=(1.5,1,0.4)$. The corresponding density is depicted in Figure 2. In practice we do not know the true model, so it might be preferable to forget seasonal variability and use instead a usual one-component GEV model $G_{\theta}$. Interestingly enough, by computing the minimizer

[TABLE]

of the Kullback-Leibler distance $d_{KL}$ between $F$ and $G_{\theta}$, we obtain a GEV depicted in Figure 2 that looks almost identical to the density of $F$. Since we are mainly interested in high quantiles, we also compare these two models by their return level curves in Figure 3. From this second point of view we can clearly identify that these two models deviate reasonably from each other in the relevant tail region. The GEV model $G_{\theta^{*}}$ underestimates the true EVI of $\gamma=\max\{0.2,0.4\}=0.4$ by $\xi^{*}=0.305$, which results in underestimating risk of extreme realizations. A similar conclusion is drawn in our simulation study in Section 5.

3.2 Accuracy of quantile estimates

Let $F=G_{\theta_{w}}\cdot G_{\theta_{s}}$ be a two-component GEV distribution, where $\theta_{w},\theta_{s}$ are two unknown parameter vectors and where $G_{\theta}$ is defined in (2). For fixed $p\in(0,1)$ we consider the problem of estimating

[TABLE]

where $\hat{\theta}_{w}$ and $\hat{\theta}_{s}$ are estimators of $\theta_{w}$ and $\theta_{s}$, respectively. It should be noted that the equation $q=\varphi(p)=\left(G_{\theta_{w}}\cdot G_{\theta_{s}}\right)^{-1}(p)$ does not have an explicit solution, except in case of $\theta_{w}=\theta_{s}$. Since $\varphi$ is monotone increasing, it is easy to obtain an accurate numerical solution. However, a non-trivial difficulty arises when one is not only interested in a point-estimate $\hat{q}_{p}$ but also in its precision or in a confidence interval.

Suppose that $\hat{\theta}_{w}=\hat{\theta}_{w,n}$ and $\hat{\theta}_{s}=\hat{\theta}_{s,n}$ are two independent, asymptotically normal and consistent estimators with asymptotic covariance matrices

[TABLE]

In our application we consider the regional TL-moment estimators from Appendix B.1, with $\hat{\theta}_{w}$ computed from winter and $\hat{\theta}_{s}$ from summer maximal flows. The corresponding theory, the asymptotic normality of the regional TL-moment approaches, has been proved in a recent paper by Lilienthal et al. (2016).

In this paper we prove that $\hat{q}_{p}$ is an asymptotically normal and consistent estimator of $q_{p}$ (see Appendix B.2) with limiting variance

[TABLE]

where $g_{\theta}$ is the density of $G_{\theta}$ and $J_{\theta}(x)=\frac{\partial}{\partial\theta}G_{\theta}(x)\in\mathbb{R}^{1\times 3}$ is the Jacobi matrix of $G_{\theta}(x)$ in $\theta$.

In practice all unknown quantities on the right-hand side of (12) can be replaced by their sample counterparts in order to obtain an estimate $\hat{\sigma}^{2}_{p}$ of $\sigma^{2}_{p}$. This, in turn, is used to calculate an asymptotically valid $(1-\alpha)$-confidence interval

[TABLE]

for $\hat{q}_{p}$ and $\alpha\in(0,1)$, where $z_{1-\alpha/2}$ denotes the $(1-\alpha/2)$-quantile of $\mathcal{N}(0,1)$.

4 Semi-parametric estimation of heavy tails

Contrary to the previous section, we only require annual observations $X$ for the computation of estimates in this section. In addition, we will not assume that the annual distribution $F$ of $X$ is determined by a finite number of parameters. Instead we apply a semi-parametric framework under a certain heavy-tail assumption: A distribution function $F$ with right-unlimited support is called a Pareto-type distribution function with extreme value index $\gamma>0$, if

[TABLE]

for some (measurable) function $L:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ satisfying $L(tx)/L(t)\rightarrow 1$ for $t\rightarrow\infty$ and all $x>0$. The characterization in (14) is called semi-parametric, since the parameter $\gamma>0$ is of core interest, while the distribution cannot be characterized by a finite-dimensional parameter vector due to the presence of the function $L$. Typical textbook examples of $L$ are $L(x)=\log(1+x)$ or simply a constant function $L\equiv c>0$. Roughly speaking, slowly varying functions can be considered as almost constant $L(x)\approx c$ for values $x>u$ above some large threshold $u$, which implies that the right tail of $F$ behaves like a power function, $\bar{F}(x)\approx c\cdot x^{-1/\gamma}$ for $x>u$. Some examples of Pareto-type distributions with extreme value index $\gamma>0$ are Student’s $t_{\nu}$ with $\gamma=1/\nu$, Fisher’s $F_{m,k}$ with $\gamma=2/k$, the $\log$-Gamma$(\lambda,\alpha)$ with $\gamma=1/\lambda$, the GEV$(\mu,\sigma,\xi)$ and the generalized Pareto GP$(\sigma,\xi)$ with $\gamma=\xi>0$.

Recall that our interest is in estimation of some high quantile $q_{p}=F^{-1}(p)$, where we often have to deal with $p\geq 1-1/n$. The Pareto-type framework allows us to estimate $F$ beyond the range of observations without relying on a parametric model: Let $U(t)=F^{-1}(1-1/t)$, $t>1$, denote the $t$-return level of $F$. It is known that Pareto-type distributions $F$ with extreme value index $\gamma>0$ satisfy

[TABLE]

We thus have $U(tx)\approx U(t)\cdot x^{\gamma}$ for large $t$, which means that a high quantile $q_{p}=U(tx)$ is approximately determined by a moderate quantile $U(t)$ and the index $\gamma$.

In practice one usually sets $t<n$, with $n$ denoting the sample length, such that $U(t)$ is estimated with satisfactory precision from a simple order statistic. The difficulty remains particularly with the estimation of $\gamma$.

4.1 Semi-parametric Index Flood

Recently, Dematteo and Clémençon (2015) and Kinsvater et al. (2016) studied the Hill estimator of $\gamma$ (Hill, 1975) under a regional tail homogeneity assumption. Their theory allows us to introduce a new approach what we call a Semi-parametric Index Flood procedure. We will use the regional estimator of $\gamma$ from the previous references and plug it in into the so-called Weissman extrapolation formula (Weissman, 1978) for semi-parametric estimation of high quantiles. For readers who are not familiar with semi-parametric heavy-tail analysis, Hill’s estimator and Weissman’s extrapolation formula, we recommend to read the introduction into these topics provided in Appendix B.3 first.

Suppose that our data consists of annual maximal river flows from $d$ different sites as in (10) and that we are interested in the estimation of, say, $F^{-1}_{1}(p)$. We assume that each marginal distribution is of Pareto-type,

[TABLE]

with not necessarily identical slowly varying functions $L_{j}$ but with the same unknown extreme value index $\gamma>0$, that is,

[TABLE]

We call assumption (17) heavy-tail homogeneity. Note that $\mathcal{H}_{0,\text{evi}}$ in combination with the Pareto-type framework (16) is much weaker than the Index Flood assumption stated in (1). We do not require the distributions $F_{j}$ to be parametric and we only need equality in EVI $\gamma$, which in the GEV framework means that we only need that the shape parameters $\xi_{1}=\ldots=\xi_{d}=\gamma>0$ are identical and positive.

Let $\hat{u}_{j}$ denote the $(n_{j}-k_{j})$-th largest observation and $\hat{\gamma}_{j}$ denote Hill’s estimator computed from the $k_{j}$ relative excesses above the threshold $u_{j}$ from sample $X_{a_{j}+1,j},X_{a_{j}+2,j},\ldots,X_{n,j}$, $k_{j}\in\{2,\ldots,n_{j}\}$. Taking advantage of assumption (17), we set

[TABLE]

where $\bm{w}=(w_{1},\ldots,w_{d})^{\prime}$ is a vector of weights summing up to 1. The left-hand side of (18), also known as Weissman’s extrapolation formula (Weissman, 1978), is motivated by (15) with $t=n_{1}/k_{1}$ and $x=k_{1}/(n_{1}(1-p))$.

In conformity with equation (1), the estimator (18) can be viewed as a semi-parametric Index Flood approach, with site-specific scales $\hat{u}_{j}$ and regional tail behavior determined by $\hat{\gamma}$. The probability of large events $x>\hat{u}_{1}$ is estimated by the inversion of (18),

[TABLE]

4.2 Selection of weights and tail samples

We recap the rules considered in Kinsvater et al. (2016) and begin with the selection of weights $\bm{w}$. For the moment, let $\mathbf{k}=(k_{1},\ldots,k_{d})^{\prime}$ be a fixed vector of integers. The joint estimator $\hat{\gamma}(\bm{w})$ on the right-hand side of (18) is asymptotically normal with

[TABLE]

and $\Sigma\in\mathbb{R}^{d\times d}$ defined in Appendix B.3.2. Using Lagrange multipliers it is straightforward to show that

[TABLE]

where $\mathbf{1}=(1,\ldots,1)^{\prime}\in\mathbb{R}^{d}$ and $\mathbf{W}$ is the set of all vectors in $\mathbb{R}^{d}$ with components summing up to 1. Thanks to Proposition 5, $\bm{w}_{\text{opt}}$ is also the minimizer of the limiting variance of the high quantile estimator from (19).

Let us now turn to the selection of $\mathbf{k}$. The numbers $k_{j}$, $j=1,\ldots,d$, represent the local tail sample lengths that are effectively used in the estimation of the common EVI $\gamma$. These numbers typically are small relative to the full sample lengths so that the Hill estimator and related procedures should be applied only in data-rich situations. Since our local sample lengths $n_{j}$ are rather short, we counteract this local drawback by combining many different stations from a region.

The fact that we deal with block maxima from $d$ different stations motivates us to set $k_{j}=\lfloor 2n_{j}^{2/3}/d^{1/3}\rfloor$ (Kinsvater et al., 2016, Sec. 3.2). Note that the ratio $k_{j}/n_{j}$ is getting smaller with increasing local sample size $n_{j}$. This is needed for the consistency of Hill’s estimator. Accounting for the dimension $d$ in $k_{j}$ turns out to be important in our data applications. We are able to reduce a typically dominant bias by taking advantage of a large dimension $d$.

5 Simulation study

We choose two numbers of sites $d\in\{5,10\}$ in order to consider common regional settings. For simplicity, we assume that all $d$ sites have the same observation period and hence the same sample lengths $n\in\{50,100\}$, which means that $a_{j}=0$ in (10) for all $j=1,\ldots,d$, even though our theory is able to deal with the general case of different sample lengths. The focus lies on the (regional) estimation of a high quantile $q_{p}=F_{j}^{-1}(p)$ for probability $p\in\{0.99,0.999\}$ at some site, say, $j=1$.

We apply Khoudraji’s device (Khoudraji, 1995) in order to generate an inter-site dependence model. We set

[TABLE]

where $C_{\theta}$ denotes the $d$-dimensional Gumbel-Hougaard copula, $\bm{\theta}=(1.5,2.5)^{\prime}$ accounts for strength of dependence and $\mathbf{c}=1/d\cdot(0,1,2,\ldots,d-1)^{\prime}$ for the amount of asymmetry. This model is simple to handle but still allows for some sort of dependence asymmetry. The latter being a typical feature of discharge data.

We compare the following high quantile estimators computed from annual maxima $X$:

Weissman’s estimator (18) with $\mathbf{k}$ and $\bm{w}$ according to Section 4.2, (W) 2. -

$L$-moment estimator from Appendix A and (L) 3. -

Trimmed $L$-moment estimator from Appendix A. (TL)

(L) and (TL) are improved regional estimators based on the method of probability weighted moments proposed by Lilienthal et al. (2016). After combining the estimation methods, the sample lengths and the number of sites, we get different scenarios to look at. Each scenario is replicated 5000 times.

Recall that in the introduction we identified two potential sources of deviations of $F$ from a GEV distribution: Finite block length $b$ and seasonal variability in general. The next two subsections are devoted to each of these sources.

5.1 Block maxima observations

Annual maximal flows are computed by maximization over finite blocks of, say, monthly observations and thus, it is common to assume that annual maxima are exactly GEV distributed. This section analyzes the behavior of the estimators when observations are in fact generated by block maxima distributions. More precisely, let $T_{\nu}$ denote the $t$-distribution function with $\nu$ degrees of freedom. We choose

[TABLE]

where in general $b\geq 2$. It can be shown that $F_{j}^{(b)}$ is the distribution function of a block maximum

[TABLE]

over independent and identically $t_{1/\xi}$-distributed random variables $Z_{1},\ldots,Z_{b}$. The standardization is used in order to obtain $\lim_{b\rightarrow\infty}F_{j}^{(b)}=G_{\mu_{j},\sigma_{j},\xi}$, which is the GEV distribution function with location $\mu_{j}$, scale $\sigma_{j}$ and shape $\xi$.

Here, we have a look at annual maximal flows, so we set $b=12$. For the parameters we choose $(\mu_{j},\sigma_{j},\xi)^{\prime}=(1.75,1,0.3)^{\prime}$ for all $j=1,\ldots,d.$

For the different scenarios, i.e., dimensions $d$, sample lengths $n$ and probabilities $p$, we came to the same conclusion, so just one particular scenario with $d=10$, $n=50$ and $p=0.99$ is depicted in Figure 4. The true 99%-quantile of this scenario is $q_{0.99}=14.151\text{ m}^{3}/\text{s}$, which corresponds to the dashed line in the figure. All competitors show a similar bias, caused by finite $b$ for the parametric estimators (model misspecification) and by finite $n$ for the Weissman approach. We complete the analysis of finite blocks $b$ and continue with seasonal observations under the GEV assumption.

5.2 Seasonal observations

Now, according to (3) and the related discussion, we set $F_{j}=F_{w,j}\cdot F_{s,j}$, where each of the seasonal distributions is a GEV with parameters $\mu_{w,j},\sigma_{w,j},\xi_{w,j}$ and $\mu_{s,j},\sigma_{s,j},\xi_{s,j}$. Here, besides the methods (W), (L) and (TL) used before, we also apply seasonal estimators from Section B.2, namely

seasonal Weissman’s estimator, (sW) 2. -

seasonal $L$-moments and (sL) 3. -

seasonal Trimmed $L$-moments. (sTL)

Based on experience from real datasets we choose the parameters $(\mu_{w,j},\sigma_{w,j},\xi_{w,j})^{\prime}=(2,1,0.2)^{\prime}$ for the winter and $(\mu_{s,j},\sigma_{s,j},\xi_{s,j})^{\prime}=(1.5,1,0.4)^{\prime}$ for the summer season. The true quantile in this scenario is $q_{0.99}=15.692\text{ m}^{3}/\text{s}$.

Figure 5 shows the boxplots for the regional scenario with $d=10$, $n=50$ and $p=0.99$. We omit the results for other combinations of $(p,n,d)$ because the conclusions drawn are similar.

The annual L- and TL-estimators are based on misspecified assumptions, since here we consider observations from two-component GEV’s. In our simulation experiments this is reflected by a non-negligible bias. Annual Weissman (W) and seasonal TL-moments (sTL) outperform their competitors, probably because they do not suffer from the previous misspecification issue. However, this is also true for seasonal Weissman (sW), and still its performance suffers from a devastating bias. We are not able to explain this poor behavior and can neither prove the presumed asymptotic normality of estimator sW.

Summing up, we do not recommend sW. We neither recommend L and TL if strong seasonal variability is present. Instead, we suggest sTL if one assumes that variability is not an issue within each season. Otherwise, we advise to apply W, provided that the dimension $d$ is not very small. Note that local estimation (not reported here) with the Weissman approach can be not recommended due to the typically dominant bias.

A comment on alternative scenarios

We also studied several other scenarios that are not reported here because of qualitatively very similar conclusions. For instance, we studied the performance under a variety of plausible parameter values for $(\mu_{w,j},\sigma_{w,j},\xi_{w,j})^{\prime}$, $(\mu_{s,j},\sigma_{s,j},\xi_{s,j})^{\prime}$, $(\mu_{j},\sigma_{j},\xi)^{\prime}$ and local sample lengths $n_{j}$, where plausibility is based on our personal experience from real data applications. We also looked at scenarios where both sources of deviations of $F$ from a GEV model, seasonal variability and block maxima observations for each season, have been used. Summing up, we always drew conclusions very similar to that of the experiments reported here.

6 Application: The Mulde river basin

The Mulde river basin with a catchment area of about 7400 km2 is located in Saxony, a federal state in the East of Germany. The basin is named after a main tributary of the Elbe river, the “Vereinigte Mulde”, with an average flow at the mouth of about 73 m3/s. From a flood frequency analysts point of view, this region is particularly interesting because there were two flood disasters in the recent past. The first one was in 2002 and the second one in 2013. Both were summer events caused by heavy rainfall during some days. The damage caused amounts to several billion euros. It is thus of particular relevance to account for risk adequately.

We study a data set of monthly maximal flows from $d=8$ sites: Lichtenwalde, Golzern, Zickau-Pölbitz, Wechselburg, Nossen, Hopfgarten, Pockau and Borstendorf. Their local sample lengths $n_{j},j=1,\ldots,d,$ are varying between 83 and 102 years. A hydrological year runs from November to October and typical seasons in this region are winter from November to April and summer from May to October. We calculate winter, summer and annual maximal flows $W_{i,j}$, $S_{i,j}$ and $X_{i,j}=\max\{W_{i,j},S_{i,j}\}$ for each available combination $(i,j)$ of year and station resulting in schemes as in (10). In order to apply our regional estimators to the annual (resp. seasonal) observations, it is recommended to check first assumption $\mathcal{H}_{0,\text{tail}}$ from (17) for annual (resp. seasonal) data. The application of homogeneity test based on the main result of Lilienthal et al. (2016) considering the hypothesis of identical shape parameters resulted in $p$-values of 0.41, 0.63 and 0.54 for annual, winter and summer schemes, respectively. This indicates that there are no serious sources of tail heterogeneity. In what follows, since regional estimation works comparably well also under moderate heterogeneity (Lettenmaier et al., 1987), we consider the sites as tail homogeneous.

Figure 6 depicts a return level plot for station Lichtenwalde, which is one of the eight sites of the group. In general, a point in the plot has coordinates \big{(}T,F^{-1}(1-1/T)\big{)}, with $T>1$ called return period and $F^{-1}(1-1/T)$ called $T$-year return level. The three lines in Figure 6 correspond to three regional estimators of $F$: The Weissman estimator, the TL-moment estimator under the usual GEV assumption and the sTL-moment estimator under the two-component GEV assumption. The first two take only annual flows $X_{i,j}$ into account, while the third one is based on seasonal observations $W_{i,j}$ and $S_{i,j}$. These three procedures turn out to be superior to other considered competitors (see Section 5). The dots \big{(}T_{i:n},X_{i:n}\big{)}, $i=1,\ldots,n$, in the figure represent the empirical counterpart: The ordered observations $X_{1:n}\leq\ldots\leq X_{n:n}$ are called empirical return levels with corresponding empirical return periods $T_{i:n}=1/(1-p_{i:n})$ and plotting positions $p_{i:n}=i/(n+1)$.

In Figure 6 we observe that estimated $T$-year return levels of the three competing estimators are pretty much the same up to $T\approx 40$ years. This corresponds to almost the same quantile curves $\hat{F}^{-1}(p)$ for $p\leq 1-1/40=0.975$. Note that the largest observation stands out from the bulk of the data. The corresponding empirical return level of around 100 years is half of the sTL and quarter of both annual estimated return levels. The fact that the seasonal approach is closer to the sample is also the case in many other of our data applications (not reported here). This, indeed, should not be always interpreted as an advantage, because of possible data-overfitting. For instance, the empirical version $F_{n}$ perfectly fits the data even though this estimator is not recommended due to its horrible variability in the tails. However, note also that seasonal estimation is based on extended information: We do not simply take annual $X_{i,j}$ but use variables $W_{i,j}$, $S_{i,j}$ instead and account for possible seasonal variability. In this particular application the variability in the tails is reflected by very different estimates of the shapes, $\hat{\xi}_{w}=0.29$ for winter and $\hat{\xi}_{s}=0.41$ for the summer season. Without accounting for the variability we would estimate the shape as $\hat{\xi}_{\text{annual}}=0.25$, which surprisingly is not in-between the first two values. This counter-intuitive result is explained by the fact that mixing two very different distributions can result in misleading overlapping.

In the end, practitioners are interested in the estimation of, say, the $p=99\%$-quantile or, equivalently, the $T=100$ year return level of the annual maximal flow distribution $F$ of Lichtenwalde. Here we focus on the regional two-component GEV approach from Appendix B.1 based on TL-moments. Our theory allows us to calculate a point estimate $\hat{q}_{p}^{(\text{sTL})}=962.5$ and its corresponding estimated 95%-confidence interval of $\hat{I}_{0.95}^{(\text{sTL})}=[680.9,~{}1244.1]$ by plugging in estimates of all unknown quantities in formula (13). For comparison, the point estimate and 95%-confidence interval for the local sTL (resp. regional W) procedure is given by $\hat{q}_{p}^{(\text{loc sTL})}=1089.6$ and $\hat{I}_{0.95}^{(\text{loc sTL})}=[568.0,~{}1611.2]$ (resp. $\hat{q}_{p}^{(\text{W})}=823.2$ and $\hat{I}_{0.95}^{(\text{W})}=[521.0,~{}1125.4]$), where local means that only observations from Lichtenwalde are used for estimation of the two-component GEV model.

7 Conclusion

We consider asymptotic methods for regional frequency analysis that are able to deal with seasonal variability and inter-site dependence. Our framework is flexible in the sense that we do allow for very different local sample lengths and we do not restrict ourselves to parametric dependence models.

The methods are investigated via simulations and illustrated by a real data application. Taking into account possible seasonal variability can have an enormous impact on the results. For instance, the return period of the severe flood in 2002 estimated by the two-component approach was only half of that by the classical one-component GEV approach. We recommend seasonal models whenever the corresponding observations are available in the given data set.

It even might be useful to consider more than $K=2$ seasons. A corresponding theoretical result, the asymptotic normality for $K$-component GEV quantiles, can be easily deduced along the same lines as the proof of case $K=2$.

Acknowledgements

We are grateful to Professor Andreas Schumann from the Department of Civil Engineering, Ruhr-University Bochum, Germany, for providing us hydrological data and for helpful discussions. The financial support of the Deutsche Forschungsgemeinschaft (SFB 823, “Statistical modelling of nonlinear dynamic processes”) is gratefully acknowledged.

Appendix A GEV parameterizations by probability weighted moments

Let $G_{\theta}$ with parameters $\theta=(\mu,\sigma,\xi)^{\prime}$ denote a GEV distribution function defined in (2). Throughout this section we will assume that $\xi<1$. The probability weighted moment (PWM) of $G_{\theta}$ of order $k\in\mathbb{N}_{0}$ is defined by

[TABLE]

It is a well-known fact that every distribution with finite expectation is uniquely determined by its sequence of PWM’s (Hosking, 2007). It is thus little surprising that the parameter $\theta$ already is determined by a finite number of PWM’s. Two such re-parameterizations of the GEV family are recapped below.

1) Method of L-moments

The most popular re-parameterization of the GEV based on the first three PWM’s $\bm{\beta}=(\beta_{0},\beta_{1},\beta_{2})^{\prime}$ is implicitly given by the equation system

[TABLE]

where $\Gamma$ denotes the gamma function and with $\lambda_{1}=\beta_{0}$, $\lambda_{2}=2\beta_{1}-\beta_{0}$, $\lambda_{3}=6\beta_{2}-6\beta_{1}+\beta_{0}$ denoting the first three so-called L-moments of $G_{\theta}$ (Hosking, 1990). For some reasons to be found in the previous reference, practitioners usually prefer to work with L-moments instead of PWM’s.

However, an explicit solution $\phi$ with $\theta=\phi(\bm{\beta})$ of the previous equation system does not exist. Instead of relying on a numerical solver, it is convenient to replace the first equation in (27) by the approximation

[TABLE]

in order to obtain an explicit solution. The error in $\xi$ caused by the approximation is known to be smaller than 0.0009 for $-1/2<\xi</1/2$ (Hosking et al., 1985) and is therefore negligible in practice.

2) Method of TL-moments

Elamir and Seheult (2003) proposed a generalization of L-moments denoted by $\lambda_{k}^{(r,s)}$, where $(r,s)\in\mathbb{N}_{0}^{2}$ is a trimming parameter: Larger values for $r$ (resp. $s$) decrease the influence of the smallest (resp. largest) observations, with $\lambda_{k}^{(0,0)}=\lambda_{k}$ being the usual (untrimmed) L-moments. For further insight we refer to Hosking (2007).

For the reminder we focus on the trimming $(r,s)=(0,1)$, which means that only the influence of the largest observations is slightly reduced and which has proven to be a good choice in our simulation experiments. Again, the procedure for the derivation of GEV parameters will be based only on the first three trimmed L-moments and, thanks to Hosking (2007, Sec. 2.4), we have that $\lambda_{1}^{(0,1)}=2\beta_{0}-2\beta_{1}$, $\lambda_{2}^{(0,1)}=3/2(4\beta_{1}-\beta_{0}-3\beta_{2})$ and $\lambda_{3}^{(0,1)}=2/3(36\beta_{2}-18\beta_{1}+2\beta_{0}-20\beta_{3})$. It can be verified that the equation system

[TABLE]

determines the GEV parameters. Similar to the first procedure, we replace the first equation in (32) by the approximation

[TABLE]

in order to obtain an explicit solution, where now $\bm{\beta}=(\beta_{0},\beta_{1},\beta_{2},\beta_{3})^{\prime}$ denotes the vector of the first four PWM’s of $G_{\theta}$.

Appendix B Asymptotic statistics

B.1 A CLT for PWMs with applications to GEV parameter estimation

We recap the main result from Lilienthal et al. (2016). These authors verify the joint asymptotic normality of sample PWMs

[TABLE]

computed from scheme (10).

Theorem 1**.**

Suppose that $\mathbf{X}_{i}$, $i\geq 1$, is a sequence of independent copies of $\mathbf{X}=(X_{1},\ldots,X_{d})^{\prime}$ with marginal distribution functions $F_{j}(x)=\mathbb{P}(X_{j}\leq x)$ and with marginal PWMs

[TABLE]

We further assume that

[TABLE]

and that $\sup_{x\in\mathbb{R}}\left|x\{F_{j}(x)(1-F_{j}(x))\}^{w}\right|<\infty$ for all $j=1,\ldots,d$ and some $w\in[0,1/2)$. Let $\hat{\bm{\beta}}_{\bm{r},n}$ denote the sample counterpart of $\bm{\beta}$ computed from scheme (10) with $n_{j}/n\rightarrow r_{j}\in\mathbb{R}_{+}$ for $n\rightarrow\infty$. Then, for $n\rightarrow\infty$, we have that

[TABLE]

where the limiting variance matrix $\Sigma_{\bm{r}}\in\mathbb{R}^{dK\times dK}$ is defined below.

Lilienthal et al. (2016) proved that the matrix $\Sigma_{\bm{r}}=\lim_{n\rightarrow\infty}{\rm Var}\left(\sqrt{n}\left(\hat{\bm{\beta}}_{\bm{r},n}-\bm{\beta}\right)\right)$ is defined block-wise by

[TABLE]

where $\bm{\beta}_{j}=\left(\beta_{0}(F_{j}),\ldots,\beta_{K}(F_{j})\right)^{\prime}$ with sample counterpart $\hat{\bm{\beta}}_{j,r_{j},n}$ computed from the $j$-th row of scheme 10. The random vectors $\mathbf{Z}_{j}=(Z_{0,j},Z_{1,j},\ldots,Z_{K-1},j)^{\prime}$, $j=1,\ldots,d$, are defined componentwise by

[TABLE]

The matrices ${\rm Cov}(\mathbf{Z}_{j},\mathbf{Z}_{\ell})$ are consistently estimated by their sample analogues: Let

[TABLE]

and $\hat{\mathbf{Z}}_{i,j}=(Z_{i,0,j},Z_{i,1,j},\ldots,Z_{i,K-1,j})^{\prime}$, $i=a_{j}+1,\ldots,n$. For $1\leq j,\ell,\leq d$, the covariance matrix ${\rm Cov}(\mathbf{Z}_{j},\mathbf{Z}_{\ell})$ is estimated by the empirical covariance matrix of the sample

[TABLE]

This empirical estimator of $\Sigma_{\bm{r}}$ is denoted by $\hat{\Sigma}_{\bm{r},n}$.

Note that the preceding theorem is purely non-parametric, i.e., except for continuity and the existence of some moments no further assumptions on the distribution are made. For the reminder of this section we consider a parametric problem, the estimation of GEV parameters under a regional homogeneity assumption: Suppose now that $F_{j}=G_{\mu_{j},\sigma_{j},\xi_{j}}$ is a GEV distribution function for all $j=1,\ldots,d$ and that the tail homogeneity assumption

[TABLE]

holds. Let $\hat{\theta}_{j,r_{j},n}=(\hat{\mu}_{j},\hat{\sigma}_{j},\hat{\xi}_{j})^{\prime}$ denote either the L- or TL-moment estimator of $\theta_{j}=(\mu_{j},\sigma_{j},\xi_{j})^{\prime}$ computed from the $j$-th row of scheme (10), with computation based on sample PWMs and the corresponding representation from Appendix A. Let $\bm{w}\in\mathbb{R}^{d}$ denote any vector of weights summing up to 1. Thanks to assumption (36), we can estimate $\theta_{j}$ consistently by the regional version \big{(}\hat{\mu}_{j},\hat{\sigma}_{j},\bm{w}^{\prime}\hat{\bm{\xi}}\big{)}. The asymptotic normality and the corresponding limiting variance can be easily derived from the previous theorem and by the delta method.

Just like for the regional Weissman estimator it is meaningful to set

[TABLE]

where $\Sigma_{\text{tail}}=\lim_{n\rightarrow\infty}{\rm Var}\left[\sqrt{n}(\hat{\bm{\xi}}-\bm{\xi})\right]$ with $\hat{\bm{\xi}}=(\hat{\xi}_{1},\ldots,\hat{\xi}_{d})^{\prime}$. Again, $\Sigma_{\text{tail}}$ is obtained by Theorem 1 and the delta method. Finally, we set $\hat{\xi}=\hat{\bm{w}}_{\text{tail}}^{\prime}\hat{\bm{\xi}}$ as our regional (T)L-moment estimator of the regional shape $\xi$, where $\Sigma_{\text{tail}}$ involved in $\bm{w}_{\text{tail}}$ is replaced by a consistent estimator $\hat{\Sigma}_{\text{tail}}$.

Lilienthal et al. (2016) introduced also a check version $\check{\Sigma}_{\bm{r},n}$ for the estimation of $\Sigma_{\bm{r}}$ working under the GEV assumption. $\check{\Sigma}_{\bm{r},n}$ is a hybrid version, consisting of parametric $3\times 3$-blocks on the diagonal and non-parametric ones outside the diagonal. Throughout the simulation section we applied the check-version because of its superior efficiency. However, in finite samples this approach sometimes leads to invalid covariance matrices reflected by negative eigenvalues. In such cases we simply replace $\bm{w}_{\text{tail}}$ by weights $\bm{w}_{\text{ind}}=(n_{1},n_{2},\ldots,n_{d})^{\prime}/\sum_{j=1}^{d}n_{j}$ proportional to local sample lengths $n_{j}$, which is the optimal choice under the assumption of spatial independence.

B.2 Quantiles of two-component GEV distributions

Let $F(x)=G_{\theta_{w}}(x)\cdot G_{\theta_{s}}(x)$ be a two-component GEV distribution, where $\theta_{w},\theta_{s}$ are two unknown parameter vectors and where $G_{\theta}$ is defined in (2). Let $p\in(0,1)$ be fixed. Suppose that our interest is in the estimation of

[TABLE]

and in assessing the estimation uncertainty of $\hat{q}_{p}$. We will further assume that $\hat{\theta}_{w}=\hat{\theta}_{w,n}$ and $\hat{\theta}_{s}=\hat{\theta}_{s,n}$ are independent satisfying

[TABLE]

for $n\rightarrow\infty$.

Theorem 2**.**

Let $g_{\theta}$ denote the density of $G_{\theta}$ and $J_{\theta}(x)=\frac{\partial}{\partial\theta}G_{\theta}(x)\in\mathbb{R}^{1\times 3}$ the Jacobi matrix of $G_{\theta}(x)$ in $\theta$. Then, for $n\rightarrow\infty$, we have that

[TABLE]

where the limiting variance is given by

[TABLE]

Proof.

Let $p\in(0,1)$ be fixed and let $\varphi:\Theta\times\Theta\subset\mathbb{R}^{6}\rightarrow\mathbb{R}$ denote the map $(\theta_{1},\theta_{2})\mapsto\left(G_{\theta_{1}}\cdot G_{\theta_{2}}\right)^{-1}(p)$. We are going to prove that $\varphi$ is a differentiable map in $(\theta_{w},\theta_{s})$.

Let $J$ denote the support of $F$ and let $I\subset J$ be some compact interval with interior point $q_{p}\in I$. Note that $\varphi$ can be decomposed into $\varphi=\psi\circ\eta$ with

[TABLE]

By Lemma 3 below and a simple multivariate extension of Kosorok (2008, Lem. 12.2), it is straightforward to verify that $\eta$ is Hadamard differentiable in $(\theta_{w},\theta_{s})$ with derivative $D\eta_{(\theta_{w},\theta_{s})}:\mathbb{R}^{6}\rightarrow\ell^{\infty}(I)$,

[TABLE]

where $J(\theta,x)=\frac{\partial}{\partial\theta}G_{\theta}(x)\in\mathbb{R}^{1\times 3}$ is the Jacobi matrix of $G_{\theta}(x)$ in $\theta$ and where $\times$ denotes matrix multiplication. By Van der Vaart (2000, Lem. 3.9.20), $\psi$ is Hadamard differentiable in $F=\eta(\theta_{w},\theta_{s})\in\ell^{\infty}(I)$ with derivative $D\psi_{F}:\ell^{\infty}(I)\rightarrow\mathbb{R}$,

[TABLE]

where $f\in\ell^{\infty}(I)$ is the derivative of $F$.

Since both, $\eta$ and $\psi$, are Hadamard differentiable, so is the composition $\varphi=\psi\circ\eta$. But this is equivalent to $\varphi$ being differentiable in the usual sense (Van der Vaart, 2000). The chain rule is also valid for Hadamard differentiability, which implies that the Jacobi matrix $D\varphi_{(\theta_{w},\theta_{s})}$ is given by

[TABLE]

Recall that $F=\eta(\theta_{w},\theta_{s})$, $q_{p}=F^{-1}(p)$ and $f(x)=g_{\theta_{w}}(x)\cdot G_{\theta_{s}}(x)+G_{\theta_{w}}(x)\cdot g_{\theta_{s}}(x)$. By the delta method and for $n\rightarrow\infty$, we finally conclude that

[TABLE]

where

[TABLE]

∎

Lemma 3**.**

Let $I_{0}$ denote the support of $G_{\theta_{0}}$ and let $I\subset I_{0}$ be some compact interval. The map $\kappa:\Theta\rightarrow\ell^{\infty}(I)$, $\theta\mapsto G_{\theta}$, is Hadamard-differentiable in $\theta_{0}$ with derivative $D\kappa_{\theta_{0}}:\Theta\rightarrow\ell^{\infty}(I)$, $\mathbf{a}\mapsto J_{\theta_{0}}\mathbf{a}$, where $J_{\theta}(x)=\frac{\partial}{\partial\theta}G_{\theta}(x)$, $x\in I$.

Proof.

Hadamard differentiability is defined as a particular kind of uniform convergence of the difference quotient of $\kappa$. We have to prove that

[TABLE]

and any $\theta_{0}+t\mathbf{a}_{t}\in\Theta$ with $\mathbf{a}_{t}\rightarrow\mathbf{a}$ for $t\rightarrow 0$. We consider the norm $\|f\|=\sup_{x\in I}|f(x)|$ on $\ell^{\infty}(I)$. It is a known fact that uniform convergence on compact sets follows by point-wise convergence and by equicontinuity.

Note that point-wise convergence, i.e., convergence in (40) for fixed $x\in I$ follows by the differentiability of the map $\theta\mapsto G_{\theta}(x)$ in $\theta_{0}$. Since $I$ is a compact set, it suffices to prove equicontinuity of the family $\{f_{t}\in\ell^{\infty}(I):~{}t\in[-1,0)\cup(0,1]\}$ with

[TABLE]

Proof of equicontinuity: Let $\varepsilon>0$ and $x_{1},x_{2}\in I$ with $|x_{1}-x_{2}|<\delta$ ($\delta$ to be defined later on). By the mean value theorem there exists a real number $x_{0}$ between $x_{1}$ and $x_{2}$, such that

[TABLE]

where $f^{\prime}_{t}(x)=\frac{\partial}{\partial x}f_{t}(x)=t^{-1}\left[g_{\theta_{0}+t\mathbf{a}_{t}}(x)-g_{\theta_{0}}(x)\right]$ and $g_{\theta}$ is the density of $G_{\theta}$. If we show that $\left|f^{\prime}_{t}(x_{0})\right|\leq c$ is bounded uniformly for all $x_{0}\in I$ and $t\in[-1,1]\backslash\{0\}$, we can set $\delta=\varepsilon/c$ and the proof of the equicontinuity is complete.

Proof of boundedness: Note first that the density $g_{\theta}$ with $\theta=(\mu,\sigma,\xi)\in\Theta$ is bounded by $|g_{\theta}(x)|\leq\sigma^{-1}|1+\xi|$. Without loss of generality we may assume that $\theta_{0}+t\mathbf{a}_{t}$ is bounded uniformly in $t\in[-1,1]\backslash\{0\}$, since $\mathbf{a}_{t}\rightarrow\mathbf{a}$ for $t\rightarrow 0$. This allows us to find a bound for $|f^{\prime}_{t}(x)|$ uniformly in $t\in[-1,-a]\cup[a,1]$ and $x\in I$, where $0<a<1$ is some constant. On the other hand (for $t$ close to zero), note that $\lim_{t\rightarrow 0}f^{\prime}_{t}(x)=\frac{\partial}{\partial\theta}g_{\theta}(x)|_{\theta=\theta_{0}}\mathbf{a}$ and that the limit is uniformly bounded for $x\in I$. Summing up, we conclude that $|f^{\prime}_{t}(x)|$ is bounded uniformly for $t\in(-a,a)\backslash\{0\}$. This completes the proof of the lemma. ∎

B.3 Semi-parametric inference

B.3.1 Introduction: Hill’s estimator and Weissman’s extrapolation formula

Let $X$ be a random variable with distribution function $F(x)=\mathbb{P}(X\leq x)$ and let $u>0$ be a real positive number. The random variable $X/u$ satisfying $X>u$ is called relative excess over the threshold $u$. Then, from (14), we obtain

[TABLE]

which means that relative excesses over $u$ approximately follow a parametric distribution for large $u$. The limit $P_{\gamma}(x)=1-x^{-1/\gamma}$, $x\geq 1$, is called Pareto distribution function with $\gamma>0$ being a shape parameter. It is straightforward to verify that a random variable $Z$ with distribution function $P_{\gamma}$ fulfills $\gamma=\mathbb{E}[\log Z]$, which suggests to estimate $\gamma$ from arithmetic means of $\log$-transformed excesses:

Let $X_{1},\ldots,X_{n}$ be independent copies of $X$ with ordered values denoted by $X_{(1)}\leq\ldots\leq X_{(n)}$. Hill’s estimator of $\gamma$ proposed in Hill (1975) is defined by

[TABLE]

for integers $1\leq k<n$ representing the number of excesses over the random threshold $u_{k,n}$. The popularity of Hill’s estimator can be explained by its computational simplicity and by the fact that it can be considered as a maximum pseudo-likelihood estimator.

Now, we choose $\hat{\gamma}=H_{k,n}$ as Hill’s estimator. In the next step, we estimate quantiles $F^{-1}(p)$ with an extrapolation formula of Weissman (1978), which is motivated by relation (15) and given by

[TABLE]

The first term $u_{k,n}=X_{(n-k)}$ is interpreted as a non-parametric estimator of a moderate quantile $F^{-1}(1-k/n)=U(n/k)$, while the remainder is used for the extrapolation into the tail region $p>1-k/n$, with $\hat{\gamma}$ controlling the extrapolation width.

B.3.2 Limit theorems for regional semi-parametric estimation

The following results are built upon technical assumptions, that, however, can be checked only if detailed information on the tail region of the joint distribution is available. Since our observation lengths are very limited, we omit these details and refer to the assumptions stated in Kinsvater et al. (2016, Section 2) and de Haan and Ferreira (2006, Theorem 4.3.8).

The following statistics are computed from observations of scheme (10), where $a_{j}=\lfloor n(1-r_{j})\rfloor$ such that $n_{j}/n\rightarrow r_{j}>0$ for $n\rightarrow\infty$ and with random vectors $\mathbf{X}_{i}=(X_{i,1},\ldots,X_{i,d})^{\prime}$ having marginal distribution functions $F_{j}=\mathbb{P}(X_{i,j}\leq x)$ of Pareto-type. The random vector $\mathbf{H}_{\mathbf{k},\bm{r},n}$ is a collection of local Hill-estimators, with $j$-th component denoted by $H_{k_{j},r_{j},n}$ computed from the $k_{j}$ largest order statistics of the sample $X_{a_{j}+1,j},X_{a_{j}+2,j},\ldots,X_{n,j}$.

Proposition 4**.**

(Proposition 1 in Kinsvater et al. (2016))*

Suppose that $\gamma_{1}=\ldots=\gamma_{d}=\gamma$, $k_{1}\rightarrow\infty$, $k_{1}/n\rightarrow 0$ and $k_{1}/k_{\ell}\rightarrow c_{\ell}$ hold for $n\rightarrow\infty$ and finite values $c_{\ell}>0$, $\ell=2,\ldots,d$. Let $\mathbf{1}=(1,\ldots,1)^{\prime}\in\mathbb{R}^{d}$. Then, under additional technical assumptions (Kinsvater et al., 2016, Sec. 2) and for $n\rightarrow\infty$, we have that*

[TABLE]

holds, where $\Sigma\in\mathbb{R}^{d\times d}$ is defined componentwise by

[TABLE]

for $1\leq l,m\leq d$ and where $x\wedge y=\min(x,y)$. For $l=m$ this reduces to $\Sigma_{l,l}=c_{l}$.

For the estimation of $\Sigma$ we plug in $\hat{c}_{\ell}=\frac{k_{1}}{k_{\ell}}$ and a consistent estimator $\hat{\Lambda}_{\ell,m}$ of the upper tail dependence copula $\Lambda_{\ell,m}$ between components $\ell$ and $m$. For the latter, we apply two different solutions: The empirical estimator studied in Schmidt and Stadtmüller (2006) or, under the extreme value dependence assumption, we use representation

[TABLE]

with corresponding Pickands function $A_{\ell,m}$ and plug in the corrected CFG-estimator $\hat{A}_{\ell,m}$ studied in Genest and Segers (2009). As a direct consequence of Proposition 4, we have

[TABLE]

where $\hat{\gamma}_{\mathbf{k},\bm{r},n}(\mathbf{w})=\mathbf{w}^{\prime}\mathbf{H}_{\mathbf{k},\bm{r},n}$ is called regional estimator of $\gamma$.

Let $\hat{F}_{j}^{-1}(p)$ denote the regional estimator of $F_{j}^{-1}(p)$ defined in (18).

Proposition 5**.**

(we refer to Kinsvater (2016, Sec. 2.7))*

Let $1\leq j\leq d$ be fixed. Let $p=p_{n}\in(0,1)$ be a sequence of numbers such that $p\rightarrow 1$, $n(1-p)/k_{1}\rightarrow 0$ and $\log(np)/\sqrt{k_{1}}\rightarrow 0$ hold for $n\rightarrow\infty$. Then, under the assumptions from Proposition 4, the same technical assumptions as in Proposition 4 and for $n\rightarrow\infty$, we have that*

[TABLE]

where $\hat{F}_{j}^{-1}(p)$ is defined in (42).

Proposition 5 allows us to derive an asymptotic $(1-\alpha)$-confidence interval

[TABLE]

where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution and with $a\cdot(1\pm b)$ denoting the interval $[a\cdot(1-b),\ a\cdot(1+b)]$ for $a\in\mathbb{R}$, $b>0$.