Comparing the Simple to Complex Automatic Methods with the Ensemble Approach in Forecasting Electrical Time Series Data

Abstract

The importance of forecasting in the energy sector as part of electrical power equipment maintenance encourages researchers to obtain accurate electrical forecasting models. This study investigates simple to complex automatic methods and proposes two weighted ensemble approaches. The automated methods are the autoregressive integrated moving average; the exponential smoothing error–trend–seasonal method; the double seasonal Holt–Winter method; the trigonometric Box–Cox transformation, autoregressive, error, trend, and seasonal model; Prophet and neural networks. All accommodate trend and seasonal patterns commonly found in monthly, daily, hourly, or half-hourly electricity data. In comparison, the proposed ensemble approaches combine linearly (EnL) or nonlinearly (EnNL) the forecasting values obtained from all the single automatic methods by considering each model component’s weight. In this work, four electrical time series with different characteristics are examined, to demonstrate the effectiveness and applicability of the proposed ensemble approach—the model performances are compared based on root mean square error (RMSE) and absolute percentage errors (MAPEs). The experimental results show that compared to the existing average weighted ensemble approach, the proposed nonlinear weighted ensemble approach successfully reduces the RMSE and MAPE of the testing data by between 28% and 82%.

1. Introduction

Electricity has an important role in human life. Therefore, research in this field is continually developing, from the development of renewable energy sources, and energy efficiency to mathematical modeling. Time series data related to the energy sector, such as electricity load, supply, and consumption, have a unique pattern and challenge researchers to produce methods with accurate forecasting values. There has been extensive discussion in [1,2,3,4,5,6,7,8,9,10,11] of the conventional model approach for forecasting electrical time series data, including the ARIMA (autoregressive integrated moving average) and exponential smoothing. Meanwhile, numerous authors [12,13,14,15,16,17,18,19,20,21] explored machine learning approaches, including Prophet and neural networks (NN). The combination of the two approaches has also been widely explored (see [22,23,24,25,26]).

In electric load forecasting, regression and ARIMA hybrid approaches have succeeded in representing the behavior of linear trends and the influence of holidays [27]. Meanwhile, for more complex patterns, the study in [28] introduces the TLSNN (two-level seasonal neural network) and the TLCSNN (two-level complex seasonal neural network) models, which take advantage of SSA (singular spectrum analysis) to decompose trend and seasonal patterns so that a more appropriate trend function can be chosen to represent other nonlinear trend patterns. Furthermore, in the TLCSNN, the seasonal function is approximated by a combination of amplitude-modulated sinusoidal functions to capture non-integer seasonal periods [29,30]. At the same time, NN (neural networks) can be used to handle nonlinearity relationships in the data.

Holt and Winter [31,32] introduced exponential smoothing as a simple method for capturing trend and seasonality patterns in time series data. This method has been a competitive alternative to the seasonal ARIMA and is widely used in applications (see [18,33,34,35,36,37,38,39,40]) because of its robustness and accuracy [41]. An overview of the development of the exponential smoothing model for time series with trend and seasonal patterns can be seen in [3]. This method has attracted researchers’ attention in electrical load data forecasting. Taylor [42] developed exponential smoothing for accommodating the trend and multiple seasonal patterns in load time series. Another study [43] developed an exponential smoothing method with a state–space approach that can model high-frequency electrical load data well. Later, another study [44] proposed a hybrid exponential smoothing-neural network and applied the method to the high-frequency Indonesian load time series. Exponential smoothing has become more popular since the study in [45] combined this method with a neural network and this took first place in the M4 competition [46,47].

The development of machine learning algorithms has encouraged researchers to propose automatic methods that make them easy to use. Users do not need to understand in-depth knowledge of statistics to be able to implement and run them. Many of these methods are old methods with a new twist. For example, the Prophet model developed by the Facebook (FB) team [48] has fundamentally the same idea as the simple model built by combining trend, seasonal, and irregular components discussed in [49]. The difference is in the more flexible trend and seasonality functions.

Several packages in R and Python offer automated models for time series forecasting. These include the ARIMA, exponential smoothing, NN, and Prophet models. Even though the model is built automatically based on specific criteria, such as minimizing AIC, RMSE, or other error evaluation tools, in some cases, the resulting model may not provide forecasting accuracy as expected. To address this issue, several researchers have developed ensemble models. An ensemble model is built from a pool of models [50] or the same model but with different parameters [51,52]. An ensemble model can also be considered to solve the problems resulting from the model, data, and parameter uncertainty.

As stated in [9,51,52], the resulting forecast of ensemble models can be calculated using the mean, median, mode, and trimmed mean techniques. These methods are simple but ignore the possibility of differences in the weights of each base models. The following year, the studies in [53] and [54] proposed the weighted ensemble method based on a nonlinear optimization and the artificial bee colony algorithm, respectively. Both [53,54] considered absolute error to obtain the weights of the base models. According to [55,56], absolute error is less sensitive to the outlier than the square error. Meanwhile, in electricity load time series, not all the extreme values can be ignored or be considered as outliers that represent corrupted parts of the data. This could be related to the community habits that occur periodically, such as the “mudik” culture in every Idul Fitri in Indonesia, which causes the load demand to decrease significantly compared to other days [28]. On the other hand, in defining an ensemble model, it is necessary to pay attention to the selection of both the base model’s appropriateness and the weighting of the base model [52,54]. Inspired by [52,53,54], this work proposes alternative weighting strategies in the ensemble model which are expected to improve forecasting accuracy by considering a pool of single methods appropriate to the time series, with trend and seasonal patterns found in electrical time series.

The motivation and contribution of this study are summarized as follows. First, considering the development of the automatic models, i.e., the auto ARIMA, the ETS method, the DSHW method, the TBATS model, NN, and Prophet models, this study evaluates the implementation of those automatic models to forecast the electricity time series data, which usually show trend and seasonal patterns. Second, in handling the problem in model selection for the historical data of the electricity time series as the source of model uncertainty, this work proposed two weighted ensemble approaches that combine the six single constituent automatic models linearly (called EnL) and nonlinearly (called EnNL), respectively, by minimizing square errors. To show the effectiveness and applicability of the proposed ensemble approach, four time series datasets with different behaviors and different characteristics, i.e., the monthly electricity net generation in the US, hourly household electricity demand in Ontario, Canada, half-hourly electricity demand in England and Wales, and half-hourly electricity demand in Victoria, Australia are considered.

The organization of this paper is as follows. Section 2 provides a brief overview of the methods used to model the time series data and also the procedure of the proposed ensemble approaches. The quantitative analysis results and the findings are summarized and discussed in Section 3. Finally, the conclusion and future directions are presented in Section 4.

2. Materials and Methods

This section provides a brief overview of the ARIMA, exponential smoothing including the ETS (error–trend–seasonal) method, the DSHW (double seasonal Holt–Winter) method, the TBATS (trigonometric, Box–Cox transformation, ARIMA error, trend, and seasonal) model, and two machine learning methods, Prophet and NN.

2.1. The ARIMA

The ARIMA is probably the most classic and popular method for time series analysis and forecasting. For the seasonal time series, the model can be represented as in Equation (1)

$${\varphi }_P\left(B^s\right){\psi }_p\left(B\right){\left(1-B\right)}^d{\left(1-B^s\right)}^D{Y}_t={\theta }_q\left(B\right){\eta }_Q\left(B^s\right){ϵ}_t,$$

(1)

and denoted as the ARIMA($p,d,q$)x${\left(P,D,Q\right)}_s$, where $s$ is the seasonal period [49,57]. The regular autoregressive and moving average factors are notated as ${\psi }_p\left(B\right)$ and ${\theta }_q\left(B\right)$, while the seasonal autoregressive and moving average are represented by ${\varphi }_P\left(B^s\right)$ and ${\eta }_Q\left(B^s\right)$, where $p, q, P$, and $Q$ indicate the order of the model. The number of regular and seasonal differencing needed to stationarize the time series is notated by $d$ and $D$, respectively. When $P=D=Q=s=0$, then the model presented in Equation (1) becomes a regular ARIMA (with no seasonal part).

Hyndman and Khandakar [57] proposed an algorithm named auto.arima, which estimates the ARIMA model parameters by considering unit root tests, MLE, and AICc. The default of this algorithm includes the determination of $d$ (between 0 and 2) using the repeated KPSS test, the stipulation of the order of the model by minimizing the AICc, and the variation in the defined model by adding or reducing the order of the model until no lower AICc is found. The chosen model is the model with the smallest AICc.

2.2. Exponential Smoothing

This work considered three automatic exponential smoothing methods provided in the R software package forecast version 8.20, i.e., the ETS method, the DSHW method, and the TBATS model. Exponential smoothing for trend and seasonal components was initially developed in [31,32], named the Holt–Winter (HW) method. The trend and seasonal patterns can be taken into account in the model either as additive or multiplicative components. The additive HW method assumes that the trend and seasonal component have a linear relationship in the model, while the multiplicative HW method considers that the input, i.e., the past observations, which are presented as a trend, and seasonal components follow a nonlinear relationship to the output. This means that the additive HW is recommended when the seasonal variation tends to be constant, while the multiplicative HW is considered when there is a change in the seasonal variation proportional to the level of the series.

The ETS method is exponential smoothing with the state–space approach, which has an automatic forecasting procedure proposed in [58]. In this case, the level, trend, and seasonal components as states change over time. The model notated by ETS (“Error”, “Trend”, “Seasonal”) method with the possibility for “Error” is additive (A) or multiplicative (M). The “Trend” is None or No trend (N), Additive (A), or Additive with a Damped parameter, while “Seasonal” is None or No seasonal (N), Additive (A), or Multiplicative (M). These state–space models are more general for the exponential smoothing model. In the analysis, the error, trend, and season components are automatically selected and the best model is chosen based on AICc. Further discussion about this method can be found in [59].

The DSHW method was developed by [42] to accommodate two seasonality that cannot be handled by the standard HW. This model includes two seasonal periods where the longer seasonal period is a multiple of a short one or in other words, there is a repetition of a short seasonal effect within the longer one [42]. The best model is fitted by minimizing the mean square error using dshw() function in the R software.

The TBATS model is more general than the DSHW method. It is also useful for modeling high-frequency time series with trend and double seasonal patterns. The seasonal period is not only limited to the integer value but can also be the non-integer value [43].

2.3. Prophet

The Prophet model developed by the Facebook team consists of trend, seasonality, and holiday components that can be written as

$$Y_t=T_t+S_t+H_t+{ϵ}_{t.}$$

(2)

where $T_t$, $S_t$, and $H_t$ represent the trend, seasonal, and holiday effect, respectively [48]. The error term, ${ϵ}_t,$ is usually assumed to be normally distributed. It accommodates any changes that are not taken into account in the model. The flexibility of the Prophet model represented in Equation (2) makes it work well for time series with trends, strong seasonal effects, and several seasons of historical data. Prophet is robust to missing data and shifts in the trend and typically handles outliers well [48].

In the analysis, a linear function and the Fourier series are employed to model the trend component and to represent the periodic component, respectively. The piecewise linear model is represented as

$$T_t={\beta }_0+{\beta }_{1}t+{\beta }_2{\left(t-c\right)}^{+}$$

(3)

where $t$ shows the time steps and $c$ is the value of break point while ${\beta }_j, j=0, 1, 2$ are unknown parameters. Equation (3) can also be written as Equation (4).

$$T_t=\left\{\begin{array}{ll}{\beta }_0+{\beta }_t\hspace{1em}\hspace{1em}& t\le c\\ {\beta }_0-{\beta }_{2}c+\left({\beta }_1+{\beta }_2\right)t& t>c\end{array}\right.$$

(4)

Meanwhile, the seasonal effect is approximated by Fourier series as in Equation (5).

$$S_t={\displaystyle \sum }_{i=1}^{n_F}{a}_i\cos \left(\frac{2\pi it}{s}\right)+b_i\sin \left(\frac{2\pi it}{s}\right)$$

(5)

where $a_i$ and $b_i$ for $i=1, 2,\dots , n_i$ are unknown parameters, $s$ is the seasonal period and $n_F$ is the order of Fourier series. In this work, the break points are determined automatically by the Prophet algorithm with the change prior scale is considered to be 0.05 and the holiday effect ($H_t$) is ignored, as in other alternative models used in this study.

2.4. Neural Networks

The neural network algorithm discussed in this work is a feedforward neural network whose inputs are lagged values of time series and the output is the forecast values. Differently from the ARIMA, it enables modeling the nonlinear relationship between input and output variables in the model. The output of the NN can be obtained by Equation (6).

$$Y_t=v_o+{\displaystyle \sum }_{j=1}^{n_h}{v}_{j}f\left(ne{t}_j\right)$$

(6)

where $v_j, j=0, 1, \dots , n_h$ are the unknown parameters and $f$ is a sigmoid function defined as

$$f\left(ne{t}_j\right)=\frac{1}{1+\exp \left(ne{t}_j\right)}$$

(7)

and

$$ne{t}_j=w_{0j}+{\displaystyle \sum }_{k=1}^{n_k}{w}_{kj}{Y}_{t-k}+{\displaystyle \sum }_{l=1}^{n_l}{w}_{lj}{Y}_{t-ls}.$$

(8)

$w_{0j} \left(j=1, \dots , n_h\right)$ are the bias feed into the $j^{th}$ node of the hidden layer, $w_{kj} \left(k=1, 2, \dots , n_k; j=1, 2, \dots , n_h\right)$ are the weights connected the input $Y_{t-k}$ to the node j in the hidden layer, and $w_{lj} \left(l=1, 2, \dots , n_l; j=1, 2, \dots , n_h\right)$ are the weights connected the input $Y_{t-ls}$ to the node j in the hidden layer, where s is the seasonal period. The terms $n_k, n_l,$ and $n_h$ represent the number of regular lagged input, seasonal lagged input, and hidden nodes, respectively. The values of $n_k$ and $n_l$ in Equation (8) are selected automatically according to the AIC value, while $n_h$ is set to the nearest integer of $\left(n_k+n_l+1\right)/2$ [60] using nnetar() function in the R software.

2.5. Proposed Ensemble Methods

The proposed ensemble methods include six steps.

Step 1:

Divide the time series data into two parts, training and testing data. Let $\left\{Y_t, t=1, \dots , N\right\}$, then the training dataset is $\left\{Y_1, Y_2, \dots ., Y_{N-24}\right\}$ and the testing dataset is {$Y_{N-24+1}, Y_{N-24+2}, \dots , Y_N\}$, where $N$ is the sample size.

Step 2:

Model the train data using automatic models, i.e., the ARIMA, the ETS method, the DSHW method, the TBATS model, NN, and Prophet.

Step 3:

Calculate the forecast values up step h ahead, where, in this work, $h=1, 2, \dots , 24$ and is the time horizon for the testing data. The forecast values, ${\widehat{Y}}_{N-24+h}^{\left(i\right)} \mathrm{for} i=1, 2, \dots , 6,$ are obtained by the models listed in step 2.

Step 4:

Assemble the models with two weighting strategies

• based on the linear relationship function, named EnL,

  $${\widehat{Y}}_{N-24+h}={\displaystyle \sum }_{i=1}^6{w}_i{\widehat{Y}}_{N-24+h}^{\left(i\right)}$$

  (9)

  where $w_i$, $i=1, 2, \dots , 6$, are obtained by minimizing square error, which can be written as in Equation (10).

  $$\min g\left(w\right)={\displaystyle \sum }_{t=1}^{N-24}{\left({\displaystyle \sum }_{i=1}^6{w}_i{\widehat{Y}}_t^{\left(i\right)}-Y_t\right)}^2$$

  (10)

  with linear constraint ${\displaystyle \sum }_{i=1}^6{w}_i=1$.
• based on the nonlinear relationship function, named EnNL

  $${\widehat{Y}}_{N-24+h}=v_0+{\displaystyle \sum }_{j=1}^{n_h}{v}_{j}f\left(w_{0j}+{\displaystyle \sum }_{i=1}^6{w}_{ij}{\widehat{Y}}_{N-24+h}^{\left(i\right)}\right)$$

  (11)

  where $f$ is a sigmoid function as written in Equation (7). The weights $v_0, v_j, w_{0j},$ and $w_{ij}$ for $i=1, 2, \dots , 6; j=1, 2, \dots , n_h$ are estimated by an NN approach. The network is trained for input–target pairs, i.e., $\left({\widehat{Y}}_t^{\left(1\right)}, {\widehat{Y}}_t^{\left(2\right)}, \dots , {\widehat{Y}}_t^{\left(6\right)}\right)$ and $Y_t, t=1, 2, \dots , N-24,$ using a backpropagation algorithm with the Levenberg–Marquart optimization. The number of $n_h$ was set between 1 and 5 and selected based on the minimum RMSE while the maximum number of the epoch was 5000.

Step 5:

Calculate the forecast values using Equation (9) for the EnL method and Equation (11) for the EnNL method.

Step 6:

Evaluate the model based on RMSE and MAPE using Equation (12) and Equation (13), respectively.

$$\mathrm{RMSE}={\left(\frac{{{\displaystyle \sum }}_{t=1}^{n_{obs}}{\left(Y_t-{\widehat{Y}}_t\right)}^2}{n_{obs}}\right)}^{1/2},$$

(12)

$$\mathrm{MAPE}=\frac{1}{n_{obs}}{\displaystyle \sum }_{t=1}^{n_{obs}}\left|\frac{Y_t-{\widehat{Y}}_t}{Y_t}\right|\times 100\%,$$

(13)

where $n_{obs}$ is the number of observations included in the calculation. In this study, the two proposed weighted ensemble approaches are compared with an existing average weighted ensemble and six single models including the ARIMA, the ETS method, the DSHW method, the TBATS model, NN and Prophet. Those six single models are considered capable of modeling time series with trend and seasonal data as in electrical time series data. Clearly, the research stages are presented in Figure 1.

The experiment consists of modeling four datasets by six single automated and three ensemble methods and computing 1 to 24 steps ahead, forecasting using all constructed models. All models are fitted on a Hewlett-Packard laptop manufactured from China with intel^® core i-7 4.70 GHz processor with 32 GB of RAM and running R 4.2.2 for the single automated models and Matlab 2023a for the ensemble models.

3. Results and Discussion

In this section, four time series data with different seasonal periods, i.e., monthly, hourly, and half-hourly, are presented as illustrative examples. The characteristic of each dataset is described in

Table 1. Characteristics of the four datasets discussed in this study.

Data	Name	Frequency	Number of Observations in the Training Data	Characteristic
1	usmelec	Monthly	317	Linear trend, increasing seasonal variation
2	ontario	Hourly	3648	Nonlinear trend, relative constant seasonal variation
3	taylor	Half-hourly	2354	No trend, constant seasonal variation
4	vic_elec	Half-hourly	2354	Fluctuate in level and seasonal variation

.

These four datasets are considered in order to show the effectiveness and applicability of the proposed ensemble approach for electricity load forecasting.

3.1. Data

The four time series datasets are described below.

Data 1:

US Monthly Electricity Total Net Generation

The first time series is the monthly total net generation of electricity from January 1973 to May 2001 presented in Figure 2. These data can be accessed freely from the library “fpp” in the R software and are named “usmelec” [61]. The data are split into two parts, i.e., the first 317 observations from January 1973 to May 1999 as the training data and the last 24 observations as the testing data.

Figure 2 shows that the data tend to increase linearly from time to time, with a change in seasonal variation.

Data 2:

Hourly Electricity Demand in Ontario

The second time series is hourly electricity demand in Ontario, Canada, from 1 August 2016, 00:00 to 31 December 2016, 23:00. These data were obtained from [62]. In the experimental study, the first 3648 observations (1 August 2016, 00:00 to 30 December 2016, 23:00) are considered as the training dataset and the rest (31 December 2016, from hour 00:00 to 23:00) as the testing set. The pattern of the data can be seen in Figure 3. Differently from Data 1, these data have a nonlinear trend. The change in seasonal variation is bigger in the beginning and more stable in the rest.

Data 3:

Half-Hourly Electricity Demand in England and Wales.

Data 3 were also discussed in [42] and can be accessed from the library “forecast” in the R software with the name “taylor”. We consider 2354 observations as the training data, i.e., data from 5 June 2000, 00:00:00 to 24 July 2000, 00:30:00. The test data consist of the last 24 observations (24 July 2000, from hour 01:00:00 to 12:30:00) [63]. Figure 4 shows that the data have no trend, with a stable seasonal variation.

Data 4:

Half-Hourly Electricity Demand in Victoria, Australia

A half-hourly electricity demand in Victoria, Australia, time series can be accessed freely from the library “tsibbledata” in the R software with the name “vic_elec” [64]. As shown in Figure 5, these data have a more complex pattern than Data 3, presented in Figure 4. In the analysis, we estimate the models using the 2354 first observations (12 November 2014, 00:00:00 to 31 December 2014, 00:30:00); and in the test, the last 24 last observations (31 December 2014, from 01:00:00 to 12:30:00).

3.2. Experimental Results

Firstly, all training data for each of the four time series are modeled using the ARIMA, the ETS method, the DSHW method, the TBATS model, NN and Prophet. The forecast values obtained from those six single models are then combined using the existing average weighted ensemble approach, and the two proposed weighted ensemble approaches, based on the linear relationship and nonlinear relationship function. The error evaluation in terms of RMSE and MAPE can be seen in

Table 2. Comparison of RMSEs and MAPEs for the training data obtained from the nine models.

	RMSE				MAPE
Model	Data 1	Data 2	Data 3	Data 4	Data 1	Data 2	Data 3	Data 4
ARIMA	6.57	237.78	1010.63	464.98	2.13%	0.96% 2	2.39%	7.35%
ETS	6.38	273.51	481.95	89.21	2.02%	1.29% 2	1.22%2	1.35% 2
DSHW	6.63	253.71	251.33	40.46	2.16%	1.13% 2	0.65% 2	0.67% 2
TBATS	6.40	246.49	229.59	43.99	2.03%	1.05% 2	0.60% 2	0.75% 2
NN	8.37	158.20	171.36	35.04	2.92%	0.69% 2	0.44% 2	0.61% 2
Prophet	8.45	1118.17	1467.91	333.53	2.93%	5.58%	4.16%	5.96%
EnA	6.17	261.69	359.83	119.79	1.98% 2	1.27% 2	0.95% 2	2.04%
EnL	6.04	157.78	177.46	29.98	1.94% 2	0.70% 2	0.46% 2	0.52% 2
EnNL	5.601	149.17 1	156.55 1	28.42 1	1.89% 1,2	0.67% 1,2	0.40% 1,2	0.49% 1,2

¹ The smallest values in each column. ² MAPE is less than or equal to 2%.

.

Table 2 shows that for the four datasets, the EnNL shows the smallest RMSE and MAPE compared to EnA and EnL, and all other single models. Meanwhile, EnL outperforms EnA. The ensemble model is worth recommending from MAPE’s viewpoint, even though EnA is slightly higher than 2% for Data 4, likewise for other automatic single models, apart from the ARIMA and Prophet. However, in selecting a model, further evaluation of the model’s performance on testing data is needed.

The performance of models in terms of RMSE and MAPE for the testing data are summarized in

Table 3. Comparison of RMSEs and MAPEs for the testing data obtained from the nine models.

	RMSE				MAPE
Model	Data 1	Data 2	Data 3	Data 4	Data 1	Data 2	Data 3	Data 4
ARIMA	9.61	859.75	5767.47	180.96	2.80%	4.85%	14.39%	3.32%
ETS	11.19	986.61	11,724.40	186.51	2.66%	5.59%	27.76%	3.56%
DSHW	9.87	1141.28	2094.31	590.19	2.74%	6.67%	5.35%	12.94%
TBATS	11.66	918.67	2903.67	103.29	2.73%	5.17%	7.14%	2.00% 2
NN	12.29	1160.34	2357.52	95.41	3.05%	6.52%	6.96%	1.98% 2
Prophet	16.27	619.86 1	732.69	146.06	4.08%	3.60% 1	1.95% 1,2	2.85%
EnA	10.30	904.52	3488.10	107.37	2.74%	5.16%	8.01%	2.22%
EnL	10.92	1138.94	644.24	228.03	2.75%	6.41%	2.03%	4.99%
EnNL	4.61 1	1166.62	627.41 1	76.55 1	1.11% 1,2	6.56%	2.12%	1.54% 1,2

¹ The smallest values in each column. ² MAPE is less than or equal to 2%.

. For Data 1 and Data 4, EnNL consistently provides the lowest RMSE and MAPE as in the training data (see Table 2). EnNL successfully modeled and forecasted Data 1, which have trend and seasonal patterns with changing variations. Supported by Figure 6, the MAPEs obtained from EnNL for all lead times up to twenty-four steps ahead are less than 2%. Within the discussion of Data 1, EnA and EnL are recommended for one-step-ahead forecasting.

As for Data 2, no model produces a MAPE value of less than 2% for up to twenty-four-steps-ahead forecast values (see Table 3). For this case, further development is needed, such as involving the holiday effect and selecting a more appropriate trend and seasonal function [28,29,65]. However, based on

Table 4. Comparison of MAPEs on the testing data of Data 2 with lead time forecast values up to seven steps ahead obtained from the nine models.

Model	MAPE
	h = 1	h = 2	h = 3	h = 4	h = 5	h = 6	h = 7
ARIMA	0.50%	0.63%	0.44%	0.53%	0.85%	1.43%	2.37%
ETS	1.03%	1.39%	1.43%	1.30%	1.10%	1.30%	1.91%
DSHW	0.13%	0.13%	0.42%	0.87%	1.48%	2.30%	3.39%
TBATS	0.55%	0.84%	0.81%	0.72%	0.93%	1.52%	2.66%
NN	0.76%	0.89%	0.74%	0.56%	0.65%	1.09%	2.07%
Prophet	2.50%	2.91%	3.11%	3.02%	2.49%	2.86%	3.66%
EnA	0.87%	1.09%	1.02%	0.77%	0.93%	1.49%	2.45%
EnL	0.86%	1.02%	0.90%	0.72%	0.75%	1.14%	2.07%
EnNL	0.75%	0.87%	0.72%	0.55%	0.64%	1.08%	2.05%

and Figure 7, all models except Prophet can be considered to forecast up to seven steps ahead. In Figure 7, it can also be seen that EnNL produces smaller MAPEs than those obtained from the other two ensemble models.

For Data 3, EnNL is the most recommended model based on the RMSE and MAPE values of the training data. Interestingly, based on Table 2, Prophet produces the smallest MAPE for the testing dataset and has a larger RMSE than EnNL. For further analysis, the MAPE behavior up to the next twenty-four half-hourly forecasts for all models is presented in Figure 8. The comparison of the actual values to forecast values up to 24 steps ahead using the Prophet method and the three ensemble approaches is also presented in Figure 9. Further, based on MAPE behavior (Figure 8) and forecasting values (Figure 9) for one to twenty-four steps ahead and considering RMSE and MAPE analysis on training data (Table 2) and testing data (Table 3), EnL and EnNL are more recommended compared to others including the Prophet. In this case, the accuracy of EnL and EnNL may still be improved by considering the influence of holidays on each constituent model.

Finally, Figure 10 shows that for Data 4, the EnNL, NN, and TBATS models produce relatively stable MAPE values of no more than 2%. Meanwhile, EnA tends to increase with increasing time ahead to be forecast. Even though the DSHW method provides relatively small RMSE and MAPE values on training data (Table 1), it cannot provide accurate forecasting values for the testing data (Table 2 and Figure 10). The ARIMA provides the largest RMSE and MAPE values on training data but produces the most accurate forecast for two to eight lead times and is less accurate after that (Figure 10).

According to the experiments, it can be summarized that most of the single automated models produced good-fitted values for the four different time series datasets. However, these models are better combined to minimize the forecast error because of the problem arising from the model uncertainty, that is mitigated by combining the results of different models in an ensemble framework. Furthermore, the experimental results show that the proposed EnNL approach significantly reduced RMSE and MAPE obtained from Data 1, Data 3, and Data 4 for up to twenty-four-steps-ahead forecasts compared to the EnA. It has reduced the RMSE and MAPE obtained from EnA by approximately 55%, 82%, 28%, 59%, 73%, and 30%, respectively. Further, for Data 2, the EnNL model forecasted up to seven steps ahead more accurately than EnA and EnL.

4. Conclusions

High-frequency electricity forecasting always become a challenging area to study. This study applied six single automatic models, i.e., the ARIMA, the ETS method, the DSHW method, the TBATS model, NN, and Prophet, and their combinations for modeling and forecasting time series data with trend and seasonal patterns encountered in electricity data. Automatic models are considered easily implemented, especially for users who do not have in-depth knowledge about the modeling procedure. Meanwhile, the ensemble models are regarded as handling the source of uncertainty both in the data and the model. In this case, we apply the automatic models provided by default by the package in R software without considering the influence of holidays. The forecast values obtained by the single automatic models were then combined based on three weighting strategies, i.e., averaging (EnA) and two proposed ensemble approaches considered the linearity (EnL) and nonlinearity relationships (EnNL), respectively. We compared the forecasting accuracy of each model based on RMSE and MAPEs. Experimental results of the four electricity time series data, i.e., US monthly electricity net generation, hourly household electricity demand in Ontario, half-hourly electricity demand in England and Wales, and half-hourly electricity demand in Victoria, show that the proposed ensemble approach, especially EnNL, provides better multistep-ahead forecast accuracy than those obtained from EnA, EnL, and other single models. In each case, it may differ for the number of future periods that can be predicted accurately.

In energy demand forecasting, understanding the character of data and various factors affecting the patterns is needed to develop the model. Involving the variable of the holiday effect and selecting the suitable weighted combination of the most suitable constituent models may considerably provide better performance accuracy.

Although the proposed ensemble approaches effectively enhance forecast accuracy, improvement can still be made. Finally, further investigations of this model development will become the focus of our future research.