- Email: [email protected]
GEFCom2014 probabilistic electric load forecasting using time series and semi-parametric regression models
International Journal of Forecasting (
)
–
Contents lists available at ScienceDirect
International Journal of Forecasting journal homepage: www.elsevier.com/locate/ijforecast
GEFCom2014 probabilistic electric load forecasting using time series and semi-parametric regression models V. Dordonnat, A. Pichavant, A. Pierrot ∗ EDF R&D, Clamart, France
article
info
Keywords: Electric load forecasting Semi-parametric regression Temperature simulation Probabilistic forecasting
abstract This paper describes the forecasting methodology used by the team ADADA in the load forecasting track of GEFCom2014. The approach includes a semi-parametric regression model for point load forecasting, a multivariate time series simulation model for temperature forecasting, and a resampling strategy for building probabilistic forecasts from point load forecast residuals. © 2015 International Institute of Forecasters. Published by Elsevier B.V. All rights reserved.
1. Introduction Load forecasting is a crucial activity for various utility applications, such as unit commitment, market operations, and the development of investment strategies. In recent years, most of the research work has been dedicated to determining the most accurate point forecasting model, but the interest in building a valid probabilistic forecasting model is increasing. Indeed, the growth in the proportion of the electricity generation that is provided by renewables (wind and solar power), with the associated uncertainty, means that more flexibility is required in demand-side management and production from thermal units. There are several different ways of producing probabilistic forecasts: building a quantile model directly, integrating weather information (ensemble forecasts in the short run, historical data or a weather generator for longer forecasting horizons), or modeling the residuals from point forecasts. This paper describes the methodology adopted by the ADADA team in the load forecasting track of GEFCom2014. The aim was to produce one-month-ahead probabilistic forecasts (as the 1%–99% quantiles of a forecast distribution) during 15 weeks (including three trial weeks which did not
∗
Corresponding author. E-mail address: [email protected] (A. Pierrot).
count towards the final score). For a complete description of GEFCom2014, see Hong et al. (this issue). The team strategy was to build a point forecasting model and simulate the uncertainty in both temperature and point load forecast errors. The paper is organized as follows. Section 2 describes the dataset and the main variables in it. Section 3 describes our forecasting methodology for both temperature and load. Our load model relies on generalized additive models (GAM), which have been proven to be effective for the short-term forecasting of French load, as was showed by Pierrot and Goude (2011). Section 4 provides our forecasting results. Section 5 concludes the paper. 2. Data description Yt ,h denotes the electric load for day t and hour h, from 1:00, January 1st, 2005, until 00:00, October 1st, 2010. Historical temperature data (in ◦ F) are also used, starting on January 1st 1:00, 2001; 25 stations are provided, denoted Tt ,h,1 , . . . , Tt ,h,25 .1 Both datasets are extended with one month of data after each week of the competition. A calendar of US Federal holidays is also used, as authorized by GEFCom2014 organizers.
1 The indexes t and h may not be written, so as to improve the readability.
http://dx.doi.org/10.1016/j.ijforecast.2015.11.010 0169-2070/© 2015 International Institute of Forecasters. Published by Elsevier B.V. All rights reserved.
2
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Fig. 1. Time series of the daily mean load, from 2005 to 2010. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. Observed load vs. the observed mean of temperatures, Tt ,h,6 , Tt ,h,10 and Tt ,h,13 , from 2005 to 2010 (left), week vs. weekend profiles (center), and winter vs. summer profiles (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The following subsections present the main features of the load and temperature data, some of which are considered in our load forecasting methodology. 2.1. Description of load data The load data exhibit a clear yearly pattern, with higher consumption levels during the summer and winter months (see Fig. 1). We also notice in Fig. 1 that the load fluctuates more in winter than in summer. This yearly pattern can be explained by the strong dependence of the load upon temperature, with both heating and cooling effects, as the left panel of Fig. 2 shows. The mean daily profiles of the load for weekdays and weekends are given in the center graph of Fig. 2. The evening peak is larger for weekdays, while the morning peak is later but larger on weekends. The daily patterns in summer and winter are different, as the right panel of Fig. 2 shows. In winter, the two peaks at 8AM and 8PM may be explained by lighting needs.
variables. The observations below are based on temperature Tt ,h,6 , but also apply to all other temperature variables. The left panel of Fig. 3 shows the estimated daily normal temperature for Tt ,h,6 . The yearly pattern is clearly visible, with the mean temperature in the summer twice that for winter. The trend on this relatively short dataset has been estimated to be zero. The center panel of Fig. 3 shows the average daily pattern of Tt ,h,6 calculated over the whole year (black line) or by month (colored lines). The minimum temperature is reached at around 6AM, and the maximum at around 3–4PM. The right panel of Fig. 3 shows the yearly pattern of the deviation from the estimated normal temperature. The winter variability is clearly higher than the mean level, whereas the summer variability is the lowest. Next, we propose a multivariate forecasting methodology for dealing with all weather variables simultaneously, so that any combination of temperatures that improves the modelling of the weather–load relationship can be used for probabilistic forecasting. 3. Methodology
2.2. Description of temperature data As was shown in the previous section, there is a clear dependence between load and temperature. Thus, it is necessary to build a proper forecasting model for temperature
We took the following steps in order to derive the forecasting quantiles of electricity demand:
• build a temperature-based deterministic load model, as described in Section 3.1.
V. Dordonnat et al. / International Journal of Forecasting (
)
–
3
Fig. 3. Estimated normal temperature Tt ,h,6 (left), daily profiles (center), and the yearly pattern of the deviation from the normal temperature (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
• build a simulation model for temperature, as detailed in Section 3.2, and generate sample paths for each month to be forecast (for example, N ∗ = 1000). • use the temperature sample paths as inputs for the deterministic load model, to obtain load sample paths. • combine the results with quantiles based on draws from deterministic load model errors on the estimation dataset in order to include the model uncertainty, as described in Section 3.3; N final paths are obtained. • derive quantiles from the N final paths for each time step. 3.1. Deterministic load model Statistical framework. The deterministic load model is built in the generalized additive model (GAM) framework, see Pierrot and Goude (2011) and Wood (2006). The general structure is as follows: Yi = Xi θ + f1 (x1i ) + f2 (x2i ) + · · ·
+ fp (xpi ) + εi ,
∀ i = 1, . . . , n,
(1)
where
• Yi , the response variable (here, electric load), belongs to • •
•
•
some exponential family distribution (here, Gaussian), E[εi ] = 0, V[εi ] = σ 2 and cov(εi , εj ) = 0, ∀ i ̸= j, Xi is a row of the model matrix for any strictly parametric model components (the day type effect, for instance) and θ is the parameter vector of the linear part, (x1,i , x2,i , . . . , xp,i ) are the covariates of the nonparametric part (such as meteorological variables), with the associated non-linear functions (f1 , f2 , . . . , fp ) being assumed to be smooth, and n is the number of observations.
The unknown parameters θ and functions fj are estimated by penalized regression methods. We control the model’s smoothness by adding a ‘‘wiggliness’’ penalty to the least squares fitting objective. The model is fitted by minimizing:
2 1 p p fj (xji ) + λj [fj′′ (x)]2 dx. Y − X θ − 0 j =1 j=1
(2)
The compromise between model fitting and the smoothness of the functions is controlled using λ: if λ is too high, there is a risk of oversmoothing, but if λ is too small, there is a risk of overfitting. Following (Wood, 2006), the GCV (generalized cross-validation) score is used to choose λ: GCV (λ) =
n 1 (yi − yi )2
n i=1 tr (I − A)2
,
(3)
where yi is the prediction of yi obtained by fitting to all data, and A is the hat matrix for the model fitted to all data: Y = AY . In practice, the models have been estimated using the mgc v package supplied in R, see Wood (2013). Weather station selection. Based on the data analysis in Section 2, we aimed to build a model that included the temperature as an explanatory variable. As no geographical information was provided, we had to determine how to associate one or more weather stations with the load. We first used the GCV score as a criterion for choosing the best temperature for our model, with lower values being better. We chose to fit hourly GAMs in order to model the strong variation in the effects of the load, depending on the current hour. The hourly model used for temperature selection is as follows: Yt ,h = fh (posant ) + θDTk,h .DayTypek,t + θh .t
+ gh (Tt ,h,s ) + εt ,h ,
t = 1, . . . , n, h = 0, . . . , 23,
(4)
where Y is the load, posan is the day of the year (from 0 for January 1st to 1 for December 31st), fh is a nonlinear function requiring estimation, DayType is the day type (k = 2 for week days, k = 1 for weekends, and hence θDTk takes only two values), t is the trend, and gh is a nonlinear function for estimating the temperature at station s. Using Eq. (4), we built 25 models, one for each of the 25 available temperatures. As our models are hourly models, we obtained 24 GCV scores for each temperature variable. Fig. 4 shows the lowest GCV scores. We observed that the models’ performances for the deterministic and probabilistic forecasts were quite different; thus, the GCV criterion could not be used, as it consists of choosing the model with the lowest GCV score. We then tried to calculate a weighted average. We used an EWA algorithm (exponentially weighted average, see
4
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Tt ,i denotes the daily mean temperature for day t at station s, MAT t ,s is the corresponding moving-average (the window is set to 31 days), and Td∗,s is the historical mean temperature for day-of-the-year d at station s. Observations on February 29th are removed from the dataset. We have MAT t ,s =
15 1 T t − j ,s , 31 j=−15
t = 1, . . . , n, s = 1, . . . , 25 Td∗,s =
Fig. 4. The lowest GCV scores. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(Cesa-Bianchi & Lugosi, 2006)), which makes predictions by combining a finite set of provided forecasts. We used this algorithm on the 25 deterministic GAM forecasts (one for each weather station). We denote the deterministic forecast for the temperature Ts by F(Ts ). The EWA algorithm converges to the selection of three models with equal weights: 1 3
F (T6 ) +
1 3
F (T10 ) +
1 3
N t =1
MAT t ,s δd(t )=d ,
(6)
d = 1, . . . , 365,
(7)
where d(·) gives the day of the year, N is the number of years in our sample, and n is the total number of observations. Note that we first consider all years to have ∗ length 365, and take T60 , which is February 29th , to be a simple mean of the preceding and following days, in order to obtain normal temperatures on 366 days. After estimating these reference temperature variables, we model the deviation from them using autoregressive (AR) models. Let 1Tt ,s be the deviation from the normal temperature on day t at station s:
1Tt ,s = Tt ,s − Td∗(t ),s 1Tt ,s =
F (T13 ).
Ps
ρp,s 1Tt −p,s + εt ,s ,
(8) s = 1, . . . , 25.
(9)
p=1
We note that the three temperatures selected by EWA had already been among the six best temperatures selected by the GCV score. We therefore define a new temperature variable to be the mean temperature of stations 6, 10 and 13: Tt ,h,26 = 31 Tt ,h,6 + 31 Tt ,h,10 + 31 Tt ,h,13 . This new temperature variable improved both the deterministic and probabilistic forecasts. Load model. Using the new temperature variable T26 , we derived other potential explanatory variables. Of all of the models we have tested, the model that presents the best average score over the year 2011 is Yt ,h = fh (posant ) + θDTk,h · DayTypek,t + θh .t
+ gh (Tt ,h,26 ) + kh (STt ,h,26 ) + εt ,h , t = 1, . . . , n, h = 0, . . . , 23,
n 1
(5)
where STt ,h,26 is an exponentially smoothed version of Tt ,h,26 (with the smoothing coefficient α being estimated by cross-validation), and fh , gh , kh are non-linear functions that require estimation. Other models have also been tested during the competition, for example using different exponential smoothing parameters, the daily maximum temperature, and updating forecasts based on the consumption of the previous day (referred to as update1day in Table 2). The results are shown in Section 4. 3.2. Temperature simulation model A multivariate simulation model for the 25 weather stations is now set. First, the daily normal temperature for all weather stations is estimated using a moving average.
The order of the AR Ps can differ between variables, and the parameters are estimated via maximum likelihood. The residuals from the AR models are white noise, but do not need to be Gaussian. Using this method, we were able to simulate N ∗ paths for each weather variable:
• For each day to be forecast, T + 1, . . . , T + K (K depends on the month to be forecast), draw N ∗ vectors from εˆ t = (ˆεt ,1 , . . . , εˆ t ,25 ), denoted (εT(k+) 1 , . . . , εT(k+) k )k=1,...,N ∗ .
For a specific month to be forecast, we use only historical residuals from the same month and the preceding and following ones for sampling. • Simulate deviations from the normal temperature using Eq. (9). • Add corresponding normal temperatures (Td∗(T +k),s )s=1,...,25; k=1,...,K so as to obtain simulations of the daily weather variables. • To get to the hourly level, apply a daily profile. This profile can be either constant or month-dependent, as is shown in the center panel of Fig. 3. Tt ,h,26 can be derived easily from this simulation step. 3.3. Probabilistic forecasts We used the simulation model presented in Section 3.2 to produce 1000 paths for the input temperature of the load model in Eq. (5). Thus, we obtained 1000 paths for forecasting the load of the month. The number of paths is set arbitrarily. Once we had taken into account the temperature uncertainty, we wanted to model the uncertainty for our load forecasting. One way to do this was first to take a large
V. Dordonnat et al. / International Journal of Forecasting (
)
–
5
Fig. 5. Raw temperature (left) and smoothed temperature (right) effects for model (5).
Fig. 6. For T26 residuals of the AR model, the left panel of the figure shows the empirical ACF for the N observations based on kernel-based density estimation (black) and Gaussian distribution fit (red), while the center and right panels show the specified bandwidth and the time series, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
sample (e.g., 1 million) for each time-step to be forecast, by sampling with replacement from the fitting residuals εˆ t ,h of our load model in Eq. (5). From this sample, we then derived sequences of quantiles 1%–99%. The temperature and load uncertainties are combined by adding each sequence of quantiles of the load uncertainty to each simulated load path, based on a temperature path. We therefore obtained 99 000 final load paths from which to compute the final 99 forecasting quantiles. We tested several other ways of modeling the load uncertainty. Instead of taking a sample of the absolute fitting residuals (referred to as loaduncert1 in Table 2), we took a sample of the relative fitting residuals (referred to as loaduncert2 in Table 2). We also tried unsuccessfully to take fitting residuals by instant or by month, one-month forecasting residuals and residuals from 1- to 31-dayahead forecasting. 4. Results 4.1. Estimation results Load model. Fig. 5 presents the estimated temperaturebased effects of the load model in Eq. (5): the effect of
the raw temperature on the left is larger than that of the smoothed temperature on the right. The shape of the raw temperature effect is similar for all hours of the day, while the smoothed temperature effect is smaller during business hours (possibly due to residential heating). The smoothed temperature for cooling remains large, possibly due to cooling in offices. Temperature model. We present the simulation results for the temperature variable T26 that was defined earlier. According to Fig. 6 (left) and to the Ljung–Box test, it seems that there is no autocorrelation left in the AR model residuals. However, the distribution of the residuals is too spiked to be Gaussian. This motivated us to draw from the empirical residuals in order to simulate the temperature. As we can see from the right panel of Fig. 6, there is a clear yearly seasonality in the residuals. This is dealt with in the simulation step when we draw from the residuals. We only use residuals from the current month and the preceding and following months for forecasting. An alternative would have been to normalize the data using a moving average on squared data. This has been tried, but was no more efficient in our tests. The daily profile used for the simulation is constant over the months.
6
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Table 1 Extract of forecast MAPE values. Model
Jan. 11
Feb. 11
March 11
Apr. 11
May 11
June 11
July 11
Aug. 11
Sept. 11
Oct. 11
Nov. 11
Mean
T6 T6_Tmax6 T6_ST6_0.992 T6_ST6_0.992_Tmax6 T6_ST6_0.996 T6_ST6_0.996_Tmax6 T6_ST6_ST10_0.992 T6_ST6_ST10_0.996 T26 T26_Tmax26 T26_ST26_0.992 T26_ST26_0.992_Tmax26 T26_ST26_0.996 T26_ST26_0.996_Tmax26
5.79 5.22 6.52 6.10 6.58 6.15 6.02 7.32 4.56 4.23 4.94 4.74 5.51 5.34
8.86 8.56 6.85 6.89 7.62 7.49 6.70 8.11 6.94 6.77 5.23 5.26 5.65 5.66
9.41 8.91 7.60 7.36 7.85 7.50 7.66 8.16 7.78 7.35 5.80 5.54 6.00 5.68
8.78 8.27 7.43 6.97 7.42 7.06 7.13 8.07 8.06 7.35 6.68 6.12 7.03 6.51
10.14 10.18 10.17 10.14 9.59 9.65 7.39 9.97 6.78 6.36 5.73 5.43 5.63 5.35
10.65 11.71 14.09 14.69 14.15 14.96 10.02 14.53 6.67 6.63 8.63 8.38 8.82 8.55
5.93 5.78 6.15 5.91 5.98 5.72 5.28 6.06 5.38 4.97 5.42 5.01 5.46 4.99
35.67 35.33 35.64 35.07 35.13 34.48 33.43 34.18 32.84 32.07 32.56 31.76 32.58 31.68
14.63 14.43 13.86 14.03 14.26 14.28 10.74 13.40 12.86 12.30 11.23 11.06 11.52 11.23
10.41 9.03 8.06 7.26 8.64 7.70 7.42 8.03 10.06 9.50 6.72 6.61 6.98 6.84
9.82 8.59 7.95 7.34 8.49 7.70 8.24 8.23 8.03 7.41 6.50 6.27 6.68 6.38
11.83 11.46 11.30 11.07 11.43 11.15 10.00 11.46 10.00 9.54 9.04 8.74 9.26 8.93
Note: The best forecast is shown in italics.
Table 2 Extract of forecasting results. Model
Task 4
Task 5
Task 6
Task 7
Task 8
Task 9
T6_T25_Tmax25_loaduncert1 T6_Tmax6_loaduncert1 T6_Tmax6_loaduncert2 T6_ST6_0.996_loaduncert1 T6_ST6_0.996_Tmax6_loaduncert1 T6_ST6_0.996_Tmax6_loaduncert2 T6_ST6_ST10_0.996_loaduncert2 T6_ST6_STmin6_0.996_loaduncert1 T6_ST6_STmin6_0.996_loaduncert2 T_PCA_loaduncert1 T26_ST26_0.992_Tmax26_loaduncert2 T26_ST26_0.996_loaduncert2 T26_ST26_0.996_Tmax26_loaduncert1 T26_ST26_0.996_Tmax26_loaduncert2 T26_ST26_0.996_Tmax26_update1day_loaduncert1
10.51 10.00 10.06 10.17 9.82 9.86 10.70 10.09 10.12
10.08 10.04 9.31 9.65 9.61 9.03 10.16 10.11
7.88 7.88 7.62 7.69 7.67 7.56 7.84 7.83
5.51 5.50 5.60 5.56 5.54 5.42 5.96 5.92
6.47 7.86 6.52 7.82 6.63 8.32 6.62 8.01 6.69 7.97 7.00 8.93 10.25 12.44 10.44 12.42
10.49 10.61 10.28 10.32
9.39 9.25 9.44 9.42
7.70 7.65 7.66 7.66
4.83 4.77 4.68 4.65 4.58 4.56 4.83 4.75 4.73 4.43 4.41 4.46 4.41
5.36 5.40 5.39 5.37
6.92 6.89 6.86 6.90
Task 10
8.55 8.71 8.67 8.63
Task 11
Task 12
Task 13
Task 14
Task 15
Mean
11.14 11.03 11.01 11.01 10.92 11.21 12.70 12.72
6.62 6.57 6.03 6.30 6.26 5.40 7.51 7.47
4.52 4.46 4.05 3.99 3.91 3.74 4.47 4.44
7.18 7.16 6.81 7.03 7.01 6.73 6.16 6.16
9.45 9.38 9.28 9.67 9.59 9.19 8.05 8.07
7.63 7.60 7.46 7.50 7.47 7.46 8.37 8.37
10.90 10.93 10.90 10.86
5.77 5.50 5.68 5.67 5.78
3.97 3.62 3.71 3.64
6.64 6.60 6.76 6.75
8.89 8.91 9.15 9.11
7.42 7.37 7.41 7.39
Notes: Effective forecasts are in bold, the best forecasts are in italics.
4.2. Forecast evaluation When trying to build an accurate point load forecasting model, the mean absolute percentage error (MAPE) has been used to select the best deterministic load model, by excluding years 2009, 2010 and 2011 from the estimation dataset in turn, and forecasting these years using historical temperatures. Table 1 gives the 2011 results (December results are not provided, due to the lack of temperature data). Thus, the estimation dataset was four years long. The pinball loss based quantile score was also used during the competition in order to test different probabilistic models for backcasting on the years 2009, 2010 and 2011. The submissions correspond to different models, with each submission being selected by looking at the best past scores for the same month. A summary of the scores is given in Table 2. Note that the best model when using a MAPE criterion was not necessarily the best one when using simulations and quantile forecasts. Tables 1 and 2 present the results obtained using the different models tested, where the name of the model indicates the variables used to model the weather dependent part. Model T26_ST26_0.992_Tmax26 is estimated with temperature T26 , the smoothed temperature ST26 with the
smoothing coefficient α = 0.992, and the daily maximum T26 temperature T max26 . loaduncert1 and loaduncert2 refer to the different approaches tested for the load model uncertainty. 4.3. Discussion Section 3.1 has presented what was our best deterministic load model on average. Other input variables have been tested as well:
• calendar-based variables: we also considered variables based on the year or month, bank holidays, and different day types. • temperature-based variables: we also considered lagged temperatures (one day, one hour, etc.), the daily maximum, minimum, mean and median, and a PCA (referred to as T_PCA in Table 2) of the 25 temperature variables. We selected variables using a MAPE criterion on validation datasets, by looking at the values of the estimated parameters and the graphs of the estimated splines; we therefore chose the day of the year instead of the month, a simple day type for week/weekend, and a linear trend. Our holiday variable did not provide any
V. Dordonnat et al. / International Journal of Forecasting (
improvement. As the MAPEs on October 2nd and 3rd, 2006, were very high, more than 30%, we removed these days from our fitting dataset. During the competition, we used past monthly scores (deterministic and probabilistic) to select one model for forecast submission. We recalculated the scores of the candidate models when new observations were released, and found that our best forecast was not always the one we had submitted. Since we found no rule for determining the best model to submit, a good way to improve would be by using model averaging. 5. Conclusion The proposed approach allowed the ADADA team to reach the top five of the load forecasting track. A deterministic load model has been developed, accounting for the day of the week as well as for the influence of weather. For providing probabilistic forecasts, we have used a resample of historical point forecast errors, together with a simulation model for temperature. Building a good probabilistic forecasting model for the electric load is a difficult task, especially when including the uncertainty on weather variables. Indeed, the best point forecasting model will not necessarily be the best when it is used with temperature simulation. It is also difficult to compare different modelling strategies, since only one observation for each hour of the day can be used for evaluating quantile forecasts.
)
–
7
Further ways to improve the modelling approach could be the inclusion of forecast averaging and the use of a better short-term forecasting model for temperature. References Cesa-Bianchi, N., & Lugosi, G. (2006). Prediction, learning, and games. Cambridge University Press. Hong, T., Pinson, P., Fan, S., Zareipour, H., Troccoli, A., & Hyndman, R. J. (2016). Probabilistic energy forecasting: state-of-the-art 2015. International Journal of Forecasting, this issue. Pierrot, A., & Goude, Y. (2011). Short-term electricity load forecasting with generalized additive models. In Proceedings of the 16th Intelligent System Applications to Power Systems Conference, ISAP 2011. Wood, S. N. (2006). Generalized additive models, an introduction with R. Chapman and Hall. Wood, S.N. (2013). Package mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness estimation.
V. Dordonnat holds an M.Sc. degree in statistics from the University of Versailles-Saint Quentin, France, and a Ph.D. degree from the VU University, Amsterdam, The Netherlands. From 2009 to 2015, she worked as a research engineer at EDF R&D, France. Since April 2015, she has been working at RTE, the French TSO, as a data scientist.
A. Pichavant holds a Master’s-level engineering degree from ENSAI, Rennes, France. She has been working at EDF R&D as a data scientist since 2012, in the load forecasting team.
A. Pierrot holds an M.Sc. degree in statistics from the University of Paris Dauphine, France. She has been working at EDF R&D as a data scientist since 2008, in the load forecasting team.
)
–
Contents lists available at ScienceDirect
International Journal of Forecasting journal homepage: www.elsevier.com/locate/ijforecast
GEFCom2014 probabilistic electric load forecasting using time series and semi-parametric regression models V. Dordonnat, A. Pichavant, A. Pierrot ∗ EDF R&D, Clamart, France
article
info
Keywords: Electric load forecasting Semi-parametric regression Temperature simulation Probabilistic forecasting
abstract This paper describes the forecasting methodology used by the team ADADA in the load forecasting track of GEFCom2014. The approach includes a semi-parametric regression model for point load forecasting, a multivariate time series simulation model for temperature forecasting, and a resampling strategy for building probabilistic forecasts from point load forecast residuals. © 2015 International Institute of Forecasters. Published by Elsevier B.V. All rights reserved.
1. Introduction Load forecasting is a crucial activity for various utility applications, such as unit commitment, market operations, and the development of investment strategies. In recent years, most of the research work has been dedicated to determining the most accurate point forecasting model, but the interest in building a valid probabilistic forecasting model is increasing. Indeed, the growth in the proportion of the electricity generation that is provided by renewables (wind and solar power), with the associated uncertainty, means that more flexibility is required in demand-side management and production from thermal units. There are several different ways of producing probabilistic forecasts: building a quantile model directly, integrating weather information (ensemble forecasts in the short run, historical data or a weather generator for longer forecasting horizons), or modeling the residuals from point forecasts. This paper describes the methodology adopted by the ADADA team in the load forecasting track of GEFCom2014. The aim was to produce one-month-ahead probabilistic forecasts (as the 1%–99% quantiles of a forecast distribution) during 15 weeks (including three trial weeks which did not
∗
Corresponding author. E-mail address: [email protected] (A. Pierrot).
count towards the final score). For a complete description of GEFCom2014, see Hong et al. (this issue). The team strategy was to build a point forecasting model and simulate the uncertainty in both temperature and point load forecast errors. The paper is organized as follows. Section 2 describes the dataset and the main variables in it. Section 3 describes our forecasting methodology for both temperature and load. Our load model relies on generalized additive models (GAM), which have been proven to be effective for the short-term forecasting of French load, as was showed by Pierrot and Goude (2011). Section 4 provides our forecasting results. Section 5 concludes the paper. 2. Data description Yt ,h denotes the electric load for day t and hour h, from 1:00, January 1st, 2005, until 00:00, October 1st, 2010. Historical temperature data (in ◦ F) are also used, starting on January 1st 1:00, 2001; 25 stations are provided, denoted Tt ,h,1 , . . . , Tt ,h,25 .1 Both datasets are extended with one month of data after each week of the competition. A calendar of US Federal holidays is also used, as authorized by GEFCom2014 organizers.
1 The indexes t and h may not be written, so as to improve the readability.
http://dx.doi.org/10.1016/j.ijforecast.2015.11.010 0169-2070/© 2015 International Institute of Forecasters. Published by Elsevier B.V. All rights reserved.
2
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Fig. 1. Time series of the daily mean load, from 2005 to 2010. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. Observed load vs. the observed mean of temperatures, Tt ,h,6 , Tt ,h,10 and Tt ,h,13 , from 2005 to 2010 (left), week vs. weekend profiles (center), and winter vs. summer profiles (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The following subsections present the main features of the load and temperature data, some of which are considered in our load forecasting methodology. 2.1. Description of load data The load data exhibit a clear yearly pattern, with higher consumption levels during the summer and winter months (see Fig. 1). We also notice in Fig. 1 that the load fluctuates more in winter than in summer. This yearly pattern can be explained by the strong dependence of the load upon temperature, with both heating and cooling effects, as the left panel of Fig. 2 shows. The mean daily profiles of the load for weekdays and weekends are given in the center graph of Fig. 2. The evening peak is larger for weekdays, while the morning peak is later but larger on weekends. The daily patterns in summer and winter are different, as the right panel of Fig. 2 shows. In winter, the two peaks at 8AM and 8PM may be explained by lighting needs.
variables. The observations below are based on temperature Tt ,h,6 , but also apply to all other temperature variables. The left panel of Fig. 3 shows the estimated daily normal temperature for Tt ,h,6 . The yearly pattern is clearly visible, with the mean temperature in the summer twice that for winter. The trend on this relatively short dataset has been estimated to be zero. The center panel of Fig. 3 shows the average daily pattern of Tt ,h,6 calculated over the whole year (black line) or by month (colored lines). The minimum temperature is reached at around 6AM, and the maximum at around 3–4PM. The right panel of Fig. 3 shows the yearly pattern of the deviation from the estimated normal temperature. The winter variability is clearly higher than the mean level, whereas the summer variability is the lowest. Next, we propose a multivariate forecasting methodology for dealing with all weather variables simultaneously, so that any combination of temperatures that improves the modelling of the weather–load relationship can be used for probabilistic forecasting. 3. Methodology
2.2. Description of temperature data As was shown in the previous section, there is a clear dependence between load and temperature. Thus, it is necessary to build a proper forecasting model for temperature
We took the following steps in order to derive the forecasting quantiles of electricity demand:
• build a temperature-based deterministic load model, as described in Section 3.1.
V. Dordonnat et al. / International Journal of Forecasting (
)
–
3
Fig. 3. Estimated normal temperature Tt ,h,6 (left), daily profiles (center), and the yearly pattern of the deviation from the normal temperature (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
• build a simulation model for temperature, as detailed in Section 3.2, and generate sample paths for each month to be forecast (for example, N ∗ = 1000). • use the temperature sample paths as inputs for the deterministic load model, to obtain load sample paths. • combine the results with quantiles based on draws from deterministic load model errors on the estimation dataset in order to include the model uncertainty, as described in Section 3.3; N final paths are obtained. • derive quantiles from the N final paths for each time step. 3.1. Deterministic load model Statistical framework. The deterministic load model is built in the generalized additive model (GAM) framework, see Pierrot and Goude (2011) and Wood (2006). The general structure is as follows: Yi = Xi θ + f1 (x1i ) + f2 (x2i ) + · · ·
+ fp (xpi ) + εi ,
∀ i = 1, . . . , n,
(1)
where
• Yi , the response variable (here, electric load), belongs to • •
•
•
some exponential family distribution (here, Gaussian), E[εi ] = 0, V[εi ] = σ 2 and cov(εi , εj ) = 0, ∀ i ̸= j, Xi is a row of the model matrix for any strictly parametric model components (the day type effect, for instance) and θ is the parameter vector of the linear part, (x1,i , x2,i , . . . , xp,i ) are the covariates of the nonparametric part (such as meteorological variables), with the associated non-linear functions (f1 , f2 , . . . , fp ) being assumed to be smooth, and n is the number of observations.
The unknown parameters θ and functions fj are estimated by penalized regression methods. We control the model’s smoothness by adding a ‘‘wiggliness’’ penalty to the least squares fitting objective. The model is fitted by minimizing:
2 1 p p fj (xji ) + λj [fj′′ (x)]2 dx. Y − X θ − 0 j =1 j=1
(2)
The compromise between model fitting and the smoothness of the functions is controlled using λ: if λ is too high, there is a risk of oversmoothing, but if λ is too small, there is a risk of overfitting. Following (Wood, 2006), the GCV (generalized cross-validation) score is used to choose λ: GCV (λ) =
n 1 (yi − yi )2
n i=1 tr (I − A)2
,
(3)
where yi is the prediction of yi obtained by fitting to all data, and A is the hat matrix for the model fitted to all data: Y = AY . In practice, the models have been estimated using the mgc v package supplied in R, see Wood (2013). Weather station selection. Based on the data analysis in Section 2, we aimed to build a model that included the temperature as an explanatory variable. As no geographical information was provided, we had to determine how to associate one or more weather stations with the load. We first used the GCV score as a criterion for choosing the best temperature for our model, with lower values being better. We chose to fit hourly GAMs in order to model the strong variation in the effects of the load, depending on the current hour. The hourly model used for temperature selection is as follows: Yt ,h = fh (posant ) + θDTk,h .DayTypek,t + θh .t
+ gh (Tt ,h,s ) + εt ,h ,
t = 1, . . . , n, h = 0, . . . , 23,
(4)
where Y is the load, posan is the day of the year (from 0 for January 1st to 1 for December 31st), fh is a nonlinear function requiring estimation, DayType is the day type (k = 2 for week days, k = 1 for weekends, and hence θDTk takes only two values), t is the trend, and gh is a nonlinear function for estimating the temperature at station s. Using Eq. (4), we built 25 models, one for each of the 25 available temperatures. As our models are hourly models, we obtained 24 GCV scores for each temperature variable. Fig. 4 shows the lowest GCV scores. We observed that the models’ performances for the deterministic and probabilistic forecasts were quite different; thus, the GCV criterion could not be used, as it consists of choosing the model with the lowest GCV score. We then tried to calculate a weighted average. We used an EWA algorithm (exponentially weighted average, see
4
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Tt ,i denotes the daily mean temperature for day t at station s, MAT t ,s is the corresponding moving-average (the window is set to 31 days), and Td∗,s is the historical mean temperature for day-of-the-year d at station s. Observations on February 29th are removed from the dataset. We have MAT t ,s =
15 1 T t − j ,s , 31 j=−15
t = 1, . . . , n, s = 1, . . . , 25 Td∗,s =
Fig. 4. The lowest GCV scores. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(Cesa-Bianchi & Lugosi, 2006)), which makes predictions by combining a finite set of provided forecasts. We used this algorithm on the 25 deterministic GAM forecasts (one for each weather station). We denote the deterministic forecast for the temperature Ts by F(Ts ). The EWA algorithm converges to the selection of three models with equal weights: 1 3
F (T6 ) +
1 3
F (T10 ) +
1 3
N t =1
MAT t ,s δd(t )=d ,
(6)
d = 1, . . . , 365,
(7)
where d(·) gives the day of the year, N is the number of years in our sample, and n is the total number of observations. Note that we first consider all years to have ∗ length 365, and take T60 , which is February 29th , to be a simple mean of the preceding and following days, in order to obtain normal temperatures on 366 days. After estimating these reference temperature variables, we model the deviation from them using autoregressive (AR) models. Let 1Tt ,s be the deviation from the normal temperature on day t at station s:
1Tt ,s = Tt ,s − Td∗(t ),s 1Tt ,s =
F (T13 ).
Ps
ρp,s 1Tt −p,s + εt ,s ,
(8) s = 1, . . . , 25.
(9)
p=1
We note that the three temperatures selected by EWA had already been among the six best temperatures selected by the GCV score. We therefore define a new temperature variable to be the mean temperature of stations 6, 10 and 13: Tt ,h,26 = 31 Tt ,h,6 + 31 Tt ,h,10 + 31 Tt ,h,13 . This new temperature variable improved both the deterministic and probabilistic forecasts. Load model. Using the new temperature variable T26 , we derived other potential explanatory variables. Of all of the models we have tested, the model that presents the best average score over the year 2011 is Yt ,h = fh (posant ) + θDTk,h · DayTypek,t + θh .t
+ gh (Tt ,h,26 ) + kh (STt ,h,26 ) + εt ,h , t = 1, . . . , n, h = 0, . . . , 23,
n 1
(5)
where STt ,h,26 is an exponentially smoothed version of Tt ,h,26 (with the smoothing coefficient α being estimated by cross-validation), and fh , gh , kh are non-linear functions that require estimation. Other models have also been tested during the competition, for example using different exponential smoothing parameters, the daily maximum temperature, and updating forecasts based on the consumption of the previous day (referred to as update1day in Table 2). The results are shown in Section 4. 3.2. Temperature simulation model A multivariate simulation model for the 25 weather stations is now set. First, the daily normal temperature for all weather stations is estimated using a moving average.
The order of the AR Ps can differ between variables, and the parameters are estimated via maximum likelihood. The residuals from the AR models are white noise, but do not need to be Gaussian. Using this method, we were able to simulate N ∗ paths for each weather variable:
• For each day to be forecast, T + 1, . . . , T + K (K depends on the month to be forecast), draw N ∗ vectors from εˆ t = (ˆεt ,1 , . . . , εˆ t ,25 ), denoted (εT(k+) 1 , . . . , εT(k+) k )k=1,...,N ∗ .
For a specific month to be forecast, we use only historical residuals from the same month and the preceding and following ones for sampling. • Simulate deviations from the normal temperature using Eq. (9). • Add corresponding normal temperatures (Td∗(T +k),s )s=1,...,25; k=1,...,K so as to obtain simulations of the daily weather variables. • To get to the hourly level, apply a daily profile. This profile can be either constant or month-dependent, as is shown in the center panel of Fig. 3. Tt ,h,26 can be derived easily from this simulation step. 3.3. Probabilistic forecasts We used the simulation model presented in Section 3.2 to produce 1000 paths for the input temperature of the load model in Eq. (5). Thus, we obtained 1000 paths for forecasting the load of the month. The number of paths is set arbitrarily. Once we had taken into account the temperature uncertainty, we wanted to model the uncertainty for our load forecasting. One way to do this was first to take a large
V. Dordonnat et al. / International Journal of Forecasting (
)
–
5
Fig. 5. Raw temperature (left) and smoothed temperature (right) effects for model (5).
Fig. 6. For T26 residuals of the AR model, the left panel of the figure shows the empirical ACF for the N observations based on kernel-based density estimation (black) and Gaussian distribution fit (red), while the center and right panels show the specified bandwidth and the time series, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
sample (e.g., 1 million) for each time-step to be forecast, by sampling with replacement from the fitting residuals εˆ t ,h of our load model in Eq. (5). From this sample, we then derived sequences of quantiles 1%–99%. The temperature and load uncertainties are combined by adding each sequence of quantiles of the load uncertainty to each simulated load path, based on a temperature path. We therefore obtained 99 000 final load paths from which to compute the final 99 forecasting quantiles. We tested several other ways of modeling the load uncertainty. Instead of taking a sample of the absolute fitting residuals (referred to as loaduncert1 in Table 2), we took a sample of the relative fitting residuals (referred to as loaduncert2 in Table 2). We also tried unsuccessfully to take fitting residuals by instant or by month, one-month forecasting residuals and residuals from 1- to 31-dayahead forecasting. 4. Results 4.1. Estimation results Load model. Fig. 5 presents the estimated temperaturebased effects of the load model in Eq. (5): the effect of
the raw temperature on the left is larger than that of the smoothed temperature on the right. The shape of the raw temperature effect is similar for all hours of the day, while the smoothed temperature effect is smaller during business hours (possibly due to residential heating). The smoothed temperature for cooling remains large, possibly due to cooling in offices. Temperature model. We present the simulation results for the temperature variable T26 that was defined earlier. According to Fig. 6 (left) and to the Ljung–Box test, it seems that there is no autocorrelation left in the AR model residuals. However, the distribution of the residuals is too spiked to be Gaussian. This motivated us to draw from the empirical residuals in order to simulate the temperature. As we can see from the right panel of Fig. 6, there is a clear yearly seasonality in the residuals. This is dealt with in the simulation step when we draw from the residuals. We only use residuals from the current month and the preceding and following months for forecasting. An alternative would have been to normalize the data using a moving average on squared data. This has been tried, but was no more efficient in our tests. The daily profile used for the simulation is constant over the months.
6
V. Dordonnat et al. / International Journal of Forecasting (
)
–
Table 1 Extract of forecast MAPE values. Model
Jan. 11
Feb. 11
March 11
Apr. 11
May 11
June 11
July 11
Aug. 11
Sept. 11
Oct. 11
Nov. 11
Mean
T6 T6_Tmax6 T6_ST6_0.992 T6_ST6_0.992_Tmax6 T6_ST6_0.996 T6_ST6_0.996_Tmax6 T6_ST6_ST10_0.992 T6_ST6_ST10_0.996 T26 T26_Tmax26 T26_ST26_0.992 T26_ST26_0.992_Tmax26 T26_ST26_0.996 T26_ST26_0.996_Tmax26
5.79 5.22 6.52 6.10 6.58 6.15 6.02 7.32 4.56 4.23 4.94 4.74 5.51 5.34
8.86 8.56 6.85 6.89 7.62 7.49 6.70 8.11 6.94 6.77 5.23 5.26 5.65 5.66
9.41 8.91 7.60 7.36 7.85 7.50 7.66 8.16 7.78 7.35 5.80 5.54 6.00 5.68
8.78 8.27 7.43 6.97 7.42 7.06 7.13 8.07 8.06 7.35 6.68 6.12 7.03 6.51
10.14 10.18 10.17 10.14 9.59 9.65 7.39 9.97 6.78 6.36 5.73 5.43 5.63 5.35
10.65 11.71 14.09 14.69 14.15 14.96 10.02 14.53 6.67 6.63 8.63 8.38 8.82 8.55
5.93 5.78 6.15 5.91 5.98 5.72 5.28 6.06 5.38 4.97 5.42 5.01 5.46 4.99
35.67 35.33 35.64 35.07 35.13 34.48 33.43 34.18 32.84 32.07 32.56 31.76 32.58 31.68
14.63 14.43 13.86 14.03 14.26 14.28 10.74 13.40 12.86 12.30 11.23 11.06 11.52 11.23
10.41 9.03 8.06 7.26 8.64 7.70 7.42 8.03 10.06 9.50 6.72 6.61 6.98 6.84
9.82 8.59 7.95 7.34 8.49 7.70 8.24 8.23 8.03 7.41 6.50 6.27 6.68 6.38
11.83 11.46 11.30 11.07 11.43 11.15 10.00 11.46 10.00 9.54 9.04 8.74 9.26 8.93
Note: The best forecast is shown in italics.
Table 2 Extract of forecasting results. Model
Task 4
Task 5
Task 6
Task 7
Task 8
Task 9
T6_T25_Tmax25_loaduncert1 T6_Tmax6_loaduncert1 T6_Tmax6_loaduncert2 T6_ST6_0.996_loaduncert1 T6_ST6_0.996_Tmax6_loaduncert1 T6_ST6_0.996_Tmax6_loaduncert2 T6_ST6_ST10_0.996_loaduncert2 T6_ST6_STmin6_0.996_loaduncert1 T6_ST6_STmin6_0.996_loaduncert2 T_PCA_loaduncert1 T26_ST26_0.992_Tmax26_loaduncert2 T26_ST26_0.996_loaduncert2 T26_ST26_0.996_Tmax26_loaduncert1 T26_ST26_0.996_Tmax26_loaduncert2 T26_ST26_0.996_Tmax26_update1day_loaduncert1
10.51 10.00 10.06 10.17 9.82 9.86 10.70 10.09 10.12
10.08 10.04 9.31 9.65 9.61 9.03 10.16 10.11
7.88 7.88 7.62 7.69 7.67 7.56 7.84 7.83
5.51 5.50 5.60 5.56 5.54 5.42 5.96 5.92
6.47 7.86 6.52 7.82 6.63 8.32 6.62 8.01 6.69 7.97 7.00 8.93 10.25 12.44 10.44 12.42
10.49 10.61 10.28 10.32
9.39 9.25 9.44 9.42
7.70 7.65 7.66 7.66
4.83 4.77 4.68 4.65 4.58 4.56 4.83 4.75 4.73 4.43 4.41 4.46 4.41
5.36 5.40 5.39 5.37
6.92 6.89 6.86 6.90
Task 10
8.55 8.71 8.67 8.63
Task 11
Task 12
Task 13
Task 14
Task 15
Mean
11.14 11.03 11.01 11.01 10.92 11.21 12.70 12.72
6.62 6.57 6.03 6.30 6.26 5.40 7.51 7.47
4.52 4.46 4.05 3.99 3.91 3.74 4.47 4.44
7.18 7.16 6.81 7.03 7.01 6.73 6.16 6.16
9.45 9.38 9.28 9.67 9.59 9.19 8.05 8.07
7.63 7.60 7.46 7.50 7.47 7.46 8.37 8.37
10.90 10.93 10.90 10.86
5.77 5.50 5.68 5.67 5.78
3.97 3.62 3.71 3.64
6.64 6.60 6.76 6.75
8.89 8.91 9.15 9.11
7.42 7.37 7.41 7.39
Notes: Effective forecasts are in bold, the best forecasts are in italics.
4.2. Forecast evaluation When trying to build an accurate point load forecasting model, the mean absolute percentage error (MAPE) has been used to select the best deterministic load model, by excluding years 2009, 2010 and 2011 from the estimation dataset in turn, and forecasting these years using historical temperatures. Table 1 gives the 2011 results (December results are not provided, due to the lack of temperature data). Thus, the estimation dataset was four years long. The pinball loss based quantile score was also used during the competition in order to test different probabilistic models for backcasting on the years 2009, 2010 and 2011. The submissions correspond to different models, with each submission being selected by looking at the best past scores for the same month. A summary of the scores is given in Table 2. Note that the best model when using a MAPE criterion was not necessarily the best one when using simulations and quantile forecasts. Tables 1 and 2 present the results obtained using the different models tested, where the name of the model indicates the variables used to model the weather dependent part. Model T26_ST26_0.992_Tmax26 is estimated with temperature T26 , the smoothed temperature ST26 with the
smoothing coefficient α = 0.992, and the daily maximum T26 temperature T max26 . loaduncert1 and loaduncert2 refer to the different approaches tested for the load model uncertainty. 4.3. Discussion Section 3.1 has presented what was our best deterministic load model on average. Other input variables have been tested as well:
• calendar-based variables: we also considered variables based on the year or month, bank holidays, and different day types. • temperature-based variables: we also considered lagged temperatures (one day, one hour, etc.), the daily maximum, minimum, mean and median, and a PCA (referred to as T_PCA in Table 2) of the 25 temperature variables. We selected variables using a MAPE criterion on validation datasets, by looking at the values of the estimated parameters and the graphs of the estimated splines; we therefore chose the day of the year instead of the month, a simple day type for week/weekend, and a linear trend. Our holiday variable did not provide any
V. Dordonnat et al. / International Journal of Forecasting (
improvement. As the MAPEs on October 2nd and 3rd, 2006, were very high, more than 30%, we removed these days from our fitting dataset. During the competition, we used past monthly scores (deterministic and probabilistic) to select one model for forecast submission. We recalculated the scores of the candidate models when new observations were released, and found that our best forecast was not always the one we had submitted. Since we found no rule for determining the best model to submit, a good way to improve would be by using model averaging. 5. Conclusion The proposed approach allowed the ADADA team to reach the top five of the load forecasting track. A deterministic load model has been developed, accounting for the day of the week as well as for the influence of weather. For providing probabilistic forecasts, we have used a resample of historical point forecast errors, together with a simulation model for temperature. Building a good probabilistic forecasting model for the electric load is a difficult task, especially when including the uncertainty on weather variables. Indeed, the best point forecasting model will not necessarily be the best when it is used with temperature simulation. It is also difficult to compare different modelling strategies, since only one observation for each hour of the day can be used for evaluating quantile forecasts.
)
–
7
Further ways to improve the modelling approach could be the inclusion of forecast averaging and the use of a better short-term forecasting model for temperature. References Cesa-Bianchi, N., & Lugosi, G. (2006). Prediction, learning, and games. Cambridge University Press. Hong, T., Pinson, P., Fan, S., Zareipour, H., Troccoli, A., & Hyndman, R. J. (2016). Probabilistic energy forecasting: state-of-the-art 2015. International Journal of Forecasting, this issue. Pierrot, A., & Goude, Y. (2011). Short-term electricity load forecasting with generalized additive models. In Proceedings of the 16th Intelligent System Applications to Power Systems Conference, ISAP 2011. Wood, S. N. (2006). Generalized additive models, an introduction with R. Chapman and Hall. Wood, S.N. (2013). Package mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness estimation.
V. Dordonnat holds an M.Sc. degree in statistics from the University of Versailles-Saint Quentin, France, and a Ph.D. degree from the VU University, Amsterdam, The Netherlands. From 2009 to 2015, she worked as a research engineer at EDF R&D, France. Since April 2015, she has been working at RTE, the French TSO, as a data scientist.
A. Pichavant holds a Master’s-level engineering degree from ENSAI, Rennes, France. She has been working at EDF R&D as a data scientist since 2012, in the load forecasting team.
A. Pierrot holds an M.Sc. degree in statistics from the University of Paris Dauphine, France. She has been working at EDF R&D as a data scientist since 2008, in the load forecasting team.