Seasonal and Station Effects Modelling to Extreme Temperature Data in South Africa

Extreme value theory has been widely applied to weather variables, and rigorous approaches have also been employed to investigate the seasonality and dependencies to extreme values of weather variables. To investigate the seasonal and station effects of daily maximum and minimum temperatures data, station and season specific effects model have been introduced in the parameters of general Pareto distribution. Then, the seasonality and station variations that are inherent in the data under consideration were assessed applying mainly the Bayesian approach. Non-informative and informative priors were used for estimation of the parameters. The seasonal and station effects parameters of the general Pareto distribution were estimated through the introduced models, allowing the sharing of information between stations and seasons. Simulation study was also carried out to investigate the precision of estimators for the GPD parameters with and without the effects, station and seasonal, to simulated data. The models employed improved precision of the station and seasonal effects parameter estimators at individual stations and in individual seasons. The study also depicted the significance of introducing seasonal and station variabilities when modelling extreme values using univariate method, which allows information to be pooled across stations and seasons. Results obtained in this study have essential scientific and practical applications. In an extreme temperature setting, designing a level without taking the station and seasonal effects into account could lead to significant under-protection. Hence, it is important to consider what is expected to be colder or warmer than usual by identifying the effects of stations and seasons in the analysis. This would benefit greatly local governments, researchers and farmers, which they can use to suggest adaptation and mitigation steps to improve resilience.


Background
A significant body of literatures have demonstrated that air and water temperatures increases globally over the past decades [1,2]. Various studies have also showed similar trends in South Africa [3][4][5][6][7][8]. Much evidence exists for an increase in the variability of rainfall, as well as the number and size of extreme wet and dry spells in several parts of the southern African region [5,[9][10][11][12]. Climate change and global warming with its associated extreme temperatures characterised by heat and cold waves pose serious economic and health challenges. Increasing temperatures and changing weather conditions affect the heating and cooling systems consequently leading to increased electricity demand for air conditioning in summer and heating systems in winter [8].
Historical data have shown that temperature has been increasing over the past decades. According to Sunday times [13], Vioolsdrif, a village in the Northern Cape of South Africa, has recorded the highest daily maximum temperature of 46 • C in 1995 and then it broke its own new record again for the highest temperature in the country reaching 53.2 • C in late 2019. Temperature extremes, which have been known to be caused by an increase in the concentration of greenhouse gases, are natural phenomena that affect our socio-economic activities. Extremely high temperatures also cause drought in the agricultural sector resulting in economic hardships and loss of human lives and livestock [14].

3
The weather conditions like extreme temperature have direct consequences on agricultural and energy generation sectors in South Africa. Moreover, it is highly expected that global temperatures will continue to rise for decades to come [15]. The rise could be largely due to greenhouse gases produced by human activities and the temperature rise has not been, and will not be, uniform or smooth across the country or over time [16]. For this reason, a realistic assessment of future behaviour of extreme temperatures is fundamental in order to understand the challenges ahead. Most importantly, building future situations will provide important input for farmers, researchers and local governments, who can use this information to propose adaptation and mitigation measures to increase resilience. In other words, appropriate policies and plans can be drawn to prepare the stake holders for changes due to extreme temperatures.

Related Literature Reviews
Numerous studies have been conducted to investigate the non-stationarity and temporal dependence of extreme temperature data in South Africa and other parts of the world e.g. [17][18][19]. Trends in extreme temperature over Nigeria from percentile-based threshold indices were analysed by [18]. The study focuses on spatial and temporal trends in extreme temperature for selected 21 weather stations in Nigeria over the years 1971-2012. It was concluded in [18] that although a majority of the stations have significant trends on warm days and warm nights, the annual trend is greatest on warm nights. The study done in Malaysia, using extreme value distribution to maximum temperature data of Penang weather station, suggested that the best model was one with a location parameter that increased linearly with time [17].
Extreme temperature was considered by [20] to see the influence on daily electricity demand in South Africa. The study established temperature as an important variable in explaining electricity demand; however, the impact varied over seasons. The empirical results showed that electricity demand in South Africa is highly sensitive to cold temperatures. Nemukula and Sigauke [19] investigated the general Pareto distribution (GPD) modelling of average minimum daily temperature in South Africa and obtained that temperature data exhibit evidence of short-range dependence and high seasonality. More recently, the point process modelling approach is considered by [21] to model both the frequency and intensity rates of the occurrence of extreme temperature data. The results of the study based on the data for the years 2000-2010 established temperature tends to be very high on consecutive days, but very hot spells tend not to last longer than 1 or 2 days.
Although temperature is an important weather variable in the energy and agricultural sector, it is important to study the impact of other climate variables such as rainfall and wind speed among others. For instance, Mason et al. [9] found evidence of significant increases in the annual extreme rainfall between 1931 and 1960 and between 1961 and 1990. Focusing on the eastern part of South Africa, Groisman et al. [12] found no significant change in the annual and summer rainfall totals for the region during the period 1906 to 1997, but they found a significant increase in the annual frequency of very heavy rainfall. Du Plessis and Burger [22] explored the magnitude and frequency of short duration rainfall events in Western Cape using extreme value distributions to nonstationary sequences for both event maxima and peak over threshold (POT) modelling. They concluded that there was no evidence of trends or indications of changes in rainfall intensities. Recently, De Waal et al. [23] applied GPD and POT methods and their results indicated that the general assumption of stationarity in design rainfall assessments should be strongly questioned. A study closely related to the present study in EVT is that of [24] who used hierarchical model to hourly extreme wind speed data using GPD with non-informative prior. Their findings, which regard station and seasonal effects as random, showed that pooling of information improved the precision of GPD model parameters fitted to the data.
On the other hand, various authors have reported that extreme value distributions have a wide range of applications to deal with extremal behaviour of weather data (see [20,[25][26][27][28][29][30]). For example, Stein [30] obtained negative shape parameter estimates for the GEV model fitted to annual temperature data at Belvedere Tower station in Central Park, New York City. Similarly, Huang et al. [29] also obtained negative shape parameter estimates for annual temperature data at the USA and some surrounding areas. However, negative shape parameter estimates for extreme value distributions may lead to underestimation of the most extreme quantiles [30]. Hence, the implications of negative shape parameter estimate for extreme temperatures should be interpreted carefully [29,30].

Motivation
Various literatures used a pragmatic way of bypassing the problem of variations and dependence within the extreme datasets of environmental variables [8,21]. More recent application to extreme temperature in South Africa is conducted by [8] to model dependences of high extreme temperature and they found strong positive extremal dependences between weather stations considered. Strong extreme dependencies were obtained, but the performance gains between the inference approaches are conflated with possible differences due to the inference paradigm used. The occurrences of extreme observations lead to the use of extreme value theory (EVT) which comprises a rich family of techniques that are suitable for modelling tail behaviour using POT approach that is modelled using the GPD [31]. Extreme data time series may indicate that the observations show seasonal patterns and variations. Dealing with seasonality is crucial to avoid the limiting argument for an extreme value analysis. This is because it can affect dependence on extreme events. The principal motivation of this paper is therefore to investigate and compare the GPD parameters and return levels which are estimated with and without seasonal and station effects, and to show how the effects are applied to model dependences of extreme temperatures data in South Africa.
The interest in inferences combining the classical extreme value modelling with the Bayesian inference is another aspect of the current study that is worth emphasising. The posterior standard deviations for the parameters provided by a Bayesian framework are much easier to interpret than the standard deviations in a frequentist approach [32]. Extreme temperatures data were analysed applying the Bayesian method by various authors (see, e.g., [20,28,33]) to evaluate the posterior distributions of the quantiles and the parameters of extreme value distributions. These authors showed that a Bayesian approach provided a complete representation over the frequentist approach, giving a full description of uncertainty in parameters and quantiles. However, they have not considered the seasonal and station effects on the prior formulation for their approaches linking with the performances of the estimators. This motivates the authors to derive a full description of uncertainty in parameters and quantiles of EVDs.

Contribution
From practical point of view, the authors attempt to quantify observable changes in extreme temperature for seasonally clustered months across selected weather stations considered in South Africa. Although various studies have been done across South Africa as a whole [6][7][8]21], no similar studies have focused on the estimation of seasonal and station effects on GPD parameters, allowing the weather stations to vary seasonally, in order to determine similarities and differences in changes within and across stations from the Bayesian perspective. A hierarchical model to hourly extreme wind speed data was used to model GPD by [24] using noninformative prior only. The present study will adopt the GPD modelling approach of [24] to cater for seasonal and station effects in the minimum and maximum temperatures data. Also, the authors extend their approach to create clustering based on seasonal dependence to the GPD using both non-informative and informative priors where the later were formulated based on seasonally clustered historical data of surrounding weather stations. Therefore, the current study focuses on introducing seasonal and station effects to an extreme value distribution using the POT approach to understand the behaviour pattern of extreme weather events. This quantifies changes due to extreme temperature events and helps to detect the changes that have already occurred for better preparedness.

Paper Organisation
The rest of the paper is organised as follows. In Sect. 2, the materials and methods for the analyses are introduced. Some exploratory data analyses were also done in this section. Models for threshold exceedances are explained briefly in Sects. 2.2 and 2.3. In Sect. 2.4, models for seasonal and station effects parameters for the GPD are briefly explained. The Bayesian estimation approach, which was employed using informative and non-informative priors, is explained in Sect. 2.5. The results and discussions for the seasonal and station effects modelling to daily maximum and minimum temperatures data at selected weather stations across South Africa are reported in Sect. 3. In Sect. 3.3, simulated data was used to study the properties of estimators for the GPD parameters obtained with seasonal and station effects, in comparison without the effects. Concluding remarks and recommendations for future studies are given in Sect. 4.

Methodology
For this study, the temperature data that correspond to a 65-year spanning the period from 1949 to 2019 series of daily maximum and minimum temperatures measured from selected weather stations in South Africa were employed. The spatial plots of network of weather stations with their respective average minimum and maximum temperatures value are given in Fig. 1. Across South Africa, a wide range of activities are influenced by differences in seasonality and there is also little consensus on the timing of seasonal boundaries [34]. For instance, the boxplots in Fig. 2 indicate that the observations of extreme daily maximum and minimum temperatures data of the weather stations may exhibit seasonal patterns and variations. It is therefore important to deal with seasonality, as this may affect the dependence of extreme events. The temperature data exhibit strong seasonality, which changes periodically through time, and appropriate measures should be taken to account for this seasonal variation in the analyses of extreme weather datasets. The most widely adopted technique to deal with data that vary seasonally is to partition the data into seasons, within which the data are assumed to be homogeneous. These seasons might for example be 'winter' and 'summer', or 'dry' and 'wet', where the seasonal variation is clearly understood 1 3 [34][35][36]. Therefore, in this paper the daily temperature data are statistically classified to create seasonal divisions of the daily maximum and minimum temperatures data.
A Euclidean cluster analysis was performed applying the Ward's D method using cluster packages in R and the seasonal classifications in months for the weather stations     Table 1. For the seasonal classifications, Euclidean cluster analysis was initially supervised at four seasonal divides. The seasonal divides are then validated by using average silhouette width (ASW) computation. The ASW value gauges the strength of within-group homogeneity and the degree of confidence in between-group distances. If the cluster was not significant, two, three, five, and six seasonal divides were utilised sequentially until it was significant [34]. Two seasonal periods were selected for analysis, winter months and summer months for the minimum temperature and maximum temperature data, respectively. Statistically classified seasonal periods are summer months that include November, December, January, February, March and April while winter months that include June, July, August and September for most of the weather stations. However, the seasonal classifications for the majority of the stations include from December to March for the maximum temperature and from June to August for the minimum temperature data.

Methods for Extreme Events
To analyse the extreme values of temperature data statistically the POT method was used, where daily values above a pre-determined threshold value were modelled by GPD [37]. The GPD and estimation methods of its parameters are presented next.

The Peak over Threshold Approach
Commonly used limit distributions for modelling extremal events can be summarised using the GPD to model high threshold exceedances [37,38]. Suppose y 1 , y 2 , … is a sequence of independent and identically distributed (IID) random variables with a common distribution F(.). Let y denote an arbitrary term of the sequence. Then, for sufficiently high threshold u, the distribution function of y condition on u, i.e. F(y|u) = P(y ≤ y − u|y > u) , follows approximately a GPD whose distribution function is of the form where (y − u)∕ > 0 provided > 0 and are the scale and shape parameters, respectively, and y ∈ [0, ∞) for ≤ 0 while y ∈ [0, ∕ ) for > 0 [39,40]. The number of observations that exceed the threshold y − u is referred to as exceedances. The shape parameter determines the nature of the tail of GPD. The GPD in expression (1) for = 0 is an exponential distribution with parameter 1∕ (for limit → 0 ). Heavy tailed Pareto distribution is attained when > 0 , while the Beta distribution with upper bound is attained when < 0 . The asymptotic results for the GPD depend on the choice of threshold value. The threshold should be high enough so that the assumptions for GPD are valid, but it also needs to be low enough so that there are sufficient data to make a good fit [41].

Estimation of the Parameters and the Quantiles
The scale and shape parameter of the GPD in expression (1) can be estimated from a time series of the exceedances. Several estimators are available; the most popular one is the maximum likelihood estimator (MLE) [42], which can be straightforwardly extended to the non-stationary case and thus is the method of choice for parameter estimation in this paper. The product of the probability density functions g(y; , ) can be obtained by differentiating the expression (1) with respect to y, and given as The joint likelihood, is therefore given by As a function of the unknown parameters and with given observations y i , the above joint likelihood is not a probability density function itself, but proportional to the probability of observing a set of exceedances y i , i = 1, ..., m given the parameters and . The MLE for the unknown scale and shape parameters can now be obtained by maximising the likelihood L( , ; y i , ..., y m ): In general, expression (4) can be solved by numerical optimisation, e.g. using a quasi-Newton method, to estimate the parameters. The standard errors of MLEs can be approximated asymptotically using the inverse of information matrix [43].
Once the parameter of the GPDs is obtained, the interest is on the application of the fitted model to estimate other quantities. For example, prediction of future extreme temperature is important to make a proper plan to minimise its negative impact. The return level for GPD associated with the return period of 1/p, denoted by y p where 1 − G(y; , ) = 1∕p and 0 < p < 1 , is obtained by inverting expression (1) [37] and is given by The return level d p is therefore a quantile of the GPD corresponding to the upper tail probability p.

Station and Seasonal Effects for the Parameters of Generalised Pareto Distribution
To allow the GPD scale, and shape, parameters vary across stations and seasons, station and seasonal effects were introduced in expression (1). Here, also the assumption was made that the extremes between stations and between seasons are independent. Let the station and seasonal effects for the GPD scale and shape parameters be denoted by j,i and j,i , respectively, for j = 1, 2, ..., 8 and i = 1, 2, ..., 6 , where j and i represent stations and seasons, respectively. The seasonality was also partitioned into different cycles between January and December from 4 to 6 months for daily maximum temperature data while from 3 to 4 months for daily minimum temperature data. Based on these assumptions, the following station and seasonal effects models were specified for GPD parameters: for j = 1, 2, ..., 8; i = 1, 2..., 4, ..., 6 where j and i represent the j th station and i th seasonal effects for scale parameter, respectively, and j and i represent the j th station and i th seasonal effects for shape parameter, respectively. Furthermore, for the combined seasonal and station effects scale parameter, log( j,i ) was used to retain the positivity of the scale parameters and for computational convenience.

Threshold Selection
For the GPD, the first step involved selection of appropriate threshold value u j,i , for j = 1, 2, ..., 8; i = 1, 2..., 4, ..., 6 and then fitting the GPD to the exceedances over u j,i values. It was also necessary to allow the threshold u j,i to vary for each combination of season and station. Therefore, it was important to choose appropriate thresholds for each season and station. For daily maximum and minimum temperatures data, the threshold values, u j,i , for each station j and season i, were chosen using the mean excess plot approach. The adequacy of selected threshold was also checked using plots of MLEs of the shape and modified scale parameters versus threshold values.

Bayesian Inference for Station and Seasonal Effects Parameters
As in the likelihood approach, suppose y 1 , y 2 , … , y m are IID and their distribution fall within the GPD family. However, the parameters j , i , j , i are now treated as random variables for which we specify prior distributions. The specification of priors enables us to supplement the information provided by the data. Let = ( j , i , j , i ) and suppose f Θ ( ) denotes density of a prior distribution for with no reference to the actual data. Then using the Bayes' theorem the density of a posterior distribution for has the form where L( | ) is the likelihood function of GPD given in expression (3) and the integral is taken over the parameter space Θ . Both non-informative and informative priors were used in this paper. The non-informative priors were constructed by assuming there was no additional information available about the parameters, apart from the data, whereas the informative priors were elicited in terms of extreme quantiles. The technique employed is briefly discussed in Sect. 2.5.1.

Prior Formulation for the Station and Seasonal Effects GPD Parameters
The informative priors were formulated according to the approach used in [44][45][46][47], that is, eliciting prior information in terms of extreme quantiles. Therefore, because of the natural ordering of the quantiles d p l , l = 1, 2 in expression (5), i.e. d p 1 < d p 2 , assumption of independent priors on d p l , l = 1, 2 would not be valid, therefore the authors used the quantile differences d In the model proposed, only two quantiles are needed to specify the GPD parameters j , i j , j and i . Therefore, the quantile differences are assumed to be independent Gamma distributions with parameters ( l , l ), l = 1, 2 , that is Based on the Gamma distribution we can also assume Based on expression (9) and the quantile differences d p l , the joint prior for ( d p 1 , d p 2 ), assuming d p 0 = 0 , has the form provided that d p 1 < d p 2 . Then using expression (5) in expression (10) and multiplying by absolute value of the Jacobian of the transformation from ( d p 1 , d p 2 ) to = ( j , i j , j , i ) leads to an expression for the prior in terms of the GPD station and seasonal effects parameters, , presented in Sect. 2.5.2.

Formulation of Informative Priors for the Station and Seasonal Effects Parameters
For the GPD station and seasonal effects, the informative priors were formulated using historical records of daily maximum and minimum temperatures data, using each season and station, of surrounding weather stations. The prior distributions were elicited using the method discussed in Sect. 2.5.1 with quantiles p l = 10 −l , for l = 1, 2 . The return level takes the form d p lj and d p li for station and seasonal effects, respectively, and the Jacobian ( ) is obtained for each station and seasonal effects parameters as and yield the following expressions and respectively, where j and i are indices of station and seasonal effects, respectively. The prior for the GPD station and seasonal effects parameters, , has the form where d p ks is the quantile differences for GPD station and seasonal effects model. Similarly, the posterior density f ( | ) is obtained by substituting the likelihood (where the GPD station and seasonal effects parameters in expression (6) are replaced in to the GPD likelihood in expression (3)) and prior density given in expression (11) in to expression (7). Then the MCMC technique employed with Metropolis-Hastings (MH) algorithm to estimate the features of the posterior distribution f ( | ) of . The joint prior for the informative GPD station and seasonal effects parameters has the form where log( j,i ) = j + i and f Θ ( GPD ) is obtained using expression (11).

Non-informative Priors for the Station and Seasonal Effects Parameters
The non-informative priors for GPD station and seasonal effects parameters, were constructed by assuming there is no external information available to formulate prior distributions apart from the data. Therefore, flat marginal independent priors were employed. Similarly to the informative prior, the joint non-informative prior density for the GPD station and seasonal effects model was then assumed to be where log( j,i ) = j + i .

Posterior Distribution
The extreme value model for seasonal and station effects extremal behaviour is specified by a set of GPD parameters for each component of the parameters. Therefore, the model for the seasonal and station effects would be specified by the vector ( j , i , j , i ) in which each parameter with in the vector determines either the seasonal or the station effects that are inherent in the extreme environmental dataset. From Bayes' theorem in expression (7), the posterior density of the seasonal and station effects GPD model can be formed by the product of the prior f Θ ( GPD ) in expressions (12) and (13), and the likelihood function in expression (3). That is where g(y 1 , ..., y m ) is the density in expression (2) for the seasonal and station effects GPD model and m is the number of observations. The prior is most naturally specified as a product of the joint priors for ( j , i , j , i ) . Then the joint prior for is expressed as in the form in expressions (12) and (13) for the informative and non-informative priors, respectively.

Results and Discussion
For the POT application on daily maximum and minimum temperatures data, a threshold has been chosen using the mean excess plot approach. The adequacy of the chosen threshold can be checked using parameter stability plots, i.e. plots of the MLEs for the seasonal and station effects shape and modified scale parameters ( j , i , j , i ) against a number of different thresholds. If the GPD is a reasonable model for the exceedances of a threshold u, then estimates of shape and modified scale parameters ( j , i , j , i ) should be approximately constant to all threshold greater than u. The MLEs of the shape and modified scale parameters versus threshold values indicate that the selected thresholds are adequate, because the estimates of the seasonal and station effects shape and scale parameters are approximately constant for all thresholds greater than the threshold value for daily maximum and minimum temperatures data. For the maximum temperature data, the exceedances above the threshold value whereas for the minimum temperature data values below the threshold value were considered as extreme observations to fit the GPD.
The prior specifications for the complete parameters vector and estimation methods for the GPD using both informative and non-informative priors are given in Sects. 2.5.2 and 2.5.3, respectively. For the informative prior, priors were elicited based on seasonally clustered data of nearby weather stations, and the precision of the marginal parameters investigated was compared with the non-informative prior.

Estimation of Model Parameters for the Stations and Seasonally Varied Temperatures Data
Here the GPD was fitted to deal with inferences about the station and seasonal effects on daily minimum and maximum temperatures data for selected weather stations in South Africa. First, the dependence was investigated visually using various plots to clarify the data series and introduce appropriately to the extreme value distributions. The seasonally clustered daily maximum and minimum temperatures data are then investigated through the Bayesian analysis using the non-informative and informative priors. For the non-informative prior, the following prior distributions are employed for the GPD station and seasonal effects parameters, namely The Gaussian distributions having zero mean and larger variances in expression (14) enforce the flat independent marginal priors, also known as diffuse or vague priors, which exhibit the lack of external information [48,49].
The MCMC method with MH algorithm was applied for sampling from approximate posterior distributions for the station and seasonal effects parameters. The values generated by the iterations of the chains were checked for convergence for all station and seasonal effects parameters of daily maximum and minimum temperatures data, and convergence was achieved. Similarly, MCMC yielded samples from the approximate posterior distributions for the 8 station effects parameters and 6 1 seasonal effects parameters, that is, for each of j , i , j and i parameters. For simulation, different starting values were considered to distinguish the chain convergence. Hence, all the chains had mixed well with the original data employed in the analysis. The posterior means associated with the standard deviations (in brackets) for the (14) f ( j ) ∼ N(0, 10000); f ( i ) ∼ N(0, 10000); N(0, 100).
GPD station and seasonal effects parameters using noninformative priors are presented in Tables 2 and 3, respectively, for the minimum and maximum temperatures data. The GPD shape parameter estimates that are considered using the sampling variations across seasons and stations were less than zero except for daily maximum temperature data of George and Plettenberg stations effects. The negative result for the shape parameter implies that the model, GPD station and seasonal effects, fitted to daily maximum temperature data had bounded upper tails. Also for the minimum temperature data, the shape parameter estimates, from both seasonal and station effects, suggest bounded upper tails to the distributions of minimum temperature data.
On the other hand, the informative priors were formulated based on seasonally clustered data of nearby daily minimum and maximum temperatures data. The posterior means and standard deviations for the GPD station and seasonal effects parameters formulated from these priors are provided in Tables 2 and 3, respectively. The results show that the posterior means of daily maximum temperature data for the GPD station and seasonal effects parameters obtained using the informative priors are close to the posterior means of non-informative priors. The posterior standard deviations of the station and seasonal effects parameters are lower for the informative priors compared with that of the noninformative priors. These show the benefit of using informative prior for the station and seasonal effects parameters. The station and seasonal effects parameter estimates with their standard deviations for the minimum temperature data are also given in Tables 2 and 3, respectively. In the same way, the posterior means for the GPD station and seasonal effects parameters from the informative priors are close to the posterior means of non-informative priors. All the posterior standard deviations of the station and seasonal effects parameters for the informative priors are lower compared with the non-informative priors.
To investigate how GPD seasonal and station effects parameters were affected by the informative priors formulated based on seasonally clustered historical data of nearby weather stations, the posterior densities of the seasonal and station effects parameters found through the elicited priors and the flat priors were compared. The estimated posterior densities of GPD seasonal and station effects parameters ( j , i , j , i ) for stations and seasonally clustered daily maximum temperature data based on George station during summer period in March are plotted and displayed in Fig. 3. Notice that the distributions of seasonal and station effects parameters obtained based on the non-informative and informative priors are symmetric for this specific station and season given in Fig. 3. On the other hand, the posterior densities of GPD seasonal and station effects parameters ( j , i , j , i ) for the informative priors that were built based on information from seasonally clustered data of nearby weather stations resulted in high peaks compared with that of the non-informative priors. The posterior densities of the remaining seasonal and station effects parameters of the model also depicted high peaks for the informative priors compared with that of the non-informative priors.
Similarly, the estimated posterior densities of GPD seasonal and station effects parameters ( j , i , j , i ) for stations and seasonally clustered daily minimum temperature data based on Durban station during winter period in June are plotted and displayed in Fig. 4. The distributions of seasonal and station effects parameters obtained based on the non-informative and informative priors are symmetric. Furthermore, the posterior densities of GPD seasonal and station effects parameters ( j , i , j , i ) for the informative priors produced high peaks compared with that of the noninformative priors. The fluctuations of the peaks for the posterior densities indicate that the posterior distributions are sensitive to the informative prior from which the prior information was shaped. In the current paper it is crucial to note that the informative priors were constructed through borrowing prior information from seasonally clustered data of nearby weather stations by considering the combinations of their means and quantiles of the data. The weather characteristics of the minimum and maximum temperatures data of those stations used for model fitting as well as the nearby weather stations from which the data were acquired are assumed to be relatively homogeneous. For instance, the East London station was selected by assuming that the weather characteristics at Bisho is similar. Moreover, according to [44], the estimates of GPD parameters based on informative priors are susceptible to the choice of weather stations used to formulate the priors. Therefore, selecting appropriate nearby weather stations for the formulations of informative priors is also an important task in order to obtain valid and robust estimates of seasonal and station effects of the GPD parameters.

Estimation of Combined Seasonal and Station Effects of the GPD Parameters
So far, the seasonal and station effects were investigated for the GPD parameters and the results are provided for each effect. To see the advantages of introduction of stations and seasonal effects in the GPD ( j,i , j,i ) models over the maximum likelihood (ML), each seasonal and stations daily maximum and minimum temperatures data were analysed separately using the ML, that is assuming the model parameters are not random. Comparison of the posterior means and standard deviations from the station and seasonal effects model with MLEs, and standard errors obtained from a maximum likelihood analysis, applied separately to each station and season, are provided in Tables 4, 5, 6, 9, 10, and 11 in Appendix for seasonally clustered daily maximum and minimum temperatures data. The combined seasonal and station effects for the Bayesian GPD parameters are obtained using both non-informative and informative priors. For the combined effect of station and seasonal effects GPD parameters, the trace plots of the chain, for example, for daily maximum temperature data of East London station during summer period in January and February are given in Fig. 11 in Appendix. The plots depicted the parameters converged well. Similarly, for all other seasons Table 3 Posterior means (standard deviations) for the seasonal effect parameters using non-informative and informative priors and stations, convergence was checked and was achieved for all seasonal and station effects GPD parameters. As can be seen from the results, the posterior standard deviations for the scale and shape parameters from informative priors are lower than the corresponding standard errors of MLEs and non-informative priors for both minimum and maximum temperatures data. The effect of using station and seasonal effects on the parameters of GPD can be assessed using the values of the standard deviations obtained from the three methods. This might illustrate one apparent advantage of the introduction of station and seasonal effects into the parameters of GPD over the analysis employed using the standard likelihood approach. Hence, it appears that the pooling of information is seen to add the precision of the analysis for the parameters of GPD model considerably. Most importantly, the results revealed that the posterior Fig. 3 Posterior densities of the generalised Pareto distribution for station and seasonal effects parameters to daily maximum temperature data using non-informative (Non-inf) and informative priors (Inf) Fig. 4 Posterior densities of the generalised Pareto distribution for station and seasonal effects parameters to daily minimum temperature data of Durban during winter period June using non-informative (Non-inf) and informative (Inf) priors standard deviations of both priors are smaller than the standard errors of the corresponding MLEs. The results are also consistent over all stations and seasons. However, the precision of the parameter estimates are obtained fully when MCMC standard errors are compared with the standard errors of the parameter estimates from MLEs. Note that the naive standard errors for the Bayesian estimates are given by dividing the actual standard deviations by the number of iterations [50]. This is because the Bayesian approach allows for an additional source of variation, which implies that the parameters have probability distributions with hyper parameters. Furthermore, it is also vital to visualise that all of the comparisons that were drawn relate to Bayesian inference applied to the station and seasonal effects model versus standard likelihood, MLEs, inferences are estimated from separate stations and seasons data. Any gain in performance and difference in the estimates that were attributed to the station and seasonal effects model were confounded with possible differences for the GPD model fitted to separate seasons and station for the likelihood method. Therefore, the results of the present study support the findings of [24] that the introduction of station and seasonal effects in the parameters of extreme value distributions increases the precision of statistical inferences. Additionally, to see the overall effects of the informative priors formulated based on seasonally clustered historical data of nearby weather stations on the combined GPD seasonal and station effects parameters ( j,i , j,i ), the posterior densities of the combined seasonal and station effects parameters found through the elicited priors and the flat priors were compared. The estimated posterior densities of GPD ( j,i , j,i ) combined seasonal and station effects parameters for East London during summer period in January obtained from daily maximum temperature data and for Durban during winter period in June obtained from daily minimum temperature data are plotted and presented in Fig. 5. The posterior distributions of combined seasonal and station effects parameters obtained based on the non-informative and informative priors are symmetric for this specific station and seasonal period given in Fig. 5. On the other hand, the posterior densities of GPD ( j,i , j,i ) combined seasonal and station effects parameters for the informative priors that were built based on information from nearby weather stations of seasonally clustered data resulted in high peaks compared with that of the non-informative priors. Similar effects on the peaks of the posterior densities were observed for the GPD seasonal and station effects parameters across all stations and seasons for both separate and combined effects to daily maximum and minimum temperatures data. These results corroborate [44,51], for which the effect assessment made between the informative prior and non-informative priors.

Estimation of Return Levels
The seasonal and station effects GPD parameters estimated in Sect. 3.1 by themselves provide little information on the occurrence probability of extreme events. A more meaningful and also more relevant quantity for risk assessment is the probability of the observed variable (here, seasonally clustered daily maximum and minimum temperatures data) exceeding a certain level. With the estimates of the parameter vector Θ = ( j , i , j , i ) for every stations and seasons, we can use the seasonal and station effects GPD scale and shape parameters to derive return levels. Hence, the vectors of observations from the posterior densities of GPD seasonal and station effects parameters could be substituted into the return level d p in expression (5) to obtain realisations from the posterior distributions of return levels.
For instance, Fig. 6 shows a map of 50-year and 100-year return levels for seasonally clustered summer months across all stations. The return level for each station is represented in a colour scale, i.e. the continuous colour-scaled background. We show the interpolation only to illustrate regional patterns. The observed 50-year return levels range from about 30 • C of daily maximum temperature to almost 35 to 48 • C of daily maximum temperature. The highest 50-year and 100-year return levels were observed for stations in the Eastern and Western Cape province of South Africa, namely George, Plettenberg and Paarl during the summer months in December, January, February and March. Specifically, the highest daily maximum temperature that is exceeded on average once in 50 and 100 years period is more than 45 • C and it is observed during summer period in March for George and Plettenberg weather stations.  Similarly, map of 10-and 100-year return levels for seasonally clustered winter months across all stations for daily minimum temperature data is provided in Fig. 7. The return levels of daily minimum temperature data for each station are represented in a colour scale, i.e. the continuous colourscaled background. The observed 10-and 100-year return levels range from about 6 • C of daily minimum temperature to 30 • C of daily minimum temperature. The highest 10-and 100-year return levels were observed for Plettenberg station throughout all winter periods from June to August. The lowest daily minimum temperature that is exceeded on average once in 10 and 100 years period is below 8 • C and 10 • C, respectively, and it is observed during winter period in June for Paarl weather station. Furthermore, Figs. 8 and 9 display boxplots of the posterior 100-year return levels for seasonally clustered daily maximum and minimum temperatures data, respectively. The return level predictions estimated for each station during summer and winter periods. The estimated return levels are visible, and the difference observed between stations and seasons appeared to be considerable. In general, the results observed from the plots indicated that the seasonal and station effects for the GPD model predicted the return level well for daily maximum and minimum temperatures data and provided important inferences involving the return levels across spaces over seasonally clustered data. The station and seasonal effects obtained from fitting GPD model in the current study also supported Plettenberg stations. On the other hand, the minimum posterior 100-year return levels for daily minimum temperature data are observed during winter period in June across all weather stations considered in the current paper except for Hermanus and George stations. The current finding is in accordance with the findings of [24,44], in which the seasonal and station effects GPD model to extreme temperature data enabled increased precision for inferences capturing station and seasonal variability. Inferences for the quantities that are of most interest, return levels, obtained in this study revealed the frequency of occurrence of extremely high or low temperatures observed during summer and winter periods. The current results obtained may be useful to agricultural and energy sectors based on the fact that the frequency analysis of the occurrence of extremely high or low temperatures is directly relevant to the modelling of electricity demand as well as agricultural productivity in which the frequency and the duration of hot and cold temperatures are important aspects from the viewpoint of the disaster risk analysts. Hence, understanding seasonal characteristics of the temperature and its extremes forms part of preparedness and mitigation measures that can be put in place to reduce the impacts of extreme events. Over the majority of land regions, the projected intensity increases, and relative frequency increases tend to be larger for more extreme hot temperature than for weaker events [52]. This is also the case for the current paper in which the increase in extreme temperature is higher for maximum temperature than for the minimum temperature data. It is also important to emphasise that extreme events occurring during the summer period would have the most dramatic impact on plant productivity [53].

Simulation to Study the Seasonal and Station Effects
The analysis of the extreme temperature data in Sect. 3.1 showed clear discrepancies in precision for the estimation of GPD parameters between the estimation methods which consider seasonal and station effects, and GPD fitted to extremes without considering the station and seasonal effects. The aim of this section is to use simulated data to compare and investigate the precision of estimators for the GPD parameters obtained under the two estimation approaches.

Simulating Stations with Various Temporal Dependence Series
To generate a dataset for the GPD that consider station and seasonal effects, different simulated stations with various temporal dependence structure between the series of observations were simulated. Then, the simulated data series were partitioned into different cycles from 1 to 12 to form seasonality for all simulated stations. The seasonality was formed ensuring that the cycles from 1 to 12 are correctly identified as months and this was done using R packages Zoo, hydroTSM and SpatialExtremes. The real dataset that was considered in Sect. 3.1 exhibits strong seasonality, and only selected months were employed in the analysis. Similarly, selected consecutive cycles among the 12 seasonally partitioned simulated data are used to demonstrate the estimation approaches. Therefore, to simulate temporally dependent extremes, a reasonable assumption for modelling is to form the series of observations x 1 , ..., x n a stationary first-order Markov chain [24]. This assumption provides a simple model for many serially correlated environmental variables [24]. The stochastic properties of such a chain are completely determined by the joint distribution of successive pairs of observations. Given a model S(x i , x i + 1 ; ); i = 1, 2, ..., n − 1 , where is the vector of GPD parameters, the joint density function for x 1 , ..., x n is given by where S(x i ; ) is the limiting distribution of exceedances, GPD, given in expression (2), whereas the numerator, S(x i , x i + 1 ; ) , in expression (15) is obtained by the bivariate extreme value distribution (BEVD), thereby capturing the first-order Markov structure. The BEVD, using logistic model dependence function, is used to simulate temporally dependent series of observations. Logistic model is one of the most flexible and accessible dependence model of BEVD. The dependence parameter denoted by ∈ (0, 1) under the logistic model measures the strength of dependence between consecutive extremes. Independence and complete dependence for the pair (x i , x i+1 ) are obtained when = 1 and in the limit as → 0 , respectively (see [37,44,54], for details on BEVD and logistic model).
To generate (x i , x i + 1 ) from the limiting BEVD, various dependence structure between x i and x i+1 were contemplated to create simulated weather stations with different temporal dependence structures between a series of observations. Therefore, at the beginning was set to have values between 0 and 1, i.e. = 0.05, 0.15, 0.3, 0.4, 0.5, 0.6, 0.8 . These values provide various simulated weather stations which consist of lower to higher degree of dependence between the series. The following procedure was used to simulate (x i , x i + 1 ) ; 1. For = 0.05 , simulate the first observation, x 1 , from the GPD marginal distribution with known = 0.3 and = −0.45.

2.
To Simulate x 2 , form the conditional distribution of x 2 |x 1 using the BEVD logistic model and the simulated observation, x 1 , in step 1. 3. Use the simulated observation, x 2 , to form the conditional distribution of x 3 |x 2 again using the logistic dependence model, to simulate x 3 from this distribution. 4. Use the simulated observation, x 3 , to form the conditional distribution of x 4 |x 3 using the logistic dependence model, and simulate x 4 from this distribution.
, 5. Continue forming and sampling from new conditionals until the required number of observations (length n) have been simulated. 6. Repeat steps (2)(3)(4)(5)  The different values of the GPD shape parameter, , were used to reflect the various tails which might be observed in a real-life dataset, whereas the GPD scale parameter was held unit constant. Also, since the different criteria about what constitutes an extreme value for each combination of season and station are varying, the threshold value was allowed to vary across seasons and stations.

Simulation Study Result
For the simulated chain (of length n), the GPD parameters were estimated with and without seasonal and station effects. The estimation procedure for the simulation was repeated N = 1000 times, with chains of length n = 40, 000 , to create sampling distributions for the GPD parameters. These sampling distributions are then compared to the true parameter values used to simulate the chains in the first place.
The estimated values associated with the standard deviations (in brackets) for the GPD parameters estimated with and without station and seasonal effects of the simulated data are presented in Tables 7 and 8. All the results from the simulated data in Tables 7 and 8 clearly show the pitfalls of simply considering station and seasonal effects. The estimated standard deviations for each parameter are considerably smaller than the standard deviation of the corresponding parameters obtained without the effects. All the analyses using simulated data that were fitted considering the effects consistently underestimated the standard deviations as compared to the models without the effects.
Similar to Table 7, Fig. 10 depicts the actual sampling distributions for the GPD parameters and , which were estimated with and without the effects when is 0.15 and 0.3, for simulated stations 2 and 3, respectively. The sampling distributions obtained without the effect are located away from the true values, in contrast to the analysis that considered the effect, where the distributions are much more close to their target values. The 95% confidence intervals (not shown here) for some of the GPD parameters that were estimated without the effect even fail to include the true parameter values, particularly for stations which were simulated with low temporal dependence. These indicate the inadequacy of the estimations approach without the effect relative to the approach which uses the effect. The inadequacy of the estimations approach without the effect could have major practical implications.
In an extreme temperature setting, designing a level without specifying the station and seasonal effects in to the analysis could result in considerable under-protection. Average daily minimum temperatures are very rare in South Africa and only occur about 8 times in a year, and result in huge increases in electricity demand [20]. Specific to the current study, the lowest and highest daily minimum and maximum temperatures that exceeded on average once in 10 and 50 years period are below 8 • C and more than 45 • C, respectively. These levels are observed during winter period in June for Paarl and during summer period in March for George and Plettenberg weather stations (see Figs. 6 and 7). Hence, an investigation of expected cooler or warmer than typical years by specifying the station and seasonal effects in to the analysis is important. All our findings support our main conclusion that the introduction of station and seasonal effects in the parameters of extreme value distributions increases the precision of statistical inferences. This helps planners and decision-makers in the electricity sector.

Conclusions and Recommendations
In this paper, the authors introduced the seasonal and station effects on the distributions of extreme values, i.e. GPD, for modelling extreme temperature data at selected weather stations in South Africa, applying mainly the Bayesian approach. The modelling of GPD to daily maximum and minimum temperatures data has been employed here as a natural framework for modelling the station and seasonal effects on the parameters of GPD in terms of seasonality and station variations that are inherent in the data under consideration. The GPD models were fitted to the complex structures of the data and the parameters were estimated through the models, allowing the sharing of information between stations using seasonally clustered data. Non-informative and informative priors were used for estimation of the parameters. This estimation resulted in apparent advantages in terms of improved precision obtained in the estimation of the parameters for the GPD station and seasonal effects models. However, the posterior means for these parameters, which were obtained using informative priors, are very close to the non-informative priors. The posterior standard deviations of the station and seasonal effects parameter estimates were smaller for the informative priors compared with those from the non-informative priors.
Furthermore, the combined station and seasonal effects parameter estimates resulted in smaller standard deviations compared with the standard errors of MLEs fitted separately to stations and seasons. This reduction in standard deviations reflects the decrease in uncertainty owing to the inclusion of the station and seasonal effects parameters in the model as well as the informative priors. In general, the expected benefit of Bayesian analysis and the station and seasonal effects approach is the improvement in precision of parameter estimates over the MLEs for daily maximum and minimum temperatures data. Overall, the results depict that the expected gain in precision was fully obtained for the station and seasonal effects parameters of GPD. This provided apparent benefits in estimating the predictions for the quantities that are of most interest, return levels, which incorporate both the randomness of the process and model uncertainty. All the results obtained in this study were found to be analogous with literature on the topic.
In general, the results of this work are particularly important to quantify the seasonal and station effects of daily minimum and maximum temperatures extremes among the various meteorological stations. This was demonstrated using simulated data to compare and investigate the precision of estimators for the GPD parameters obtained with the effects relative to the estimation approach without the effect. All our findings support our main conclusion that the introduction of station and seasonal effects in the parameters of extreme value distributions increases the precision of statistical inferences. Several developments in the extreme value theory were motivated by applying univariate method to environmental datasets. Introducing seasonal and station variabilities to the univariate method, which allows information to be pooled across stations and seasons, should be preferred over the standard method resulting in improved precision of the estimates. The recommendation here is thus to introduce seasonal and station effects to the method in order to make proper inference when modelling the extremal process. The problems related to dependencies are relatively easy to overcome using the method employed in the current paper. Also, the approach appears to work extremely well when estimating return levels for stations with contrasting patterns of extremes. The modelling approach proposed in the current paper will help to reveal some useful information needed for planning by climatologists, meteorologists, agriculturalists, decision-makers and planners. Moreover, a number of recent developments in extremal methods have not yet been fully exploited in such applications. This information could also help to focus our attention on the multivariate and spatial modelling of extreme temperature data. The inference on the extremal behaviour of extreme temperature over several weather stations could be improved with these methods and this is the subject of future research.    Funding Open access funding provided by University of South Africa.

Data Availability
The data that support the findings of this study are available from South African Weather Service (SAWS) but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of SAWS.
Code Availability R Code is available and can be obtained on reasonable request from the authors.

Declarations
Ethics Approval Ethics approval was granted by the University of South Africa School of Science Ethics Review Committee (Date: July 01, 2019; No: 2019/SSR-ERC/028).

Conflict of Interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.