2.1 The Mann–Kendall test

Hydrological time series are often composed by non-normally independent realizations of phenomena, and this characteristic makes the use of nonparametric trend tests very attractive (Kundzewicz and Robson, 2004). The Mann–Kendall test is a widely used rank-based tool for detecting monotonic, and not necessarily linear, trends. Given a random variable z, and assigned a sample of L independent data z=(z₁, …, z_L), the Kendall S statistic (Kendall, 1975) can be defined as follows:

where “sgn” is the sign function.

The null hypothesis of this test is the absence of any statistically significant trend in the sample, whereas the presence of a trend represents an alternative hypothesis. Yilmaz and Perera (2014) reported that serial dependence can lead to a more frequent rejection of the null hypothesis. For L≥8, Mann (1945) reported that Eq. (1) is approximatively a normal variable with a zero mean and variance which, in the presence of t_m ties of length m, can be expressed as

In practice, the Mann–Kendall test is performed using the Z statistic

which follows a standard normal distribution. Using this approach, it is simple to evaluate the p value and compare it with an assigned level of significance or, equivalently, to calculate the Z_α threshold value and compare it with Z, where Z_α is the (1−α) quantile of a standard normal distribution.

Yue et al. (2002b) observed that autocorrelation in time series can influence the ability of the MK test to detect trends. To avoid this problem, a correct approach with respect to trend analysis should contemplate a preliminary check for autocorrelation and, if necessary, the application of pre-whitening procedures.

A nonparametric tool for a reliable estimation of a trend in a time series with N pairs of data is the Sen's slope estimator (Sen, 1968), which is defined as the median of the set of slopes δ_j:

where i>k.

2.2 Likelihood ratio test

The likelihood ratio statistical test allows for the comparison of two candidate models. As its name suggests, it is based on the evaluation of the likelihood function of different models.

The LR test has been used multiple times (Tramblay et al., 2013; Cheng et al., 2014; Yilmaz et al., 2014) to select between stationary and nonstationary models in the context of nested models. Given a stationary model characterized by a parameter set θ_st and a nonstationary model, with parameter set θ_ns, if $l\left({\widehat{\mathit{\theta }}}_{\mathrm{st}}\right)$ and $l\left({\widehat{\mathit{\theta }}}_{\mathrm{ns}}\right)$ are their respective maximized log likelihoods, the likelihood ratio test can be defined using the deviance statistic

D is (for large L) approximately ${\mathit{\chi }}_m^{\mathrm{2}}$ distributed, with $m=\dim \left({\mathit{\theta }}_{\mathrm{ns}}\right)-\dim \left({\mathit{\theta }}_{\mathrm{st}}\right)$ degrees of freedom. The null hypothesis of stationarity is rejected if D>C_α, where C_α is the (1−α) quantile of the ${\mathit{\chi }}_m^{\mathrm{2}}$ distribution (Coles, 2001).

Besides the analysis of power, we also checked (in Sect. 3.3) the approximation $D\sim {\mathit{\chi }}_m^{\mathrm{2}}$ as a function of the sample size L for the evaluation of the level of significance.

2.3 Akaike information criterion ratio test

Information criteria are useful tools for model selection. It is reasonable to retain that the Akaike information criterion (AIC; Akaike, 1974) is the most famous among these tools. Based on the Kullback–Leibler discrepancy measure, if θ is the parameter set of a k-dimensional model (k=dim(θ)), AIC is defined as

The model that best fits the data has the lowest value of the AIC between candidates. It is useful to observe that the term proportional to the number of model parameters allows one to account for the increase of the estimator variance as the number of model parameters increases.

Figure 1An empirical distribution of AIC_R and the rejection threshold AIC_R,α of the null hypothesis (stationary GEV parent).

Download

Sugiura (1978) observed that the AIC can lead to misleading results for small samples; thus, he proposed a new measure for the AIC:

where a second-order bias correction is introduced. Burnham and Anderson (2004) suggested only using this version when $L/k_{\max }<\mathrm{40}$, with k_max being the maximum number of parameters between the models compared. However, for larger L, AIC_c converges to AIC. For a quantitative comparison between the AIC and AIC_c in the extreme value stationary model selection framework, the reader is referred to Laio et al. (2009).

In order to select between stationary and nonstationary candidate models, we use the ratio

where the subscripts indicate the AIC value obtained for a stationary (st) and a nonstationary (ns) model, both fitted with maximum likelihood to the same data series.

Considering that the better fitting model has a lower AIC, if the time series arises from a nonstationary process, the AIC_R should be less than 1; the opposite is true if the process is stationary.

In order to provide a rigorous comparison between the use of the MK, LR, and AIC_R, we evaluated the AIC_R,α threshold value corresponding to the significance level α using numerical experiments.

More in detail, we adopted the following procedure:

1. N=10 000 samples are generated from a stationary GEV parent distribution, with known parameters;

2. for each of these samples the AIC_R is evaluated by fitting the stationary and nonstationary GEV models described in Sect. 2.4, thus providing its empirical distribution (see probability density function, pdf, in Fig. 1);

3. exploiting the empirical distribution of AIC_R, the threshold associated with a significance level of α=0.05 is numerically evaluated. This value, AIC_R,α, represents the threshold for rejecting the null hypothesis of stationarity (which in these generations is true) in 5 % of the synthetic samples.

This procedure was applied both for the AIC and AIC_c. The experiment was repeated for a few selected sets of the GEV parameters, including different trend values, and different sample lengths, as detailed in Sect. 3.

2.4 The GEV parent distribution

The cumulative distribution function of the generalized extreme value (GEV) distribution (Jenkinson, 1955) can be expressed as follows:

where ζ, σ, and ε are known as the position, scale, and shape parameters, respectively; θ_st=[ζ, σ, ε] is a general and comprehensive way to express the parameter set in the stationary case. The flexibility of the GEV, which accounts for the Gumbel, Fréchet, and Weibull distributions as special cases (for ε=0, ε>0 and ε<0 respectively) makes it eligible for a more general discussion about the implications of nonstationarity.

Traditional extreme value distributions can be used in a nonstationary framework, modeling their parameters as function of time or other covariates (Coles, 2001), producing ${\mathit{\theta }}_{\mathrm{st}}⟶{\mathit{\theta }}_{\mathrm{ns}}=[{\mathit{\zeta }}_t$, σ_t, ε_t].

In this study, only a deterministic linear dependence on the time t of the position parameter ζ has been introduced, leading Eq. (7) to be expressed as follows:

with

and θ_ns=[ζ₀, ζ₁, σ, ε].

It is important to note that Eq. (8) is a more general way of defining the GEV and has the property of degenerating into Eq. (7) for ζ₁=0; in other words, Eq. (7) represents a nested model of Eq. (8) which would confirm the suitability of the likelihood ratio test for model selection.

According to Muraleedharan et al. (2010), the first three moments of the GEV distribution are as follows:

Here, $g_k=\mathrm{\Gamma }(\mathrm{1}-k\mathit{\epsilon })$, with $k\in Z^{+}$ and Γ(⋅), is the gamma function. It is worth noting that, following Eqs. (10)–(12), the trend in the position parameter only affects the mean, while the variance and skewness remain constant.

In this work, we used the maximum likelihood method (ML) to estimate the GEV parameters from sample data. The ML allows one to treat ζ₁ as an independent parameter, as well as ζ₀, σ and ε. For this purpose, we exploited the “extRemes” R package (Gilleland and Katz, 2016).

2.5 Numerical evaluation of test power and significance level

The power of a test is related to the type II error and is the probability of correctly rejecting the null hypothesis when it is false. In particular, defining α (level of significance), the probability of a type I error, and β, the probability of a type II error, we have a power value of 1−β. The maximum value of power is 1, which correspond to β=0, i.e., no probability of a type II error. In most applications, the conventional values are α=0.05 and β=0.2, meaning that a 1-to-4 trade-off between α and β is accepted. Thus, in our experiment we always assumed a significance level of 0.05, and, for the following description of results and discussion, we considered a power level of less than 0.8 to be too low and, hence, unacceptable. In Sect. 4, we report further considerations regarding this choice. For each of the tests described in Sect. 2.1, 2.2, and 2.3, the power was numerically evaluated according to the following procedure:

1. N=2000 Monte Carlo synthetic series, each of length L, are generated using the nonstationary GEV in Eqs. (8) and (9) as a parent distribution with a fixed parameter set θ_ns=[ζ₀, ζ₁, σ, ε] with ζ₁≠0.

2. The threshold AIC_R,α associated with a significance level of α=0.05 is numerically evaluated, as described in Sect. 2.3, using the corresponding parameter set θ_st=[ζ₀, σ, ε] of the GEV parent distribution.

3. From these synthetic series, the power of the test is estimated as

   $$\mathrm{rejection}\mathrm{rate}={\displaystyle \frac{N_{\mathrm{rej}}}{N}},$$

   where N_rej is the number of series for which the null hypothesis is rejected, as in Yue et al. (2002a).

The same procedure, with N=10 000, was used in order to check the actual significance level of the test, which is the probability of type I error, i.e., the probability of rejecting the null hypothesis of stationarity when it is true. The task was performed by following the abovementioned steps 1 to 3 while replacing θ_ns with θ_st in step 1); in such a case, the rejection rate N_rej∕N represents the actual level of significance α.

We used a reduced number of generations (N=2000) for the evaluation of power as a good compromise between the quality of the results and computational time. N=2000 was also used by Yue et al. (2002a).

3.1 Evaluation of AIC_R, with the AIC and AIC_c

Considering the nonstationary GEV four-parameter model, in order to satisfy the relation $L/k_{\max }<\mathrm{40}$ suggested by Burnham and Anderson (2004), a time series with a record length no less than 160 should be available. Following this simple reasoning, the AIC should be considered not to be applicable to any annual maximum series showing a changing point between the 1970s and 1980s (e.g., Kiely, 1999). In our numerical experiment, the second-order bias correction of Sugiura (1978) should always be used, as we have $L/k_{\max }=\mathrm{70}/\mathrm{4}=\mathrm{17.5}$ for the nonstationary GEV for a maximum sample length of L=70. Nevertheless, we checked if using the AIC or AIC_c may affect results. For this purpose, we evaluated the percentage differences between the power of the AIC_R obtained by means of the AIC and AIC_c from synthetic series. In Fig. 2, the empirical probability density functions of such percentage differences, grouped according to sample length, are plotted for generations with ε=0.4 and different values of σ. It is interesting to note that the error distribution only shows a regular and unbiased bell-shaped distribution for L=70. We then observe a small negative bias (about −0.02 %) for L=50, while a bias of −0.08 with a multi-peak and negatively skewed pdf is noted for L=30. The latter pdf also has a higher variance than the others. The purpose of this figure is to show that the difference between the power obtained with the AIC and the power obtained with the AIC_c is negligible. Different peaks in one curve (L=30) can be explained by the merging of sample errors obtained for different values of σ. Similar results were obtained for all values of ε, which always provided very low differences and allow for the conclusion to be reached that the use of the AIC or AIC_c does not significantly affect the power of AIC_R for the cases examined. This follows the combined effect of the sample size (whose minimum value considered here is 30) and the limited difference in the number of parameters in the selected models. In the following, we will refer to and show only the plots obtained for the AIC_R in Eq. (6) with the AIC evaluated as in Eq. (4).

3.2 Dependence of the power on the parent distribution parameters and sample size

The effect of the parent distribution parameters and the sample size on the numerical evaluation of the power and significance level of the MK, LR, and AIC_R tests for different values of ε, σ, and ζ₁ is shown in Fig. 3. The curves represent both the significance level, which is shown for ζ₁=0 (true parent is the stationary GEV), and the power, which is shown for all other values ζ₁≠0 (true parent is the nonstationary GEV). Each panel in Fig. 3 shows the dependence of the power and significance level of MK, LR and AIC_R on the trend coefficient for one set of parameter values and different sample sizes. In all panels, the test power strongly depends on the trend coefficient and sample size. This dependence is also affected by the parent parameter values. In all cases, the power reaches 1 for a strong trend and approaches 0.05 (the chosen level of significance) for a weak trend (ζ₁ close to 0). In all combinations of the shape and scale parameters (and especially for short samples) for a wide range of trend values, the power exhibits values well below the conventional value of 0.8. The curves' slope between 0.05 and 1 is sharp for long samples and gentle for short samples. It also depends on the parameter set, with slopes generally being gentler for higher values of the scale (σ) and shape (ε) parameters of the parent distribution. A significant difference in the power between the MK, LR, and AIC_R tests is observable when the sample size is smaller and even more so when the parent distribution is heavy-tailed ($\mathit{\epsilon }=+\mathrm{0.4}$).

Figure 3Dependence of test power on the trend coefficient, sample size, scale, and shape of the parent parameters.

Download

In particular, for ε=0, −0.4 and L=50, 70, it is possible to report a slightly larger power of LR with respect to the AIC_R and MK, but values are very close to each other. However, the reciprocal position of MK and AIC_R power curves is interesting; in fact, the AIC_R power is always larger than that of the MK, except when $\mathit{\epsilon }=-\mathrm{0.4}$, for all values of the scale parameter.

A higher difference is found for a heavy-tailed parent distribution ($\mathit{\epsilon }=+\mathrm{0.4}$). While LR still has the largest power value, the difference with respect to AIC_R remains small and the MK power value almost always collapses to values smaller than 0.5.

The practical consequences of such patterns are very important and are discussed in Sect. 4.

3.3 Sensitivity and evaluation of the actual significance level

We evaluated the threshold values (corresponding to a significance level of 0.05) for accepting/rejecting the null hypothesis of stationarity according to the methodologies recalled in Sect. 2.1 and 2.2 for the MK and LR tests and introduced in Sect. 2.3 for AIC_R. Based on these thresholds, we exploited the generation of series from a stationary model (ζ₁=0) in order to numerically evaluate the rate of rejection of the null hypothesis, i.e., the actual significance level of the tests considered in the numerical experiment, following the procedure described in Sect. 2.5.

Table 1 shows the numerical values of the actual level of significance, obtained numerically, to be compared with the theoretical value of 0.05 for all of the sets of parameters and sample sizes considered. Among the three measures for trend detection, the LR shows the worst performance. The results in Table 1 show that the rejection rate of the (true) null hypothesis is systematically higher than it should be, and it is also dependent on parent parameter values. This effect is exalted when the parent distribution has an upper boundary ($\mathit{\epsilon }=-\mathrm{0.4}$) and for shorter series (L=30). In practice, this implies that when using the LR test, as described in Sect. 2.2, there is a higher probability of rejecting the null hypothesis of stationarity (if it is true) than expected or designed.

Table 1The actual level of significance of the tests for different sample sizes, scales, and shapes of the parent parameters.

Download Print Version | Download XLSX

Conversely, the performance of MK with respect to the designed level of significance is less biased and is independent of the parameter set. Similar good performance is trivially obtained for the AIC_R, whose rejection threshold is numerically evaluated.

Figure 4Enlargement of the power test curves in the case (σ=15, $\mathit{\epsilon }=-\mathrm{0.4}$), with focus on the actual level of significance (ζ₁=0).

Download

The plot in Fig. 4 is displayed in order to focus on the actual value of the level of significance and, in particular, on the LR approximation $D\sim {\mathit{\chi }}_m^{\mathrm{2}}$ as a function of the sample length L. The difference between the theoretical and numerical values of the significance level is represented by the distance between the bottom value of the curve (obtained for ζ₁=0, i.e., the stationary GEV model) and the chosen level of significance 0.05, represented by the horizontal dotted line. In particular, in Fig. 4, results for the parameter set (σ=15, $\mathit{\epsilon }=-\mathrm{0.4}$) show that the actual rate of rejection is always higher than the theoretical one and changes significantly with the sample size; this means that the ${\mathit{\chi }}_m^{\mathrm{2}}$ approximation leads to a significant underestimation of the rejection threshold of the D statistic. Moreover, it seems that the LR power curves (in red) are shifted toward higher values as a consequence of the significance level overestimation, meaning that the LR test power is also overestimated due to the approximation $D\sim {\mathit{\chi }}_m^{\mathrm{2}}$. These results suggest the use of a numerical procedure for the LR test (such as that introduced for AIC_R in Sect. 2.3) for evaluating the D distribution and the rejection threshold.

Table 2Actual level of significance of the AIC_R test for AIC${}_{R,\mathit{\alpha }}=\mathrm{1}$.

Download Print Version | Download XLSX

Figure 5AIC_R,α thresholds for different parameter sets vs. sample size.

Download

Other considerations can be made regarding the use of AIC_R. As explained in Sect. 2.3, we empirically evaluated the AIC_R,α threshold value using numerical generations with a significance level 0.05 for each of the parameter sets and sample sizes considered. Similar results were obtained using the AIC_c, which are not shown for brevity. We found a significant dependence of AIC_R,α on the sample size. Figure 5 shows the AIC_R,α curves obtained for each of the parameter sets vs. sample size. It is also worth noting that all curves asymptotically trend to 1 as L increases. This property is due to the structure of the AIC and the peculiarity of the nested models used in this paper: while using a sample generated with weak nonstationarity (i.e., when ζ₁→0 in Eq. 9), the maximum likelihood of the model shown in Eq. 7, $l\left({\widehat{\mathit{\theta }}}_{\mathrm{st}}\right)$, tends toward $l\left({\widehat{\mathit{\theta }}}_{\mathrm{ns}}\right)$ of the model shown in Eq. 8, leaving only the bias correction (2k in equation 4) to discriminate between competing AIC values in model selection applications. As a consequence, the AIC_R,α should always be lower than 1; however, when increasing the sample size, both the likelihood terms $-\mathrm{2}l\left({\widehat{\mathit{\theta }}}_{\mathrm{st}}\right)$ and $-\mathrm{2}l\left({\widehat{\mathit{\theta }}}_{\mathrm{ns}}\right)$ in Eq. (4) will also increase, pushing AIC_R toward the limit 1. Conversely, Fig. 5 shows that the threshold value AIC_R,α is significantly smaller than 1 up to L values well beyond the length usually available in this kind of analysis. Hence, the numerical evaluation of the threshold has to be considered as a required task in order to provide an assigned significance level to model selection. In contrast, the simple adoption of the selection criteria AIC_R<1 (i.e., AIC${}_{R,\mathit{\alpha }}=\mathrm{1}$) would correspond to an unknown significance level that is dependent on the parent distribution and sample size. In order to highlight this point, we evaluated the significance level α corresponding to AIC${}_{R,\mathit{\alpha }}=\mathrm{1}$, following the procedure described in Sect. 2.5, by generating N=10 000 synthetic series (from a stationary model) for any parameter set and sample length. The results, provided in Table 2, show that α ranges between 0.16 and 0.26 in the explored GEV parameter domain and mainly depends on the sample length and the shape parameter of the parent distribution.

Figure 6Sample variability of ML-ζ₁ and δ vs. the trend coefficient ζ₁.

Download

3.4 Sample variability of parent distribution parameters

In our opinion, the results shown above, with respect to the performance of parametric and nonparametric tests, are quite surprising and important. It is proved that the preference widely accorded to nonparametric tests, due to the fact that their statistics are allegedly independent from the parent distribution, is not well founded. Conversely, the use of parametric procedures raises the problem of correctly estimating the parent distribution and, for the purpose of this paper, its parameters. Moreover, as the trend coefficient ζ₁ is a parameter of the parent distribution under nonstationary conditions, the proposed parametric approach provides a maximum likelihood-based estimation of the same trend coefficient, which is hereafter referred to as ML-ζ₁. For a comparison with nonparametric approaches, we also evaluated the sample variability of the Sen's slope measure (δ) of the imposed linear trend. Furthermore, in order to provide insights into these issues, we analyzed the sample variability of the maximum likelihood estimates ML-ε and ML-σ (from the same sets of generations exploited above) for different parameter sets and sample lengths.

We evaluated sample variability s[⋅], as the standard deviation of the ML estimates of parameter values obtained from synthetic series. In the upper panels of Fig. 6, we show s[ML-ζ₁], and in the lower panels, we show the Sen's slope median s[δ]. In both cases, the sample variability of the linear trend is strongly dependent on sample size and independent from the true ζ₁ value in the range examined [−1, 1]. It reaches high values for short samples and, in such cases, its dependence on the scale and shape parent parameters is also relevant. The ML estimation of the trend coefficient is always more efficient than Sen's slope, and this is observed for heavy-tailed distributions in particular.

In Fig. 7, we show the empirical distributions of the Sen's slope δ and ML-ζ₁ estimates obtained from samples of size L=30 with a parent distribution characterized by σ=15 and $\mathit{\epsilon }=[-\mathrm{0.4}$, 0, 0.4], providing visual information about the range of trend values that may result from a local evaluation. Similar results, characterized by smaller sample variability, as shown in Fig. 6, are obtained for L=50 and L=70 and are not shown for brevity.

Figure 7Empirical distributions of δ and ML-ζ₁ evaluated from samples with L=30 and σ=15 vs. the trend coefficient ζ₁.

Download

Figure 8 shows the sample variability of ML-ε and ML-σ, which is still independent of the true ζ₁ for values of ε=0 and 0.4, whereas for the upper-bounded GEV distributions ($\mathit{\epsilon }=-\mathrm{0.4}$) it shows a significant increase for higher values of σ and high trend coefficients ($\left|{\mathit{\zeta }}_{\mathrm{1}}\right|>\mathrm{0.5}$). The randomness of results for L=30 and σ=[15, 20] is probably due to the reduced efficiency of the algorithm that maximizes the log-likelihood function for heavy-tailed distributions.

Figure 8Sample variability of ML-ε and ML-σ vs. the trend coefficient ζ₁.

Download

Figure 9Empirical distributions of ML-ε evaluated for σ=15 from samples with L=30 and L=70 vs. the trend coefficient ζ₁.

Download

Figure 10Empirical distributions of ML-σ evaluated for σ=15 from samples with L=30 and L=70 vs. the trend coefficient ζ₁.

Download

In order to better analyze such patterns, for the scale and shape parent parameters we also report the distribution of their empirical ML estimates for different parameter sets vs. the true ζ₁ value used in generation. The sample distribution of ML-ε for σ=15 is shown in Fig. 9 for L=30 and L=70. The sample distribution of ML-σ for σ=15 is shown in Fig. 10 for L=30 and L=70. The panels show that the presence of a strong trend coefficient may produce significant loss in the estimator efficiency, which is probably due to deviation from the normal distribution of the sample estimates for long samples. This suggests the need for more robust estimation procedures that provide higher efficiency for estimates of ε and σ in the case of a strong observed trend. It should be highlighted that efficiency in the parameter estimation increases with sample size for ε=[0, 0.4], whereas it decreases for both ε and σ in the case of $\mathit{\epsilon }=-\mathrm{0.4},$ where the trend of the location parameter implies a shift in time of the distribution upper bound.