Lionel Lenôtre
email: lionel.lenotre@polytechnique.edu

Statistical Post-Processing of Weather Forecasts for the Renewable Energy Production by Eolian Devices

Keywords: Classification. 62-02, 62P12. Keyword. Statistical Post-Processing, Survey, Ensemble Forecast, Numerical Weather Predictions, Renewable Energy Production

This paper reviews different statistical methods dedicated to the post-processing of Numerical Weather Predictions and Ensemble Forecast. We focus on the application of the post-processing to problems linked to the production of electricity by eolian devices. The basic idea is to give a concise panorama of the methods commonly used nowadays. We pay a particular attention to the mathematics involved in the methods. We do not compare the methods and do not provide some preferences.

Introduction

Problematic. The renewable energy generation providers use various mathematical tools to manage operations on the spot (daily and intraday) electricity market [START_REF] Bathurst | Trading wind generation in short term energy markets[END_REF][START_REF] Fabbri | Assessment of the cost associated with wind generation prediction errors in a liberalized electricity market[END_REF][START_REF] Lund | Large-scale integration of wind power into different energy systems[END_REF]. Such tools must take into account the uncertainty of the future hourly production rates of renewable energy generation devices with an horizon of approximately 48 hours. For that purpose, the quantification of the uncertainty in weather predictions has to be solidly managed. The required forecast horizon for the spot market is called short-term in meteorology. For such an horizon, the future climatic conditions can be furnished by Numerical Weather Predictions and Ensemble Forecast [START_REF] Costa | A review on the young history of the wind power short-term prediction[END_REF] issued by meteorological centers. Both of these methods are mainly based on physic. Now other methods can be used to obtain short-term weather forecasts such as time series trained on historical data [START_REF] Brown | Time Series Models to Simulate and Forecast Wind Speeds and Power[END_REF] or neural networks [START_REF] Soman | A review of wind power and wind speed forecasting methods with different time horizons[END_REF]. These last two type of methods are statistical and rely on historical data. Hence, pure statistical methods has to be actualized regularly compared to Numerical Weather Predictions. Neural Networks had proved to be very satisfying since it perform online correction of itself. Inclusion of physical concerns can be performed in statistical methods. They often provide improvement and reduce the necessity of actualization. As explained in [START_REF] Costa | A review on the young history of the wind power short-term prediction[END_REF], it is very difficult to determine which methods are the most effective and must be preferred. Here we decide to only deal with numerical weather predictions and Ensemble Forecast.

A Numerical Weather Prediction is a scenario computed with a physical model. It contains forecast values of different weather variables at some time steps whose order is from minutes to few hours. This scenario is purely deterministic and it is not rare that it is completely wrong. A common manifestation of this fact is the lack of trust some people has on weather forecast given by news channel on television. The technique of Ensemble Forecast tries to render the variability of numerical weather predictions by producing a set of different possible scenarios using methods numerical weather predictions. Nevertheless, having different scenarios does not really measure the uncertainty and surely does not provides a clear quantification like a cumulative distribution function or a probability density function do it. Hence, numerical weather predictions and Ensemble Forecast cannot be incorporated directly into useful mathematical tools for operational management such as optimal stochastic control problems. This means that numerical weather predictions and Ensemble Forecast has to be post-processed.

Content. In this survey, we present a few suitable methods for the post-processing of Numerical Weather Predictions and Ensemble Forecast with the aim of solving problems linked to the generation of electricity using eolian devices. Some of them use simple Numerical Weather Predictions while others use Ensemble Forecast which are an improvement of Numerical Weather Predictions. However all the presented methods rely on the same fundamental idea which is to use the forecasts as hints. In other words, the methods infer statistical quantities or functionals of the conditional distribution of a weather variable of interest X knowing some forecasts F = (F 1 , . . . , F N ) coming from an Ensemble Forecast containing N scenarios or from N weather variables correlated to the one of interest in the case of simple Numerical Weather Predictions. Among the methods, a few directly produces results from an Ensemble Forecast while most of them use historical data containing couples of forecasts and realizations. For pretty clear reasons, the methods for the post-processing of Numerical Weather Predictions are only of this last type.

Some the methods presented in this survey were not designed for weather variables directly linked to eolian devices. Since we decide to present the original article which introduce a method, it explain the presence of estimation of statistical quantities or functionals for weather variables such as sea level or temperature in spite of wind speed or wind direction. However, we provides references of the applications or modifications of these methods to the case of wind. At last but not least, we assume weather variables as being in one dimension all along this survey. The extension of some methods to weather variables in high dimension would misleadingly appeared to be very straightforward. Hence, we invite the reader to be very careful.

Two distinctive parts compose this survey. The first one exposed in §2 can be read independently. It gives a small account on models of numerical weather predictions and point out some sources of uncertainty which are encountered in these model. It also present some basic facts about Ensemble Forecast. Precisely, it details the idea behind this method and explain the rendering of the uncertainty of numerical weather predictions by Ensemble Forecasts. In §3, we give defined some scoring rules for the assessment of statistical methods which are used in the sequel to measure the performances of the reviewed statistical methods. The second part is given in §4, §5 and §6. Precisely, we review methods for the computation of quantiles and confidence intervals and the estimation of cumulative distribution functions and probability density functions. We principally detail the mathematical part of these methods and the rationales. We make few comments on the performances of the methods, but we do not enter into details since a large number of articles such as [START_REF] Wilks | Comparison of ensemble-MOS methods in the Lorenz ' 96 setting[END_REF][START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF] provide benchmarks and comparisons.

Numerical Weather Predictions and Ensemble

In various fields of physics such as orbital mechanics, the future state of a system is completely determined by the present state of the system. The evolution of the system is completely governed by causal relationships, that is the laws of physics. This is the path of determinism. In meteorology, it is globally admitted that a fully deterministic method fails to predict future states of the atmosphere. For instance, it is often the case that weather predictions performed are not realized or even that some events cannot be predicted. In this section, we propose to provide some rationales behind the failure of determinism and to present a commonly used workaround solution.

Uncertainty in Numerical Weather Predictions. As we already said above, Numerical Weather Predictions are forecast values of different weather variables at some time steps whose order is from minutes to few hours and whose sum is usually around 48 hours. These forecasts are technically obtained using a numerical model which solves in silico a system of N differential equations coming from a spatial discretization of the laws governing the behavior of the atmosphere [START_REF] Leutbecher | Ensemble forecasting[END_REF]. Such a numerical model is constructed as follows. The physical laws representing the dynamic of the atmosphere [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF] which are

• the conservation of momentum,

• the conservation of mass,

• the equation of state for ideal gases,

• the conservation of energy,

• the conservation equation for water mass. are combined into a physical model depicting a certain interpretation of the behavior of the atmosphere. A filtering approximations [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF] is then performed to allow for the discretization process. The current approximation in use is mainly the hydrostatic one. It has progressively replaced the quasi-geostrophic approximation [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF][START_REF] Lewis | Roots of Ensemble Forecasting[END_REF]. Once these two steps are performed, a mathematical model of the dynamic of the atmosphere is obtained. In order to be solved in silico, the mathematical model is transformed into a numerical model. This step basically consist in choosing a set of numerical methods such as Finite Element or Finite Difference and to set consistent boundary conditions [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF]. These choices are usually made depending on the available computing resources, the desired computational time, etc... On balance, Numerical Weather Predictions comes from the derivation of three nested models, each of which resulting from a sequence of choices whose validity can always be impeached to some extent. In other words, the capability of the numerical model to truthfully render the dynamic of the atmosphere is actually completely uncertain. Moreover, truncation errors and numerical artefacts intrinsic to the algorithms within the numerical methods has to be added.

Once a numerical model has been derived, then it has to be initialized in order to obtain Numerical Weather Predictions In few words, some initial conditions derived from an estimation of the present state of the atmosphere has to be plugged. The estimation procedure must be performed very carefully since the system of differential equations involved in the mathematical model is very chaotic and thus is the numerical model too. For instance, very different results can be obtained from a numerical model while it is simply initialized with multiple rounded estimates of the same initial conditions [START_REF] Lorenz | Atmospheric predictability experiments with a large numerical model[END_REF]. In practice, meteorological centers compute the initial conditions using statistical methods such as Kalman filtering on a combination of observations and very short-term forecast also called first guess [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF]. Such a process is usually named data assimilation in the literature. It follows easily from the description above that the approximated initial conditions are subject to uncertainty coming from the errors in the measurement of observations and the error in the numerical model which is here used to provide the first guess. The iterative use of very short-term forecast can even be seen as the initiator of a particular propagation of the uncertainty through a sort of a loop [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF]Fig. 5.1.2]. The uncertain initial conditions coupled with the amplification property of the numerical model are actually the most important source of uncertainty in the field of Numerical Weather Predictions [START_REF] Leutbecher | Ensemble forecasting[END_REF].

Ensemble Forecast and Uncertainty Quantification. Monte Carlo methods are usually very practical to solve sets of differential equations with uncertain initial conditions. In meteorology, this approach is unfeasible. In fact, the number D of degrees of freedom of the numerical models used is rather too high (D ≈ 10 8 -10 6 ) [START_REF] Leutbecher | Ensemble forecasting[END_REF]. Besides, the probability density function which is used to reflect the uncertainties in the initial conditions and to initialize the Monte Carlo methods is itself uncertain. It also has to be approximated as we are dealing with a problem possessing a high number of degrees of freedom.

Ensemble forecast was designed as a workaround solution in order to replace Monte Carlo methods. It recycles the idea of using random initial conditions. In a nutshell, a numerical model used to provide Numerical Weather Predictions is run several times using small variations of initial conditions obtained through the procedure of data assimilation. The computation of such small variations is performed using deterministic methods [START_REF] Buizza | A Comparison of the ECMWF, MSC, and NCEP Global Ensemble Prediction Systems[END_REF][START_REF] Leutbecher | Ensemble forecasting[END_REF]. Conceptually, the small variations represent the samples of the initial probability density which can not be clearly determined. Each run thus provides samples of each of the different weather variables at some time steps for a certain time lapse. Ensemble forecast can also deal with the uncertainty coming from choices performed in the physical, mathematical and numerical modeling. In fact, small variations can be plugged into different models [START_REF] Lewis | Roots of Ensemble Forecasting[END_REF]. On average, N ≈ 20 to 50 ensemble members are produced due to the computational cost.

An Ensemble Forecast provides raw data which accounts for the uncertainty in the future weather condition. Hence, functionals representing the repartition of the uncertainty such as a cumulative or a probability density function has to be recovered. As multiple outputs from a numerical model, the data produced by an Ensemble Forecast are simply values that can be possibly taken by a weather variable at regularly spaced consecutive times for the next 24 or 48 hours. These data are usually separated into Ensemble members, that is a sequence of values indexed by the time or a path coming from a particular run for each variable. Actually, the most common practice is to infer functionals of the marginal distribution of a weather variable (distribution at a given time step), that is to manipulate probability distributions on R d where d ∈ N. Rarely, the law of the whole path are examined. In other word, statistical methods dealing with probability distribution on some space of functions or sequences are not really in the actual trend. Among the methods for the inference of marginal distributions, some use only the predicted paths of the variable of interest while others add the paths from other weather variables to account for potential correlation. However, almost all of them rely on the same idea of using the paths of Ensemble Forecast as predictors of weather variables. Basically, the common tactic is to infer the conditional marginal distribution of the weather variables knowing the samples coming from Ensemble Forecast. Every methods presented below used this tactic. As an illustration of this principle, the method detailed in §6 all aims at estimating the conditional density f (x|F) of X knowing the N Ensemble members (F j ) 1≤j≤N .

3 Measures of Performance.

In this section, we define the Rank-histogram which is tools designed to examine the qualities of an Ensemble Forecast. Indeed, even with the most advanced post-processing, the quantification of uncertainty can not be well done from an Ensemble Forecast if it is utterly biased or under-dispersive [START_REF] Hamill | Interpretation of Rank Histograms for Verifying Ensemble Forecasts[END_REF]. Then we provide the definition of the two most used scoring rules for the assessment of statistical post-processing. We do not go into detail. Some references are provided for more informations and interpretations of these two scoring rules.

Rank-Based Histogram. In [START_REF] Hamill | Evaluation of Eta-RSM Ensemble Probabilistic Precipitation Forecasts[END_REF], the author define Rank-Histograms and use it as a tool in order to assess the presence of some required probabilistic properties in an Ensemble Forecast. Rank-Histograms are very easy to derive from the data. Let N be the number of members of an Ensemble Forecast and (F j , R j ) 1≤j≤M be M historical pairs of forecast/realization where, for each j = 1, . . . , M , F j = (F 1 , . . . , F n ) is the vector formed with the Ensemble members. For each j = 1, . . . , M , incorporate R j to the vector R j to form a new vector Fj of dimension N + 1. Then sort the elements of Fj in the ascending order and capture the rank of R j . Repeat the procedure M times and form histogram using the relative frequency of apparition of each rank. According to [START_REF] Hamill | Evaluation of Eta-RSM Ensemble Probabilistic Precipitation Forecasts[END_REF], the closest the Histogram is to the one of a uniformly distributed random variable, the better the Ensemble Forecast is. Indeed, a rather high under-dispersitivity or a too large bias would directly overcharge some ranks/beans. The Brier Score. For binary variables, the Brier Score (BS) is a practical equivalent of the squared errors [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF]. It is defined as

BS = M j=1 (p j -o j ) 2
where, for each j = 1, . . . , M , p j is the forecast probability of a boolean event whose realizations are o j . In the literature, this score is sometimes called the half-Brier score [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF] by some authors. The original Brier score together with an extension to categorical variables is given in [START_REF] Gneiting | Strictly Proper Scoring Rules, Prediction, and Estimation[END_REF].

The Brier Score is the most common score in use to measure the accuracy of the estimation of the quantiles of a predictive distribution. Basically, let q be the estimation of the θ-quantile q θ of a variable X where θ ∈ [0, 1] using the method M , then the Brier score of M is

BS M = M j=1 (p -1(R J ≤ q θ )) 2
where R J are observations of X. The lesser the Brier score of M is, the better M is.

The Continuous Ranked Probability Score. The Continuous Ranked Probability Score (CRPS) is defined as

CRPS(F, R) = - R (F (y) -1(y ≥ R)) 2 dy
where F is the cumulative probability function of the distribution used to forecast and R is a realisation [START_REF] Gneiting | Strictly Proper Scoring Rules, Prediction, and Estimation[END_REF]. By looking at the definition of BS M above, one clearly see that the CRPS correspond to the integral of the Brier score extended to continuous random variable. Due to the difficulty to obtain a closed form for the CRPS, it is mostly numerically computed using the formula

CRPS(F, R) = 1 2 E[|X -X |] -E[|X -R|]
where X and X are independent random variables distributed according to the predictive distribution. Indeed, a numerical approximation using a Monte Carlo method can easily be obtain with such a formula. The Continuous Ranked Probability Score offers the advantage of generalizing the absolute error [START_REF] Gneiting | Strictly Proper Scoring Rules, Prediction, and Estimation[END_REF]. Hence, the average of the CRPS of different pairs of forecasts/observations can be compared to the absolute error of a deterministic forecast. For more details on this last point, we refer to [START_REF] Hersbach | Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems[END_REF] which provides a clear presentation of the CRPS and its links with the Brier score and the absolute error.

Computation of Quantiles and Confident Intervals.

Let X be a random variable. In the context of the post-processing of Numerical Weather Predictions or Ensemble Forecast, X is either a weather variable subject to uncertainty or the error between a forecast and a realization considered as a noise. The θ-quantiles of the distribution of X where θ ∈ [0, 1] is defined as the value x such that

P{X ≤ x} = θ.
Quantile forecasts of weather variables related to wind are used in various operational problems on the electricity market such as the quantification of the cost of a suboptimal decision [START_REF] Morales | Integrating renewables in electricity markets -Operational problems[END_REF]. In fact, optimal decision always relates on quantiles with a given coverage rate [START_REF] Gneiting | Quantiles as optimal point forecasts[END_REF]. On the other side, quantiles of the error between forecasts and observations help to provide confidence intervals around some predicted values and thus assess the reliability of the forecast weather condition. Low and stable wind speed are very important for maintenance of the production devices, while very high wind speed has to be taken into account since the devices stop in order not to break.

Numerous methods were developed to obtain quantile forecasts. Here we decide to focus on direct methods, that is those who do not involve cumulative or probability density functions. In fact, some weather variables cannot be handled using standard parametric distributions. For example, the distribution of the error in the case of the wind power is not Gaussian [START_REF] Pinson | On-line assessment of prediction risk for wind power production forecasts[END_REF]. Two very different works are detailed here, one on classical quantiles, the other on confidence intervals. Each of these work presents minimal assumptions. Direct applications are also given in the paper where they were introduced. At last, both of these methods deal with wind power which can be naturally considered has a weather variable since the production rate of an eolian device mainly relates wind speed and direction, temperature and humidity.

Quantiles Estimation using Local Regression. In [START_REF] Bremnes | Probabilistic wind power forecasts using local quantile regression[END_REF], the author proposes to estimate the θ-quantile q θ (F) of the wind power E p produces by eolian devices knowing the numerical weather predictions F = (F 1 , . . . , F K ) of K weather variables which are relevant to predict wind power and thus serve as predictors. The author decides to use the local linear quantile regression [START_REF] Yu | Local Linear Quantile Regression[END_REF] which approximate q θ (F) by a linear functional of the predictors

F ∈ R K , that is q θ (F) = α 0 + α T F,
where the regression parameters α 0 ∈ R and α ∈ M K R, the space of square matrices of size K, are the solutions of the minimization problem arg min

(b 0 ,b) M i=1 ρ(e p,i -(b 0 + b T (F i -F))) w F i -F 2 h λ (F) ,
coming from an empirical version of

q θ (x) = arg min a∈R E[ρ θ (E p -a)|X = x]
ρ θ (u) = uθ, u ≥ 0, u(θ -1), otherwise,
and

w(u) = (1 -u 3 ) 3 , u ∈ [0, 1], 0, otherwise,
• 2 is the Euclidian norm and and h λ is the distance of F to λM -nearest predictor value [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF] where λ ∈ (0, 1] and λM ∈ N.

The above estimation procedure of θ-quantiles is implemented by the author using numerical weather predictions given by Hirlam10. Numerical experiments are performed to show the importance of the selection of predictors. However, the author does not provide a comparison with an alternative method for the numerical computation of quantiles. Such a work is actually performed in an other article by the same author [START_REF] Bremnes | A comparison of a few statistical models for making quantile wind power forecasts[END_REF] where the local linear quantile regression is compared to a local Gaussian model and to the Nadaraya-Watson estimator for cumulative distribution function. Overall, the benchmark does not really discriminate any of the method. We can retain that the local linear quantile regression is a serious direct non-parametric method compared to other two which required the derivation of quantiles from the cumulative distribution function and has parametric assumptions.

Confident Intervals using Resampling. In [START_REF] Pinson | On-line assessment of prediction risk for wind power production forecasts[END_REF], the authors propose a method to compute confidence intervals of wind power forecasts. One of the strength of this method is the possibility to perform On-line update of the confidence intervals. This property is of very high value for decision problems [START_REF] Morales | Integrating renewables in electricity markets -Operational problems[END_REF]. In this method, numerical weather predictions from an Ensemble Forecast are transformed into wind power forecasts which serve as predictors and can be considered as a set of scenario for wind power. This transformation is performed by plugging the relevant weather variables into a power curve of an eolian devices. To compute the confidence intervals, the authors estimate the quantiles of the distribution of the prediction errors

e t+k,t = R t+k -F t+k,t
where F t+k,t is the wind power forecast for the look-ahead time issued at time t and R t+k is the measured wind power at time t + k. The first step of the estimation is to collect the predictions errors made by the model in the past using a sliding windows defined on the index t. The resulting stack of errors is then sorted using fuzzy sets. In a nutshell, the range of values taken by the variable of interest is divided into sets and the errors coming from the previous forecast are kept only if the forecasts belong to the same set as the current one. Two types of sets are defined here. The first concerns the wind speed whose range is divided into two zones depending on the cut-off level, that is the threshold value of wind power where the eolian device stops in order no to break. The second one concerns the range of the wind power which is divided into three zones in order to capture the amplification of the error caused by the power curve. The author choose a particular slicing detailed in [25, pp. 122]. The second step of the method is a resampling procedure which must be performed multiple times. Basically, N values in the stack of errors are drawn randomly with replacement where N is the size of the stack. Then the values are sorted in the ascending order and the λ percent lowest and (1 -λ) percent biggest values where λ represents the percentiles of interest are picked and recorded. Once the resampling step is performed a large number of times, the percentiles are simply given by the averages of the respectively lowest and biggest recorded values.

The method described above is tested by the authors on data coming from Irish wind farms. They provide measures of confidence on the computed intervals [START_REF] Pinson | On-line assessment of prediction risk for wind power production forecasts[END_REF]Table II,pp. 130] which actually speak from themselves and provide a clear assessment of the method. In [START_REF] Pinson | Properties of quantile and interval forecasts of wind generation and their evaluation[END_REF], they more deeply investigate the performance of the method. More precisely, they construct a measuring tool for the assessment of methods which aims at computing quantile forecasts and confidence intervals. They apply it to the present method. Now in [START_REF] Pinson | On-line assessment of prediction risk for wind power production forecasts[END_REF], the authors also construct a method to narrow their confidence intervals. They introduce a so-called meteo-risk index which measures the spread of the forecasts. They use it to define a scale factor which allows to reduce confidence intervals in case of stable weather conditions. Overall, it is the great simplicity of this method which is very seducing compared the previous one we presented. We do no perform a benchmark to compared the methods, but it is clearly an interesting idea for a future work.

Cumulative Distribution Functions and Thresholds

The estimation of cumulative distribution functions is an essential step to numerically solve decision problem such as the Two-Price Market [START_REF] Morales | Integrating renewables in electricity markets -Operational problems[END_REF]. This functional is also a cornerstone to evaluate the probability of given threshold values, that is P{X ≤ x}, P{X > x} or P{x 1 ≤ X ≤ x 2 } where X denotes a weather variable. Two very interesting threshold for the electricity production is the minimum required wind speed for an eolian device to work and the maximum wind speed before the cut-off of the device. Moreover, the cumulative distribution function can be used to compute the quantiles through inversion. This can be performed numerically with classical algorithm such as the Newton-Raphson methods. We would not present examples of such a process nor assessing any performance of this method versus the direct computation of quantiles.

Here we propose to review two different methods for the computation of cumulative distribution functions. The first one is deterministic while the second one is statistical. Both of them rely on very classical technology. Albeit, the choice of certain parameters or forms of the involved functional was performed carefully by the authors in order to deal with the characteristics of weather variables.

Direct Model Output. We present here one of the first method historically introduced to post-process Ensemble Forecast. It is thus worth to take a look at it. Moreover, it is a very simple method which can be implemented directly after the reading of this paragraph. The bad point is that this method is based on a plotting formula which is controversial since the validity and the reliability is not clearly assess mathematically.

Direct model output tries to estimate the quantities P{X ≤ x} using a formula called the Tukey Plotting Position [37, (3.18) pp. 40 and Table 3.2 pp. 41]. The advantage of this formula which derive from the classical democratic voting method is that it correct deficiencies such as the assignment of zero probability to very small quantiles. It arises, for example, when the classical Gumbel formula is used [37, pp. 285].

Let us now explain the Direct model output. The N Ensemble members are regrouped with the value of the x-quantile which is here estimated empirically from the Ensemble members using any classical estimators. Then P{X ≤ x} is simply estimated using the formula

P{X ≤ x} = Rank(x) -1 3 (N + 1) + 1 3
where Rank(x) is the position of x in the context of the Ensemble members.

The problem with such an estimator is that it does not correct any of the issues such as under dispersion. It basically leads to inaccurate and overconfident results without a very large ensemble [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF].

Estimation using Logistic Regression. Among the first methods for the estimation of P{X ≤ x} where X denotes as usual a weather variable, regression-based methods were very popular [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF]. One of the most successful of such methods is certainly Logistic Regression [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF]. Basically, when a weather variable X is measured, we either have {X ≤ x} or {X > x}. In other words, we have the realization of a Bernoulli law. The estimation of P{X ≤ x} thus simply become an approximation of the parameter p of this Bernoulli law. The idea behind the Logistic Regression is to estimate the logit of p with a linear functional of the predictors F = (F 1 , . . . , F N ) which are the values of the N Ensemble members. It follows

Logit(p) = log p 1 -p = b 0 + b 1 F 1 + • • • + b N F N
where the b i 's are the regression's parameters to estimates and the F i 's are the aforementioned predictors. In practice, the estimation of the b i 's is performed using the maximum likelihood estimator with

Λ(F 1 , . . . , F N ) = M i=1 1(R i ≤ q) exp(b 0 + b 1 F i 1 + • • • + b N F i N ) + (1 -1(R i ≤ q)) 1 + exp(b 0 + b 1 F i 1 + • • • + b N F i N )
where the (R i ) 1≤i≤M are observed outcomes of X and (F i j ) 1≤i≤M, 1≤j≤N are historical Ensemble members. Here Λ(b 1 , . . . , b N ) is the joint likelihood function for the multiple occurrence of the Bernoulli experiments with p is replaced using the Logit formula by

exp(b 0 + b 1 F 1 + • • • + b N F N ) 1 + exp(b 0 + b 1 F 1 + • • • + b N F N )
In [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF], the authors propose a two-predictors logistic regression applied to daily temperature, medium-range temperature and precipitation forecasts. They choose to take the Ensemble's mean F and the Ensemble standard deviation S as predictors. With such a setting,

P{X ≤ x} = exp(α 0 (x) + α 1 (x) F + α 2 (x) F S) 1 + exp(α 0 (x) + α 1 (x) F + α 2 (x) F S)
where α 0 (x), α 1 (x) and α 2 (x) are the parameters that must be estimated for each desired x. The authors called this method LR(2) and compared it to a logistic regression with a single predictor as well as other methods such as Nonhomogeneous Gaussian Regression [START_REF] Gneiting | Calibrated Probabilistic Forecasting Using Ensemble Model Output Statistics and Minimum CRPS Estimation[END_REF]. However, we do not know the algorithm used to performed the Maximum Likelihood Estimation. The experiments can therefore not be reproduce exactly. The comparison shows that LR(2) performed quite well particularly in terms of rank probability score [38, Fig. 1 and2] and Brier score [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF]Fig. 3]. It also points out that LR(2) is a serious challenger to the Nonhomogeneous Gaussian Regression. Overall, the authors acknowledge that their method must be tested on the different weather variables before drawing completely reliable conclusions.

For the electricity market, the method we presented here is very relevant but actually not for direct weather variables. In fact, it was assessed to be a very powerful method for two categorical variables involved in the analysis of the impact of renewable energy generation on electricity markets and which depends on weather forecasts [START_REF] Morales | Integrating renewables in electricity markets -Operational problems[END_REF].

Estimations of Probability Density Functions

Probability density functions are certainly the most well-known probabilistic functionals. It provides the complete repartition of the mass of a random variable and provide thus more information than quantiles or cumulative distribution functions. Moreover, it has a strong physical meaning. Indeed, as we explain above when we describe the models of numerical weather prediction, a probability density function is the natural output of a Monte Carlo methods designed to numerically solve a model.

For the electricity marker, the probability density functions are of major interest since they are involved in the design of optimal offering strategies or in the optimal quantification of reserve requirements [START_REF] Morales | Integrating renewables in electricity markets -Operational problems[END_REF]. The problem of the two-price market can also be solved using probability density functions and the problem of the one-price market requires it.

Here we present a few methods of statistical post-processing of an Ensemble Forecast which provide a probability density function as an output.

Gaussian Kernel Dressing. In [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF], the authors compare their method to compute cumulative distribution function to the integration of a density estimated using a kernel estimator combined with Gaussian ensemble dressing [START_REF] Roulston | Combining dynamical and statistical ensembles[END_REF][START_REF] Wang | Improvement of ensemble reliability with a new dressing kernel[END_REF]. The latter method is often called Gaussian kernel dressing [START_REF] Bröcker | From ensemble forecasts to predictive distribution functions[END_REF]. In few words, this method estimate the conditional density f σ (x|F) of E[X|F] where X is a weather variable and F = (F 1 , . . . , F N ) is an Ensemble forecast of X with N members using [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]Chap. 2.4]

f σ (x|F) := 1 σN N i=1 G x -F i σ
where G denotes a Gaussian and σ is a bandwidth estimator determined by a method of Gaussian ensemble dressing. The density f σ (v|X) is basically estimated by a sequence of bumps, each of which is centered at X i and has a length proportional to σ.

Originally, the methods of Gaussian ensemble dressing were designed to cheaply extend the number of members in an Ensemble Forecast and also remove the bias/errors in the members using historical pairs of forecasts-validations [START_REF] Roulston | Combining dynamical and statistical ensembles[END_REF]. Here we examine the method of Gaussian ensemble dressing developed in [START_REF] Wang | Improvement of ensemble reliability with a new dressing kernel[END_REF] and combined with kernel dressing in [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF]. In [START_REF] Wang | Improvement of ensemble reliability with a new dressing kernel[END_REF], the authors use the second moment to dress the Ensemble Forecast. This choice is made to priorly treat the problem of underdispersivity of Ensemble Forecast. It also overcomes some of the limitations of the dressing method developed in [START_REF] Roulston | Combining dynamical and statistical ensembles[END_REF]. In a nutshell, the authors assume each ensemble members to be of the form

F i = Fi + F i
where Fi denotes the mean of X and F i follows the recentered distribution of V . Hence, they consider that a dressed ensemble member has to be

Ψ i = Fi + F i + ε
where ε is Gaussian independent of (F i ) 1≤i≤N . As a matter of fact, Ψ ∼ N ( F , σ) with a σ which is computed using the second moment constraint. According to the computation performed in appendix of [START_REF] Wang | Improvement of ensemble reliability with a new dressing kernel[END_REF] and aggregated in the formula (9a),

σ = M k=1 ( F k -R k ) 2 -(1 + 1 M ) M k=1 σ k
where M denotes the number of pairs of available historical data belonging to the same season, (R k ) 1≤k≤M are the observed outcomes and F k 1≤k≤M and σ k 1≤k≤M are respectively the empirical means and variance of the historical Ensemble Forecasts.

Part of the experiments on this methods were already commented right above. The weather variables use in those test are still daily temperature, medium-range temperature and precipitation. In term of results, the Brier scores of the Gaussian kernel dressing is given in [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF]Fig. 6]. The most interesting part is certainly the good rank probability scores on a short training period obtained by Gaussian kernel dressing [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF]Table 1]. Overall, it appears to be a very appealing methods in the case of a lack of historical data.

A method based on Gaussian kernel dressing were also developed in [START_REF] Pinson | Ensemble-based probabilistic forecasting at Horns[END_REF] to forecast the wind power using an Ensemble Forecast. The authors of this article model the variance of the gaussian with a Logistic function. The two parameters of the function are assumed to be the same for every kernels. For the estimation of the parameters, the authors build an adaptative estimation using a recursive maximum likehood estimation and the logit of the parameters. The authors compare forecasts produced with their method to raw Ensemble Forecast. Clearly, their Gaussian kernel dressing provide tighter and sharper forecast. However, a comparison of their method to [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF] would have been very interesting. Now this work shows that Gaussian kernel dressing can be used on wind power with success. This is not the case for the different methods to compute densities presented above.

The Screening Regression Procedure. The procedure of multiple linear regression provides an estimate of a weather variable X by relating it to a set of N outputs F 1 , . . . , F n of a model of numerical weather predictions. As an equation,

X = a 0 + a 1 F 1 + • • • + a N F N ,
where the parameters a 0 , . . . , a N minimize the sum of the squares of the estimation errors on the historical data, that is (a 0 , . . . , a N ) = arg min

(b 0 ,...,b N ) M j=1 (R j -b 0 + b 1 F j 1 + • • • + b N F j N ) 2 .
Here R j 1≤j≤M and (F j k ) 1≤j≤M, k=1,...,N represent the historical data of respectively the true realization of X and the outputs F 1 , . . . , F N .

The main problem with multiple linear regression is the selection of the predictors. Indeed, too few or to many of them can provide an estimation which performed poorly. In [START_REF] Glahn | The Use of Model Output Statistics (MOS) in Objective Weather Forecasting[END_REF], the authors choose a method of selection called the screening procedure. First, the variable which is the most highly correlated to V is captured. Then, the others predictors are selected iteratively regarding their contribution to the variance reduction, that is their ability to minimize

1 n M j=1 (R j -R) 2 -1 n M j=1 (R j -F j ) 2 1 n M j=1 (R j -R) 2 (1)
where R is the empirical mean of the true realizations R j 1≤j≤M of X and Fj are the estimate of X using the procedure of multiple linear regression with the predictors F 1 , . . . , F N and the historical predictions F j 1 , . . . , F j N , j = 1, . . . , M . The multiple linear regression with the screening procedure was experimented in [START_REF] Glahn | The Use of Model Output Statistics (MOS) in Objective Weather Forecasting[END_REF] on various weather variables such as surface wind, maximum temperature or cloud amount. The point of the tests is to prove the efficiency of the post-processing of numerical weather predictions. For the probability of precipitation, the authors provides a comparison of the Brier scores of their method versus basic forecasts at local station [10, Fig. 2]. On average, the multiple linear regression performed better. This fact is confirmed on maximum temperature when the method is compared to the so-called Klein-Lewis [START_REF] Glahn | The Use of Model Output Statistics (MOS) in Objective Weather Forecasting[END_REF]Fig. 4]. Finally, we present this method in order to provide one which deals with Numerical Weather Predictions only, not Ensemble Forecast. In addition, the article which introduced this method provide good guidance on the selection of variables for weather predictions.

Nonhomogeneous Regression Techniques. Nonhomogeneous Regression Methods has received several attention as a method to post-process Ensemble Forecasts. Multiple variants of the basic concepts exposed in [START_REF] Gneiting | Calibrated Probabilistic Forecasting Using Ensemble Model Output Statistics and Minimum CRPS Estimation[END_REF] and called EMOS were developed [START_REF] Lerch | Comparison of nonhomogeneous regression models for probabilistic wind speed forecasting[END_REF]. The classical regression method where a weather variable X is approximated by

X = a + b 1 F 1 + • • • + b N F N + ε
where ε is assumed to be Gaussian and F 1 , . . . , F N are the predictors, suggests to take has predictive probability density function a Gaussian whose mean and variance respectively equal the regression estimate and the mean of the square prediction errors. This is a basic estimation by the maximum likelihood estimator of the regression coefficients. In [START_REF] Gneiting | Calibrated Probabilistic Forecasting Using Ensemble Model Output Statistics and Minimum CRPS Estimation[END_REF], the authors modify the above regression equation by assuming

Var(ε) = c + dS 2
where S is the empirical variance of the Ensemble Forecast and c and d has to be computed. They also decide to replace the classical estimation of the coefficients a, b 1 , . . . , b N , c and d by the maximum likelihood estimator. They choose to use take the parameters a, b 1 , . . . , b N , c and d which minimize the CRPS which can be here computed analytically, that is

(a, b 1 , . . . , b N , c, d) = arg min (A,B 1 ,...,B N ,C,D) 1 M M k=1 (C+DS 2 k ) 1 2 Z k [2Φ(Z k ) -1] + 2φ(Z k ) - 1 √ π where Z k = R k -(A + B 1 F k 1 + • • • + B N F k N ) (C + DS 2 k ) 1 2
.

Here Φ and φ respectively denote the cumulative distribution function and the probability density function of a standard Gaussian, (F k , F k 1 , . . . , F k N , S k ), k = 1, ..., M stand for the historical data of realizations, Ensemble members and empirical variance of the Ensemble members. . The estimation of the parameters is performed numerically by the authors using the Broyden-Fletcher-Goldfarb-Shanno algorithm which is implemented in the R language.

The method is tested using the University of Washington ensemble over the Pacific Northwest to forecast the sea level and the surface temperature. The basic scoring rules such as the Brier score and the CRPS are used to assess the performance of the method. The author also compare forecasts to raw ensemble, de-biased ensemble and climatological forecasts. The regression method provided by the authors clearly outperformed the climatological forecast. It more lightly beats the bias-corrected Ensemble Forecast [13, Table 4]. Some comparisons with other statistical methods are given in [START_REF] Wilks | Comparison of Ensemble-MOS Methods Using GFS Reforecasts[END_REF] using the rank probability score and the Brier score. In term of applications, this method was proven to be very relevant to forecast low and high wind regime when the kernel is respectively replaced by a truncated Gaussian kernel and a so-called GEV distribution [START_REF] Lerch | Comparison of nonhomogeneous regression models for probabilistic wind speed forecasting[END_REF]. Nonhomogeneous Regression can be easily combined with regime-switching. Hence, multiple wind regimes can be thus treated with the right kernel. This is clearly a great strength and thus a good encouragement to use this method on wind farm locate on places with irregular wind regimes.

Estimation using Bayesian Model Averaging. In [START_REF] Raftery | Using Bayesian Model Averaging to Calibrate Forecast Ensembles[END_REF], the authors propose to estimate the conditional density E[X|F] using a simplification of Bayesian model averaging. As usual, X and F = (F 1 , . . . , F N ) respectively denotes a weather variable and a vector of N Ensemble members [START_REF] Hoeting | Bayesian Model Averaging: A Tutorial[END_REF]. The method can be seen as an improved kernel dressing. Boiled down to its essence, the idea is to estimate the conditional density f (x|F) as a mixture, that is

f (v|F) = N k=1 w k g k (x|F k )
Here, for each k = 1, ...N , g k (x|F k ) is a conditional density which represents the information provided by F k and w k is a weight of the contribution of F k . The weights w 1 , . . . , w n can be understood as the probability that F k is the member of the ensemble which provides the most accurate forecast. In classical Bayesian model averaging, w k is replaced by P(g k (x|F k )|X D ) where X D represents the historical data, that is the posterior probability of the estimation of X using the kernel g k (x|F k ). For the estimation of temperature and sea level, the authors choose to take E[X|F k ] ∼ N (a k + b k F k , σ 2 ). They directly estimates a k and b k where k = 1, . . . , N using linear regression with paired historical data of true realizations and Ensemble Forecasts. It can be considered as a bias-correction process.

parameter σ and the weight are, on the other side, estimated using the maximum likelihood estimator and the same pairs of historical data. It has actually become usual to replace classical Bayesian estimation procedures by the maximum likelihood estimator in Bayesian model averaging to avoid expensive computations. In some instances, it provides a justified large-sample approximation [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF].

The authors performed various benchmarks of their method where the EM algorithm is use to estimate the maximum likelihood. They compare their method to various raw Ensemble Forecasts using different measurement tools such as the CRPS or the RMSE. Overall, it is showed that the method performed very well [28, Table 4 to 6 and 8 to 12]. The authors also investigated the possibility to reduce the number of ensemble members by disqualifying those with a poor weight [START_REF] Raftery | Using Bayesian Model Averaging to Calibrate Forecast Ensembles[END_REF]Table 7]. Due to its performance, the present method was adapted to the prediction of wind speed [START_REF] Sloughter | Probabilistic Wind Speed Forecasting Using Ensembles and Bayesian Model Averaging[END_REF]. In this article, several modification were performed since the kernels are chosen to be Gamma distributed. Different techniques of the estimation of the maximum likelihood were tried. However, it appears that the different methods perform similarly [START_REF] Sloughter | Probabilistic Wind Speed Forecasting Using Ensembles and Bayesian Model Averaging[END_REF]Figure 4]. Some benchmarks which compared the method to basic Ensemble Forecast and climatology are also provided. The result are similar [START_REF] Sloughter | Probabilistic Wind Speed Forecasting Using Ensembles and Bayesian Model Averaging[END_REF]Table 2] to those obtained in [START_REF] Raftery | Using Bayesian Model Averaging to Calibrate Forecast Ensembles[END_REF].

Bayesian Model Averaging is actually one of the leading method for the post-processing of Ensemble Forecast. It really takes into account the model uncertainty arising in Ensemble Forecast. However, several non-trivial improvements has to be performed in order to used this method to treat problems linked to the electricity market [START_REF] Sloughter | Probabilistic Wind Vector Forecasting using Ensembles and Bayesian Model Averaging[END_REF].

  using a local linear kernel weighted version of E[ρ θ (E p -a)|X = x] [30, 35]. Here F i and e p,i i = 1, . . . , M are respectively M historical values of the predictors F and the wind power E p , ρ θ and w are respectively the "check" function and the tricube weight function given by

Acknowledgement. This work was supported and funded by the Agence National pour la Recherche through the research project CAESARS. We thank E. Gobet and C. Alasseur for their various comments and their help to improve this paper.