Bayesian Generalised Additive Extreme Value Modelling of Maximum Wind Speeds in a Semi-Arid Region

1. Introduction

Climate change is a major concern facing the world today, with extreme events increasingly having impacts on property, infrastructure, economy, and agriculture. These events are broadly defined as the time and location where weather, climate, or environmental conditions exceed the usual range of what has occurred in the past [1]. Within South Africa, the Northern Cape province stands out as a semi-desert province where strong winds frequently affect the communities, threaten infrastructure and agricultural productivity. Considering these significant impacts of extreme events, it is therefore imperative to study their patterns to help better preparedness.

Extreme event studies have increasingly become important due to increased changes in the world climate. Different extreme value and additive frameworks which have been mostly useful in modelling and forecasting extreme weather events include the Generalised Extreme Value (GEV) [2], Generalised Additive Models (GAMs) and Bayesian Generalised Extreme Value (GEV) models, among others [3]. Recent works has employed a hybridisation of the framework called Generalised Additive Extreme Value (GAEV) to model extreme events by allowing extreme value parameters to depend on covariates through spline functions [4]. However, such applications remain limited and scattered across hydrology and environmental science. Furthermore, while Bayesian methods have been increasingly recognised for their ability to incorporate prior information and quantify uncertainty [5,6,7], comparative studies among Bayesian and frequentist approaches in general and extreme wind speed modelling context remain scarce, particularly in semi-arid regions. In South Africa, and specifically in semi-arid regions such as the Northern Cape, extreme wind events pose significant risks to infrastructure, agriculture, and livelihoods. Yet, empirical studies applying advanced statistical frameworks to localised wind extremes are limited. From the limited studies, very few have used hybrid GAEV which offer better predictive performance by adequately capturing the underling properties of the extreme events, but the focus has been limited to temperature and rainfall extremes. In addition, comparative studies of the performance of the hybrid GAEV models estimated via maximum likelihood and Bayesian estimation approaches with applications to extreme wind speed remain scarce, particularly within semi-arid regions of South Africa. Therefore, there is a need to fill this gap, owing to the need for accurate modelling and forecasting of wind speed in semi-arid regions.

This study addresses these gaps by extending and applying the hybrid GAEV framework to a new and practically important problem; extreme wind speed modelling in semi-arid South Africa and then provide a comparison between Bayesian and frequentist estimation within this framework. The extension of the GAEV framework allows the parameters of the GEV distribution (location, scale, and shape) to vary smoothly with covariates such as seasonality, temperature, humidity, and cloud cover via spline functions. The method is then used to fit wind speed data from Kimberley and Upington stations, to the frequentist and Bayesian versions of the hybrid model. Subsequently, the performance of the two approaches is compared to identify the most accurate method.

The originality of this study lies in its formalization or extension of the GAEV framework to explicitly use splines to allow the location, scale, and shape parameters to vary smoothly with covariates such as seasonality, temperature, humidity, and cloud cover to capture non-linear influences on extreme wind behaviour more effectively. A further novel contribution is the explicit comparison between frequentist and Bayesian estimation of GAEV models. By employing the maximum likelihood estimation alongside Bayesian inference via Markov Chain Monte Carlo (MCMC). We then compare the performances of the frequentist and Bayesian approaches. This comparative evaluation is relatively underexplored in African climate modelling literature and demonstrates the practical advantages of Bayesian approaches for forecasting extremes.

Finally, the application of this framework to extreme wind speed events in Kimberley and Upington, represents a significant regional innovation. The Northern Cape is a semi-arid, climate sensitive region where infrastructure, agriculture, and livelihoods are highly vulnerable to wind extremes. By employing advanced statistical modelling within this context, the study contributes to region specific insights that can inform disaster preparedness, infrastructure planning, and climate adaptation strategies. For the rest of the paper, the materials and methods are presented in Section 2, the results and discussions in Section 3 while the conclusion is presented in Section 4.

2. Materials and Methods

2.1. Study Source and Area

The data used in this paper consist of the maximum daily temperature, maximum relative humidity, maximum cloud cover, and maximum wind speed from 1 January 2019 to 31 December 2024. The data was specifically sourced from the Open-Meteo API which is available at https://open-meteo.com/en/docs/historical-weather-api. The specific areas of interest are Kimberley and Upington, located in the Northern Cape Province of South Africa. These areas have predominantly semi-arid climates, see Figure 1, which makes them appropriate for this study.

The dataset consisted of hourly observations, which were aggregated to daily maximum wind speed to have a block maxima data for the extreme value analysis. The aggregation of hourly observations to daily maximum resulted in a loss of information regarding the "within day" variability and the timing and duration of extreme wind events.

2.2. Test for Stationary

As highlighted by [8], assessing stationarity and the existence of unit roots remains a fundamental step in time series analysis because these properties directly affect model specification, estimation, and forecasting accuracy. Since many statistical time series techniques are developed under the assumption that the underlying process is stationary, it is important to verify this assumption before proceeding with further analyses in the present study. To evaluate stationarity, three widely used procedures are employed, namely the Augmented Dickey–Fuller (ADF), Phillips–Perron (PP), and Kwiatkowski–Phillips–Schmidt–Shin (KPSS) tests. The Augmented Dickey–Fuller (ADF) test is a commonly applied method for examining whether a time series contains a unit root and is therefore non-stationary [9]. By accounting for higher-order serial correlation through lagged difference terms, the ADF test provides a robust framework for determining the stationarity properties of a series. The hypotheses considered in the test are given as follows:

• $H_0$: the time series contains a unit root, non-stationary;
• $H_a$: the time series does not contain a unit root, stationary.

The ADF test uses a regression model:

$$\Delta y_t=\alpha +\beta t+\gamma y_{t-1}+{\displaystyle \sum _{i=1}^p}{\theta }_i\Delta y_{t-i}+{ϵ}_t,$$

(1)

where $\alpha$ is a constant, $\beta$ is the coefficient on time trend, and $ϵ$ represents the white noise process. The null hypothesis is $\gamma =1$, and the alternative hypothesis is $\gamma \ne 1$.

The PP test is closely related to the ADF test. It serves as an alternative to the ADF test and applies non-parametric corrections to account for serial correlation and heteroskedasticity in the error terms [10].The hypothesis for the PP test is formulated as indicated below:

• $H_0:$ time series has a unit root, non-stationary;
• $H_a:$ time series does not have a unit root, stationary.

The test involves the estimation of the regression equation:

$$y_t=\alpha +\beta t+\rho y_{t_1}+{ϵ}_t,$$

(2)

where $y_t$ represents the time series, $\alpha$ is the intercept, $\beta t$ is a deterministic time trend, $\rho$ is the autogressive coefficient, and ${ϵ}_t$ is a white noise error term.

The KPSS test models a time series $y_t$ as:

$$y_t=r_t+{\beta }_t+{ϵ}_t,$$

(3)

where ${\beta }_t$ is a deterministic trend, ${ϵ}_t$ is a stationary error term, and $r_t$ is a random walk which is then expressed as

$$r_t=r_{t-1}+u_t,$$

(4)

where $u_t\sim IID(0,{\sigma }_u^2)$.

The KPSS test assumes the time series data is stationary under the null hypothesis and non-stationary under the alternative hypothesis, i.e.,

• $H_0:$ time series is stationarity;
• $H_a:$ time series not stationarity.

If ${\sigma }_u^2=0$, then $r_t$ is constant and time series is stationary. If ${\sigma }_u^2>0$, then time series is non-stationary and $r_t$ follows a stochastic trend.

The KPSS test statistic is calculated as:

$$KPSS={\displaystyle \frac{1}{T^2{\widehat{\sigma }}^2}}{\displaystyle \sum _{t=1}^T}{S}_t^2,$$

(5)

A high KPSS statistic value suggests the presence of a unit root and a low KPSS value means the null hypothesis of stationarity canoot be rejected.

2.3. Normality Tests

The normality tests are usually used for checking normality, that is, if the data are normally distributed or not. This paper used the Shapiro-Wilk (SW) and Jarque-Bera (JB).

SW tests the null hypothesis that the data are normally distributed against the alternative that it is not. According to [11], the null hypothesis is rejected when then value of W statistic is small and the p-value is less that the specified level of significance which indicates a significant deviation from normality.

The SW test is given as:

$$W={\displaystyle \frac{{\left({\sum }_{i=1}^n{a}_i{X}_i\right)}^2}{{\sum }_{i=1}^n{(X_i-\overline{X})}^2}},$$

(6)

where $x_i$ is the ith order statistic and $\overline{x}={\displaystyle \frac{x_1+...+x_n}{n}}$ is the sample mean.

The coefficient $a_i$ are given by:

$$(a_1,...,a_n)={\displaystyle \frac{m^T{V}^{-1}}{{\left(m^T{V}^{-1}{V}^{-1}m\right)}^{\frac{1}{2}}}},$$

(7)

where

$$m={(m_1,...,m_n)}^T$$

(8)

$m_1,...,m_n$ are the expected values of the order statistic of independent and identically distributed $\left(IID\right)$ random variables from the standard normal distribution and V is the covariance matrix of the order statistic.

The JB is also a normality test that is widely used due to its simplicity and satisfactory performance [12]. The test is defined as:

$$JB={\displaystyle \frac{n}{6}}(S^2+{\displaystyle \frac{1}{4}}{(K-3)}^2),$$

(9)

where n is the number of observations , S is the sample skewness, and K is the sample kurtosis:

$$S={\displaystyle \frac{\frac{1}{n}{\sum }_{i=1}^n{(x_i-\overline{x})}^3}{{\left(\frac{1}{n}{\sum }_{i=1}^n{(x_i-\overline{x})}^2\right)}^{\frac{3}{2}}}}.$$

(10)

The hypothesis test:

• $H_0:$ the series follows a normal distribution;
• $H_a:$ the series does not follow a normal distribution.

2.4. Goodness-of-Fit

The Goddness-of-Fit (GoF) tests statistics are usually used for checking the validity of an assume probability distribution model, i.e to evaluate how well the GEVD explains extreme event. This paper used the Kolmogorov-Smirnov (KS) test the Anderson-Darling (AD). According to [13] the KS test is a non-parametric test applied to test whether the sample under consideration is from a hypothesised distribution or to compare whether two samples come from identical distribution. The hypotheses are formulated as:

• $H_0:F\left(x\right)=F^{*}\left(x\right)$ for all x from $-\infty$, i.e., the series follows a specified distribution;
• $H_a:F\left(x\right)\ne F^{*}$ the series does not follow the specified distribution.

The test statistic is also defined as:

$$D_{m}ax=Ma{x}_x|F^{*}\left(x\right)-F\left(x\right)|$$

(11)

The null hypothesis for this test is rejected at 5% significance level if $D_{m}ax$ calculated is greater than the tabulated value $D_{0.05}={\displaystyle \frac{1.36}{\sqrt{n}}}.$

Meanwhile, AD test is used to assess whether data comes from a specified distribution with the hypotheses formulated as

• $H_0:$ the series follows a specified distribution;
• $H_a:$ the series does not follow a specified distribution.

The test statistic is defined as:

$$A^2=-n-{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}(2i-1)[lnF\left(x\right)+ln(1-F_n\left(x\right))],$$

(12)

where the data x is the ordered sample of size n, and $F\left(x\right)$ is the theoretical CDF to which the sample is compared. The A-D test assign more weight to the tail of the distribution than the K-S test.

2.5. Generalised Additive Models

Generalised Additive Models (GAMs) are an extension of the generalised linear models (GLM)

$$g\left(E\left[Y\right]\right)={\beta }_0+{\beta }_1{X}_1+...+{\beta }_k{X}_k,$$

(13)

where Y is the response variable, $E\left[Y\right]$ is the expected value of the response, $g(.)$ is the link function (connects the mean of Y to the linear predictor), ${\beta }_0,{\beta }_1,...,{\beta }_k$ are the model coefficients and $X_1,...,X_k$ are the independent variables. The GAMs allow non-linear relationships between predictors and the response through the use of smooth functions [14]. GAMs extend Generalised Linear Models (GLMs) by incorporating smooth functions of the predictors.

Let Y be a response random variable and $X=(X_1,...,X_k)$ be a set of k independent variables. GAMs consists of a random component, an additive component, and a link function relating to the two component,

$$g\left(E\left[Y\right]\right)={\beta }_0+s_1\left(x_1\right)+...+s_k\left(x_k\right)={\beta }_0+{\displaystyle \sum _{i=1}^k}{s}_i\left(x_i\right),$$

(14)

where $g(.)$ is a link function, $E\left[Y\right]$ is the expected value of the response variable, ${\beta }_0$ is the constant, $x_i$ are the independent variables and $s_i\left(x_i\right),i=1,...,k$ are smooth functions of the predictor variable. The smooth functions $s_i\left(x_i\right)$ are usually fit using non-parametric curve estimating methods such as the local-regression, splines, or smooth splines.

GAMs focuses on exploring data non-parametrically and are more suitable for exploring the data set and visualising the relationship between the dependent variable and the independent variables [15].

2.6. Extreme Value Theory

Extreme Value Theory (EVT) is known as the study of extremal properties of random processes and its objective is to quantify the stochastic behaviour of a process as usually high or low levels. According to [16], the intention to use the EVT is to predict the occurrence of extreme or rare events using historical data. Block maxima and the peaks-over-threshold (POT) are the two fundamental approaches in EVT.

2.6.1. Generalised Extreme Value Distribution

The Generalised Extreme Value Distribution (GEVD) arises from the theory of extreme values and provides a unified framework for modelling rare and extreme events [17]. It is typically applied using the block maxima approach, where the largest observation within each predefined time block is extracted for analysis. Under the extreme value theorem, the distribution of suitably standardised block maxima converges to the GEVD as the sample size becomes large. Consequently, the GEVD is well suited for analysing maxima obtained from independent and identically distributed ($i.i.d$.) observations collected over fixed intervals. As noted by [18], the GEVD is frequently adopted as an approximate probability model for the largest values observed within sufficiently long sequences of random variables. A notable feature of the GEVD is that it encompasses three classical extreme value distributions, namely the Gumbel, Fréchet, and Weibull distributions, commonly referred to as Type I, Type II, and Type III extreme value distributions, respectively. By incorporating these distributions into a single family, the GEVD is capable of representing a wide variety of tail behaviours and distributional shapes through a common modelling framework.

$$G_{\xi }(x,\mu ,\sigma ,\xi )=exp\left\{-{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]}^{-1/\xi }\right\},\mathrm{if}1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)>0\mathrm{and}\xi \ne 0,$$

(15)

where $\mu$, $\sigma$ and $\xi$ are the location, scale and shape parameters, respectively. For $-\infty <\mu <\infty$, $\sigma >0$ and $-\infty <\xi <\infty$. For $\xi =0$ Equation (15) is the Gumbel class of distribution which is given by

$$G_{\xi }(x,\mu ,\sigma ,\xi )=exp\left\{-exp\left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\right\},\mathrm{if}\xi =0.$$

(16)

For $\xi >0$ and $\xi <0$, the distribution in Equation (15) corresponds to the Fréchet and negative Weibull distributions, respectively. The location parameter specifies where the distribution is centred, the scale parameter reflects the degree of variability in the data, and the shape parameter influences the nature of the distribution’s tails.

The GEVD’s extreme quantiles are estimated as follows:

$$z_p=\left\{\begin{array}{c}\mu -{\displaystyle \frac{\sigma }{\xi }}\left[1-{\left\{-log(1-p)\right\}}^{-\xi }\right],\xi \ne 0\hfill \\ \mu -\sigma log\left[-log(1-p)\right],\xi =0,\hfill \end{array}\right.$$

(17)

Here, ${\displaystyle \frac{1}{p}}$ represents the return period and $z_p$ is known as the return level associated with the return period ${\displaystyle \frac{1}{p}}$. This is the level which is expected to be exceeded on average once every ${\displaystyle \frac{1}{p}}$ years. The GEVD allows us to model the extreme values of block maxima.

2.6.2. Generalised Pareto Distribution

Within the Peaks-Over-Threshold (POT) framework, the Generalized Pareto Distribution (GPD) is commonly employed to describe observations that exceed a suitably chosen threshold. In extreme value analysis, it is particularly useful for characterizing the tail behaviour of an underlying distribution. The GPD has three parameters, the threshold parameter $\mu$, the shape parameter $\xi$ and the scale parameter $\sigma$. The GPD is considered the most efficient distribution compared to the GEVD, the equation is given by:

$$G_{\xi }(x,\mu ,\sigma ,\xi )=\left\{\begin{array}{c}1-{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]}^{-1/\xi },\mathrm{for}1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)>0\mathrm{and}\xi \ne 0\hfill \\ 1-exp\left(-{\displaystyle \frac{x-\mu }{\sigma }}\right),\xi =0,\hfill \end{array}\right.$$

(18)

where $\sigma >0$ and $x>\mu$. Equation (18) becomes an exponential distribution when $\xi =0$ and uniform when $\xi =-1$. When $\xi <0$, the distribution has a finite end point and is referred to as short-tailed. When $\xi >0$, the distribution is referred to as heavy-tailed.

2.6.3. Parameter Estimation

By employing maximum likelihood and Bayesian method to estimate the parameters $\mu ,\sigma ,\mathrm{and}\xi$, we will gain a valuable insight into the behaviour of the extreme events, under study.

The Maximum Likelihood Estimation Method

Parameter estimation under the maximum likelihood framework is achieved by identifying the set of parameter values that maximizes the likelihood function associated with the assumed probability distribution [19]. It is a widely used method for estimating parameters in statistical models [20] and it is an effective method when dealing with large samples. However, it is know to performs badly when the sample size is small.

Consider a sample $X_1,...,X_n$ of independent and identically distributed observations drawn from a GEVD. For $\xi \ne 0$, the log-likelihood funct is given by:

$$\ell (\mu ,\sigma ,\xi )=-nlog\sigma -\left(1+{\displaystyle \frac{1}{\xi }}\right){\displaystyle \sum _{i=1}^n}log\left(1+\xi \left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)\right)-{\displaystyle \sum _{i=1}^n}{\left(1+\xi \left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)\right)}^{-1/\xi },$$

(19)

where $i=1,...,n$ and $1+\xi \left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)>0$. When $\xi =0$, the log-likelihood function becomes:

$$\ell (\mu ,\sigma )=-nlog\sigma -{\displaystyle \sum _{i=1}^n}\left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)-{\displaystyle \sum _{i=1}^n}exp\left(-{\displaystyle \frac{x_i-\mu }{\sigma }}\right).$$

(20)

Let $X_1,...,X_n$ be a random sample of n with extreme value $X_{(}{n}_{)}$ above a threshold u, the log-likelihood function of the GPD be expressed as

$$\ell (\sigma ,\xi )=-nlog\sigma -\left(1+{\displaystyle \frac{1}{\xi }}\right){\displaystyle \sum _{i=1}^n}log\left(1+\xi {\displaystyle \frac{x_i}{\sigma }}\right),\xi \ne 0,$$

(21)

provided $\left(1+\xi {\displaystyle \frac{x_i}{\sigma }}\right)>0$ for $i=1,...,n$. The likelihood function may be rewritten as

$$\ell (\sigma ,\xi )=-nlog\sigma -{\displaystyle \frac{1}{\sigma }}{\displaystyle \sum _{i=1}^n}{x}_i,\mathrm{when}\xi =0.$$

(22)

Bayesian Method

Bayesian estimation of GEVD parameters involves updating prior beliefs using the likelihood derived from observed data. This process yields a posterior distribution that reflects both prior assumptions and sample information, resulting in a more comprehensive understanding of the parameters [19]. The model parameters are not usually random, but they are, however, treated as random by assigning probability distributions to describe the uncertainty of their values. The Bayesian method uses a computational algorithm called Markov Chain Monte Carlo (MCMC) to estimate the posterior distributions of model parameters. The parameter estimation is made through the posterior distribution which is computed using Bayes’ Theorem:

$$p(\theta \mid y)\propto p(y\mid \theta )·p\left(\theta \right)p(\theta \mid y)={\displaystyle \frac{p(y\mid \theta )·p\left(\theta \right)}{p\left(y\right)}}={\displaystyle \frac{p(y\mid \theta )·p\left(\theta \right)}{{\int }_{\theta }p(y\mid \theta )p\left(\theta \right)d\theta }},$$

(23)

where $p\left(\theta \right|y)$ is the posterior distribution of $\theta$, $p\left(y\right|\theta )$ is the likelihood function, $p\left(\theta \right)$ is the prior distribution which represents prior knowledge about the parameters, and $p\left(y\right)$ is the marginal likelihood. The Bayesian GAMs estimation:

$$p({\beta }_0,\left\{f_i\right\}\mid \mathbf{y})\propto p(\mathbf{y}\mid {\beta }_0,\left\{f_i\right\})·p\left({\beta }_0\right)·\prod _{i}p\left(f_i\right),$$

(24)

where y is the independent variable being observed, $p({\beta }_0,f_i|y)$ is the posterior distribution of the intercept and smooth functions, $p\left(y\right|{\beta }_0,f_i)$ is the likelihood which describes how the data relate to the model, $p\left({\beta }_0\right)$ is the prior on the intercept, $p\left(f_i\right)$ is the prior on the smooth function, and ${\prod }_{i}p\left(f_i\right)$ is the product over all smooth function priors, assuming they are independent. Suppose $y_1,y_2,...,y_n$ are iid and their distribution fall within the GEV or GPD family. The parameters are treated as random variables.

$$p(\mu ,\sigma ,\xi \mid \mathbf{y})\propto p(\mathbf{y}\mid \mu ,\sigma ,\xi )·p\left(\mu \right)·p\left(\sigma \right)·p\left(\xi \right),$$

(25)

where $p(\mu ,\sigma ,\xi |y)$ is the posterior distribution of the GEVD/GPD parameter, $p\left(y\right|\mu ,\sigma ,\xi )$ is the likelihood, and $p\left(\mu \right)·p\left(\sigma \right)·p\left(\xi \right)$ are the prior distributions.

2.7. Generalised Additive Extreme Value

Rather than using the GAMs and EVT separately, this paper integrates them through the Generalised Additive Extreme Value (GAEV) model allowing the GEV parameters to vary smoothly with covariates. The GAEV model is a hybrid model that combines the flexibility of GAMs with the principles from EVT [4]. The GAMs framework allows us to model the non-linear effects of covariates on the GEVD parameters.

Let $Y\left(x\right)$ be a random variable that depends on the covariate x, such that

$$Y\left(x\right)\sim GEV\left(\mu \right(x),\sigma (x),\xi (x\left)\right),$$

(26)

where parameters $\mu \left(x\right),\sigma \left(x\right)>0,\mathrm{and}\xi \left(x\right)$ are modelled as smooth functions of covariates through GAMs.

The location model:

$$\mu \left(x\right)={\beta }_0+{\displaystyle \sum _{i=1}^k}{s}_i\left(x_i\right)$$

(27)

the log scale model:

$$log\left(\sigma \left(x\right)\right)={\alpha }_0+{\displaystyle \sum _{j=1}^k}{p}_j\left(x_j\right)$$

(28)

and the shape model:

$$\xi \left(x\right)={\rho }_0+{\displaystyle \sum _{m=1}^k}{h}_m\left(x_m\right),$$

(29)

where ${\beta }_0,{\alpha }_0,$ and ${\rho }_0$ are intercepts, $s_i(·),p_j(·)$ and $h_m(·)$ and $x_i,x_j$ and $x_m$ are non-linear functions that include seasonality (day of year), maximum temperature, maximum relative humidity, and maximum cloud cover.

The GAEV model is usually applied to model the distribution of block maxima, allowing the estimation of return levels.

2.8. Model Selection

Model selection involves choosing the best model from a list of candidates based on data analysis [21]. Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Leave-One-Out Cross Validation (LOO-CV) are selection models that are used to compare how well statistical models fit data.

AIC aims to find the model that best explains the data with the fewest number of parameters, hence balancing goodness of fit and model complexity. Given the number of parameters in the model as k and L as the estimated likelihood, the AIC is formulated as. A model with the lowest AIC value is preferred.

$$AIC=2K-2log\left(L\right),$$

(30)

Using similar definitions for k and L, the BIC is also fromulated as,

$$BIC=klog\left(n\right)-2log\left(\widehat{L}\right),$$

(31)

where n is the number of data points. BIC selects a model that maximises the posterior model probability. Unlike the AIC and the BIC which are used to compare the performance of any two models, the LOO-CV is used to assess the fitness of different estimated Bayesian models [7].

2.9. Model Diagnostics

Model diagnostics includes using various tools and techniques to check model validity and assumptions, they assess how well a chosen model fits the data [22] and whether it satisfies the necessary assumptions of the modeling framework. Another way of assessing the goodness-of-fit is by utilizing the diagnostic plots; Probability-Probability (P-P), Quantile-Quantile (Q-Q), Residuals vs. Seasonality, and Density Histogram.

The P-P plot compares the empirical probabilities to the theoretical probabilities. For a good fitted model, the empirical probabilities will overlap with the theoretical probabilities to form a straight diagonal line. plots [23]. Checks the overall distribution fit.

The Q-Q plot compares the empirical quantiles and the theoretical quantiles from the fitted model, deviations from the diagonal line shows a mismatch in tail estimation [23]. It checks the tail behaviour and overall fit of the GEV.

The residual vs seasonality plot examines whether seasonal structure has been adequately captured by the smooth terms.

The density histogram compares the empirical and fitted density of the maxima.

2.10. Metrics for Evaluating Forecasts

Forecast evaluation assesses how well the fitted model predicts unseen or future data. It helps with comparing and assessing model performance. Evaluating the effectiveness of our models, the following metrics are to be applied; Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE). The values of RMSE, MAE, and MAPE can measure how how good the prediction results produced by a model are [24]. Given n as the length of the sample, $y_i$ as the actual value and ${\widehat{y}}_i$ as the predicted value, the metrics are calculated as follows:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2},$$

(32)

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\left|y_i-{\widehat{y}}_i\right|,$$

(33)

$$\mathrm{MAPE}={\displaystyle \frac{100}{n}}{\displaystyle \sum _{i=1}^n}\left|{\displaystyle \frac{y_i-{\widehat{y}}_i}{y_i}}\right|,$$

(34)

3. Results and Discussion

3.1. Exploratory Data Analysis

3.1.1. Summary Statistics

Table 1 shows the descriptive statistics of both Kimberley and Upington stations. With Kimberly showing the maximum wind speed of $35.30\mathrm{km}/\mathrm{h}$ and average of $16.90\mathrm{km}/\mathrm{h}$. And Upington showing the maximum wind speed of $47.70\mathrm{km}/\mathrm{h}$ and an average of $22.30\mathrm{km}/\mathrm{h}$. The distribution of Kimberley a positive skewness since the value is positive and a platykurtic given the kurtosis value is more than 3 which implies more extreme wind speed events. While Upington also shows a positive distribution and leptokurtic with a kurtosis value of less than 3 which means fewer extreme wind speed events.

3.1.2. Stationarity Tests

Table 2 presents the results of the stationary test for the daily maximum wind speed for Kimberly, and Upington.

Evidence from the table suggests that the time series data are stationary for all the stations. This is because the p-values for the ADF, PP and KPSS tests are all less than the $(\alpha =0.05)$ level of significance. Overall, the ADF and PP results show strong evidence that the daily maximum wind speed data at all stations are stationary.

3.1.3. Test for Normality Results

Table 3 presents the normality test for the daily maximum wind speed. S-W and J-B test results show that the null hypothesis of normality is rejected at $0.05$ level of significance for the daily maximum wind speed across all stations. These results are consistent with the kurtosis values from Table 1. In conclusion, evidence from the table indicate that the daily maximum wind speed for both stations do not follow a normal distribution.

3.2. Model Fitting Results

3.2.1. Estimated Frequentist Parameters

The table below, Table 4, presents the parametric intercept estimates and their corresponding standard error (in parantheses) for the GEV parameters; the location, scale, and shape.

From the table, there is clear evidence that Upington has the highest location intercept, under the Location model, implying higher extremes, with 19.50, than Kimberley. The shape intercepts are negative for both stations, suggesting a Weibull distribution, upper bounded tails and finite maximum wind speed.

The tables below, Table 5 and Table 6, present smooth terms for Kimberley and Upington. The estimated degrees of freedom (edf) determines the flexibility of the curve; an edf value closer to 1 means the term behaves like a straight line (there is a linear effect) and an edf greater than 1 means the term captures non-linear patterns. Furthermore, the importance of these smooths are checked through the p-values, [4].

The result shows that day of year (seasonality trend) and temperature have highly significant smooth effects in the Location model. The day of year smooth term shows high significance across all models in both stations which indicates strong evidence of seasonal effects on the daily maximum wind speed.

In Kimberley, Table 5, all the smooth terms are significant across all models except for maximum cloudcover under the Scale model. The smooth terms showing non-linear relationships except the maximum cloudcover under the Shape model with edf = 1.01.

Table 6, the Location and Scale models have significant smooth terms. The smooth terms show approximately linear effects on the Shape model, however, maximum temperature seems insignificant with p-value more than 0.05.

Overall, the Location model seems to produce more highly significant terms over the Scale and Shape models. (see Appendix A).

3.2.2. Goodness-of-Fit Tests for Extreme GAMs

Table 7 presents the Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) goodness-of-fit test results for Kimberley and Upington stations using the mle parameters.

From Table 7, the p-values for both tests are higher than 5% thus the tests are statistically insignificant for the location model across both stations. This implies that the daily maximum wind speed follows the specified GAEVD model, that is, the Location model fits well. As for the Scale model Upington show highly significant p-value, leading to the rejection of the null hypothesis and the conclusion that the station does not follow the specified GAEVD model while Kimberley does follow the specified GAEVD model. The p-values for Kimberly shows insignificance under the Shape model which implies that the daily maximum wind speed follows the specified GAEVD model while Upington does not follow the GAEVD model since the p-values are all less than 5%. Similarly, the table provides clear evidence that the maximum daily wind speed follows the specified GAEVD model in Kimberly stations.

Figure 2 and Figure 3 below display four standard diagnostic plots for the best fitted model, Location model. From each figure, the top-left is the Probability-Probability (P-P) plot, top-right is the Quantile-Quantile (Q-Q) plot, bottom-left is the Residual vs Seasonality plot, and bottom-right is the Density Histogram of Residuals plot. The figures show that the P-P and Q-Q plots are almost linear illustrating a good fit of the GEV to the daily maximum wind speed. The density histogram shows a positive skewness across both stations.

3.2.3. Model Selection

Table 8 presents model selection. The best model is seelected using Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The lowest values of the AIC and BIC represent the best model.

The Location model produced lower values across both stations, this implies that Location model is the best fitted model, with Kimberley having the lowest AIC and BIC values. Since the Location model outperforms the Scale and Shape models, this indicates that changes in the central tendency are the best explanation for extremes, rather than variations in dispersion or tail heaviness.

3.2.4. Metrics Evaluation

Table 9 presents the RMSE, MAE, and MAPE evaluation metrics used to assess the forecasting performance of the Location, Scale, and Shape models, for both in-sample and out-of-sample. The model with the smallest RMSE, MAE, and MAPE is considered the best model. In the in-sample metrics evaluations performance, the Location model demonstrated better performance than the Scale and Shape models across all metrics. The results show that the Location model gives the most precise representation of extremes (better predictive accuracy), with lower RMSE, MAE, and MAPE values, while Scale and Shape models show weaker performances across both stations.

Kimberley produced the lowest metrics value under the Scale model, in the out-of-sample performance, indicating that the Scale model gives better predictive accuracy. Similarly, the Scale model performs better than the Location and Shape for Upington. Overall, Scale model produces lower metrics values across all stations, in the out-of-sample performance.

3.2.5. Return Level

The return level estimates for daily maximum wind speed at different return periods are given in Table 10, and the accompanying 95% confidence intervals are shown in parentheses.

The results in the table shows the return levels increase as the periods increase. Kimberley has the lowest 5 year return level of 19.51 km/h with a narrower confidence interval ranging from 19.30 km/h to 1973 km/h and Upington has 27.14 km/h return level with a bit wider confidence interval ranging from 26.76 km/h to 27.54 km/h. This suggests that a maximum wind speed of 19.51 km/h is expected to increase atleast once in 5 years. 41.71 km/h maximum wind speeds are expected to increase atleast once in 200 years in Upington and 27.68 km/h in Kimberley. From the table, it is clear that Upinton will experience more increase in the wind speed along the years as compared to Kimberly.

The figure, Figure 4, below depicts the return levels plots for the stations Kimberley and Upington.

3.3. Diagnostics of Estimated Bayesian Extreme GAMs

The performance and convergence of the Bayesian GAEV models were assessed using posterior probability plots, the Gelman–Rubin statistic (Rhat), Bulk ESS, Tail ESS, and LOO criteria [7]. Table 11 and Table 12 present the corresponding diagnostic results for the location, scale, and shape models. Since all estimated parameters, namely day of year, temperature, relative humidity, and cloud cover, yielded Rhat values of 1, there is strong evidence that the MCMC chains converged successfully and produced reliable posterior estimates.

3.3.1. Goodness-of-Fit Test for Extreme GAMs

Using the posterior parameter estimates derived from the Bayesian analysis, goodness-of-fit diagnostics were performed to evaluate the appropriateness of the fitted GAEVD model. In particular, the K–S and A–D tests were applied to examine whether the distribution of daily maximum wind speed observations is consistent with the assumed GAEVD model.

The goodness-of-fit statistics reported in Table 13 provide evidence regarding the suitability of the Bayesian GAEVD models for the two study locations. For the model governing the location parameter, both the K–S and A–D procedures produced probability values above the conventional 0.05 threshold for Kimberley and Upington. As a result, the assumption that the observed daily maximum wind speed data originate from the fitted GAEVD cannot be rejected, indicating satisfactory model performance at both stations. Different conclusions emerge for the scale component. While the Kimberley results remain compatible with the fitted distribution, the corresponding outcomes for Upington suggest a lack of agreement between the observed data and the assumed model. This is reflected by statistically significant test results, with reported probability values of approximately 0.005. Consequently, the hypothesis that the scale model adequately represents the Upington data is not supported. With respect to the shape component, the diagnostic tests again reveal contrasting behaviour between the two sites. The Kimberley station exhibits probability values exceeding the 5% significance level, providing no indication of a departure from the fitted GAEVD. In contrast, the corresponding values for Upington fall below the specified significance threshold, suggesting that the fitted shape model does not adequately capture the characteristics of the daily maximum wind speed observations at that location.

3.3.2. Model Performance

The table below, Table 14, presents the LOO-CV tests results. The prediction accuracy is measured by the expected logarithmic predictive density (elpd_loo) and leave-one-out information criterion (looic) compares the model(s). The highest elpd_loo and lowest looic values indicate better model (superior model) performance over other models, and the pareto estimates evaluates the reliability of the model where the values less than 0.7 indicate the results are stable and reliable.

From the table, there is clear evidence that the Location model outperforms the Scale and Shape model across both stations. Kimberley has a higher elpd_loo of -5605. and lower looic of 9234.9 value for the Location model, making the Location model superior over Scale and Shape models. Overall, the Location model performs better than the Scale and Shape models due to higher elpd_loo and lower looic values. Additionally, since the pareto estimates for all the models and across the stations are less than 0.7, that is an indication that the models are reliable and stable.

3.3.3. Metrics Evaluation

Table 15 presents the Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE) evaluation metrics used to assess the forecasting performance of the Location, Scale, and Shape models, for both in-sample and out-of-sample. The model with the smallest RMSE, MAE, and MAPE is considered the best model.

The in-sample performance shows that the location model performs better than the Scale and Shape models across all metrics. These results show that the Location model gives the most precise representation of extremes (better predictive accuracy), with lower RMSE, MAE, and MAPE values, while Scale and Shape models show weaker performances across both stations.

Similar to the in-sample performance, the Location model outperforms the Scale and Shape model across all metrics and stations in the out-of-sample performance. Overall, the Location model is superior over the Scale and Shape model.

3.3.4. Posterior Plots

Posterior predctive plots are a way of checking whether the Bayesian models can reproduce the features of the observed data. In the plots, the dark blue lines (y) represents the posterior distribution of the maximum wind speed, while the lighter blue lines ($y_{\mathrm{rep}}$) represents the distribution of 10 random draws.

Posterior predictive checks for Kimberley shows weaker fit while Upington shows good fit with Location captured well. Upington seems to be more reliable compared to Kimberley.

Figure 5. Kimberley Posterior Predictive checks.

Figure 5. Kimberley Posterior Predictive checks.

3.3.5. Return Levels

Table 16 displays the predicted return levels of daily maximum wind speed across selected return periods and their corresponding 95% confidence intervals indicated in parentheses.

The Bayesian method provides slightly higher estimates, the results in the table shows the return levels increase as the periods increase. Kimberley has the lowest 10 year return level of 21.44 km/h with a narrow confidence interval ranging from 21.18 km/h to 21.72 km/h and Upington has 30.63 km/h return level with a confidence interval ranging from 30.12 km/h to 31.19 km/h. This suggests that a maximum wind speed of 21.44 km/h is expected to increase atleast once in 10 years. 39.85 km/h maximum wind speeds are expected to increase atleast once in 100 years in Upington and 26.52 km/h in Kimberley. From the table, it is clear that Upinton will experience more increase in the wind speed along the years as compared to Kimberly.

The Figure 7, below depicts the return levels plots for the stations.

3.4. Frequentist and Bayesian Performance Comparison

The tables below, Table 17 and Table 18, present the comparison of frequentist and Bayesian metrics (in parantheses). The in-sample performance comparison between the frequentist and Bayesian models show that both approaches deliver similar results with slight differences in the metrics and stations. The Bayesian metrics show lower values compared to the frequentist. For RMSE and MAE, the Bayesian models achieved lower values across all models and stations both for the in-sample and out-of-sample performances. However, for the MAPE, the frequetist achieved lower values in the in-sample performance compared to the out-of-sample.

The out-of-sample performance across all metrics show the Bayesian models to be more powerful than the frequestist. Overall, the results of this paper shows the Bayesian method performs better than the frequentist method for the data used to train the models and also to better to test the models.

Figure 6. Upington Posterior Predictive checks.

Figure 6. Upington Posterior Predictive checks.

Figure 7. Return Level plots.

Figure 7. Return Level plots.

4. Conclusion

Extreme wind speed events have significant impacts on infrastructure, human livelihoods, the economy and agriculture in South Africa. This paper presented a statistical analysis of daily maximum wind speed in a semi-arid region using both frequentist and Bayesian Generalised Additive Extreme Value (GAEV) models. The data were first tested for stationarity using Augmented Dickey-Fuller (ADF), Phillips-Perron (PP) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests, all of which revealed that the series was stationary, which confirmed that the raw series were suitable for extreme value analysis. Although the raw series showed stationary, the inclusion of covariates (seasonality, maximum temperature, maximum relative humidity, and maximum cloud cover) in the GAEV framework explicitly allowed for non-stationary behaviour in the extremes, therefore, capturing climatological drivers of variability. The normality was then assessed using the SW and JB tests where the results revealed that the data were not normally distributed.

Model selection was evaluated using the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Leave-One-Out Cross-Validation (LOO-CV). For forecasting, Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE) were applied to check the performance of the models. Across the considered model specifications, the Location model provided the best overall fit and predictive performance according to these criteria, outperforming the Scale and Shape models.

For the return level estimates, the Bayesian framework produced narrower credible intervals, which account for the uncertainties in the return estimates. This suggest that the return estimates are reliable for decision making. The Location model incorporated day of year (seasonality), maximum temperature, maximum relative humidity, and maximum cloud cover as covariates. The results showed that seasonal cycles and temperature are the most influential drivers of extreme wind speed behaviour, while humidity and cloud cover had less significant effects. These results show the significance of seasonal and thermal dynamics in influencing wind extremes and offer useful climatological insights for semi-arid regions.

Overall, the findings demonstrated that the Bayesian GAEV model outperformed the frequentist GAEV model, particularly in terms of predictive accuracy. Therefore, the study concludes that the Bayesian GAEV framework provides a more effective approach to modelling and forecasting daily maximum wind speed extremes at the Kimberley and Upington stations. Despite the findings of this paper, few limitations were encountered, $\left(1\right)$ using secondary historical weather data from 2019–2024, which provided short observations and limited the ability to capture long-term climate variability and rare extreme wind speed events, $\left(2\right)$ there is limited literature on the modelling of extreme events using frequentist and Bayesian GAEV approaches, making direct comparisons with previous studies challenging.

Future research should $\left(1\right)$ extend the study period beyond 2019–2024 and include additional weather stations in all South African provinces to better capture long-term variability and rare extreme events, $\left(2\right)$ incorporate additional atmospheric and land-surface variables to improve model performance. Finally, integrating the Generalised Pareto Distribution (GPD) within a Bayesian additive modelling framework using threshold exceedance methods would provide further insight into high-frequency extreme wind speed events. Furthermore, future works should consider probabilistic scoring rules or tail-focused measures such as the quantile score for the 95th percentile, to evaluate tail performance.