The GEV distribution is discussed by Fisher and Tippett (1928) as the most suitable approximation for modelling the maxima or minima of a long sequence of finite variables. Suppose that the independent and identically distributed (i.i.d.) finite sequence X₁,…,X_n constitutes a random sample of size n that is chosen from the random variable X whose marginal distribution function is F. Let M_n = max{X₁,…,X_n}, then the extremal types theorem states that for the suitable normalising constants {a_n > 0} and {b_n}, the distribution function of the re-scaled block maxima Mn* is given by:

for a non-degenerate distribution function G, then G belongs to a unified GEV distribution family that is given by:

defined for {x:1+γ(x−μσ)>0},

where − ∞ < μ < ∞ is the location parameter, σ > 0 is the scale parameter and − ∞ < γ < ∞ is the shape parameter.

The ordinary GEV distribution in equation 2 converges to one of the three families of extreme value distributions, depending on the rate of decay of the tail that is indexed by the shape parameter γ. When γ = 0, G_γ(x) reduces to a type I or a short-tailed unbounded Gumbel family of distributions (Beirlant et al. 2004). When γ < 0, G_γ(x) is thin-tailed and we get a type II family which is the Weibull class of distributions with an upper bound given by (Beirlant et al. 2004):

If γ > 0, then G_γ(x) belongs to a type III family which is a heavy-tailed Frèchet class of distributions that is bounded below by (Beirlant et al. 2004):

The survival distribution of the GEV distribution is given by:

defined for {x:1+γ(x−μσ)>0} and y = 0.

By letting p = Pr(X > x) and rearranging equation 3, we get the quantile function that is given by:

As p → 0 and γ<0,xp=μ−σγ.

The quantile function given in equation 4 is used in estimating high quantiles and predicting the probability of exceedance levels in section ‘Empirical results and discussion’. The GEV distribution in equation 2 is then extended to give an asymptotic model for the joint distribution of the r largest order statistics within annual blocks, for fixed values of r. We initially define:

and then identify the limiting behaviour of this variable, for fixed k, as n → ∞. If equation 1 holds, then, for fixed k:

defined on {x:1+γ(x−μσ)>0}.

where:

with:

The GEV distribution for the r largest order statistics:

is the joint probability density function given in Coles (2001) as:

valid for − ∞ < μ < ∞, σ > 0 and − ∞ < γ < ∞; x^(r) ≤ x^(r–1) ≤ …≤ x⁽¹⁾; and:

x(k):1+γ(x(k)−μσ)>0 fir k = 1,2,…,r.

The likelihood function for the r largest order statistics model when γ = 0 is given by:

defined for 1+γ(xi(k)−μσ)>0,k=1,…,ri,i=1,…m.

However, the validity of the GEV distribution for r largest order statistics depends on the choice of the number of order statistics r, which is critical. If r is too large, the asymptotic support for the model is likely to be violated, leading to the occurrence of bias (Bader et al. 2015; Coles 2001). Similarly, if r is too small, few observations will be generated and the variance of the estimator can be high (Bader et al. 2015; Coles 2001). The details of the asymptotic theory of order statistics that are highlighted in this study are obtained in various parts of the literature, including Falk and Wisheckel (2016).

A plot of the AMDT in degrees Celsius using data for the non-winter season (September–April of each year) is given in Figure 1 for the period 2000–2010. The results of the r largest asymptotic order statistics model that is fitted to the 10 largest annual temperatures over a period of 11 years are summarised in Table 1. Several values of the maximised negative log-likelihood estimate (λ_i), for i = 1,…,10, have led to the maximum likelihood estimates of the target parameters as illustrated in Table 1 with standard errors in parentheses. However, the attention is limited to r ≤ 6 order statistics as a result of the reasonable doubt on the validity of the model for r ≥ 7. The maximised negative log-likelihood estimates (λ_i) seem to be stable for r ≤ 6 and start decreasing significantly when r ≥ 7. The standard errors of estimates increase significantly when r ≥ 7, implying poor fit of the model when r ≥ 7. According to Soares and Scotto (2004), the above-mentioned observations also support the doubt of the fitted GEV_r distribution when r ≥ 7. For r ≤ 6, the standard errors of the estimates (μ^,σ^, and γ^) decrease as the values of r increase, implying an increase in precision of the model (Bader et al. 2015; Coles 2001; Soares & Scotto 2004).

The sufficient evidence of the validity of the Weibull distribution family towards modelling of AMDT in South Africa is revealed by the estimates of the shape parameter for all values of r. A negative value of the shape parameter in classical EVT implies that the GEV distribution converges towards the Weibull class of distributions. All estimates of γ^ in Table 1 are negative values and the Weibull class of distributions is the most appropriate. The graphical diagnostic tools (probability–probability [P–P] plot and quantile–quantile [Q–Q] plot) for assessing accuracy in the fit of the annual maximum temperature to the r largest asymptotic GEV distribution are given for each value of r in Figures 2, Table 1-A1 and Table 2-A1. However, the paramount attention in checking the fit of the models in this study is limited to the Q–Q plots based on the fact that the linear pattern of the Q–Q plot is unlikely to be affected by the shifts in the location, scale and symmetry of the distribution, an advantage which is not the case in the P–P plot. Looking at Figures 2 and Table 1-A1, the Q–Q plots for r = 2, 3 and 4 display the dots that are near linear. There is doubt on the fit of the model at r ≥ 6. The model fitted for r = 4 with λ₄ = − 8.5826 seems to be the one that possesses a reasonably good fit. The graphical diagnostics which are the P–P, Q–Q, return level and the density plots are given for r = 4 in Figure 3-A1. To justify the appropriateness of the Weibull class of distributions as a proper model for the AMDT in South Africa, the confidence interval for γ is estimated for all values of r. All the point estimates of the shape parameter γ are found to be negative values. It is important to compare the outcomes of point estimation with the results of interval estimation. If indeed the point estimates of the shape parameter are negative values, then they must be bounded below and above by the negative confidence limits. For example, the confidence interval for γ considering r = 1 is given by:

This leads us to conclude that at a 95% level of confidence, the value of the shape parameter γ is expected to be enclosed within the interval (−1.1505; −0.2191), which makes sense because γ^=−0.6848 when r = 1. The confidence interval for γ is then estimated for all r values, and it is noted that all the upper limits are negative values, and hence enclose the point estimates γ^. This justifies that the Weibull class of distributions can be used for modelling AMDT in South Africa. These results are included in the last column of Table 1.

There are various tests that are considered in the literature as relevant to assess the goodness of fit of the models to the data. The Anderson–Darling test is suitable for assessing the goodness of fit of the heavy-tailed distributions. In this article, we use the EDT that is discussed in Bader et al. (2015) as a specification test for assessing the goodness of fit of GEV_r for r largest order statistics (r > 1). This is derived as the difference in entropy of GEV_r and GEV_r-1 models. Considering a continuous random variable X with density function f, the entropy difference is given by:

which is the mean of the negative log-likelihood that is estimable by the sample average of the contribution to the log-likelihood from the observed data. Suppose that r − 1 top order statistics fit the GEV_r-1 distribution. Then the difference in the log-likelihood between GEV_r-1 and GEV_r provides a measure of deviation from the null hypothesis H0r. However, large deviation from the expected difference under H0r suggests a possible misspecification of the null hypothesis (Bader et al. 2015). From the log-likelihood contribution, the difference in log-likelihood for the i th block:

is given by:

It is critical to test the goodness of fit of the GEV_r distribution with a sequence of null hypotheses (Bader et al. 2015). This study desires to test the null hypothesis H0r: the GEV_r distribution fits the sample of the r largest order statistics well for r = 1,…,R, where R is the maximum, predetermined number of top order statistics to test. The extreme value analysis (EVA) R package of Bader and Yan (2016) is used for conducting the EDT in this article.

The results of the test are summarised in Table 2, which shows in column 1 the values of order statistics r to be tested. The p-values from the individual tests at each value of r in column 2 establish r = 4 as the most appropriate GEV_r model to the data. This is based on p = 0.7834 which is the highest probability value. Basically, the rule is to reject H0r if p-value ≤ α, where α is the level of significance. At 10% level of significance, we fail to reject H0r when r = 4, implying that GEV_r=4 is the most appropriate model for AMDT in South Africa. Looking at the test statistics in column 4, the critical value Z_α/2 = 1.645 also supports that we fail to reject H0r when r = 4. The p-values of the Forward and the Strong Stops are also provided together with parameter estimates for each value of r in Table 2.

Estimates of extreme quantiles and probabilities of exceedance levels of the annual maximum distribution are obtained using the quantile function given in equation 4. These are calculated with the use of parameter estimates of model r₄ which is established as the best fitting model to the data. The quantity x_p in equation 4 is the return level associated with the return period 1/p. The level x_p is expected to be exceeded on average once every 1/p years (Coles 2001). The deviance statistic is used in assessing the goodness of fit of the models. This is given by:

where λ(r_i) and λ(r_j) are the maximum likelihood functions of r_i and r_j, respectively. A test of the validity of model based on r_i relative to r_j at the α level of significance is to reject r_i if D_(i,j) > Cα, where Cα is the (1 – α) quantile of the χ12 distribution (Coles 2001; Smith 1989; Soares & Scotto 2004).

The point estimates x^p for the 5, 10, 25, 50 and 100 years return levels corresponding to the return period 1/p are calculated using equation 4 and the outcomes are summarised in Table 3. Looking at 25 years return level, for example, the level x_0.04 = 32.15 is expected to be exceeded on average once every 25 years. For the greater accuracy in the estimation of the return level and parameters in the model for r = 4, profile log-likelihoods for µ, σ and γ are given in Figures 4-A1, 5-A1 and 6-A1, respectively. Table 4 shows the comparison of the models using the deviance statistic given in equation 15. For example, the deviance statistic for comparing λ(r₁) and λ(r₂) is calculated as D_(1,2) = 2 (−15.9670 – (−12.5255)) = −6.88. Similar to Soares and Scotto (2004), the deviance statistics in this article are significant when comparing the log-likelihoods among D_(2,3), D_(3,4) and D_(4,5). The comparison is invalid for D_(1,2) and D_(5,6). The test for selecting the best model is conducted at 1% level of significance for which χ12=6.64. For D_(2,3) and D_(4,5), λ(r₃) and λ(r₅) are rejected as D_(2,3) > 6.64 and D_(4,5) > 6.64. It is therefore clear that D_(3,4) = 5.49 < 6.64 reveals failure to reject λ(r₄), which implies the validity of r = 4 as an order statistic that possesses the reasonably good fit as suggested by the graphical diagnostic tools. The Q–Q and the P–P plots given in Figure 1-A1 show that the Weibull family is the appropriate distribution for the maximum temperatures and also that r = 4 possesses the most reasonable fit of the model out of the six order statistics. Equation 4 is also used for estimating the future return levels (extreme quantiles) for the different return periods as illustrated in Table 5. For example, the 90th quantile is calculated as follows:

The number of observations that are above the estimated tail quantile x_0.1 = 32.0 are then counted and found to be 2. For the observed number of exceedances, we get 0.1 × 44 = 4.4 ≈ 4, where 44 is the number of observations available in 11 years by four order statistics.

Table 5 presents a summary of the estimated tail quantiles at different tail probabilities. The tail quantiles (temperatures) are given in column 2. The observed number of temperatures that is larger than the estimated tail quantiles is shown in column 3, while column 4 shows the corresponding number estimated using the GEV distribution. The estimation of the tail quantiles in this article establishes maximum temperatures that are expected to be exceeded. Considering the 90th quantile, for example, x_p = 32 °C, is the maximum temperature that is expected to be exceeded on average once every 10 years. This is the temperature above which we expect the occurrence of extreme heat that may lead to natural hazards such as hot spells or heat waves. From the viewpoint of energy demand forecasting, this is the maximum temperature above which the demand of electricity does not increase significantly.