The transformed-stationary approach: a generic and simplified methodology for non-stationary extreme value analysis

Abstract

Statistical approaches to study extreme events require, by definition, long time series of data. In many scientific disciplines, these series are often subject to variations at different temporal scales that affect the frequency and intensity of their extremes. Therefore, the assumption of stationarity is violated and alternative methods to conventional stationary extreme value analysis (EVA) must be adopted. Using the example of environmental variables subject to climate change, in this study we introduce the transformed-stationary (TS) methodology for non-stationary EVA. This approach consists of (i) transforming a non-stationary time series into a stationary one, to which the stationary EVA theory can be applied, and (ii) reverse transforming the result into a non-stationary extreme value distribution. As a transformation, we propose and discuss a simple time-varying normalization of the signal and show that it enables a comprehensive formulation of non-stationary generalized extreme value (GEV) and generalized Pareto distribution (GPD) models with a constant shape parameter. A validation of the methodology is carried out on time series of significant wave height, residual water level, and river discharge, which show varying degrees of long-term and seasonal variability. The results from the proposed approach are comparable with the results from (a) a stationary EVA on quasi-stationary slices of non-stationary series and (b) the established method for non-stationary EVA. However, the proposed technique comes with advantages in both cases. For example, in contrast to (a), the proposed technique uses the whole time horizon of the series for the estimation of the extremes, allowing for a more accurate estimation of large return levels. Furthermore, with respect to (b), it decouples the detection of non-stationary patterns from the fitting of the extreme value distribution. As a result, the steps of the analysis are simplified and intermediate diagnostics are possible. In particular, the transformation can be carried out by means of simple statistical techniques such as low-pass filters based on the running mean and the standard deviation, and the fitting procedure is a stationary one with a few degrees of freedom and is easy to implement and control. An open-source MATLAB toolbox has been developed to cover this methodology, which is available at https://github.com/menta78/tsEva/ (Mentaschi et al., 2016).

Introduction

Extreme value analysis (EVA) attains a great importance in several applied sciences, particularly in earth science, because it is a fundamental tool to study the magnitude and frequency of extreme events and their changes (e.g., Alfieri et al., 2015; Forzieri et al., 2014; Jongman et al., 2014; Resio and Irish, 2015; Vousdoukas et al., 2016a). Climatic extreme events are usually associated with disasters and damages with significant social and economic costs. A correct statistical evaluation of the strength of extreme events related to their average return period is crucial for impact assessment, for the evaluation of the risks affecting human lives and activities, and for planning actions regarding risk management and prevention (e.g., Hirsch and Archfield, 2015; Jongman et al., 2014).

Often it is necessary to apply EVA to non-stationary time series, i.e., series with statistical properties that vary in time due to changes in the dynamic system. In particular, climate change can induce variations in the statistical properties of time series of climatic variables. For example, an intensification of the meridional thermal gradient at middle latitudes on a global scale would lead to an increase of the climatic variability (e.g., Brierley and Fedorov, 2010), resulting in a reduction of the average return period of storms with a given strength. Consequently, in the study of climate change, an accurate statistical estimation of middle to long-term extremes is inherently connected to the application of non-stationary methodologies.

While a general theory about non-stationary EVA has not yet been formulated (Coles, 2001), there are several studies describing methodologies for the estimation of time-varying extreme value distributions on non-stationary time series, which rely on the pragmatic approach of using the standard extreme value theory as a basic model that can be further enhanced with statistical techniques (e.g., Coles, 2001; Davison and Smith, 1990; Hüsler, 1984; Leadbetter, 1983; Méndez et al., 2006).

An established technique consists in expressing the parameters of an extreme value distribution as time-varying parametric functions ($M$) of time for some custom parameters (${\mathit{\alpha }}_i$, ${\mathit{\beta }}_i$, ${\mathit{\gamma }}_i$ …). By means of a fitting process such as the maximum likelihood estimator (MLE), it is then possible to fit the values of (${\mathit{\alpha }}_i$, ${\mathit{\beta }}_i$, ${\mathit{\gamma }}_i$ …) to model the extremes of the non-stationary series. Appropriate implementations of such a methodology, hereinafter referred to as the established method (EM), produce meaningful results, as proved by a number of contributions (e.g., Cheng et al., 2014; Gilleland and Katz, 2016; Izaguirre et al., 2011; Méndez et al., 2006; Menéndez et al., 2009; Mudersbach and Jensen, 2010; Russo et al., 2014; Sartini et al., 2015; Serafin and Ruggiero, 2014).

A drawback of this approach is that there is no general indication on how to formulate the function $M$. As a rule the model should be as simple as possible. For this reason, typically several formulations of $M$ are tested, and then the best model is chosen through a balance between high likelihood and low degrees of freedom, for example by means of the Akaike criterion (Akaike, 1973). Furthermore, the choice of $M$ depends on the statistical model chosen for the extreme value analysis: for example, for the same series the $M$ used for the generalized extreme value (GEV) model is different from the $M$ used for the generalized Pareto distribution (GPD) model. Moreover, the EM requires non-stationary statistical fitting techniques that are relatively complex to implement and control, because the detection of the time-varying properties of the series is incorporated into the fitting of the extreme value distribution.

Another commonly used approach for dealing with non-stationary series is to divide them into quasi-stationary slices and apply the stationary theory to each slice (e.g., Vousdoukas et al., 2016a). This technique is referred to in the text as “stationary on slice” (SS). Although this technique enables the detection of meaningful trends for short return periods, it has the drawback of reducing the size of the sample used for the EVA, implying larger uncertainty in the estimation of long return periods.

This study aims to contribute to the field of non-stationary EVA by introducing the transformed-stationary (TS) extreme value methodology, which decouples the analysis of the non-stationary behavior of the series from the fitting of the extreme value distribution. For this purpose, it introduces a standard methodology to model the variations of the statistical properties of the series.

The remainder of the paper is structured as follows. In Sect. 2, the TS methodology is described and discussed in a general and theoretic way and implementation details are outlined. In Sect. 3, the validation of the methodology is presented. Section 4 illustrates a comparison with other common approaches for the EVA of non-stationary series, such as EM and SS for modeling time series characterized by seasonal cycles and time series showing long-term trends. In Sect. 5, the results are discussed and in Sect. 6, the most important conclusions are drawn.

Methods and data

Theoretical background

The TS methodology consists of three steps: transforming a non-stationary time series $y\left(t\right)$ into a stationary series $x\left(t\right)$, performing a stationary EVA, and back-transforming the resulting extreme value distribution into a time-dependent one.

The transformation $y\left(t\right)$ $\to$ $x\left(t\right)$ we propose is $$x\left(t\right)=f(y,t)={\displaystyle \frac{y\left(t\right)-T_y\left(t\right)}{C_y\left(t\right)}},$$ where $T_y\left(t\right)$ is the trend of the series, i.e., a curve representing the long-term, slowly varying tendency of the series, and $C_y\left(t\right)$ is the long-term, slowly varying amplitude of a confidence interval that represents the amplitude of the distribution of $y\left(t\right)$. In particular, if $C_y\left(t\right)$ equals the long-term varying standard deviation $S_y\left(t\right)$ of the series $y\left(t\right)$, Eq. (1) reduces to a simple time-varying renormalization of the signal: $$x\left(t\right)=f(y,t)={\displaystyle \frac{y\left(t\right)-T_y\left(t\right)}{S_y\left(t\right)}}.$$ For simplicity, in the remainder of this paper we will limit our analysis to Eq. (2), knowing that all the considerations can be easily extended to any time-varying confidence interval $C_y\left(t\right)$.

Equation (2) guarantees that the average of $x\left(t\right)$ and its standard deviation are uniform in time, which is a necessary condition for $x\left(t\right)$ to be stationary. In particular, the transformed signal $x\left(t\right)$ has a mean equal to 0 and a variance equal to 1. It is worth noting that the transformed series $x\left(t\right)$ is not necessarily stationary: a series with a constant trend and a uniform standard deviation may still have a time-dependent auto-covariance that would invalidate the hypothesis of stationarity (i.e., the condition of a series with statistical moments constant in time). Before proceeding with the analysis, therefore, a stationarity test should be carried out to ensure that $x\left(t\right)$ is stationary and that its annual maxima can be fitted by a stationary extreme value distribution. For example, a simple test can be performed to ensure that higher order statistics such as skewness and kurtosis are roughly constant along the series.

Once the hypothesis of stationarity of $x\left(t\right)$ is verified, we can estimate the distribution GEV${}_X\left(x\right)$ that best fits its extremes, for example through MLE. GEV${}_X\left(x\right)$ is then given by $${\mathrm{\text{GEV}}}_X\left(x\right)=\text{Pr}(X<x)=\exp \left\{-,{\left[\mathrm{1},+,{\mathit{\epsilon }}_x,\left({\displaystyle \frac{x-{\mathit{\mu }}_x}{{\mathit{\sigma }}_x}}\right)\right]}^{-\mathrm{1}/{\mathit{\epsilon }}_x}\right\},$$ where the shape (${\mathit{\epsilon }}_x$), scale (${\mathit{\sigma }}_x$), and location (${\mathit{\mu }}_x$) parameters do not depend on the time. To find the time-dependent distribution GEV${}_Y(y$, $t)$ that fits the non-stationary time series $y\left(t\right)$ we can write that $$\begin{array}{cc}{\text{GEV}}_Y\left(y\right)& =\text{Pr}\left[Y\right(t)<y]=\text{Pr}\left[f^{-\mathrm{1}},(,X,,,t,),<,y\right]\\ & & =\text{Pr}[X<f(y,t\left)\right]={\text{GEV}}_X\left[f\right(y,t\left)\right],\end{array}$$ where $f(y$, $t)$ is the transformation from $y$ to $x$ given by Eq. (1), and $f^{-\mathrm{1}}(x$, $t)$ is its inverse: $$f^{-\mathrm{1}}(x,t)=y\left(t\right)=S_y\left(t\right)\cdot x+T_y\left(t\right).$$ It is always possible to compute GEV${}_Y(y$, $t)$ from GEV${}_X\left(x\right)$ because $f(y$, $t)$ is a monotonically increasing function of $y$ for every time $t$, because the standard deviation $S_y\left(t\right)$ is always positive.

Using Eqs. (3) and (5) in Eq. (4), we find $$\begin{array}{cc}{\text{GEV}}_Y& (y,t)={\text{GEV}}_X\left[f\right(y,t\left)\right]\\ & =\exp \left\{-,{\left[\mathrm{1},+,{\mathit{\epsilon }}_x,\left({\displaystyle \frac{\frac{y-T_y\left(t\right)}{S_y\left(t\right)}-{\mathit{\mu }}_x}{{\mathit{\sigma }}_x}}\right)\right]}^{-\mathrm{1}/{\mathit{\epsilon }}_x}\right\}\\ & & =\exp \left\{-,\left[\mathrm{1},+,{\mathit{\epsilon }}_x,{\left({\displaystyle \frac{y-T_y\left(t\right)-{\mathit{\mu }}_x\cdot S_y\left(t\right)}{{\mathit{\sigma }}_x\cdot S_y\left(t\right)}}\right)}^{-\mathrm{1}/{\mathit{\epsilon }}_x}\right]\right\}.\end{array}$$ Therefore, if $x\left(t\right)$ is fitted by the stationary distribution GEV${}_X\left(x\right)$ then $y\left(t\right)$ is fitted by the time-dependent distribution GEV${}_Y(y$, $t)$ with shape, scale, and location parameters given by $$\begin{array}{}& & {\mathit{\epsilon }}_y={\mathit{\epsilon }}_x,& & {\mathit{\sigma }}_y\left(t\right)=S_y\left(t\right)\cdot {\mathit{\sigma }}_x,& & {\mathit{\mu }}_y\left(t\right)=S_y\left(t\right)\cdot {\mathit{\mu }}_x+T_y\left(t\right).\end{array}$$ It can be shown that the time-dependent GEV parameters given by Eqs. (7)–(9) are the same as the time-varying parameters ${\mathit{\epsilon }}_{\text{ns}}$, ${\mathit{\sigma }}_{\text{ns}}$, and ${\mathit{\mu }}_{\text{ns}}$ of a non-stationary distribution GEV${}_{\text{ns}}$ that would be obtained from a non-stationary MLE on the series $y\left(t\right)$, and which are given by $$\begin{array}{}& & {\mathit{\epsilon }}_{\text{ns}}=\text{const}.,& & {\mathit{\sigma }}_{\text{ns}}=S_y\left(t\right)\cdot a,& & {\mathit{\mu }}_{\text{ns}}=S_y\left(t\right)\cdot b+T_y\left(t\right),\end{array}$$ for varying parameters $a$ and $b$. In fact, if $p_{GX}\left(x\right)$ is the probability density function (PDF) associated with the distribution GEV${}_X\left(x\right)$, then the MLE for GEV${}_X\left(x\right)$ is estimated so that $$\sum \log \left[p_{GX},\left(,x,\right)\right]=\text{max},$$ which involves vanishing derivatives of Eq. (13) on the GEV parameters ${\mathit{\epsilon }}_x$, ${\mathit{\sigma }}_x$, and ${\mathit{\mu }}_x$. For example, considering the scale parameter ${\mathit{\sigma }}_x$, we see that $$\sum {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}_x}}\log \left[p_{GX},\left(x,,,{\mathit{\sigma }}_x\right)\right]=0.$$ The non-stationary MLE maximizes the log-likelihood of the non-stationary PDF $p_{\text{Gns}}(y$, $t)$ associated with GEV${}_{\text{ns}}$, in function of the parameters $a$ and $b$. For example, considering the parameter $a$, we impose $$\sum {\displaystyle \frac{\partial }{\partial a}}\log \left[p_{\text{Gns}},(,y,,,a,,,t,)\right]=0.$$ Let us assume that $p_{\text{Gns}}(y$, $t)$ coincides with the PDF $p_{GY}(y$, $t)$ associated with the distribution GEV${}_Y(y$, $t)$ given by Eq. (6) and that $a$ $=$ ${\mathit{\sigma }}_x$. Considering that $$p_{GY}(y,t)={\displaystyle \frac{\partial }{\partial y}}{\text{GEV}}_Y(y,t)=p_{GX}\left(x\right){\displaystyle \frac{\partial }{\partial y}}f(y,t)={\displaystyle \frac{p_{GX}\left(x\right)}{S_y\left(t\right)}},$$ we obtain $$\begin{array}{cc}\sum {\displaystyle \frac{\partial }{\partial a}}& \log \left[p_{\text{Gns}},(,y,,,a,,,t,)\right]=\sum {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}_x}}\log \left[p_{GY},\left(y,,,{\mathit{\sigma }}_x,,,t\right)\right]\\ & =\sum {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}_x}}\log \left[{\displaystyle \frac{p_{GX}\left(x,,,{\mathit{\sigma }}_x\right)}{S_y\left(t\right)}}\right]\\ & =\sum {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}_x}}\left\{\log ,\left[p_{GX},\left(x,,,{\mathit{\sigma }}_x\right)\right],-,\log ,\left[S_y,\left(,t,\right)\right]\right\}\\ & & =\sum {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}_x}}\log \left[p_{GX},\left(x,,,{\mathit{\sigma }}_x\right)\right]=\mathrm{0},\end{array}$$ where the last step is possible because $S_y\left(t\right)$ does not depend on ${\mathit{\sigma }}_x$.

The same principle can be applied differentiating $\sum \log \left[p_{GY}\right(x$, ${\mathit{\mu }}_x$, $t\left)\right]$ $=$ 0 on the location parameter ${\mathit{\mu }}_x$ to maximize the log likelihood, finding the condition $$\begin{array}{cc}\sum & {\displaystyle \frac{\partial }{\partial {\mathit{\mu }}_x}}\log \left[p_{GY},\left(x,,,{\mathit{\mu }}_x,,,t\right)\right]\\ & & =\sum {\displaystyle \frac{\partial }{\partial {\mathit{\mu }}_x}}\log \left[p_{GX},\left(x,,,{\mathit{\mu }}_x\right)\right]=0.\end{array}$$ Therefore, if $x$ is stationary, the condition of maximum likelihood for $p_{GX}\left(x\right)$ coincides with the condition of maximum likelihood for $p_{GY}(y$, $t)$, and applying MLE for fitting the stationary parameters (${\mathit{\sigma }}_x$, ${\mathit{\mu }}_x$) is equivalent to fitting the parameters ($a$, $b$) of by Eqs. (10)–(12) by non-stationary MLE. The equivalence between the two methodologies suggests that the TS approach is dual to the EM approach, meaning that any implementation of EM is equivalent to an implementation of the TS approach for some transformation $f(y$, $t)$ $:$ $y\left(t\right)$ $\to$ $x\left(t\right)$ (see Appendix A for a more detailed discussion). One can also prove that Eq. (1) allows a general TS formulation with a constant shape parameter, i.e., all the TS models with a constant ${\mathit{\epsilon }}_y$ can be connected to Eq. (1) (see Appendix A). This last result is remarkable, because it shows that Eq. (1) is exhaustive for all the TS models with a constant shape parameter.

The findings drawn above are general and can be applied also to peak over threshold (POT) methodologies, because the GPD is formally derived from the GEV as the conditional probability that an observation beyond a given threshold $u$ is greater than $x$. In particular, the POT/GPD parameters are given by $$\begin{array}{}& & u_y\left(t\right)=S_y\left(t\right)\cdot u_x+T_y\left(t\right),& & {\mathit{\epsilon }}_y={\mathit{\epsilon }}_x=\text{const}.,& & {\mathit{\sigma }}_{\text{GPD}y}\left(t\right)={\mathit{\sigma }}_y\left(t\right)+{\mathit{\epsilon }}_y\left[u_y,\left(,t,\right),-,{\mathit{\mu }}_y,\left(,t,\right)\right]=S_y\left(t\right)\cdot {\mathit{\sigma }}_{\text{GPD}x},\end{array}$$ where $u_x\left(t\right)$ and $u_y\left(t\right)$ are the thresholds of the $x$ and $y$ time series, ${\mathit{\epsilon }}_y$ $=$ ${\mathit{\epsilon }}_x$ is the shape parameter, ${\mathit{\sigma }}_{\text{GPD}x}$ and ${\mathit{\sigma }}_{\text{GPD}y}\left(t\right)$ are the GPD scale parameters of $x$, and $y$, ${\mathit{\sigma }}_y$, and ${\mathit{\mu }}_y$ are the scale and location parameters of a GEV associated with the GPD, which have been included in Eq. (19) to make it clear how the parameter ${\mathit{\sigma }}_{\text{GPD}y}\left(t\right)$ can be derived.

It is worth noting that the TS methodology is neutral for a stationary series, i.e., the application of this methodology to a stationary series leads to the same results as a stationary EVA with the same underlying statistical model. That is because in such case $T_y$ and $S_y$ are constant, and Eq. (2) reduces to a constant translation and scaling.

Modeling seasonality

Often, we would like to model extreme events that show seasonality, for example with local winter extremes that differ in magnitude from summer extremes. A simple way to add the seasonal cycle to Eqs. (7)–(9) is by expressing the trend $T_y\left(t\right)$ and the standard deviation $S_y\left(t\right)$ as $$\begin{array}{}& & T_y\left(t\right)=T_{\mathrm{0}y}\left(t\right)+s_{\text{T}}\left(t\right),& & S_y\left(t\right)=S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right),\end{array}$$ where $T_{\mathrm{0}y}\left(t\right)$ and $s_{\text{T}}\left(t\right)$ are, respectively, the long-term varying and seasonal components of the trend, $S_{\mathrm{0}y}\left(t\right)$ is the long-term varying standard deviation and $s_{\text{S}}\left(t\right)$ is the seasonality factor of the standard deviation. In the notation, the subscript 0 denotes the long-term varying components. Applying Eqs. (22)–(24) to Eq. (2), we obtain $$x\left(t\right)={\displaystyle \frac{y\left(t\right)-T_{\mathrm{0}y}\left(t\right)-s_{\text{T}}\left(t\right)}{S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right)}}.$$ The time-varying GEV parameters can be expressed as $$\begin{array}{}& & {\mathit{\epsilon }}_y={\mathit{\epsilon }}_x=\text{const}.,& & {\mathit{\sigma }}_y\left(t\right)=S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right)\cdot {\mathit{\sigma }}_x,& & {\mathit{\mu }}_y\left(t\right)=S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right)\cdot {\mathit{\mu }}_x+T_{\mathrm{0}y}\left(t\right)+s_{\text{T}}\left(t\right),\end{array}$$ and the time-varying POT/GPD parameters can be expressed as $$\begin{array}{}& & u_y\left(t\right)=S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right)\cdot u_x+T_{\mathrm{0}y}\left(t\right)+s_{\text{T}}\left(t\right),& & {\mathit{\epsilon }}_y={\mathit{\epsilon }}_x=\text{const},& & {\mathit{\sigma }}_{\text{GPD}y}\left(t\right)=S_{\mathrm{0}y}\left(t\right)\cdot s_{\text{S}}\left(t\right)\cdot {\mathit{\sigma }}_{\text{GPD}}.\end{array}$$

Implementation

The implementation of the TS methodology is illustrated in Fig. 1. The fundamental input is represented by the series itself, and the core of the implementation consists of a set of algorithms for the elaboration of the time-varying trend $T_{\mathrm{0}y}\left(t\right)$, standard deviation $S_{\mathrm{0}y}\left(t\right)$, and seasonality terms $s_{\text{T}}\left(t\right)$ and $s_{\text{S}}\left(t\right)$.

TS methodology: block diagram.

[Figure omitted. See PDF]

In this study, we propose algorithms based on running means and running statistics (see Sect. 2.2.1). Hence, an important aspect is the definition of a time window $W$ for the estimation of the long-term statistics $T_{\mathrm{0}y}\left(t\right)$ and $S_{\mathrm{0}y}\left(t\right)$ and of a time window $W_{\text{sn}}$ for the estimation of the seasonality. The computation of $T_{\mathrm{0}y}\left(t\right)$ and $S_{\mathrm{0}y}\left(t\right)$ acts as a low-pass filter removing the variability within $W$. Therefore, $W$ should be chosen short enough to incorporate in the analysis the variability above the desired timescale but long enough to exclude noise, short-term variability, and sharp variations in the statistical properties of the transformed series. For example, in studies of long-term climate changes a reasonable choice is to impose $W$ $=$ 30 years, because this is the generally accepted time horizon for observing significant variations in climate (e.g., Arguez and Vose, 2011; Hirabayashi et al., 2013). It is worth stressing that the chosen value of $W$ should be verified a posteriori to ensure that the transformed series is stationary. The time window $W_{\text{sn}}$ is used to estimate the intra-annual variability of the standard deviation (see Sect. 2.2.1). In Fig. 1, the input corresponding to the seasonal time window $W_{\text{sn}}$ is drawn in a dashed box because its value is easier to choose than the value of $W$. For the examined case studies, a value of 2 months for $W_{\text{sn}}$ always resulted in a satisfactory estimation of the seasonal cycle.

In this implementation of the TS methodology, the estimation of the long-term statistics is separated from the estimation of the seasonality. This allows to study the long-term variability of the extreme values as is typically done when studying extremes on an annual basis, as well as the combination of long-term and seasonal variability to evaluate extremes on a monthly basis.

After the estimation of $T_{\mathrm{0}y}\left(t\right)$, $S_{\mathrm{0}y}\left(t\right)$, $s_{\text{T}}\left(t\right)$, and $s_{\text{S}}\left(t\right)$ we can apply Eq. (2) and perform a stationary EVA on the transformed series. It is important to stress that the stationary EVA is performed on the whole time horizon. The stationarity of the transformed signal allows us to apply different techniques for the EVA. In this study, we illustrate the GEV and GPD approaches, but an interesting development would be the elaboration of non-stationary techniques for other approaches such as those described by Goda (1988) or Boccotti (2000), based on the TS methodology.

The final step of the implementation is the back-transformation of the fitted extreme value distribution into a non-stationary one as given by Eqs. (10)–(12) and (25)–(27) for GEV and by Eqs. (19)–(21) and (28)–(30) for GPD.

Estimation of trend, standard deviation, and seasonality

There are several possible ways of estimating the slowly varying trend and standard deviation and their seasonality. We propose here a simple methodology based on a running mean and standard deviation. We formulate the trend $T_{\mathrm{0}y}\left(t\right)$ as a running mean of the signal $y\left(t\right)$ on a multi-yearly time window $W$, $$T_{\mathrm{0}y}\left(t\right)=\sum _{tt=t-W/\mathrm{2}}^{tt=t+W/\mathrm{2}}y\left(tt\right)/N_t,$$ where $N_t$ is the number of observations available during the time interval [$t$ $-$ $W/\mathrm{2}$, $t$ $+$ $W/\mathrm{2}$]. The seasonality of the trend relative to a given month of the year can be estimated as the average monthly anomaly of the de-trended series. For a given month of the year the seasonality is then $$s_{\text{T}}\left(\text{month}\right[t\left]\right)=\sum _{\text{years}}{\displaystyle \frac{\left[y,\left(,t,t,\right),-,T_{\mathrm{0}y},\left(,t,t,\right)\right]{|}_{tt\in \text{month}\left(t\right)}}{N_{\text{month}}}},$$ where the subscript $tt$ $\in$ month[$t$] indicates that the averaging operation is limited to time intervals within each considered month of the year. For example, the seasonality of January is computed as the average for all months of January of the detrended signal. To estimate the slowly varying standard deviation, we execute a running standard deviation with the same time window used to estimate $T_{\mathrm{0}y}\left(t\right)$: $$\begin{array}{cc}{S}_{\mathrm{0}y}& \left(t\right){|}_{\text{ROUGH}}=\sum _{tt=t-W/\mathrm{2}}^{tt=t+W/\mathrm{2}}\\ & & \sqrt{\left[y\right(tt)-\stackrel{\mathrm{‾}}{y}(tt\in [t-W/\mathrm{2},t+W/\mathrm{2}]){]}^{\mathrm{2}}/N_{\text{Wsn}}},\end{array}$$ where the subscript ROUGH stresses the fact that this expression is sensitive to outliers and that its direct employment leads to a relevant statistical error, as explained in Sect. 2.2.2. To overcome this problem, we smooth $S_{\mathrm{0}y}\left(t\right){|}_{\text{ROUGH}}$ with a moving average on a time window smaller than $W$, for example $W/L$ with $L$ $=$ 2: $$S_{\mathrm{0}y}\left(t\right)=\sum _{tt=t-W/\mathrm{2}L}^{tt=t+W/\mathrm{2}L}L{S}_{\mathrm{0}y}\left(tt\right){|}_{\text{ROUGH}}/N_t.$$ It is worth stressing that, in general, a further smoothing of the results of running means and standard deviations is appropriate if it reduces the error and improves the detection of the slowly varying statistical behavior of the time series. This is because the estimation of $T_{\mathrm{0}y}\left(t\right)$ and $S_{\mathrm{0}y}\left(t\right)$ involves a low-pass filter to smooth the signal on timescales lower than $W$ and remove high-frequency variability.

To estimate the seasonality we perform another running standard deviation $S_{\text{sn}}\left(t\right)$ on a time window $W_{\text{sn}}$ much shorter than 1 year, in the order of the month: $$\begin{array}{cc}{S}_{\text{sn}}& \left(t\right)=\sum _{tt=t-W_{\text{sn}}/\mathrm{2}}^{tt=t+W_{\text{sn}}/\mathrm{2}}\\ & & \sqrt{{\left[y,\left(,t,t,\right),-,\stackrel{\mathrm{‾}}{y},\left(t,t,\in ,\left[t,-,W_{\text{sn}},/,\mathrm{2},,,t,+,W_{\text{sn}},/,\mathrm{2}\right]\right)\right]}^{\mathrm{2}}/N_t}.\end{array}$$ The seasonality of the standard deviation can then be computed as the monthly average of the ratio between $S_{\text{sn}}\left(t\right)$ and $S_{\mathrm{0}y}\left(t\right)$: $$s_{\text{S}}\left(\text{month}\right[t\left]\right)=\sum _{\text{years}}{\displaystyle \frac{\left[S_{\text{sn}},\left(,t,t,\right),/,S_{\mathrm{0}y},\left(,t,t,\right)\right]{|}_{tt\in \text{month}\left(t\right)}}{N_{tt\in \text{month}\left(t\right)}}}.$$ The estimated seasonality terms $s_{\text{T}}$ and $s_{\text{S}}$ are periodic with a period of 1 year. In order to smooth them and remove any possible noise in the signal, we take into account only their first three Fourier components computed in a period of 1 year, corresponding to components with a periodicity of 1 year, 6 months, and 3 months.

Statistical error

Since there is an inherent error in the estimation of the trend, standard deviation and seasonality given by Eqs. (32)–(36), we need to estimate this error and propagate it to the statistical error of the parameters of the non-stationary GEV and GPD distributions. In general, given a sample $d$ of data with size $N$, average $\stackrel{\mathrm{‾}}{s}$, variance var($s$), and standard deviation $S\left(s\right)$ we have $$\begin{array}{}& & \text{var}\left(\stackrel{\mathrm{‾}}{d}\right)=\text{var}\left(\stackrel{\mathrm{‾}}{d}\right)/N\Rightarrow \text{Err}\left[\stackrel{\mathrm{‾}}{d}\right]=S\left(d\right)/\sqrt{N},& & \mathrm{\text{var}}\left[\text{var}\right(d\left)\right]\approx \mathrm{2}\text{var}(d{)}^{\mathrm{2}}/N\Rightarrow \text{Err}[S\left(d\right)]\approx S(d)\cdot \sqrt[\mathrm{4}]{\mathrm{2}/N}.\end{array}$$ Equation (37) represents the error on the average and can be obtained by propagating the intrinsic error of each observation, given by the standard deviation $S\left(s\right)$, to expression $\stackrel{\mathrm{‾}}{s}$ $=$ $\sum s_i/N$. Equation (38) represents the error on the standard deviation and can be evaluated considering that with a Gaussian approximation quantity $S$ $=$ $\sum _N{s}_i^{\mathrm{2}}/$var($s$) follows a chi-squared distribution with standard deviation 2$N$.

Using Eqs. (37) and (38), we can estimate the error on $T_{\mathrm{0}y}\left(t\right)$ and $S_{\mathrm{0}y}\left(t\right){|}_{\text{ROUGH}}$ as $$\begin{array}{}& & \text{Err}\left[T_{\mathrm{0}y}\right]\approx S_{\mathrm{0}y}/\sqrt{N_t},& & \text{Err}\left[S_{\mathrm{0}y}\right]{|}_{\text{ROUGH}}\approx S_{\mathrm{0}y}\cdot \sqrt[\mathrm{4}]{\mathrm{2}/N_t}.\end{array}$$ As mentioned in Sect. 2.2.1, Eq. (40) tends to return rather high values of the error relative to $S_{\mathrm{0}y}\left(t\right)$. For example, if we are considering a time window of 20 years with an observation every 3 h, we have $$N_t\approx 59000\Rightarrow {\displaystyle \frac{\text{Err}\left[S_{\mathrm{0}y}\right]{|}_{\text{ROUGH}}}{S_{\mathrm{0}y}}}\approx 7.6\mathit{\%}.$$ Using expression Eq. (34) for the estimation of $S_{\mathrm{0}y}\left(t\right)$ overcomes this issue because we can estimate the uncertainty in $S_{\mathrm{0}y}\left(t\right)$ as the error of the standard deviation averaged over the time window $W/L$, which is significantly lower than the error given by Eq. (41). Using Eq. (37), we find $$E\left[S_{\mathrm{0}y}\right]\approx {\displaystyle \frac{\text{Err}\left[S_{\mathrm{0}y}\right]{|}_{\text{ROUGH}}}{\sqrt{N_t/L}}}=S_{\mathrm{0}y}\cdot \sqrt[\mathrm{4}]{{\displaystyle \frac{\mathrm{2}{L}^{\mathrm{2}}}{N_t^{\mathrm{3}}}}}.$$ We can estimate the error on the seasonality of the trend $s_{\text{T}}$ by adding the error estimated for $T_{\mathrm{0}y}\left(t\right)$ to that of the monthly mean. As the statistical error of independent Gaussian variables sum vectorially, we obtain $$\text{Err}\left[s_{\text{T}}\right]=\sqrt{{\text{Err}}^{\mathrm{2}}\left[\text{mntmean}\right(y\left)\right]+{\text{Err}}^{\mathrm{2}}\left[T_{\mathrm{0}y}\right]},$$ where the mntmean($y$) operator represents the monthly average of $y$. If, for example, one considers the month of January, it is the average computed on all months of January in the time series. Assuming the error on mntmean($y$) as approximately constant within the year, it follows that $$\text{Err}\left[\text{mntmean}\right(y\left)\right]\approx S_{\mathrm{0}y}/\sqrt{N_{\text{month}}}\approx S_{\mathrm{0}y}\cdot \sqrt{12/N_{\text{tot}}},$$ where $N_{\text{month}}$ is the number of observations corresponding to the considered month, $N_{\text{tot}}$ is the total number of elements of the series $y\left(t\right)$, $N_{\text{month}}$ $\approx$ $N_{\text{tot}}/12$. Therefore, Eq. (43) can be rewritten as $$\text{Err}\left[s_{\text{T}}\right]\approx S_{\mathrm{0}y}\sqrt{12/N_{\text{tot}}+\mathrm{1}/N_t}.$$ The error on $s_{\text{S}}$ can be estimated as the error of the average ratio $s_{\text{S}}/S_{\mathrm{0}y}$. Using Eq. (38), the error of the ratio $s_{\text{S}}/S_{\mathrm{0}y}$ is given by $$\begin{array}{cc}\text{Err}& \left[{\displaystyle \frac{S_{\text{sn}}}{S_{\mathrm{0}y}}}\right]\approx \sqrt{{\left({\displaystyle \frac{\text{Err}\left[S_{\text{sn}}\right]}{S_{\mathrm{0}y}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{S_{\text{sn}}}{S_{\mathrm{0}y}^{\mathrm{2}}}},\text{Err},\left[S_{\mathrm{0}y}\right]\right)}^{\mathrm{2}}}\\ & & \approx {\displaystyle \frac{S_{\text{sn}}}{S_{\mathrm{0}y}}}\sqrt{\sqrt{{\displaystyle \frac{\mathrm{2}}{N_{\text{sn}}}}}+\sqrt{{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{N_{\text{t}}^{\mathrm{3}}}}}}\approx s_{\text{S}}\sqrt[\mathrm{4}]{{\displaystyle \frac{\mathrm{2}}{N_{\text{sn}}}}},\end{array}$$ where $N_{\text{sn}}$ is the average number of observations within the time window $W_{\text{sn}}$ and assuming $N_t$ $\gg$ $N_{\text{sn}}$. We can then estimate the error on $s_{\text{S}}$ as the error of the monthly average of $s_{\text{S}}/S_{\mathrm{0}y}$: $$\begin{array}{cc}\text{Err}\left[s_{\text{S}}\right]& \approx \text{Err}\left[{\displaystyle \frac{S_{\text{sn}}}{S_{\mathrm{0}y}}}\right]/\sqrt{N_{\text{month}}}\approx s_{\text{S}}\sqrt{{\displaystyle \frac{12}{N_{\text{tot}}}}}\sqrt[\mathrm{4}]{{\displaystyle \frac{\mathrm{2}}{N_{\text{sn}}}}}\\ & & =s_{\text{S}}\sqrt[\mathrm{4}]{{\displaystyle \frac{288}{N_{\text{tot}}^{\mathrm{2}}{N}_{\text{sn}}}}}.\end{array}$$ Using Eqs. (40), (45), and (47) we can estimate the error on the time-varying GEV parameters as $$\begin{array}{}& & \text{Err}\left[{\mathit{\epsilon }}_y\right]=\text{Err}\left[{\mathit{\epsilon }}_x\right],& & \mathrm{\text{Err}}\left[{\mathit{\sigma }}_y\right]=\sqrt{{\left(S_{\mathrm{0}y},\cdot ,s_{\text{S}},\cdot ,\text{Err},\left[{\mathit{\sigma }}_x\right]\right)}^{\mathrm{2}}+{\left(S_{\mathrm{0}y},\cdot ,\text{Err},\left[s_{\text{S}}\right],\cdot ,{\mathit{\sigma }}_x\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_{\mathrm{0}y}\right],\cdot ,s_{\text{S}},\cdot ,{\mathit{\sigma }}_x\right)}^{\mathrm{2}}},& & \mathrm{\text{Err}}\left[{\mathit{\mu }}_y\right]=\sqrt{{\left(S_{\mathrm{0}y},\cdot ,s_{\text{S}},\cdot ,\text{Err},\left[{\mathit{\mu }}_x\right]\right)}^{\mathrm{2}}+{\left(S_{\mathrm{0}y},\cdot ,\text{Err},\left[s_{\text{S}}\right],\cdot ,{\mathit{\mu }}_x\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_{\mathrm{0}y}\right],\cdot ,s_{\text{S}},\cdot ,{\mathit{\mu }}_x\right)}^{\mathrm{2}}+{\text{Err}}_{\text{T}}^{\mathrm{2}}},\end{array}$$ and the error on the time-varying GPD parameters as $$\begin{array}{}& & \mathrm{\text{Err}}\left[u_y\right]=\sqrt{{\left(S_{\mathrm{0}y},\cdot ,s_{\text{S}},\cdot ,\text{Err},\left[u_x\right]\right)}^{\mathrm{2}}+{\left(S_{\mathrm{0}y},\cdot ,\text{Err},\left[s_{\text{S}}\right],\cdot ,u_x\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_{\mathrm{0}y}\right],\cdot ,s_{\text{S}},\cdot ,u_x\right)}^{\mathrm{2}}+{\text{Err}}_{\text{T}}^{\mathrm{2}}},& & \text{Err}\left[{\mathit{\epsilon }}_y\right]=\text{Err}\left[{\mathit{\epsilon }}_x\right],& & \mathrm{\text{Err}}\left[{\mathit{\sigma }}_{\text{GPD}y}\right]=\sqrt{{\left(S_{\mathrm{0}y},\cdot ,s_{\text{S}},\cdot ,\text{Err},\left[{\mathit{\sigma }}_{\text{GPD}x}\right]\right)}^{\mathrm{2}}+{\left(S_{\mathrm{0}y},\cdot ,\text{Err},\left[s_{\text{S}}\right],\cdot ,{\mathit{\sigma }}_{\text{GPD}x}\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_{\mathrm{0}y}\right],\cdot ,s_{\text{S}},\cdot ,{\mathit{\sigma }}_{\text{GPD}x}\right)}^{\mathrm{2}}},\end{array}$$ where $${\text{Err}}_{\text{T}}^{\mathrm{2}}={\text{Err}}^{\mathrm{2}}\left[T_{\mathrm{0}y}\right]+{\text{Err}}^{\mathrm{2}}\left[s_{\text{T}}\right].$$

Data and validation

To assess the generality of the approach, the TS methodology has been validated on time series of different variables, from different sources and with different statistical properties.

The analysis of annual and monthly maxima has been carried out on time series of significant wave height at two locations: the first located in the Atlantic Ocean, west of Ireland ($-$10.533${}^{\circ }$ E, 55.366${}^{\circ }$ N), and the second close to Cape Horn (60.237${}^{\circ }$ E, $-$57.397${}^{\circ }$ N). The data have been obtained by means of wave simulations performed with the spectral model Wavewatch III^® (Tolman, 2014) forced by the wind data projections of the RCP8.5 scenario (van Vuuren et al., 2011) of the CMIP5 model GFDL-ESM2M (Dunne et al., 2012) on a time horizon spanning from 1970 to 2100. This data set is referred to from now on as GWWIII. Here, the TS methodology is used in order to examine its applicability to climate change studies. The annual and monthly analyses have been repeated on a series of water level residuals offshore of the Hebrides Islands (Scotland, $-$7.9${}^{\circ }$ E, 57.3${}^{\circ }$ N) obtained from a 35-year hindcast of storm surges at European scale (Vousdoukas et al., 2016a, b) forced by the ERA-Interim reanalysis data (Dee et al., 2011). This data set is further referred to as JRCSURGES.

For the annual maxima of the considered series, we further compare the TS methodology with the SS technique as implemented by Alfieri et al. (2015) and Vousdoukas et al. (2016a). For this purpose, we extracted time series from projections of streamflow in the Rhine and Po rivers covering a time horizon from 1970 to 2100 (Alfieri et al., 2015), from now on referred to as JRCRIVER. Also, the two series of significant wave height of west Ireland and Cape Horn extracted from the GWWIII data set have been used in this comparison.

Finally, we compare the TS methodology and the EM for monthly maxima using time series of significant wave height extracted from a 35-year wave hindcast database (Mentaschi et al., 2015) near the locations of La Spezia and Ortona. The analysis of this data set, further referred to as WWIII_MED, focuses on a comparison between seasonal cycles modeled by the two approaches.

Results

Waves: annual extremes

The validation of the TS methodology was performed first on the time series of significant wave height of west Ireland and Cape Horn from the GWWIII data set. We verified first the non-seasonal transformation given by Eq. (2) and the time-dependent GEV and GPD given by Eqs. (7)–(9) and (19)–(21), respectively. By ignoring the seasonality, this formulation is suitable for finding extremes and peaks on an annual basis. For technical reasons the two series do not have data in two time intervals, from 2005 to 2010 and from 2092 to 2095. The impact of the missing data on the analysis is small, however, especially if we choose a time window $W$ large enough for the estimation of the trend and standard deviation using Eqs. (31) and (33). In particular, for this analysis we chose a time window of 20 years, which is long enough to ensure the accuracy of the results and short enough to include the multi-decadal variability of a 130-year time series.

The results of the analysis for the two time series are illustrated in Figs. 2 and 3. Panel (a) of each figure shows the original time series and its slowly varying trend and standard deviation. Panel (b) illustrates the normalized series obtained through the transformation given by Eq. (1), allowing an evaluation at a glance of the stationarity of the normalized series. The mean and the standard deviation of the normalized series plotted in panel (b) are 0 and 1, respectively. Higher order statistics such as skewness and kurtosis are included in the graphics to support the assumption of stationarity of the normalized series. From the normalized time series we extracted the annual maxima and estimated the corresponding non-stationary GEV as given by Eqs. (7)–(9) (see Figs. 2c and 3c). Moreover, we performed a POT selection of the extreme events on the normalized series. The threshold was defined in order to have on average five events per year, following Ruggiero et al. (2010), corresponding for both of the series to the 97th percentile. From the resultant POT sample we estimated the corresponding non-stationary GPD as given by Eqs. (19)–(21) (see Figs. 2d and 3d). In Fig. 2c and d and Fig. 3c and d, the shape parameters $\mathit{\epsilon }$ estimated by the MLE for the GEV and the GPD are also reported. Inter-decadal oscillations in the annual maxima are modeled for both of the series, though they are more pronounced for the west Ireland time series. Moreover, for both series there is a tendency for the annual maxima to increase. This is more pronounced for the Cape Horn series, where the increase in the annual maxima of significant wave height estimated by GWWIII is about 2 m.

Long-term analysis of the projections of significant wave height in Cape Horn: (a) series, its trend, and standard deviation; (b) the normalized series with higher order statistical indicators; (c) non-stationary GEV of annual maxima; (d) non-stationary GPD of annual peaks. In (c) and (d) the values of the shape parameter $\mathit{\epsilon }$ best fitted for the GEV and GPD distributions are reported.

[Figure omitted. See PDF]

Long-term analysis of the projections of significant wave height in Cape Horn: (a) series, its trend, and standard deviation; (b) the normalized series with higher order statistical indicators; (c) non-stationary GEV of annual maxima; (d) non-stationary GPD of annual peaks. In (c) and (d) the values of the shape parameter $\mathit{\epsilon }$ best fitted for the GEV and GPD distributions are reported.

[Figure omitted. See PDF]

It is worth noting that for both of the considered series, the statistical mode of GEV and GPD grows faster in time than the slowly varying trend $T_y\left(t\right)$. This is due to the fact that the growth of the location parameter ${\mathit{\mu }}_y\left(t\right)$ of the non-stationary GEV (Eq. 7) and of the threshold $u_y\left(t\right)$ of the non-stationary GPD (Eq. 19) are related not only to the growth of $T_y\left(t\right)$ but also to the growth of $S_y\left(t\right)$. The upper tail of the distributions grows even faster because the scale parameter is also proportional to $S_y\left(t\right)$.

Average error components for the long-term analysis of the projections of significant wave height extracted at west Ireland and Cape Horn, for non-stationary GEV and GPD. The error is dominated by the component due to the stationary MLE.

Yearly maxima: trend-only analysis
Error	West	Cape Horn	%
component	Ireland	error (m)	(err${}^{\mathrm{2}}$)
(average)	error (m)
Non-stationary GEV
$S_y$ $\cdot$ Err$\left[{\mathit{\sigma }}_x\right]$	0.0371	0.0372	100 %
Err$\left[S_y\right]$ $\cdot$ ${\mathit{\sigma }}_x$	5.876 $\times$ 10${}^{-\mathrm{4}}$	5.818 $\times$ 10${}^{-\mathrm{4}}$	$<$ 0.1 %
Err$\left[{\mathit{\sigma }}_y\right]$	0.0371	0.0372	100 %
$S_y$ $\cdot$ Err$\left[{\mathit{\mu }}_x\right]$	0.0538	0.0536	97.7 %
Err$\left[S_y\right]$ $\cdot$ ${\mathit{\mu }}_x$	3.6 $\times$ 10${}^{-\mathrm{3}}$	3.4 $\times$ 10${}^{-\mathrm{3}}$	0.4 %
Err$\left[T_y\right]$	7.4 $\times$ 10${}^{-\mathrm{3}}$	7.0 $\times$ 10${}^{-\mathrm{3}}$	1.85 %
Err$\left[{\mathit{\mu }}_y\right]$	0.053	0.054	100 %
Non-stationary GPD
$S_y$ $\cdot$ Err$\left[{\mathit{\sigma }}_{\text{GPD}x}\right]$	0.0418	0.0310	100 %
Err$\left[S_y\right]$ $\cdot$ ${\mathit{\sigma }}_{\text{GPD}x}$	1.12 $\times$ 10${}^{-\mathrm{3}}$	8.9 $\times$ 10${}^{-\mathrm{4}}$	$<$ 0.1 %
Err$\left[{\mathit{\sigma }}_{\text{GPD}y}\right]$	0.0418	0.0310	100 %
$S_y$ $\cdot$ Err$\left[u_x\right]$	0.1489	0.1376	100 %
Err$\left[S_y\right]$ $\cdot$ $u_x$	1.9 $\times$ 10${}^{-\mathrm{3}}$	1.7 $\times$ 10${}^{-\mathrm{3}}$	$<$ 0.1 %
Err$\left[u_y\right]$	0.1491	0.1278	100 %

Average error components for the seasonal analysis of the projections of significant wave height extracted at west Ireland and Cape Horn, for non-stationary GEV and GPD. The error is dominated by the component due to the stationary MLE.

Monthly maxima: seasonal analysis
Error	West Ireland	Cape Horn	%
component	error (m)	error (m)	(err${}^{\mathrm{2}}$)
(average)
Non-stationary GEV
$S_{\mathrm{0}y}$ $\cdot$ $s_{\text{S}}$ $\cdot$ Err$\left[{\mathit{\sigma }}_x\right]$	0.0135	0.0138	99.7 %
$S_{\mathrm{0}y}$ $\cdot$ Err$\left[s_{\text{S}}\right]$ $\cdot$ ${\mathit{\sigma }}_x$	7.2 $\times$ 10${}^{-\mathrm{4}}$	7.6 $\times$ 10${}^{-\mathrm{4}}$	0.3 %
Err$\left[{\mathit{\sigma }}_y\right]$	0.0135	0.0138	100 %
$S_{\mathrm{0}y}$ $\cdot$ $s_{\text{S}}$ $\cdot$ Err$\left[{\mathit{\mu }}_x\right]$	0.019	0.020	96.6 %
$S_{\mathrm{0}y}$ $\cdot$ Err$\left[s_{\text{S}}\right]$ $\cdot$ ${\mathit{\mu }}_x$	0.0014	0.0017	0.7 %
Err$\left[s_{\text{T}}\right]$	4.86 $\times$ 10${}^{-\mathrm{6}}$	5.25 $\times$ 10${}^{-\mathrm{6}}$	$<$ 0.1 %
Err$\left[{\mathit{\mu }}_y\right]$	0.0204	0.0214	100 %
Non-stationary GPD
$S_{\mathrm{0}y}$ $\cdot$ $s_{\text{S}}$ $\cdot$ Err$\left[{\mathit{\sigma }}_{\text{GPD}x}\right]$	0.025	0.029	100 %
$S_{\mathrm{0}y}$ $\cdot$ Err$\left[s_{\text{S}}\right]$ $\cdot$ ${\mathit{\sigma }}_{\text{GPD}x}$	9.4 $\times$ 10${}^{-\mathrm{4}}$	9.9 $\times$ 10${}^{-\mathrm{4}}$	$<$ 0.1 %
Err$\left[{\mathit{\sigma }}_{\text{GPD}y}\right]$	0.0253	0.0293	100 %
$S_{\mathrm{0}y}$ $\cdot$ $s_{\text{S}}$ $\cdot$ Err$\left[u_x\right]$	0.1061	0.1205	100 %
$S_{\mathrm{0}y}$ $\cdot$ Err$\left[s_{\text{S}}\right]$ $\cdot$ $u_x$	0.0011	0.0014	$<$ 0.1 %
Err$\left[u_y\right]$	0.1063	0.1207	100 %

Seasonal analysis of the projections of significant wave height in west Ireland: (a) series, its trend, and standard deviation; (b) the normalized series with higher order statistical indicators; (c) non-stationary GEV of annual maxima; (d) non-stationary GPD of annual peaks. In (c) and (d) the values of the shape parameter $\mathit{\epsilon }$ best fitted for the GEV and GPD distributions are reported. For the sake of clarity, only a 5-year time slice is reported.

[Figure omitted. See PDF]

Seasonal analysis of the projections of significant wave height in Cape Horn: (a) series, its trend, and standard deviation; (b) the normalized series with higher order statistical indicators; (c) non-stationary GEV of annual maxima; (d) non-stationary GPD of annual peaks. In (c) and (d) the values of the shape parameter $\mathit{\epsilon }$ best fitted for the GEV and GPD distributions are reported. For the sake of clarity, only a 5-year time slice is reported.

[Figure omitted. See PDF]

The impact of the statistical error in the slowly varying trend and the standard deviation on the uncertainty of the distribution parameters have been examined using Eqs. (48)–(50) and (51)–(53), which, for the non-seasonal analysis, reduce to $$\begin{array}{}& & \text{Err}\left[{\mathit{\epsilon }}_y\right]=\text{Err}\left[{\mathit{\epsilon }}_x\right],& & \text{Err}\left[{\mathit{\sigma }}_y\right]=\sqrt{{\left(S_y,\cdot ,\text{Err},\left[{\mathit{\sigma }}_x\right]\right)}^{\mathrm{2}}+{\left(\text{Err},\left[s_y\right],\cdot ,{\mathit{\sigma }}_x\right)}^{\mathrm{2}}},& & \mathrm{\text{Err}}\left[{\mathit{\mu }}_y\right]=\sqrt{{\left(S_y,\cdot ,\text{Err},\left[{\mathit{\mu }}_x\right]\right)}^{\mathrm{2}}+{\left(\text{Err},\left[T_y\right],\cdot ,{\mathit{\mu }}_x\right)}^{\mathrm{2}}+{\text{Err}}^{\mathrm{2}}\left[T_y\right]},\end{array}$$ for the GEV, and to $$\begin{array}{}& & \mathrm{\text{Err}}\left[u_y\right]=\sqrt{{\left(S_y,\cdot ,\text{Err},\left[u_x\right]\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_y\right],\cdot ,u_x\right)}^{\mathrm{2}}+{\text{Err}}^{\mathrm{2}}\left[T_y\right]},& & \text{Err}\left[{\mathit{\epsilon }}_y\right]=\text{Err}\left[{\mathit{\epsilon }}_x\right],& & \mathrm{\text{Err}}\left[{\mathit{\sigma }}_{\text{GPD}y}\right]=\sqrt{{\left(S_y,\cdot ,\text{Err},\left[{\mathit{\sigma }}_{\text{GPD}x}\right]\right)}^{\mathrm{2}}+{\left(\text{Err},\left[S_y\right],\cdot ,{\mathit{\sigma }}_{\text{GPD}x}\right)}^{\mathrm{2}}},\end{array}$$ for the GPD. The result is that for the non-seasonal analysis the error due to the estimation of the trend and standard deviation is negligible with respect to the error associated with the stationary MLE. In Table 1, the values of the different components of the compared error in Eqs. (55)–(57) and (58)–(60) are reported together with the total error estimated for each parameter of the non-stationary GEV and GPD. Since the threshold $u_x$ of the stationary GPD was selected to have on average five events per year, the error has been computed as the uncertainty related to this definition. The percentage contribution to the squared error is also reported in Table 1 in a single column because the percentages estimated for the two series are roughly equal. The error for both GEV and GPD and for the two series is clearly dominated by the error associated with the estimation of the parameters of the stationary distributions ([$S_y$ $\cdot$ Err[${\mathit{\sigma }}_x$]] and [$S_y$ $\cdot$ Err[${\mathit{\mu }}_x$]] for the GEV and [$S_y$ $\cdot$ Err[${\mathit{\sigma }}_{\text{GPD}x}$]] and [$S_y$ $\cdot$ Err[$u_x]$] for the GPD).

Waves: monthly extremes

The seasonal formulation of the approach is suitable to estimate extreme value distributions on a monthly basis. Hence, we applied Eq. (24) to estimate the normalized series, then fitted a stationary GEV of monthly maxima by means of a MLE that was back-transformed into a non-stationary GEV through Eqs. (25)–(27). It is worth stressing that for the stationary MLE, the entire normalized series was used, covering a time horizon of 130 years. For the GPD, we selected the threshold in order to have on average 12 events per year, corresponding to the 93rd percentile for both series. Results are displayed in Fig. 4 for the location of west Ireland and in Fig. 5 for Cape Horn. To make the seasonal cycle distinguishable in these figures, we plotted only a slice of 5 years from 2085 to 2090. The meaning of the four panels in Figs. 4 and 5 is the same as in Figs. 2 and 3. The non-stationary extreme value distribution estimated for the location of west Ireland presents a strong seasonal cycle with higher and more broad-banded extremes during winter. For Cape Horn, the seasonal cycle is weaker, with the extremes of significant wave height slightly lower during the local summer. The estimated PDF for the seasonal GEV and GPD are significantly lower than those estimated for the non-seasonal analysis because in the seasonal analysis we consider monthly extremes, while in the non-seasonal one we consider annual extremes.

It is worth stressing that in the study of the monthly maxima, the long-term trend is also estimated even if it cannot be appreciated in Figs. 4 and 5 due to the short time horizon represented.

Table 2 reports the components of the statistical error due to the uncertainty in the estimation of the seasonality, together with the components of the stationary MLE. The error components relating to the uncertainty in the estimation of $T_{\mathrm{0}y}$ and $S_{\mathrm{0}y}$ were omitted as they are negligible compared with the error associated with the fitting of the stationary extreme value distribution (see Sect. 3.1). In Table 2, we can see that, as for the non-seasonal analysis, the error for both GEV and GPD and for the two series is clearly dominated by the uncertainty associated with the estimation of the parameters of the stationary distributions, though in this case the error related to the stationary MLE is significantly smaller than that found for the non-seasonal analysis due to the larger sample of data.

Residual water levels

To verify the performance of the TS methodology on a series from a different source, of a different size, and with different statistical characteristics, we tested it on a series of water level residuals extracted from the JRCSURGES data set for an off-shore location of the Hebrides Islands, Scotland ($-$7.9${}^{\circ }$ E, 57.3${}^{\circ }$ N). This series is characterized by a flat trend $T_y\left(t\right)$ because the model results are approximately constant averaged. Therefore, almost all the variability is modeled by the TS methodology in the standard deviation $S_y\left(t\right)$. Since the time horizon of this series is shorter than that of the GWWIII projections, a time window of 6 years was adopted for the computation of the trend to better identify its inter-annual variability. The results of the TS analysis of the yearly maxima are shown in Fig. 6. The series displays also a strong seasonal behavior with annual maxima usually occurring during the local winter (for brevity, the seasonal analysis is not illustrated).

An interesting aspect is that the estimated standard deviation $S_y\left(t\right)$ presents a strong correlation ($\mathit{\rho }$ $=$ 0.79) with the annual means of the North Atlantic Oscillation (NAO) index. This is illustrated in Fig. 7, where the scatter plot of $S_y\left(t\right)$ vs. the annual means of the NAO index (Fig. 7a) and the two time series (Fig. 7b) are represented. As a consequence, the estimated annual maxima are also correlated with the NAO index.

Comparison with other approaches

Stationary methodology on time slices for long trend estimation

A comparison was carried out between the TS methodology and the SS technique, consisting of a stationary analysis on quasi-stationary slices of data. This analysis was carried out on river discharge projections for the Po and the Rhine extracted from the JRCRIVER data set and on the projections of significant wave height extracted from the GWWIII data set for the locations of west Ireland and Cape Horn. The TS methodology was applied with a time window of 30 years to estimate a non-stationary GPD of annual maxima. The SS technique was carried out using a GPD approach on time slices of 30 years from 1970 to 2000, 2020 to 2050, and 2070 to 2100. For both methodologies, the threshold was selected to have on average five peaks per year.

Long-term analysis of the residual water levels modeled at the Hebrides Islands: (a) series, its trend, and standard deviation; (b) the normalized series with higher order statistical indicators; (c) non-stationary GEV of annual maxima; (d) non-stationary GPD of annual peaks. In (c) and (d) the values of the shape parameter $\mathit{\epsilon }$ best fitted for the GEV and GPD distributions are reported.

[Figure omitted. See PDF]

Time-varying standard deviation $S_y\left(t\right)$ estimated by means of the TS methodology vs. the yearly average of the NAO index, indicated by scatter plot (a) and time series (b).

[Figure omitted. See PDF]

Return level plots for the discharge of the Rhine River at its mouth, TS (black continuous line), 95 % confidence interval for the TS methodology (green band), and SS methodology (black dashed line) for the time slices 1970–2000, 2020–2050, and 2070–2100.

[Figure omitted. See PDF]

Return levels modeled by the TS methodology ($x$ axis) vs. those modeled by the SS methodology ($y$ axis) for the discharge of the Rhine and Po rivers and the significant wave height in west Ireland and Cape Horn. The three series of dots represent the three time slices. Dots' color represents the return period. The blue lines represent the maximum 30-year return level.

[Figure omitted. See PDF]

Long-term variations of the extremes of projected river discharge for the Rhine and Po rivers, and of projected significant wave height for west Ireland and Cape Horn: normalized bias (NBI) and maximum difference (max diff) between the return levels estimated with the TS methodology and the SS approach, and mean 95 % confidence interval amplitude expressed as percentage of the return level, for return periods of 5, 10, 30, 100, and 300 years.

Return period		5 years	10 years	30 years	100 years	300 years
Rhine	NBI	$-$1.07 %	$-$1.51 %	$-$2.35 %	$-$3.43 %	$-$4.53 %
(River dis.)	Max diff	$-$3.58 %	$-$4.40 %	$-$5.92 %	$-$7.81 %	$-$9.69 %
Mean conf. int. (TS)		4.90 %	5.54 %	6.68 %	8.01 %	9.27 %
Mean conf. int. (SS)		17.99 %	21.34 %	26.87 %	33.16 %	39.04 %
Po	NBI	1.47 %	2.06 %	2.92 %	3.69 %	4.25 %
(River dis.)	Max diff	5.87 %	4.88 %	5.60 %	9.57 %	13.06 %
Mean conf. int. (TS)		5.08 %	5.77 %	7.00 %	8.46 %	9.84 %
Mean conf. int. (SS)		16.77 %	20.07 %	25.45 %	31.47 %	36.99 %
W. Ireland	NBI	$-$0.28 %	$-$0.14 %	0.07 %	0.27 %	0.43 %
(Waves Hs)	Max diff	$-$0.91 %	$-$1.14 %	$-$1.48 %	2.06 %	2.51 %
Mean conf. int. (TS)		1.97 %	2.22 %	2.63 %	3.05 %	3.41 %
Mean conf. int. (SS)		7.73 %	9.01 %	10.95 %	12.91 %	14.54 %
Cape Horn	NBI	$-$1.07 %	$-$1.13 %	$-$1.17 %	$-$1.18 %	$-$1.18 %
(Waves Hs)	Max diff	$-$1.87 %	$-$2.36 %	$-$3.12 %	$-$3.92 %	$-$4.59 %
Mean conf. int. (TS)		1.74 %	2.03 %	2.52 %	3.07 %	3.57 %
Mean conf. int. (SS)		6.40 %	7.70 %	9.80 %	12.09 %	14.15 %

Results are illustrated in Fig. 8, where the return levels of the projected discharge of the Rhine are shown for three time slices. In Fig. 8 the continuous black line and the green band represent the return levels and the 95 % confidence interval estimated by the TS methodology, where the dashed black line represents the return levels estimated by the stationary EVA on the considered slice (labeled in the legend as SS). The return levels estimated for short return periods by the two methodologies are close, while they tend to spread for high return periods. This fact is also evident from Fig. 9, where the return levels estimated by the two methodologies are plotted against each other for the river discharge of the Rhine and the Po and for the significant wave height of west Ireland and Cape Horn. We can see that for the analyzed time series the two methodologies are in good agreement for return periods below 30 years while they spread for larger return periods. Some quantitative data about this fact are shown in Table 3, which reports the normalized bias NBI of the return levels of the two methodologies, defined as $$\text{NBI}=\stackrel{\mathrm{‾}}{\left({\text{RL}}_{\text{TS}},-,{\text{RL}}_{\text{cmp}}\right)}/\stackrel{\mathrm{‾}}{{\text{RL}}_{\text{cmp}}},$$ where RL${}_{\text{TS}}$ and RL${}_{\text{cmp}}$ are the return levels obtained by the TS and the SS methodology, respectively. Table 3 also includes the maximum deviation between the return levels estimated by the TS and by the SS methodology, as well as the 95 % confidence interval amplitude expressed as a percentage of the return level. The NBI and the maximum deviations were obtained by comparing results of the two techniques on the three 30-year time windows. From Table 3 we can see that the maximum deviation for return periods up to 30 years is always below 6 %, while for higher return periods it increases up to 13 % for the discharge of the Po. Moreover, the confidence intervals estimated for SS are always larger than those for TS, especially for large return periods. This is mainly due to the fact that for the stationary analysis on the quasi-stationary time slices we consider a sample of only 30 years, which leads to wider uncertainty ranges especially in the estimation of large return periods such as 100 and 300 years. This also explains the sharp variations of high return levels that we find between the three time windows using the SS approach. These variations are likely more related to the uncertainty in estimating the levels associated with long return periods rather than to climatic changes. The TS methodology allows a more accurate estimation of high return levels because it uses the whole sample of 130 years, and this represents one of the strengths of the TS methodology vs. SS. It is finally worth noting that the relative confidence interval estimated by both methodologies for the series of river discharge is larger than that estimated for the series of significant wave height. This is because for wave height data the minimum distance between two peaks has been set to at least 3 days, while for river discharge it has been set to 7 days.

Established non-stationary method for seasonal variability

Section 3 shows that the TS methodology is mathematically equivalent to a particular implementation of the EM methodology as described for example by Coles (2001), Izaguirre et al. (2011), Menéndez et al. (2009), and Sartini et al. (2015). For the sake of completeness, we show here the results of a comparison between the performances of TS and of a different formulation of the EM methodology. In its formulation, the parameters of the non-stationary GEV of the monthly maxima are expressed as $$\begin{array}{}& & \mathit{\mu }\left(t\right)={\mathit{\beta }}_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{N_{\mathit{\mu }}}\left[{\mathit{\beta }}_{\mathrm{2}i-\mathrm{1}},\cos ,\left(,i,\mathit{\omega },t,\right),+,{\mathit{\beta }}_{\mathrm{2}i-\mathrm{1}},\sin ,\left(,i,\mathit{\omega },t,\right)\right],& & \mathit{\sigma }\left(t\right)={\mathit{\alpha }}_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{N_{\mathit{\sigma }}}\left[{\mathit{\alpha }}_{\mathrm{2}i-\mathrm{1}},\cos ,\left(,i,\mathit{\omega },t,\right),+,{\mathit{\alpha }}_{\mathrm{2}i-\mathrm{1}},\sin ,\left(,i,\mathit{\omega },t,\right)\right],& & \mathit{\epsilon }\left(t\right)={\mathit{\gamma }}_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{N_{\mathit{\epsilon }}}\left[{\mathit{\gamma }}_{\mathrm{2}i-\mathrm{1}},\cos ,\left(,i,\mathit{\omega },t,\right),+,{\mathit{\gamma }}_{\mathrm{2}i-\mathrm{1}},\sin ,\left(,i,\mathit{\omega },t,\right)\right],\end{array}$$ where ${\mathit{\beta }}_{\mathrm{0}}$, ${\mathit{\alpha }}_{\mathrm{0}}$, and ${\mathit{\gamma }}_{\mathrm{0}}$ are the stationary components; ${\mathit{\beta }}_i$, ${\mathit{\alpha }}_i$, and ${\mathit{\gamma }}_i$ are the harmonics' amplitudes; $\mathit{\omega }$ $=$ 2$\mathit{\pi }{T}^{-\mathrm{1}}$ is the angular frequency, with $T$ corresponding to 1 year; $N_{\mathit{\mu }}$, $N_{\mathit{\sigma }}$, and $N_{\mathit{\epsilon }}$ are the number of harmonics; and $t$ is expressed in years. Therefore, the parameters ${\mathit{\beta }}_i$, ${\mathit{\alpha }}_i$, and ${\mathit{\gamma }}_i$ have been optimized through a non-stationary MLE in order to fit the monthly maxima of the non-stationary series. Different combinations of $N_{\mathit{\mu }}$, $N_{\mathit{\psi }}$, and $N_{\mathit{\epsilon }}$ have been tested and the best model was chosen as the one presenting the lowest value of the Akaike criterion (Akaike, 1973) given by $$\text{AIC}=\mathrm{2}k-\mathrm{2}\log \left(L\right),$$ where $k$ is the number of degree of freedoms of the model and $L$ is the likelihood. In particular, the maximum value tested for $N_{\mathit{\mu }}$ and $N_{\mathit{\psi }}$ is 3, while the maximum considered value of $N_{\mathit{\epsilon }}$ is 2. In general, this model can be extended to incorporate long-term trends, but the two series examined in this test display flat trends. Hence, Eqs. (62)–(64) are adequate to model them.

Seasonal cycle estimated by TS methodology and by the EM for the series of significant wave height of La Spezia and Ortona. The red continuous (dashed) line represents the location parameter $\mathit{\mu }$ estimated by TS (EM). The green continuous (dashed) line represents the sum between the location parameter $\mathit{\mu }$ and the shape parameter $\mathit{\sigma }$ estimated by TS (EM). The dots represent the monthly maxima. The shape parameters ${\mathit{\epsilon }}_{\text{TS}}$ and ${\mathit{\epsilon }}_{\text{EA}}$ estimated by the two methodologies have been also reported for the two series.

[Figure omitted. See PDF]

Return levels for La Spezia and Ortona for the month of January, estimated by the TS methodology (black continuous line) and by the EM (black dashed line). The green area represents the 95 % confidence interval estimated by the TS approach.

[Figure omitted. See PDF]

In the comparison, the EM and the seasonal TS methodology (GEV only) were applied to the same series of significant wave heights relative to the WWIII_MED data set described in Sect. 2.3. For the transformed-stationary approach, a 10-year time window was used for the computation of the long-term trend. The results of the two methodologies are similar, with a roughly flat trend and strong seasonal pattern. The comparison of the seasonal cycles estimated by the two techniques is represented in Fig. 10 for the two series. Here, the continuous red and green lines are the location and scale parameters ($\mathit{\mu }$ and $\mathit{\sigma }$, respectively) as estimated by the TS approach. The dashed red and green lines are the location and scale parameters estimated through the EM. The blue dots represent the monthly maxima, while the color scale represents the time-varying probability density estimated by the transformed-stationary methodology. Since for both of the series the models selected based on the Akaike criterion have a constant shape parameter $\mathit{\epsilon }$, these are reported together with those estimated by the TS methodology.

The GEV parameters estimated by the two approaches are in good agreement. The small differences have relatively small impact on the return levels as one can see in Fig. 11, where the return levels estimated by the two methodologies for the month of January are plotted. For both series, the return levels estimated by EM lie within the 95 % confidence interval estimated by TS. Table 4 reports the values of normalized bias (NBI) between the return levels estimated by TS and EM, defined as in Eq. (61), and the mean 95 % confidence interval amplitude expressed as a percentage of the return level. In Table 4 the values of NBI are reported for the four seasons for return periods of 5, 10, 30, 50, and 100 years, for both La Spezia and Ortona. In the definition of seasons that is used, winter starts on 1 December, spring on 1 March, summer on 1 June, and autumn on 1 September. We did not report return levels of periods greater than 100 years because the extension of the data covers only 35 years, hence the estimates for such periods are inaccurate for both methodologies. The average deviation between RL${}_{\text{TS}}$ and RL${}_{\text{cmp}}$ for the considered time series is rather small and remains below 7 % for all seasons. The confidence intervals estimated for TS are smaller than those estimated for EM, because the stationary MLE of TS has fewer degrees of freedom than the non-stationary one of EM, and is therefore affected by smaller uncertainty.

Normalized bias between the return levels estimated by the TS methodology and the EM methodology for the estimation of the seasonal variations, and mean 95 % confidence interval amplitude expressed as percentage of the return level, for return periods of 5, 10, 30, 50, and 100 years, for the four seasons, for significant wave height in La Spezia and Ortona.

Return period		5 years	10 years	30 years	50 years	100 years
La Spezia	NBI winter	1.19 %	1.51 %	1.95 %	2.14 %	2.39 %
(Waves Hs)	NBI spring	0.59 %	0.55 %	0.59 %	0.64 %	0.71 %
	NBI summer	4.75 %	5.28 %	5.99 %	6.27 %	6.62 %
	NBI autumn	$-$1.17 %	$-$1.03 %	$-$0.78 %	$-$0.66 %	$-$0.50 %
Mean conf. int. (TS)		2.68 %	3.05 %	3.63 %	3.90 %	4.25 %
Mean conf. int. (EM)		5.90 %	6.72 %	8.01 %	8.59 %	9.35 %
Ortona	NBI winter	3.74 %	4.23 %	4.91 %	5.20 %	5.57 %
(Waves Hs)	NBI spring	4.26 %	4.39 %	4.62 %	4.74 %	4.91 %
	NBI summer	$-$3.66 %	$-$3.44 %	$-$3.07 %	$-$2.90 %	$-$2.66 %
	NBI autumn	1.41 %	1.45 %	1.59 %	1.68 %	1.81 %
Mean conf. int. (TS)		3.18 %	3.75 %	4.70 %	5.15 %	5.78 %
Mean conf. int. (EM)		5.21 %	5.92 %	7.10 %	7.67 %	8.45 %

Discussion

Extreme value analysis is a subject of broad interest not only for earth science but also for other disciplines such as economy and finance (e.g., Gençay and Selçuk, 2004; Russo et al., 2015), sociology (e.g., Feuerverger and Hall, 1999), geology (e.g., Caers et al., 1996), and biology (e.g., Williams, 1995), among others. As a consequence, non-stationarity of signals is a common problem (e.g., Gilleland and Ribatet, 2014). In this respect, it is important to stress that the TS methodology is general, and its applicability only requires the stationarity of the transformed signal. Therefore, even if in this study the technique was applied only to series related to earth science, it can be employed in all disciplines dealing with extremes.

Given that the extreme value statistical model is an important component of applications such as those discussed here (e.g., Coles, 2001; Hamdi et al., 2014), it is important to stress that the theory was formulated in a way that is not restricted to GEV and GPD, but can be extended to any statistical model for extreme values. In particular, since the GEV distribution is a generalization of the Gumbel, Frechet, and Weibull statistics, TS can be reformulated separately for these three distributions, as well as for the commonly used $r$-largest approach statistics (e.g., Coles, 2001; Hamdi et al., 2014). Finally, an extension of TS to statistical models not based on the GEV theory (e.g., Boccotti, 2000; Goda, 1988) may open the way to their non-stationary generalization and could be an interesting direction for future research.

The transformation consists in simple, time-varying normalization of the signal through the estimation of trend, slowly varying standard deviation and seasonality, and allows different types of analysis. The first product of the methodology is its capability to estimate the extreme values of the signal. Next, the TS approach enables the analysis of long-term variability. As an example, it was shown to be useful in relating the long-term trend of the signal with the NAO climatic index (see Sect. 3.3). Finding correlations of natural parameters with climatic indices is a theme of common interest in earth science, especially in view of climate change (e.g., Barnard et al., 2015; Dodet et al., 2010; Plomaritis et al., 2015). If a time series is correlated to a climatic index in the long term, an advantage of the TS methodology is that it can model extremes correlated to the index without considering it explicitly in the computation. Finally, the TS methodology allows to describing the seasonal variability of extremes, which is also critical for climate studies (e.g., Sartini et al., 2015; Menendez et al., 2009; Méndez et al., 2006).

As shown in Sect. 4, the TS methodology has advantages over SS (e.g., Vousdoukas et al., 2016a) and EM (e.g., Cheng et al., 2014; Gilleland and Katz, 2016; Izaguirre et al., 2011; Méndez et al., 2006; Menéndez et al., 2009; Mudersbach and Jensen, 2010; Russo et al., 2014; Sartini et al., 2015), both in terms of accuracy of the results and its conceptual and implementation simplicity. In particular, in the comparison with the SS methodology for long-term variability, the return levels estimated by the two techniques are similar for return periods for which the SS is accurate. The use of the whole time horizon of the series represents a major advantage of TS over SS because it allows more accurate estimations of the return levels associated with long return periods. A conceptual advantage of the TS methodology over EM is that it decouples the detection of the non-stationary behavior of the series from the fitting of the extreme value distribution. The study of the time-varying statistical features of the series is delegated to the transformation, and takes place before the fitting of the extreme value distribution. This fact provides a simple diagnostic tool to evaluate the validity of the model applied to a particular series: the model is valid if the transformed series is stationary. This is useful for validating the output of the approach. Moreover, the decoupling simplifies both the detection of non-stationary patterns and the fitting of the extreme value distribution. In particular, the detection of non-stationary patterns can be accomplished by means of simple statistical techniques such as low-pass filters based on the running mean and standard deviation, and the fitting of the extreme value distribution can be obtained through a stationary MLE with a small number of degrees of freedom that is easier to implement and control. Moreover, unlike many implementations of EM (e.g., Cheng et al., 2014; Gilleland and Katz, 2016; Izaguirre et al., 2011; Méndez et al., 2006; Menéndez et al., 2009; Sartini et al., 2015; Serafin and Ruggiero, 2014), the detection of non-stationary patterns described in this paper does not require an input parametric function $M$ for the variability. This makes the TS methodology well suited for massive applications with the simultaneous evaluation of many time series, for which a common definition of $M$ would be difficult (e.g., Vousdoukas et al., 2016a).

It is worth remarking that the EM implemented, for example, using Eq. (62), is able to model a shape parameter varying in time, unlike the TS using the transformation given by Eq. (1). While in principle this is a weak point of the TS methodology described here, assuming a constant shape parameter is a reasonable assumption for most cases, because in general simple models should be preferred to complex ones (e.g., Coles, 2001). In particular, using EM the Akaike criterion (Akaike, 1973), that favors simple models with fewer degrees of freedoms, often selects models with a fixed shape parameter (e.g., Sartini et al., 2015; Menendez et al., 2009). Moreover, the finding that a non-stationary GEV always corresponds to a transformation of the non-stationary time series into a stationary one, shown in Appendix A, suggests that a generalization of the TS methodology is possible in order to include models with time-varying shape parameters.

Conclusions

This paper describes the TS methodology for non-stationary extreme value analysis. The main assumption underlying this approach is that if a non-stationary time series can be transformed into a stationary one to which the stationary EVA theory can be applied, then the result can be back-transformed into a non-stationary extreme value distribution through the inverse transformation. The proposed methodology is general and, even if in this study we applied it only to series related to earth science, it can be employed in all disciplines dealing with EVA. Moreover, though we discussed it only for GEV and GPD, it can be extended to any other statistical model for extremes.

As a transformation we proposed a simple time-varying normalization of the signal estimated by means of a time-varying mean and standard deviation. This simple transformation was also adapted to describe the seasonal variability of the extremes. In addition, it was proven to provide a comprehensive model for non-stationary GEV and GPD distributions with a constant shape parameter, which means that it can be applied to a wide range of non-stationary processes. The formal duality between the TS and more established approaches has also been proven, suggesting that a complete generalization of the TS approach would allow including models with a time-varying shape parameter.

The methodology was tested on time series of different variables, sizes, and statistical properties. An evaluation of the statistical error associated with the transformation showed that, for the examined series, this is negligible with respect to the error associated with the stationary MLE (the squared error is 2 orders of magnitude smaller) and to that related to the estimation of the threshold for GPD.

The TS methodology was compared with a stationary EVA applied on quasi-stationary slices of non-stationary series (i.e., SS) for the estimation of the long-term variability of extremes, and with the EM to non-stationary EVA. The return levels estimated by TS are shown to be comparable to those obtained by these two methodologies. However, the TS approach has advantages over both SS and EM. With respect to SS, the TS uses the whole time series for fitting the extreme value distribution, guaranteeing a more accurate estimation at larger return periods. With respect to EM, the TS decouples the detection of the non-stationarity of the series from the fit of the extreme value distribution, involving a simplification of both steps of the analysis. In particular, the fit of the distribution can be accomplished using a simple MLE with a few degrees of freedom and is easy to implement and control. The detection of non-stationarity can be performed by means of easily implemented and fast low-pass filters, which do not require as input any parametric function for the variability. This makes the methodology well suited for massive applications where the simultaneous evaluation of several time series is required.

An implementation of the TS methodology has been developed in an open-source MATLAB toolbox (tsEva), which is available at https://github.com/menta78/tsEva/.

Data availability

The source code used to implement the case studies presented in this work is available at https://github.com/menta78/tsEva/ (Mentaschi et al., 2016).

Duality between the established method and the TS methodology

Here, we show that if the extremes of a time series $y\left(t\right)$ are fitted by a non-stationary distribution GEV${}_Y(y$, $t)$ then there is a family of transformations $f(y$, $t)$ : $y\left(t\right)$ $\to$ $x\left(t\right)$ such that GEV${}_Y(y$, $t)$ $=$ $GE{V}_X\left[f^{-\mathrm{1}}\right(x$, $t\left)\right]$, where GEV${}_X\left(x\right)$ is a stationary GEV fitting the extremes of a supposed stationary series $x\left(t\right)$.

To prove this, we expand relationship GEV${}_Y(y$, $t)$ $=$ GEV${}_X\left[f^{-\mathrm{1}}\right(x$, $t\left)\right]$, finding $${\left\{\mathrm{1},+,{\mathit{\epsilon }}_x,\left[{\displaystyle \frac{f(y,t)-{\mathit{\mu }}_x}{{\mathit{\sigma }}_x}}\right]\right\}}^{\mathrm{1}/{\mathit{\epsilon }}_x}={\left\{\mathrm{1},+,{\mathit{\epsilon }}_y,\left(,t,\right),\left[{\displaystyle \frac{y-{\mathit{\mu }}_y\left(t\right)}{{\mathit{\sigma }}_y\left(t\right)}}\right]\right\}}^{\mathrm{1}/{\mathit{\epsilon }}_y\left(t\right)},$$ where [${\mathit{\epsilon }}_y\left(t\right)$, ${\mathit{\sigma }}_y\left(t\right)$, ${\mathit{\mu }}_y\left(t\right)]$ are the time-varying parameters of GEV${}_Y(y$, $t)$ and [${\mathit{\epsilon }}_x$, ${\mathit{\sigma }}_x$, ${\mathit{\mu }}_x$] are the constant parameters of GEV${}_X\left(x\right)$. Solving for $f(y$, $t)$, we find $$f(y,t)={\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}_x}}\left\{{\mathit{\sigma }}_x,{\left[\mathrm{1},+,{\mathit{\epsilon }}_y,\left(,t,\right),\left({\displaystyle \frac{y-{\mathit{\mu }}_{y}t}{{\mathit{\sigma }}_y\left(t\right)}}\right)\right]}^{\frac{{\mathit{\epsilon }}_x}{{\mathit{\epsilon }}_y\left(t\right)}},-,{\mathit{\sigma }}_x,+,{\mathit{\epsilon }}_x,{\mathit{\mu }}_x\right\}.$$ Equation (A2) defines a family of functions because the values of the stationary GEV parameters [${\mathit{\epsilon }}_x$, ${\mathit{\sigma }}_x$, ${\mathit{\mu }}_x$] can be assigned arbitrarily. Furthermore, if we chose ${\mathit{\epsilon }}_x$ $\ne$ 0 then $f(y$, $t)$ is monotonic in $y$ for every time $t$ and can therefore be inverted, while for ${\mathit{\epsilon }}_x$ $=$ 0 a Gumbel-specialized formulation can be derived from Eq. (A1).

In the particular case of ${\mathit{\epsilon }}_y$ $=$ const. $=$ ${\mathit{\epsilon }}_x$ function $f(y$, $t)$ reduces to $$f(y,t)={\displaystyle \frac{y-{\mathit{\mu }}_y\left(t\right)+{\mathit{\mu }}_x/{\mathit{\sigma }}_x\cdot {\mathit{\sigma }}_y\left(t\right)}{{\mathit{\sigma }}_y\left(t\right)/{\mathit{\sigma }}_x}},$$ which is equivalent to Eq. (1) provided that $T_y$ $=$ ${\mathit{\mu }}_y$ $-$ ${\mathit{\mu }}_x/$ ${\mathit{\sigma }}_x$ $\cdot$ ${\mathit{\sigma }}_y$ and $C_y$ $=$ ${\mathit{\sigma }}_y/{\mathit{\sigma }}_x$. Hence, we can say that Eq. (1) allows a general TS formulation for models with a constant shape parameter, because we can arbitrarily impose ${\mathit{\epsilon }}_x$ $=$ ${\mathit{\epsilon }}_y$ in Eq. (A2) if we assume a constant ${\mathit{\epsilon }}_y$. This finding is remarkable because it proves that any non-stationary GEV model with constant ${\mathit{\epsilon }}_y$ can be connected to Eq. (1).

Equation (A2) alone is not enough to formulate a fully generalized TS approach, because in Eq. (A2) the non-stationary GEV parameters [${\mathit{\epsilon }}_y\left(t\right)$, ${\mathit{\sigma }}_y\left(t\right)$, ${\mathit{\mu }}_y\left(t\right)$] are regarded as known variables, which is an incorrect assumption in practical applications. But it is enough to say that any implementation of the non-stationary established method is equivalent to a transformation into a supposed stationary series $x\left(t\right)$. Therefore, Eq. (A2) could be used as a diagnostic tool for implementations of the established method: a condition for the validity of the non-stationary model is that the transformed $x\left(t\right)$ series is stationary.

Acknowledgements

The authors would like to thank Simone Russo of the JRC and Francesco Fedele of the GIT for the precious suggestions, and Niall McCormick of the JRC for the careful review of this manuscript. This work was co-funded by the JRC exploratory research project Coastalrisk and by the European Union Seventh Framework Programme FP7/2007-2013 under grant agreement no. 603864 (HELIX: High-End cLimate Impacts and eXtremes; http://www.helixclimate.eu). Edited by: S. Archfield Reviewed by: two anonymous referees