Modelling Long-Term Monthly Rainfall Variability in Selected Provinces of South Africa: Trend and Extreme Value Analysis Approaches

: Extreme rainfall events have made signiﬁcant damages to properties, public infrastructure and agriculture in some provinces of South Africa notably in KwaZulu-Natal and Gauteng among others. The general global increase in the frequency and intensity of extreme precipitation events in recent years is raising a concern that human activities might be heavily disturbed. This study attempts to model long-term monthly rainfall variability in the selected provinces of South Africa using various statistical techniques. The study investigates the normality and stationarity of the underlying distribution of the whole body of rainfall data for each selected province, the long-term trends of the rainfall data and the extreme value distributions which model the tails of the rainfall distribution data. These approaches were meant to help achieve the broader purpose of this study of investigating the long-term rainfall trends, stationarity of the rainfall distributions and extreme value distributions of monthly rainfall records in the selected provinces of South Africa in this era of climate change. The ﬁve provinces considered in this study are Eastern Cape, Gauteng, KwaZulu-Natal, Limpopo and Mpumalanga. The ﬁndings revealed that the long-term rainfall distribution for all the selected provinces does not come from a normal distribution. Furthermore, the monthly rainfall data distribution for the majority of the provinces is not stationary. The paper discusses the modelling of monthly rainfall extremes using the non-stationary generalised extreme value distribution (GEVD) which falls under the block maxima extreme value theory (EVT) approach. The maximum likelihood estimation method was used to obtain the estimates of the parameters. The stationary GEVD was found as the best distribution model for Eastern Cape, Gauteng, and KwaZulu-Natal provinces. Furthermore, model ﬁtting supported non-stationary GEVD model for maximum monthly rainfall with nonlinear quadratic trend in the location parameter and a linear trend in the scale parameter for Limpopo, while in Mpumalanga the non-stationary GEVD model with a nonlinear quadratic trend in the scale parameter and no variation in the location parameter ﬁtted well to the monthly rainfall data. The negative values of the shape parameters for Eastern Cape and Mpumalanga suggest that the data follow the Weibull distribution class, while the positive values of the shape parameters for Gauteng, KwaZulu-Natal and Limpopo suggest that the data follow the Fréchet distribution class. The ﬁndings from this paper could give information that can assist decision makers establish strategies for proper planning of agriculture, infrastructure, drainage system and other water resource applications in the South African provinces. The study analysed the long-term trends of monthly rainfall data in the ﬁve selected provinces of South Africa from 1900 to 2017. Two trend analysis techniques were applied in this study, the Mann-Kendall test and Sen’s slope. Findings from the Mann-Kendall test revealed statistically signiﬁcant monotonic decreasing trends in Eastern Cape, Gauteng and Kwazulu-Natal provinces, while in Limpopo and Mpumalanga provinces the trends were also revealed to be monotonically decreasing, but insigniﬁcant. The Mann-Kendall test statistic ﬁndings for Eastern Cape, Gauteng, KwaZulu-Natal and Limpopo were in agreement with the ﬁndings from Sen’s slope estimator method. However, the Mann-Kendall test ﬁndings for Mpumalanga slightly differed from Sen’s slope estimator ﬁndings. This slight difference can be interpreted as a conﬁrmation of the insigniﬁcance of long-term trend and slope for Mpumalanga monthly rainfall. The study further analysed the modelling of monthly rainfall extremes using the non-stationary GEVD approach which belongs to the block maxima realisation The maximum likelihood estimation method was used to obtain the estimates of the parameters. The stationary GEVD was found as the best distribution model for Eastern Cape, Gauteng and KwaZulu-Natal provinces. Furthermore, model ﬁtting supported non-stationary GEVD models for Limpopo and Mpumalanga maximum monthly rainfall, with nonlinear quadratic trend in the location parameter and a linear trend in the scale parameter for Limpopo, while for Mpumalanga the non-stationary GEVD model with a nonlinear quadratic trend in the scale parameter and no variation in the location parameter ﬁtted well to the maximum monthly rainfall data. The study further revealed that the maximum monthly rainfall for Eastern Cape and Mpumalanga can be modelled by distributions in the negative-Weibull domain, while maximum monthly rainfall data for Gauteng, KwaZulu-Natal and Limpopo follow distributions in the Fréchet distribution class.


Introduction
According to [1], climate change is possibly the biggest environmental problem facing the globe. Masereka et al. [2] stated that flood risks are caused by extreme rainfall events that have resulted in flood disasters that accounted for about 47% of all weather-related

Test for Stationarity
Statistical theory offers a wide range of unit root tests, with the most commonly used being augmented Dickey-Fuller (ADF) test, Phillips-Perron (PP) test and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test [16]. In this study ADF, PP and KPSS are used to test whether the monthly rainfall data for selected provinces of South Africa are stationary.

Augmented Dickey-Fuller (ADF) Test
The ADF test was employed in this study to check whether the monthly rainfall data for selected provinces of South Africa are stationary.
The ADF test is assessed under the following hypotheses: There exists a unit root and the time series is non-stationary, H 1 : The time series is stationary.
The ADF test consists of estimating the following regression model: where β is a constant, β 1 is the coefficient on time trend. The null hypothesis is δ = 1, and the alternative hypothesis is δ = 1, while t is a pure white noise error term and the ADF follows an asymptotic distribution [17].

Phillips-Perron (PP) Unit Root Test
The Phillips-Perron test is a more developed test, introduced in 1988 and it has the same null hypothesis with ADF test and also uses the same critical values with it [18,19] . The PP test makes a non-parametric correction to the t-statistic. The PP test involves the equation coming from Dickey-Fuller test: where t is I(0) and it can be heteroscedastic. For this reason, the test estimates the equation: The PP method estimates the non-augmented DF test equation and modifies the t-ratio of the coefficient, so that serial correlation does not affect the asymptotic distribution of the test statistic. The PP test is based on the statistic: The PP test is assessed under the following hypotheses: There is a unit root, H 1 : There is no unit root.

Kwiatkowski-Phillips-Schmidt-Shin (KPSS)
The ADF and PP test mentioned in the previous sections are testing the null hypothesis that the time series y t is integrated of order one, I(1). The opposite case, that is, testing the null hypothesis that the time series y t is I(0) is described by the KPSS test [20]. KPSS builts on the idea that the time series is stationary around a deterministic trend and is calculated as the sum of deterministic trend, random walk and stationary random error. It is based on the model: where d t = ∑ p i=o β i t i , for p = 0, 1, contains deterministic parts of the model constant or deterministic trend, t are independent and identically distributed (iid) error terms ∼ N(0, σ 2 ), and r t is a random walk with variance σ 2 u and u t . KPSS test is based on the likelihood method test of the hypothesis that random walk has a zero variance, i.e., H 0 :σ 2 u = 0, which means that r t is a constant, against the alternative H 1 :σ 2 u > 0. The test statistic is written as: where s t = ∑ T t=1ˆ , t = 1, 2, ..., T, andσ 2 is the estimate of variance σ 2 of process t from equation (5). Critical values are derived by a simulation method and are listed in [20]. The advantage of the KPSS is that to some extent KPSS alleviates the problem that is present with the ADF test [21].
Kwiatkowski et al. [20] argue that KPSS test can differentiate a series that appears to be stationary, series that appears to have a unit root, and series for which the data are not sufficiently informative to be sure whether they are stationary or integrated.
The KPSS test is assessed under the following hypotheses: H 0 : The series does not have a unit root or is stationary, H 1 : The series has a unit root or is not stationary.

Trend Test
This study used non-parametric Mann-Kendall (M-K) test statistic, Sen's slope estimator and time series plots to investigate the long-term trend of the monthly rainfall and its variability across the selected provinces.

Non-Parametric Mann-Kendall (M-K) Test Statistic
The non-parametric Mann-Kendall (M-K) test statistic is frequently used to quantify the significance of monotonic trend in hydrometeorological time series [22,23]. The M-K test statistic is defined as sgn(e i − e j ), (7) where n is the number of extreme values. If S is positive, then there is an increasing trend, but if S is negative, then there is a decreasing trend, and sgn(e i − e j ) is a sign function given by: Under the null hypothesis of no trend, the theoretical mean of S is 0 and its variance is given by Var(S) = n(n − 1)(2n + 5) − g ∑ p=1 t p (t p − 1)(2t p + 5) /18, (9) where g is the number of tied groups (a tied group is a set of sample data having the same value), and t p is the number of data points in the pth tied group. If no tied group exist, this process can be ignored [23]. In cases where the sample size n > 30, the normalised test statistic Z can be used to statistically quantify the significance of the trend. Z is calculated using the following equation: Positive values of Z indicate an increasing trend, while negative Z values show decreasing trends. In a one-tailed test at a significance level of α, the null hypothesis of no trend is rejected if | Z | > z α , where z is the standard normal variable. In this study, the significance level was set to be 5%.

Sen's Slope Estimator
Sen's slope estimator non-parametric method was used to estimate the magnitude of trends in the time series data [24]. The slope of "n" pairs of data can be first estimated by using the following equation: In this equation, X j and X k denote data values at time j and k, respectively, and time j is after time k (k ≤ j). The median of "n" values of β i is the Sen's slope estimator test. A negative β i value represents a decreasing trend, a positive β i value represents an increasing trend over time.
If "n" is an even number, then the Sen's slope estimator is computed by using the following equation: If "n" is an odd number, then the estimated slope by using the Sen's slope method can be computed as follows: β med = (β [n+/2] ). (13) β med is tested by a two tailed test at 100 (1−α) % confidence level, and the true slope of monotonic trend can be estimated by using a non-parametric test [25,26].

Time Series Plots
A time series plot is simply a graph in which the data values are arranged sequentially in time. It is commonly used to give a pictorial view of the data series over time.

Test for Normality
According to [27], there are several parametric and non-parametric methods of assessing whether data are normally distributed or not. These methods can be split into two groups: graphical and statistical. The most frequently used methods include: Quantilequantile (Q-Q) plots, density plots, probability-probability (P-P) plots, Anderson-Darling test (AD), Shapiro-Wilk (SW) test, D'Agostino-Pearson K2 (DPK) test, chi-square test, Jarque-Bera (JB) test, kurtosis test, Shapiro-Francia (SF), skewness test, robust Jarque-Bera (RJB) test among others. In this study the JB, SW and chi-square methods are employed to check whether the monthly rainfall data are normally distributed. The SW test is one of the most popular tests for normality assumption diagnostics, and has good properties of power based on correlation within given observations and associated normal scores [28]. The JB and chi-square tests are among the most widely used techniques for testing normality of the data.

Jarque-Bera (JB) Test
The JB test statistic is expressed as: where √ b 1 and b 2 are the skewness and kurtosis measures and are given by m 3 (m 2 ) 3/2 and m 4 (m 2 ) 3 , respectively; and m 2 , m 3 and m 4 , are second, third and fourth central moments, respectively. The JB test statistic is chi-square distributed with two degrees of freedom. The hypothesis test for the JB test procedure is H 0 : The monthly rainfall data is normally distributed, versus, H 1 : The monthly rainfall data do not come from a normal distribution.

Shapiro-Wilk (SW) Test
The SW test is of the form: where m = n 2 if n is even, while m = (n−1) 2 and x (i) represent the i th order statistic of a sample. The constants a i are given by: (a 1 , a 2 , ..., a n ) = and m is given by m = (m 1 , m 2 , ..., m n ) T , where m 1 , m 2 , ..., m n are the expected values of order statistics of iid random variables sampled from the standard normal distribution, and V is the covariance matrix of those order statistics [27].
The SW test is assessed under the following hypotheses: The monthly rainfall data is normally distributed, H 1 : The monthly rainfall data does not come from a normal distribution.

Chi-Square Test
The chi-square goodness-of-fit test is defined as: where (O i ) and (E i ) refer to the ith observed and expected frequencies, respectively, and n is the number of groups. When the null hypothesis is true, the above test statistic follows a chi-square distribution with k − 1 degrees of freedom [27].
The chi-square test is assessed under the following hypotheses: The monthly rainfall data are sampled from a normal distribution, H 1 : The monthly rainfall data are not sampled from a normal distribution.

Extreme Value Theory Techniques
In extreme value theory (EVT) two approaches exist: the block maxima (BM) and the peaks-over-threshold (POT) methods. According to [29], the BM is an approach in EVT that consists of dividing the observation period into non-overlapping periods of equal sizes. The current study utilises the BM approach in a changing climate to model monthly rainfall of the five selected provinces of South Africa.

Stationary Generalised Extreme Value Distribution
Generalised extreme value distribution (GEVD) is the family of asymptotic distribution that describes the behaviour of extreme conditions. The GEVD consists of three extreme value distributions namely: Gumbel, Fréchet and Weibull families which are also referred to as type I, II and III extreme value distributions [30][31][32]. The cumulative probability distribution for GEVD is of the form: where x are the extreme values from the blocks, µ, σ and ξ are the location, scale and shape parameters, respectively. For ξ > 0, we obtain the Fréchet distribution, for ξ = 0, we get the Gumbel distribution and for ξ < 0, we get the Weibull distribution.

Non-Stationary Generalised Extreme Value Distribution
The non-stationary GEVD model is the fundamental modification of the stationary GEVD model [30]. To account for non-stationary GEVD, the location parameter µ and the scale parameter σ are assumed to vary with time t and possibly other covariates [32,33]. The non-stationary GEVD is given by: In the simplest case, the following regression structures could be examined for the location and scale parameters: allowing up to quadratic dependence on time t and keeping the shape parameter constant [34].

Parameter Estimation of Non-Stationary GEVD
Parameters of the non-stationary GEVD are estimated using the method of maximum likelihood (ML). For a sample of N observations, the ML of the time-dependent GEVD in (18) was determined by maximising the log-likelihood function, expressed with timevarying parameters: where N is the number of years of observation. To obtain the GEVD parameter estimators that maximise Equation (18) we use the interior algorithm based nonlinear optimisation in the MATLAB Optimisation Toolbox [22].

Goodness-of-Fit
Goodness-of-fit test statistics are used for checking the validity of a specified or assumed probability distribution model. In this study, Kolmogorov-Smirnov (K-S) test, Anderson-Darling (A-D) and graphical methods, were applied to identify the best model.

Kolmogorov-Smirnov (K-S) Test
The K-S test, based on the empirical cumulative distribution function is used to decide if a sample comes from a hypothesised continuous distribution [35][36][37]. The K-S statistic D is defined as the largest vertical distance between theoretical and the empirical cumulative distribution (CDF) and is formulated as follows: where X i are random samples, i = 1, 2, ..., n, and the CDF is The K-S test is estimated under the following hypotheses: The monthly rainfall data follow a specified distribution, H 1 : The monthly rainfall data do not follow the specified distribution.

Anderson-Darling (A-D) Test
The A-D test statistic (A 2 ) is defined as: The A-D test is used to compare the fit of an observed CDF to an expected CDF. This test gives more weight to the tails of the distribution than the K-S test [36,37].
The A-D test is estimated under the following hypotheses: H 0 : The monthly rainfall data follow a specified distribution, H 1 : The monthly rainfall data do not follow the specified distribution.

Graphical Test
Alam et al. [35] stated that graphical test is one of the most simple powerful techniques for selecting the best-fit model. To check if the time-dependent GEVD fit well to the monthly rainfall data, the following graphical tests were used.
Quantile-quantile (Q-Q) plots Quantile-quantile (Q-Q) plot, is a comparison of an empirical form for estimating the exceedance and the inverse of fitted distribution function. Any departure from linearity indicates model failure in perfectly fitting the data [38].
Probability-probability (P-P) plots Probability-probability (P-P) plot is a comparison of an empirical (usually percentage rank) and the fitted distribution function. In case of perfect fit, the data would line up on the diagonal of the probability plots [35,38].

Return level plots
In these plots the empirical estimates of the return level functions are added. If there is an agreement between the model-based curve and empirical estimates, then the model is suitable for the data [35,38].

Choice of Preferred Model
When time-dependent GEVD is considered with covariates, there are a number of possible models to select from [39]. In order to select between model fits, a test of the likelihood ratio test also known as the deviance (D) statistic is used. For models M 0 ⊂ M i , we define the D statistic as: where l 0 (M 0 ) and l i (M i ) are the maximised log-likelihood under models M 0 and M i , respectively. The asymptotic distribution of D is given by χ 2 k distribution with k degrees of freedom, where k is the difference in dimensionality of M i and M 0 . The calculated deviance statistic, D, is compared to critical values from χ 2 k at α level of significance. Large values of D suggest that M i explains substantially more of the variation in the data than M 0 [32] and [39].

Exploratory Data Analysis
This section is divided into three sections: descriptive statistics, stationarity tests and normality tests.

Descriptive Statistics
The descriptive statistics evaluated are the mean, standard deviation, median, kurtosis, skewness, minimum and the maximum monthly rainfall amount for each province. The summary of the descriptive statistics for each province is presented in Table 1.
From Table 1, the monthly rainfall data for each province has a mean valueX > Q 2 (Median), indicating that the monthly rainfall data is positively skewed and this is confirmed by the positive values of skewness. Eastern Cape, Gauteng, KwaZulu-Natal provinces have kurtosis greater than three which indicate heavy tails than a normal distribution, while Limpopo and Mpumalanga have kurtosis less than three which indicate lighter tails than a normal distribution.
The standard deviation for all the five provinces ranges from 31.28 to 57.23 mm per month. KwaZulu-Natal province has the highest standard deviation with the value of 57.23 mm per month which indicates a large variation in the monthly rainfall series, while Mpumalanga province has the lowest standard deviation of 31.28 mm per month which implies a small variation in the monthly rainfall series.
The minimum monthly rainfall ranges between 0.01 mm and 0.50 mm per month where Eastern Cape receives the highest minimum rainfall of 0.50 mm per month, while Gauteng and KwaZulu-Natal receive the lowest minimum rainfall of 0.01 mm per month.
The maximum monthly rainfall lie between 111.00 mm and 478.80 mm per month where KwaZulu-Natal receives the highest maximum monthly rainfall of 478.80 mm per month followed by Gauteng with the maximum rainfall of 438.10 mm per month. Mpumalanga receives the lowest maximum rainfall of 111.00 mm per month.

Test for Stationarity Results
The augmented Dickey-Fuller(ADF), Phillips-Perron (PP) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests were used to check for stationarity of monthly rainfall data for selected provinces of South Africa. Table 2 shows the results of the ADF, PP and KPSS tests.
The ADF and PP tests were assessed under the following hypotheses: There exists a unit root and the time series is non-stationary, H 1 : The time series is stationary.
The KPSS test was tested under the following hypotheses: The series does not have a unit root test (or series is stationary). H 1 : The series has a unit root (or series is not stationary).
From Table 2 the p-values of the ADF test statistics for Eastern Cape, Limpopo and Mpumalanga are significant (p < 0.05), suggesting that the monthly rainfall data for these three provinces are stationary. The ADF p-values for Gauteng and KwaZulu-Natal are insignificant (p > 0.05), suggesting that the monthly rainfall data for these two provinces are not stationary at 5% level of significance.
Also, from Table 2 the p-values of the KPSS test for all five provinces are significant (p < 0.05), suggesting that the monthly rainfall data are not stationary. Furthermore, from Table 2 the p-values of the PP test for all five provinces are significant (p < 0.05), implying that the monthly rainfall data are stationary.
Overall, based on all the stationarity test findings, we conclude that the monthly rainfall data are not stationary for the majority of the provinces.

Test for Normality Results
In this study we formally tested for normality of the monthly rainfall data using the Jarque-Bera (JB), Shapiro-Wilk (SW) and chi-square tests. Table 3, presents the results of the normality tests.
The JB, SW and chi-square tests are assessed under the following hypotheses: The monthly rainfall data are normally distributed, versus, H 1 : The monthly rainfall data do not come from a normal distribution.
From Table 3, the results for all the three normality tests are significant (p < 0.05), which suggests that the monthly rainfall data for all the five provinces do not come from a normal distribution.

Results and Discussion
This section is divided into two sections namely; trend analysis and model fitting.

Trend Analysis Results
Mann-Kendall test statistic, Sen's slope estimator and time series plots were used to analyse the long-term trends of the monthly rainfall data for the five provinces. The Mann-Kendall test statistic and Sen'slope results are presented in Table 4. The outcome of the Mann-Kendall test results revealed that in Eastern Cape, Gauteng and KwaZulu-Natal provinces there were significant monotonic decreasing long-term trends (p < 0.05 and τ negative), while in Limpopo and Mpumalanga there were no significant monotonic decreasing long-term trends (p > 0.05 and τ negative). Sen's slope values for Eastern Cape, Gauteng and KwaZulu-Natal showed significant decreasing magnitudes of trends, which were corresponding with the Mann-Kendall test results. While in Limpopo, Sens' slope value revealed an insignificant decreasing magnitude of trend, which also supports the findings from Mann-Kendall test. On the other hand, Sen's slope value for Mpumalanga showed no magnitude of trend, which slightly differs from the results of Mann-Kendall test. The latter findings illustrate the insignificance of the decreasing monotonic trend for Mpumalanga. Figure 1 illustrate the monthly rainfall data time series plots for Eastern Cape, Gauteng, KwaZulu-Natal, Limpopo and Mpumalanga provinces. The time series plots in Figure 1 do not exhibit any significant discernible long-term trends for all the provinces. This justifies the use of Mann-Kendall test to help uncover the hidden long-term trends in the monthly rainfall series in Table 4.

Non-Stationary GEVD Modelling of Annual Block Maxima Rainfall Data
The time series plots of the annual block maxima rainfall series are shown in Figure 2. There seems to be some strong evidence for a positive long-term trend over the years, for all the provinces. A substantial part of the variability in the data can probably be explained by a systematic variation in rainfall over the years. One way of capturing this trend is by allowing the GEVD location and scale parameters to vary with time [40]. From Figure 2 a simple linear trend in time seems plausible for our annual maximum rainfall X t , and we can use the model where µ(t) and σ(t) are the time-dependent location and scale parameters, respectively. In the present study, eight models are proposed for the non-stationary GEVD: M 1 , M 2 , M 3 , M 4 , M 5 , M 6 , M 7 and M 8 . The reference model is denoted by M 0 and is the stationary GEVD [40]. Model M 1 has a linear trend in the location parameter such that µ(t) = µ 0 + µ 1 t, σ(t) = σ and ξ(t) = ξ; Model M 2 has a linear trend in the scale parameter such that µ(t) = µ, log σ(t) = exp(σ 0 + σ 1 t) and ξ(t) = ξ; Model M 3 has a linear trend in both location and scale parameters such that µ(t) = µ 0 + µ 1 t, log σ(t) = exp(σ 0 + σ 1 t) and ξ(t) = ξ; Model M 4 has a nonlinear quadratic trend in the location parameter and a linear trend in scale parameter such that µ(t) = µ 0 + µ 1 t + µ 2 t 2 , log σ(t) = exp(σ + σ 1 t) and ξ(t) = ξ; Model M 5 has a linear trend in the location parameter and a nonlinear quadratic trend in the scale parameter such that µ(t) = µ 0 + µ 1 t, log σ(t) = exp(σ 0 + σ 1 t + σ 2 t 2 ) and ξ(t) = ξ; Model M 6 has a nonlinear quadratic trend in both location and scale parameters such that µ(t) = µ 0 + µ 1 t + µ 2 t 2 , log σ(t) = exp(σ 0 + σ 1 t + σ 2 t 2 ) and ξ(t) = ξ; Model M 7 has a nonlinear quadratic trend in the location parameter with no variation in scale such that µ(t) = µ 0 + µ 1 t + µ 2 t 2 , σ(t) = σ and ξ(t) = ξ; Model M 8 has a nonlinear quadratic trend in the scale parameter with no variation in the location parameter such that µ(t) = µ, log σ(t) = exp(σ 0 + σ 1 t + σ 2 t 2 ) and ξ(t) = ξ.

Eastern Cape
The stationary GEVD model for Eastern Cape data (i.e., model M 0 ) has a maximum negative log-likelihood (NLLH) of 556.765 (see Table 5). A GEVD model with linear trend in the location parameter (i.e., M 1 ) has a maximum NLLH of 555.820. The deviance statistic for comparing these two models is therefore, D = 2(556.769 − 555.820) = 1.898, which is small compared to χ 2 1 (0.05) = 3.841. Thus, allowing for a linear trend in the location parameter does not improve on our stationary GEVD model, M 0 . Therefore, M 1 is not a worth model to consider.
Consider the pair of models (M 0 , M 2 ) from Table 5. The deviance statistic is 2(556.769 − 555.724) = 2.090, which is small compared to χ 2 1 (0.05) = 3.841. Thus, allowing for a linear trend in the scale parameter does not improve on our stationary GEVD model, therefore, we reject model M 2 and conclude that is not worthwhile to allow for a linear trend in the scale parameter.
From Table 5 Table 5, the model pair (M 0 , M 8 ), which allows for nonlinear quadratic trend in scale parameter with no variation in location parameter, has a deviance statistic of 0.354, which is too small compared to the critical value of 5.991 with 2 degrees of freedom. Thus, allowing for a quadratic trend in the scale parameter with no variation in the location parameter does not improve on the stationary GEVD. Overall, the final model for Eastern Cape is the stationary GEVD model, M 0 . The general model for Eastern Cape is given by The shape parameter (−0.012) for the model, M 0 , in (28) indicates that the rainfall data for Eastern Cape can be modelled by the Weibull class of distributions since the shape parameter ξ < 0. The diagnostic plots for the stationary GEVD model in (28) are presented in Figure 3. The diagnostic plot results in Figure 3 show that the stationary GEVD model, M 0 , is the best fit for the Eastern Cape monthly rainfall data.

Goodness-of-fit test for Eastern Cape GEVD model
The goodness-of-fit test based on Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests were performed in order to check if the maximum monthly rainfall data for Eastern Cape follow a stationary GEVD model. Table 6 presents the results of the K-S and A-D goodness-of-fit tests for the selected stationary GEVD model for Eastern Cape.
The hypotheses are formulated as follows H 0 : The monthly rainfall data follow a specified distribution, and H 1 : The monthly rainfall data do not follow the specified distribution.
Since the p-values for both the K-S and A-D tests are greater than the 5% level of significance, α = 0.05, we conclude that the maximum monthly rainfall for Eastern Cape follow the specified stationary GEVD. The model pairs (M 0 , M 1 ) and (M 0 , M 2 ) from Table 7 have the same critical value of χ 2 1 (0.05) = 3.841 with the deviance statistic values of 0.022 and 0.250 for M 1 and M 2 , respectively. Since the values of the deviance statistics for M 1 (0.022) and M 2 (0.250) are smaller than the critical value of 3.841, we conclude that both models do not provide any improvement in fit over the stationary GEVD model. From Table 7 The shape parameter (0.117) for the stationary GEVD model, M 0 , in (29) indicates that the rainfall data for Gauteng can be modelled using Fréchet distribution class since the shape parameter ξ > 0. The diagnostic plots for the stationary GEVD model in (29) are presented in Figure 4. The diagnostic plot results in Figure 4 reveal that the stationary GEVD model, M 0 , in the Fréchet domain of attraction is the best fit for Gauteng monthly rainfall data.

Goodness-of-fit test for Gauteng GEVD model
Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests were used to determine whether maximum monthly rainfall data for Gauteng follow a stationary GEVD. Table 8 presents the results of the K-S and A-D goodness-of-fit tests for Gauteng stationary GEVD model.
The results from Table 8 show that the p-values for both the K-S and A-D tests are not significant (p > 0.05). Therefore, we conclude that the maximum monthly rainfall for Gauteng province follow the specified stationary GEVD.   Table 9 share a critical value of χ 2 2 (0.05) = 5.991 with deviance statistic values of 2.248 and −0.176 for M 7 and M 8 , respectively. Since the values of the deviance statistics are smaller than the critical value of 5.991 with 2 degrees of freedom, it implies that both models do not provide any improvement in fit over the stationary GEVD model.
The model pair (M 0 , M 6 ), which allows for nonlinear quadratic trend in both the location and scale parameters in Table 9, has a deviance statistic of 0.598 which is too small compared to the critical value of 9.488 with 4 degrees of freedom. Thus, allowing for a quadratic trend in both the location and scale parameters does not improve on the stationary GEVD model. Overall, the final best model for KwaZulu-Natal is the stationary GEVD model, M 0 . The general model for KwaZulu-Natal is given by The shape parameter (0.070) for the model M 0 , in (30) suggests that the rainfall data for KwaZulu-Natal can be modelled using Fréchet class of distributions since the shape parameter ξ > 0. The diagnostic plots for the stationary GEVD model in (30) are presented in Figure 5. The results in Figure 5 show that the stationary GEVD model, M 0 , is the best fit for KwaZulu-Natal maximum monthly rainfall data since all the four diagnostic plots suggest a reasonable good fit.

Goodness-of-fit test for KwaZulu-Natal GEVD model
Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests were used to determine whether maximum monthly rainfall data for KwaZulu-Natal follow a stationary GEVD model. Table 10 presents the K-S and A-D goodness-of-fit tests results for KwaZulu-Natal GEVD.
From Table 10, the p-values for both K-S and A-D tests are insignificant (p > 0.05) at 5% level of significance. Thus, we conclude that the maximum monthly rainfall for KwaZulu-Natal follow the specified stationary GEVD model. The stationary GEVD model for Limpopo data (i.e., model M 0 ) has a maximum NLLH of 669.707. A GEVD model with linear trend in the location parameter (i.e., M 1 ) has a maximum NLLH of 666.705 (see Table 11). The deviance statistic for comparing these two models is therefore D = 2(669.707 − 666.705) = 6.004, which is greater than the critical value of 3.841 with 1 degree of freedom. Therefore, model M 1 provides an improvement in fit over the stationary GEVD model. The likelihood ratio test for µ 1 = 0 has p-value = 0.005, which is significant at 5% level of significance (p < 0.05). This clearly shows that the non-stationary GEVD model is worthwhile and provides an improvement in fit over the stationary GEVD model.
Consider the pair of models (M 0 , M 2 ) from Table 11. The deviance statistic is 2(669.707 − 665.327) = 8.760, which is large compared to χ 2 1 (0.05) = 3.841. Thus, allowing for a linear trend in the scale parameter improves on the stationary GEVD model. The likelihood ratio test for σ 1 = 0 has p-value of 0.001, implying that the linear trend in the scale parameter is significant at 5% level of significance (p < 0.05). This indicates that model M 2 is important and does provide an improvement in fit over the stationary GEVD model. From Table 11, the pair of models (M 0 , M 3 ), has the deviance statistic of 11.014, which is greater than the critical value of 5.991 with 2 degrees of freedom, implying that model M 3 provides an improvement in fit over the stationary GEVD model. The likelihood ratio test for µ 1 = 0 has p-value = 0.067, which indicates that the likelihood ratio test is not significant at 5% level of significance (p > 0.05), while the likelihood ratio test for σ 1 = 0 has p-value = 0.013, which suggests that the likelihood ratio test is significant at 5% level of significance (p < 0.05). The other pairs from Table 11, (M 0 , M 4 ) and (M 0 , M 5 ), have deviance statistic values of 19.040 and 7.900, respectively. These results suggest that model M 4 , which allows for nonlinear quadratic trend in the location parameter and linear trend in the scale parameter, provides an improvement in fit over the stationary GEVD model since the value of the deviance statistic (19.040) is larger as compared to the value of χ 2 3 (0.05) = 7.815. The likelihood ratio test for µ 1 = 0 has p-value = 0.001, for µ 2 = 0 it has p-value of 0.002, and for σ 1 = 0 it has p-value = 0.034, which are all significant at 5% level of significance (p < 0.05). Also, model M 5 which allows for linear trend in the location parameter and a nonlinear quadratic trend in the scale parameter, provides an improvement in the stationary GEVD model since the value of the deviance statistic is greater than the value of χ 2 3 (0.05) = 7.815. The likelihood ratio test for µ 1 = 0 has p-value = 0.236, which is not significant at 5% level of significance (p > 0.05), while the likelihood ratio tests for σ 1 = 0, and σ 2 = 0, all have p-values < 0.001, which are both significant at 5% level of significance (p < 0.05).
The model pair (M 0 , M 6 ), which allows for nonlinear quadratic trend in both the location and scale parameters in Table 11, has a deviance statistic of 9.046 which is small compared to the critical value of 9.488 with 4 degrees of freedom. Thus, allowing for a quadratic trend in both the location and scale parameters is not worthwhile in fit over the stationary GEVD model M 0 . The likelihood ratio test for µ 1 = 0 has p-value = 0.145, and for µ 2 = 0 has p-value = 0.185, which is insignficant at 5% level of significance (p > 0.05), while the likelihood ratio test for σ 1 = 0, and σ 2 = 0, all have p-values < 0.001, which are both significant at 5% level of significance (p < 0.05).
Consider the model pair (M 0 , M 7 ) in Table 11 with deviance statistic of 15.820, which is greater than the critical value of χ 2 2 (0.05) = 5.991, indicating that the non-stationary GEVD model provides an improvement in fit over the stationary GEVD model. The likelihood ratio tests for µ 1 = 0, and µ 2 = 0 have p-values < 0.001, which indicate that the likelihood ratio tests are significant at 5% level of significance (p < 0.05) for the quadratic trend in the location parameter with no variation in the scale parameter. This implies that the non-stationary GEVD model, M 8 , is worthwhile and does give an improvement in fit over the stationary GEVD model.
Consider the model pair (M 0 , M 8 ) from Table 11 with χ 2 2 (0.05) = 5.991 and deviance statistic of 9.338. The likelihood ratio tests for σ 1 = 0 and σ 2 = 0 have p-values <0.001. These results show that the nonlinear quadratic trend in scale parameter with no variation in the location parameter is significant at 5% level of significance (p < 0.05). The deviance statistic (9.338) is greater than the critical value of 5.991, which implies that the nonstationary GEVD model, M 8 , is important and does provide an improvement in fit over the stationary GEVD model.
Overall, Limpopo has five competing non-stationary GEVD models: M 1 , M 2 , M 4 , M 7 and M 8 , for which only two models were considered based on their deviance statistic values as main and alternative best models. The best non-stationary GEVD model is M 4 , which has a nonlinear quadratic trend in the location parameter and a linear trend in the scale parameter, and is given by The alternative non-stationary GEVD model is M 7 , which has a nonlinear quadratic trend in the location parameter and no variation in the scale parameter, and is given by: The shape parameters in (31) and (32), that is, 0.040 and 0.047 for the models M 4 and M 7 , respectively, are positive, which suggests that the rainfall data for Limpopo can be modelled using the Fréchet distribution class since the shape parameter ξ > 0. The diagnostic plots for the non-stationary GEVD model in (31) are presented in Figure  6. The results in Figure 6 show that model M 4 is the best fit for Limpopo maximum monthly rainfall data since the two diagnostic plots indicate a reasonable good fit for the non-stationary GEVD model with a nonlinear quadratic trend in the location parameter and a linear trend in the scale parameter.

Goodness-of-fit test for Limpopo non-stationary GEVD model
Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests were used to determine whether maximum monthly rainfall data for Limpopo follows the non-stationary GEVD model, M 4 . Table 12 presents the K-S and A-D goodness-of-fit tests.
From Table 12, the p-value for the K-S test is insignificant (p > 0.05), implying that the maximum monthly rainfall for Limpopo follows the non-stationary GEVD model, while the results from the A-D test suggest that the maximum monthly rainfall for Limpopo do not follow the specified non-stationary GEVD model (p < 0.05). This contradiction in the results of the two goodness-of-fit tests may be a cause for concern, and may suggest that the selected non-stationary GEVD model, M 4 , may not model the extreme right tails of the Limpopo maximum monthly rainfall data quite well. The model pairs (M 0 , M 1 ) and (M 0 , M 2 ) in Table 13 share the critical value of From Table 13, the pair of models (M 0 , M 3 ) has a deviance statistic of 19.530, which is greater than the critical value of 5.991 with 2 degrees of freedom, implying that model M 3 provides an improvement in fit over the stationary GEVD model. The likelihood ratio tests for µ 1 = 0 and σ 1 = 0 have p-values < 0.001, which indicate that the likelihood ratio tests are significant at 5% level of significance (p < 0.05) for both the location and scale parameters, implying that the non-stationary GEVD model, M 3 , is important and does provide an improvement in fit over the stationary GEVD model.
The other model pairs from Table 13 The likelihood ratio test for µ 1 = 0 has p-value= 0.392, and for µ 2 = 0 it has p-value of 0.096, which are both not significant at 5% level of significance (p > 0.05), but the likelihood ratio test for σ 1 = 0 has p-value < 0.001, which is significant at 5% level of significance (p < 0.05). On the other hand, model M 5 which allows for linear trend in the location parameter and a nonlinear quadratic trend in the scale parameter, provides an improvement in fit over the stationary GEVD model since the value of the deviance statistic is greater than the value of χ 2 3 (0.05) = 7.815. The likelihood ratio test for µ 1 = 0, σ 1 = 0 and σ 2 = 0, all have p-values < 0.001, which are significant at 5% level of significance (p < 0.05).
The model pair (M 0 , M 6 ) in Table 13, which allows for nonlinear quadratic trend in both the location and scale parameters, has a deviance statistic of 24.512 which is greater than the critical value of 9.488 with 4 degrees of freedom. Thus, allowing for a quadratic trend in both location and the scale parameters is worthwhile in fit over the stationary GEVD model, M 0 . The likelihood ratio test for µ 1 = 0 has p-value = 0.499, and µ 2 = 0 has p-value = 0.303, which is insignificant at 5% level of significance (p > 0.05), while the likelihood ratio tests for σ 1 = 0 and σ 2 = 0 all have p-values < 0.001, which are significant at 5% level of significance (p < 0.05).
Consider the model pair (M 0 , M 7 ) in Table 13, with a deviance statistic of 6.394 which is greater than the critical value of χ 2 2 (0.05) = 5.991. These results show that the nonstationary GEVD model provides an improvement in fit over the stationary GEVD model. The likelihood ratio test for µ 1 = 0 has p-value = 0.369 and for µ 2 = 0 it has p-value = 0.449, which are both not significant at 5% level of significance (p > 0.05). This implies that model M 7 , with a quadratic trend in the scale parameter and no variation in the location parameter is not worthwhile over the stationary GEVD model.
These results show that the nonlinear quadratic trend in scale parameter with no variation in the location parameter is significant at 5% level of significance (p < 0.05). The deviance statistic (29.150) is greater than the critical value of 5.991, which implies that the nonstationary GEVD model, M 8 , is important and does provides an improvement in fit over the stationary GEVD model. In general, Mpumalanga has five competing non-stationary GEVD models: M 1 , M 2 , M 3 , M 5 and M 8 , for which only two models were considered based on their deviance statistic values as main and alternative best models. The best non-stationary GEVD model is M 8 , which has a nonlinear quadratic trend in the scale parameter and no variation in the location parameter, and is given by The alternative non-stationary GEVD model, is M 5 , which has a linear trend in location parameter and nonlinear quadratic trend in scale parameter and is given by: The shape parameters in (33) and (34), that is, −0.161 and −0.006 for the respective models M 8 and M 5 are negative, which indicate that the rainfall data for Mpumalanga can be modelled using Weibull distribution class since the shape parameter ξ < 0. The diagnostic plots for the non-stationary GEVD model in (33) are presented in Figure 7. The results in Figure 7 show that the non-stationary GEVD model, M 8 , is the best fit for Mpumalanga maximum monthly rainfall data since the two diagnostic plots suggest a reasonable good fit for the non-stationary GEVD model with a quadratic trend in the scale parameter and no variation in other parameters.

Goodness-of-fit test for Mpumalanga non-stationary GEVD model
Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests were used to determine whether maximum monthly rainfall data for Mpumalanga follow the non-stationary GEVD model, M 8 . Table 14 presents the K-S and A-D goodness-of-fit test results for Mpumalanga non-stationary GEVD model, M 8 .
From Table 14, the p-value for the K-S test is insignificant (p > 0.05), implying that the maximum monthly rainfall for Mpumalanga follows the specified non-stationary GEVD model. On the other hand, the results from the A-D test contradict the results from the K-S test. The explanation for this contradiction is similar to that given for the Limpopo province best model.

Conclusions
In this paper, stationarity test, which included the augmented Dickey-Fuller (ADF), Phillips-Perron (PP) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests, was done. The findings from the KPSS test suggest that monthly rainfall data for all the five provinces are not stationary, while the findings from the PP test contradict those from KPSS test. On the other hand, the findings for the ADF stationarity test for Eastern Cape, Limpopo and Mpumalanga suggest that the monthly rainfall data are stationary, which is a further contradiction to the KPSS findings. However, ADF stationarity test findings for Gauteng and KwaZulu-Natal provinces concur with those from KPSS test. The study also employed Jarque-Bera (JB), Shapiro-Wilk (SW) and chi-square test methods to check whether the monthly rainfall data were normally distributed. Findings from the JB, SW and chi-square normality tests revealed that the monthly rainfall data for all the five provinces do not come from a normal distribution.
The study analysed the long-term trends of monthly rainfall data in the five selected provinces of South Africa from 1900 to 2017. Two trend analysis techniques were applied in this study, the Mann-Kendall test and Sen's slope. Findings from the Mann-Kendall test revealed statistically significant monotonic decreasing trends in Eastern Cape, Gauteng and Kwazulu-Natal provinces, while in Limpopo and Mpumalanga provinces the trends were also revealed to be monotonically decreasing, but insignificant. The Mann-Kendall test statistic findings for Eastern Cape, Gauteng, KwaZulu-Natal and Limpopo were in agreement with the findings from Sen's slope estimator method. However, the Mann-Kendall test findings for Mpumalanga slightly differed from Sen's slope estimator findings. This slight difference can be interpreted as a confirmation of the insignificance of long-term trend and slope for Mpumalanga monthly rainfall.
The study further analysed and discussed in detail the modelling of monthly rainfall extremes using the non-stationary GEVD approach which belongs to the block maxima realisation [40]. The maximum likelihood estimation method was used to obtain the estimates of the parameters. The stationary GEVD was found as the best distribution model for Eastern Cape, Gauteng and KwaZulu-Natal provinces. Furthermore, model fitting supported non-stationary GEVD models for Limpopo and Mpumalanga maximum monthly rainfall, with nonlinear quadratic trend in the location parameter and a linear trend in the scale parameter for Limpopo, while for Mpumalanga the non-stationary GEVD model with a nonlinear quadratic trend in the scale parameter and no variation in the location parameter fitted well to the maximum monthly rainfall data. The study further revealed that the maximum monthly rainfall for Eastern Cape and Mpumalanga can be modelled by distributions in the negative-Weibull domain, while maximum monthly rainfall data for Gauteng, KwaZulu-Natal and Limpopo follow distributions in the Fréchet distribution class.
Model diagnostics, which included the Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests among others, further confirmed that the maximum monthly rainfall for Eastern Cape, Gauteng and KwaZulu-Natal follow the stationary GEVD, while for Limpopo and Mpumalanga the K-S findings showed that the maximum monthly rainfall for these two provinces follow the non-stationary GEVD model. The latter findings could not be confirmed by the A-D goodness-of-fit test.
Findings from this study can help us with information necessary for decision makers to establish strategies for proper planning of agriculture, infrastructure, drainage system and other water resource applications in South Africa. These findings may also assist South African government agencies to improve the socio-economic conditions of the country under the changing rainfall patterns and impending global warming.
This study will form a benchmark for monthly rainfall studies of this kind in these provinces of South Africa. Further studies may look to extend this study into spatial extremes, copula and conditional extremes modelling, as well as Bayesian extreme value modelling approaches.
Author Contributions: All authors contributed equally to the production of the article. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors would like to acknowledge South African Weather Service (SAWS) for providing the data that was used for this research study. We are also very thankful to the University of Limpopo where this study was carried out.

Conflicts of Interest:
The authors declare no conflict of interest.