DRY MONTHS IN THE AGRICULTURAL REGION OF RIBEIRÃO PRETO, STATE OF SÃO PAULO-BRAZIL: AN STUDY BASED ON THE EXTREME VALUE THEORY

The application of the Extreme Value Theory (EVT) to model the probability of occurrence of extreme low Standardized Precipitation Index (SPI) values leads to an increase of the knowledge related to the occurrence of extreme dry months. This sort of analysis can be carried out by means of two approaches: the block maxima (BM; associated with the General Extreme Value distribution) and the peaks-over-threshold (POT; associated with the Generalized Pareto distribution). Each of these procedures has its own advantages and drawbacks. Thus, the main goal of this study is to compare the performance of BM and POT in characterizing the probability of occurrence of extreme dry SPI values obtained from the weather station of Ribeirão Preto-SP (1937-2012). According to the goodness-of-fit tests, both BM and POT can be used to assess the probability of occurrence of the aforementioned extreme dry SPI monthly values. However, the scalar measures of accuracy and the return level plots indicate that POT provides the best fit distribution. The study also indicated that the uncertainties in the parameters estimates of a probabilistic model should be taken into account when the probability associated with a severe/extreme dry event is under analysis.


INTRODUCTION
It is well known that an extreme dry month can cause great damage to crop yield production.Strongly related to a deficit in the expected rainfall amounts, this departure from a climate normal condition frequently leads to a shortage of available water for plant growth.The Standardized Precipitation Index (SPI1 ) is a probability-based procedure developed for characterizing rainfall departures with regard to an expected local climate condition (GUTTMAN, 1999;BORDI et al., 2007;SANTOS et al., 2011 andBLAIN, 2012a, among others).Consequently, the SPI is suitable for detecting a shortage of water supply especially during those months in which high precipitation amounts are expected (BLAIN 2012a,b).As can be noted from the SPI calculation algorithm (BLAIN, 2012a), for those locations with a distinctly dry season, the lowest (or driest) SPI values tend to be observed during the rainy season.Moreover, according to BORDI et al. (2007) andHAYES et al. (2011), we may assume that the characterization of a dry event does not rely entirely on low precipitation amounts; it also depends on a negative anomaly of these amounts with respect to an expected local climate condition (BORDI et al., 2007).This last feature may be seen as one of the reasons that lead the SPI to be widely used by governmental agricultural institutions, such as Empresa Brasileira de Pesquisa Agropecuária (EMBRAPA), Instituto Agronômico (IAC) and Instituto Nacional de Meteorologia (INMET).
According to SANTOS et al. (2011) an important task in drought risk evaluation is to assess the rarity of severe dry events.Hence, assessing the probability of occurrence of extreme low SPI values leads to an increase of the knowledge related to the occurrence of extreme dry months.After BORDI et al. (2007) have stated that the SPI is more suitable than the rainfall for charactering extreme dry monthly events they used the Extreme Value Theory (EVT) to model the probability of extreme dry periods in Sicily-Italy.
A question that arises when an extreme value analysis is being carried out is related to what should be regarded as an extreme event (WILKS, 2011).A practical decision is to adopt the socalled block maxima or minima approach (BM), in which the blocks are considered equal to one year.The BM selects the largest (or the smallest) single value in each year or block (WILKS, 2011).
A desirable feature of the use of the BM is that it removes the influence of serial correlations on the estimates.However, this procedure also leads to a large loss of information (or data) by sampling only one value from each block (the average sample rate of the BM is equal to 1 event per year; λ=1).To overcome this undesirable feature one may adopt the procedure called peaks-overthreshold (POT).According to SUGAHARA et al. (2009) the abovementioned large loss of data has been an important reason for the increasing use of the POT.This last approach selects any value larger than (or smaller than) a threshold regardless of their year of occurrence (WILKS, 2011).The asymptotic distribution of a dataset that was correctly sampled by the POT is the Generalized Pareto Distribution (GPD) (COLES, 2001, SUGAHARA et al., 2009, KATZ, 2010and WILKS, 2011).Naturally, a dataset can only be correctly sampled by the POT if an appropriate threshold is chosen.
As described by COLES (2001), a too low threshold may violate the asymptotic basis of the function, leading to bias.On the other hand, a too high threshold may lead to a great loss of data resulting in high variance of the parameters estimations.SUGAHARA et al. (2009) state that choosing an appropriate threshold is a fundamental step for using the POT.At this point, we may pose the same question presented in WILKS (2011): Which of these approaches (BM or POT) is more useful for a given application?According to authors such as MADSEN et al. (1997) this question should be empirically evaluated for each location under analysis.
The aforementioned arguments, associated with the fact that a severe/extreme rainfall deficit frequently leads to crop yield reductions (BACK & UGGIONI, 2010), have motivated us to verify (i) whether the BM and/or POT can be used to assess the probability of occurrence of extreme dry SPI values in an agricultural region of the State of São Paulo, Brazil and (ii) which of these approaches provides the best characterization of the probability of occurrence of extreme dry months.Thus, the goals of this study were to apply both BM and POT to sample the extreme SPI data, obtained from the weather station of Ribeirão Preto , and to empirically compare the performance of these both approaches in providing the best fit extreme value distribution.It is expected that this study should provide a better assessment of the probability of occurrence of extreme dry months that may be linked to crop yield reductions.

MATERIAL AND METHODS
Monthly rainfall totals were obtained from the weather station of Ribeirão Preto (21º11'S, 47º48 'W, 620m -1937to 2012;State of São Paulo, Brazil).As can be noted from studies such as BLAIN (2009) the shape of these monthly probability functions vary from a strongly skewed distribution similar to those observed in arid climate (July) to a bell-shaped distribution similar to those observed in equatorial climate (January).All statistical tests of this study were performed at the 5% significance level.
As described by GUTTMAN (1999) the SPI calculation starts by determining a parametric distribution (cdf[x]) capable of properly describing the long-term observed precipitation.According to BLAIN (2011BLAIN ( , 2012b)), the rainfall monthly series of Ribeirão Preto can be considered as coming from a 2-parameter gamma distribution or a 3-parameter Pearson Type-III distribution (PE3).The cumulative probability [H(PRE)], associated with a given rainfall amount, is obtained from the following mixed distribution (equation 1) where m is the number of zeros observed in a dataset composed by N observations.
The final step of the SPI calculation algorithm is to subject H(PRE) to an equiprobabilistic transformation.This last step aims to transform H(PRE) into a normally distributed process with zero mean and unit variance.The SPI (dry) categories are shown in Table 1.Based on distinct locations of the State of São Paulo, BLAIN (2011) has indicated that the SPI series obtained from the PE3 meets more frequently the normality assumption (intrinsic to the use of this standardized index) than the SPI series obtained from the gamma.Therefore, we estimated the SPI monthly values by using the PE3 distribution (GUTTMAN, 1999, SANTOS et al., 2010, BLAIN, 2011and BLAIN, 2012b).
Finally, it is worth emphasizing that the SPI can be calculated for several time scales.However, as described in BORDI et al. (2007), time scales greater than the monthly scale introduces temporal persistence or serial correlations into its series.These authors state that the presence of a significant serial correlation is an undesirable feature for applying the EVT.Thus, we decided to estimate the SPI only at the monthly scale.The presence of serial correlation in the SPI monthly series was evaluated by using the autocorrelation function (acf; WILKS, 2011).The acf coefficients were estimated up to lag 24 (months).Although the number of lags (#Lags) is an arbitrary choice, the upper limit #Lags ≤ N/4 was respected.

Block maxima and Peaks-over-threshold
The Extremal Types Theorem implies that when the highest values of a process are stabilized with suitable sequences, they have a limiting distribution that (theoretically) must be one of the three types of the extreme value distributions, i.e.Gumbel (Type I), Fréchet (Type II) or Weibull (Type III).According to WILKS (2011) this theorem is also applicable to the distributions of extreme minima.However, in this case the sign of the data must be reversed.As described in COLES (2001), each one of these three different types of distributions gives distinct representations of extreme value behavior.Therefore, in some applications of the EVT, it is usual to adopt, at least, one of these types, and then to estimate the parameters of the adopted function.Naturally, a technique should be applied to select the most appropriate type for each data sample (COLES, 2001).Fortunately, these three types can be combined into a single function; the General Extreme Value Distribution (GEV; equation 2 in its cumulative form).By using the GEV function, the data themselves specify the most appropriate type of tail behavior and there is no need to make a priori judgments about which individual type (Gumbel, Fréchet or Weibull) should be adopted (COLES, 2001).
The parameters of the GEV distribution were estimated by using the method of maximum likelihood (ML) as described in several studies such as Coles (2001), WILKS (2011), FURIÓ & MENEU (2010) and BLAIN & CAMARGO (2012).Hereafter, these estimative will be represented by the following characters: μ (location), σ (scale), ξ (shape) and, for the sake of brevity, will be referred to as parameters.As described in studies such as DELGADO et al. (2010), while σ defines the spread of the distributions, μ defines the position of the GEV function with regard to the origin.The Gumbel distribution corresponds to the case ξ=0.The Fréchet and Weibull distributions correspond, respectively, to the cases ξ>0 and ξ<0 (COLES, 2001).According to CANNON (2010) ξ is the main determinant of the tail behavior.SUGAHARA et al. (2009) observed that several strategies have been proposed for selecting an appropriate threshold for using the GPD (equation 3, in its cumulative form).However, according to these authors, the scientific literature still lacks of a general and objective method to define such a threshold.In this study, the threshold was selected by means of two different procedures described in COLES ( 2001).The first one, called mean residual life plot, is an exploratory method that is based on the mean of the GPD in the sense that the expectation E(SPI-threshold|SPI>threshold) is a linear function of all valid thresholds.The second one is based on the search of stability of the parameters estimations obtained by adopting different thresholds.Further information regarding these two procedures can be found in COLES (2001).The parameters of the GPD distribution were also estimated by using the ML, which is frequently regarded as a robust parameter estimation procedure (HADDAD & RAHMAN, 2010).

Goodness-of-fit
The chi-square test (χ 2 ) and the Kolmogorov-Smirnov test (KS) are often used to check the fit of a given data sample to a parametric distribution (WILKS, 2011).The χ 2 is suitable for discrete random variables because to calculate it, the range of the data must be divided into discrete classes (WILKS, 2011).On the other hand, the KS test compares the theoretical and the empirical cumulative distribution.Hence, for continuous data such as the SPI, the KS test tends to be more powerful than the χ 2 (WILKS, 2011).The null hypothesis (H0) of the KS test states that the data was drawn from the theoretical distribution under analysis.However, according to authors such as WILKS (2011) and VLČEK & HUTH (2009) the original KS calculation algorithm is applicable if, and only if, the parameters of the theoretical distribution have not been estimated from the same sample used to evaluate the fit of the parametric distribution.
Given that the parameters of the extreme value distributions were fitted using all available data, the KS test had to be modified.This adapted method is frequently mentioned as Kolmogorov-Smirnov/Lilliefors test (KS-L; WILKS 2011).The statistical simulations required for its calculation were based on the method called "non-uniform random number generation by inversion" as described by WILKS (2011).The null distribution of the KS-L test were generated from Ns=10000 synthetic data samples.
Regarding extreme-value statistics, the Anderson-Darling (AD) test (equation 4) is based on the sum of the squares of the differences between theoretical and empirical distribution with a weight function [Ψ(SPI)] which emphasizes the discrepancies in both upper and lower tails (SHIN et al., 2012).This last feature cannot be observed in the KS-L algorithm.
According to SHIN et al. (2012) the AD test is a powerful test because it emphasizes the tails differences.However, one may correctly argue that equation 4 gives equal weights to both upper and lower tails of the distribution (SHIN et al., 2012).So, studies addressing drought events are naturally interested in the lower tail of the (SPI) distributions rather than in its upper tails.Thus, the use of a Ψ(SPI) function that, for the present case, gives emphasis to discrepancies at the lower tails of the distributions becomes an interesting choice.Considering this background, a modified form of the AD test (equation 5; AU) that gives emphasis to the lower tails of the distributions should also be used.As for the KS-L, the null hypothesis of the AU and AD were obtained by using the procedure called "non-uniform random number generation by inversion" (Ns=10000).
After analyzing equations 4 and 5, one may correctly argue that a weakness of the AD and AU tests is that both F(SPI) and G(SPI) tend to 1 as the SPI absolute values increase.Naturally, the same weakness can be associated with the use of the KS-L test.According to COLES (2001), this deficiency can be avoided by analyzing the quantile estimation, which can be obtained by inverting equations 2 (QGEV) and 3 (QGPD).The mean absolute error (MAE) and the mean squared error (MSE) were used to evaluate the agreement between these estimates (QGEV and QGPD) and the SPI values.Further information regarding the MAE and the MSE can be found in WILKS (2011).Finally, the return level plot, as described in several studies such as COLES ( 2001), was used as a qualitative assessment of the goodness-of-fit of the two parametric models.

RESULTS AND DISCUSSION
The acf function has indicated that the SPI monthly series can be seen as free from significant serial correlations.All acf coefficients fell within the white noise limits (Figure 1).This last result concurs with those found by BLAIN (2012b) for distinct weather stations of the State of São Paulo.This lack of serial correlation is also consistent with the great temporal variability of the rainfall totals observed in the State of São Paulo (BLAIN, 2009 andBLAIN &KAYANO, 2011) in the sense that a rainy month can be preceded and/or followed by a dry month.

Block maxima and Peaks-over-threshold
The mean residual life plot (not shown here for the sake of brevity) has indicated that the expectation E(SPI-threshold|SPI>threshold) may be seen as a linear function of threshold values equal to or smaller than -0.85.So, by adopting several thresholds values (raging from -2.0 to -0.5) it was observed that the parameters of equation 2 remains approximately constant for values equal to or smaller than -0.85.At this point, it is worth emphasizing that a threshold ≤ -0.85 leads to an average sample rate λ=2.10 events per year which is, naturally, greater than the sample rate obtained from the BM approach.Thus, by following the idea that a threshold must not result in the violation of the asymptotic basis of the GPD function and it also must avoid the great loss of data associated with the BM approach (COLES, 2001;BORDI et al., 2007;SUGAHARA et al., 2009 andWILKS, 2011), we adopted a threshold equal to -0.85, which is also virtually equal to the threshold (-0.84) used by SANTOS el al. (2011).
The parameters of the two distributions are: μ =1.213 [0.0688], σ=0.416 [0.0384], ξ=0.089 [0.0524] (GEV); σ*=0.459[0.0506], ξ=0.0636 [0.0762] (GPD); [.] represents the standard error associated with each estimate.Positive values of the shape parameter imply that both distributions have no upper limit (COLES, 2001;WILKS, 2011).Considering that the rainfall distributions, from which the SPI values were calculated, have a physical lower limit associated with the zero value, these positive values of ξ should be regarded as an indication that the real lower limit of the SPI distributions could not be specified.It also should be remembered that all SPI data were multiplied by -1.
Goodness-of-fit tests such as the KS-L, AD and AU provide empirical evidences for rejecting or accepting the null hypothesis that a given parametric distribution may be used to describe the probability of occurrence of a given variable.This last statement is based on the fact that these tests quantify the agreement between the empirical cumulative distribution and the theoretical cumulative function (WILKS 2011 andSHIN et al., 2012).
The abovementioned tests have indicated a good agreement between the empirical cumulative distributions and both theoretical distributions.The p-values associated with each outcome were always greater than the 5% critical level (Table 2).Thus, we may assume that the GEV as well as the GPD may be used to evaluate the probability of occurrence of the extreme dry SPI values obtained from the weather station of Ribeirão Preto.When looking at Table 2, one may also note that the lowest values of the MAE and MSE are obtained by using the GPD.At this point it is worth emphasizing that the accuracy with which the quantiles are estimated is a fundamental criterion for selecting a distribution (MADSEN et al., 1997).By considering that the MAE and MSE are scalar measures of agreement (WILKS, 2011) that evaluates the reliability of quantile estimates made by a given model (HADDAD & RAHMAN, 2011), we may assume that the GPD fits better than the GEV to the extreme dry SPI values obtained from the weather station of Ribeirão Preto.Considering the 95% confidence interval, the visual analysis of Figure 2 clearly indicates that the GEV model is inconsistent with the most extremes SPI values obtained from the weather station of Ribeirão Preto.However, the same analysis indicates that the GPD model is consistent with all observed data.Thus, similar to what was indicated by the MAE and MSE, the results depicted in Figure 2 indicate that the GPD is better than the GEV for describing the probability associated with extreme dries SPI values obtained from the abovementioned location.Moreover, the analysis of Figure 2 also indicates that parameter uncertainties, associated with a probabilistic model, should be taken into account when the probability of an extreme event is being assessed.According to COLES & PERICHI (2003) uncertainties in extreme value are substantial, so the sample variability has to be taken into account in any extreme value analysis.  .The black dots represent observed data.The dashed lines represent the 95% confidence interval.
As can be noted from Table 3 the most extreme dry monthly event occurred during the Month of April 2000.This event has a cumulative probability equal to 99,87%.No extreme dry SPI value was observed during the 1980's.The years of 1956 and 1963 present three months with extreme low SPI values.By considering the analyzed period , one may note that approximately 10% of the 76 years presents at least one extreme dry month that may cause great damage to the crop yield production.In this view, by using equation 3 and the parameters estimated in this study, it is possible to verify that the (average) cumulative probability of occurring at least one dry month Pr(x≥SPI=-2) during a given year in the location of Ribeirão Preto is approximately 92%.

CONCLUSSION
The peaks-over-threshold approach can be used to assess the probability of occurrence of extreme dry SPI monthly values obtained from the Weather Station of Ribeirão Preto -State of São Paulo, Brazil.
When compared to the Block Maxima, the peaks-over-threshold provides a better characterization of the probability of occurrence of extreme low SPI values that describe a shortage of available water for plant growth.Accordingly, we recommend the use of this latter approach to assess the probability of occurrence of an extreme dry month that can cause great damage to crop yield production.
The uncertainties in the parameters estimates should be taken into account when the probability associated with a severe/extreme dry event is being assessed.
FIGURE 1. Auto-correlation function (acf) applied to the monthly SPI series of the location of Ribeirão Preto (1937-2012), SP.

FIGURE 2 .
FIGURE 2. Return level plots for the Generalized Pareto Distribution (a) and Generalized Extreme Value Distribution (b) fitted to extreme and severe dry SPI values.Ribeirão Preto-SP.The black dots represent observed data.The dashed lines represent the 95% confidence interval.

TABLE 1 .
The SPI classification system.

TABLE 2 .
Kolmogorov-Smirnov/Lilliefors (KS-L) test, Anderson-Darling test (AD), Modified Anderson-Darling test (AU), mean absolute error (MAE) and mean squared error (MSE) used to evaluate the fit of the General Extreme Value distribution (GEV) and the Generalized Pareto distribution (GPD) to the extreme dry SPI values obtained from the weather station of Ribeirão Preto-SP.

TABLE 3 .
Probability of occurrence associated with extreme dry months observed for the weather station of Ribeirão Preto, State of SãoPaulo-Brazil (1937-2012)