**1. Introduction**

The programming of the production of electrical energy is a complex task carried out by the Transmission System Operators (TSO). Their main objective is the supply of electricity, assuring the distribution. In addition, the energy production at the lowest possible cost to the consumer, without incurring losses is a key point [1]. Any mismatch between the programmed and the really demanded energy produces huge losses. Hobbs [2] determined that the losses produced by a 1% gap between planning and reality can cost up to millions of dollars. Hong et al. [3] determined that for a maximum peak demand central of 1 GW, a 1% error in the prediction can involve costs of \$600,000 per year. The TSOs use for their predictions time series forecasting [4–6]. One of the most common used technique is the exponential smoothing methods, and especially the Holt-Winters models, due to their easy to understand and implement features [7,8]. Exponential smoothing methods use the data observed in the past to make predictions, assigning an exponentially decreasing weight to the older information against the newer. The way to assign the relevance of newer data against the older is performed by the weight assigned to smoothing parameters. The smoothing parameters are bounded in the range [0,1]. The closer to 0 the more important the older data is, and the contrary when closer to 1. The Holt-Winters model is described in Equations (1)–(4)

$$S\_{t-} = \alpha \left(\frac{X\_t}{I\_{t-s}}\right) + (1-\alpha)(S\_{t-1} + T\_{t-1})\tag{1}$$

$$T\_{t-1} = \gamma (S\_{t-1} S\_{t-1}) + (1 - \gamma) T\_{t-1} \tag{2}$$

$$I\_t = -\delta \left(\frac{X\_t}{S\_t}\right) + (1 - \delta)I\_{t-s} \tag{3}$$

$$\hat{X}\_{t+h} = \begin{array}{c} (\mathbb{S}\_t + hT\_t) I\_{t-s+h} \end{array} \tag{4}$$

where *St*, and *Tt* are the equations for the level and additive trend, with smoothing parameters *α* and *γ*. *It* are the seasonal indices of length *s*, with smoothing parameters *δ*. *Xt* is the observed data. Finally, *<sup>X</sup>t*+*<sup>h</sup>* is the equation to forecast h time instants ahead. It collects the information contained in the model and makes predictions.

The introduction of double and triple seasonal models (HWT), described in [9,10], improved the accuracy of these methods, and their use is being generalized. Additionally, these methods behave very accurate for such type of series, with a strong seasonal effect [11,12]. The generalization to include up to *n* seasonal patterns is known as Multiple Seasonal Holt-Winters (nHWT), developed in [13]. The model is defined as in Equations (5)–(8),

$$S\_{t-} = \alpha \left(\frac{X\_t}{\prod\_{i=1}^{n\_s} I\_{t-s\_i}^{(i)}}\right) + (1-\alpha)(S\_{t-1} + T\_{t-1})\tag{5}$$

$$T\_t = \gamma \dot{(S\_t - S\_{t-1})} + (1 - \gamma)T\_{t-1} \tag{6}$$

$$I\_t^{(i)} = \\\delta^{(i)}\left(\frac{X\_t}{S\_t \prod\_{j=1, j\neq i}^{n\_s} I\_{t-s\_j}^{(j)}}\right) + (1 - \delta^{(i)})I\_{t-s\_i}^{(i)} \quad i = 1, \ldots, n\_s \tag{7}$$

$$\hat{X}\_{t+h} = \quad (S\_t + hT\_t) \prod\_{i=1}^{n\_t} I\_{t-s\_i + h}^{(i)} + \varphi\_{AR} \varepsilon\_t \tag{8}$$

where *I* (*i*) *t* are the seasonal indices of length *si*, with smoothing parameters *δ*(*i*). There are as many equations as seasonal patterns allocated in the time series, *ns*. The *ϕAR* parameter is entered to correct the model including effectively the first order autocorrelation error (*<sup>ε</sup>t*), known as AR(1) adjustment.

Depending on the way the equations are defined, additive or multiplicative methods can be used. The method of combining the smoothing equations will determine the Holt-Winters model to be used. The combination provides 30 different models, according to the trend and seasonality combination, as well as AR(1) adjustment, as described in [13]. The nomenclature for these models is as follows: Three letters to describe the combination, where the first one describes the trend method, the second one the seasonality method and the last one whether the adjustment was applied. This combination is shown in Table 1. Furthermore, the model is described by adding as subscript the different seasonal patterns considered, using the length of the cycle. As an example, the *AMC*24,168 one is a model with additive trend and multiplicative seasonal, the model is adjusted using first order autocorrelation error and it has two seasonalities, one daily (subscript 24) and one weekly (subscript 168).

The smoothing parameters are obtained by adjusting the model to data observed in the past, and try to reproduce the same pattern for the near future. The adjustment is done by solving a non-linear problem, in which the 1-hour-ahead forecasting error is minimised. This topic is better developed in Section 3. The smoothing parameters obtained are assumed to be optimal, and the model is then exploited. The optimality of these parameters is locally obtained, valid only for the dataset and the model selected. Moreover, the same parameter can highly vary from one method to another using the same dataset. Thus, it is a big deal for forecasters to understand the behaviour of the model considering the locally optimised parameters. The future predictions accuracy will be closely related to the smoothing parameters. The optimisation procedure is using an optimisation algorithm, that do not take care concerning the reality of the series. The forecaster needs to make use of the own experience

to validate and approve the obtained parameters. These actions are crucial for the entire process. Some factors have an influence on the parameter values. The method to obtain initial values for seeding the model impact on the forecasting accuracy [14] although after an optimisation of the parameters, accuracy differences are also minimised [15]. Climate conditions are not a big deal in terms of short-term forecasting accuracy [12]. However, it seems reasonable the parameters ought to be influenced. If only two seasonal patterns are considered, seasons modify the trend, whereas if three seasonalities are considered, the intra-year seasonal parameter should be influenced. The calendar also has a huge influence on the predictions and parameters, that must deal with some irregularities of the series [16–18]. Therefore, forecasters should always pay close attention to the parameters values.

**Table 1.** Multiple seasonal Holt-Winters models' nomenclature according to trend and seasonal method. The first letter determines the trend, N: none, A: additive, d: damped additive, M: multiplicative and D: damped multiplicative. The second letter is used for Seasonality, with only N (none), A (additive) and M (multiplicative) option. The third letter is used for the AR(1) adjustment: L (no adjustment) and C (adjusted).


The forecasting procedure performed by TSOs must ensure continuously accurate forecasts for the electricity system. These predictions are used by the TSO for the operational planning and unit assignment, while they are also used by the market for the spot price settlement [19,20]. The Spanish electricity market (OMIE) operates similarly [21,22], where the unique Spanish TSO, Red Eléctrica de España (REE), supplies forecasts of demand for the next week every Wednesday, for the next 24 h every day, and a revision of the daily demand every six hour [23]. The market uses this information for the bidding process, as well as REE itself to carry out operational planning. Thus, there is an enormous responsibility in such eagerly awaited forecasts. In fact, REE uses a complex algorithm [24] to provide forecasts.

With such high frequency forecasting, there is a need for the model to be updated continuously. The forecaster must deal with the varying smoothing parameters and take decisions on the engaged forecasts. In the Holt-Winters literature there is a strong discussion about whether the parameters should be continuously readjusted as the series progresses, or on the contrary, they must be immobile and try to exploit them for as long as possible. The values of the parameters are also discussed, since high values of the parameters indicate a grea<sup>t</sup> damping to adapt to the new observed values, compared to low values, where the initial values are given greater weight. Before this doubt, the following question is added: how do the parameters respond to a specific time series, such as the hourly electricity demand in Spain? Make it sense the obtained values? What is their stability? All these questions have not been answered in the previous literature.

This article describes how we have performed a stability analysis on the parameters of the nHWT models applied to the electricity demand series in Spain. The main objective was to understand the behaviour of smoothing parameters against several different datasets with changes in climate conditions, time series components such as trend or seasonal, calendar effects among others. This analysis will help the forecasters to face the optimised model before making forecasts, considering the values of the parameters. It is also important to check whether it is necessary or not to update the smoothing parameters of the Holt-Winters models in order to obtain a good accuracy over time. On the other hand, the results of this study will show the influence of the calendar effect on these parameters, and therefore the need to develop new Holt-Winters models that take it into account [16].

The article is organized as follows: Section 2 reviews the existing literature related to the smoothing parameters of Holt-Winters models. Section 3 presents the methodology followed in order to predict the Spanish electricity demand time series and to analyze the variability of the parameters; in Section 4 the results obtained for a given double and triple seasonal model are shown and in Section 5 the results are discussed. Finally, the conclusions reached in this article are shown.

### **2. Related Work**

Short and medium term electricity demand time series forecasting were extensively studied in the literature applying both statistical methods and machine learning techniques. Recently, the machine learning methods have focused on big data [25], especially deep learning [26,27], and ensemble methodologies [28,29].

Within classical methods, the accuracy of the forecasts in the Holt-Winters models has been widely studied, especially how to select the best method and adjust the parameters [30]. As the Holt–Winters models are recursive, an initialization value is needed to feed the model. Thus, the initialization methods for Holt-Winters were also another intense field of research especially with the emergence of double and triple seasonal Holt-Winters models [31,32]. In fact, new methods to initialize the level, trend and seasonality in multiple seasonal Holt–Winters models were recently developed in [14]. These seed values have influence on the smoothing parameters' values as well as the forecasting accuracy [33]. However, when the series adjustment is made, the initial values lose influence [15]. After fitting the Holt-Winters model, the obtained smoothing parameters are assumed to be optimal for the data set. Thus, the study of their values or their stability has not been worked in depth in the literature. This is due, in large part, to the fact that the models do not allow a theoretical mathematical analysis.

Archibald [34] used 406 monthly series from the M competitions. He analyzed the models with additive seasonality, and demonstrated that the values of the parameters within the range [0,1] are not always invertible. Thus, it is necessary to use only a set range within the region of invertibility.

Lawton [35] used state spaces to analyze Holt-Winters models, and although his work focused on the normalization of the seasonal component, he also analyzed the stability of the parameters. From previous work on filters [36,37], it determined that Holt-Winters models are not asymptotically stable. The values of the eigenvectors of the stability matrices depend on the *α*, *γ* and *δ* parameters.

Some authors [38] worked on obtaining the limits of the smoothing parameters. They analyzed state space parameters, and stated the parameters ought to meet 0 < *α* < *γ*. Finally, the authors in [39] established a series of criteria in the parameters so that the models can be "predictable", a term that mints whether a series will be able to make forecasts with constant mean and variance. Osman et al. [40] obtained the values of the main vectors and confirmed that the models proposed in [39] produce stable predictions, but if regressors are used, they must be invariant over time.

Another interesting study related to this paper is the study of the minimum size required for the sample. Hyndman et al. [41] analyzed this situation and concluded that there is no clear answer, but that there should be many more observations than parameters. In this regard, García-Díaz and Trull [13] verified that for a double seasonal model, a data set length of one year produces more stable parameters than 8–10 weeks.

#### **3. Materials and Methods**

The time series used in this study is the hourly electricity demand in Spain, within the years from 2008 to 2017. It was provided by REE through its website www.ree.es. This time series shows clearly several seasonalities: the intraday which is repeated every 24 h, and the intra-weekly, which is repeated every 168 h. A third seasonality has also been included, which relates the data of the series with the seasonal periods of the year. Various strategies were proposed to deal with the 3rd seasonality, since its length depends largely on the concept to be covered. Taylor uses the calendar year, so that the years comprise 52 weeks, although to adjust the process in leap years, it uses 53 weeks [10]. Most authors have chosen to use the solar year, which includes 365.25 days, and seasonally adjusted [42]. In our case, it was considered to be 365.25 days, or 8766 h.

This series was used in former works [43,44] and especially using the nHWT methods [13,14,16]. This experience showed the authors that many factors could have an impact on the model parameters, as well as its forecast accuracy, key point for the proper exploitation of the nHWT models. The influence of the initial values was analysed in [14]. However, there is still a lack of analysis regarding the meaning of the parameter space and their relation with forecasts.

The method described in this paper analyses the smoothing parameters of the nHWT forecasting procedure when working with the hourly electricity demand in Spain. It performs an analysis of the parameter values, by checking their behaviour depending on the method chosen from Table 1. That implies to analyse the stability against some influential factors derived from our experience such as seasons, years, calendar, and inner variability.

To match the previously mentioned issues, raw data must be used, avoiding any smoothing or modification of the series. Double seasonal nHWT used a dataset size of 8–10 weeks, it was chosen as this seems to be long enough for this purpose. Alternatively, an analysis of the dataset size is also performed. Triple seasonal nHWT used a dataset of three years.

Randomness in the series selection must be assured. To observe the maximum possible variability, randomly selected data sets were used, but always complying with the following premises:


As a result, for the seasonal double analysis, 250 random samples of 10 weeks were extracted, of which, the first 8 weeks were used for adjustment, while the other two weeks were used for validation. Because the climatic effect present in the series must be assessed, the samples were selected so that they belong to different climatic periods and several years. Figure 1 shows the working scheme for this analysis.

The parameters optimization was done using a minimization algorithm, Nelder-Mead's simplex [45], where the minimization function was the Root Mean of the Squared Error (RMSE), described in Equation (9).

$$RMSE = \sqrt{\frac{1}{h} \sum\_{t=1}^{h} (\hat{X}\_t - X\_t)^2} \tag{9}$$

where *h* is the number of samples to be predicted.

Once the model was fit and the parameters obtained, 24-hours ahead forecast during two weeks were performed, and compared against the real data. The forecast accuracy was measured by using the Mean Average Percentage Error (MAPE), described in Equation (10).

$$MAPE(\%) = 100 \frac{1}{h} \sum\_{t=1}^{h} \frac{|\hat{X}\_t - X\_t|}{|X\_t|} \tag{10}$$

The use of RMSE for the model optimisation is preferred [46]. On the contrary, the most commonly used indicator to compare forecast accuracy is MAPE [47], in special when comparing demand forecasts. The triple seasonal analysis was performed similarly, but it was needed a higher dataset to fit the model and obtain the parameters. Thus, 80 random samples from three full years and two weeks were sampled. The first set of three years was used to fit the model, whereas the two week set was used for validation purposes. In this case, only the *AMC*24,168,8766 model was used for the analysis.

**Figure 1.** Demand split in years and seasons for stability analysis.
