A hybrid of artificial neural network, exponential smoothing, and ARIMA models for COVID-19 time series forecasting

2.1Description of the datasets

There are currently 85,040,066 confirmed cases and 1,844,920 deaths from the coronavirus COVID-19 outbreak as of January 03, 2021, 10:37 GMT.

We consider total confirmed cases as of close of business on January 22, 2020 through June 19, 2020. In addition, we considered the data for the outcomes of all closed cases, which resulted in death, from February 2, 2020 through June 19, 2020 (Worldometer, 2020).

Figure 1 shows the total confirmed COVID-19 cases and reflects the fact that there were 1,201,483 confirmed cases of COVID-19 worldwide at the end of the day on June 19, 2020. This figure includes deaths and recovered or discharged patients (i.e., cases with an outcome). Total confirmed cases ranged between 580 and 8,750,990 with a mean of 2,418,848.85. This summary indicates that the spread of total confirmed cases varies among the various countries.

Figure 1.

Time series plot of total confirmed cases and death rate of COVID-19.

In addition, Fig. 1 shows the outcome of all closed cases (death rate) from February 2, 2020 through June 19, 2020. This figure illustrates the cumulative number of total deaths over the cumulative number of closed cases. The death rate ranges between 5.64% and 41.8% with a mean of 15.01 and a median of 14.32.

In this study, 10% of the sample size is used as the testing sample. A training sample is used for the model building, and the testing sample is used for the model validation at the end of analysis. We considered the first 135 observations as training sample for the total confirmed cases and 125 observations as training sample for the death rate of COVID-19 data.

2.2Assessment of forecasting accuracy

The accuracy of forecasts can be determined by considering how well a model performs on new data that were not used when fitting the model. Accuracy is an important issue in forecasting; therefore, researchers tend to add more and more variables to their proposed model. Safi and White (2017) considered the issue of whether a complex model does a better job than a simple one.

Many measures of forecasting accuracy have been developed and several authors have discussed the fundamental usage of these measurements and compared the accuracy of forecasting methods with univariate time series data (see, i.e., Cryer & Chan, 2008; Hyndman & Athanasopoulos, 2018; Wei, 2006). The best forecasting models selected will be compared using three different forecasting accuracy measuring criteria: mean absolute error (MAE), the root mean squared error (RMSE), and mean absolute percentage error (MAPE).

Definition 1. Suppose y1,y2,…,yn denotes the dataset. The MAE is then defined as:

Definition 2. The RMSE is defined as:

The RMSE has been popular, largely because of its theoretical relevance in statistical modeling (Hyndman et al., 2002). However, this measure is more sensitive to outliers than the MAE, which has led some authors (e.g., see Armstrong, 2001) to recommend using other measures of forecasting accuracy. The MAE is preferable in cases of the existence of outliers. Using the MAE or RMSE is recommended when comparing forecasting methods on a single data set. This means, the MAE and RMSE should be used if all forecasts are measured on the same scale.

Definition 3. The MAPE is defined as:

The MAPE presents the forecast error in terms of percentage and hence it is scale invariant and unit-free (Lyhagen et al., 2015). The MAPE is a simple average of absolute percentage errors. It is recommended to use MAPE when comparing the accuracy of the same or different methods on different time series data with different scales, unless the data contain zeros or small values (Hyndman & Koehler, 2006).

The evaluation criterion for these measures of forecasting accuracy is that the smaller the value obtained, the better the model’s forecasting ability (McKenzie, 2011).

Definition 4. The efficiency of the proposed forecasting method relative to that of the benchmark method in terms of the RMSE, ψ, is defined as:

where RMSEp and RMSEb represent the RMSE from the proposed and the benchmark methods, respectively. Usually the benchmark method is the most naïve method (Hyndman, 2006). A ratio of less than one indicates that the forecasting performance of the proposed method is more efficient than the benchmark method and if ψ is close to one, then the proposed forecasting method is nearly as efficient as the benchmark forecasting method. Otherwise, the proposed method performed poorly (White & Safi, 2016).

3.1The ARIMA model

The general ARIMA (p,d,q) model is given by Box et al. (2015)

where d⩾1 is the degree of differencing, ∇=1-B is the differencing operator, and the lag operator B is defined as B⁢Yt=Yt-1, the operator that provides the previous value of the series. ϕ⁢(B) and θ⁢(B) are polynomials of degrees p and q in B,

and

The best ARIMA model is chosen according to its Akaike information criterion (AIC), AICc, or BIC value.

3.2The exponential smoothing model

The exponential smoothing method has been in use since the 1950s and is still the most popular forecasting method used in business and industry. The field of exponential smoothing has undergone a substantial evolution over the past 10 years.

Let Y1,Y2,…,Yn denote the time series of interest and let 𝒙t=(lt,bt,st,st-1,…,st-(m-1)), where lt denotes the level, bt denotes the trend, and st denotes the seasonal component, all at time t. Then, the state space models are:

where {εt} is the residual at time t and εt∼NID⁢(0,σ2). The methodology of the exponential smoothing state space model is fully automated using R-software (Hyndman et al., 2008).

3.3The ANN model

The nnetar function is used to fit neural networks. This function is described as a feed-forward neural network with a single hidden layer and lagged inputs for forecasting univariate time series. The nnetar function fits a neural network autoregressive 𝑁𝑁𝐴𝑅⁢(p,P,k) model. For a non-seasonal time series, the default is the optimal number of lags, according to the AIC value, for a linear autoregressive 𝐴𝑅⁢(p) model. For a seasonal time series, the default values is P=1 where p is chosen from the optimal linear model fitted to the seasonally adjusted data and k=12⁢(p+P+1) (rounded to the nearest integer). By default, 25 networks with random starting values are trained and their predictions are averaged (Hyndman, 2004).

3.4The hybrid model

The hybrid model fits multiple individual model specifications to enable the easy creation of ensemble forecasts. The hybrid model consists of a combination of three models: the ARIMA, the exponential smoothing, and the ANN models. Looking at a time series composed of autocorrelated linear and nonlinear components, we have:

Fitting L^t using the ARIMA model with et as the residual yields:

The error term consists of nonlinear relationships with previous errors. The nonlinear relationships can be modeled from the past residuals as follows:

Then, using an ANN model to predict et as an estimate for Nt, we can calculate the forecast:

These models were built using hybrid Model command in the forecast Hybrid package. This package fits multiple models and combines them using either equal weights or weights based on in-sample errors. There are six models: ARIMA, exponential smoothing, theta, neural network autoregression (NNAR), seasonal and trend decomposition, and the trigonometric seasonal + exponential smoothing method + Box-Cox Transformation (TBATS) model for heterogeneity and the autoregressive moving average model for residuals + trends + seasonal (including multiples and non-integer periods).

4.1Test for normality

We carried out the Anderson-Darling (AD) normality test to determine if the residuals after approximation for the four models followed a normal distribution (Thode, 2002).

For the total confirmed cases’ residuals of COVID-19 data, normality test yielded the following AD values: 8.8922, 7.5803, 6.22, and 9.645, with corresponding p-values of < 2.2e-16, < 2.2e-16, 2.333e-15, p-value < 2.2e-16 for ARIMA, ETS, ANN, and Hybrid models, respectively.

For the death rates residuals’, the normality test yielded the following AD values: 4.976, 0.535, 2.368, and 1.029, with corresponding p-values of 2.25e-12, 0.168, 5.002e-06, and 0.01006, for ARIMA, ETS, ANN, and Hybrid models, respectively. Both the test for normality and the probability plots for the residuals of the total confirmed cases and the death rate indicate that the normality assumption was not satisfied for all cases except for the residuals of the death rate for ETS model.

Since the data are not normally distributed (without loss of generality), to compare the performance of the models through the two datasets, we used the forecasting accuracy measure, the MAE, over the forecasting period for each model. The smaller values of MAE indicate higher forecasting accuracy. Therefore, the ratios of the MAE of the hybrid model to those of the ARIMA, ETS, and ANN models were calculated for analysis.

4.2Relative efficiencies of the empirical results

Table 1 lists the complete empirical results for the ratios of the hybrid model’s RMSE, MAE, and MAPE of residuals after approximation to those of the ARIMA, ETS, and ANN models for total confirmed cases and the death rate. Since the data are not normally distributed, we discuss these empirical results based on the ratios of the hybrid model’s MAE to those of the ARIMA, ETS, and ANN models.

4.3Total confirmed cases data

Applying the ANN model with average of 10,000 networks, each of which is a 1-25-1 network, with 76 weights and estimated noise variance of 1.33e+08. This indicates that 25 networks were trained and that their predictions were averaged. For the hybrid model, the weights for the ARIMA, ETS, and ANN models were chosen to be equal. For the ARIMA model, the result shows that the best fit model was the ARIMA (0, 2, 1) with estimated moving average coefficient -0.6791 and estimated noise variance of 1.77e+08 (with AIC = 2907.24, AICc = 2907.33, and BIC = 2913.02).

For the exponential smoothing model, simple exponential smoothing with additive errors is fitted. The estimated values of the model smoothing parameters are α^= 0.8779, β^= 0.3243, with initial level l0=-2273.75 and initial trend b0= 1507.71 and the estimated noise variance equals 13271.8 (with AIC = 3231.369, AICc = 3231.834, and BIC = 3245.895).

The relative efficiencies based on residuals after approximation of the hybrid model to the ARIMA, ETS, and ANN models equal ψ= 1.5089, 1.5638, and 0.5890, respectively. Equivalently, the relative efficiencies of the ETS model to the ARIMA, ANN, and Hybrid models equal ψ= 0.9649, 0.3766, and 0.6395, respectively. This result indicates that the ETS MAEs equal 96.49%, 37.66%, and 63.95% of that of the ARIMA, ANN, and Hybrid models, respectively.

Therefore, the ETS model perfroms more efficient as compared to ARIMA, and was superior than the ANN and Hybrid models for total confirmed cases of COVID-19 data. However, as a second choice in this case, we could use the ARIMA model since it is almost as efficient as the ETS model, keeping in mind that it is not a perfect substitute.

4.4Death rate data

The ANN model with an average of 10,000 networks, each of which was a 1-25-1 network, with 76 weights and estimated noise variance of 0.2494. This indicates that 25 networks were trained and that their predictions were averaged. For the hybrid model, the weights for the ARIMA, ETS, and ANN models were chosen to be equal. For the ARIMA model, the result shows that the best fit model is ARIMA (2,0,2) with zero mean and estimated noise variance equaling 0.09114 (with AIC = 66.23, AICc = 66.73, and BIC = 80.37). For the exponential smoothing model, simple exponential smoothing with multiplicative errors is fitted. The estimated values of the model smoothing parameters are, α^=β^= 0.711, ϕ^= 0.925. The initial level l0= 45.4072 and initial trend b0=-2.6433 with an estimated variance of 0.0156 (with AIC = 232.8629, AICc = 233.5748, and BIC = 249.8328).

Figure 2.

Forecasts and Actual values in total confirmed cases and death rate.

Figure 3.

Time Series Plot of percentage change* from the forecast in total confirmed cases and death rate. *The percentage change is (Xt-Xt-1)/Xt-1×100%, where Xt and Xt-1 represent the total confirmed cases or death rate at time t and t-1, respectively.

The relative efficiencies based on residuals after approximation of the hybrid to the ARIMA, ETS, and ANN models equal ψ= 1.1682, 1.1417, and 0.7593, respectively. Equivalently, the relative efficiencies of the ARIMA model to the ETS, ANN, and Hybrid models equal ψ= 0.9773, 0.3766, and 0.6395, respectively. This result indicates that the ETS MAEs equal 97.73%, 37.66%, and 63.95% of that of the ETS, ANN, and Hybrid models, respectively. Therefore, the ARIMA model perfroms more efficient as compared to ETS, and was superior than the ANN and Hybrid models for death rate of COVID-19 data. However, as a second choice in this case, we could use the ETS model since it is almost as efficient as the ARIMA model.

Therefore, the hybrid model performed poorly as compared to the ARIMA and ETS models and was more efficient than the ANN model for forecasting the total confirmed cases and the death rate of COVID-19 data. Therefore, the ANN model is not the best choice for forecasting COVID-19 death rates, and ARIMA outperforms the other models.

Figure 4.

Actual data in total confirmed cases and death rate.

Figure 5.

Time Series Plot of percentage change in total confirmed cases.