Application of time series models for heating degree day forecasting

In this study, the monthly HDD data in France were obtained from the official website of the Statistical Office of the European Union (2019). A total of 528 data covering the period of January 1974 and December 2017 were used to develop the forecasting model. The unit of data is°C * day. Table 1 presents a part of the HDD data used in the analysis.

A time series consists of four different components, namely trend, seasonal, cyclical, and random components. The individual component analysis of the time series is the basis for the decomposition method. There are two decomposition approaches, which are additive and multiplicative decomposition. In the analysis of consecutive data, the multiplicative decomposition approach is preferred. Therefore, in this study, this approach was used, in which time series can be expressed as the product of the four components of the series and are Y = T.S.C.I, where Y represents the original data, whereas T, S, C, and I represent the trend, seasonal, cyclical, and random components of the series, respectively (Calis et al. 2017). In this method, first, the effect of the random change is reduced by applying an exponential correction to the data. Then, the seasonal and trend components of the data are eliminated by calculating the seasonal index and trend equation, respectively. Finally, the model was created assuming that the seasonal and trend components are in multiplication.

The Box–Jenkins method, also known as the autoregressive integrated moving average (ARIMA) model, was used for time-based analysis and time-based modeling of HDD data. This method applies the autoregressive moving average (ARMA) or ARIMA models to find the best fit of a time series model to the historical values of a time series (Baldigara and Koic 2015; Box et al. 2015). The seasonal autoregressive integrated moving average model (SARIMA) is an expanded form of ARIMA. SARIMA processes are designed to model trends, seasonal patterns, and correlated time series and have proven to be successful in estimating short-term fluctuations (Kam et al. 2010; Baldigara and Koic 2015). The SARIMA model consists of (1) automatic regression, (2) difference, and (3) moving average. SARIMA model is represented as SARIMA (p, d, q) (P, D, Q)_S in which p, d, q represent the degree of nonseasonal linear autoregressive model (AR), non-seasonal difference, and degree of nonseasonal moving average (MA) model, respectively, whereas P, D, Q represent the degree of seasonal AR, seasonal difference, and the degree of seasonal MA. In addition, the seasonality length is represented by S. The d and D parameters are considered when the series is not stationary. In this study, to select the most appropriate SARIMA (p, d, q) (P, D, Q)_S model for HDD series, the significance of the coefficients of the models was checked by the Ljung-Box-Pierce Chi-squared statistics and t-test. In addition, the corrected R² values, the sum of the squares of the residuals, the Akaike Information Criteria (AIC) and the Schwarz Information Criteria (SIC) were considered. Furthermore, the residuals analysis including the autocorrelation function (ACF) and the partial autocorrelation function (PACF) of the residuals were conducted. Minitab 18.0 and EViews 10.0 packet programs were used to calculate all the statistics and the model selection criteria.

To evaluate the performances of the models developed by the decomposition and Box–Jenkins methods, the accuracy of the fitted values was determined by calculating mean absolute percentage error (MAPE), mean absolute deviation (MAD), and mean squared deviation (MSD), which are formulated as follows: MAPE=∑t=1n|(yt−y^t)/yt|n×100 (yt≠0)MAD=∑t=1n|(yt−y^t)|nMSD=∑t=1n|(yt−y^t)|2n\begin{array}{*{20}{c}} {{MAPE} = \frac{{\sum\nolimits_{t = 1}^n {\left| {\left( {{y_t} - {{\hat y}_t}} \right)/{y_t}} \right|} }}{n} \times 100\,\,\,\,\,\,\,\,\,\,\,\,\left( {{y_t} \ne 0} \right)} \\ {{MAD} = \frac{{\sum\nolimits_{t = 1}^n {\left| {\left( {{y_t} - {{\hat y}_t}} \right)} \right|} }}{n}} \\ {{MSD} = \frac{{\sum\nolimits_{t = 1}^n {{{\left| {\left( {{y_t} - {{\hat y}_t}} \right)} \right|}^2}} }}{n}} \end{array} MAPE corresponds to the relative scale of the prediction error between the forecasted value and the actual value. It should be noted that the smaller the error is in MAPE, the more accurate the prediction is (Kam et al. 2010). Similarly, the smaller the MAD and MSD are, the more accurate the prediction is.

In this study, the ratio-to-moving-average method was used for calculating the seasonal index. In this method, first, the moving averages were calculated and the effects of seasonal and random changes were eliminated. Since the period numbers were equal to the season type numbers and the months of the year represented the seasons in the HDD data set, 12-months moving averages were calculated. Then, using the calculated moving averages, the central moving averages and the ratios of the actual data to these averages were calculated. Most of the random changes were eliminated by calculating the arithmetic mean of these ratios.

It should be noted that the sum of the calculated averages of each month should be 12.00. However, the fact that the sum of the averages was obtained as 11.99, the correction factor (1.0008) was calculated and applied to the calculated averages of the months. The corrected averages, in other words, seasonal indices, were obtained by multiplying the average ratio value of each month by this correction factor (Table 2). It can be said that the number of HDD in January is 109% higher than the typical index which is 100%, whereas in July and August, it is 94% less.

Next, the actual HDD data were divided into the seasonal indices to obtain the seasonally adjusted data. Actual data and the seasonally adjusted data are shown in Figure 3.

Fig. 3

After the seasonal effect was eliminated from the series, the trend effect of the data was determined and eliminated to develop the forecasting model. For this purpose, the trend equation was determined for the seasonally adjusted data. It should be noted that the time data are numbered as 1, 2, 3..., 528 for the trend equation. Accordingly, the trend line and trend equation were obtained, as shown in Figure 4.

Fig. 4

The trend values were calculated by using the trend equation. The trend component was removed from the data set by dividing the seasonally adjusted data by the trend values. The graph of trend corrected data is presented in Figure 5.

Fig. 5

After eliminating the trend effect, the random component effect in each time series was removed by the exponential correction method. In this method, the corrected time series data at time t was calculated by substituting the numbers from 0.1 to 0.9 with intervals of 0.1 instead of α in the equation below. Yt=αyt+α(1−α)yt−1+α(1−α)2yt−2+…+(1−α)t−1y1(0≤α≤1;t=1,2,3,…n)\begin{array}{*{20}{c}} {{Y_t} = \alpha {y_t} + \alpha \left( {1 - \alpha } \right){y_{t - 1}} + \alpha {{\left( {1 - \alpha } \right)}^2}{y_{t - 2}} + \ldots + {{\left( {1 - \alpha } \right)}^{t - 1}}{y_1}} \\ {\left( {0 \leqslant \alpha \leqslant 1;t = 1,2,3, \ldots n} \right)} \end{array} In addition, the sum of the squares of the differences between the actual data and the corrected data was calculated for each α value. The least squares sum was obtained for α = 0.9. Therefore, this α value was chosen for forecasting future observations. As a result of the first-order exponential correction by taking α = 0.9, the random component is removed and the graph of the obtained cyclic indices is given in Figure 6.

Fig. 6

As a result, since the random component is eliminated, the forecast equation is obtained as follows: Yc=Tc*S*T{{\text{Y}}_{\text{c}}} = {{\text{T}}_{\text{c}}}*{\text{S}}*{\text{T}} where Y_c represents the forecasted HDD values; T_c represents the trend and the seasonality effect eliminated HDD values; S is the seasonal index and T is the trend index. The actual HDD values and the forecasted values calculated by using the equation are presented in Figure 7.

Fig. 7

In this study, 79 tentative models were generated. Table 3 presents 44 of these models. It should be noted that 44 tentative models shown in the table were generated without constant. Tentative models other than those with a star are also tested by including constant.

As a first step, the models which comply with both of the following conditions were selected: (1) the models with correlations that are statistically significant; (2) the models with the Chi-squared statistics of the residuals of k=12, 24, 36 months delay are statistically insignificant, in other words, the residuals of the models are not correlated with each other. As a result, nine models were selected for further analysis. In the next step, the adjusted R² value, the residual sum of squares, the AIC, and the SIC were calculated to select the appropriate model. The selected models and the selection criteria values are presented in Table 4.

Regarding the selection of the most appropriate forecasting model among the models selected in the first step, the criteria shown in Table 5 were considered. The higher the model’s adjusted R² value, the lower the residual sum of squares, the lower the AIC and the SIC are, the better the model is. Accordingly, the model that addresses most of the criteria at the same time is selected. Table 3 shows that the adjusted R² value of the models varies between 0.464422 and 0.936911 and the highest values are observed in the first and the seventh models. Then, the residuals sum of squares and AIC and SIC values of these models were compared. The values of the aforementioned criteria for the seventh model are smaller than those of the 1^st model. As a result, the seventh model, namely SARIMA (2,0,1)(1,0,1)₁₂, was selected as the suitable and the final forecasting model. However, it should be noted that the graph of the residuals should be examined to check the appropriateness of the final forecasting model. The graphs obtained through Minitab are presented in Figure 8.

Fig. 8

It should be noted that the residuals should be normally distributed with zero mean and constant variance. The normal probability plot and histogram in Figure 5 show that residuals except for the last residuals approximately suit to normal distribution. In addition, “versus fits” graph shows that the variance tends to increase slightly. It should be noted that this situation can be caused by the large variation of the data in the series. In addition to the plots in Figure 3, ACF and PACF plots of the residuals, which are the most examined plots in the model selection, are presented in Figures 9 and 10.

Fig. 9

Fig. 10

The plots indicate that most of the autocorrelations of the residuals are zero, in other words, they are not correlated. Thus, it can be concluded that the chosen model is appropriate.

The monthly actual HDD values, as well as the monthly forecasted HDD values and the confidence intervals obtained through the SARIMA (2,0,1) (1,0,1)12 model, are shown in Figure 11. The plot shows that the forecasted HDD values are close to the actual HDD values. In addition, the actual and forecasted HDD values of 2017 are presented in Table 4, which also indicates that the forecasted HDD values are close to the actual HDD values.

Fig. 11