Evaluating the uncertainty and reliability of standardized indices

The spatio-temporal variability of hydrological variables is often described by defining and quantifying specific standardized indices (Mishra & Singh 2010). Meteorological standardized indices, such as the Standardized Precipitation Index (SPI) (Mckee et al. 1993) and the Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al. 2010), represent the most common applications, but there are also examples of indices relating to more complex hydrological variables such as the streamflow (Standardized Discharge Index (SDI), Nalbantis & Tsakiris 2009) or the runoff (Standardized Runoff Index, Shukla & Wood 2008). Recently, multivariate standardized indices have also been proposed (Hao & AghaKouchak 2013). Standardized indices express any observed value of the reference hydrological variable (i.e., precipitation, climatic water balance, streamflow, etc.) in terms of standard normal deviate, thus allowing a rapid and effective identification and quantification of the anomalies with respect to the ‘normal’ condition. The main advantage of standardization lies in the possibility of comparing the anomalies both in space and time, independently of any climatic or seasonal difference.

In the traditional parametric approach, the fundamental initial step in the quantification of a standardized index is the choice of a probability distribution that is suitable for describing the specific reference variable. Sometimes this step does not have a clear and univocal solution (Guttman 1999), thus introducing a first problematic issue: indeed, if different distributions are used, a single value of the reference variable corresponds to different standardized values (with a consequent reduction of the comparability feature). Several studies have investigated this aspect, mainly in relation to the SPI (Guttman 1999; Kumar et al. 2009; Blain 2011; Angelidis et al. 2012) and SPEI (Stagge et al. 2015).

To avoid these problems and minimize the uncertainty associated with the selection and estimation of parametric distribution functions, it is possible to derive the probability distribution of a certain hydrological variable by a nonparametric approach. This is typically obtained by using an empirical distribution function (Farahmand & AghaKouchak 2015) or a nonparametric kernel density estimator (Kumar et al. 2016).

Whatever the approach is, the values of the standardized indices are affected by uncertainty because the corresponding cumulative probabilities are estimated on the basis of a limited sample of observations. The quantification of such type of uncertainty can be very useful for practical applications (water resources management, prevention, etc.), because it makes it possible to define the reliability (or reasonableness) of the index estimate and of the dependent actions. However, this topic has rarely been addressed in the literature. No studies have discussed the uncertainty associated with nonparametric standardized indices. For the parametric approach, no studies are available for SPEI or RDI, and the only example for SDI can be found in Hong et al. (2015). More studies can be found for the SPI. For example, Naumann et al. (2012) and Hu et al. (2015) have applied bootstrap-based procedures to describe and quantify the confidence intervals (CI) associated with SPI estimates. Vergni et al. (2015) performed a more comprehensive bootstrap-based analysis to evaluate, for central Italy, the variability of the mean size of the CI of parametric SPI varying the time scale, the length of the time series, and the underlying distribution. In that analysis, the size (ΔS) of the CI was simply expressed in standard deviate units. Of course, this method is useful for comparing the contribution of different factors to the overall uncertainty, but it is not very effective for understanding the reliability of the estimate.

For this reason, a first objective of the present paper is to improve the methodology illustrated in Vergni et al. (2015) by proposing some more effective measures of both uncertainty and reliability. These measures are defined on the basis of CI expressed in years of return period. Second, the paper illustrates a possible method to extend this type of uncertainty analysis to nonparametric standardized indices.

The illustrative case study refers to SPI and SPEI indices, whose time series were calculated for a single station in central Italy. The results are discussed in order to analyze the variability of the proposed measures of uncertainty in relation to: sample size, type of approach (parametric or nonparametric), time scale, type of standardized index, and type of anomaly (excess or deficit).

The first step in the calculation of a generic standardized index S consists in the quantification, for each year i (i = 1, …, N), of the values,⁠, i.e., the values describing the reference variable X in relation to a particular period of interest j (typically the month), and to k-1 past consecutive months (k is called time scale). Then, a parametric or nonparametric probability function is fitted to homogenous (by j and k) samples of the values ⁠. The fitting is performed for each calendar month, in order to take into account the climatic differences due to seasonality. After this step, the cumulative probability of an observed value can be estimated and, by an equiprobability transformation (Abramowitz & Stegun 1965), it is possible to derive the corresponding standard normal deviate, which represents the standardized index S.

The standardized indices considered in this paper are the SPI (McKee et al. 1993) and the SPEI (Vicente-Serrano et al. 2010). For SPI, the reference hydrological variable X is the cumulative precipitation, while for SPEI, X is the climatic water balance (difference between precipitation and reference evapotranspiration ET0).

In the parametric calculation, the underlying distributions assumed for the cumulative precipitation (SPI) are both the two-parameter gamma (GAM) and the Pearson type III (PE3) distributions, since they are those most used in the literature (Blain 2011). For the climatic balance (SPEI), the analysis is based on both the generalized logistic (GLO) distribution and the generalized extreme value (GEV) distribution. The GLO was considered in the original algorithm proposed by Vicente-Serrano et al. (2010), but Stagge et al. (2015), in a recent analysis for Europe, showed that the goodness-of-fit of the GEV can be better than that of the GLO. The parameters of the distributions selected were estimated from L-moments (Hosking 1990). For both indices, each computation routine was performed assuming a unique probability distribution with parameters depending on the month of the year. As an alternative, it would be possible to assume a probability distribution dependent on the month of the year.

This type of nonparametric approach was preferred over that based on the empirical cumulative distribution functions (ECDF) (Farahmand & AghaKouchak 2015), which are less suitable for the calculation of the bootstrap CI. Concerning this, more details will be provided in the specific section.

Both SPI and the SPEI were analyzed at the 1-, 3-, 6- and 12-month time scales.

A second criterion was based on the normality test (Wu et al. 2007) applied to the resulting time series of the standardized indices. According to Wu et al. (2007), the series of a generic standardized index S are considered non-normal when the following three criteria are met: (1) Shapiro-Wilk statistic, W, is less than 0.96; (2) p-value associated to the W statistic is less than 0.1; and (3) the absolute value of the median S is greater than 0.05. The rationale behind this test is that standardized indices were designed to represent drought and wet anomalies in a similar way, and this requirement can be achieved only by normally distributed series. Non-normality frequently occurs when the distribution selected exhibits a poor fitting (Vergni et al. 2015).

Recently, Stagge et al. (2015) have applied the normality test to SPI and SPEI series computed on the basis of reanalysis data in Europe to assess the goodness-of-fit of several parametric distributions. They demonstrated that the normality test produces results similar to (or even more restrictive than) those obtained by the more traditional Kolmogorov–Smirnov and Anderson–Darling tests.

The CI for a generic standardized index S can be determined by using a bootstrap resampling technique. The theoretical details of the bootstrap approach can be found in Efron & Tibshirani (1993), and an example of its application in the estimation of the CI for parametric SPI is available in Hu et al. (2015). A brief description of the steps required for a generic index S is provided here in the following.

Input data consist of N values ⁠, i (i = 1,…, N) of the reference hydrological variable X for a given month j, time scale k, and location.

Then, M bootstrapped samples of size N with replacement from the original N values are generated. The further sequence of operations varies according to the type of calculation of the index S (parametric or nonparametric).

For parametric standardized indices, the steps are the following:

• (1) Choice of the probability distribution to be fitted to and estimation of the distribution parameters for each of the M bootstrapped samples, to obtain a family of M possible populations.

• (2) Calculation of M values of cumulative probability for each single value ⁠.

• (3) Transformation of the M values of in M values of the corresponding standardized index S through an equiprobability transformation (Abramowitz & Stegun 1965).

• (4) Analysis of the ECDF of the M values of S to identify the (α/2)^th and (1-α/2)^th percentiles, which are assumed, respectively, as the lower (S_min) and upper (S_max) limits of the CI for S with probability P= (1-α).

For nonparametric standardized indices, steps 1 and 2 are replaced by the following ones:

• (1) Choice of a nonparametric method for the estimation of the cumulative probability of the variable ⁠.

• (2) Calculation of M values of the nonparametric cumulative probability for each single value ⁠.

The other steps are analogous to those of the parametric approach but it is necessary to replace with ⁠.

For this type of application, we found that a kernel density function (Equation (1)) works better than an ECDF, in particular for the determination of the CI associated with the extreme observations. In fact, using an ECDF, the minimum and maximum observations have constant cumulative probabilities in all the bootstrapped samples, thus not allowing construction of their CI. Instead, the kernel density function provides a smoothed and continuous probability distribution whose characteristics (e.g., the bandwidth) are different in each bootstrapped sample.

All the analyses in this paper were based on CI with P = 90% and on M = 1,000.

The sizes of the CI associated with the point estimates of a standardized index have been used in this paper to derive some measures of uncertainty and reliability, whose definitions and calculation details are provided in the next section.

ΔT is not calculated for the events characterized by opposite upper and lower limits of the CI (i.e., S_min < 0 and S_max > 0). In fact, under these circumstances, the return periods associated with the lower and upper limits are not comparable (the former is the return period for the non-exceedance probability, the latter is the return period for the exceedance probability). However, the events having S_min < 0 and S_max > 0 correspond to near-normal conditions (S ≈ 0), for which the uncertainty is lower and the practical interest is limited.

As for ΔS, the greater the ΔT, the greater is the uncertainty. However, thanks to its unit of measure, ΔT is more informative than ΔS. Of course, it is not possible to define absolute ΔT thresholds, beyond which the estimate should be considered unacceptable or unreliable. Indeed, this depends on several factors, such as the return period associated with the point estimate S and the aims of the monitoring (i.e., the maximum admissible risk), as well as subjective considerations.

As for ΔT, T_ratio is also not calculated for events characterized by opposite upper and lower limits of the CI. The more T_ratio approaches 1, the lower is the uncertainty.

T_ratio is particularly suitable for evaluating the reliability of the index estimate. It can, in fact, be assumed that a certain estimate of the standardized index S is unreliable if the corresponding T_ratio is greater than a subjective threshold τ. In the case study illustrated in this paper, a τ = 3 was adopted. This relatively low (i.e., severe) threshold was chosen to obtain a significant number of unreliable estimates also under conditions characterized by the lowest overall uncertainty (i.e., for analysis based on very long time series of observations).

The original daily data were used in input to the LARS-WG 5.0 weather generator (Semenov & Barrow 2002) in order to obtain daily time series of different lengths (30, 60, 90, 120, 150, and 180 years), while maintaining the statistical properties of the original data. This was done in order to analyze the effect of the length of the available time series on the uncertainty associated with the standardized indices.

All the daily precipitation and temperature data were aggregated at the monthly time scale. Monthly estimates of reference evapotranspiration (required for SPEI calculation) were then obtained from monthly precipitation and monthly mean minimum and maximum temperature, by applying a modified version of the Hargreaves and Samani equation (Droogers & Allen 2002).

Months with zero rainfall are not very common (1.4%), and in most cases they occur in summer (34% and 21% in August and July, respectively). Zero rainfall values are not present at time scales greater than 2 months. As indicated by some authors (e.g., Wu et al. 2007), the presence of zero values is the primary reason for non-normal SPI series. However, in the authors' opinion, the low occurrence of zero values in the case study considered should not have significant consequences on the results.

Taking into account the number of indices (2), time scales (4), sample sizes (6) and calculation methods (3, of which, 1 nonparametric and 2 parametric with different underlying distributions), 144 different series of standardized indices were considered in the case study.

Table 1 summarizes the goodness-of-fit (in terms of mean RMSE) of nonparametric and parametric distributions for the case study, varying the hydrological variable, the time scale, and the sample size. For the sake of simplicity, only the two extreme sample sizes (30 and 180 years) are reported. The fitting was performed for each calendar month, therefore Table 1 shows average RMSE results.

From the results of Table 1, it is evident that the nonparametric approach provides a better goodness-of-fit than parametric approaches. This result is expected since the kernel density function is directly fitted to the sample observations, while in the parametric method, the sample characteristics are ‘filtered’ by the estimated distribution parameters. Of course, for sample sizes smaller than those here considered, the nonparametric approach could also lead to misleading probabilities estimations.

The comparison between the parametric distributions shows that the PE3 performs better than GAM for precipitation, probably due to the fact that PE3 has one more parameter than GAM. For the climatic balance, the performance of GEV is slightly better than that of GLO.

The 144 series of standardized indices were also analyzed by applying the normality test (Wu et al. 2007) for each calendar month. The normality assumption is not rejected for all the series derived from the nonparametric calculation, confirming the good performance of this type of approach. For the series derived from the parametric distributions, the overall percentages of rejection of the normality assumption are: 5.9% for GAM, 1.7% for PE3, 1.4% for GLO, and 1.7% for GEV. These values are consistent with those obtained by similar studies (Stagge et al. 2015) and confirm that all the selected distributions provide an adequate description of the reference variables underlying the SPI and SPEI indices. A slightly worse performance is obtained by the adoption of the GAM for SPI. For this distribution, a detailed analysis of the results by time scale (not shown) revealed that most of the rejections are associated with the 1-month time scale. By excluding this time scale, the performance of the GAM is similar to that obtained by the PE3 distribution. This lower performance of the GAM distribution for short accumulation periods is consistent with Stagge et al. (2015).

No relationship was found between the results of the normality test and the length of the time series or the calendar month.

Analyzing Figure 2, the following general comments can be made:

• (1) The point estimates (full lines) of the standardized indices obtained from different approaches are rather similar, particularly for the larger dataset.

• (2) As expected, whatever the calculation approach is, the uncertainty is much lower when the analysis is based on relatively larger datasets.

• (3) The uncertainty bounds obtained from different calculation approaches are almost similar for normal conditions and for moderate drought or wet conditions. As the severity of the events increases (e.g., values around the 7th and 61st months of the series), the uncertainty increases and also the differences among different approaches become more evident (particularly for the smaller dataset). This behavior has been already observed in Vergni et al. (2015) and it is due to the fact that the extreme values tend to intercept the tails of the bootstrapped probability functions, with a consequent explosion of the uncertainty.

This last result underlines the importance of an adequate sample size to obtain a reliable analysis of extremes. For such type of events, the inconsistency issues deriving from the choice of different suitable parametric distributions, affect, not only in the point estimates, but also in the corresponding measures of uncertainty.

The proposed uncertainty measures ΔT and T_ratio were calculated for each single value of the 144 series considered in the case study.

Both ΔT and T_ratio are calculated on the basis of the CI of standardized indices, of which an illustrative example is provided in Figure 2. As a consequence, also ΔT and T_ratio exhibit the already-discussed tendency to increase as the absolute value of the index S increases.

Figure 3 also indicates the advantage of both ΔT and T_ratio in comparison to ΔS. This last measure, in the range −1.5 < S < 1.5, exhibits a more or less constant and irreducible limit of about 0.4 (for all the combinations considered in the figure). The specific value (0.4) is mainly dependent on the sample size (90 years) to which Figure 3 is related. However, it can be observed that in the same range −1.5 < S < 1.5, the corresponding ΔT and T_ratio are close to 0 and 1, respectively, thus indicating a good reliability of the estimates of the indices. Therefore, a ΔS ≈ 0.4 (based on a 90-year time series) can be considered a low-uncertainty condition, but this interpretation has required the quantification of the corresponding ΔT or T_ratio. Moreover, apparently negligible increases of ΔS (for |S| > 2) can conceal dramatic increases of uncertainty when they are analyzed by the corresponding ΔT or T_ratio. This could lead to misleading interpretations of the uncertainty assessment based on ΔS.

These favorable practical aspects of ΔT and T_ratio in comparison with ΔS hold in general, irrespective of the location, calculation approach, sample size, distribution, and index type.

In the following sections a detailed analysis of the behavior of the measures ΔT and T_ratio, varying the calculation approach, the sample size, the time scale, the types of index, and anomaly, is carried out in relation to the case study considered.

Figure 4 clearly shows that the length of the time series is one of the most important factors in the regulation of the uncertainty (here measured in terms of ΔT). As expected, the uncertainty decreases as the length of the time series increases. However, this influence is evident up to a length of 90 years. For larger samples, this effect is less evident, and the variability of ΔT due to other factors becomes negligible. The most interesting differences arise from the analysis of ΔT for the smaller sample sizes (30 and 60 years).

First of all, it can be observed that the nonparametric approach usually leads to ΔT smaller and less variable (with the time scale and the index) than those obtained from parametric calculations. Moreover, the ΔT values obtained from the nonparametric approach are systematically lower for the drought than for the wet anomalies. This asymmetry in the CI derives from the fact that the variables considered are usually positively (right) skewed, thus a greater uncertainty is expected on the skew side (i.e., wet anomalies). This aspect is correctly and consistently captured only by the nonparametric approach, that, as also shown in Table 1, provides the best fitting to the data. The parametric distributions (of both SPI and SPEI) provide a correct description of this asymmetric uncertainty only for the 6- and 12-month time scales. Furthermore, the uncertainty of the parametric estimates related to the 30-year sample is much more variable than that related to the nonparametric estimates (this is particularly evident for both the 1-month and 12-month time scales).

As regards the comparison between parametric distributions, for SPI, the ΔT obtained from the PE3 and GAM are usually aligned, apart fron the 1-month time scale, for which the GAM (likely due to the lower goodness-of-fit) presents much higher uncertainty than PE3 for drought conditions. For SPEI, the GEV has, in most cases, slightly lower ΔT values than GLO for both drought and wet anomalies.

The comparative analysis illustrated in Figure 4 was also carried out on the basis of ΔS, instead of ΔT (not shown), with practically identical conclusions. Indeed, as previously explained, the flaw in ΔS does not lie in its capacity to quantify uncertainty, but rather in its practical interpretability.

As observed for ΔT, the percentage of UE also decreases as the length of the time series increases. This decrease is particularly marked when passing from the 30-year to the 60-year lengths, while for longer time series the influence of this factor is less evident (Figure 5). This result clearly demonstrates that the analyses based on only 30 years of observations (which is usually considered the minimum requirement for the calculation of standardized indices) are particularly critical from a reliability viewpoint. On the other hand, in this region of data availability, small increments in the number of observed data can lead to significant improvement of the estimate reliability.

The detailed analysis of the UE values obtained for SPI leads to comments similar to those presented in relation to ΔT. For SPEI, the results obtained for the nonparametric approach are aligned with those obtained by considering ΔT. For the parametric distributions, instead, the results are quite different from those obtained by analyzing ΔT: first of all it can be observed that also the parametric distributions often present the expected asymmetry (i.e., less reliability for the wet anomalies). Moreover, in most cases, the GLO has lower percentages of UE than GEV (with the sole exception of the wet anomalies for the 30-year time series). In this regard, it can be observed that from an applicative viewpoint, it is preferable to have lower UE percentages than lower mean ΔT values. Therefore, the analysis reveals that, for the case study examined, GLO is a better choice than GEV for SPEI calculation. This result is consistent with the recent study by Vicente-Serrano & Beguería (2016), which is in response to the analysis of Stagge et al. (2015).

Therefore, it can be observed that the uncertainty measures proposed in this paper (particularly the percentage of UE) can be a useful tool also in the selection of the most suitable distribution, which, with other characteristics being equal, should guarantee the lowest percentage of unreliable estimates.

In this paper, some bootstrap-based measures of uncertainty and reliability applicable to the standardized indices have been defined. These measures are based on the quantification of the size of the CI of each estimate in terms of years of return period. This allows one to obtain a clear and effective evaluation of the reliability of the index estimate. The application of the measures proposed has been illustrated in relation to the SPI and the SPEI, pointing out the contribution of different factors to the overall uncertainty. In this context, the uncertainty analysis can also be considered as a method, complementary to the goodness-of-fit, for selecting the calculation method most suitable for quantifying a standardized index. Indeed, in this choice, methods able to reduce the percentage of unreliable estimates (other factors being equal) should be preferred. Although the case study is based on a single weather station, some results can be considered valid in general: in particular, nonparametric methods, in the typical regions of data availability, allow one to obtain lower (and more accurate) uncertainty levels than parametric methods.

In the paper, only univariate standardized indices have been considered. The next step might be the extension of this type of uncertainty analysis to bivariate indices. Concerning this, it will be important to adopt a correct resampling technique of the multivariate data (in order to maintain the dependence structure among variables).

The present computational times of bootstrap technique are reasonable, and therefore this type of uncertainty analysis is advisable as a part of any research or application based on standardized indices.

This research was financially supported by Fondazione Cassa Risparmio Perugia, project code 2015.0350 021.