Bayesian bootstrap quantile regression for probabilistic photovoltaic power forecasting

Photovoltaic (PV) systems are widely spread across MV and LV distribution systems and the penetration of PV generation is solidly growing. Because of the uncertain nature of the solar energy resource, PV power forecasting models are crucial in any energy management system for smart distribution networks. Although point forecasts can suit many scopes, probabilistic forecasts add further flexibility to an energy management system and are recommended to enable a wider range of decision making and optimization strategies. This paper proposes methodology towards probabilistic PV power forecasting based on a Bayesian bootstrap quantile regression model, in which a Bayesian bootstrap is applied to estimate the parameters of a quantile regression model. A novel procedure is presented to optimize the extraction of the predictive quantiles from the bootstrapped estimation of the related coefficients, raising the predictive ability of the final forecasts. Numerical experiments based on actual data quantify an enhancement of the performance of up to 2.2% when compared to relevant benchmarks.


Introduction
The recent growth in distributed Photovoltaic (PV) power generation systems fosters the exploitation of renewable energy and adds further flexibility to electrical distribution systems, which now experience a significant amount of generation in close proximity to load, with obvious advantages in terms of reduced line congestion and losses. PV power is however nondeterministic, as it strictly depends on weather and environmental conditions [1]. This brings significant challenges for transmission and distribution system operators and for market operators. They continually have to deal with local production that may not respond coherently to day-ahead prior generation programs [2]. The impact of such uncertainty on the network operation is therefore severe. For example, large margins of energy reserve are necessary to be able to deal with the deviations from prior generation programs, and the distributed energy resources of smart grids must be managed under such uncertainty [3]. In order to optimally accomplish such tasks, forecasts of PV power for a specific time horizon must be available up to few days before the actual energy generation.
The relevant literature on PV power forecasting is quite heterogeneous [4]. Most researchers and practitioners use deterministic PV power forecasting, i.e., the disposal of a single value that predicts the actual PV power at a specific time horizon [5]. However, PV power forecasting is best addressed within probabilistic frameworks, because of the intrinsic randomness of the physical phenomena. Literature reviews and competition surveys [5][6][7][8] indicate a varied state-of-the-art, although the number of contributions devoted to probabilistic PV power forecasting is much smaller than that in the deterministic framework. Considering that a single spot value can always be extracted from probabilistic forecasts while the opposite is clearly unfeasible, research and contributions in the latter are highly recommended.
Probabilistic PV power forecasting systems range from pure statistical models to hybrid physical-statistical models. High-performance solutions are based on Quantile Regression (QR) models [9][10][11], machine learning approaches (such as gradient boosting [12], quantile regression forests [10,13,14] and k-nearest neighbors [15]). It is worth noting that, although the analytic formulation of QR models is much simpler than machine learning approaches, QR predictions still are somehow competitive in most cases.
Generally speaking, the integration of Numeric Weather Predictions (NWPs) into PV power forecasting models is quite mandatory to improve predictions, particularly for large time horizons [6]. Usually several weather variables are available to forecasters, and thus model selection (i.e., the selection of the most informative predictors for the final model) is typically applied to discard uninformative inputs. To further maximize the exploitation of the available input data, ensemble approaches [7] have been applied with success in probabilistic energy forecasting [10,[16][17][18][19]. These approaches are commonly divided into boosting, stacking, and bagging. Boosting consists of building a "strong" model by combining several weaker ones, which are trained iteratively. Stacking consists of combining the outcomes of several models which are solved in parallel or in cascade, in order to form the final prediction. Bagging consists of building the final prediction as a combination of the output of the same model trained multiple times, using input data resampled with replacement (i.e., bootstrap aggregated) [7].
This paper provides a contribution to probabilistic ensemble PV power forecasting within the bagging framework, based on the interaction between a QR model and a Bayesian bootstrap. A Bayesian bootstrap is the Bayesian analogue for bootstrapping, originally presented in [20]. Bayesian models have been applied with success in probabilistic energy forecasting [16,[21][22][23], although applications of Bayesian bootstrap are still very rare. Thus further contributions are worthy of attention.
Like other bootstrapping techniques, the Bayesian bootstrap can improve the probabilistic forecasts by using resampled data with replacement, which allows for differentiating the output predictions. In the forecasting system presented in this paper, the Bayesian bootstrap works analytically to find the posterior distributions of the QR model parameters, thus differentiating itself from the traditional bootstrap which instead relies on random picks among the available input data. In this paper, a procedure to extract an optimal point from the posterior distributions of the QR model parameters is specifically developed, and this procedure is added to the forecasting system, in order to generate the final PV power forecasts for the target forecast horizon.
The proposal is validated using actual data collected at a real PV installation in Switzerland and actual data published in the framework of the Global Energy Forecasting Competition 2014 [8]. Extensive numerical experiments, based on these data and NWPs provided by an external source [24], are presented in this paper. Several related benchmarks are also presented to validate the accuracy of the forecasts under comparative analyses.
This paper is organized as follows. Section 2 presents the forecasting system and the analytic models used in each stage of the system. Section 3 presents the benchmarks and the error indices used to validate the proposal. Section 4 present the numerical experiments, and the paper is concluded in Section 5.

Methods
The proposed PV power forecasting system based on Bayesian Bootstrap Quantile Regression (BBQR) is illustrated in Fig. 1. The inputs of the system are NWPs NW and historical measured PV power data P. The proposed forecasting system consists of three stages.
The first stage is model selection, i.e., the selection of the most informative predictors among the available pool of predictors. This is performed by evaluating the performance of multiple QR models having different combinations of predictors, and by picking the model which returns the smallest error.
Inputs are therefore pooled in order to form predictor data X (i.e., independent variables in the QR model) which is informative for the PV power (i.e., the dependent variable in the QR model). Data are then partitioned into training datasets P tra and X tra (a 1 × T vector and a T × M matrix, respectively), and into validation datasets P val and X val (a 1 × V vector and a V × M matrix, respectively).
Training data thus contains T occurrences which are used only to train models, whereas validation data contains V occurrences for model selection to develop and refine the forecasting system. M is the number of predictors contemplated in the generic QR model, which therefore has M + 1 parameters.
Multiple QR models are trained, and predictions for the validation period are issued with each model. Since predictions are given in terms of predictive quantiles, the Pinball Score (PS) [25] is considered in order to select the best model. In particular, the QR model returning the smallest PS for the validation period is considered as the most skilled, and is selected as the underlying QR model for the remainder of the system. For notation, the underlying QR model selected in this first stage of the system has M * predictors and M * + 1 parameters.
The second stage consists of applying Bayesian bootstrapping over the selected underlying QR model, in order to estimate the posterior distribution of the parameters of the QR model. Specifically, the Bayesian bootstrap returns R samples extracted from each of the M * + 1 posterior distributions of the M * + 1 parameters of the QR model. As will be shown later, these samples are extracted from a multivariate Dirichlet distribution. A Monte Carlo sampling method then extracts R samples (P hα q i h ) of predictive α q -quantiles of PV power for the target horizon h. The models and the stages of the forecasting system are discussed in the following sub-Sections.

Quantile regression modeling
A QR model links the target variable (i.e., the predictive α q -quantileP hα q i h of PV power at the target time horizon h) to predictors x h ¼ fx 1 h ; …; x M h g related to the time horizon h but available at the forecast origin h − k. The forecast lead time is indicated by k. This paper focuses on day-ahead forecasting with hourly time resolution, assuming that forecasts are issued at midnight of day D-1 for the entire day D (i.e., k = 1, …, 24), although the proposal can also be applied to other short-term PV power forecasting frameworks. Note that theoretically 24 models should be developed (i.e., one for each hour of the day), but since the PV power production is deterministically zero during the night, only the 16 models corresponding to lead times k = 5, …, 20 are considered. In order to lighten the notation in the entire paper, we will not make reference to the forecast lead time k in the symbols although the methodological section is related to a specific lead time k.
The link imposed by the generic QR model for the PV power α q -quantile is: Note that (1) is linear with the parameters, although some predictors can be obtained as multiplicative terms between two or more variables (this allows the introduction of interaction effects among variables [26]).
Parameters are estimated in the training step by minimizing an error score on known data (i.e., supervised training). The PS fits this purpose well, since it can be applied directly on predictive quantiles and, for this reason, it is applied in this paper in order to evaluate the accuracy of PV power forecasts. The minimization problem is: where PSðP tra ;P hα q i tra Þ is the PS of the T forecastsP hα q i tra issued for the training period of length T, calculated with respect to the actual occurrences of PV power in the training set P tra ¼ fP t 1 ; …; P t T g.
Although it is not directly explained in (2), the fore-castsP hα q i tra are obtained from (1), and thus they are functions of β hα q i and are dependent on the T × M matrix X tra which contains the corresponding predictors for the training period. It is:P and therefore (2) can be rewritten in compact form as: where G(•) is a function obtained by combining (2) and (3).

First stage: model selection
In the first stage of the proposed forecasting system, the optimal model is selected among a pool of candidates, which differ in the predictors used to generate the predictions. In this paper, NWPs and 1-day lagged PV power, together with their coupled interactions, form the pool of candidate predictors. The considered NWPs are: total cloud coverage, clear-sky irradiance, total irradiance, total precipitation, pressure and air temperature [24]. Two hypotheses are added to reduce the search dimension for the optimal model: i) if a coupled interaction is a predictor of the model, the two individual variables are forced to occur in the model; ii) only models containing NWPs of total cloud coverage, clear-sky irradiance and total irradiance are considered because of their recognized importance in PV power forecasting. The underlying QR model selected under these hypotheses is the one which minimizes the PS across the Q quantile coverages, i.e., the same optimal combination of M * predictors is selected for the Q considered quantile coverages. To avoid overfitting, the minimum PS is evaluated on the validation dataset P val ¼ fP v 1 ; …; P v V g , which is not used for training the model.

Second stage: Bayesian bootstrap quantile regression
Like traditional bootstrapping techniques, the Bayesian bootstrap can improve the probabilistic forecasts by using resampled data with replacement, which allows for differentiating the output predictions. The Bayesian bootstrap is specifically applied on the underlying QR model selected in the previous stage, in order to evaluate the posterior distribution of the M * + 1 parameters. As shown in the remainder of this subsection, BBQR consists of extracting weights from a Dirichlet distribution R times (once for each bootstrap replicate), building R multinomial distributions using the occurrences and the weights, sampling with replacement from these R distributions and calculatingβ hα q i from (2)-(4) on the bootstrapped data. Therefore, the posterior distribution of the QR model parameters is given by R samplesβ  [27,28] is provided.
The T × (M * + 1) occurrence matrix Y tra ¼½P 0 tra X tra is initially obtained from the transposed training set P 0 tra and the corresponding T × M * matrix X tra . The n th row vector y t n ¼ fP t n ; x 1 tn ; …; x M Ã tn g , taken from the occurrence matrix Y tra , contains the target variable and the predictors at the time step t n . It may be viewed as an item coming from some generic, unknown multinomial distribution F(y), with T available realizations (i.e., past occurrences) y t 1 ; …; y t T .
As shown earlier, the estimated parametersβ hα q i of the QR model come from (4), and therefore they can be viewed as a function of G[F(y)]: In bootstrap (either traditional or Bayesian [20,27,28]), the unknown distribution F(y) is searched for among distributions of the type F T (y): where δ y tn is a degenerate probability measure for the n th vector y t n of occurrences, and ω t n is an assigned weight. For consistency, the weights must satisfy the following conditions: In a traditional bootstrap, the function G[F(y)] is esti- With reference to the generic r th replicate, the functional G½ F hri T ðyÞ is calculated using the weights ω hri ¼ fω hri t 1 ; … ; ω hri t T g , that are obtained by a random extraction from the multinomial distribution: and normalizing by T.
The Bayesian bootstrap differs from the traditional bootstrap since the bootstrapped weights ω hri are not obtained by random extraction from distribution (8). Instead, the vector ω ¼ fω t 1 ; …; ω t T g is the object of Bayesian analysis, which aims at evaluating a posterior distribution p(ω| Y tra ) of this vector of weights, given occurrence data Y tra . A prior distribution p(ω) should be imposed upon the parameters ω to start the Bayesian inference [20,27,28]. A convenient choice is to select a Dirichlet distribution, which is a conjugate prior for the multinomial distribution of y [20,27,28]. In such a case, the posterior distribution p(ω| Y tra ) is itself a Dirichlet distribution, f Dir (1, …, 1; 1, …, 1) [20,27]. This allows applying a Monte Carlo sampling method to get the Bayesian bootstrapped samplesβ

Third stage: extraction of a single value from the Bayesian bootstrapped predictive PV power
The R samplesP hα q i h of the predictive α q -quantile of PV power can be interpreted as probabilistic predictions for the predictive quantile. Sample quantiles and confidence intervals of the predictive α q -quantile of PV power can therefore be estimated fromP Since it will be of use later in the paper, the generic sample τ q -quantile estimated fromP Probabilistic PV power forecasts are usually given in terms of predictive distribution or a set of predictive quantiles at different coverage levels, and the redundancy given by multiple samples for each quantile level can lead to misinterpretation of the results in practical utilization of forecasts. A dedicated procedure is developed in this paper to reduce the redundancy of the forecasts by extracting a single value from the samples of the predictive α q -quantile of PV power. This single value is treated as the final predictive α q -quantileP fromP hα q i h as: where the specific coverage τ Ã q of this sample quantile is the object of an optimization problem aimed at minimizing the PS of the final forecasts calculated on the validation set P val , i.e.: Note that this procedure is made independent for each quantile coverage α 1 , …, α Q , for simplicity. Therefore, possible quantile crossing in the final PV power forecasts is corrected by post-processing the results with simple sorting in ascending order across the Q coverages.

Benchmarks and error indices
This section presents the benchmarks and the error indices used to assess the validity of the proposal and to compare the probabilistic forecasts.

Benchmarks
Five benchmarks are considered to provide a fair comparison of the results. They are listed below.
Simple QR (SQR): the first benchmark [10,16] is introduced to be used as a reference in which each predictive quantile is directly provided as a single value, rather than by passing through the bootstrap. This allows assessing whether the bootstrap is effective or not in improving forecasts. To provide a fair comparison, the same model selection procedure presented in Section 2.2 within the framework of the proposed forecasting system based on BBQR is also adopted for SQR. Traditional Bootstrap QR (TBQR): the second benchmark [29] is introduced in order to evaluate if the Bayesian bootstrap is more effective than the traditional bootstrap in increasing the skill of the final forecasts. For fair comparison, the same QR model selection procedure, the same Monte Carlo sampling method and the same procedure to extract a single value from the R samples, presented respectively in Sections 2.2, 2.3 and 2.4, are applied to TBQR. Quantile Regression Neural Network (QRNN) and Gradient Boosting Regression Tree (GBRT): the third and fourth benchmarks of QRNN [30] and GBRT [31] are introduced to provide independent references that do not come from QR-based models.
QRNN is formulated to simultaneously predict several PV power quantiles to reduce the quantile crossing effect [30]. The basic neural network architecture used to develop QRNN is the multilayer perceptron with a single hidden layer. Hyperparameter optimization is performed through a validation procedure on the validation dataset P val , to maintain statistical fairness with the other models. QRNN is implemented using the qrnn package in R [32].
GBRT is developed individually for each considered quantile, and the predictions are post-processed in a sorting procedure to avoid quantile crossing. Also in this case, the hyperparameter optimization is performed through a validation procedure on the validation dataset P val , to maintain statistical fairness with the other models. GBRT is implemented using the gbm package in R [33].
Seasonal Persistence Model (SPM): the fifth benchmark is based on the underlying daily periodicity of the PV power pattern driven by the rotation of the Earth around its own axis [16]. In practice, each predictive quantile is the PV power observed the day before: This benchmark is added in order to provide a naive, unbiased reference for comparison.

Error indices
Two error indices are used to quantify the accuracy of the forecasts. The first index is the abovementioned PS, which is a strictly proper score [25] that simultaneously addresses the reliability and the sharpness of forecasts [16,25]. Its formulation is: where the indicator function IðP h ;P hα q i h Þ is: A comprehensive PS can be obtained averaging across multiple forecast issues (e.g., the V issues in the validation set) and summing over the Q quantiles. PS is negatively oriented, so a smaller PS indicates better forecasts. The Normalized PS (NPS) is provided in our numerical experiments to get scaleindependent results. It is: where P rated is the rated power of the PV system. The second error index is the Average Absolute Coverage Error (AACE), and it addresses the reliability of the forecasts, i.e., the correspondence between the estimate and the nominal coverages of the predictive quantiles [34]. Because of its intrinsic properties, it can only be formulated for multiple forecast issues. For sake of clarity it is referred in (16) to the validation set (although it can be easily adapted to other data sets). In such a case, the estimated α qcoverageα q is: and the Absolute Coverage Error (ACE) on the nominal α q -quantile is: The AACE across the Q coverages can be obtained as a percentage value as: The AACE is negatively oriented, so a smaller AACE indicates more reliable forecasts.

Results and discussion
The proposed forecasting system is assessed using two actual PV power datasets.
Dataset #1 consists of data collected at a PV installation within the ReIne smart grid laboratory in Switzerland [35]. Data span February 1, 2016 to November 30, 2018, and are initially pre-processed to correct potential outliers, missing and bad data, and to average the measurements in order to obtain hourly PV power data.
Dataset #2 consists of zone-1 data published in the framework of the Global Energy Forecasting Competition 2014 [8] at an unspecified location in Australia, and data span April 1, 2012 to June 30, 2014. Since these data were already pre-processed by the competition organizers and arranged with an hourly resolution, no further processing is applied on dataset #2.
NWPs are provided in both cases by the European Centre for Medium-range Weather Forecast (ECMWF) [24] for the locations at which the PV systems are installed, and for the same time intervals related to datasets #1 and #2. All the weather forecasts refer to the midnight run, i.e., NWPs are available at 24:00 of day D-1 for the entire day D, in order to suit the desired forecasting framework [24].
Data are normalized in the 0-1 range to accommodate for the very different intervals spanned by the variables. The normalized valuez h of the occurrence of the generic variable z at time h is: where z min and z max are the minimum and maximum values registered for the variable, respectively. The available datasets are partitioned in three subsets. Dataset #1 is split into a training set spanning February 1, 2016 to December 31, 2017 (i.e., T = 16800), a validation set covering first 6 months of 2018 (i.e., V = 4344), and a test set only used to assess the accuracy of forecasts, covering the remaining 5 months of 2018 (i.e., 3672 forecast issues). Similarly, dataset #2 is split into a training set from April 1, 2012 to October 31, 2013 (i.e., T = 13896), a validation set covering the following 5 months (i.e., V = 3624), and a test set covering the remaining 3 months of 2014 (i.e., 2184 forecast issues).
The number R of bootstrapped is searched for in the range 1000-10,000, considering the size of the two datasets. Four tests with R = 1000, 2000, 5000, 10000 were run, and the performances in these four cases were checked on the validation datasets P val . A good compromise was found by selecting R = 5000 for both datasets. 5000 is the value used eventually to predict the PV power in the test periods.

Results for dataset #1
The outcome of the first stage of the proposed forecasting system determines the model selected for the specific PV system. For dataset #1, the selected model is: where tcc h , csi h and ti h are NWPs of total cloud coverage, clear-sky irradiance and total irradiance, respectively. The number of predictors of the selected model is therefore M * = 7. The forecast results for the test set of dataset #1 are shown in Table 1 via NPS (summed across the Q = 19 quantiles and averaged through the test set) and AACE. BBQR returns a NPS smaller than SQR, TBQR, QRNN, GBRT and SPM benchmarks by 2.2%, 0.6%, 5.4%, 1.4% and 51.0%, respectively. Bootstrapping increases the accuracy of forecasts, since both the bootstrapped methods (BBQR and TBQR) outperform SQR, although the Bayesian-based procedure slightly outperforms the traditional bootstrapping procedure in terms of NPS.
Further details on the skill of forecasts can be evaluated from the results of the experiments. Figure 2 shows the NPS, averaged through the test set of BBQR, SQR and TBQR forecasts for each nominal quantile level. Figure 3 shows the NPS (summed across the Q = 19 quantiles and averaged through the test set) of BBQR, SQR and TBQR forecasts versus the forecast lead time. The similar patterns illustrated in these two figures are determined by the same underlying QR model used in all three forecasting methods. It can be determined that peak NPS occurs for middle coverage levels, whereas the NPS changes with the lead time for two reasons: i) forecasts inevitably tend to lose accuracy as lead time increases, and ii) the "bell-shaped" PV power patterns have small errors in proximity to dawn and dusk. BBQR forecasts are also the most reliable, as the AACE is reduced by 59%, 6.7%, 56.3% and 61.2%, with respect to SQR, TBQR, QRNN and GBRT, respectively. SPM AACE is not presented, since SPM forecasts are the same for each quantile coverage. To compare the probabilistic QR-based forecasts in detail, Fig. 4 shows the reliability diagrams, outlining the estimated coverages versus nominal, of BBQR, SQR and TBQR. Both BBQR and TBQR show similar patterns, with slightly overestimated coverages in the range 0.5 to 0.8. However, the SQR coverages are overestimate for all nominal levels.
In order to provide a graphical interpretation of the PV power forecasts versus time, Fig. 5 shows the BBQR prediction intervals for 1 week of the test period. Prediction intervals are given for rates 90%, 50% and 10%, and they are plotted together with the actual PV power.

Results for dataset #2
The results for dataset #2 are presented in a more compact form, to avoid duplication. Table 2 shows the NPS and AACE for BBQR and the benchmarks. As seen, BBQR returns an NPS that is smaller than the benchmarks for this dataset, too. The respective reductions are around 4.7%, 2.4%, 6.3%, 2.0% and 53.4%, compared to SQR, TBQR, QRNN, GBRT and SPM. Bootstrapping is again proved to increase the accuracy of forecasts since both bootstrapped methods (BBQR and TBQR) outperform the SQR. The Bayesian-based procedure slightly outperforms the traditional bootstrapping procedure in terms of NPS.
The AACE of BBQR forecasts is again the smallest, accounting for reductions of 41.3%, 7.7%, 30.4% and 34.9% with respect to SQR, TBQR, QRNN and GBRT, respectively. Figure 6 shows the reliability diagrams of the probabilistic QR-based forecasts. BBQR and TBQR show similar patterns, with slightly overestimated coverages in the range 0.5 to 0.8, whereas the SQR coverages are underestimated for all nominal levels.

Discussion
The results obtained in the experiments denote the ability of BBQR to slightly increase the performance of the probabilistic predictions. Compared to traditional bootstrap approaches, the NPS is reduced in the range from 0.6% to 2.4%, and the overall reliability is slightly increased. The proposed method also performs well when compared with the state-of-theart consolidated probabilistic models in energy forecasting, such as QRNN and GBRT.
These promising results indicate the applicability of Bayesian bootstrap techniques for estimating the parameters of different models, thus not limiting the  analysis to QR-based models, in order to consolidate the technique in probabilistic energy forecasting. Some limitations apply to the type of prior and posterior distributions of the parameter. As shown in the methodology section, although conjugate priors ease the process of Bayesian bootstrap sampling, numerical methods (e.g., Metropolis-Hastings or Gibbs sampling) can be applied to draw samples from the posterior distributions of parameters even if the prior is not conjugate, thus allowing for generalizing the approach under different assumptions.

Conclusions
This paper provides a new insight on bootstrapping in energy forecasting. A novel forecasting system based on Bayesian bootstrap is applied to day-ahead prediction of PV power. A Bayesian bootstrap is used to provide information regarding the parameters of an underlying QR model, and a Monte Carlo sampling method is specifically developed in order to generate the final forecasts. The proposed forecasting system is compared to several benchmarks, to assess the efficacy of the Bayesian bootstrap with respect to traditional bootstrap and others not using a bootstrapping technique. Numerical experiments based on actual PV power data and on NWPs confirm that BBQR is a valid approach with slightly increased accuracy and reliability of PV power forecasts.
This paper opens new research perspectives in further applications of Bayesian approaches and the Bayesian bootstrap in probabilistic energy forecasting, and in the development of techniques to extract final forecasts from bootstrapped results. : Generic m th element of the vector of predictors related to the target horizon h; X tra : Predictor data of the training period; X val : Predictor data of the validation period; y tn : Generic n th row vector taken from the occurrence matrix Y tra ; Y tra : Occurrence matrix, containing PV power data and predictor data of the training period; z h : Variable z at the target horizon h;z h : Normalized variable z at the target horizon h; z min : Minimum value registered for variable z; z max : Maximum value registered for variable z; α q : Generic q th quantile coverage level;α q : Generic q th estimated coverage; β hαqi : Vector of parameters of the QR model for the α q -quantile; β hαqi m : Generic m th parameter of the QR model for the α q -quantile;β hαqi : Vector of estimated parameters of the QR model for the α q -quantile;β hαqi m : Bootstrapped samples of the generic m th parameter of the QR model for the α q -quantile;β hαqi m : Generic m th estimated parameter of the QR model for the α q -quantile; δ y tn : Degenerate probability measure for the generic n th vector y tn of occurrences; τ Ã q : Optimized quantile coverage level of the sample quantile estimated from the baggedP hαqi h ; ω: Vector of weights assigned to occurrences; ω tn : Generic n th weight assigned to the vector y tn of occurrences; ω 〈r〉 : Vector of weights assigned to the occurrences in the r th bootstrap replicate; ω hri tn : Generic n th weight assigned to the vector y tn of occurrences in the r th bootstrap replicate