Novel approaches to the energy load unbalance forecasting in the Italian electricity market

Abstract

In the present paper we study the statistical properties of the Italian daily electricity load market, by mean of different statistical methods, such, e.g., the exponential smoothing model, the ARMA-ARIMA model and the ARIMA-GARCH model, also providing results about the goodness of each of the proposed approaches. Moreover, we show how the aforementioned models behave if exogenous regressors, as the day of the week or the temperature, are additionally taken into account. Analysed methods are then exploited to perform the one-day ahead energy load prediction, where the main focus is on guessing the right sign of the energy load unbalance.

1 Introduction

Energy Markets are financial markets that deal with the demand and supply of energy as well as to what concern financial derivatives structured on them. One of the most relevant characteristics of energy market frameworks relies on the fact that energy in general cannot be stored efficiently. Latter fact is just one of the main reasons why such particular type markets turn to be highly complicated both from theoretical, empirical and numerical point of view, see, e.g., [1] and references therein. During last years an increasing attention has been devoted particularly with the aim to obtain robust forecast methods capable to provide effective estimates about the production and consumption of, e.g., oil, gas, electricity and performances of the methods related to their productions, e.g., solar cells, oil plants, wind turbines, etc. Latter goals have to deal with the need to take into account both natural and social variables, such, e.g., weather conditions, firm needs, urban energy consumption, energy transportation, possible ways to storage electricity, as in the case of batteries, etc. Last but not least, almost every country has its particular type of energy market, with its own regulations, controlling bodies, etc. The sum of the aforementioned factors results in a financial framework which is also strongly characterized by high volatility levels, whose behaviour can be effectively studied exploiting the approach used in, e.g., [2], and references therein, or using more sophisticated approximation techniques, as those highlighted in [3], and references therein.

The latter implies that a general theory to treat energy markets, even if we restrict ourselves just to particular forecasting settings, is not a realistic goal and ad hoc studies have to be done in order to obtain effective results.

The present work aims at addressing a first fundamental task in energy trading in the Italian energy market, that is the forecast of energy load. Latter challenge is of crucial relevance particularly from the financial point of view by the side of agencies that produce and sell energy on the market under the supervision of the Italian Power Exchange (IPEX), managed by the Gestore del Mercato Elettrico (GME in Italian), which is the exchange for electricity and natural gas spot trading in Italy. In particular IPEX comprises two spot markets, namely the Day Ahead (in Italian: Mercato del Giorno Prima, or MGP) and the Intra Day (in Italian: Mercato Infragiornaliero, or MI). The latter, which since the 10th of February, 2015, is divided into five components, or sessions, which are called MI1, MI2, MI3, MI4, MI5.

Because there are not any efficient and economical sustainable ways to store electricity, power systems need to be constantly balanced between production and consumption, see, e.g., [4]. The latter implies that an accurate next day imbalance forecast has to be derived if one wants to obtain profit as well as to avoid losses caused by wrong imbalance in sign, whence the need to have concrete methods able to predict next days such a sign.

To the best of our knowledge, just few results have been already developed concerning last issue, moreover most of them, if not all, do not address the Italian market, but rather USA markets, as in the case of activities related to the California Energy Commission, or the Germany market, particularly with respect to the recently launched plan called Energiewende or energy turnaround.

We therefore intend to address this crucial topic in a series of paper, where we aim at giving a characterization as extensive as possible of the problem of energy load forecasting within the Italian energy market framework.

The present paper is so structured as follows: in Section 2 we give an overview of the Italian energy markets functioning and of its main peculiarities, addressing also the main problem that motivates our study; in Section 3 we provide a quick overview of the theoretical foundation of the method we later exploit in our analysis; eventually, in Section 4 we study the aforementioned problem with respect to a concrete real case.

2 The Italian energy markets functioning

The Italian power exchange market (IPEX) (Mercato elettrico italiano) is a free system that allows producers and consumers to enters into hourly contracts for buying or selling electricity. Such a market is divided into two main markets: future market and spot market, which are themselves divided into different sessions. We would like to underline that such a market implies a particular treatment of its financial basis since its nature is rather different from, e.g., the one characterising the usual derivatives/options/assets scenarios, see, e.g., [5, 6], and references therein.

In the future markets (FM) participants buy and sell bilateral contracts for delivery of energy on a specified future date. In such a contracts the parties are obliged to sell and buy the agreed amount of energy.

In the spot market (SM) any market operator has as a counterpart the transport system operator (TSO) (Gestore del Mercato Elettrico). Spot market is divided into the day ahead market (mercato del giorno prima) (MGP), the intra day market (mercato infragiornaliero) (MI), which since February 2015 is composed of five sessions (MI1, MI2, MI3, MI4 and MI5), and the Mercato per il Servizio di Dispacciamento (MSD).

The day ahead market

In the MGP negotiations for energy trading for any hour of the next day take place. This market is based upon an implicit auction mechanism, each player submits his bids composed by a quantity and a price representing the maximum price at which he is willing to buy energy or the minimum price at which he is willing to sell energy.

At closure of each session, one for every hour of the day, all the offers are processed and either accepted or refused according to the System Marginal Price, that is, bids, resp. asks, are ordered from the lowest to the highest price, resp. from the highest to the lowest price, the equilibrium price and total exchanged energy are determined by the intersection of the two curves.

The intra day market

The MI is composed of five sessions, namely MI1, MI2, MI3, MI4 and MI5. MI1 and MI2 take place the day before the actual delivery, whereas MI3, MI4 and MI5 occur the same day of the actual delivery. In each session every operator can modify his program of injection or withdrawal of energy. Also MI sessions follows the same exact rules of price formation of MGP.

The dispatching services market

Dispatching guarantees the overall equilibrium between production and loading and thus ensures the correct functioning of the national electric grid. In Italy the dispatching system is managed by Terna S.p.A., which is the owner of the high voltage national transmission network. Terna, in order to guarantee the proper functioning of the electrical network, has to deal with the congestion resolution activities between the different market areas, the creation of the reserve of energy and the real-time balance between production and consumption.

In MSD Terna obtains the necessary reserves to the dispatching service by acting as central counterparty in negotiations with operators enabled to the dispatching service. In this market all accepted bids are remunerated at the price presented, according to the paid as bid method. The offers of purchase in MSD are also called downward, meaning that such offers will be accepted if it is necessary to reduce the amount of energy generated, while offers to sell are called upward.

2.1 The problem of the unbalance forecasting

As already mentioned, the continuous balance between production and consumption of energy is a fundamental task in order to guarantee the correct functioning of the whole electrical network. This balance is guaranteed by Terna through MSD. In order to create the necessary energy reserves to balance the grid, Terna needs to know as precisely as possible the production of different plants.

While for traditional sources plants, such as coal, gas and other fossil fuels, is a relatively simple task to predict the next day production, for non-programmable renewable sources plants, such as wind and solar, this forecast is very difficult task, with many factors affecting the final outcome.

In order to ease the work of Terna to balance the network, all the producers from traditional sources and programmable renewable sources, such as for instance some types of hydroelectric energy, are required to provide to Terna the exact production plan for the next day; in the event that these programs are not observed, then the actor has to pay a penalty. Producers of non-programmable renewable sources that do not meet the scheduled production, however, incur in penalties or rewards depending on the relative sign between their unbalance and the unbalance of the macro-zone in which the plant is located.

In Italy there are two different macro-zones for balancing purposes: the northern macro-zone consists of all the regions of Northern Italy, including Emilia Romagna, whereas the southern macro-one consists of all other regions. The aggregate zonal unbalance is the algebraic sum, changed in sign, of the amount of energy procured by Terna in MSD at a given time in a given macro-zone. When the aggregate unbalance is positive means that the energy produced is greater than the energy scheduled and then most of the offers accepted in MSD were downward; when the aggregate unbalance is negative the opposite happens.

A producer of non-programmable renewable sources which has produced more than declared in an area with positive aggregate unbalance will be required to pay a penalty, as it has helped to increase the unbalance in the area. Similar thing happens in the case in which the actor has produced less than declared and the macro-zone has negative aggregate imbalance. Conversely, if a producer is unbalanced in the opposite direction to the macrozonal aggregate imbalance he is rewarded as he is helping to balance the market. Table 1 gives reward and penalties for non-programmable renewable plants.

Table 1 Imbalance price settlement mechanism

Energy that was scheduled the day before is paid at the price determined in MGP $P_{\mathrm{MGP}}$, whereas unbalanced energy is paid according to Table 1. Above we have defined by $P_{\mathrm{MB}}^{\downarrow}$, resp. $P_{\mathrm{MB}}^{\uparrow}$, the average price used by Terna to decrease, resp. increase, the generation of energy. Moreover the following relation holds:

$$P_{\mathrm{MB}}^{\downarrow} \leq P_{\mathrm{MGP}} \leq P_{\mathrm{MB}}^{\uparrow} . $$

3 Statistical methods

In the present section we analyse the theoretical mathematical foundation of the method that will be later used to concrete address our energy balance task. In particular, in order to make the present work as much as self-contained as possible, we provide an overview of standard approaches that can be used in similar frameworks, referring the interested reader to, e.g., to [7, 8] for a deeper introduction to them.

The exponential smoothing model :

Let us first consider the exponential smoothing (ES), see, e.g., [7–9], which is mainly based on predictive procedures built starting from an exponentially weighted average of past observations. The general model involves a state vector $y_{t} = (l_{t},b_{t},s_{t},\dots,s_{t-m+1})$, where $l_{t}$ represents the level of the series, $b_{t}$ represents the growth and $s_{t}$ is the seasonal component, coupled with a state space equation of the form

$$\textstyle\begin{cases} X_{t} = w(y_{t-1}) + r(X_{t-1}) \epsilon_{t}, \\ y_{t} = f(y_{t-1}) + g(X_{t-1}) \epsilon_{t}; \end{cases} $$

$\epsilon_{t}$ being a centred Gaussian noise with finite variation $\sigma^{2}_{t}$, and for some suitable functions, to be choose according to the model one wishes to fit, w, r, f and g, see, e.g., [10] or also [8]. We remark that in what follows we will not choose a particular ES model to be fitted to data, but rather we run a routine in order to choose the best performing ES model.

The ARMA-ARIMA model :

When one concerns the study of time series, the AutoRegressive Moving Average (ARMA) models play a central role because they are capable of describe weakly stationary stochastic processes with a rather restricted set of assumptions, being mainly based on the use of two polynomials: the first one takes into account the autoregressive character of the data set, while the second takes into account the moving average. It is worth to mention that such a method results as a combination of the Moving Average method (MA) together with an AutoRegressive (AR) one. In particular, denoting by $X_{t}$ the unknown value of the series of interest at time t, which is in fact treated as random variable, a pth order Auto Regressive method ($\operatorname{AR}(p)$) is defined as follows

$$X_{t} - \sum_{k=1}^{p} \phi_{k} X_{t-k} = \epsilon_{t} , $$

where $\epsilon_{t}$ is a general random noise, while the coefficients $\phi_{1},\dots,\phi_{p}$ are the AR (or regression) coefficients. In the most simple case the noise $\epsilon_{t}$ is assumed to be Gaussian, however more general type of random disturbance can be considered.

Concerning the moving average component of the ARMA model, it is defined by the Moving Average method of order q, which will be indicated by $\operatorname{MA}(q)$, and it is defined as follows

$$X_{t} = \epsilon_{t} + \sum_{k=1}^{q} \theta_{k} \epsilon_{t-k} , $$

$\epsilon_{t}$ still being the noise, or uncertainty, affecting our observations, or elements of the time series we are studying, which is not necessarily of Gaussian type. As before X is the process that we would like to forecast, on the basis of previous observations, while $\theta_{1},\dots,\theta_{q}$ are the moving average parameters.

Merging the methods already introduced, it is possible to define the AutoRegressive Moving Average method of order $(p,q)$, indicated by $\operatorname{ARMA}(p,q)$, and defined as follows

$$ X_{t} - \sum_{k=1}^{p} \phi_{k} X_{t-k} = \epsilon_{t} + \sum _{k=1}^{q} \theta_{k} \epsilon_{t-k} . $$

(1)

A further step is represented by the AutoRegressive Integrated Moving Average (ARIMA) method, which allows to take into account time series which are not of stationary type. In particular an ARIMA model of order $(p,d,q)$, indicated by $(\operatorname{ARIMA}(p,d,q))$ where d is the degree of differencing, namely it represents how much the time series we are dealing with is far from being stationary, is defined as follows

$$\phi(B) \nabla^{d}X_{t} = \theta(B) \epsilon_{t} , $$

where $\nabla X_{t} = (1-B)X_{t}$ is the lag 1 differencing operator, and B is the backward shift operator defined by $B^{n} X_{t} := X_{t-n}$, considering the use of the following short notations

$$\begin{aligned} &\phi(B) = 1- \phi_{1} B^{1} - \cdots- \phi_{p} B^{p} , \\ &\theta(B) = 1 + \theta_{1} B^{1} + \cdots+ \theta_{q} B^{q} . \end{aligned} $$

We refer the interested reader to, e.g., [8], Section 3.4, for further details.

The ARMA-GARCH model :

The ARMA-GARCH model models the mean equation via an $\operatorname{ARMA}(p,q)$ model, see Eq. (1), whereas the random noise components, represented by $\epsilon_{t}$, are modelled with a Generalized Autoregressive Conditional Heteroskedastic (GARCH), model. In particular the GARCH model of order $(p,q)$, indicated by $(\operatorname{GARCH}(p,q))$, defines the residual $\epsilon_{t}$ appearing in equation (1) as follows

$$\begin{aligned} &h_{t} = \epsilon_{t} \sigma_{t} , \\ &\sigma^{2}_{t} = \alpha_{0} + \sum _{k=1}^{q} \alpha_{k} h_{t-k}^{2} + \sum_{k=1}^{p} \beta_{k}^{2} \sigma^{2}_{t-k} , \end{aligned} $$

for some positive coefficients $\alpha_{k}, \beta_{k} \geq0$. The latter approach allows to overcome one of the main issue affecting the ARMA process, namely the fact that the mean equation cannot take into account for heteroskedastic effects of the time series process, as, e.g., happens for the so called fat tails.

ARMA-GARCH method with exogenous variables :

The ARMA-GARCH models can be enriched by considering also the role played by the so called exogenous variables. In particular exogenous variables can be added to the ARMA component, as well as to the GARCH model one. However, for the sake of brevity, here we only consider the case of exogenous variables added to the ARMA process, the case related to the GARCH one being analogous.

We thus define the AutoRegressive Moving Average method with exogenous variables ($\operatorname{ARMAX}(p,q,r_{1},\dots,r_{n})$), as follows

$$\phi(B)X_{t} = \theta(B)\epsilon_{t} + \sum _{k=1}^{n} \psi^{k}(B) Y^{k}_{t} , $$

where $r_{i}$’s are the orders of the exogenous variables $Y^{1},\dots ,Y^{n}$ and

$$\psi^{k}(B) = \psi^{k}_{0} + \psi^{k}_{1} B + \cdots+ \psi^{k}_{r_{i}} B^{r_{i}} , $$

and we refer to [8], Section 3.4.9, for a detailed introduction to modelling with exogenous variables.

4 Estimating the model

In what follows we will apply the methods introduced in Section 3 to analyse the net hourly energy load imbalance in the Italian northern macro-zone, where the positive, resp. negative, sign is to be intended as explained in Section 1. We recall that the Italian northern macro-zone is composed by all regions from northern Italy, including Emilia Romagna.

First we study the hourly time series of energy load, applying time series method recalled in Section 3 to outline which of them better perform according to a concrete criterion that will be specified later on. Then, we will consider the task of forecasting the next day imbalance sign. We would like to underline that we are mainly interested not in the prediction of the exact amount of the zonal imbalance, but rather in predicting the right sign of imbalance, since the latter is the main factor affecting the energy trading in Italy, nowadays. The latter characteristic is due to the imbalance mechanics on which the MSD is based, see Table 1 and Section 2.1, for details. Last but not least, the forecast of the next day imbalance sign for the Italian energy market constitutes the main novelty of the present work.

As often have been pointed out in literature, see, e.g., [8], instead of considering the hourly times, it is more effective to deal with 24 different daily time series, one for every single hour of the day. In fact, the hourly time series results to be highly intractable from a statistical point of view, with several trend and seasonal components, as it is clearly shown in Figure 1, where the auto correlation function, Figure 1(b), and the partial auto correlation function, Figure 1(c), exhibit a clear daily component with a significance correlation at lag 24.

Figure 1

Hourly energy load time series.

Therefore, motivated by previous facts, we perform our study on the time series of the effective unbalance at a given hour. In particular, in what follows we focus our attention on one of the most challenging hour of the day, namely we consider the time series at 2 p.m. Nevertheless, it is worth to mention that different hours, albeit being characterized by their specific peculiarities, can be analogously treated exploiting our approach. It is worth to mention that our analysis will be performed with respect to different time windows, showing how changing the number of days taken into consideration may slightly affect the overall fit of the model. We stress that the problem of choosing the right time window is really a hard task in energy markets, mainly because it is affected by a large number of seasonal components. The latter implies that, going too far in the past, may only lead to an increase of the overall instability of the model. Furthermore, since the national regulation of the Italian energy market has been changed several times during last years, one has to take past values with particular care on the chosen period in order to avoid to treat non homogeneous numbers resulting as the output of different regulatory settings.

On the basis of previous considerations, our study is based on the time series of the net energy imbalance in the northern macro-zone, that starts from the 1st of January 2015, see Figure 2, resulting in 509 daily observations.

Figure 2

Time series $\pmb{h:14:00}$ .

A first choice on the best method to use can be done looking at both the sample autocorrelation function (ACF) and the partial autocorrelation function (PACF), see Figure 2 bottom panel. In fact, also taking into consideration the graphs in Figure 2, the considered time series seems to be non-stationary, with clear evidences of seasonal behaviours. Nevertheless, in what follows we will choose the best model looking at the bias-corrected Akaike’s Information Criterion (AIC), defined as follows

$$\mathrm{AIC} := - 2 \log\mathcal{L} + \frac{2 d n}{n-d-1} , $$

being d the size of the mode, while $\log\mathcal{L}$ represents the log-likelihood function defined as

$$\log\mathcal{L} =- \frac{n}{2} \log\bigl(2 \pi\sigma^{2}\bigr) - \frac{1}{2} \log\Gamma- \frac{1}{2 \sigma^{2}}L^{T} \Gamma^{-1}L , $$

being $L = (L_{1},\dots,L_{n})^{T}$ the observations of a stationary time series with $L^{T}$ its transpose, $\sigma^{2}$ the noisy variance and Γ is the auto-covariance matrix of L, see, e.g., [8], p.84. We stress that a lower AIC means a better fit. Other widely used criterion are the Akaike’s Final Prediction Error (FPE) and the Bayesian Information Criterion (BIC), sometimes also referred to as Schwarz Information Criterion (SIC), we refer the interested reader to such criterion to [8], Section 3.4.6. All of aforementioned criteria provide insights about the goodness of statistical fit, hence giving precious informations on how to select the most performing model. In particular they are characterized by the peculiarity of penalizing artificial overfit which is caused, in most of the cases, by an excessive complexity, e.g., too many parameters have been considered, with respect to the chosen model.

The model selection procedure consists of several steps. Concerning the ES model, after having chosen a time window, we estimate it considering the associated AIC value. Afterwards, $\operatorname{ARIMA}(p,d, q)$ processes are calibrated to the stochastic component for different values of parameters p and q, then the AIC criteria is used to choose the best $(p,q)$ combination. Analogously, we proceed for the simple ARMA-GARCH model, resp. for the ARMA-GARCH enriched by exogenous variables. As regard the exogenous variables, several choices can be made. As an example, Figure 2 enlightens that the considered time series exhibits a weekly seasonality, so that the first exogenous variable we take into account is the day of the week. The second exogenous parameter that has to be taken into account is temperature. In particular Figure 3 represents the ACF of the net load time series plotted against temperature, there is a clear evidence of a positive correlation between energy load and the temperature values, see also, e.g., [8], Section 3.2.3, for a deeper treatment of the topic.

Figure 3

ACF function of the load time series against temperature.

Once we have calibrated all the aforementioned models, we exploit the associated AIC values to chose the one that has better performed. Eventually, the same procedure is repeated over different time windows. It is worth to mention, as briefly said above, that choosing the right time windows to be used to perform the statistical analysis turns to be a rather difficult task. In the framework of energy market, latter problem is even more complicated, playing, at the same time, a more relevant role. In fact, energy related time series often exhibit different seasonal components, being the energy load correlated to the day of the week as much as the season on is considering. Besides, it can happen that some exogenous variables may affect the data only in a given season, being irrelevant during the others. Concerning the latter issues, we refer the interested reader to [8], Section 3.2.4.

We list the resulting AIC values for all the exploited models, with respect to particularly relevant time windows.

A first immediate consequence that can be drawn from Table 2 is that the ARMA(X)-GARCH performs considerably better than the other models, with a slightly increase in obtained accuracy when the exogenous variables are taken into account. Nevertheless, Table 2 also implies that the main improvement is given by considering the GARCH model, rather than the predictor variables. In fact, it can be seen that also the ARIMAX model performs better than the ARIMA one, but its fit is considerably worst than the ARMA-GARCH.

Table 2 Resulting AIC values for considered models over different time windows

A further comparison between the ARMA-GARCH type models with and without exogenous variables is shown in Figure 4. In fact it can be seen how both the empirical distribution of standardized residuals, Figure 4 top panel, and the 99% confidence interval, Figure 4 bottom panel, for the ARMA-GARCH and the ARMAX-GARCH model behave similarly.

Figure 4

Empirical distribution of standardized residuals (top panels) and 99% VaR (bottom pannels).

Finally, in Table 3 ,we have reported the AIC criterion for each model, whit respect to some relevant hours of the day, such as the net load at $h: 03:00$, $h: 10:00$, $h: 18:00$ and $h: 22:00$. An immediate glimpse at Table 3 shows how the same conclusions drawn above can be done also if one consider different hourly time series.

Table 3 AIC value for each statistical model at different time of the day

4.1 On the forecast for the next day unbalance

As mentioned at the beginning of the current Section, we are mainly interested in the forecast of the next day energy load, with particular attention to the overall sign of the imbalance in the macro-zone, rather than to what concerns the exact quantity of imbalanced energy. Moreover, our main goal is to obtain accurate short time previsions of the right sign of next days imbalance. The latter is due to the specific mechanism behind the trading strategies explained in Section 2.1, see in particular Table 1. For a general treatment of forecasting within different energy markets and/or with respect to different scopes, we refer to, e.g., [2, 11–13], while we refer to, e.g., [11, 13] for the study of the long-time horizon forecast and to [14], and references therein, for the treatment of related computational issues.

In particular, our forecasting procedure is structured as follows: we chose an appropriate positive threshold and, if the forecasted value is in absolute value higher than the threshold, then we enter the market according to the predicted unbalance sign, otherwise, namely, as it is most likely to happen, the outcome is too close to zero implying a high probability to unbalance in the wrong direction, we do not enter the market. We underline that the aforementioned threshold can be chosen according to different parameters, for instance we considered the estimated volatility of the time series.

In Figure 5 we show the next day predictions obtained exploiting different models, but with a fixed time window, for the time series at $h: 14:00$. The top panel in Figure 5 represents the worst performing methods according to AIC value, that is ES model, Figure 5(a), and ARIMA model, Figure 5(b); grey bounds represent the 99% and 95% prediction interval. The bottom panel instead is concerned with the two best performing model, that is the ARMA-GARCH model, Figure 5(c), and the ARMAX-GARCH model, Figure 5(d); the yellow bound in the bottom panel represents the 95% prediction interval. From Figure 5 it can be seen how the fact that the models have different goodness of fit implies that the forecast for next days energy load is different. In fact both ES and ARIMA models have almost constant forecasted values, whereas the ARMA-GARCH and ARMAX-GARCH models seems to be more sensitive to oscillations in the energy load imbalance.

Figure 5

Next days predicted values.

Eventually, in Table 4, we show the performances of each of the used methods. The evaluation has been obtained according to the following criterion: we fix a time window and we fit all the models, then we check the volatility of the times series to decide the best threshold, using it to predict the next day imbalance. Then, if the forecasted value is above the threshold in absolute value, we enter the market. Afterwards, we check with the actual datum if the predicted imbalance sign is correct, and we shift the time series to repeat the procedure. As a result, Table 4 shows in the left part the number of times we have entered the market, along with number of times we would have guessed the right, resp. wrong, imbalance sign. Latter results have been obtained using a fixed time window of 100 days and still considering the daily times series at $h: 14:00$. The right hand side of Table 4, reports the overall performance reached by each method over the whole 509 days composing our global data set, considering all the associated daily time series.

Table 4 Different methods: overall performance

Table 4 shows on one side the expected results, while on the other it highlights some peculiar features. In fact, as expected from results reported in Section 4, the ARIMA model is poorly performing, being indeed the worst when dealing with the daily time series as well as considering the overall performance. Viceversa, again as we expect, the ARMA-GARCH model performs rather well, with a slight improvement when one also considers the exogenous variables. What turns to be an unexpected result concerns the ES model performance, since, according to the AIC criterion, see Table 2, it is by far the worst performing one. Nevertheless, if one addresses the problem of forecasting it can be seen that the ES model outperforms the ARMA-GARCH model without exogenous variables, being also, even if by a few percentage, better than the ARMA-GARCH with exogenous variables. In particular, it appears that the ES model and the ARMAX-GARCH model perform similarly in predicting the next day sign, being the main difference among the two represented by the number of times one enters the market which is greater for the ES, most probably because it overestimates the next day outcome, leading to a higher number of plays.

5 Conclusion and further developments

The present work constitutes a first step towards the solution of the highly difficult task of next day energy imbalance forecast. Such an ambitious goal is affected by several issues due to many different reasons, such, e.g., the difficulty to find statistically good data, which means that time series are rather often dirtied by regulatory changes, exogenous noises, sensors faults, etc. Besides, a rather important issue, which is typically not considered, will play a fundamental role in the next future, namely the one concerning the problem of optimal allocation/transportation of energy resources/products, a particularly difficult task that is intrinsically linked to the solution of stochastic optimal problems stated on networks, and whose solution has been the subject of an increasing number of researches during recent years, not only within the energy market framework, see, e.g., [3, 15–17] and references therein.

Moreover, different factors can affect energy loads, and we have considered just the day of the week and temperature. Even if the latter appears, from standard literature, to be the most relevant, nevertheless recent studies have shown that also different factors could play a crucial role, as for the case of renewable energies. In particular, considering renewable energies lead to at least two non-trivial problems: exact values for such an exogenous parameter are not always available, moreover the exact value of its next day production is rather tricky to be predicted with enough accuracy. Previous reasons suggest a very careful and detailed study, that goes beyond the aims of the present work, but whose results will be of great relevance since the production of renewable energy plays a fundamental role in the zonal unbalance, that is why we will address this key task in a future work. The second being how to chose the most relevant renewable energy with respect to its impact on energy loads, since the effect of different renewable energies may vary a lot from region to region.

We outline that the latter point cannot be neglected in order to develop a solid method to predict future energy loads, indeed such a subject will be the main focus of our future works. Last but not least, we would like to underline that the most of the computational part which has been developed so far with respect to the type of problems we have analysed in the present paper, is mainly based on Monte Carlo type techniques. Such a type of numerical approach is particularly ineffective for our purposes, because of its slow rate of convergence and poor accuracy, at least compared to more sophisticated methods as the ones based on the Polynomial Chaos Expansion approach, see, e.g., [14] and references therein.