Estimation of Extreme Quantiles of Global Horizontal Irradiance: A Comparative Analysis Using an Extremal Mixture Model and a Generalised Additive Extreme Value Model

Solar power poses challenges to the management of grid energy due to its intermittency. To have an optimal integration of solar power on the electricity grid it is important to have accurate forecasts. This study discusses the comparative analysis of semi-parametric extremal mixture (SPEM), generalised additive extreme value (GAEV) or quantile regression via asymmetric Laplace distribution (QR-ALD), additive quantile regression (AQR-1), additive quantile regression with temperature variable (AQR-2) and penalised cubic regression smoothing spline (benchmark) models for probabilistic forecasting of hourly global horizontal irradiance (GHI) at extremely high quantiles (τ = 0.95, 0.97, 0.99, 0.999 and 0.9999). The data used are from the University of Venda radiometric in South Africa and are from the period 1 January 2020 to 31 December 2020. Empirical results from the study showed that the AQR-2 is the best fitting model and gives the most accurate prediction of quantiles at τ = 0.95, 0.97, 0.99 and 0.999, while at 0.9999-quantile the GAEV model has the most accurate predictions. Based on these results it is recommended that the AQR-2 and GAEV models be used for predicting extremely high quantiles of hourly GHI in South Africa. The predictions from this study are valuable to power utility decision-makers and system operators when making highrisk decisions and regulatory frameworks that require high-security levels. This is the first application to conduct a comparative analysis of the proposed models using South African solar irradiance data, to the best of our knowledge.


Background
The generation of power from clean energy sources makes an important contribution to sustainable development. Energy policymakers are establishing more penetration of clean energy onto the grid, such as solar and wind power [1]. Variability of accurate renewable energy forecasts and lack of consistency on renewable power generation are the limitations for using renewable energy sources (RES) to produce power for supplying peak energy demand [1]. Accurate forecasting models are useful to system operators and market players in solving problems created by the integration of RES and also viable to handle their variability and uncertainty [2]. However, this study focuses on modelling and forecasting global horizontal irradiance (GHI) at extremely high quantiles. Solar power is one of the most preferable RES to be integrated onto the grid because its contract price is cheaper since it depends on primary sources [3]. The energy sector needs to take the direction of exploring solar energy as it is one of the cleanest energy sources [4]. Accurate solar power forecasts are useful for an economic operation dispatch, optimal unit commitment and ensuring the stability of a national grid. It also reduces the uncertainties of solar energy sources and results in assuring safety and easier grid management. Lack of accurate forecasts of solar power generation creates problems for grid energy management [3,5].
Solar irradiance forecasting is required to guide the grid operators to apply the relevant procedure to optimise electricity production and reduce production costs. Solar irradiance forecasts are used in an electrical grid for the integration of photovoltaic (PV) power [5]. Power system operators can have problems because of inaccurate extreme quantiles of GHI forecasts. Accurate extreme forecasts of GHI quantiles help decision-makers in making high-risk decisions and regulatory frameworks that require high-security levels. However, accurate forecasts are important for the effective operation of the electrical grid when solar energy is integrated. It has been shown that solar power forecasting is essential in different areas of planning operations in the energy sector such as unit dispatch and renders schedules of production of power from renewable energy sources for the next hours or days. These forecasts are also used to know in advance the amount of solar power that will be integrated into the grid in the following hours or days. The purpose of producing new forecasting models for solar power is to provide accurate forecasts [6].

A review of the solar irradiance forecasting literature
Accurate forecasts are very important in the energy utility because operators rely on them to maintain the uncertainty created when the penetration of solar energy increases [5]. According to [7], accurate PV power forecasting is required for electric power system operation during the installation of grid-connected PV generation on a huge scale. Application of PV power forecasting is not only used in decision making for grid dispatch, it is also used to complement the control of other multiple power sources. Potential statistical modelling errors including changing weather patterns pause challenges in forecasting accurately solar power [8]. To reduce the intermittent solar power and operational expenses when integrating it on the electricity grid it is important to have accurate forecasts and also capture the uncertainty surrounding the forecasts [9].
Reliable accurate forecasts are required to maintain the balance of fluctuating PV energy generation and also make the integration to render better quality service [10]. Solar power forecasts are mainly classified periodically, ranging from very short-term, short-term, medium-term and long-term [11]. According to [11], very short-term forecasts consist of 15 minutes to a few hours ahead with a granularity of 30 seconds to 5 minutes, which are used for events ramping and variability related to operations. For short-term forecasts, twenty-four up to 168 hours ahead with a granularity of hourly. These forecasts are useful for transmission scheduling including unit commitment. Medium-term forecasts which are important for planning and asset optimisation range from a week to a year ahead. Long-term forecasts are generated from a year up to several years ahead and are used for planning of the extension of the electricity grid including the building of new power plants [11].
Many researchers have been focusing on short-term forecasting with an hourly cycle [7]. However, the focus of this study is on hourly GHI forecasting at extremely high quantiles. According to references [12] and [13], solar irradiance short-term forecasts are essential for operational planning including scheduling of power systems. PV power forecasting consists of models such as physical, statistical, machine learning, and hybrid integration models [7]. The forecasting techniques of solar irradiance are typically separated into physical and statistical approaches including machine learning algorithms [10,14]. Physical approaches use numerical weather predictions (NWP) and statistical approaches use historical time solar time series data [3,4,10,14]. There are three types of solar irradiance which are direct normal irradiance (DNI), diffuse horizontal irradiance (DHI) and global horizontal irradiance (GHI), respectively [12]. GHI is the sum of DNI and DHI.
Solar energy plays an important role in electricity markets [5]. Reference [5] compiled a review of solar power forecasting focusing on the current improvements and future projections of solar energy. The idea of the study was to present a precise forecast with its analysis of the economic implications. The study revealed that most of the research has been done focusing on day-ahead forecast horizon [5]. Accurate solar power forecasts are required by utilities and independent system operators, as well as solar power producers and energy traders for core operations [11]. Reference [11] proposed a set of standards for performance evaluation of solar power forecasts. The standards consist of sound methodologies and a lot of field experience, which provide accurate inter-agency comparisons of forecast performance [11]. A hybrid short-term forecasting model was proposed by reference [7] to solve the problem of forecasting accuracy using solar data from Yunnan PV power generating plant. Several forecast evaluation metrics were used which include normalised absolute average error, normalised root-mean-square 118 Estimation of Extreme Quantiles of Global Horizontal Irradiance: A Comparative Analysis Using an Extremal Mixture Model and a Generalised Additive Extreme Value Model error, among others. Results from the study showed that the proposed model was the best fitting model [7].
Solar irradiance forecasting using South African data is discussed in the literature. In a recent study, reference [14] predicted hourly global horizontal irradiance using data from the University of Pretoria radiometric station. A comparative analysis of long short-term memory networks (LSTM), support vector regression (SVR) and feed-forward neural networks (FFNN) models was done with results showing that the FFNN model gives the most accurate forecasts. The forecasts from the individual models were then combined to get an ensemble LSTM-SVR-FFNN model via quantile regression averaging (QRA) and convex combination (CC). The ensemble LSTM-SVR-FFNN model via QRA was the best fitting model with the most accurate forecasts. The study suggested that the sufficient and thorough metrics for evaluation, as well as statistical tests, render more perspective into the proposed forecasting methods [14]. Using data from Tellerie radiometric station in the Northern Cape province of South Africa, reference [3] developed partially linear additive quantile regression models which were then used in forecasting day-ahead GHI. The selection of variables was done using the least absolute shrinkage and selection operator (Lasso) via hierarchical interactions. The forecasts from the individual models were then combined using QRA and CC. The study revealed that QRA yields the most accurate forecasts [3].
Reference [15] used a multiple linear regression model to forecast global solar irradiance in South Africa. Covariates used in the study were weather variables, traditional extraterrestrial irradiance including sunshine hours. The study showed that the solar irradiance models performance and accuracy in other areas improved by the inclusion of weather parameters [15]. Reference [16] solved the issue of long-range dependence ingrained in the solar irradiance data in South Africa using three models which are seasonal autoregressive fractionally integrated moving average (SARFIMA), harmonically Coupled SARIMA and regression model with SARFIMA error terms. In the study, additive quantile regression was used as a benchmark model. Empirical results from the study suggest that long memory is anti-persistent in all the developed models [16].
Reference [2] estimated extremely low (below 0.05-quantiles) and extremely high (above 0.95-quantiles) quantiles of wind power data from a wind power plant located in Galicia, Spain. The study applied conditional extreme value theory estimators by integrating gradient boosting trees and a truncated generalized Pareto distribution. Based on pinball losses, the proposed method outperformed the benchmark models. The study did not include the weather variables as additional covariates to improve the accuracy of the proposed model. However, this study will include temperature variables as a way to improve the precision of models and also use different models for probabilistic forecasting of distribution tails.

Research highlights and contributions
The main contribution of this study is that it compares semi-parametric extremal mixture (SPEM), generalised additive extreme value (GAEV) or quantile regression via asymmetric Laplace distribution (QR-ALD), additive quantile regression (AQR-1), additive quantile regression with temperature variable (AQR-2) and penalised cubic regression smoothing spline (benchmark) models in predicting extremely high hourly GHI data. To the best of our knowledge, this is the first study to do a comparative analysis of the proposed models using data from South Africa. The following is an overview of the study's highlights and findings: • The forecasts from SPEM, GAEV, BM, AQR-1 and AQR-2 models were integrated using simple average (AVG) and median (MED) methods to develop AVG and MED forecasts, • Based on both the continuous rank probability score (CRPS) and pinball loss function (PLF), the AQR-2 model was identified to be the most accurate model for quantiles τ = 0.95, 0.97, 0.99 and 0.999 and the GAEV model was also identified to be the most accurate model for quantile τ = 0.9999, • The results from this study yield improved forecasts for extremely high quantiles of hourly global horizontal irradiance (GHI), • This is the first study to compare the SPEM, GAEV, BM, AQR-1 and AQR-2 models in estimating extremely high quantiles of hourly GHI using data from South Africa, • Providing accurate solar power forecasts which are needed in the integration of solar power in the electricity grid.
The sections of the study are arranged as follows: Section 2 discusses the proposed models. Section 3 presented the empirical results and discussion while conclusion of the study presented Section 4. The proposed modelling framework is based on the work of [17].
The excess distribution function over a sufficiently high threshold µ is then given by [17,18]: where the cumulative distribution function (CDF) F t is assumed to be in the domain of attraction of the Fréchet class of distributions. The CDF F t,µ (x) can be approximated by a Pareto distribution, which is defined by: where the parameter θ > 0 denotes the conditional tail index (extremal index) and an unknown threshold µ ≥ x 0 , which depends on t. The distribution function F t can be approximated by the empirical distribution function. The semi-parametric extremal mixture model (SPEM) is defined by [17,18]: where y t,µ,θ is the global horizontal irradiance (GHI). From equation (3), F t (x t ) is a bulk model, in which we fit a non-parametric model (kernel density) and is the tail model, in which we fit the Pareto distribution. In the study we use the kernel estimator to estimate the parameter θ. A point-wise data driven procedure is used to choose a threshold µ. The bandwidth denoted by h will be determined by minimising the cross-validation function. Let K(.) denote a kernel function which is assumed to be non-negative continuous and symmetric to the real line such that K(x) ≤ 1. The set of weights for i = 1, ..., n defined by the kernel function is given by: where bandwidth parameter h > 0.

Threshold selection
In reference, [19], extremal mixture models for estimating sufficiently high thresholds are proposed to overcome the drawback of fixed threshold models which do not capture threshold uncertainty. In their paper, [19] grouped the threshold methods into four classes which are the classical fixed threshold, re-sampling based, tail fraction estimation and extremal mixture models. Reference [20] developed an extremal mixture model which combines a non-parametric kernel density for the bulk distribution and GDP for the tail distribution. Uncertainty around threshold selection and parameters estimation is addressed by Bayesian inference [21]. In a recent study, reference [22] discussed an automatic threshold selection method. Reference [21] developed a time-varying threshold which is coupled with a positive shift factor. In this present study, the selection of the threshold is done by testing a goodness-of-fit test discussed by reference [18], for the parametric part of the model 3. The k upper statistics are used to test the tail adjustment of Pareto distribution. If the test is not rejected, the number k of upper statistics is then increased. The tail adjustment is tested again until it is rejected.

Inference
The maximisation of the weighted quasi-log-likelihood estimator is defined by [17,18]: wheren t,h,µ = n i=1 W t,h (t i )I {Xti>µ} represents the weighted number of observations above the threshold µ. For any τ ∈ [0.99, 1) the estimator of the τ quantile of X t is given by: andτ µ =F t,h (µ).

Generalised additive extreme value model
The high threshold µ(x) is estimated using quantile regression (QR) based on a pre-specified quantile τ . The likelihood-based inference is achieved by comparing the tilted loss function used in QR and the asymmetric Laplace distribution (ALD). Let Y (x) be a random variable indexed by a covariate x [23,24]. Then: and the density function of ALD is given by [23,24]: where µ(x) corresponds to quantile 0 < τ < 1 at location x, σ > 0 and ρ τ (y) = y(x) τ (x) − I{y(x) < 0} denotes the check function, for indicator function I. Assume generalised additive model (GAM) forms in covariate x defined by [23,24]: where y t is a GHI, β ij denotes basis coefficients and b ij is the basis functions.

Inference: Restricted maximum likelihood
In order to make inference for β, penalised log-likelihood function is defined by [23,24]: where λ = (λ 1 , ..., λ I ) represents the smoothing parameters, S λ denotes a penalty matrix with elements determined by the basis function b ij that has been chosen. And penalty matrix is defined by S λ = I i=1 λ i S i , where S i denotes rows and columns of the matrix corresponds to b i j , i = i and consists of zeros. The parameters of extreme value distribution are estimated by a restricted maximum likelihood defined by [23,24]: where for given λ,β λ miximises (β λ , λ), H(β T λ ) = −∇∇ T (β, λ)| β=β λ and |S λ | + is the product of positive eigenvalues of matrix S λ .

Extremal index
The censored log-likelihood model is used for estimation of an extremal index, which is given by [24]: where w u = −1 log τ .

Additive quantile regression model
An additive quantile regression (AQR) model is one which combines GAM and QR models. The AQR model is defined by [25,26,27]: where y t,τ is a GHI at time t = 1, .., n at quantile τ , s i,τ denote the smooth functions and ξ t,τ is the error term. In this study we are going to consider two scenarios: one AQR model (16) (AQR-1), which has one linear covariate t and two AQR model (17) (AQR-2) which has one linear covariate t and temperature variable. The two models are written as:

Benchmark model for estimating extreme conditional quantiles
We first fit a penalised cubic regression smoothing spline given by [21]: where y t denotes our GHI, λ is a smoothing parameter and ε t is the error term. We then extract residuals ε t = π 1 (x) −π 1 (x) and then estimate the shift factor µ τ ∈ , τ ∈ [0:95; 1). The shift factor has to be sufficiently large to satisfy the asymptotic conditions when we fit the GPD [21].
To estimate the shift factor µ τ , threshold stability plots and extremal mixture models will be used. Table 1 presents a summary of some of the advantages and disadvantages of the proposed models, which are SPEM, GAEV, AQR and BM.

Evaluation of probabilistic forecasting methods and error measures
This section presents the methods for evaluation and comparison of probabilistic forecasts from proposed models of extreme quantiles. In this study two error measures will be used, are the continuous rank probability score (CRPS) and pinball loss function (PLF). When the quantile score is low, the forecasting model is more accurate.

Continuous ranked probability score
Continuous ranked probability score (CRPS) can handle both calibration and sharpness at the same time [28]. The CRPS is given by [2,28,29]: where F is the forecast distribution and QS τ represents the quantile score defined by: where I denotes an indicator function.

Pinball loss function
Even if CRPS is useful for measuring the quality of cumulative distribution function forecast, the study by reference [30] indicated that the CRPS is more applicable to distributions with smaller uncertainty intervals. However, a more suitable scoring rule is the PLF or quantile loss [28]. The PLF is defined by: if yt > qτ , where q τ is a quantile forecast and y t denotes the observed value of hourly GHI.

Combining the extreme quantiles
In this study, simple average and median methods are used to combine the extreme quantiles from the proposed models.

Simple average method
A simple average method performs better than other sophisticated methods when combining forecasts. This method calculates a simple average of the endpoints of the intervals. Simple averages are used in combining point forecasts due to their robustness, simplicity and providing accurate forecasts. The lower and upper averages models are defined by

Median method
The median method is a known method used in summarising data because is less sensitive to extreme values than the mean. The

Empirical results and discussion
This study used GHI data for the period, 1 January 2020 to 31 December 2020, giving a sample size of 4083 observations. The data is obtained from the Southern African Universities Radiometric Network (SAURAN) website (http://www.sauran.net). In this study, the station of focus is the University of Venda radiometric station found in an inland location at Vuwani Science Research Centre, in the Limpopo province of South Africa. It is on latitude -23.13100052, longitude 30.42399979 and elevation 628 m. Figure 1 shows a picture of the pyranometer at the University of Venda radiometric station. The study seeks to predict extremely high quantiles of GHI using R statistical software.

Exploratory data analysis
The summary statistics of hourly GHI are given in Table 2. Since the mean and median of the GHI data are not equal, the distribution is not normally distributed. The non-normality of the distribution of GHI data is confirmed by the skewness and kurtosis values which are 0.5734277 and -0.9057017, respectively. Figure 2 presents the time series plot of GHI together with histogram and box plots. These plots suggest that the GHI data does not follow a Gaussian distribution. A plot of hourly GHI superimposed with a non-linear trend is shown in Figure 3 and this best fit was determined based on "mgcv" developed by reference [32].

Forecasting results
The models considered in this study are the SPEM, GAEV, BM, AQR-1 and AQR-2. The R-packages such as "extremefit" developed by reference [18] is used for the SPEM, "evgam" developed by reference [23] is used for the GAEV model and "qgam" developed by reference [27] is used for AQR models. The quantiles considered in this study are τ = 0.95, 0.97, 0.99, 0.999 and 0.9999. All the data is used for predictions. The plots of actual hourly GHI and forecasts of GHI at quantiles (τ = 0.95, 0.97, 0.99, 0.999 and 0.9999) using SPEM, GAEV, BM, AQR-1 and AQR-2 models are given in Figure 4 -6 respectively. It is shown that the quantile predictions for each model follow the actual hourly GHI data remarkably well at a high quantile. The simple average (AVG) and median (MED) methods are used to combine forecasts for the SPEM, GAEV, BM, AQR-1 and AQR-2 models at each quantile. The plots of actual hourly GHI and forecasts of GHI using AVG and MED methods at quantiles (τ = 0.95, 0.97, 0.99, 0.999 and 0.9999) are given in Figure 7. Table 3 gives a summary of the values of the evaluation metrics (PLF and CRPS) for the proposed models. At quantile levels, i.e τ = 0.95, 0.97, 0.99 and 0.999 AQR-2 is the most accurate model based on the PLF and at τ = 0.9999 GAEV is the most accurate model based on the PLF and CRPS. Figure 8 shows the plots of actual hourly GHI and forecasts of GHI from the best-fitting model AQR-2 at quantiles τ = 0.95, 0.97, 0.99 and 0.999. Figure 9 also shows the plots of actual hourly GHI and forecasts of GHI from the best-fitting model GAEV at quantile τ = 0.9999.

Discussion of results
This study was motivated by previous studies such as [3,14], among others. The focus of the study was on a comparative analysis of SPEM, GAEV, BM, AQR-1 and AQR-2 models in the prediction of extremely high hourly GHI using South African data. The forecasts from SPEM, GAEV, BM, AQR-1 and AQR-2 models were combined using AVG and MED methods to develop AVG and MED forecasts. The PLF and CRPS were used to measure the accuracy of models at each quantile including AVG and MED forecasts. Based on both the PLF and CRPS, the best-fitting model for quantiles τ = 0.95, 0.97, 0.99 and 0.999 was found to be AQR-2, while for the 0.9999-quantile the best model is GAEV. The most accurate model was selected based on the smallest value of PLF and CRPS. It was revealed that fitting a model (i.e AQR-2) including temperature variable improves the forecasts because the temperature is the main predictor of GHI. Prediction of extremely high quantiles of GHI solar power could help system operators to know the possible largest solar power which can be integrated on the national grid. Accurate forecasts of

130
Estimation of Extreme Quantiles of Global Horizontal Irradiance: A Comparative Analysis Using an Extremal Mixture Model and a Generalised Additive Extreme Value Model Figure 9. Plots of actual hourly GHI and forecasts of GHI at high quantiles using GAEV model at τ = 0.9999 GHI produced by this study are useful by decision-makers and system operators in effectively balancing demand and supply of electricity that is environmentally friendly and also secures the future economic prosperity of the country.

Conclusion
This study carries out a comparative analysis of SPEM, GAEV, BM, AQR-1, AQR-2, AVG and MED models in the prediction of extremely high hourly GHI using data from the University of Venda radiometric station. The thrust of the study was on hourly GHI at extremely high quantiles i.e τ = 0.95, 0.97, 0.99, 0.999 and 0.9999. The AQR-2 was discovered to be the most accurate model for quantiles (τ = 0.95, 0.97, 0.99 and 0.999) and the GAEV model was discovered to be the most accurate model for quantile (τ = 0.9999). The best-fitting models were also used to forecast the GHI for each quantile. The findings of this study could be useful to decision-makers from power utilities in the optimal integration onto the electricity grid of the energy generated from solar power plants. In future research inclusion of additional covariates such as cloud cover, wind speed, seasons, among others can improve the predictive abilities of the models. Furthermore, it will also be more interesting to use the proposed models on the prediction of other renewable energy sources, such as wind power generation.