TEXT SIZE

search for



CrossRef (0)
Forecasting with a combined model of ETS and ARIMA
Communications for Statistical Applications and Methods 2024;31:143-154
Published online January 31, 2024
© 2024 Korean Statistical Society.

Jiu Oha, Byeongchan Seong1,a

aDepartment of Applied Statistics, Chung-Ang University, Korea
Correspondence to: 1 Department of Applied Statistics, Chung-Ang University, 84 Heukseok-ro, Dongjak-gu, Seoul 06974, Korea. E-mail: bcseong@cau.ac.kr
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2022R1F1A1074340).
Received September 26, 2023; Revised November 13, 2023; Accepted November 27, 2023.
 Abstract
This paper considers a combined model of exponential smoothing (ETS) and autoregressive integrated moving average (ARIMA) models that are commonly used to forecast time series data. The combined model is constructed through an innovational state space model based on the level variable instead of the differenced variable, and the identifiability of the model is investigated. We consider the maximum likelihood estimation for the model parameters and suggest the model selection steps. The forecasting performance of the model is evaluated by two real time series data. We consider the three competing models; ETS, ARIMA and the trigonometric Box-Cox autoregressive and moving average trend seasonal (TBATS) models, and compare and evaluate their root mean squared errors and mean absolute percentage errors for accuracy. The results show that the combined model outperforms the competing models.
Keywords : ETS, ARIMA, hybrid models, state space models, forecasting performance
1. Introduction

The most prevalent methods for analyzing and forecasting time series data are exponential smoothing (ETS) and autoregressive integrated moving average (ARIMA) models. Since the ETS model was combined with the state space model (SSM), its development has begun to accelerate further. For details of the existing ETS models, see Gardner (1985, 2006) and Hyndman et al. (2008). The Box-Jenkins’ ARIMA model has become the most standard time series analysis model since Box and Jenkins (1976).

Recently, a lot of efforts has been made to upgrade the forecasting performances of the ARIMA and ETS models. For example, there are extensions to vector forms or the introduction of time-varying coefficients such as Seong (2020) and Campagnoli et al. (2009). The most successful effort may be the ensemble or combined models that can be classified into three types.

The first type is to independently fit each individual model and combine them only in the forecast step by using a combination technique such as a mean or median. Among others, Lee and Seong (2022) consider univariate time series models, which are well known in the field of forecasting, to study the forecasting performances of their simple combinations.

The second type is to complement some models with others, what is called the sequential approach. For example, we can consider the ETS model with ARIMA residuals (errors). Taylor (2003) showed that the forecasting performance is improved by fitting the AR(1) model to the residuals in the ETS model. As an extended form, De Livera et al. (2011) developed a Trigonometric Box–Cox ARMA trend seasonal (TBATS) model which consists of ETS components, a Fourier term and ARMA errors based on the Box-Cox transformed time series.

The last type considers a sort of model combination from the model building stage. A representative tool that makes it possible is the SSM. Further, the SSM with a single source of error (SSOE) is more convenient to use than with multiple sources of error (MSOE). See Hyndman et al. (2002, 2008) for details. In this context, we investigate, the ETS+ARIMA model, which is based on the SSM with the SSOE, that combines the two models and estimates them simultaneously. We explicitly describe the model setup and compare its forecasting performance with that of the existing methods by using two real data sets.

The structure of this paper is as follows. In Section 2, we introduce the innovations SSM with SSOE and represent ETS and ARIMA models with the SSM. In Section 3, we define the ETS+ARIMA model based on the SSM. In Section 4, we compare the forecasting performance of the ETS+ARIMA model with that of the existing forecast models. Section 5 presents the conclusion.

2. ETS and ARIMA by the SSM

The SSM consists of measurement and transition equations as follows:

yt=w(vt-l)+r(vt-l)ɛt,vt=f(vt-l)+g(vt-l)ɛt,

where vt is the state vector, containing the time series components, such as level, trend and seasonal; l is a lags vector; w(·) is a measurement function; r(·) is an error function; f (·) is a transition function; g(·) is a persistence function; and the error term εt follows the normal distribution with mean zero and variance σ2. Now we represent ETS and ARIMA models in the SSM form of Equation (2.1).

2.1. ETS model

The ETS model analyzes time series data using the three time series components: Error, trend, and seasonal. The fundamental models can be classified into the 15 forms shown in Table 1, depending on the combination of trend and seasonal components. As the error component has two types of additive or multiplicative errors, we can consider a total of 30 models.

In Table 1, the alphabetic symbols, A, Ad, M, and Md, indicate the additive, additive-damped, multiplicative, and multiplicative-damped trend components, respectively. N denotes the case where the trend or seasonal component is absent. Usually, we label the ETS model as ETS(*,*,*), which includes the alphabetic symbols between the parentheses that indicate error, trend, and seasonal components in order. For example, ETS (A,A,A) denotes that all three components are additively linked and its model equation for time series yt is expressed as follows:

yt=t-1+bt-1+st-m+ɛt,t=t-1+bt-1+αɛt,bt=bt-1+βɛt,st=st-m+γɛt,

where ℓt, bt, and st denote the level, trend, and seasonal components, respectively; m is the seasonal period; and α, β, and γ are the smoothing constants. φbt may be used instead of bt to dampen the trend where φ is a damping parameter. These equations can now be expressed in the SSM form of equation (2.1):

yt=wvt-l+ɛt,vt=Fvt-l+gɛt,

where

w=(111),         F=(110010001),         g=(αβγ),         vt=(tbtst),         vt-l=(t-1bt-1st-m),         l=(11m).

ETS(M,M,M), the pure multiplicative versions of ETS(A,A,A), can be formulated using natural logarithms in the following way:

log yt=wlog(vt-l)+log(1+ɛt),log vt=Flog(vt-l)+log(lk+gɛt),

where 1k = (1, 1, 1)′ is the vector of ones, containing k elements (number of components in the model). For more details, see Hyndman et al. (2008) and Svetunkov (2022a).

2.2. ARIMA model

The SSM representation for the ARIMA model is constructed through the level variable yt but not the differenced variable (1 − B)yt for a lag operator B with Bkyt = ytk.

First of all, the variable yt for ARIMA(p, d, q)(P, D, Q)m is expressed as follows:

yt=j=1Kηjyt-j+j=1Kθjɛt-j+ɛt,

where ηj and θj are coefficients for AR and MA polynomials, respectively. K is the highest polynomial order calculated by K = max(p + d + (P + D)m, q + Qm). For example, if the MA order is greater than the AR order in the maximum arguments, ηj = 0 when j > p + d + (P + D)m. In the opposite situation, θj = 0 when j > q + Qm. Then, we obtain the measurement equation:

yt=j=1Kvj,t-j+ɛt,vi,t-i=ηiyt-i+θiɛt-i.

Replacing yti of (2.8) by yt-i=Σj=1Kvj,t-j-i+ɛt-i, we can obtain the transition equation for i = 1, 2, …, K as:

vi,t-1=ηij=1Kvj,t-j-i+(ηi+θi)ɛt-i.

For example, ARIMA(1,1,2) for a nonseasonal time series yt with P = D = Q = 0 can be expressed in the SSM form:

yt=v1,t-1+v2,t-2+ɛt,v1,t=(1+η1)(v1,t-1+v2,t-2)+(1+η1+θ1)ɛt,v2,t=-η1(v1,t-1+v2,t-2)+(-η1+θ2)ɛt.

The matrix form, F, w, g, vt, and the lags vector l can be defined as follows:

F=(1+ηt-ηt),         w=(11),         g=(1+η1+θ1-η1+θ2),         vt=(v1,tv2,t),         l=(12).

Similar to the ETS model, the SSM representation for a multiplicative ARIMA model can be expressed by using natural logarithms as follows:

log yt=j=1Klog vj,t-j+log(1+ɛt),log vi,t=ηij=1Klog vj,t-j+(ηi+θi)log(1+ɛt)         fori=1,2,,K.

This multiplicative model is called the log ARIMA model.

3. The ETS+ARIMA model

3.1. The SSM representation

We now formulate the ETS+ARIMA model in the SSM. When the error is additive, the model can be expressed as follows:

yt=wEvE,t-lE+wAvA,t-lA+ɛt,vE,t=FEvE,t-lE+gEɛt,vA,t=FAvA,t-lA+gAɛt,

where the subscripts ‘E’ and ‘A’ correspond to the ETS or ARIMA models, respectively. This model simultaneously estimates the two model equations in contrast to the existing sequential approach such as the ETS model with ARIMA errors.

We examine the ETS(A,N,A)+ARIMA(2,0,0) model as an example of an ETS+ARIMA model.

The model equation can be expressed as follows:

yt=t-1+st-m+v1,t-1+v2,t-2+ɛt,t=t-1+αɛt,st=st-m+γɛt,v1,t=φ1v1,t-1+φ1v2,t-2+φ1ɛt,v2,t=φ2v1,t-1+φ2v2,t-2+φ2ɛt.

It is rewritten in the matrix form of Equation (2.3) where

w=(1111),         F=(1000010000φ1φ100φ2φ2),         g=(αγφ1φ2),         vt=(tstv1,tv2,t),         vt-l=(t-1st-mv1,t-1v2,t-2),         l=(1m12).

If the error is multiplicative, the model is expressed as follows:

logyt=wElog vE,t-lE+wAlog vA,t-lA+log(1+ɛt),logvE,t=FElog vE,t-lE+log(1k+gEɛt),logvA,t=FAlog vA,t-lA+gAlog(1+ɛt).

We note that these SSM representations make it easy to combine ETS and ARIMA. This is also true for the estimation of parameters, which will be discussed in Section 3.3. However, if we consider an ETS model with ARIMA errors, which, for example, has the following model equation, it will make model identification and estimation much more difficult.

yt=t-1+bt-1+st-m+ηt,t=t-1+bt-1+αɛt,bt=bt-1+βɛt,st=st-m+γɛt,

where ηt denotes ARIMA(p, d, q)(P, D, Q)m errors:

φ(B)ΦP(Bm)(1-B)d(1-Bm)Dηt=θ(B)ΘQ(Bm)ɛt,

and where φ(B) and ΦP(Bm) denote nonseasonal and seasonal AR polynomials, respectively, and θ(B) and ΘQ(Bm) denote nonseasonal and seasonal MA polynomials, respectively.

3.2. Identifiability of ETS+ARIMA

The ETS+ARIMA models should be treated with care in that the ETS and ARIMA models may duplicate each other in some cases. For example, when we look at the combination of ETS(A,N,N) and ARIMA(0,1,1), the model is defined as follows:

yt=t-1+v1,t-1+ɛt,t=t-1+αɛt,v1,t=v1,t-1+(1+θ1)ɛt.

However, the model does not have unique parameters and thus is unidentifiable because it is indistinguishable from a model where α* = αc and θ1*=θ1+c for arbitrary constant c are substituted for α and θ1, respectively.

As a general recommendation, Svetunkov (2022a) proposes a short list of guidelines to avoid unreasonable models that cannot be identified.

  • Use ARIMA(0, 1, q) if q > 1 or use ETS(A, N, N) if q ≤ 1 in modeling ETS(A, N, N)+ARIMA(0, 1, q).

  • Use ARIMA(0, 2, q) if q > 2 or use ETS(A, A, N) if q ≤ 2 in modeling ETS(A, A, N)+ARIMA(0, 2, q).

  • Use ARIMA(p, 1, q) when either p > 1 or q > 2, or use ETS(A, Ad, N) when p ≤ 1 and q ≤ 2 in modeling ETS(A, Ad, N)+ARIMA(p, 1, q).

  • Avoid using ETS(M, N, N)+log ARIMA(0, 1, 1).

  • Avoid using ETS(M, M, N)+log ARIMA(0, 2, 2).

  • Avoid using ETS(M, Md, N)+log ARIMA(1, 1, 2).

The three models mentioned in points (4) through (6) above are referred to as potentially unidentifiable combinations. It is also known that mixing additive ETS with multiplicative ARIMA or multiplicative ETS with additive ARIMA does not make sense from the modeling point of view.

To increase the identifiability of the model, the roles of ETS and ARIMA can be separated. In other words, a model (i.e., ETS with ARMA errors) could be constructed so that the ETS model part explains the non-stationarity and the ARMA part takes care of the stationarity of the remainder. For example, De Livera et al. (2011)’s TBATS model utilizes this concept. The ETS+ARIMA model considered in this study does not separate the roles for non-stationarity and stationarity. This is because the ETS and ARIMA models treat non-stationarity differently. Therefore, it would be better for ETS and ARIMA to complement each other in their roles for non-stationarity, especially if the overlap mentioned in Section 3.2 can be avoided. In other words, it is reasonable to view the meaning of “+” in the ETS+ARIMA model as an ensemble of the two models rather than the separation of roles.

3.3. Estimation and model selection

We can consider maximum likelihood estimation (MLE) for parameters depending on the distributional assumption. When the error is additive with εt ~ N(0, σ2), the log likelihood for y1, …, yT is expressed by

log=-T2log (2πσ2)-12t=1Tɛt2σ2.

Then, σ^2=(1/T)Σt=1Tet2 for residuals et. However, when the error is multiplicative, the log likelihood is expressed by

log=-T2log (2πσ2)-12t=1Tɛt2σ2-t=1Tlog|μy,t|,

where yt = μy,t(1 + εt).

The basic idea underlying ETS+ARIMA components selection is to regard ARIMA as an addition to ETS in the framework. This means that we can do model selection in the following steps:

  • Step (1): Select ETS components by using y1, …, yT,

  • Step (2): Extract Akaike information criterion (AIC) and residuals of the best ETS model in Step (1),

  • Step (3): Calculate the autocorrelation function (ACF) and the partial autocorrelation function (PACF) based on the residuals of Step (2), and create plausible candidate ARIMA models by using the ACF and PACF. For automatic modelling, we can use the Hyndman-Khandakar algorithm in Hyndman et al. (2002).

  • Step (4): Fit one suitable model among the candidates in Step (3) to the residuals in Step (2),

  • Step (5): If AIC of the current ETS+ARIMA model is lower than that in Step (2), go to Step (4). Otherwise, go to Step (6),

  • Step (6): As the final model, estimate the model with the lowest AIC among the models fitted from Steps (4) to (5).

Since the above model selection steps are to fit ETS first, and then consider ARIMA for the residuals, the initial values of the ETS components (level, trend and seasonal) required in Step (1) can be obtained in the same way as in the pure ETS model. For details, see Hyndman et al. (2008). Therefore, Steps (4) to (5) can be understood as a procedure of iteratively finding the order of the appropriate ARIMA model. Still, there may be various ways to construct an ETS+ARIMA model using with different sequences between ETS and ARIMA selections. Svetunkov (2022a) suggests that starting with ETS and then selecting ARIMA orders produce a robust forecasting model.

4. Performance evaluation for ETS+ARIMA

We fit the ETS+ARIMA models to two real data sets to evaluate its forecasting performance. The first data is the monthly particulate matter concentration (PM10) series of Jung-gu, Seoul obtained from the national statistical portal (https://kosis.kr/index/index.do). The data period is from April 2001 to June 2022 with 267 observations, and the training data and test data are divided on the basis of before and after June 2022. The second data is the annual Canadian lynx capture series from 1821 to 1934, with 114 observations. The test data is from 1925. We note that this series is non-seasonal time series but has a cycle of about 9 or 10 years. This data is available from R package ‘forecast’.

4.1. Monthly PM10

The time series plot for PM10 is shown in Figure 1 where there was no trend of increasing or decreasing to a certain level, except for the rise from October 2001 to June 2002. In the figure, we observe that there is a distinct seasonal period with m = 12 and the PM10 concentrations are high in every March and April.

4.2. Canadian annual lynx capture

Figure 2 shows the time series plot for the annual lynx capture in Canada. From the ranges 1820–1830, 1860–1870, and 1900–1910, the peak (highest value) of the capture was approximately 7,000; however, there was no tendency to increase or decrease to a specific level. It is interesting that, instead of seasonal period, there is a cycle repeating approximately every 9 to 10 years.

4.3. Evaluation of forecasting performance

In this section, we compare the forecasting performances of ETS+ARIMA models with that of the existing models: ETS, ARIMA, and TBATS. ETS+ARIMA model is fitted by using the adam() function in the smooth package (Svetunkov, 2017; 2022b) of R, ETS by ets(), ARIMA by auto.arima(), and TBATS by tbats() functions in the forecast package of R (Hyndman et al., 2023). To compute forecasting accuracies, we use two common measures. The first is the mean absolute percentage error (MAPE) shown below:

MAPE=1hj=1h|yT+j-y^T+j||yT+j|,

where T is the forecast origin and h is the forecast horizon. MAPE is the most widely used measure in evaluating the forecasting accuracy. It is an efficient measure when all the time series data are far from zero, but otherwise, it is likely to distort the overall error rate. Thus, the second measure is the root-mean-square error (RMSE) shown below:

RMSE=j=1h(yT+j-y^T+j)2h.

Using RMSE makes interpretation easier by converting error indicators into units similar to the actual values. However, it is sensitive to outliers because of its scale-dependency.

4.3.1. Performances for PM10 series

To estimate the order of the ETS+ARIMA model from the PM10 series in Jung-gu, Seoul, ETS(A,N,M) is suitable for the ETS component and ARIMA(0,0,16) or ARIMA(20,0,0) is suitable through the PACF and ACF of the residuals. We fitted two ETS+ARIMA models to examine the robustness of the model. The competing models, ETS and ARIMA, are fitted into ETS(A,N,M) and ARIMA(0, 1, 1)(0, 1, 1)12, respectively. The TBATS model is fitted with MA(1) errors and two Fourier-like seasonal terms for the monthly period without a Box-Cox transformation.

Figure 3 shows the true values (bold solid line) and 12-month-ahead forecasts of the three competing models and two ETS+ARIMA models. In Table 2, we compute two accuracy measures for the forecasting performance. The forecasting superiority may be indistinguishable through the figure, but when comparing RMSE and MAPE from the table, the superiority of ETS+ARIMA is clearly visible.

It seems that all kinds of forecasts overall do not fit the true values well; especially from mid-2020 to mid-2021. The forecast (two-dashed line) by the ARIMA model is relatively under-estimated. On the other hand, four models (two ETS+ARIMAs, ETS and TBATS) with elements of ETS show forecasts with very similar values. In conclusion, looking at the RMSE and MAPE values in Table 2, it is clear that two ETS+ARIMAs show the best forecasting performance.

4.3.2. Performances for lynx series

To estimate the order of the ETS+ARIMA model in Canada’s annual lynx capture series, ETS(M,N,N) is fitted for the ETS component and ARIMA(8,0,0) or ARIMA(3,0,0) is fitted through the residuals. The competing models used were ETS(M, N, N), ARIMA(8,0,0) and TBATS without seasonality (i.e., BATS model). Note that all three competing models were fitted as non-seasonal models because the lynx series has no seasonality.

Figure 4 shows the true values and 10-year-ahead forecasts of the three competing models and two ETS+ARIMA models. In Table 3, we compute two accuracy measures for forecasting performance.

The forecasting superiority is distinguishable from both the figure and the table where the superiority of ETS+ARIMA is clear, even if the ARIMA model shows the second-best performance. Interestingly, the ETS and TBATS models do not follow the characteristics of the cycle shown by the true values and only show horizontal forecasts. The ETS+ARIMA models fit the cycle of the time series well. This is a result that was not observed in the seasonal time series, such as PM10 data. In other words, it can be seen as evidence showing the potential of the ETS+ARIMA model’s excellent forecasting performance by specializing in stationary time series with cycles.

5. Conclusion

In this paper, we introduce the ETS+ARIMA model, which combines the two models simultaneously but not combining them sequentially. In addition, we evaluate the forecasting performance of the model by applying it to real data sets: Seasonal and cyclic time series. The results of the evaluation show the ETS+ARIMA model having better forecasting performances than the competing models. In particular, the model shows excellent forecasting performance in stationary time series with cycles. For further studies, it is necessary to investigate forecasting performances for purely non-seasonal or multivariate time series.

Figures
Fig. 1. Time series plot for monthly PM10 in Jung-gu, Seoul.
Fig. 2. Time series plot for Canadian annual lynx capture.
Fig. 3. Forecasts plot for PM10.
Fig. 4. Forecasts plot for lynx data.
TABLES

Table 1

The main exponential smoothing models

Trend componentSeasonal component

N (None)A (Additive)M (Multiplicative)
N (None)N, NN, AN, M
A (Additive)A, NA, AA, M
Ad (Additive damped)Ad, NAd, AAd, M
M (Multiplicative)M, NM, AM, M
Md (Multiplicative damped)Md, NMd, AMd, M

Table 2

Comparison of forecasting accuracies for PM10

ModelsRMSEMAPE
ETS(A,N,M)12.96723.803
ARIMA(0,1,1)(0,1,1)1213.52424.208
TBATS12.89624.219
ETS(A,N,M)+ARIMA(0,0,16)12.71123.670
ETS(A,N,M)+ARIMA(20,0,0)12.45423.944

Table 3

Comparison of forecasting accuracies for lynx series

ModelsRMSEMAPE
ETS1290.011136.880
ARIMA645.58235.045
TBATS1433.702154.509
ETS(M,N,N)+ARIMA(8,0,0)545.62425.600
ETS(M,N,N)+ARIMA(3,0,0)945.22841.946

References
  1. Box GEP and Jenkins GM (1976). Time Series Analysis: Forecasting and Control, San Francisco, Holden-Day.
  2. Campagnoli P, Petrone S, and Petris G (2009). Dynamic Linear Models with R, New York, New York, Springer.
    CrossRef
  3. De Livera AM, Hyndman RJ, and Snyder RD (2011). Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association, 106, 1513-1527.
    CrossRef
  4. Gardner ES (1985). Exponential smoothing: The state of the art. Journal of Forecasting, 4, 1-28.
    CrossRef
  5. Gardner ES (2006). Exponential smoothing: The state of the art-Part II. Journal of Forecasting, 22, 637-666.
    CrossRef
  6. Hyndman RJ, Athanasopoulos G, and Bergmeir C, et al. (2023) Forecast: Forecasting Functions for Time Series and Linear Models R package version 8.21.1 .
  7. Hyndman RJ, Koehler AB, Ord JK, and Snyder RD (2008). Forecasting with Exponential Smoothing: The State Space Approach, Berlin, Heidelberg, Springer.
    CrossRef
  8. Hyndman RJ, Koehler AB, Snyder RD, and Grose S (2002). A state space framework for automatic forecasting using exponential smoothing methods. International Journal of Forecasting, 18, 439-454.
    CrossRef
  9. Lee S and Seong B (2022). Performance for simple combinations of univariate forecasting models. The Korean Journal of Applied Statistics, 35, 385-393.
  10. Seong B (2020). Smoothing and forecasting mixed-frequency time series with vector exponential smoothing models. Economic Modelling, 91, 463-468.
    CrossRef
  11. Svetunkov I (Array). Statistical models underlying functions of ‘smooth’ package for R. Working Paper of Department of Management Science, Lancaster University 2017:1, 1-52.
  12. Svetunkov I (2022b) Smooth: Forecasting Using State Space Models, R package version 3.2.0 .
  13. Svetunkov I (2023). Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM), Florida, Chapman and Hall/CRC.
    CrossRef
  14. Taylor JW (2003). Short-term electricity demand forecasting using double seasonal exponential smoothing. Journal of the Operational Research Society, 54, 799-805.
    CrossRef