Dennis Van Der Meer

Pierre Pinson

Simon Camal

Georges Kariniotakis

CRPS-based online learning for nonlinear probabilistic forecast combination

Keywords:

published or not. The documents may come

Introduction

The combination of probability distributions issued by experts has a long history that can be traced back to at least [START_REF] Stone | The opinion pool[END_REF]. The linear opinion pool, as labelled by [START_REF] Stone | The opinion pool[END_REF], is the convex combination of component probability distributions. In probabilistic forecasting, forecasters aim to maximize the sharpness of the forecasts, subject to calibration [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF]. Calibration refers to the agreement between the forecasts and observed probabilities; for instance, when a forecaster predicts daily overcast conditions with 80% probability, cloudiness should actually occur on 80 of the 100 days with such conditions. Ensemble forecasts from numerical weather prediction (NWP) models tend to exhibit underdispersion, which implies that the forecasts are overconfident [START_REF] Wilks | Chapter 3 -univariate ensemble postprocessing[END_REF]. Conversely, probabilistic forecasts can be overdispersed, which means that the forecaster is underconfident and issues forecasts with too much variance. These types of miscalibration can negatively affect decision-making based on such forecasts and care needs to be taken to ensure proper calibration. On this point, [START_REF] Hora | Probability judgments for continuous quantities: Linear combinations and calibration[END_REF] notes that "there are theoretical reasons for questioning the use of linear combinations of experts' probabilities. This concern stems from understanding that well-calibrated experts cannot be combined without introducing miscalibration." Nevertheless, linear forecast combination is commonly used, especially in combination with batch learning. Following the above, it is important to note that the context in this manuscript is different from that of deterministic, or point, forecast combination. In the context of the latter, the observation that the simple average of experts' forecasts tends to be more accurate than the optimally weighted combination is referred to as the "forecast combination puzzle" and is attributed to the combination variance that is affected by weight estimation [START_REF] Claeskens | The forecast combination puzzle: A simple theoretical explanation[END_REF].

Batch learning

Batch learning requires a separate data set that accurately represents the test set to optimize the combination weights. Consequently, the optimally weighted forecast combination does not necessarily outperform the best component forecast if the separate data sets are not comparable (e.g., [START_REF] Hall | Combining density forecasts[END_REF]. Nevertheless, there is value to optimizing the weights as shown by [START_REF] Martin | Optimal probabilistic forecasts: When do they work[END_REF] who consistently outperform the naive forecast combination with equal weights. A notable example frequently applied to NWP forecasts is Bayesian Model Averaging (BMA), where ensemble members are dressed with a probability density function (PDF) and the weights optimized using the logarithmic score [START_REF] Raftery | Using Bayesian Model Averaging to Calibrate Forecast Ensembles[END_REF]. In their application, BMA is an efficient method to calibrate the ensemble because the ensemble under investigation is underdispersed. To reduce overdispersion, [START_REF] Bracale | A probabilistic competitive ensemble method for short-term photovoltaic power forecasting[END_REF] minimize a weighted sum of the continuous ranked probability score (CRPS) and deviation from calibration. Similarly, Jose, Grushka-Cockayne and [START_REF] Jose | Trimmed opinion pools and the crowd's calibration problem[END_REF] propose the exterior-trimmed opinion pool heuristic method that effectively removes expert forecasts with low and high means or cumulative distribution function (CDF) values. Additionally, they introduce the interior-trimmed opinion pool that removes expert forecasts with moderate means or CDF values to increase the dispersion in case the experts are overconfident.

Averaging probabilities, which is also referred to as vertical averaging, can be shown to be at least as accurate in terms of CRPS as the average CRPS of the experts [START_REF] Lichtendahl | Is it better to average probabilities or quantiles?[END_REF]. Besides averaging probabilities, it is also possible to average quantiles, which can be referred to as horizontal averaging. In their comparative study, [START_REF] Lichtendahl | Is it better to average probabilities or quantiles?[END_REF] show that the average quantile forecast is always sharper than the average probability forecast and that the former is therefore better suited in case the component forecasts are well calibrated. In the area of load forecasting, [START_REF] Wang | Combining probabilistic load forecasts[END_REF] linearly combine quantiles and show that their approach does not always outperform the component models although they do not evaluate forecast calibration. In a similar fashion, [START_REF] Bracale | Developing and comparing different strategies for combining probabilistic photovoltaic power forecasts in an ensemble method[END_REF] minimize the pinball loss to optimally combine quantile forecasts of photovoltaic (PV) power and improve accuracy. However, similar to [START_REF] Wang | Combining probabilistic load forecasts[END_REF], the authors do not evaluate the calibration of the component nor the combined forecasts. [START_REF] Taylor | Combining probabilistic forecasts of covid-19 mortality in the united states[END_REF] forecast the cumulative COVID-19 mortality and face interesting challenges such as a lack of historical forecasts, and therefore use trimming techniques as well as the simple average, the median forecast and weights based on the inverse quantile score to find that the latter performs at least as well as the simple average.

Besides trimming and averaging quantiles, one can also apply nonlinear transformations to the linear opinion pool to improve probabilistic calibration. [START_REF] Gneiting | Combining predictive distributions[END_REF] describe two such methods, namely the spread-adjusted linear pool (SLP) and beta-transformed linear pool (BLP). SLP adjusts the spread of the component forecasts and can consequently mitigate-to a certain extent-overdispersion caused by linearly combining calibrated forecasts. [START_REF] Möller | Probabilistic temperature forecasting with a heteroscedastic autoregressive ensemble postprocessing model[END_REF] apply SLP to post-processed temperature forecasts issued by the European Center for Medium-range Weather Forecasts (ECMWF) ensemble prediction system and show that it effectively lowers CRPS compared to the component forecasts. However, a limitation of SLP is that the method fails to be flexibly dispersive, which is to say that it is unable to sufficiently adjust the spread to produce neutrally dispersed forecasts, especially when the component forecasts are neutrally dispersed or underdispersed [START_REF] Gneiting | Combining predictive distributions[END_REF]. In contrast, BLP is exchangeably flexibly dispersive and as such is able to transform the predictive distributions such that the second moment of the resulting probability integral transform (PIT) can attain any value in the open interval (0, 1∕4), with 1∕12 indicating neutral dispersion [START_REF] Gneiting | Combining predictive distributions[END_REF]. [START_REF] Van Der Meer | Generalizing renewable energy forecasting using automatic feature selection and combination[END_REF] apply both SLP and BLP to combine PV power forecasts and show that SLP outperforms BLP, which is caused by a lack of representativeness of the training data that affects parameter learning more in case of BLP than in case of SLP. Finally, [START_REF] Bassetti | Bayesian nonparametric calibration and combination of predictive distributions[END_REF] develop a Bayesian nonparametric approach that extends the parametric class of calibration functions, i.e., BLP, using a possibly unknown number of beta mixtures, which can be interpreted as a mixture of local combination models.

Online learning

Nonstationary data, extensive training times for complex machine learning models, and data storage present challenges in batch learning. To illustrate this further, consider batch gradient descent. In this method, the objective is to update the model parameter 𝜃 by calculating the loss gradient ∇𝐿(𝑓 (𝑥 𝑖 , 𝜃), 𝑦 𝑖 ) across a training set indexed by 𝑖 ∈ 1, … , 𝑁, which consists of input-output pairs (𝑥 𝑖 , 𝑦 𝑖 ). The parameter update, after processing a batch containing 𝑁 data points, is computed as follows:

𝜃 new = 𝜃 old -𝜂 ⋅ 1 𝑁 ∑ 𝑁 𝑖=1 ∇𝐿(𝑓 (𝑥 𝑖 , 𝜃 old ), 𝑦 𝑖 ),
with 𝜂 representing the learning rate. The aforementioned expression highlights that batch learning necessitates the storage of all 𝑁 data points and the consecutive computation of the gradient 𝑁 times. Unlike batch learning, online learning is computationally less expensive and does not require to store historical data. Moreover, it integrates real-time data and therefore adapts to trends and seasonalities. Returning to the example of gradient descent, the online version instead updates 𝜃 as soon as a new input-output pair (𝑥 𝑖 , 𝑦 𝑖 ) becomes available as 𝜃 new = 𝜃 old -𝜂 ⋅ ∇𝐿(𝑓 (𝑥 𝑖 , 𝜃 old ), 𝑦 𝑖 ). In the context of forecast combination, such real-time adaptation would be relevant when one of the experts is better at forecasting in a declining market rather than a rising market [START_REF] Winkler | Probability forecasts and their combination: A research perspective[END_REF].

Recent years have seen a stark increase in the number of publications on online forecast combination. For instance, [START_REF] Thorey | Online learning with the continuous ranked probability score for ensemble forecasting[END_REF] consider online ridge regression and exponentiated gradient to update weights based on CRPS minimization so as to linearly combine the ensemble members into a CDF. In a subsequent study, [START_REF] Thorey | Ensemble forecast of photovoltaic power with online crps learning[END_REF] extend their research by applying the ML-Poly algorithm introduced by [START_REF] Gaillard | A second-order bound with excess losses[END_REF] to linearly combine ensemble members from ECMWF and Metéo France. This extension demonstrated enhanced performance compared to using raw ensembles, although it should be noted that their forecasts still exhibited underdispersion. In the aforementioned studies, the authors show that the regret of their algorithm is logarithmic as a function of time. It is important to note that in the field of online learning, here used interchangeably with Online Convex Optimization (OCO), and of which forecast combination is a subset, the primary objective is to minimize regret rather than loss. Regret is defined as the difference in performance between the online player and the optimal fixed strategy in hindsight [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF]. Instead of combining ensemble members, V'yugin and Trunov (2019) combine predictive CDFs using Vovk's Aggregating Algorithm and show that their method offers a time-independent upper bound on regret. [START_REF] Zamo | Sequential aggregation of probabilistic forecasts-application to wind speed ensemble forecasts[END_REF] modify the objective function to comprise the CRPS and the Jolliffe-Primo test for rank histogram flatness to improve calibration of their linearly combined forecasts. The aforementioned studies concern linear forecast combination, which, as previously discussed, leads to overdispersion.

Regarding quantile forecasting, [START_REF] Berrisch | CRPS learning[END_REF] observe potential differences in accuracy across various segments of the experts' predictive distributions. To address this, they propose a linear pointwise combination algorithm that aggregates quantiles based on CRPS minimization. Their research reveals that this approach yields a more uniformly distributed loss throughout the predictive distribution. In a related study, [START_REF] Krannichfeldt | Online ensemble approach for probabilistic wind power forecasting[END_REF] adapt the pinball loss to remain "passive" when the loss is below a threshold but "aggresively" adjusts the weights when a new sample causes the loss to exceed the threshold, resulting in improved CRPS and pinball loss although their method appears to be outperformed by benchmark models in terms of calibration.

Contributions

In some applicative fields like wind power forecasting, the value of forecast combination was recognized already 20 years ago, with the first operational models based on spot forecast combination set-up by the Spanish Transmission System Operator (see, e.g., [START_REF] Sánchez | Adaptive combination of forecasts with application to wind energy[END_REF]. Today, this is considered as a mainstream approach in business practices in renewable energy forecasting. Though when it comes to probabilistic forecasting, several research challenges remain. As outlined in the previous sections, the literature on probabilistic forecast combination-both batch and online-is expanding rapidly, which has been predicted based on the popularity of forecasting competitions and developments in machine learning and expert forecasting [START_REF] Winkler | Probability forecasts and their combination: A research perspective[END_REF]. In this work, we extend the beta-transformed linear pool proposed by [START_REF] Gneiting | Combining predictive distributions[END_REF] to the online setting to mitigate miscalibration caused by linear forecast combination. The beta-transformed linear pool is exchangeably flexibly dispersive, meaning that there exists a set of parameters that ensures that the combined forecasts are probabilistically calibrated while the component forecasts are exchangeable [START_REF] Gneiting | Combining predictive distributions[END_REF]. The method that we develop is able to adapt in nonstationary contexts and is integrated into the Online Newton Step (ONS) algorithm that moves in the direction of an approximate Hessian and the gradient, the former of which is additionally used to project the weights back onto the simplex. The CRPS is employed to guide the learning process as it is exponentially-concave, an attribute of the CRPS that permits logarithmic regret when used as a cost function in the ONS algorithm. To summarize, we contribute to the state of the art of online probabilistic forecast combination in the following ways:

• We develop a nonlinear and online method to combine probabilistic forecasts that is exchangeably flexibly dispersive.

• The proposed method relies on the CRPS, an exponentially-concave function, in conjunction with the Online Newton Step algorithm. This combination permits logarithmic regret and accommodates the most comprehensive scenario, wherein experts provide full predictive distributions.

• We demonstrate the effectiveness of our method through two simulation studies from the literature and a realworld wind power forecasting study, and make the code publicly available to facilitate its uptake.

The remainder of this work is organized as follows. Section 2 describes the probabilistic forecast combination framework, i.e., the linear and beta-transformed linear pool, as well as the CRPS that is used to guide the learning process. Next, we introduce online convex optimization, develop the necessary mathematics and introduce the ONS algorithm in Section 3. We present the results of two simulation studies and a real-world wind power case study in Section 4 and conclude this work in Section 5.

Probabilistic forecast combination framework

In this section, we provide a comprehensive explanation of forecast combination. Initially, we delve into the concept of linear combination, followed by an exploration of how the Beta distribution is employed to transform this linear combination. Throughout the remainder of this work, we consider a total of 𝑚 experts indexed by 𝑗 and denote a predictive distribution by F . Furthermore, we use lower case normal font for realizations of scalar variables and upper case for scalar variables, lower case bold font for vectors and upper case bold font for matrices.

Classical linear combination

The linear opinion pool is defined as (e.g., [START_REF] Gneiting | Combining predictive distributions[END_REF]:

F (𝑦) = 𝑚 ∑ 𝑗=1 𝑤 𝑗 F𝑗 (𝑦). (1)
We adopt the abbreviations of [START_REF] Gneiting | Combining predictive distributions[END_REF] in this work. Consequently, in case of equal weights, i.e., 𝑤 𝑗 = 1∕𝑚, we refer to (1) as the ordinary linear pool (OLP), which was labelled the linear opinion pool by [START_REF] Stone | The opinion pool[END_REF]. In contrast, when the weights are optimized with respect to a score, we refer to (1) as the traditional linear pool (TLP). The main motivation to combine forecasts is to harness the wisdom of the crowd, which is an effective method in point forecasting due to the expertise and diversity of the crowd [START_REF] Soll | The "wisdom of crowds" effect[END_REF]. Instead, linearly combining diverse probabilistic forecasts further increases the dispersion, which may lead to calibration issues [START_REF] Hora | Probability judgments for continuous quantities: Linear combinations and calibration[END_REF].

Beta-transformed linear pool

To mitigate miscalibration caused by linear combination, [START_REF] Gneiting | Combining predictive distributions[END_REF] propose a nonlinear transformation of the linear pool ( 1) by means of the Beta CDF. The beta-transformed linear pool (BLP) forecast encapsulates (1) and is defined as [START_REF] Gneiting | Combining predictive distributions[END_REF]:

F𝑎,𝑏 (𝑦) = 𝐼 𝑎,𝑏 ( 𝑚 ∑ 𝑗=1 𝑤 𝑗 F𝑗 (𝑦) ) , ( 2
)
where 𝐼 𝑎,𝑏 is the regularized incomplete beta function with shape parameters 𝑎 and 𝑏. Recall that BLP is able to transform the predictive distributions such that the second moment of the resulting PIT can attain any value in the open interval (0, 1∕4) with fixed weights and 𝑎 > 0 and 𝑏 > 0 [START_REF] Gneiting | Combining predictive distributions[END_REF]. Note that when 𝑎 = 𝑏 = 1, (2) equals (1). Throughout the remainder of this work, we will abbreviate the linear opinion pool as 𝑧 = ∑ 𝑚 𝑗=1 𝑤 𝑗 F𝑗 (𝑦). The regularized incomplete beta function, which is also known as the Beta CDF, is defined as:

𝐼 𝑎,𝑏 (𝑧) = 𝐵 𝑎,𝑏 (𝑧) 𝐵 𝑎,𝑏 , (3)
where

𝐵 𝑎,𝑏 (𝑧) = ∫ 𝑧 0 𝑢 𝑎-1 (1 -𝑢) 𝑏-1 d𝑢 = Γ(𝑎)𝑧 𝑎 2 F 1 (𝑎, 1 -𝑏; 𝑎 + 1; 𝑧), ( 4
)
𝐵 𝑎,𝑏 = ∫ 1 0 𝑢 𝑎-1 (1 -𝑢) 𝑏-1 d𝑢 = Γ(𝑎)Γ(𝑏) Γ(𝑎 + 𝑏) . ( 5
)
In the above equations, 𝐵 𝑎,𝑏 (𝑧) is the incomplete Beta function, 𝐵 𝑎,𝑏 is the complete Beta function, Γ is the gamma function and 2 F 1 is a regularized hypergeometric function. These identities are useful to derive the gradient in Section 3.1.

CRPS-based learning

Strictly proper scoring rules minimize in expectation under the true model. These rules are therefore recommended for forecast evaluation because they encourage truth telling from a forecaster (e.g., Gneiting and Raftery, 2007). Conversely, a forecaster may minimize such a rule to identify the true model in the learning stage. The CRPS is a strictly proper scoring rule that evaluates the entire predictive distribution and therefore permits the most general setting where experts issue complete predictive distributions. Moreover, the CRPS is exponential-concave in case of bounded support, which implies that the score is strongly convex in the direction of the gradient but not necessarily elsewhere [START_REF] Korotin | Mixability of integral losses: A key to efficient online aggregation of functional and probabilistic forecasts[END_REF][START_REF] Hazan | Introduction to Online Convex Optimization[END_REF]. Exp-concavity is an important attribute of the CRPS because-in combination with particular online learning methods-it allows for accelerated learning, as Section 3.2 will clarify. The CRPS is defined as follows:

CRPS( F𝑡 , 𝑦 𝑡 ) = ∫ 1 0 ( F𝑡 (𝑥) -1{𝑥 ≥ 𝑦 𝑡 } ) 2 d𝑥, (6)
where 1 denotes the indicator function that is 1 when the condition inside the curly brackets is true. The limits of the integral are taken to be 0 and 1, respectively, since the data are normalized. Note that ( 6) is an instantaneous value and that these are averaged over a test set to rank competing forecasts. Similarly, one can minimize the expected value of ( 6) to learn the true parameters of a forecast model or, as in our case, learn the optimal combination parameters. Henceforth we will omit the time index to simplify the notation.

Distribution-oriented forecast verification

While the CRPS can be used to rank competing forecasts, it does not reveal specific types of miscalibration. The PIT is a common verification tool to visualize the calibration of probabilistic forecasts and can be computed as 𝑍 = F (𝑌 ), where F and 𝑌 are series of predictive distributions and observations, respectively. F is said to be probabilistically calibrated if 𝑌 ∼ F , which implies that 𝑍 is a standard uniform distribution [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF]. Consequently, the variance of 𝑍 is constrained to the closed interval [ 0, 1∕4 ] and the distribution is flat when var(𝑍) = 1∕12. Nevertheless, a PIT histogram can deviate from flatness due to randomness induced by a test set of limited length. To account for randomness, [START_REF] Bröcker | Increasing the reliability of reliability diagrams[END_REF] propose consistency bands that represent the maximum estimation uncertainty that can be expected for a test set of specific length. In other words, as long as the deviation from flatness remains within the consistency bands, one cannot reject the hypothesis that the probabilistic forecasts are reliable. Herein, we use dotted lines to visualize the consistency bands.

Online convex optimization

Online convex optimization (OCO) can be seen as a game where a player repeatedly makes decisions over time. At time step 𝑡, the player chooses from the 𝑛-dimensional convex decision set 𝒙 𝑡 ⊆ ℝ 𝑛 . After the player's decision, a convex cost function 𝑓 𝑡 is revealed and the player incurs a loss 𝑓 𝑡 (𝒙 𝑡 ). The performance of OCO algorithms is measured in terms of regret, which is defined as the difference in performance between the online player and the optimal fixed strategy in hindsight. In the context of online forecast combination, the forecast aggregator is the online player who adapts the aggregation function based on the most recent performance and who regrets-in hindsight-not choosing the fixed aggregation function that optimizes performance. The reason for using regret is that OCO permits an adversary who can choose different cost functions as the game progresses. In this framework, OCO is concerned with algorithms that realize sublinear regret with increasing test length, implying that, on average, the algorithms perform as well as the best strategy in hindsight. For an in-depth treatment of OCO and its various algorithms, the reader is referred to [START_REF] Orabona | A modern introduction to online learning[END_REF] and [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF].

Although there are settings in OCO, such as Multi-Armed Bandit, where a player does not have access to the loss gradient, we focus on gradient-based online learning. One advantage of gradient-based online learning is that it generally permits tighter bounds on regret. The properties of the CRPS described in Section 2.2 make it an attractive scoring rule in an online learning context as well as for evaluating forecasts. What follows is the derivation of the CRPS gradient with respect to the shape parameters of the Beta CDF and the linear combination weights. Afterwards, we introduce the OCO algorithm that we employ.

Derivation of the CRPS gradient

Recall that we consider 𝑚 experts indexed by 𝑗; the vector of weights is subsequently defined as

𝒘 = ( 𝑤 1 , … , 𝑤 𝑚 ) ⊤
and the vector of parameters is defined as 𝒙 = ( 𝑎, 𝑏, 𝒘 ⊤ ) ⊤ . The gradient of the CRPS with respect to 𝒙 is then defined as follows:

∇CRPS( F𝑎,𝑏 , 𝑦) = ( 𝜕CRPS 𝜕𝑎 , 𝜕CRPS 𝜕𝑏 , 𝜕CRPS 𝜕𝑤 1 , … , 𝜕CRPS 𝜕𝑤 𝑚 ) ⊤ . (7)

Partial derivative with respect to the weights

In deriving the partial derivatives, we first consider the linear combination of the component forecasts as it is encapsulated in the beta-transformed linear pool. This is relevant because the linear pool can act as a benchmark for the beta-transformed linear pool.

Linear The combination weights are contained within the integral, cf. ( 6). We therefore use the Leibniz integral rule for differentiation under the integral sign to derive the partial derivatives with respect to weights 𝑤 𝑗 since they are not integration variables. To improve the numerical properties of the learning process, [START_REF] Pinson | Adaptive modelling and forecasting of offshore wind power fluctuations with markov-switching autoregressive models[END_REF] propose a logit transform w𝑗 of 𝑤 𝑗 , such that:

w𝑗 = ln ( 𝑤 𝑗 1 -𝑤 𝑗 ) , ( 8
)
which constrains 𝑤 𝑗 to the open interval (0, 1). Using the chain rule, the partial derivative of the CRPS with respect to 𝑤 𝑗 becomes:

𝜕CRPS 𝜕𝑤 𝑗 = 𝜕CRPS 𝜕 w𝑗 𝜕 w𝑗 𝜕𝑤 ,
which means that:

𝜕CRPS 𝜕 w𝑗 = 𝜕CRPS 𝜕𝑤 𝑗 𝜕𝑤 𝑗

𝜕 w𝑗 .

Notice that

𝜕 w𝑗 𝜕𝑤 𝑗 = 1 𝑤 𝑗 -𝑤 2

𝑗

. We have:

𝜕CRPS 𝜕𝑤 𝑗 = ∫ 1 0 2 ( FTLP (𝑥) -1{𝑥 ≥ 𝑦} ) 𝜕 FTLP (𝑥) 𝜕𝑤 𝑗 d𝑥 = 2 ∫ 1 0 F𝑗 (𝑥) ( FTLP (𝑥) -1{𝑥 ≥ 𝑦} ) d𝑥, (9)
and therefore:

𝜕CRPS 𝜕 w𝑗 = 2 ( 𝑤 𝑗 -𝑤 2 𝑗 ) ∫ 1 0 F𝑗 (𝑥) ( FTLP (𝑥) -1{𝑥 ≥ 𝑦} ) d𝑥. ( 10
)
Nonlinear The partial derivative in the nonlinear setting is similar to (9), except that FTLP (𝑥) is replaced with F𝑎,𝑏 (𝑧).

Notice that 𝑧 appears in the upper limit of the integral in (4) and we therefore require Leibniz's integral rule. We defer the derivation to Appendix A.1 and instead present the final result below:

𝜕CRPS 𝜕 w𝑗 = 2 ( 𝑤 𝑗 -𝑤 2 𝑗 ) 𝐵 𝑎,𝑏 ∫ 1 0 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) 𝑧 𝑎-1 (1 -𝑧) 𝑏-1 F𝑗 (𝑥) d𝑥. (11)

Partial derivative with respect to shape parameter 𝑎

Similarly, shape parameter 𝑎 is not an integration variable and we therefore apply the Leibniz integral rule again. Important to note here is that both 𝑎 and 𝑏 are strictly positive. Similar to weight 𝑤 𝑗 , we apply a change of variable ã = ln(𝑎) to improve the stability of the algorithm [START_REF] Pinson | Adaptive modelling and forecasting of offshore wind power fluctuations with markov-switching autoregressive models[END_REF]. The partial derivative with respect to shape parameter 𝑎 is presented below, whereas the derivation is deferred to Appendix A.2:

𝜕CRPS 𝜕 ã = 2𝑎 ∫ 1 0 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) × ( F𝑎,𝑏 (𝑧) (ln(𝑧) -𝜓(𝑎) + 𝜓(𝑎 + 𝑏)) - Γ(𝑎)Γ(𝑎 + 𝑏) Γ(𝑏) 𝑧 𝑎 3 F 2 (𝑎, 𝑎, 1 -𝑏; 𝑎 + 1, 𝑎 + 1; 𝑧) ) d𝑥, ( 12
)
where 𝜓 is the digamma function.

Partial derivative with respect to shape parameter 𝑏

A similar approach to derive 𝜕CRPS∕𝜕𝑏 can be used as was used in the previous section. Furthermore, an identical change of variable b = ln(𝑏) is used to improve the stability of the algorithm. Equation ( 13) presents the result, while the derivation is deferred to Appendix A.3.

𝜕CRPS 𝜕 b = 2𝑏 ∫ 1 0 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) × ( Γ(𝑎 + 𝑏)Γ(𝑏) Γ(𝑎) (1 -𝑧) 𝑏 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) + F𝑏,𝑎 (1 -𝑧) (𝜓(𝑏) -𝜓(𝑎 + 𝑏) -ln(1 -𝑧)) ) d𝑥. (13)

Online Newton step

Several algorithms have been proposed in the field of online convex optimization that are analogues to wellknown offline algorithms (see, e.g., [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF], for an overview). Online Gradient Descent (OGD), as described in Section 1.2, is the analogue to Gradient Descent in which the algorithm moves in the direction of the gradient at every iteration and is therefore a first-order algorithm [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]. Despite its simplicity, OGD is linear in time and achieves sublinear regret in the worst case and logarithmic regret for strongly convex loss functions [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF].

One difficulty associated with OGD is the need for precise tuning of the step size or learning rate to achieve the desired regret. Additionally, it may be that the requirement of strong convexity is too stringent in practice. In such instances, it is possible to utilize exp-concave loss functions that are strongly convex in the direction of the gradient but not necessarily elsewhere [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF]. The Online Newton Step (ONS) is an algorithm that guarantees logarithmic regret for exp-concave loss functions and therefore does not require an adaptive step size [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF]. ONS is analogues to the Newton-Raphson method in that it moves in the direction of an approximated Hessian and the gradient, i.e., 𝑨 -1 𝑡 ∇ 𝑡 , yet is based only on first-order information. The operation 𝑨 -1 𝑡 ∇ 𝑡 can potentially yield a set of weights that lie beyond the boundaries of the unit simplex. To avoid constraints, [START_REF] Pinson | Adaptive modelling and forecasting of offshore wind power fluctuations with markov-switching autoregressive models[END_REF] parameterize the transition probabilities on the 𝑚-dimensional unit sphere. Instead, we perform a projection onto the unit simplex with a norm induced by 𝑨 𝑡 rather than the Euclidean norm [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF]. More details on this projection will be given below.

Recall that we consider 𝑚 experts, resulting in a weight vector 𝒘 ∈ ℝ 𝑚 and a parameter vector 𝒙 ∈ ℝ 𝑛 that includes the shape parameters of the Beta distribution (3). Algorithm 1 summarizes the ONS algorithm in detail. Note that 𝛾 is a scaling factor of the starting point of the update step, i.e., 𝑨 -1 0 , whereas 𝜂 represents. The step size remains constant for two reasons: (i) ONS guarantees logarithmic regret for exp-concave losses and (ii) maintaining a fixed step size allows the algorithm to adapt effectively in non-stationary environments, whereas reducing step sizes could hinder this adaptability. To avoid inverting a potentially large matrix, [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF] recommend a recursion step for 𝑨 -1 𝑡 using the Sherman-Morrisson formula [START_REF] Sherman | Adjustment of an Inverse Matrix Corresponding to a Change in One Element a Given Matrix[END_REF]:

𝑨 -1 𝑡 = ( 𝑨 𝑡 + ∇ 𝑡 ∇ ⊤ 𝑡 ) -1
First Author et al.: Preprint submitted to Elsevier Algorithm 1: Online Newton step [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF][START_REF] Wintenberger | Lecture notes in Online Convex Optimization[END_REF])

Data: convex set , 𝑇 , 𝒙 1 ∈  ⊆ ℝ 𝑛 , parameters 𝛾, 𝜂 > 0, 𝑨 0 = 1∕𝛾 2 𝑰 𝑛 , 𝑨 -1 0 = 𝛾 2 𝑰 𝑛 for 𝑡 ← 1 to 𝑇 do
Play 𝒙 𝑡 and observe cost 𝑓 𝑡 (𝒙 𝑡 );

𝑨 𝑡 = 𝑨 𝑡-1 + ∇ 𝑡 ∇ ⊤ 𝑡 ; 𝑨 -1 𝑡 = 𝑨 -1 𝑡-1 - 𝑨 -1 𝑡-1 ∇ 𝑡 ∇ ⊤ 𝑡 𝑨 -1 𝑡-1 1+∇ ⊤ 𝑡 𝑨 -1 𝑡-1 ∇ 𝑡 ; Newton step: 𝒚 𝑡+1 = 𝒙 𝑡 -𝜂 1 𝛾 𝑨 -1 𝑡 ∇ 𝑡 ; 𝒗 𝑡+1 ∈ ℝ 𝑚 ⊂ 𝒚 𝑡+1 ∈ ℝ 𝑛 ; Projection (weights only) with weighted norm ‖⋅‖ 𝑫 𝑡 : 𝒘 𝑡+1 = 1 2 arg min 𝒘∈ ‖𝒘 -𝒗 𝑡+1 ‖ 2 𝑫 𝑡 end = 𝑨 -1 𝑡-1 - 𝑨 -1 𝑡-1 ∇ 𝑡 ∇ ⊤ 𝑡 𝑨 -1 𝑡-1 1 + ∇ ⊤ 𝑡 𝑨 -1 𝑡-1 ∇ 𝑡 . ( 14
)
In this work, the weights should sum to 1. Although not a strict requirement, we additionally ensure that the weights are positive. Whereas OGD uses a Euclidean projection step, ONS requires a generalized projection onto the unit simplex Λ that is in the vector norm induced by 𝑨 𝑡 , which is a convex program [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF]. Nonetheless, the presence of off-diagonal elements within 𝑨 𝑡 introduces complexity to the optimization problem. Therefore, it is more advantageous to approximate it by employing

𝑫 𝑡 = diag ( 𝑨 𝑡 )
. Consequently, the generalized projection onto Λ breaks down into a series of independent scalar minimization problems for each component of 𝒘. These can be resolved by solving a piecewise linear equation through sorting, as described by [START_REF] Held | Validation of subgradient optimization[END_REF]. Specifically, we want to solve the following convex optimization problem [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF])

𝒘 𝑡+1 = 1 2 arg min 𝒘∈ ‖𝒘 -𝒗 𝑡+1 ‖ 2 𝑫 𝑡 , ( 15
)
where 𝒗 𝑡+1 ∈ ℝ 𝑚 comprises only the weights and ‖𝒗 𝑡+1 -𝒘‖ 2

𝑫 𝑡 ∶= ( 𝒘 -𝒗 𝑡+1 ) ⊤ 𝑫 𝑡 ( 𝒘 -𝒗 𝑡+1 )
. The piecewise linear equation can be derived by solving a system of linear equations obtained from the Karush-Kuhn-Tucker (KKT) optimality conditions (see Appendix B), which results in the weighted soft-threshold [START_REF] Wintenberger | Lecture notes in Online Convex Optimization[END_REF])

𝒘 * = 𝑫 -1 𝑡 SoftThreshold ( 𝑫 𝑡 𝒗 𝑡+1 , 𝜈 * ) . ( 16
)
Algorithm 2 summarizes the aforementioned steps. Note that in our application, the convex set  ⊆ ℝ 𝑛 is unbounded and open since the shape parameters 𝑎 and 𝑏 are strictly larger than 0. We therefore forego a regret analysis and instead empirically show the efficacy of the proposed method.

Results

We illustrate the efficacy of the proposed method through three case studies. The first two case studies are based on synthetic data sets adapted from relevant studies, a time-invariant (Section 4.2.1) and a time-varying process (Section 4.2.2). The synthetic case studies enable us to assess the performance of the proposed method within a controlled environment. To facilitate uptake of the proposed method, we make the code to generate the synthetic data and run the experiments available on GitHub1 . Finally, we apply the proposed method to a real-world wind power forecast case study in Section 4.3.

Benchmarks

In order to evaluate the effectiveness of the proposed online combination method, referred to as BLP, we incorporate several benchmark models. The most obvious benchmarks are the component models that BLP is designed to combine, Algorithm 2: Simplex projection with weighted norm ‖⋅‖ 𝑫 𝑡 [START_REF] Wintenberger | Lecture notes in Online Convex Optimization[END_REF])

Data: 𝒘 ∈ ℝ 𝑚 and 𝑫 𝑡 = diag ( 𝑨 𝑡 ) if 𝒘 ∈ Λ then Return 𝒘; else Sort (𝑑 𝑡 𝑣 𝑡+1 ) 1 ≥ … ≥ (𝑑 𝑡 𝑣 𝑡+1 ) 𝑚 ; Find 𝑑 0 = max { 1 ≤ 𝑗 ≤ 𝑚; (𝑑 𝑡 𝑣 𝑡+1 ) 𝑖 - 1 ∑𝑑 0 𝑖=1 𝑑 -1 𝑡,𝑖 ( ∑ 𝑑 0 𝑖=1 𝑣 𝑡+1,𝑖 -1 ) } ; Define 𝜈 * = 1 ∑𝑑 0 𝑖=1 𝑑 -1 𝑡,𝑖 ( ∑ 𝑑 0 𝑖=1 𝑣 𝑡+1,𝑖 - 1
) ;

Return 𝒘 * = 𝑫 -1 𝑡 SoftThreshold ( 𝑫 𝑡 𝒗 𝑡+1 , 𝜈 * ) end and these will be elaborated upon in the respective sections. Furthermore, there are three combination models characterized by fixed weights, meaning that their weights remain constant throughout the entire duration. These combination models are OLP with equal weights, TLP * (cf. ( 1)), and BLP * (cf. ( 2)). Note that the asterisk indicates that the model parameters are optimized in hindsight, meaning that on average they are optimal for the testing data set. Finally, we include two benchmarks based on online learning: (i) TLP, and (ii) the CRPS learning approach by [START_REF] Berrisch | CRPS learning[END_REF], available on CRAN2 , with all options set to the default values and which we refer to as PFC. Note that the parameters required by the online learning methods are updated at every time step.

Synthetic data

For both synthetic case studies, we initialize the weights as 1∕𝑚 and the shape parameters 𝑎 = 𝑏 = 1, which effectively is the ordinary linear pool (OLP). In total, we run 150 simulations for 11,500 time steps and report the median as well as the 90% confidence interval of the parameter estimates. Regarding the hyperparameters of the ONS algorithm, we perform an exhaustive grid search for 𝛾 and 𝜂 and select the hyperparameters that minimize CRPS. Because of the different nature of the case studies, these parameter grids are not identical.

Time-invariant process

We adapt a simulation study from [START_REF] Gneiting | Combining predictive distributions[END_REF], in which the data generating process is the combination of standard normal random variables 𝑋 0 , 𝑋 1 , 𝑋 2 , 𝑋 3 and 𝜖 without a temporal relation:

𝑌 = 𝑋 0 + 𝑎 1 𝑋 1 + 𝑎 2 𝑋 2 + 𝑎 3 𝑋 3 + 𝜖, ( 17
)
where 𝜖 represents an error term. [START_REF] Gneiting | Combining predictive distributions[END_REF] note that 𝑋 0 may represent public information, such as weather forecasts from a meteorological office, while 𝑋 1 , 𝑋 2 and 𝑋 3 represent private information measured by forecasters 1, 2 and 3, respectively. These experts issue the following probabilistic forecasts:

𝑓 1 =  ( 𝑋 0 + 𝑎 1 𝑋 1 , 1 + 𝑎 2 2 + 𝑎 2 3 ) , ( 18
)
𝑓 2 =  ( 𝑋 0 + 𝑎 2 𝑋 2 , 1 + 𝑎 2 1 + 𝑎 2 3 ) , ( 19
)
𝑓 3 =  ( 𝑋 0 + 𝑎 3 𝑋 3 , 1 + 𝑎 2 1 + 𝑎 2 2 ) , ( 20
)
where 𝑎 1 = 𝑎 2 = 1 and 𝑎 3 = 1.1. Consequently, the component forecasts are probabilistically calibrated by design and their linear combination leads to overdispersion [START_REF] Gneiting | Combining predictive distributions[END_REF]. As mentioned, we compare our method with linear and nonlinear combinations whose weights are static over time, denoted with an asterisk. To that end, the parameters are optimized over the test data by minimizing the logarithmic score as in [START_REF] Gneiting | Combining predictive distributions[END_REF].

Simulation results

Given the stationary and time-invariant nature of this case study, the hyperparameter grid includes low 𝛾 values for increased stability and relatively high 𝜂 values to learn quickly. Specifically, the grid comprises the 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.00 Figure 1a presents the evolution of the parameters over time during linear and nonlinear forecast combination. The shaded area represents the 90% confidence interval associated to the parameter, computed over 150 experiments. When considering linear forecast combination, it becomes apparent that the third expert typically receives the highest weight. This pattern is expected because the third expert consistently provides slightly sharper predictive distributions. The weights of experts 1 and 2 are similar, which can be explained by the fact that their forecasts are identical in expectation. We observe a similar pattern in the case of nonlinear forecast combination. However, both shape parameters are larger than 1, indicating that the Beta transformation sharpens the predictive distributions. Crucially, the shape parameters are very similar and this is expected as their primary purpose is to enhance sharpness.

Figure 1b presents the PIT histograms for the experts and the combination methods. Recall that departures from flatness may arise from a limited-sized test set, which is why Fig. 1b also includes consistency bars. This figure demonstrates that the experts exhibit the intended probabilistic calibration, highlighting that the linear combination methods (namely, OLP, TLP, and TLP * ) result in overdispersed forecasts. In contrast, both the online BLP and the fixed BLP * effectively transform the combined forecasts to generate calibrated forecasts, whereas the PFC model tends to be slightly miscalibrated at the outer edge of the distributions.

Considering CRPS, Table 1 shows that all combination methods substantially outperform the component models. Furthermore, it can be seen that OLP, TLP and TLP * perform quite similar. BLP and BLP * perform identically, which can partly be attributed to the fact that the process is time-invariant. In other words, the ONS algorithm converges to the true weights, which implies that logarithmic regret is hereby empirically shown. Surprisingly, PFC underperforms on this synthetic data set although it is important to recall that all settings have been left to their default values.

Table 1

Realized CRPS for the component models and combination methods. The results are presented as the mean and standard error of 150 experiments based on synthetic data adapted from [START_REF] Gneiting | Combining predictive distributions[END_REF]. Note that the asterisk indicates a combination method with optimal parameters in hindsight. Finally, we evaluate the computational complexity of Algorithm 1, which is dominated by the computation of the integrals (11), ( 12) and ( 13). On a 2020 M1 MacBook Pro running R 4.2.2, it takes approximately 10 ms (per component model), 33 ms and 44 ms to evaluate the aforementioned integrals, respectively. One iteration of Algorithm 1 requires, on average, 110 ms when considering 3 component models.

Time-varying process

We additionally illustrate the efficacy of the proposed method on a simulation study adapted from [START_REF] Berrisch | CRPS learning[END_REF]. Specifically, the data generating process is defined as:

𝑌 𝑡 =  ( 0.15asinh ( 𝜇 𝑡 ) , 1 ) , ( 21
)
where 𝜇 𝑡 = 0.99𝜇 𝑡-1 + 𝜖 𝑡 and 𝜖 𝑡 a standard normal random variable. In this study, two experts provide constant probabilistic forecasts:

𝑓 1 =  (-1, 1) , ( 22
)
𝑓 2 =  (3, 4) . ( 23
)
Similar to the previous case study, we compare our method with linear and nonlinear combination methods whose parameters are optimal in hindsight, learned by minimizing the logarithmic score over the test set and indicated by an asterisk.

Simulation results

Here, we select a grid that includes higher 𝛾 values than in the previous case study to increase adaptation speed, while retaining the same 𝜂 values. Specifically, the grid comprises the combinations of 𝛾 ∈ {0.175, 0.200, 0.225, 0.250, 0.275, 0.300, 0.325} and 𝜂 ∈ {0.150, 0.175, 0.200, 0.225, 0.250}. The combination that minimizes CRPS is 𝛾 = 0.275 and 𝜂 = 0.175.

The forecasts by the two experts in this case study are both biased; however, the first expert does issue forecasts that are correctly dispersed, whereas the second expert issues forecasts that are overdispersed. Consequently, Fig. 2a shows that both linear and nonlinear combination weigh the first expert more than the second (recall that the shaded area represents the 90% confidence interval associated with the parameters). To optimize the calibration of the resulting combination, TLP (the linear combination method) distributes the weights to include both experts to a degree. In contrast, BLP quickly disregards the second expert and instead uses the shape parameters to modify the predictive distributions. Note that the shape parameters are quite different from each other; this is because they have to account for the bias of the first expert.

Figure 2b reveals the miscalibration of the two experts. Their miscalibration is such that linearly combining the forecasts does not result in calibrated forecasts, which can be expected because of the overdispersion of the second expert and their opposing biases. Nonlinear combination significantly improves calibration and the online method is preferred over the static method (indicated with an asterisk). Nevertheless, these combined forecasts are not perfectly calibrated at the edges of the distributions. In this case study, the PFC generates forecasts that are closest to perfect calibration although also here some deviation can be observed at the edges.

In terms of CRPS, Table 2 shows that BLP * outperforms the online method, i.e., BLP. This is not entirely unexpected since the data generating process is stationary, meaning that the best strategy in hindsight can be expected to be highly competitive. The PFC method outperforms others, primarily because it combines quantiles, allowing it to incorporate the strengths of component forecasts exhibiting biases in opposing directions. In terms of linear combination, we observe that TLP outperforms TLP * , indicating that online learning does benefit linear forecast combination. 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.0 0.1 0.2 0.3

Probability integral transform Relative frequency (-)

Figure 2: In (a), the evolution of the weights over time organized by linear and nonlinear combination (rows), where the shaded area represents the 90% confidence interval computed over 150 experiments. In (b), the PIT histograms of the experts and combination methods again with the asterisk indicating the optimal parameters in hindsight.

Table 2

Realized CRPS for the component models and combination methods. The results are presented as the mean and standard error of 150 experiments based on synthetic data adapted from [START_REF] Berrisch | CRPS learning[END_REF]. Note that the asterisk indicates a combination method with optimal parameters in hindsight.

F1

Real-word time series

We apply the proposed combination method on probabilistic power forecasts issued for a wind farm located in midwest France. The wind farm has a nominal capacity of 16,000 kW and the data set ranges from 2018-09-30 until 2020-09-30 at 15 min resolution. To flag anomalies, we use the OpenOA library in Python and linearly interpolate flagged instances [START_REF] Perr-Sauer | OpenOA: An Open-Source Codebase For Operational Analysis of Wind Farms[END_REF].

For brevity, we consider only 15 min ahead, 3 h, 6 h and 24 h ahead forecasts. Especially further into the future, NWP forecasts are essential. Furthermore, [START_REF] Winkler | Probability forecasts and their combination: A research perspective[END_REF] argue that weighted forecast combination of experts that are highly correlated can become unstable. We therefore include forecasts from three NWP models, namely (i) the High-Resolution Forecast (HRES) from the European Centre for Medium-Range Weather Forecasts (ECMWF), (ii) the Global Forecast System (GFS) from the National Centers for Environmental Prediction (NCEP), and (iii) the global Arpège model from Metéo France. The average of the four NWP grid points closest to the wind farm, as well as the last observed power and wind speed, serve as input to three machine learning models that we describe next. The last stage before forecast combination is to use post-processing models to convert the wind speed forecasts into probabilistic wind power forecasts that are to be combined. For this, we employ three models, specifically (i) Generalized Boosted Regression Models (GBM), (ii) Quantile Regression Forests (QRF), and (iii) Quantile Regression (QR). GBM is a technique where simple models, here regression trees, are repeatedly combined in a stage-wise fashion to minimize the prediction loss [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF]. In contrast, QRF is based on random forests in which independent regression trees are fitted on subsets of the feature space and sample space, making them easily parallelizable [START_REF] Meinshausen | Quantile regression forests[END_REF]. Finally, QR is a linear model where the pinball loss is minimized to predict conditional quantiles [START_REF] Koenker | Regression Quantiles[END_REF]. A thorough treatment of these models is out of the scope of this manuscript and we refer to standard textbooks such as [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Furthermore, it is important to note that the aim of this study is not to generate the most accurate component forecasts; rather, the aim is to show the effectiveness of probabilistic forecast combination with our proposed method. The aforementioned models have been selected because they are well documented and can easily be used in R or Python. We use the R packages gbm to implement GBM [START_REF] Greenwell | gbm: Generalized Boosted Regression Models[END_REF], quantregForest to implement QRF [START_REF] Meinshausen | quantregForest: Quantile Regression Forests[END_REF] and quantreg to implement QR [START_REF] Koenker | quantreg: Quantile Regression[END_REF]. Finally, we use the function contCDF of the R package ProbCast to convert quantile forecasts into continuous CDFs with generalized Pareto distribution tails (Browell, Gilbert, McFadzean and Tawn, 2022).

We perform an exhaustive grid search to select the hyperparameters that minimize CRPS for each machine learning model. Given that there are in total 216 combinations to be validated3 , we restrict ourselves to a single forecast horizon (3 h) and a single wind park, thus resulting in 54 combinations. This strategy can be motivated by the fact that the post-processing models are expected to have a similar learning pattern across horizons, except for perhaps the first forecast horizon. For hyperparameter validation, we train each model on 5 weeks of data and validate on the following 3 weeks, which is repeated 12 times in a rolling fashion to cover the entire first year. Similarly, we perform an exhaustive grid search to determine the hyperparameters that minimize CRPS when running the ONS algorithm. Specifically, the grid comprises the combinations of 𝛾 = 𝜂 ∈ {0.10, 0.22, 0.34, 0.46, 0.58, 0.70}. The first four months of the test set are used to validate the ONS algorithm hyperparameters, whereas the remaining eight months are used to test the selected ONS hyperparameters. Table 3 presents the hyperparameter combinations for the linear and nonlinear forecast combination methods and the forecast horizons that minimize CRPS. It is interesting to note that the hyperparameters vary as a function of combination methods and the forecast horizon, indicating that careful validation is necessary. As before, we compare the linear and nonlinear combination methods with the optimal parameters in hindsight, i.e., the parameters are learned based on out-of-sample test data. These static models are denoted with an asterisk.

Measure-oriented forecast analysis

Table 4 presents forecast results across the entire test set for all combinations of weather models, static and online combination methods, as well as the post-processing methods, in terms of CRPS. In the table, the best performing combination of weather model and post-processing model, i.e., expert, is underlined while the best performing combination method is in bold font. Table 4 indicates that there is minimal distinction among the models for the initial forecast horizon. It is worth highlighting that, in this case, the standard deviation of the CRPS surpasses its mean value. An in-depth analysis (which is not presented here) unveiled the skewness in the CRPS distribution. Specifically, the 25 th percentile and median CRPS values are considerably lower than the mean CRPS, and there is also a noteworthy maximum value, approaching approximately 0.76 (varies depending on the model).

Table 4

Average and standard deviation of the CRPS in percentage of nominal capacity, computed over the entire test set. The best performing expert is underlined, whereas the best performing combination method is in bold. Note that the asterisk indicates a combination method with optimal parameters in hindsight. As expected, Table 4 shows that the CRPS increases with the forecast horizon, with the largest relative increase occurring from 15 min to 3 h ahead. The relative increase is especially noticeable for QR, which can be expected since the linear model fails to accurately represent the nonlinear relationship between wind speed and wind power.

Post

Furthermore, the table shows that weather forecasts from ECMWF generally result in more accurate wind power forecasts, whereas weather forecasts from GFS result in less accurate wind power forecasts. Of the combination methods, it can be seen that the linear static method with optimal parameters in hindsight, i.e., TLP * , achieves lower CRPS than the online method, i.e., TLP. Conversely for nonlinear combination, the proposed method BLP improves CRPS up to 4.91% compared to BLP * , indicating that online learning combined with the additional flexibility of the Beta transformation is a valuable improvement when dealing with nonstationary time series. Moreover, the proposed method improves CRPS up to 7.26% when compared to the most accurate expert. In comparison with PFC, the difference in CRPS varies per forecast horizon and is minimal.

Figure 3 presents the CRPS as well, except as cumulative and five-day moving averages for the most accurate expert and the combination methods, and only for the 24 h ahead horizon. In terms of cumulative average, the figure shows that BLP is consistently on par with PFC and that together they outperform the other combination methods, as well as the most accurate expert. Although OLP performs worst in this figure, it is worth noting that it is a competitive method when compared to the experts (cf. Table 4). Nonetheless, when examining the moving average, it becomes evident that adaptive techniques prove beneficial in instances where one or more experts underperform. This is notably the case prior to January and around the outset of April. Particularly during the latter period, it is evident that the fixed method OLP exhibits significantly poorer performance compared to the top-performing expert, along with TLP, PFC, and BLP, with the latter demonstrating the best performance during this period. Concerning the period prior to January, it is noteworthy that BLP exhibits a considerably superior performance compared to the other methods, suggesting the important influence of shape parameters. This aspect will be explored further in the upcoming section. Figure 4 presents the three-day moving average of the weights as they evolve throughout the test period, organized by linear and nonlinear weights, as well as the shape parameters used by the nonlinear forecast combination algorithm for the 24 h forecast horizon.

Combination weights

For the case of linear combination, Fig. 4 shows that expert ECMWF-GBM generally is given the most weight, which is consistent with the results from Table 4. It is worth highlighting, however, the varying weights assigned to experts during different periods. Notably, expert GFS-GBM carries the highest weight during the initial period, while expert MF-QRF assumes significant weighting in the later stages of the test set. Particularly noteworthy is the fact that GFS-GBM maintains a weight exceeding 0.3 for an extended duration, despite not performing as competitively across the entire test set (as shown in Table 4).

When we consider nonlinear combination, as illustrated in Fig. 4, it is evident that expert ECMWF-GBM maintains the highest weighting. However, the remaining experts are allocated substantial weights, with the exception of those relying on QR as a post-processing model. Notably, during the last period of the test set, 4 experts are assigned weights more than 0.1 while the shape parameters are close to being identical and larger than 1. The latter suggests that these experts are overdispersed during the final period and that the proposed algorithm leverages the shape parameters to enhance the probabilistic calibration, similar as in Section 4.2.1.

Distribution-oriented forecast analysis

ECMWF-GBM OLP TLP TLP* BLP BLP* PFC 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.00 Finally, we examine the probabilistic calibration using PIT histograms, which are presented in Fig. 5. Recall that deviations from flatness can be due to a test set of limited size, which is why Fig. 5 additionally presents the consistency bars. The figure illustrates that the most accurate expert aligns closely with probabilistic calibration. Although not shown here, it is worth noting that the majority of experts exhibit satisfactory calibration in the central portion of their predictive distributions. However, all experts do experience slight miscalibration in the outer regions of their predictive distributions.

As the experts are already approaching probabilistic calibration, employing linear combination methods such as OLP, TLP and TLP * leads to overdispersed forecasts. The optimal nonlinear combination in hindsight, BLP * , is closest to perfect calibration, again indicating that the Beta transformation is capable of generating competitive forecasts that are calibrated as well. However, BLP shows signs of miscalibration at the lower end of the predictive distributions, suggesting a positive bias-a characteristic not observed in the predictions of the experts. This discrepancy can likely be attributed to the period up to and including April, during which the shape parameter 𝑎 exceeds shape parameter 𝑏, as depicted in Fig. 4. In general, setting 𝑎 > 𝑏 shifts the predictive distribution upwards, potentially introducing a positive bias when the component forecasts are unbiased. One possible remedy for this issue could involve periodic recalibration of the hyperparameters of ONS. It is worth noting that benchmark PFC lags in terms of probabilistic calibration, exhibiting a notable deviation at the lower end of the predictive distributions.

Conclusions

In this study, we have expanded the application of the beta-transformed linear pool to the online setting by deriving the gradient of the continuous ranked probability score (CRPS) with respect to the shape parameters of the Beta distribution and the combination weights. We have selected the CRPS because it evaluates the entire predictive distribution and therefore yields the most general setting where experts issue complete predictive distributions. In addition, the CRPS is exponential-concave, which allows for accelerated learning in combination with the Online Newton Step (ONS) algorithm. The motivation for our approach is that linear combination, which is the predominant form of forecast combination, always leads to overdispersed forecasts in case the experts are probabilistically calibrated, which is a requirement in probabilistic forecasting. In a time-invariant simulation study using synthetic data, we have shown that the proposed method converges to the optimal combination strategy in hindsight, meaning that the average regret goes to zero. In another simulation study with time-varying synthetic data, we have shown that the proposed method approaches the optimal combination strategy in hindsight but does not achieve zero average regret. Finally, we employed the proposed method to combine probabilistic forecasts of nine experts, resulting from all combinations of three weather models and three post-processing methods, on a real-world wind power data set. We showed that the online and offline combination methods, except the naive version, always perform as well as the best expert. More importantly, the proposed method outperformed the most accurate expert by up to 7.26% and the optimal combination strategy in hindsight by up to 4.91% in terms of CRPS, indicating that online learning combined with the additional flexibility of the Beta transformation is a valuable improvement for nonstationary time series. In all case studies, we observed that the ONS algorithm, which uses the derived gradient to update the parameters, is sensitive to its hyperparameters. Hence, it is important to carefully validate these to optimize the performance of the algorithm.

CRediT authorship contribution statement

A. Derivation of the gradient

A.1. Partial derivative with respect to the weights

Recall that we can compute 𝜕CRPS 𝜕𝑤 𝑗 as follows:

𝜕CRPS 𝜕𝑤 𝑗 = ∫ 1 0 2 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) 𝜕 F𝑎,𝑏 (𝑧) 𝜕𝑤 𝑗 d𝑥, (24)
Also recall that the CDF of the Beta distribution is a quotient where the denominator only depends on 𝑎 and 𝑏. Therefore, we focus on the incomplete beta function 𝐵 𝑎,𝑏 (𝑧). As mentioned, the weights appear in the upper limit of the integral in (4) as 𝑧 = ∑ 𝑚 𝑗=1 𝑤 𝑗 F𝑗 (𝑦) and we therefore require Leibniz's integral rule. In general terms, it is defined as follows:

𝑑 𝑑𝑥 ( ∫ 𝑏(𝑥) 𝑎(𝑥) 𝑓 (𝑥, 𝑡)𝑑𝑡 ) = 𝑓 (𝑥, 𝑏(𝑥)) 𝑑 𝑑𝑥 𝑏(𝑥) -𝑓 (𝑥, 𝑎(𝑥)) 𝑑 𝑑𝑥 𝑎(𝑥) + ∫ 𝑏(𝑥) 𝑎(𝑥) 𝜕 𝜕𝑥 𝑓 (𝑥, 𝑡)𝑑𝑡. ( 25
)
Now, we want to find partial derivative 𝜕 𝜕𝑤 𝑗 , so we plug (4) into (25):

𝜕𝐵 𝑎,𝑏 (𝑧)

𝜕𝑤 𝑗 = 𝜕 𝜕𝑤 𝑗 ( ∫ 𝑧 0 𝑢 𝑎-1 (1 -𝑢) 𝑏-1 𝑑𝑡 ) = 𝑧 𝑎-1 (1 -𝑧) 𝑏-1 𝜕 𝜕𝑤 𝑗 𝑧 -0 𝑎-1 (1 -0) 𝑏-1 𝜕 𝜕𝑤 𝑗 0 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ =0 + ∫ 𝑧 0 𝜕 𝜕𝑤 𝑗 𝑢 𝑎-1 (1 -𝑢) 𝑏-1 𝑑𝑡 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ =0 = ( 𝑚 ∑ 𝑗=1 𝑤 𝑗 F𝑗 (𝑦) ) 𝑎-1 ( 1 - 𝑚 ∑ 𝑗=1 𝑤 𝑗 F𝑗 (𝑦) ) 𝑏-1 F𝑗 (𝑦). (26)
Accordingly, 𝜕𝐼 𝑎,𝑏 (𝑧)∕𝜕𝑤 𝑗 requires multiplication of (26) with 1∕𝐵 𝑎,𝑏 since the denominator in the latter term is not affected by the derivative:

𝜕𝐼 𝑎,𝑏 (𝑧)

𝜕𝑤 𝑗 = 1 𝐵 𝑎,𝑏 𝑧 𝑎-1 (1 -𝑧) 𝑏-1 F𝑗 (𝑦). (27)
Finally, we apply the logit transform (8) to attain the partial derivative of the CRPS with respect to the 𝑗 th weight by plugging ( 27) into (24):

𝜕CRPS 𝜕 w𝑗 = 2 ( 𝑤 𝑗 -𝑤 2 𝑗 ) 𝐵 𝑎,𝑏 ∫ 1 0 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) 𝑧 𝑎-1 (1 -𝑧) 𝑏-1 F𝑗 (𝑥) d𝑥. (28)

A.2. Partial derivative with respect to shape parameter 𝑎

The Leibniz integral rule of 𝜕CRPS∕𝜕𝑎 results in:

𝜕CRPS 𝜕𝑎 = ∫ 1 0 2 ( F𝑎,𝑏 (𝑧) -1{𝑥 ≥ 𝑦} ) 𝜕 F𝑎,𝑏 (𝑧) 𝜕𝑎 d𝑥. (29)
Consequently, it is necessary to compute the partial derivative of (3) with respect to 𝑎. Given that the regularized incomplete Beta function is the ratio of the incomplete Beta function and the complete Beta function, we use the quotient rule. First, we compute 𝜕𝐵 𝑎,𝑏 ∕𝜕𝑎 by applying the product rule to the final equality in (5):

𝜕𝐵 𝑎,𝑏 𝜕𝑎 = 𝜕 𝜕𝑎 Γ(𝑎)Γ(𝑏) Γ(𝑎 + 𝑏) = 𝜕 𝜕𝑎 Γ(𝑎)Γ(𝑏)Γ(𝑎 + 𝑏) -1 .
The partial derivative of Γ(𝑎) with respect to parameter 𝑎 is:

𝜕 𝜕𝑎 Γ(𝑎) = 𝜕 𝜕𝑎 ∫ ∞ 0 𝑢 (𝑎-1) exp(-𝑢) d𝑢 = ∫ ∞ 0 𝜕 𝜕𝑎 𝑢 (𝑎-1) exp(-𝑢) d𝑢 = ∫ ∞ 0 𝑢 (𝑎-1) ln(𝑢) exp(-𝑢) d𝑢 = 𝜓(𝑎)Γ(𝑎), (30)
The partial derivative with respect to 𝑎 of Γ(𝑎 + 𝑏) can be computed in a similar manner such that 𝜕Γ(𝑎 + 𝑏)∕𝜕𝑎 = 𝜓(𝑎 + 𝑏)Γ(𝑎 + 𝑏). Using this result and (30) in the product rule we find:

𝜕𝐵 𝑎,𝑏 𝜕𝑎 = 𝜓(𝑎)Γ(𝑎)Γ(𝑏)Γ(𝑎 + 𝑏) -1 -Γ(𝑎)Γ(𝑏)Γ(𝑎 + 𝑏) -2 𝜕 𝜕𝑎 Γ(𝑎 + 𝑏) = 𝜓(𝑎)Γ(𝑎)Γ(𝑏) Γ(𝑎 + 𝑏) - Γ(𝑎)Γ(𝑏)𝜓(𝑎 + 𝑏) Γ(𝑎 + 𝑏) = Γ(𝑎)Γ(𝑏) Γ(𝑎 + 𝑏) (𝜓(𝑎) -𝜓(𝑎 + 𝑏)) . ( 31
)
The next step is to compute the partial derivative of the incomplete Beta function in (4): ) d𝑥.

𝜕𝐵 𝑎,𝑏 (𝑧) 𝜕𝑎 = 𝜕 𝜕𝑎 ∫ 𝑧 0 ) 2 - Γ (
(42)

B. Weighted projection on the simplex

Recall that the projection of the updated weights 𝒗 𝑡+1 onto Λ with weighted norm ‖⋅‖ 𝑫 𝑡 , where 𝑫 𝑡 = diag ( 𝑨 𝑡 ) , is the following convex optimization problem [START_REF] Wintenberger | Lecture notes in Online Convex Optimization[END_REF])

𝒘 𝑡+1 = 1 2 arg min 𝒘∈ ‖𝒘 -𝒗 𝑡+1 ‖ 2 𝑫 𝑡 . ( 43
)
From this, the can be derived:

where 𝒃 = [ 𝟏 𝑚 ] ⊤ is a vector of ones, 𝜆 𝑗 is a dual variable associated with the nonnegativity constraint of the 𝑗 th component forecast weight and 𝜈 the dual variable associated with the equality constraint. The next step is to compute the gradient:

∇𝐿(𝑤, 𝜆, 𝜈) = ⎛ ⎜ ⎜ ⎝ 𝑫 𝑡 ( 𝒘 -𝒗 𝑡+1 ) -𝜆 + 𝜈𝒃 ⊤ -𝑤 𝒃 ⊤ 𝒘 -1 ⎞ ⎟ ⎟ ⎠ . ( 45
)
Setting the gradient (45) to 0, we find the Karush-Kuhn-Tucker (KKT) conditions: (48) where 𝑑 𝑡,𝑖 is the 𝑖 th diagonal element of 𝑫 𝑡 with identical ordering, which implies that

𝜈 * = 1 ∑ 𝑑 0 𝑖=1 𝑑 -1 𝑡,𝑖 ( 𝑑 0 ∑ 𝑖=1 𝑣 𝑡+1,𝑖 -1
) .

(49)

Figure 1 :

 1 Figure 1: In (a), the evolution of the weights over time organized by linear and nonlinear combination (rows), where the shaded area represents the 90% confidence interval computed over 150 experiments. In (b), the PIT histograms of the experts and combination methods again with the asterisk indicating the optimal parameters in hindsight.

Figure 3 :

 3 Figure 3: Cumulative average and moving average of the best performing expert as well as the online combination methods for the 24 h ahead forecast horizon.

Figure 4 :

 4 Figure 4: The evolution of linear and nonlinear parameters over time, color coded according to the combination of weather model and post-processing model for the 24 h forecast horizon.

Figure 5 :

 5 Figure 5: PIT histograms of the component and combination models again with the asterisk indicating the optimal combination with static weights.

  Dennis van derMeer: Conceptualization, Data Curation, Investigation, Methodology, Software, Visualization, Writing -Original Draft, Writing -Review & Editing. Pierre Pinson: Conceptualization, Funding Acquisition, Writing -Original Draft, Writing -Review & Editing. Simon Camal: Data Curation, Funding Acquisition, Project Administration, Software, Writing -Original Draft. Georges Kariniotakis: Funding Acquisition, Project Administration, Resources, Writing -Original Draft.

=== 1 -

 1 𝑎)𝑧 𝑎 2 F 1 (𝑎, 1 -𝑏; 𝑎 + 1; 𝑧) Γ(𝑎)Γ(𝑏) Γ(𝑎+𝑏) (𝜓(𝑎) -𝜓(𝑎 + 𝑏)) ( 𝑎, 𝑎, 1 -𝑏; 𝑎 + 1, 𝑎 + 1; 𝑧) 𝑏 (𝑧) (𝜓(𝑎) -𝜓(𝑎 + 𝑏)) 𝐵 𝑎,𝑏 = 𝐼 𝑎,𝑏 (𝑧) ln(𝑧) -Γ(𝑎 + 𝑏) Γ(𝑎)Γ(𝑏)𝑧 𝑎 Γ(𝑎) 2 3 F 2 (𝑎, 𝑎, 1 -𝑏; 𝑎 + 1, 𝑎 + 1; 𝑧) -𝐼 𝑎,𝑏 (𝑧) (𝜓(𝑎) -𝜓(𝑎 + 𝑏)) = 𝐼 𝑎,𝑏 (𝑧) (ln(𝑧) -𝜓(𝑎) + 𝜓(𝑎 + 𝑏)) -Γ(𝑎)Γ(𝑎 + 𝑏) Γ(𝑏) 𝑧 𝑎 3 F 2 (𝑎, 𝑎, 1 -𝑏; 𝑎 + 1, 𝑎 + 1; 𝑧). (34)As mentioned in Section 3.1.2, the change of variable ã = ln(𝑎) improves the stability of the algorithm. Given that 𝜕𝑎∕𝜕 ã = 𝑎 and plugging (34) in to (29) in combination with the definition of the beta-transformed linear pool forecast in (2), we find the partial derivative of the CRPS with respect to ã: 𝑏(𝑧) -1{𝑥 ≥ 𝑦} ) 𝜕 F𝑎,𝑏 (𝑧) (𝑧) (ln(𝑧) -𝜓(𝑎) + 𝜓(𝑎 + 𝑏)) -Γ(𝑎)Γ(𝑎 + 𝑏) Γ(𝑏) 𝑧 𝑎 3 F 2 (𝑎, 𝑎, 1 -𝑏; 𝑎 + 1, 𝑎 + 1; 𝑧)derivative of the gamma function with respect to 𝑏 is similar to that in (30), we find that: 𝑏) -𝜓(𝑎 + 𝑏))= 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) . (37)For the partial derivative of the incomplete beta function with respect to 𝑏, we use the following useful property:𝐼 𝑎,𝑏 (𝑧) = 1 -𝐼 𝑏,𝑎 (1 -𝑧) ⇔ 𝐵 𝑎,𝑏 (𝑧) 𝐵 𝑎,𝑏 = 1 -𝐵 𝑏,𝑎 (1 -𝑧) 𝐵 𝑎,𝑏 . (38)From equality (38) it follows that 𝐵 𝑎,𝑏 (𝑧) = 𝐵 𝑎,𝑏 -𝐵 𝑏,𝑎(1 -𝑧). Considering the definition of the incomplete beta function in (4), switching 𝑎 and 𝑏 is important to attain a solution that depends predominantly on 𝑏. Then:𝜕𝐵 𝑎,𝑏 (𝑧) 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) -∫ )𝑢 𝑏-1 (1 -𝑢) 𝑎-1 d𝑢. (39)Similar to the partial derivative with respect to 𝑎, we use integration by parts and some algebra to find that:𝜕𝐵 𝑎,𝑏 (𝑧) 𝜕𝑏 = 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) -1 𝑏 ( ln(𝑢)𝑢 𝑏 2 𝐹 1 (𝑏, 1 -𝑎; 𝑏 + 1; 𝑢) -∫ 𝑢 𝑏-1 2 𝐹 1 (𝑏, 1 -𝑎; 𝑏 + 1; 𝑢) d𝑢 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) -1 𝑏 ( ln(𝑢)𝑢 𝑏 2 𝐹 1 (𝑏, 1 -𝑎; 𝑏 + 1; 𝑢) -𝑢 𝑏 𝑏Γ(𝑏) 2 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 𝑢) 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) + (1 -𝑧) 𝑏 Γ(𝑏) 2 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) -ln(1 -𝑧)𝐵 𝑏,𝑎 (1 -𝑧). (40)Using (37), (40) and the quotient rule, we can compute the partial derivative of the regularized incomplete beta function with respect to 𝑏:𝜕𝐼 𝑎,𝑏 (𝑧) 𝜕𝑏 = ( 𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) + (1 -𝑧) 𝑏 Γ(𝑏) 2 3 F 2 (𝑏,𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) -ln(1 -𝑧)𝐵 𝑏,𝑎 (1 -𝑧) 𝑏 (𝑧)𝐵 𝑎,𝑏 (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) 𝐵 2 𝑎,𝑏 = Γ(𝑎 + 𝑏)Γ(𝑏) Γ(𝑎) (1 -𝑧) 𝑏 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) -ln(1 -𝑧) 𝐵 𝑏,𝑎 (1 -𝑧) 𝐵 𝑎,𝑏 + (𝜓(𝑏) -𝜓(𝑎 + 𝑏)) ( 𝐵 𝑎,𝑏 𝐵 𝑎,𝑏 -𝐵 𝑎,𝑏 (𝑧) 𝐵 𝑎,𝑏 ) First Author et al.: Preprint submitted to Elsevier = Γ(𝑎 + 𝑏)Γ(𝑏) Γ(𝑎) (1 -𝑧) 𝑏 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) + 𝐼 𝑏,𝑎 (1 -𝑧) (𝜓(𝑏) -𝜓(𝑎 + 𝑏) -ln(1 -𝑧)) . (41) The final step is to plug (41) into (36) in combination with the change of variable: 𝑧) 𝑏 3 F 2 (𝑏, 𝑏, 1 -𝑎; 𝑏 + 1, 𝑏 + 1; 1 -𝑧) + F𝑏,𝑎 (1 -𝑧) (𝜓(𝑏) -𝜓(𝑎 + 𝑏) -ln(1 -𝑧))

  𝑗 + 𝜈(𝒃 ⊤ 𝒘 -1),

,(

  ⊤ + 𝜆 * ) 𝒃 ⊤ 𝒘 * = 1 𝑤 * 𝑗 = 0 or 𝑤 * 𝑗 > 0 and 𝜆 * = 0 𝑫 𝑡 𝒗 𝑡+1 , 𝜈 * ) . (47)Subsequently, set ‖𝒘 * ‖ 0 = 𝑑 0 to establish the relation 𝑡,𝑖 𝜈 * ,

Table 3

 3 The combination of 𝛾 and 𝜂 that minimize CRPS for linear and nonlinear forecast combination, i.e., TLP and BLP respectively, and for all forecast horizons.

			TLP				BLP
		15 min	3 h	6 h	24 h 15 min	3 h	6 h	24 h
	𝛾	0.58	0.46 0.70 0.58	0.22	0.70 0.46 0.22
	𝜂	0.58	0.22 0.58 0.34	0.46	0.70 0.70 0.70

www.github.com/xyz First Author et al.: Preprint submitted to Elsevier

https://cran.r-project.org/web/packages/profoc/index.html First Author et al.: Preprint submitted to Elsevier

Four forecast horizons, one wind park, three data sources, three forecast models and six hyperparameter combinations results in 216 combinations.

Acknowledgement

The authors acknowledge ECMWF, NCEP and Metéo France for providing numerical weather predictions. The present research was carried as part of the Smart4RES Project (European Union's Horizon 2020, No. 864337).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The partial derivative can be found through integration by parts and setting 𝑣 = ln(𝑢) and d𝑤 = 𝑢 𝑎-1 (1 -𝑢) 𝑏-1 , which gives d𝑣 = 1∕𝑢 and 𝑤 = 𝑢 𝑎 2 𝐹 1 (𝑎, 1 -𝑏; 𝑎 + 1; 𝑢)∕𝑎. Plugging this in to (32) gives:

The final step to compute the partial derivative of the regularized incomplete Beta function with respect to parameter 𝑎 is, as mentioned at the beginning of this section, to use the quotient rule:


Keywords: spider toxins, dorsal root ganglia neurons, pain, voltage-gated sodium channels, electrophysiology

published or not. The documents may come

Antinociceptive targets for spider toxins in adult dorsal root ganglia neurons

* Corresponding author; E-mail: evelyne.benoit@cea.fr

Although necessary for human survival, pain may sometimes become pathologic if long-lasting and associated with alterations in its signaling pathway. Opioid painkillers are officially used to treat moderate to severe, and even mild, pain. However, the consequent strong and not so rare complications that occur, including addiction and overdose, combined with pain management costs, remain an important societal and economic concern. In this context, animal venom toxins represent an original source of peptides for therapeutic treatment to fight against pain and avoid central side effects, and the most appropriate targets seem to be those of neurons located in the dorsal root ganglia (DRG) which convey pain from the skin and tendons to the central nervous system (CNS). Hence, the DRG neurons are well-known to express various families of transmembrane proteins, including ion channels, G-protein-coupled receptors and gap junctions/pannexins. Among the ion channel family, the most studied target for pain treatment is the voltage-gated sodium (NaV) channel subtype NaV1.7. The attraction of scientists for this channel subtype is due to its validation by human genetic diseases as a pain target. Many studies have been reported in the literature to describe gating modulators or pore blockers that affect the functional properties of NaV1.7 channel subtype. The present overview will focus on the fascinating spider venom toxins which represent an original source of proteins possessing complex structures associated with specific electrophysiological effects and prone to be more selective for the NaV1.7 channel subtype mainly expressed in DRG neurons. It will detail (i) the characteristics of these primary sensory neurons, (ii) the electrophysiological properties of the different NaV channel subtypes expressed in these neurons, with a particular attention on the NaV1.7 subtype, and (iii) the features of spider venom toxins, shaped of inhibitory cysteine knot motif, that present high affinity for the NaV1.7 subtype associated with evidenced analgesic efficacy in animal models.

  Tânia C. Gonçalves, Evelyne Benoit, Michel Partiseti, Denis Servent To cite this version: Tânia C. Gonçalves, Evelyne Benoit, Michel Partiseti, Denis Servent. Antinociceptive targets for spider toxins in adult dorsal root ganglia neurons. RT25 -25th Meeting of the French Society of Toxinology -Toxins and Health: From Molecules to Organisms, Nov 2018, Paris, France. . hal-

2018 34 Antinociceptive targets for spider toxins in adult dorsal root ganglia neurons

  Toxines et Santé : des Molécules aux Organismes -Toxins and Health: from Molecules to Organisms 25 èmes Rencontres en Toxinologie -25 th Meeting on Toxinology Paris, 28 et 29 novembre 2018 -Paris, November 28 th and 29 th ,

	Tânia C. GONCALVES 1,2 , Evelyne BENOIT 2,3,* , Michel PARTISETI 1 , Denis SERVENT 2

1 Sanofi R & D, Integrated Drug Discovery -High Content Biology, F-94440 Vitry-sur-Seine, France. 2 DRF, Institut des Sciences de la Vie Frédéric Joliot, SIMOPRO, CEA de Saclay, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France. 3 Institut des Neurosciences Paris-Saclay, UMR CNRS/Université Paris-Sud 9197, F-91198 Gif-sur-Yvette, France.

This research was funded by a collaborative grant (#153114) between Sanofi Research & Development (Chilly-Mazarin, France) and the French Alternative Energies and Atomic Energy Commission (CEA, Gif sur Yvette, France). T.C.G. was supported by a doctoral CIFRE fellowship from Sanofi.


Jean-Marc Ribero

Jean-Pierre Damiano

Robert Staraj

Accurate original synthesis model for