Inadequacy of the gamma distribution to calculate the Standardized Precipitation Index

Inadequação da distribuição gama para o cálculo do Índice Padronizado de Precipitação R E S U M O O Índice Padronizado de Precipitação foi desenvolvido para ser um método probabilístico que monitora os déficits de precipitação pluvial, de forma normalizada. Assim, a confiabilidade deste índice é afetada pelo uso de distribuições inadequadas aos dados de precipitação. Com isto, o objetivo do trabalho foi: avaliar o ajuste da distribuição gama aos totais de precipitação acumulados em distintas escalas temporais (Pelotas, Rio Grande do Sul), analisar o desempenho de demais distribuições no cálculo desse índice e avaliar o pressuposto de normalidade das séries desse quantificador de seca, calculadas com base nessas distintas distribuições. Os testes Lilliefors e de normalidade indicaram que a distribuição gama não é apropriada para o cálculo do índice em diversas escalas temporais. A distribuição generalizada normal mostrou-se apropriada. O uso desta função também resultou no maior número de séries com distribuição normal. Recomenda-se, então, que os sistemas de monitoramento de seca e os estudos acadêmicos reavaliem o uso da distribuição gama no cálculo do Índice Padronizado de Precipitação. Um código computacional capaz de calcular este índice de seca a partir da distribuição generalizada normal é também apresentado.


Introduction
Drought is a natural phenomenon that occurs in virtually all parts of the globe.It may be regarded as a slow-moving hazard that presents no single or universal definition (Hayes et al., 2011;Lloyde, 2014).However, it is well known that a prolonged and abnormally dry period can severely impact virtually all human and natural processes (Blain, 2012).Accordingly, a great number of drought indices have been proposed to monitor the drought conditions in several parts of the world.Among then, the Standardized Precipitation Index (SPI) has been widely used in both academic and operational modes (Wu et al., 2007;Jain et al., 2010) to detect rainfall deficits summed over several time spans (or timescales).
In Brazil, the SPI is used by drought early warning systems of several governmental agricultural institutions such as Empresa Brasileira de Pesquisa Agropecuária (EMBRAPA), Instituto Agronômico (IAC), Instituto Nacional de Meteorologia (INMET) as well as by the Instituto Nacional de Pesquisas Espaciais (INPE).Moreover, the Lincoln Declaration on Drought Indices stated that the national hydrological and meteorological services around the Globe should use the SPI to characterize meteorological droughts (Hayes et al., 2011).This declaration also stated that the World Meteorological Organization should encourage the widespread use of the SPI in countries wanting to track meteorological droughts.
The SPI may be calculated for several timescales depending on the particular interest of the users (Dutra et al., 2013) with typically values ranging from 1 to 24 months.Thus, the SPI is able to simultaneously describe wet conditions at a given timescale and dry conditions at other timescale (Türkeş & Tatli, 2009).However, the time responses to drought of different systems such as water reservoirs and/or soil moisture are frequently not known (Vicente-Serrano et al., 2012).Therefore, for each particular application and region a detailed study should be carried out to relate the different SPI values (calculated at different timescales) to the variable of interest (e.g.available soil moisture or natural reservoirs; Dutra et al., 2013).Further information on the advantages and drawbacks of the SPI and its association with the different types of drought can be found in several studies including Fernandes et al. (2010) and Blain (2012).
One of the fundamental steps of the SPI calculation algorithm is to estimate the cumulative probability of a given rainfall amount summed over a particular timescale.Therefore, the confidence in the SPI results may be affected by the use of a distribution that does not provide a suitable goodness-of-fit for the rainfall data (Wu et al., 2005;2007;Sienz et al., 2012;Stagge et al., 2015).Although authors such as Guttman (1999) and López-Moreno & Vicente-Serrano (2008) have already provided evidences indicating that the gamma distribution may not always be the most appropriate choice; this fundamental step is frequently based on this 2-parameter distribution (Wu et al., 2007;Dutra et al. 2013; among many others).Moreover, the gamma distribution may not give the best goodness-of-fit for the rainfall data summed over different timescales because it only has two free parameters (Wu et al., 2007).
The SPI final values are obtained by performing an equiprobability transformation of the cumulative probability estimated in the above-mentioned step (Wu et al., 2007).From a mathematical point of view, the essential goal of this latter step is to transform the SPI series into normally distributed series with zero mean and unitary variance.From a meteorological point of view, the essential goal of this transformation is to ensure that the SPI is capable of representing drought and flood events in a similar probabilistic fashion (Wu et al., 2007).Thus, the parametric distribution used to calculate this probabilitybased index should also be capable of providing normally distributed SPI series.Further details on the importance of meeting the normality assumption of the SPI series may be found in Wu et al. (2007), Blain (2012) and Stagge et al. (2015).
Finally, although McKee et al. (1993) have stated that the required length of record for rainfall used in the SPI calculation is 'ideally a continuous period of at least 30 years'; Guttman (1994) stated that 70 to 80 years of records are required for parameter estimation stability in the tails of the rainfall distributions.Depending on the distribution and on the method used to estimate its parameters, the required length of record may be even longer (Wu et al., 2005).
Considering the above-mentioned SPI features the following three questions are posed: Is the 2-parameter gamma distribution suitable for fitting the rainfall series summed over several timescales in a given location?Is there an alternative distribution capable of performing such task?Is this alternative distribution capable of providing a higher number of normally distributed SPI series than the gamma distribution does?Therefore, the goals of this study were: (i) to assess the fit of the gamma distribution for the rainfall amounts summed over several timescales (1 to 24 months; obtained from one of the oldest Brazilians weather stations -Pelotas, Rio Grande do Sul, 1896-2011); (ii) to assess the goodness-of-fit of alternative distributions to such rainfall series and; (iii) to evaluate the normality assumption of the SPI series obtained from the distributions used to calculated this drought index.

Material and Methods
The rainfall daily data used in this study were obtained from the weather station of Pelotas, Rio Grande do Sul State, Brazil (21º 52' S, 52º 21 W and 13.2 m; 1896-2011; Federal University of Pelotas, RS).The choice of Pelotas is because it is one of the longest meteorological series of Brazil showing a high consistency (Blain et al., 2009).From the daily data, it was derived the rainfall series used to calculate the SPI at several timescales (1 to 24 months).It is also emphasize that Pelotas is one of the fewest Brazilians weather station capable of meeting the required length of record for parameter estimation stability in the tails of the rainfall distributions (Guttman, 1994;Wu et al., 2005).All statistical tests were performed at the 0.05 significance level.
The SPI calculation starts by determining a cumulative density function (cdf) capable of properly describing the longterm observed rainfall series (Guttman, 1999).Once a pdf is chosen, the cumulative probability of a given rainfall amount [H(PRE)] is calculated from the following mixed distribution: where: N -sample size; and m -number of zeros of the data sample.
As described in Wu et al. (2007) the final step of the SPI calculation is to apply an equiprobability transformation to the [H(PRE)] values (Abramowitz & Stegun, 1965).
The Kolmogorov-Smirnov (KS) and the chi-square tests (χ 2 ) are frequently used to assess the fit of a given parametric distribution to a given data sample (Wilks, 2011).Nevertheless, to calculate the χ 2 one needs to divide the data into discrete classes.Consequently, the χ 2 test operates more naturally for discrete random variable (Wilks, 2011).
The KS test compares the cumulative density function (cdf) to the empirical function.Thus, for continuous data such as rainfall data, the KS test is often more powerful than the χ 2 (Wilks, 2011).However, an important feature of the KS test is that it must not be applied to the same data sample used to estimate the parameters of the distribution (Vlček & Huth, 2009).Therefore, since the parameters of distributions have been estimated from all available data, a modification of the original KS test referred to as Lilliefors test was used.
The Lilliefors test also compares the cumulative density function to the empirical function.However, while the null hypothesis of the KS test does not depend upon the explicit form of the distribution under analysis, the null hypothesis of the Lilliefors depends on the values of the shape parameter.Accordingly, the critical values for the Lilliefors test are usually derived from statistical simulations procedures (Shin et al., 2012;Heo et al., 2013;Blain, 2014).Further information on the Lilliefors test can be found in studies such as Blain (2014).
Finally, it was also took into account the domain of the probability distribution functions.Therefore, if an observed rainfall amount falls beyond the domain of a given distribution (described in Table 1), this distribution will be regarded as inappropriate to assess the probabilistic structure of the series from which this rainfall amount was observed.
As previously described, the parametric distribution used to calculate the SPI should be capable of providing normally distributed SPI series.Therefore, it was applied the same algorithm proposed by Wu et al. (2007) to verify the normality The parametric distributions used in this study, 2-parameter gamma, general extreme value, generalized logistic, 3-parameter generalized normal, generalized Pareto, Pearson Type III, 4-parameter Kappa and 5-parameter Wakeby are described in Table 1.The greek characters are the parameters of the distributions that were estimated from all available data by means of the L-moments method (Guttman, 1999;Queiroz & Chaudhry, 2006;Blain & Meschiatti, 2014).
# Standardized Precipitation Index Table 1.Parametric distributions used to calculated the SPI # (2) (3) assumption of the SPI series.A SPI series was considered non-normal when the following criteria were simultaneously met: (a) absolute value of the median greater than 0.05; (b) Shapiro-Wilk statistic lower than 0.96 and p-values ≤ 0.10.Further information on the Shapiro-Wilk test may be found in several studies including Razali & Wah (2011).

Results and Discussions
Several rainfall amounts observed at the weather station of Pelotas have fall beyond the domain of both Wakeby and Pearson Type III distribution.For instance, by considering the parameters of the Pearson Type III distribution obtained from the monthly series of October (µ = 100.64;ξ = 61.43 and σ = 1.33) and the domain of this distribution (Table 1; 8.26 mm ≤ PRE ≤ +∞) the probability of the observed rainfall amount for October 1924 (6.7 mm) cannot be estimated.Naturally, this prevents the Pearson Type III distribution to be used as a statistical basis for a drought early warning system.
In spite of the fact that the Wakeby distribution is able to exhibit shapes that other 2-or 3-parameter distributions cannot mimic (Sen & Niedzielski, 2010;Morgan et al., 2011;Santos et al., 2011), this 5-parameter function has also failed to describe the probabilistic structure of a large number of rainfall series observed at the weather station of Pelotas (37% of the 264 series).For instance, the observed rainfall amount for March 1905 is 1.0 mm.However, by considering the parameters of the Wakeby distribution for such Month (ξ = 11.22,α = 103.28,β = 1.63, g = 59.24, δ = -0.01) the domain of this function is 11.22 mm ≤ PRE ≤ ξ+α/β-g/δ.Therefore, the Pearson Type III and Wakeby have been discarded as alternative distributions to calculate the SPI.
Among the 264 rainfall series (Figure 1A; 12 Months and 24 timescales for each distribution) the Lilliefors test indicated that the gamma 2-parameter distribution cannot be used to describe the probabilistic structure of approximately 10% of the series (27 series).This result is consistent with the statement that a 2-parameter distribution may not give the best goodnessof-fit for the rainfall data summed over different time scales (Wu et al., 2007).This latter finding also agrees with the studies indicating that the gamma distribution may not be the most appropriate choice for assessing the probability of occurrence of the rainfall amounts used to calculate the SPI (Guttman, 1999;Lopez-Moreno & Vicente-Serrano, 2008).
The results presented in Figure 1 also indicate that 3-parameter distributions tend to be preferred to 2-parameter distributions (Guttman, 1999).For instance, the generalized extreme value distribution performed slightly better than the gamma distribution.This 3-parameter function cannot be used to assess the probability of the rainfall amounts obtained from 6.1% of the series (Figure 1B).The generalized normal and the generalized logistic were the only two distributions that have presented rejection rates lower than the adopted significance level.By considering that the significance level is the probability of falsely to reject the null hypothesis, these low rejection levels may be regarded as evidences in favor of the use of these two distributions (Stagge et al., 2015).In this view, the results presented in Figure 1C indicate that the generalized logistic distribution cannot be used to describe the probabilistic structure of 4.1% of the 264 series.Among all the evaluated distributions, the generalized normal distribution (Figure 1D) have presented the best performance.This distribution has failed to fit only 3.0%, of the 264 rainfall series.
The Kappa and General Pareto distributions have presented performances far worse than any other distribution evaluated in Figure 1 (not show here because of the sake of brevity).These distributions cannot be used to describe the probabilistic structure of more than 50% of the series.
After evaluating the results depicted in Figure 1, it becomes worth emphasizing that the analysis of the SPI literature revels that virtually all studies use a single distribution (frequently the 2-parameter gamma) to calculate this drought index (Guttman, 1999;Wu et al., 2007;Lopez-Moreno & Vicente-Serrano, 2008;Blain 2012;Dutra et al., 2013;among many others).In other words, the use of the SPI has been based on the assumption that, for a given location, a single distribution is always capable of fitting the rainfall amounts summed over several timescales.The results of the goodness-of-fit tests (Figure 1) have provided evidences indicating that this assumption may not always be met.
The normality test indicated that the gamma distribution has failed to provide normally distributed SPI series also for 10% of the analyzed series (Figure 2A).The results of the normality test are consistent with those presented in Figure 1 in the sense that the generalized normal distribution outperforms all other distributions (including the 2-parameter gamma) in providing normally distributed SPI series.The normality tests have revealed that the median of the SPI values calculated from generalized normal distribution fall more often within the [-0.05:+0.05]interval than the median of the SPI values obtained from the other distributions.Naturally, this is an important feature for an index developed to represent floods and drought events in a similar probabilistic fashion (Wu et al., 2007;Blain, 2012).
By using the r-package 'Lmonco' (Asquith, 2013) as well as the following r-code, one is able to calculate the SPI for any rainfall series.Naturally, the rainfall data should to have been drawn from a population with generalized normal distribution.Before using this r-code, the user should first write it in the r-script and the data must be sorted by month (all evaluated years must have 12 rainfall amounts).A basic knowledge of the r-software is also required.The code is easily adapted to calculate the SPI with other distributions and/or for higher numbers of timescales.

Conclusions
1.The results of this study have provided evidences indicating that the drought early warning systems as well as the academic studies should re-evaluate the widespread use of the 2-parameter gamma distribution in the Standardized Precipitation Index (SPI) calculation algorithm.
2. Based on one of the largest rainfall records of Brazil (the weather station of Pelotas, Rio Grande do Sul -1896Sul - -2011) ) this study also demonstrated that even for a particular location the assumption that a single distribution is always able to fit the rainfall amounts summed over several timescales may not be fulfilled.
3. The 2-parameter gamma distribution is not suitable for fitting the rainfall amounts summed over several timescales (1 to 24 months) and obtained from one of the longest rainfall records of Brazil (the weather station of Pelotas, Rio Grande do Sul).
4. Among all the probabilistic models evaluated from the rainfall data of Pelotas, Rio Grande do Sul -Brazil, the generalized normal distribution has presented the best performance in fitting the rainfall series.
5. The use of this 3-parameter distribution also resulted in the highest number of normally distributed SPI series.

Figure 1 .
Figure 1.Lilliefors test applied to the rainfall series used to calculate the Standardized Precipitation Index (SPI) for the weather station of Pelotas Rio Grande do Sul State, Brazil The dark squares denote that the distribution cannot be used to assess the probability of the rainfall amounts.The open squares indicate the opposite