Ensembling methods for countrywide short term forecasting of gas demand

TL;DR

This study compares various forecasting models for short-term gas demand prediction in Italy, demonstrating that ensemble methods outperform individual models and existing TSO predictions, with potential for improved accuracy.

Contribution

It introduces and evaluates ensemble forecasting methods for gas demand, showing their superiority over individual models and current industry predictions.

Findings

Ensemble predictors outperform base models in accuracy.

Models outperform Transmission System Operator (TSO) predictions.

Ensemble methods show consistent improvements in out-of-sample validation.

Abstract

Gas demand is made of three components: Residential, Industrial, and Thermoelectric Gas Demand. Herein, the one-day-ahead prediction of each component is studied, using Italian data as a case study. Statistical properties and relationships with temperature are discussed, as a preliminary step for an effective feature selection. Nine "base forecasters" are implemented and compared: Ridge Regression, Gaussian Processes, Nearest Neighbours, Artificial Neural Networks, Torus Model, LASSO, Elastic Net, Random Forest, and Support Vector Regression (SVR). Based on them, four ensemble predictors are crafted: simple average, weighted average, subset average, and SVR aggregation. We found that ensemble predictors perform consistently better than base ones. Moreover, our models outperformed Transmission System Operator (TSO) predictions in a two-year out-of-sample validation. Such results suggest that combining predictors may lead to significant performance improvements in gas demand forecasting.

Figures

Tables

Table 1. Table 1: Features

Feature	Reference time	Type	Series
Gas demand series	t-1	continuous	RGD, IGD, TGD
Gas demand series	t-7	continuous	RGD, IGD, TGD
Gas demand series	$\mathrm{sim}\left(t\right)$	continuous	RGD, IGD, TGD
Gas demand series	$\mathrm{sim}\left(t-1\right)$	continuous	RGD, IGD, TGD
Forecasted temperature	t	continuous	RGD, IGD, TGD
Forecasted temperature	t-1	continuous	RGD, IGD, TGD
Forecasted temperature	t-7	continuous	RGD, IGD, TGD
Forecasted temperature	$\mathrm{sim}\left(t\right)$	continuous	RGD, IGD, TGD
Forecasted HDD	t	continuous	RGD, IGD
Forecasted HDD	t-1	continuous	RGD, IGD
Forecasted HDD	t-7	continuous	RGD, IGD
Forecasted HDD	$\mathrm{sim}\left(t\right)$	continuous	RGD, IGD
Forecasted HCDD	t	continuous	TGD
Forecasted HCDD	t-1	continuous	TGD
Forecasted HCDD	t-7	continuous	TGD
Forecasted HCDD	$\mathrm{sim}\left(t\right)$	continuous	TGD
Weekday	t	categorical	RGD, IGD, TGD
Holiday	t	dummy	RGD, IGD, TGD
Day after holiday	t	dummy	RGD, IGD, TGD
Bridge holiday	t	dummy	RGD, IGD, TGD

Equations

$${\mathrm{GD}}_t={\mathrm{IGD}}_t+{\mathrm{TGD}}_t+{\mathrm{RGD}}_t$$

$$\mathrm{HCDD}=\left|T_c-T\right|$$

$${\tau }^{\ast }=\arg \underset{\tau }{\min }|\mathrm{yearday}(\tau )-\mathrm{yearday}(t)|$$

$$f(\text{x})=\sum _{j=1}^p{x}_j{\beta }_j={\text{x}}^T\mathbf{\beta },\mathbf{\beta }\in {\mathbb{R}}^p$$

$$\mathbf{f}=\left[\begin{array}[]{cccc}f(\mathbf{x_{1}})&f(\mathbf{x_{2}})&\ldots&f(\mathbf{x_{n}})\end{array}\right]^{T}=\mathbf{X}\boldsymbol{\beta}$$

$${\mathbf{\beta }}^{\text{ridge}}$$

$${\mathbf{\beta }}^{\text{LASSO}}$$

$${\mathbf{\beta }}^{\text{elastic net}}$$

$$\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }^2=\sum _{i=1}^{n}L(y_i,f({\mathbf{x}}_{\mathbf{i}})),L(y,\hat{y}):=(y-\hat{y}{)}^2$$

$$L_{ϵ}(y,\hat{y}):=\{\begin{array}{ll}0,& |y-\hat{y}|<ϵ\\ |y-\hat{y}|-ϵ,& \text{otherwise}\end{array}$$

$$\Vert f{\Vert }_K=\sum _{j=1}^{\infty }\sum _{i=1}^{\infty }{\alpha }_j{\alpha }_{i}K(x_j,x_i)$$

$$f(\mathbf{x})=\sum _{j=1}^{\infty }{\alpha }_{j}K(\mathbf{x},{\mathbf{x}}_{\mathbf{j}})$$

$$f^{\text{SVR}}:=\arg \underset{f\in \mathcal{H}}{\min }\sum _{i=1}^n{L}_{ϵ}(y_i,f({\mathbf{x}}_{\mathbf{i}}))+\frac{\lambda }{2}\Vert f{\Vert }_K^2$$

$$f^{\text{SVR}}(\mathbf{x})=\sum _{i=1}^n{c}_{i}K({\mathbf{x}}_{\mathbf{i}},\mathbf{x})$$

$${\hat{f}}^{\mathrm{A}}({\mathbf{x}}_{\ast })=\frac{1}{M}\sum _{i=1}^M{\hat{f}}_i({\mathbf{x}}_{\ast })$$

$$\{\begin{array}{l}{\hat{f}}^{\mathrm{LS}}({\mathbf{x}}_{\ast })={\sum }_{i=1}^M{w}_i{\hat{f}}_i({\mathbf{x}}_{\ast })\\ w_i\ge 0\\ {\sum }_{i=1}^M{w}_i=1\end{array}$$

$$\sum _{k=2}^8\left(\genfrac{}{}{0ex}{}{9}{k}\right)=501$$

$$\mathrm{MAE}=\frac{1}{n}\sum _{t=1}^n|y_t-{\hat{y}}_t|$$

Full text

1 Introduction

Natural gas is one of the most important energy sources in Italy: it feeds domestic and industrial heating, production processes and thermoelectic power plants. Data from SNAM Rete Gas, the Italian Transmission System Operator (TSO), show that the total Gas Demand (GD) is made of three main components: Residential Gas Demand (RGD), Industrial Gas Demand (IGD), and Thermoelectric Gas Demand (TGD). In 2018, RGD accounted for 41.5% of the total consumption, IGD for 25.4% and TGD for the remaining 33.1% [1].

Accurate forecasts of the overall GD, as well as of its three main components, are of primary importance to energy providers, in order to improve pipe reservation and stock planning and also prevent financial penalties due to network unbalance. Moreover, GD is closely correlated with natural gas price, which is a key input for determining the optimal production plan of thermal power plants.

Several works addressed the forecasting of natural gas demand: comprehensive reviews are [2] and [3]. The latter proposes a classification along four dimensions: geographical area, time horizon, method, and inputs. Herein, we are interested in country-level, one-day-ahead predictions, based on statistical learning models that leverage past gas demand, temperature and calendar features as input variables.

With respect to prediction horizon, a distinction is made between long-term forecasting, featuring an horizon of months or years and short-term, with an horizon of one or few days. Focusing on country-wide predictions, a long-term model based on temperature was proposed in [4] to forecast Turkish demand. The importance of the relation between weather and gas demand was also highlighted in [5] and [6]. In [7] a statistical model was applied to forecast the long-term evolution of Slovenian demand, while different kinds of so-called "grey models" were applied in [8] and [9] to forecast Chinese demand.

In [10], short-term forecasting of UK natural gas demand was addressed using support vector regression with false neighbours filtered. According to the authors, the method performed better than Auto-Regressive Moving Average (ARMA) models and neural networks (ANN). Azadeh et al. [11] proposed an adaptive network-based fuzzy inference system (ANFIS) to predict Iranian gas demand, which improved on classical time series methods and ANN. A more advanced model, combining wavelet transform, genetic algorithm, ANFIS and ANN was applied in [12] to the Greek gas distribution network.

Long-term evolution of the Italian gas demand is investigated in [13] and [14]: macroeconomic indicators, such as gross domestic product and gas prices and climatic factors are used to build scenarios of RGD and overall GD evolution up to 2030. However, to the best of our knowledge, the daily series of Italian GD and its peculiar features have not been studied in the literature and no result about its short-term forecasting has been presented.

A previous work [15], focusing on Italian RGD, proposed and compared five prediction models: ridge regression, Gaussian Process (GP), nearest neighbours, Artificial Neural Networks (ANN), and torus model, concluding that ANN and GP provided the best results. Herein, we extend the analysis along three directions: first, also the prediction of IGD and TGD is addressed, thus enabling the prediction of the overall Italian GD; second, four additional base forecasters (LASSO, elastic net, random forest, and support vector regression) are considered; third and finally, the use of ensemble predictors, i.e. forecasters obtained by the suitable aggregation of base forecasts, is investigated. More precisely, based on the nine base models (the five discussed in [15] and the four additional ones), four ensemble predictors are considered: simple average, weighted average, subset average, and support vector regression aggregation.

Ensembling, also known as blending, is known to be an effective technique to improve overall accuracy and stability, see e.g. [16, 17]. Recently, ensemble predictors have been proven successful in forecasting electric load [18], whose series shows a periodic structure similar to the one of GD.

The contribution of this paper is thus threefold: first, we present and discuss the statistical properties of Italian IGD and TGD; on these data, we develop, apply, and compare nine machine-learning models; finally, we explore the use of ensemble predictors, assessing the consequent improvements.

The paper is organised as follows. In Section 2 we present the forecasting problem and the available data, while in Section 3 we describe the most relevant features of IGD and TGD time series. After discussing feature engineering (Section 4), in Section 5 the adopted models are concisely overviewed. Training and hyperparameter tuning are presented in Section 6. Results are reported in Section 7 and some concluding remarks (Section 8) end the paper.

2 Problem statement

In this paper, the prediction of the Italian daily GD is addressed both at aggregated and disaggregated level. Apart from minor components that can be neglected, at any day $t$, the overall GD is given by the sum of Industrial Gas Demand (IGD), Thermoelectric Gas Demand (TGD) and Residential Gas Demand (RGD):

[TABLE]

IGD includes demand by industrial plants, while TGD only accounts for the fuel required by thermoelectric power plants. The residential component and its forecasting was studied in a companion paper [15] and is therefore not further discussed herein.

Three one-day-ahead forecasting problems are considered: (i) IGD, (ii) TGD, and (iii) overall GD. The last forecast will be obtained by summing (i) and (ii) with the prediction of RGD.

The datasets for both TGD and IGD are 12 year long, ranging from 2007 to 2018, and consist of 3 fields: date ($t$), forecasted average temperature in Northern Italy ($T$)111Weather forecast were provided by an Italian specialised company, and gas demand (either IGD or TGD). Temperature in Northern Italy was considered as this region has the most rigid climate, and is thus more sensible to heating requirements. In Fig. 1 the complete series of RGD, IGD, TGD and overall GD are displayed.

3 Exploratory analysis

3.1 Industrial gas demand

Industrial Gas Demand (IGD) does not exhibit strong trends: a significant decrease is only recorded in 2009, following the financial crisis started in the previous year. The series presents weekly and yearly seasonal patterns. In particular, as most of the industrial facilities stop or slow down production during the weekend, IGD is lower on Saturdays and Sundays. In August and at the end of December, typical holiday periods, IGD drops to about half of its average value. Other holidays, such as Easter and the Labour Day, result in similar effects. During the year, IGD shows a decrease from January to August and an increase from September to December, due to the use of gas for environmental heating. All these features can be appreciated in Fig. 2, where 11 years of IGD are superimposed, aligning weekdays to better highlight periodic behaviours.

The periodogram, plotted in Fig. 3, exhibits peaks at periods of 365.25 and 7 days, while other relevant values are ascribable to multiple harmonics of the fundamental ones. Notably, differently from what happens for RGD [15], the weekly seasonality prevails on the yearly one in terms of magnitude.

Temperature is known to be a major determinant of gas demand [3, 5, 6]. In order to take into account that the need for heating ceases when temperature raises above $18\,^{\circ}{\rm C}$, it is useful to refer to the so-called Heating Degree Days (HDD), defined as $\mathrm{HDD}=\max(18-T,0)$, where $T$ is the temperature in degrees Celsius. The scatter plots of IGD against temperature and IGD against HDD are reported in Fig. 4.

3.2 Thermoelectric gas demand

Differently from IGD, Thermoelectric Gas Demand (TGD) shows a clear trend, see Fig. 1. From 2008 to 2014 TGD decreases, mostly due to the growing importance of renewable sources of electric power, while, since 2014, the trend stabilises, likely due to the decrease in subsidies to the installation of photovoltaic systems.

Moreover, TGD shows a greater variability compared to IGD and RGD, as seen from the year-over-year plot in Fig. 5. TGD is indeed influenced by several factors, including prices of electric power, gas, and European Emission Allowance (EUA) certificates, which exhibit a large volatility [19]. This explains why yearly periodicity is relatively less important in TGD than in IGD and RGD. The periodogram in Fig. 6 shows that, also for TGD, the main seasonal component is the weekly one, which is consistent with Italian power demand [20].

The scatter plot of TGD against temperature, displayed in the left panel of Fig. 7, shows a peculiar U-shaped pattern: TGD increases as weather gets colder, but also when it gets hotter. In fact, in summer more thermoelectric production is required because air conditioning pushes the demand for electric power. This U-shaped pattern justifies the introduction of a suitable feature variable, herein named Heating and Cooling Day Degrees (HCDD). More precisely

[TABLE]

We found that $T_{c}=16°C$ maximises the linear correlation between TGD and HCDD.

4 Feature extraction

Based also on the exploratory analysis, the features used in the prediction algorithms include: autoregressive terms, calendar features, temperature and its derived variables HDD and HCDD. Table 1 reports the complete list of features.

4.1 Autoregressive features

In order to capture the yearly periodicity and the concurrent effect of holidays, it is convenient to resort to the notion of similar day, see e.g. [15], here recalled for convenience. Letting $t$ and $\tau$ be calendar dates, we define:

• •

$\mathrm{year}(t)$: the year to which day $t$ belongs;

• •

$\mathrm{weekday}(t)$: the weekday of day $t$, e.g. Monday, Tuesday, etc;

• •

$\mathrm{yearday}(t)$: the day number within $\mathrm{year}(t)$ starting from January 1, whose yearday is equal to $1$.

Definition 1** (Similar Day [15]).**

If $t$ is not a holiday222According to the Italian calendar, holidays are: 1 January, 6 January, 25 April, 1 May, 2 June, 15 August, 1 November, 8, 25 and 26 December, Easter and Easter Monday., its similar day $\tau^{*}=\mathrm{sim}(t)$ is

[TABLE]

subject to

• •

$\mathrm{year}(\tau)=\mathrm{year}(t)-1$;

• •

$\mathrm{weekday}(\tau)=\mathrm{weekday}(t)$;

• •

$\tau$* is not a holiday.*

If $t$ is a holiday, its similar day $\tau^{*}=\mathrm{sim}(t)$ is the same holiday in the previous year.

To predict $y_{t},y\in\{\mathrm{IGD},\mathrm{TGD}\}$, we included as features $y_{t-1}$, $y_{t-7}$, $y_{\mathrm{sim}(t)}$ and $y_{\mathrm{sim}(t-1)}$.

4.2 Calendar features

As calendar features, we introduced binary dummy variables to account for weekdays and holidays. Dummy variables were also added to identify: (i) extended holidays, i.e. working days preceded and followed by either Saturdays, Sundays or holidays, and (ii) days after holidays, i.e. working days which immediately follow a holiday and are not extended holidays.

4.3 Temperature features

We selected as features also forecasted temperatures $T_{t}$, $T_{t-1}$, $T_{t-7}$ and $T_{\mathrm{sim}(t)}$. In view of what shown in Section 3, for IGD, also HDD values at the same times were introduced, while, for TGD, HCDD replaced HDD.

5 Predictive models

We selected nine base models, which can be grouped into three categories:

linear models: ridge regression, lasso, Torus model [20], support vector regression, and elastic net 2. 2.

non-linear models: random forest, neural networks 3. 3.

non-parametric models: Gaussian Process, nearest neighbour

Four of them, namely ridge regression, Gaussian Process (GP), Torus model, nearest neighbours and neural networks, were already applied to RGD: we thus refer to [15] for their description.

Moreover, we introduced four ensemble models, described in Section 5.2, which aggregate forecasts issued by a subset of the basic models: (i) simple average, (ii) average on an optimised subset of base forecasts (subset average), (iii) weighted average, (iv) support vector regression.

5.1 Base models

For the design of the predictors, we assume that $n$ training data pairs $(\textbf{x}_{i},y_{i})$, $i=1,\ldots,n$, are used to devise a prediction rule $f(\cdot)$. In particular, letting $(\textbf{x}_{*},y_{*})$ denote a novel input-output pair, $f(\textbf{x}_{*})$ will be used to predict $y_{*}$. Herein, $\textbf{x}_{i}\in\mathbb{R}^{p}$, $p<n$, is a vector whose entries are given by the $p$ features associated with the target $y_{i}$.

In the following, $\textbf{y}=y_{i}\in\mathbb{R}^{n}$ is the vector of the training target data, while $\textbf{X}=\{x_{ij}\}\in\mathbb{R}^{n\times p}$ is the matrix of the training input data, where $x_{ij}$ is the $j$-th feature of the $i$-th training pair $(\textbf{x}_{i},y_{i})$.

Linear models

Ridge regression [21], LASSO [22] and elastic net [23] are methods to identify the parameters $\beta_{j}$ of the linear-in-parameter predictor:

[TABLE]

where \mathbf{x}^{T}=\left[\begin{array}[]{cccc}x_{1}&x_{2}&\ldots&x_{n}\end{array}\right]. Accordingly, the vector of the predicted training targets is

[TABLE]

To prevent overfitting and improve generalization capabilities, in all the three methods the loss function includes a penalty on the magnitude of $\boldsymbol{\beta}$:

[TABLE]

The three methods share the same standard quadratic loss

[TABLE]

but use different penalties that result in specific shrinking patterns. In ridge regression, the quadratic penalty shrinks parameters toward the origin. In the LASSO the penalty on the sum of absolute values has the effect of zeroing the least relevant parameters, thus enforcing some degree of sparsity. In the elastic net an intermediate effect is achieved.

Both $\lambda$ and $\alpha$ play the role of hyperparameters: $\lambda$ controls the strength of the parameter shrinkage, while $\alpha$ tunes the balance between penalties on the $L_{1}$ and the $L_{2}$ norm of $\boldsymbol{\beta}$.

An alternative method to estimate $f$ is by means of support vector regression (SVR), see e.g. [24]. In this case, in place of the quadratic loss, an $\epsilon$-insensitive loss function is used:

[TABLE]

Moreover, the assumption is made that $f\in\mathcal{H}$, where $\mathcal{H}$ is a Reproducing Kernel Hilbert Space (RKHS) [25], whose reproducing kernel is denoted by $K(\cdot,\cdot)$. Under this assumption,

[TABLE]

where ${\alpha_{i}}$ are such that

[TABLE]

The SVR estimate is defined as

[TABLE]

The hyperparameters are the real-valued constants $\lambda$ and $\epsilon$.

Although $f^{\text{SVR}}(\mathbf{x})$ is a nonlinear function of $\mathbf{x}$, the Representer Theorem (see e.g. [26]) ensures that there exist coefficients $c_{i}$ such that the predictor function can be written as a linear combination of kernel functions centered at $\mathbf{x_{i}}$

[TABLE]

so that also the SVR predictor, though implementing a nonlinear function of the features, has a linear-in-parameter structure.

Nonlinear models: Random forest

The Random forest method (see e.g. [27]) is based on so-called Classification and Regression Trees (CART) [28]. CARTs perform a recursive feature-wise partitioning of the input space and fit local linear regressions in each region of the final partition. CARTs are known to be unstable and prone to overfitting. In order to overcome these limitations, random forest models grow multiple CARTs, resorting to so-called data and feature bagging. Bagging or bootstrap aggregating is a random selection of a subset of a dataset that is repeated multiple times (with replacement). Models are then trained on each selected subset. By applying bagging to both data and features, each tree gets trained on different samples and feature sets. Forecasts performed by all models are then averaged to get the final prediction, leading to a more stable model.

5.2 Ensemble models

Four alternative aggregation techniques were considered. Apart from the Simple average, the calibration of the other three methods requires a specific ensemble training dataset not used for training the base models.

Simple average. The most trivial aggregation is the arithmetic average of the forecasts achieved by base models. Given a test input $\mathbf{x}_{*}$, and $M$ base forecasts $\hat{f}_{i}(\mathbf{x}_{*})$, $i=1,...,M$, the ensemble forecast is

[TABLE]

Weighted average. A second option is the weighted average of base forecasts:

[TABLE]

The weights $w_{i}$ are obtained by minimizing the sum of squared residuals between the ensemble forecast and the target vector on the ensemble training dataset.

Subset average. The third ensemble method computes the average of a suitable subset of predictors. The chosen subset is obtained by a brute force search within the set of all possible subsets, choosing the subset of predictors whose average minimises the MAE computed on the ensemble training dataset. In our case, excluding the complete subset made of all the nine predictors (already considered as simple average), and the nine base models, the number of candidate subsets is

[TABLE]

SVR aggregation. The fourth ensemble method trains a SVR model on the ensemble training dataset, using base forecasts as features.

6 Implementation notes

The available data range from 2007 to 2018. Four one-year long test sets, ranging from 2015 to 2018, were used to obtain a comparative assessment of the 13 models, including 9 base models and 4 ensemble ones. Each test set was associated to a set of training data, that were organised differently depending on the nature of the considered model, either base or ensemble.

Training of base models. The training set, called base training set $\mathcal{T}_{\mathrm{base}}(\mathcal{Y})$, is made of all data previous to the test year $\mathcal{Y}$. For instance, if $\mathcal{Y}=\{2017\}$ is taken as test set, the 9 base models were trained on the base training set $\mathcal{T}_{\mathrm{base}}(\{2017\})=\{2007,\ldots,2016\}$.

Training of ensemble models. In this case two training sets were considered. The year before the test set $\mathcal{Y}$ was used as ensemble training set $\mathcal{T}_{\mathrm{ens}}(\mathcal{Y})$, while the remaining data were used to train the 9 base models that enter the aggregation. For instance, if $\mathcal{Y}=\{2017\}$ is taken as test set, the 9 base models were trained on $\mathcal{T}_{\mathrm{base}}(\{2016\})=\{2007,\ldots,2015\}$, while the ensemble models were trained on $\mathcal{T}_{\mathrm{ens}}(\{2017\})=\{2016\}$.

Hyperparameters of the Torus model were tuned by maximising AIC, those of the Gaussian Process by maximising the marginal likelihood, while for all the other base models five-fold cross validation was used.

Test of base models. For the test set $\mathcal{Y}$, base predictions were computed using the base models trained on $\mathcal{T}_{\mathrm{base}}(\mathcal{Y})$.

Test of ensemble models. Given the forecasts provided by the base models trained on $\mathcal{T}_{\mathrm{base}}(\mathcal{Y})$, the ensemble forecasts were obtained using ensemble models trained on $\mathcal{T}_{\mathrm{ens}}(\mathcal{Y})$.

Out-of-sample performances were evaluated in terms of Mean Absolute Error (MAE):

[TABLE]

where $y_{t}$ and $\hat{y}_{t}$, $y\in\{\mathrm{RGD},\mathrm{IGD},\mathrm{TGD},\mathrm{GD}\}$, are the actual value and its forecast, while $n$ is the number of samples in the considered test set. MAE was preferred over percent or relative error metrics due to the large range of values assumed by the target variable, which would give undue importance to poor performances during low-demand periods.

7 Results

Mean absolute errors for RGD, IGD, TGD, and total GD are reported, respectively, in Table 2, Table 3, Table 4 and Table 5.

For what concerns base models, GP, ANN and SVR achieved the best average MAE across all the gas demands, with differences between each other smaller than 0.10 MSCM. Notably, results achieved by such models were also stable across different test sets. On the other hand, KNN was consistently the worst performer, due to its poor capability of modelling influence of temperature and holidays.

Ensemble models consistently outperformed base ones. In particular, subset average achieved the best average MAE on all four types of gas demand: the three disaggregated demands and the total one. A possible explanation is that different models are better at capturing specific behaviours: in [15], for instance, it was shown that the ANN model achieved the best results in winter, while GP in summer, suggesting that the former is better at modelling the impact of weather, while the latter can better follow seasonal patterns. Aggregation can indeed mitigate errors committed by single models, thus increasing overall accuracy and robustness.

The improvement due to aggregation was particularly evident for RGD (Table 2), where the best base model (GP) was outperformed by the best ensemble model (subset average) by 0.25 MSCM. The gap between base and ensemble models was smaller for the other two gas demands: GP and SVR are worse than subset average by 0.07 MSCM, for IGD (Table 3); SVR is worse than subset average by 0.09 MSCM, for TGD (Table 4). Finally SVR is worse than subset average by 0.29 MSCM, for the global Italian GD (Table 5). The 2018 forecasts and the corresponding residuals provided by the best ensemble predictor, namely subset average, are displayed in Fig. 8 and Fig. 9, respectively.

To the best of our knowledge, the only term of comparison available for the task addressed in this paper is given by the forecasts of the global Italian GD issued by SNAM Rete Gas, the Italian Transmission System Operator (TSO) [29]. In 2017 and 2018, the improvement is neat: the out-of-sample MAE of SNAM predictions was 9.62 MSCM in 2017 and 8.30 MSCM in 2018, while our best model (subset average) scored 5.16 MSCM in 2017 and 5.46 MSCM in 2018, see Table 5.

8 Conclusions

We analyzed the industrial and thermoelectric components of Italian daily gas demand, completing a previous study concerning residential demand. We found that industrial and thermoelectric demand show different relationships with temperature and crafted features to properly take them into account.

Several forecasting models were investigated and compared: nine base models plus four ensemble models. Aggregated models were found to be consistently more effective than base ones. In particular, in 2017 and 2018 the best ensemble model, i.e. subset average, outperformed forecasts provided by the Italian TSO.