A Deep Learning BiLSTM Encoding-Decoding Model for COVID-19 Pandemic Spread Forecasting

Abstract

The COVID-19 pandemic has widely spread with an increasing infection rate through more than 200 countries. The governments of the world need to record the confirmed infectious, recovered, and death cases for the present state and predict the cases. In favor of future case prediction, governments can impose opening and closing procedures to save human lives by slowing down the pandemic progression spread. There are several forecasting models for pandemic time series based on statistical processing and machine learning algorithms. Deep learning has been proven as an excellent tool for time series forecasting problems. This paper proposes a deep learning time-series prediction model to forecast the confirmed, recovered, and death cases. Our proposed network is based on an encoding–decoding deep learning network. Moreover, we optimize the selection of our proposed network hyper-parameters. Our proposed forecasting model was applied in Saudi Arabia. Then, we applied the proposed model to other countries. Our study covers two categories of countries that have witnessed different spread waves this year. During our experiments, we compared our proposed model and the other time-series forecasting models, which totaled fifteen prediction models: three statistical models, three deep learning models, seven machine learning models, and one prophet model. Our proposed forecasting model accuracy was assessed using several statistical evaluation criteria. It achieved the lowest error values and achieved the highest R-squared value of 0.99. Our proposed model may help policymakers to improve the pandemic spread control, and our method can be generalized for other time series forecasting tasks.

1. Introduction

During the writing of this paper, an infectious disease from a new generation of Coronavirus (Coronavirus disease 2019, COVID-19) appeared in several countries. Transportation restrictions also underwent changes due to the COVID-19 pandemic spread.

As a precautionary procedure, several countries have stopped international air flights more than once. The first version of COVID-19 appeared in Wuhan city in China at the end of December 2019. COVID-19 was announced by the World Health Organization (WHO) to be a global pandemic on 11 March 2020 [1]. COVID-19 exponentially spread over all the world and highly affected the healthcare systems in several countries [2]. The total number of positive confirmed cases reached about 80 million people with death cases of about 1.7 million people [3]. This high death rate has led researchers and scientists in different fields to look for ways to address the challenges of this virus and work on overcoming the epidemic.

Artificial intelligence (AI) tool plays an important role to help the human survival during the COVID-19 pandemic [4,5]. Many researchers have employed AI tool in different applications [6,7,8,9,10,11,12].

In [13], an early detection and diagnosis system for COVID-19 was proposed. The patient provided his early disease symptoms, and the patient was classified as a positive or negative case. In [14], people contact tracing systems were applied in nine counties with several tracing technologies, such as Bluetooth, WiFi, and GPS. This tracing aimed to reduce the pandemic spread.

In [15], several AI algorithms were applied to CT and X-Ray imaging modalities to classify and detect COVID-19. These algorithms have succeeded in differentiating between COVID-19, Non-COVID-19, and pneumonia. These algorithms have achieved high accuracy classification reached to 99 %. In [16], an AI algorithm was introduced to recognize the potential old drugs with anti-coronavirus activities.

In [17], AI was integrated with surveillance systems that enable Non-mask violation detection. In [18], an AI algorithm was able to monitor and track social distancing. The system monitored the minimum safe distancing 2-m, which the governments and health authorities have recommended in shopping centers, streets, and schools. In [19], the AI algorithm was introduced to forecast the future spread of COVID-19 in several countries and globally. Therefore, forecasting pandemic spread represents a real important policymaker parameter that needs accurate prediction models.

Deep learning networks were introduced by LeCun [20], which represent the new era of AI algorithms. The traditional AI algorithms are based on pre-processing, feature extraction and selection, and classification. These sequences consume much time and lacking to generalization [21]. Deep learning introduced the convolutional neural network (CNN), which extracts the features without pre-processing, hand-crafted features extraction. CNN can automatically extract valuable features through consequences convolutional layers. These layers are defined as separable learnable parameters, which acquire their parameters during their training process.

The recurrent neural network (RNN) is more efficient with time series analysis problems. A special type of recurrent neural network (RNN) is long-short term memory (LSTM), which works in one-directional. Bidirectional long-short term memory (BiLSTM) is proposed to work in both past and future directions. Both LSTM and BiLSTM can extract temporal features and predict time-series [22]. Gated recurrent unit (GRU) network is also a special type of deep network with a simpler structure than LSTM, with higher performance and less training processing time.

Deep learning has been demonstrated in several applications for the COVID-19 pandemic, and it achieved better results [23]. There is a real need to employ deep network capabilities to build a high-performance pandemic spread prediction model. Policymakers are always in real need of accurately forecasting the pandemic spread. The correct forecasting guides the governments to make the correct decisions, such as (1) opening and closing procedures that stop the pandemic outbreak and save human life, (2) preparing the healthcare systems for patient flow demand, especially the intensive care units (ICU) demand and providing the medical needs, (3) applying the correct economic decision for rationalizing the government spending.

One of the most crucial COVID-19 research dimensions is forecasting the pandemic spread through confirmed, deaths, and recovered cases. Therefore, the COVID-19 pandemic spread prediction was investigated based on countries, period, and model. Each country has its own characteristics, such as population, contact rate, healthcare sector efficiency, strict opening and closing policies. Therefore, each country prediction spread model can lead to individual findings. At country level, the previous researches focused on some countries, such as Spain [24], Italy [24], China [25], India [26], Saudi Arabia [27], Brazil [28], South Africa [29], and United States of America (USA) [30].

With the start of the COVID-19 spread, scientists attempted to build models to predict the future scenario of the pandemic to estimate the peak of the confirmed cases and death cases. Most of these studies gathered the data with starting of February 2020 to the end of March 2020. Then, the scientists tried to find the answer to the pandemic spread stopping. Therefore, they included the period starting from early February to the middle of June.

There are five main frameworks for COVID-19 that were utilized to achieve the COVID-19 spread prediction task. The first framework is based on statistical models, which are auto-regressive (AR), moving average (MA), and autoregressive integrated moving average (ARIMA). The second framework is based on machine learning models. Some of these models are regression models, such as linear regression, Lasso regression, Elastic Net, Theil-Sen regression, RANSAC regression, Huber regression, support vector machine (SVM), and decision trees. The third framework is based on deep learning models, such as CNN, LSTM, BiLSTM, and GRU. The fourth framework is based on dynamic system models, which are susceptible infectious recovered (SIR), susceptible exposed infectious recovered (SEIR), immunity susceptible exposed infectious recovered (MSEIR). Finally, the fifth framework is based on recently developed prophet models, such as Prophet Model and Google trend prophet model.

The remainder of our paper is organized as follows. Section 2 discusses the literature for the time series prediction task of the COVID-19 pandemic. Section 3 discusses the material utilized in this study with our proposed method. Section 4 presents the experimental section with a detailed discussion. Finally, Section 5 concludes our work and suggests future works.

2. Review of Predictions Models

This section covers most of the time series prediction models utilized to forecast the COVID-19 pandemic spread. First, we cover the most well-known time series forecasting algorithms based on statistical, machine learning, and deep learning approaches. Second, we cover the pandemic spread in different countries during a different times. Finally, we cover the encoder-decoder deep learning models that have been applied for several time-series forecasting tasks.

Several studies in the literature covered different countries and presented models that can predict the COVID-19 pandemic spread. In [31], the authors proposed a comparative study to predict the confirmed positive cases in several countries. The study included the period from 22 January 2020 till 17 June 2020. They have compared between the CNN LSTM, stacked LSTM, and BiLSTM. They assessed their results based on several metrics. However, they have applied their model to forecast cases for 16 days. Their deep learning models hyper-parameters have not been tuned, and they did not compare their results to other time series forecasting statistical models.

In [24], the authors utilized the ARIMA model to predict the confirmed positive cases in the highest infected countries in Europe, which are Italy, Spain, and France. The study included the period from 21 February 2020 to 15 April 2020. Their study achieved RMSE value 1654 for Italy, RMSE value 2031 for Spain, RSME value 971 for France. However, the study did not cover other statistical models. In addition, the studying duration was limited. The achieved RMSE Value for the three countries’ pandemic spread prediction was very high.

In [32], the authors proposed an excellent comparative study between statistical models (ARIMA and SVR) and three deep models (LSTM, BiLSTM, and GRU). Their study included the period from 22 January 2020 to 27 June 2020. Their study covered 10 different countries. The BiLSTM network achieved the best performance with RMSE value 0.007 for China.

In [33], the authors have proposed a deep assessment methodology and fractional calculus to predict the pandemic spread in 8 different countries (China, Italy, Turkey, Spain, France, Germany, UK, and the USA). The study included the period from 13 February 2020 to 19 June 2020. They built a deep assessment algorithm based on fractional calculus and deep LSTM. However, the daily confirmed cases curve is estimated to be a Gaussian function, which is not suitable for every country as they found. In [34], the authors have proposed a machine learning method based on exponential smoothing (ES) algorithm to predict the pandemic spread in 8 different countries (China, Italy, Turkey, Spain, France, Germany, UK, and the USA).

The study included the period from 22 January 2020 to 4 June 2020. The proposed algorithm is compared to other machine learning algorithms, such as linear regression (LR), support vector machine (SVM), least absolute shrinkage, and selection operator (LASSO). Their study was applied on a global dataset, ignoring each country’s particular case study or comparing their findings with well-known statistical time series forecasting models. Their experiments showed that the ES performed its best performance when following it by LR and LASSO models. However, this strategy was more complex and required a high computational cost.

In [35], the authors proposed an SEIR model to predict COVID-19 confirmed cases. The study included the period from 22 February 2020 to 8 October 2020. Furthermore, the study included two countries, which are Egypt and Iraq. Their proposed model estimated that the peak value of the confirmed cases is 4254 cases in Iraq and 1534 cases in Egypt based on the Gaussian model, while their model estimated that the peak value is 490,900, and 105,000 in Iraq and Egypt, respectively based on the logistic model. However, their findings have not been compared with similar studies in the literature.

In [36], the authors proposed seven machine learning algorithms to predict the positive confirmed cases in Egypt. The study included the period from 15 February 2020, to 15 June 2020. They investigated the performance of different types of regression models, which are are exponential polynomial, quadratic, third-degree, fourth-degree, fifth-degree, sixth-degree, and logit growth. However, their results were good, and we cannot generalize their study.

In [37], the authors introduced a COVID-19 spread prediction model based on google trend algorithm in India, USA, and the UK. The authors applied the study in the period from 24 February 2020 to 20 May 2020. The study included predicting the daily new cases, the cumulative cases, and the death cases. They optimized their proposed system based on Grey Wolf Optimizer (GWO) and compared their findings with the ARIMA prediction model.

In [38], the authors proposed a simple statistical model to predict positive conformed cases globally. The authors applied the study in the period from 22 January 2020 to 21 May 2020. Their model has not been compared with other models. On the other hand, their forecasting results were negatively biased in the case of death cases predictions and positively biased when new spikes of death and confirmed cases happened. They mentioned that their forecasting model would perform better when the established wave remains stable.

In [39], the authors proposed a deep learning model to predict the confirmed cases in several countries, which are Brazil, India, Russia, South Africa, Mexico, Peru, Chile, Colombia, and Iran. The authors applied the study from 20 January 2020, to 3 August 2020. The deep learning models were based on multi-head attention and LSTM network. Then, the authors optimize their model hyper-parameters based on Bayesian optimization.

Each country represents a particular case according to several characteristics, such as population, contact rate, number of tourists, opening and closing policies. Saudi Arabia’s population has reached about 34 million citizens as reported by the unified national platform [40]. Furthermore, there are millions of foreign workers that are existed in Saudi Arabia. As its importance to the muslims, there are about 2 million that are visiting Saudi Arabia every year to perform their Hajj. On the other hand, millions of Muslim people are visiting Saudi Arabia continuously whole the year to perform Umrah. Therefore, there is a real need to predict the future confirmed cases of COVID-19 to help policymakers make recommendations and procedures that can stop the spread of infectious cases. There are a few articles in the literature that have covered Saudi Arabia’s pandemic spread prediction.

In [41], the authors have proposed a statistical model based on ARIMA time-series forecasting to predict the COVID19 pandemic spread in Saudi Arabia. This study has covered the duration from 2 March 2020 to 20 April 2020. Their model achieved a lower RMSE value of 21.17. However, their model was examined for only the last 15 days. The authors did not present enough comparative analysis with other methods in the literature to prove their findings. They have reported only the actual cases vs. the predicted cases with no deep investigations.

In [42], the authors have proposed a deep learning algorithm based on the LSTM network to predict the COVID19 pandemic spread in Saudi Arabia. This study has covered the duration from 2 March to 10 October 2020. The algorithm was examined at three different periods, and the authors have compared their findings with the other statistical models in the literature. However, they have achieved a high RMSE value, which reached 160.608 with confirmed cases prediction and 100.039 with death cases prediction.

In [43], the authors have proposed a system dynamic based on SIR and SIR-F models to predict the COVID-19 pandemic spread in Saudi Arabia. This study has covered the duration from 2 March to 12 June 2020. In addition, the study has investigated the COVID19 pandemic spread for susceptible, infected and recovered people. However, the authors did not compare their findings with the other statistical or machine learning models in the literature.

In [44], the authors proposed a mathematical model for predicting the confirmed cases in Saudi Arabia based on fractal-fractional derivative. However, they investigate the COVID-19 cases for a brief duration from 1 March 2020 till 22 April 2021. They have not also investigated the model’s ability to be applied to other countries.

In [45], the authors proposed a SIR model to predict the cases in Saudi Arabia, the Philippines, Singapore, and Indonesia. They expand the study from 2 March to 23 December 2020. Their model achieved an R-squared value of 0.95. However, they did not utilize any of the error evaluation metrics to prove their findings.

In [46], the authors proposed a neural network with a self-organizing map for the spatial analysis of data. On the other hand, they utilized the fuzzy fractal technique to capture the temporal trends. They applied their model to Belgium, Italy, United States, and Mexico from 21 January to 30 January 2021. However, this period is concise to measure the model’s performance.

In [47], the authors proposed a hybrid mathematical model based on the short-term forecast (STF) and long-term forecast (LTF) models. Their system was applied to Jordan only. However, their model achieved an R-squared value of 0.84, which is very low related to the similar works in the literature.

In [48], the authors proposed a fractional SEIR-AHQ model to predict COVID-19 cases in Beijing, Chongqing, Tianjin, and Heilongjiang. They investigate the period from 22 January 2020 to 5 March 2020 based on the data of the first 10 (early stage), 20 (middle-stage), and 30 (late-stage) days, respectively. They proved that the non-medical intervention criteria play an essential role in COVID-19 control.

In [49], the authors proposed an LSTM deep learning model to predict COVID-19 spread cases in Egypt from 14 February 2020 to 15 August 2020. Unfortunately, their model was limited to only one country and has not been applied to other countries.

From the previous studies, there is more a real challenge to increase the studying period, investigate the first wave of COVID-19 pandemic spread at its starting and ending, increase R-Squared value and decrease the achieved RMSE value for the prediction results. Also, there is a real need to employ recent deep learning algorithms to increase forecasting performance.

The Encoder–Decoder long short-term memory (LSTM) was introduced for natural language processing (NLP) tasks. The Encoder–Decoder architecture is based on a recurrent neural network (RNN). It achieved a successful performance versus other methods in the literature, specifically in the area of text translation [50]. Recently, the Encoder–Decoder long short-term memory (LSTM) has been applied for several time series forecasting tasks, such as power consumption [51], metal temperature [52], air pollutant [53] behaviour prediction [54], and gas concentration [55]. However, the LSTM core for Encoder–Decoder architecture needs to be developed using recent deep units. Also, there is a real need to apply Encoder–Decoder architecture for pandemic spread prediction tasks.

In this paper, our contributions are as follows: we extend the studying period to cover the pandemic spread from its start point to almost its stopping in Saudi Arabia. We investigate the pandemic spread for (confirmed positive cases, death cases, and recovery cases) in Saudi Arabia. We also investigated the pandemic spread during two periods at the first wave starting and ending. We proposed a deep network model based on encoder-decoder BiLSTM network to forecast the COVID-19 pandemic spread in Saudi Arabia. We optimized the hyper-parameters of our proposed network as follows: the utilized training optimizer, learning rate, units, and dropout ratio. We compared our proposed algorithm with more than fifteen time-series forecasting approaches; statistical, machine learning, and deep learning.

Our proposed system is compared with the other algorithms in the literature that have been applied to the same country in the literature. We examined the generalization of the capability of our proposed system to several countries.The proposed system will need to be applied to more countries to prove its superior performance. The hyper-parameters are critical parameters that need to be developed. Moreover, our proposed system has reflected its superior performance and achieved the best performance compared to the other literature methods.

3. Material and Methods

3.1. Material

In favor of the transparency of the health authorities in Saudi Arabia, there is a daily recording of the COVID-19 cases. The recording includes the confirmed cases, death cases, and recovered cases. The information of all Saudi cases exists on the Saudi ministry of health website (https://covid19.moh.gov.sa, accessed on 1 January 2021). All utilized data were acquired from this website, the existed interactive dashboards, and the available application program interface (API). We utilized the available dataset for the Saudi Arabia dataset at three dimensions of investigation; The first dimension is investigating all COVID-19 cases (confirmed, recovered, and death). The second dimension is specifying a short studying period at the COVID-19 pandemic starting from 2 March 2020 to 31 May 2020. The third dimension specifies long studying period at the COVID-19 pandemic ending from 2 March 2020 to 27 December 2020.

We split the available dataset into 80% for the training process and 20% for testing the proposed model for each studying period. For the short studying period, the training samples are from 22 January 2020 to 17 May 2020, and the testing samples are from 18 May 2020 to 31 May 2020 for about 14 days. For the long studying period, the training samples are from 22 January 2020 to 27 October 2020, and the testing samples are from 28 October 2020 to 27 December 2020 for about 60 days. This splitting percentage reflects the robustness of our proposed network more than the other works in the literature.

To extend our investigations to other countries, we utilized the dataset available at the data repository for the 2019 Novel Coronavirus Visual Dashboard operated by the Johns Hopkins University Center for Systems Science and Engineering (JHU-CSSE). We investigated both deaths and confirmed cases in the countries of interest.

3.2. Methods

In this work, we employ three kinds of methodologies: statistical models, machine learning models, and deep learning models. Our proposed system is described in Figure 1. First, we collect the available COVID-19 time series data set with it three categories (confirmed cases, death cases, and recovered cases) for Saudi Arabia is split into 80% for training and 20% for testing. Then, we utilize a data pre-processing stage, which includes reading CSV files and data standardization of the available data.

In the hyper-parameter selection stage, we select the optimal parameters from the previous studies and investigate the optimal hyper-parameters to train our proposed network. After the training process for each model, we evaluate the prediction results based on statistical performance indices and visualization plots. The Saudi COVID-19 time series data is employed to optimize our proposed network hyper-parameters. Finally, our proposed optimized network is examined by forecasting the global dataset’s death and confirmed cases for several countries.

3.2.1. Data Pre-Processing

After reading the COVID-19 time-series data from a comma -separated values (CSV) files, we apply the feature scaling stage to the input time series. There are two primary data feature scaling methodologies that are applied to the time-series data, which are normalization [56] and standardization [57]. These methodologies play a crucial role in overall system performance, especially the time series data [58]. This study investigates which feature scaling technique is more suitable with the COVID-19 data time series.

The features scaling using normalization technique can be defined as follows:

$$Y_{norm}=\frac{X-X_{min}}{X_{max}-X_{min}}$$

(1)

where $Y_{norm}$ represents the normalized data between 0 and 1 for input data X. $X_{max}$ and $X_{min}$ represents the maximum value and the minimum value of the input X, respectively. The features scaling using standardization technique can be defined as follows:

$$Y_{stand}=\frac{X-\mu }{\sigma }$$

(2)

where $Y_{stand}$ represents the standardized data for the input data X. µ represents the mean value of the given input X and $\sigma$ represents the standard deviation the given input X.

3.2.2. Our Proposed Approach

To investigate the superior performance of our proposed approach, we assess the comparison between our proposed approach and the other COVID-19 forecasting algorithms in the literature, such as machine learning, deep learning, and statistical algorithms. The machine learning regression algorithms are linear regression, support vector regression (SVR), lasso, Huber, RANSAC, TheilSen, and ElasticNet.

We set the parameters of the machine learning algorithms as in [34]. The deep learning algorithms applied previously for COVID-19 time series forecasting were LSTM, BiLSTM, GRU, and Encoder–Decoder BiLSTM. We set the hyper-parameters for each deep learning network as in [32,42]. For the statistical approach, we investigate the performance of three ARIMA models in the literature to investigate the COVID-19 spread in Saudi Arabia. We set the (P,Q,D), which controls the performance of ARIMA model as in [27,41,59].

The traditional RNNs predict each input sample corresponds to an output sample for the same time step. The Encoder–Decoder deep architecture was demonstrated to solve the sequence-to-sequence mapping models [55]. Therefore, it does not correlate the input time series correspond to the predicted time series since the inputs and the predicted samples are not correlated, and their lengths can be different. Therefore, it maps the time series samples of different lengths to each other.

The disadvantage of this model is the challenge to summarize a long sequence into a single vector, and the model often forgets the earlier parts of the input sequence when processing the last parts. Therefore, we develop the traditional Encoder–Decoder model by replacing the LSTM unit by BiLSTM unit model to overcome forgetting earlier samples in the sequence. The Encoder–Decoder deep architecture aims to extract more valuable data representation in compressed form, enabling the network to obtain the most discriminated features of the training input data [54].

An auto-encoder consists of encoding and decoding stages as shown in Figure 2. First, the encoder is utilized to read the input time series data and encode it into a fixed-length vector. The second stage is called a decoder. The decoder is utilized for decoding the fixed-length vector. Then, it produces the forecast sequence. Our proposed approach is based on Encoder-Decoder BiLSTM architecture. The proposed approach consists of a BiLSTM layer for both encoding and decoding layers, RepeatVector layer, dropout layer, and time distributed layer as shown in Figure 2. The Encoder–Decoder is based on Recurrent Neural Networks (RNNs). We employ the deep BiLSTM as a particular type of RNN for encoding and decoding procedures in our proposed approach.

RNNs are networks that are organized into successive layers, and each layer of neurons represents nodes [60]. RNNs are a vanilla neural network consisting of input layers, hidden layers, and output layers where the neurons connect. However, in RNNs, each neuron is assigned to a fixed time step. Furthermore, each neuron in the hidden layer is also connected in a time-dependent direction. Finally, each input and output neuron is attached to the hidden layer with its corresponding time step. RNNs have several advantages for time series forecasting: not significantly affected with missed values, tracking complex time series patterns, and modeling data flexibility as each sample depends on the previous one. However, the vanilla RNN cannot forecast long period time series.

Long short term memory (LSTM) is a special kind of RNN deep learning model that was developed by [61]. LSTM has the advantage of predicting long-term time series. This is because it uses a particular combination of hidden units, element-wise products, and aggregates within units to execute gates responsible for controlling memory cells. First, each cell is created to hold information without alteration for long periods. Then, the learnable weight values must be updated to forecast the next step, which demands the preservation of information from the initial steps. Unfortunately, the simple structure of RNN makes it learn only a limited number of short-term relationships, and it does not perform well during forecasting the long-term series. However, LSTM can well forecast long-term series, and it overcomes the vanishing gradient that happens in the state-of-the-art RNN.

LSTM consists of a cell current state with three gates: input, output, and forget, as shown in Figure 2. The cell state is the network memory responsible for retrieving the sample along the input sequence. The input gate determines the relevant information to add the previous time steps. The forget gate keeps the previous time step and determines how the previous memory remembers and forgets. Finally, the output gate decides the value of the current time step.

The forget gate equation is computed as follows:

$$f\left(t\right)=\sigma (x\left(t\right)U_f+h(t-1)W_f)$$

(3)

The input gate equation is computed as follows:

$$i_1\left(t\right)=\sigma (x\left(t\right)U_i+h(t-1)W_i)$$

(4)

$$i_2\left(t\right)=tanh(x\left(t\right)U_g+h(t-1)W_g)$$

(5)

$$i\left(t\right)=i_1\left(t\right)\ast i_2\left(t\right)$$

(6)

The cell state equation is computed as follows:

$$C\left(t\right)=\sigma \left(f\right(t)\ast C(t-1)+i(t\left)\right)$$

(7)

The output gate equation is computed as follows:

$$O\left(t\right)=\sigma (x\left(t\right)\ast U_O+h(t-1)W_O)$$

(8)

$$h\left(t\right)=tanh\left(C_t\right)\ast O\left(t\right)$$

(9)

where i(t), f(t), and O(t) are the input, forget, and output gates at time t, respectively; $W_i$ and $U_i$ represent the hidden layer weights that are the input of input gate; $W_f$, $U_f$ represent the weights of the hidden layer corresponding to the forget gate; $W_o$ and $U_o$ represent the weights of the hidden layer corresponding to the output gate; and $C_t$ and $h_t$ represent the outcome of the cell and the outcome of the layer, respectively [62].

One of the few limitations founded on the LSTM cell is that it cannot see the future time sample. To overcome this limitation, in [63], the authors introduced the bidirectional-LSTM (BiLSTM) as shown in Figure 2. For the input sequence X(t) with time step t, the output sequence y(t) is calculated from $\overrightarrow{h}\left(t\right)$ in the forward direction, and $\overleftarrow{h}\left(t\right)$ in the backward direction. Therefore, we replace the LSTM unit of the traditional Encoder–Decoder architecture with a BiLSTM unit in our proposed approach to take the advantages of BiLSTM over the LSTM unit.

The RepeatVector layer is essential to fit the encoder data dimension with the decoder data dimension. The RepeatVector is employed to repeat the one fixed-length data for each time step in the output sequence.

In our proposed approach, the RepeatVector layer repeats the outputs of the previous encoder stage for one time. Therefore, for the input shape with size (1, 64), its output shape after the RepeatVector will be (1, 1, 64), since the inputs were repeated only once. The input should be at least 3D, and the dimension of the index one will be considered to be the temporal dimension. Therefore, the 3D output can be processed later through the decoding stage.

After the decoding stage, we apply the dropout layer to the decoding stage output. The dropout layer was introduced by [64] to increase the performance of the neural networks. This happens through dropping some of the weights inside the network as in Equation (10). In addition, the connections of the drop weights are skipped, as shown in Figure 3. The dropout technique supports the regularization strategy, reducing the risk of co-adaptation, and reducing the over-fitting. In this work, we investigate the optimal drop percentage and its effect on our proposed approach.

$${\displaystyle \widehat{{\mathbf{w}}_j}=\left\{\begin{array}{cc}{\mathbf{w}}_j,\hfill & \mathrm{with}P\left(c\right)\hfill \\ \mathbf{0},\hfill & \mathrm{otherwise}\hfill \end{array}\right.}$$

(10)

where $\widehat{{\mathbf{w}}_j}$ represents the output dropout matrix with the dropping probability c for the input weight matrix ${\mathbf{w}}_j$.

The Time-Distributed layer performs the same dense layer function to the output of the dropout layer for one time step at a time. Time-Distributed layer processes each received input sample using the dense layer. Therefore, we apply a dense layer to every temporal slice of the input data with an index considered to be the temporal dimension. For the COVID-19 case time-series forecasting model, the dense layer output is set to be 1.

3.2.3. Prediction Results Evaluation Criteria

We assess the statistical performance in terms of three error measures, which are the mean absolute error (MAE), root mean square error (RMSE), and mean square error (MSE). The three metrics are defined in Equations (11)–(13), respectively. Moreover, we employ the coefficient of determinations (R-Squared) specified in this section for performance evaluation as in Equation (15). The low MSE, RMSE, and MAE values indicate the best forecasting performance. On the other hand, the high R-Squared value indicates the best forecasting performance. We also utilize the visualization plots for training, validation, and prediction plots. This proves the forecasting results trend stability with time-series data variation.

$$MAE=\frac{1}{N}\sum _{i=1}^N\left|C_i-{\widehat{C}}_i\right|$$

(11)

$$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N{(C_i-{\widehat{C}}_i)}^2}$$

(12)

$$MSE=\frac{1}{N}\sum _{i=1}^N{(C_i-\widehat{C_i})}^2$$

(13)

$$R-Squared=1-\frac{\frac{1}{N}{\sum }_{i=1}^N{\left(C_i-{\widehat{C}}_i\right)}^2}{\frac{1}{N}{\sum }_{i=1}^N{\left(C_i-{\overline{C}}_i\right)}^2}$$

(14)

where $C_i$ represents the true value of the COVID-19 cases, ${\widehat{C}}_i$ represents the predicted value of the COVID-19 cases at time i, ${\overline{C}}_i$ is defined as ${\overline{C}}_i=\frac{1}{N}{\sum }_{i=1}^N{C}_i$, and N represents the number of samples.

Furthermore, we utilize the mean absolute scaled error (MASE), which is a scale-free error metric [65]. The low MASE value indicates the best forecasting performance. The error is represented in the MASE metric as a ratio compared to a baseline average error as in Equation (15).

$$\mathrm{MASE}=\frac{1}{I}\sum _{i=1}^I\left|q\left(t_i\right)\right|$$

(15)

$$q\left(t_i\right)=\frac{e\left(t_i\right)}{\frac{1}{N}{\sum }_{n=2}^N\left|X\left(t_n\right)-X\left(t_{n-1}\right)\right|}$$

(16)

where $q\left(t_i\right)$ represents the scaled error value, $e\left(t_i\right)$ represents the error value at time $t_i$ is defined as $e\left(t_i\right)=(C_i-{\widehat{C}}_i)$, I are steps of forecasts, and $X\left(t_n\right)$ = {1, 2, …, N} are the existing observations that used for the training process of the forecasting model.

4. Results

During our experiment, we utilized the software Python 3.8 in the Spyder platform. The hardware system contained Quad-Core 2.9 GHz Intel i7 with 16 GB RAM. The GPU computation was performed through NVIDIA Ge-Force 840M with 4 GB built-in RAM and compute-capability 5.0.

We organized our experiments as follows: Experiment 1 aimed to find the optimal hyper-parameters for our proposed model. It also aims to find the best optimizer to train our proposed model. Experiment 2 examined the effect of features scaling strategy on the proposed model. Experiment 3 proposed a comparison between our proposed method and fifteen forecasting models from the literature. Furthermore, most of these models were employed in forecasting the Saudi Arabia COVID-19 spread. Experiments 1, 2, 3, and 4 were based on the Saudi Arabia COVID-19 dataset for confirmed, recovered, and death cases. Finally, Experiment 4 aims to prove the ability of our proposed model to be applied to other countries, including Brazil, India, South Africa, Germany, Italy, Turkey, and Spain.

4.1. Experiment 1

In Experiment 1, we investigate the optimal parameters to be utilized in the training process of our model. In addition, we investigate the effect of the optimization algorithms of our proposed model. We employ the Saudi Arabia COVID-19 dataset for this target.

Generally speaking, the pandemic wave starts with the first confirmed case reported, and it ends with about zero cases reported before the regrowth of confirmed cases again [66]. For example, in Saudi Arabia, the first case of COVID-19 was reported on 2 March 2020 [67]. With the end of December 2020, the first wave of COVID-19 spread was almost ended as reported in [68,69,70]. Therefore, during the studying period from 2 March 2020 until 27 December 2020, there was only one COVID-19 spread wave in Saudi Arabia, as reported in the literature. In our experiments, we investigate the model performance during two different time series in the start of the first wave and the ending of the first wave. Furthermore, we also investigate the confirmed, recovered, and death cases in all periods.

To optimize the hyper-parameters of our model, we investigate three model parameters: the hidden units, initial learning rate value, and drop percentage as in [42]. First, we fix the initial learning rate at 0.0005 and the drop percentage at 0.2. Then, we tune the number of hidden units with (16, 32, 64, and 128). Secondly, we fix the hidden unit at 128 and the drop percentage at 0.2. Then, we tune the initial learning rate value with (0.01, 0.001, 0.005, and 0.0005). Thirdly, we fix the units at 128 and the initial learning rate at 0.0005. Then, we tune the dropout percentage with values (0.1, 0.2, 0.3, and 0.4). We set the (Adaptive Moment Estimation Algorithm) Adam optimizer as the training optimizer for this experiment as it proved its superior performance for forecasting the same task [32,42,66]. Therefore, we can deliver the optimal parameters for our proposed model.

For the forecasting of confirmed cases of the COVID-19 first wave starting as shown in Table S1, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9929, 9.681 × 10², 8.031 × 10², and 9.371 × 10⁵, respectively. For the forecasting of confirmed cases of the COVID-19 first wave ending as shown in Table S1, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage.

It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9981, 2.001 × 10², 1.711 × 10², and 3.981 × 10⁴, respectively. The superior performance of our proposed model to forecast the confirmed cases for along 14 days during the first wave start is shown in Figure 4a. The superior performance of our proposed model to forecast the confirmed cases for along 60 days during the first wave end is shown in Figure 4b.

For the forecasting of recovered cases of COVID-19 first wave starting as shown in Table S2, our proposed model achieved the best performance with the following parameters: 16 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9657, 2.271 × 10³, 2.081 × 10³, and 5.171 × 10⁶, respectively. For the forecasting of recovered cases of COVID-19 first wave ending as shown in Table S2, our proposed model achieved the best performance with the following parameters: 16 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9984, 2.451 × 10², 2.141 × 10², and 6.021 × 10⁴, respectively.

The performance of our proposed model to forecast the recovered cases for along 14 days during the first wave starting is shown in Figure 4c. We noticed that the instability of the prediction line vs. the test line. The superior performance of our proposed model to forecast the recovered cases for along 60 days during the first wave ending is shown in Figure 4d.

For the forecasting of death cases of COVID-19 first wave starting as shown in Table S3, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9954, 3.96, 3.40, and 10.56, respectively. For the forecasting of death cases of COVID-19 first wave ending as shown in Table S3, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage.

It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9989, 8.00, 6.55, and 60.39, respectively. The superior performance of our proposed model to forecast the death cases for along 14 days during the first wave starting is shown in Figure 4e. The superior performance of our proposed model to forecast the death cases for along 60 days during the first wave end is shown in Figure 4f.

In the literature, there is no deep investigation of the optimizer effect on the forecasting model [32,42]. Therefore, we investigate the effect of the optimization algorithms of our proposed model as shown in Tables S4 and S5 through two studies.

In study 1, we investigated the performance of different optimizers with the median hyper-parameter values, which are 48 hidden units, initial learning rate 0.003, and dropout percentage 0.25. We employ the median values introduced in the first experiment to have a fair comparison between the different optimizers. In study2, we investigate the performance of our proposed optimal hyper-parameters, which are 128 hidden units, initial learning rate 0.0005, and dropout percentage 0.2 with different optimizers. These optimization algorithms are Root Mean Square Propagation (RMSprop), Stochastic Gradient Descent (SGD), Adam, Adaptive Maximum (Adamax), Adaptive Delta (Adadelta), and Adaptive Gradient (Adagrad).

We selected the confirmed cases during the first wave starting and ending to investigate the optimizer effect. In study1, Adam optimizer achieved the best forecasting performance for forecasting the confirmed cases during the first wave start. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9908, 1.101 × 10³, 9.281 × 10², and 1.221 × 10⁶, respectively. For the forecasting of the confirmed cases during the first wave ending, Adam optimizer achieved the best forecasting performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9789, 6.581 × 10², 6.411 × 10², and 4.321 × 10⁵, respectively.

In study 2, we investigated the performance of different optimizers with the optimal hyper-parameters values, which are 128 hidden units, initial learning rate 0.0005, and dropout percentage 0.2. Adam optimizer achieved the best forecasting performance for forecasting the confirmed cases during the first wave start. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9929, 9.681 × 10², 8.031 × 10², and 9.371 × 10⁵, respectively. For the forecasting of the confirmed cases during the first wave ending, Adam optimizer achieved the best forecasting performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9981, 2.001 × 10², 1.711 × 10², and 3.981 × 10⁴, respectively. Adamax optimizer achieved the second rate.

It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9783, 6.661 × 10², 6.381 × 10², and 4.431 × 10⁵, respectively. RMSprop optimizer achieved the third rate. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9783, 6.661 × 10², 6.381 × 10², and 4.431 × 10⁵, respectively. However, SGD, Adadelta, and Adagrad achieved bad performance during our forecasting task. From this experiment, we concluded the importance of the three investigated hyper-parameters to control the performance of the proposed forecasting model. We noticed that the initial learning rate value variation was a critical parameter on our proposed model’s performance. However, the variation of hidden units or dropout percentage slightly varied our proposed model’s performance. It is recommended to set the parameter for our model with 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage.

We also noticed the stability of our model between the forecasting of the short time series during the COVID-19 first wave starting and forecasting of the long time series during the COVID-19 first wave ending. The exact optimal hyper-parameters were the same for both first wave start and end.

Generally, the forecasting performance during the first wave end was better than during the first wave start. This can be explained by the high variation of all cases during the first wave start. Based on several statistical assessments, the forecasting performance was similar for the confirmed, recovered, and death cases with shallow variance errors.

We also concluded that the Adam optimizer was the best optimization algorithm to train our proposed model. A similar performance between the two studies was noticed for the performance of different optimizers. The Adamax optimizer achieved the second rank among optimizers. However, SGD, Adadelta, and Adagrad achieved a lousy performance during our forecasting task. We also noticed that the crucial role of investigating the optimization algorithms for deep learning models.

4.2. Experiment 2

In this experiment, we investigate the performance of our proposed model based on different feature scaling techniques. We investigate the features scaling techniques normalization and standardization. Based on several statistical assessment criteria, the features standardization is better than the features normalization technique as shown in Figure 5a,b. We select the confirmed cases during the first wave starting and ending to investigate the features scaling effect. We set the hyper-parameters as 128 hidden units, initial learning rate 0.0005, and dropout percentage 0.2.

We noticed that the standardization has a higher R-Squared than the normalization technique. We also noticed that the standardization has lower RMSE, MSE, and MAE than the normalization technique. Therefore, we conclude that the standardization is much better than the normalization technique for our proposed model.

4.3. Experiment 3

This experiment compares our proposed approach and four deep learning methods, three statistical ARIMA models, seven machine learning algorithms, and one prophet model.

ARIMA model is a statistical time-series method that has been extensively utilized to forecast infectious diseases. It gathers the Autoregressive (AR) model and Moving Average (MA) with integration based on the decomposition method. In which all current and historical residual series values in the present time series are expressed linearly [41]. The ARIMA model is controlled as ARIMA (p,d,q) values. p represents the auto-regressive seasonal order, d represents the autoregressive non-seasonal order, and q represents the non-seasonal moving average order. The prophet model was utilized to forecast COVID-19 cases in Saudi Arabia [71]. The prophet model utilized the Fourier spectral data to compute the seasonality impact for forecasting time series data in the medium-term and long-term [72].

This model has been successfully utilized in Saudi Arabia to predict COVID-19 cases. The supervised machine learning algorithms were successfully applied to predict COVID-19 cases based on regression models [34]. The deep learning RNN networks can automatically capture seasonal and trends characteristics of the time series [71]. In this experiment, we utilized four different deep architectures for time-series forecasting tasks: LSTM, BiLSTM, GRU, and Encoder–Decoder LSTM.

For each model parameters, we set their values as in the literature. The BiLSTM model was applied to Saudi Arabia [42]. The BiLSTM model hyper-parameters were set as 100 hidden units and 0.005 initial learning rate value. The GRU and LSTM hyper-parameters were generally applied in several countries [32]. They set the parameters as 128 hidden units, 0.001 initial learning rate value. Furthermore, we compared our proposed method based on the Encoder–Decoder BiLSTM model and the Encoder–Decoder LSTM one.

We set the optimal parameters as 128 units, 0.2 drop percentage, and 0.0005 initial learning rate for both models. In this way, we have a fair comparison to study the effect of using BiLSTM unit in the encoder decoder deep architecture. We utilized the Adam optimizer for all deep learning models training process. The three ARIMA models were applied to Saudi Arabia. Each model’s parameters were investigated previously in these studies. The (p,q,d) values were set in ARIMA Model1 [41] as (2,1,1) in ARIMA Model2 [27] as (1,1,1) in ARIMA Model3 [59] as (0,2,0). The machine learning models setup parameters were utilized as in [34] and [32].

For the forecasting of confirmed cases of COVID-19 first wave starting as shown in

Table 1. A comparative analysis of our proposed approach with the previous methods based on statistical evaluation criteria for confirmed cases in Saudi Arabia.

	COVID-19 First Wave Starting				COVID-19 First Wave Ending
Forecasting Model	MSE	MAE	RMSE	R-squared	MSE	MAE	RMSE	R-squared
Prophet Model	4.721 × 109	6.761 × 104	6.871 × 104	−31.8956	1.271 × 1011	3.561 × 105	3.561 × 105	−6090.89
ARIMA_Model1	1.301 × 107	3.251 × 103	3.611 × 103	0.909	4.711 × 106	1.541 × 103	2.171 × 103	0.7733
ARIMA_Model2	1.251 × 107	3.181 × 103	3.541 × 103	0.9126	3.131 × 106	1.301 × 103	1.771 × 103	0.8493
ARIMA_Model3	1.231 × 107	3.141 × 103	3.501 × 103	0.9143	7.141 × 106	1.851 × 103	2.671 × 103	0.6566
Lasso	2.831 × 109	4.671 × 104	5.321 × 104	−70.9598	2.121 × 1010	1.451 × 105	1.451 × 105	−1578.24
RANSACRegressor	5.581 × 109	6.741 × 104	7.471 × 104	−140.6868	8.111 × 108	2.841 × 104	2.851 × 104	−59.5539
HuberRegressor	4.771 × 109	6.191 × 104	6.901 × 104	−120.0091	5.351 × 108	2.311 × 104	2.311 × 104	−38.9581
LinearRegression	4.211 × 109	5.791 × 104	6.491 × 104	−105.8295	8.111 × 108	2.841 × 104	2.851 × 104	−59.5539
SVR_linear	4.991 × 109	6.361 × 104	7.071 × 104	−125.7476	3.391 × 108	1.831 × 104	1.841 × 104	−24.336
ElasticNet	5.871 × 108	1.941 × 104	2.421 × 104	−13.9027	8.611 × 109	9.281 × 104	9.281 × 104	−641.912
TheilSenRegressor	8.431 × 109	8.411 × 104	9.181 × 104	−212.9676	1.411 × 109	3.761 × 104	3.761 × 104	−104.638
GRU	6.661 × 106	2.041 × 103	2.581 × 103	0.9497	6.381 × 107	7.961 × 103	7.991 × 103	−2.1222
BiLSTM	3.631 × 107	5.791 × 103	6.031 × 103	0.726	2.321 × 105	4.511 × 102	4.821 × 102	0.9886
LSTM	2.951 × 107	4.711 × 103	5.431 × 103	0.5619	1.261 × 107	3.401 × 103	3.551 × 103	0.3164
Encoder–Decoder-LSTM	2.501 × 106	1.531 × 103	1.581 × 103	0.9812	8.921 × 106	2.981 × 103	2.991 × 103	0.5636
Our Proposed Model	9.371 × 105	8.031 × 102	9.681 × 102	0.9929	3.981 × 104	1.711 × 102	2.001 × 102	0.9981

, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9929, 9.681 × 10², 8.031 × 10², and 9.371 × 10⁵, respectively. The Encoder–Decoder LSTM model achieved the second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9812, 1.5811 × 10³, 1.5311 × 10³, and 2.501 × 10⁶, respectively. GRU model achieved the third rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9497, 2.041 × 10³, 6.661 × 10⁶, and 2.581 × 10³, respectively.

ARIMA Model3 achieved the third rank among all compared models and the first rank in all ARIMA Models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9143, 3.501 × 10³, 3.141 × 10³, and 1.231 × 10⁷, respectively. Among deep learning models, the LSTM model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.562, 5.431 × 10³, 4.711 × 10³, and 2.951 × 10⁷, respectively. For all machine learning models, we noticed bad forecasting results with high error values and negative coefficients of determination. TheilSen regression algorithm achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −212.97, 9.181 × 10⁴, 8.141 × 10⁴, and 8.431 × 10⁹, respectively.

For the forecasting of the confirmed cases of COVID-19 first wave ending as shown in

Table 2. COVID-19 spread forecasting results in countries with double spread wave based statistical evaluation criteria.

Country	Cases	Period	MSE	MAE	RMSE	MASE	R_Squared
Germany	Death Cases	Period 1	5.701 × 102	2.081 × 101	2.391 × 101	12.2	0.9881
		Period 2	9.891 × 105	8.151 × 102	9.941 × 102	83.27	0.9723
	Confirmed Cases	Period 1	2.931 × 105	4.541 × 102	5.411 × 102	8.79	0.9633
		Period 2	4.701 × 109	6.821 × 104	6.851 × 104	45.3	0.9595
Italy	Death Cases	Period 1	4.091 × 103	4.931 × 101	6.391 × 101	0.2	0.9904
		Period 2	1.541 × 105	3.021 × 102	3.921 × 102	4.87	0.9987
	Confirmed Cases	Period 1	2.581 × 105	4.581 × 102	5.081 × 102	0.06	0.9742
		Period 2	3.951 × 109	6.181 × 104	6.281 × 104	33.09	0.9797
Turkey	Death Cases	Period 1	3.651 × 102	1.601 × 101	1.911 × 101	2.57	0.9864
		Period 2	9.601 × 104	2.271 × 102	3.101 × 102	0.69	0.9884
	Confirmed Cases	Period 1	8.221 × 105	7.531 × 102	9.071 × 102	0.87	0.976
		Period 2	1.221 × 1010	6.591 × 104	1.101 × 105	47.88	0.9731
Spain	Death Cases	Period 1	3.321 × 103	4.911 × 101	5.761 × 101	0.96	0.9898
		Period 2	5.721 × 105	5.801 × 102	7.561 × 102	3.47	0.9707
	Confirmed Cases	Period 1	8.461 × 105	8.381 × 102	9.201 × 102	0.71	0.9141
		Period 2	2.421 × 109	4.871 × 104	4.921 × 104	7.93	0.9418

, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9981, 21 × 10², 1.711 × 10², and 3.981 × 10⁴, respectively. BiLSTM model achieved the second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9886, 4.821 × 10², 4.511 × 10², and 2.321 × 10⁵, respectively. In spite of, BiLSTM model is slightly higher than our proposed model, the BiLSTM model has higher error values.

The Encoder–Decoder LSTM model did not achieve a good performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.5636, 2.981 × 10³, 2.991 × 10³, and 8.921 × 10⁶, respectively. ARIMA Model2 achieved the third rank among all compared models and the first rank in all ARIMA Models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8493, 1.771 × 10³, 1.301 × 10³, and 3.131 × 10⁶, respectively. Among the deep learning models, the GRU model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −2.1222, 7.991 × 10³, 7.961 × 10³, and 6.381 × 10⁷, respectively. For all examined machine learning models, we noticed bad forecasting results with high error values and negative coefficients of determination.

In addition, the prophet model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −6090.89, 3.561 × 10⁵, 3.561 × 10⁵, and 1.271 × 10¹¹, respectively. The performance of our proposed model to forecast the confirmed cases for along 14 days during the first wave start is shown in Figure 6a. We noticed superior performance of our proposed model to forecast the confirmed cases for along 60 days during the first wave end as shown in Figure 6b.

For the forecasting of recovered cases of COVID-19 first wave starting are shown in

Table 3. A comparative analysis of our proposed approach with the previous methods based on statistical evaluation criteria for recovered cases in Saudi Arabia.

	COVID-19 First Wave Starting				COVID-19 First Wave Ending
Forecasting Model	MSE	MAE	RMSE	R-squared	MSE	MAE	RMSE	R-squared
Prophet Model	1.771 × 109	4.001 × 104	4.211 × 104	−9.3176	1.191 × 1011	3.451 × 105	3.451 × 105	−3067.49
ARIMA_Model1	5.541 × 107	5.951 × 103	7.441 × 103	0.6777	3.141 × 107	4.271 × 103	5.611 × 103	0.1881
ARIMA_Model2	3.681 × 107	4.731 × 103	6.061 × 103	0.7862	6.811 × 107	6.401 × 103	8.251 × 103	−0.7588
ARIMA_Model3	3.041 × 107	4.341 × 103	5.511 × 103	0.8234	1.111 × 108	8.541 × 103	1.051 × 104	−1.8675
Lasso	1.541 × 107	3.271 × 103	3.921 × 103	0.7733	2.681 × 1010	1.641 × 105	1.641 × 105	−976.727
RANSACRegressor	1.811 × 108	1.341 × 104	1.351 × 104	−1.6702	1.381 × 109	3.711 × 104	3.711 × 104	−49.2069
HuberRegressor	2.051 × 107	3.731 × 103	4.531 × 103	0.6984	1.051 × 109	3.231 × 104	3.231 × 104	−37.0974
LinearRegression	2.381 × 107	3.991 × 103	4.881 × 103	0.649	1.381 × 109	3.711 × 104	3.711 × 104	−49.2069
SVR_linear	3.831 × 107	4.961 × 103	6.191 × 103	0.4357	6.911 × 108	2.631 × 104	2.631 × 104	−24.1915
ElasticNet	1.591 × 108	1.261 × 104	1.261 × 104	−1.3486	9.541 × 109	9.761 × 104	9.771 × 104	−346.681
TheilSenRegressor	4.141 × 107	5.191 × 103	6.431 × 103	0.3905	5.421 × 109	7.361 × 104	7.361 × 104	−196.374
GRU	1.941 × 107	4.361 × 103	4.411 × 103	0.8713	1.951 × 107	4.401 × 103	4.421 × 103	0.4877
BiLSTM	1.211 × 108	1.011 × 104	1.101 × 104	0.197	8.571 × 104	2.471 × 102	2.931 × 102	0.9978
LSTM	9.451 × 107	7.871 × 103	9.721 × 103	−0.553	9.401 × 106	2.781 × 103	3.071 × 103	0.7271
Encoder–Decoder-LSTM	2.031 × 107	3.871 × 103	4.511 × 103	0.8653	5.311 × 106	2.281 × 103	2.301 × 103	0.8608
Our Proposed Model	5.171 × 106	2.081 × 103	2.271 × 103	0.9657	6.021 × 104	2.141 × 102	2.451 × 102	0.9984

, Our proposed model achieved the best performance with the following parameters: 16 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.965, 2.271 × 10³, 2.081 × 10³, and 5.171 × 10⁶, respectively. The GRU model achieved the second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8713, 4.411 × 10³, 4.361 × 10³, and 1.941 × 10⁷, respectively. The Encoder–Decoder LSTM model achieved a good performance relative to the other methods. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8653, 4.511 × 10³, 3.8711 × 10³, and 2.031 × 10⁷, respectively.

ARIMA Model3 achieved the third rank among all compared models and the first rank in all ARIMA Models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8234, 5.511 × 10³, 4.341 × 10³, and 3.041 × 10⁷, respectively. Among deep learning models, the LSTM model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −0.553, 9.721 × 10³, 7.871 × 10³, and 9.451 × 10⁷, respectively. Among machine learning models, the Lasso regression model achieved the highest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.773, 3.921 × 10³, 3.271 × 10³, and 1.541 × 10⁷, respectively. We noticed bad forecasting results with high error values and negative coefficients of determination. Prophet model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −9.3176, 4.211 × 10⁴, 4.001 × 10⁴, and 1.771 × 10⁹, respectively.

For the forecasting of recovered cases of COVID-19 first wave ending as shown in Table 3, our proposed model achieved the best performance with the following parameters: 16 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9984, 2.451 × 10², 2.141 × 10², and 6.021 × 10⁴, respectively. BiLSTM model achieved the second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9978, 2.931 × 10², 2.471 × 10², and 8.571 × 10⁵, respectively. The Encoder–Decoder LSTM model achieved the third rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8608, 2.311 × 10³, 2.2811 × 10³, and 5.311 × 10⁶, respectively.

The three ARIMA Models achieved a low performance among all compared models in forecasting of recovered cases of long time-series. The best ARIMA Model is Model1. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.1881, 5.611 × 10³, 4.271 × 10³, and 3.141 × 10⁷, respectively. The BiLSTM model achieved the highest performance among deep learning models and the second rank among all compared models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9978, 2.931 × 10², 2.471 × 10², and 8.571 × 10⁴, respectively. We noticed bad forecasting results for all machine learning models with high error values and negative coefficients of determination.

In addition, the prophet model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −3067.49, 3.451 × 10⁵, 3.451 × 10⁵, and 6.811 × 10⁷, respectively. The performance of our proposed model to forecast the recovered cases for along 14 days during the first wave starting is shown in Figure 6c. We noticed the superior performance of our proposed model to forecast the confirmed cases for along 60 days during the first wave ending as shown in Figure 6d.

For the forecasting of death cases of COVID-19 first wave starting as shown in

Table 4. A comparative analysis of our proposed approach with the previous methods based statistical evaluation criteria for death cases in Saudi Arabia.

	COVID-19 First Wave Starting				COVID-19 First Wave Ending
Forecasting Model	MSE	MAE	RMSE	R-squared	MSE	MAE	RMSE	R-squared
Prophet Model	80,698	276.6609	284.0755	−18.3656	3.281 × 107	5725.99	5731.193	−549.971
ARIMA_Model1	165.524	8.139	12.8656	0.9603	1.361 × 104	85.9144	116.7805	0.7712
ARIMA_Model2	155.3177	7.9651	12.4627	0.9627	1.111 × 104	73.9975	105.474	0.8134
ARIMA_Model3	137.216	7.1448	11.7139	0.9671	4.911 × 104	167.2461	221.6259	0.1761
Lasso	1150.8528	32.9337	33.9242	0.4582	8.651 × 106	2936.2	2941.442	−202.413
RANSACRegressor	2396.0025	47.5562	48.949	−0.128	5.061 × 105	703.9295	711.3375	−10.8963
HuberRegressor	3396.2835	56.6015	58.2776	−0.5989	5.091 × 105	705.7586	713.593	−10.9718
LinearRegression	4308.111	63.7156	65.6362	−1.0281	5.061 × 105	703.9295	711.3375	−10.8963
SVR_linear	4769.5801	67.4808	69.0621	−1.2454	5.501 × 105	734.0821	741.8875	−11.94
ElasticNet	5292.7366	71.4414	72.7512	−1.4916	3.461 × 106	1858.153	1860.545	−80.384
TheilSenRegressor	5913.6893	75.4037	76.9005	−1.784	8.541 × 105	919.0879	923.8529	−19.0662
GRU	32.7193	4.8849	5.7201	0.9904	5.921 × 104	240.6431	243.2687	−0.0133
BiLSTM	677.8149	24.9696	26.0349	0.8016	4.241 × 101	5.7041	6.512	0.9993
LSTM	1030.9777	30.1047	32.1088	0.0894	4.121 × 103	52.8893	64.1746	0.9192
Encoder–Decoder-LSTM	87.7449	9.3539	9.3672	0.97	1.431 × 104	119.71	119.73	0.754
Our Proposed Model	16.1006	3.1147	4.0126	0.9953	6.391 × 101	6.5497	7.9967	0.999

, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.995, 4.01, 3.11, and 16.1, respectively. On the other hand, the GRU model achieved second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.967, 11.7139, 7.1448, and 137.216, respectively. The Encoder–Decoder LSTM model achieved the third rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9743, 9.3672, 9.3539, and 87.7449, respectively.

ARIMA Model3 achieved the third rank among all compared models and the first rank in all ARIMA Models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9671, 11.71, 7.14, and 137.22, respectively. Among deep learning models, the LSTM model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −0.553, 9.721 × 10³, 7.871 × 10³, and 9.451 × 10⁷, respectively. In addition, we noticed bad forecasting results with high error values and negative coefficients of determination. However, the Lasso regression model achieved acceptable results. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.45, 32.93, 33.92, 8.431 × 10⁹, and 1150.8, respectively.

For the forecasting of death cases of COVID-19 first wave ending as shown in Table 4, our proposed model achieved the best performance with the following parameters: 128 hidden units, 0.0005 initial learning rate, and 0.2 dropout percentage. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9984, 2.451 × 10², 2.141 × 10², and 6.021 × 10⁴, respectively. BiLSTM model achieved the second rank. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9978, 2.931 × 10², 2.471 × 10², and 8.571 × 10⁴, respectively. In spite of, BiLSTM model is slightly higher than our proposed model, the BiLSTM model has higher error values.

ARIMA Model1 achieved the first rank in all ARIMA Models. However, it is low compared to the other models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.1881, 5.611 × 10³, 4.271 × 10³, and 3.141 × 10⁷, respectively. Among deep learning models, GRU model achieved a good performance compared to the other deep learning and all models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.4877, 4.421 × 10³, 4.401 × 10³, and 1.951 × 10⁷, respectively.

The Encoder–Decoder LSTM model achieved a bad performance relative to the other deep learning models. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.7546, 119.7271, 119.7084, and 1.431 × 10⁴, respectively. For all machine learning models, we noticed bad forecasting results with high error values and negative coefficients of determination. In addition, the prophet model achieved the lowest performance. It demonstrated an R-Squared, RMSE, MAE, and MSE of −3067.4922, 3.451 × 10⁵, 3.451 × 10⁵, and 6.811 × 10⁷, respectively. The performance of our proposed model to forecast the death cases for along 14 days during the first wave starting is shown in Figure 6e. We noticed the superior performance of our proposed model to forecast the confirmed cases for along 60 days during the first wave ending as shown in Figure 6f.

From these results, we noticed that our proposed encoding–decoding deep learning network based on BiLSTM was better than the other deep learning networks. Moreover, our proposed method achieved a higher performance than the traditional Encoder–Decoder LSTM model. We noticed that the performance of the Encoder–Decoder model during the first wave starting was better than its performance during the first wave ending. This can be explained as, in each time step, our proposed method generates a point in forecasting, and it searches for a set of time series positions in the source of the time series where the most relevant information is intensified.

Thus, the model can forecast a target time series point based on the context vectors associated with these source positions and all the previous generated target points. We noticed that the BiLSTM model performance is better than LSTM model. Despite the GRU model performing well in the forecasting of the confirmed cases for short time series, it presented bad performance with forecasting of the confirmed cases for long time series. We noticed a similar performance between forecasting models for confirmed and death cases.

We observed the variation of the LSTM performance between forecasting of the short time-series and long time-series. ARIMA models are a powerful statistical forecasting tools. However, their (p,q,d) values control their performance. For time-series forecasting evaluation criteria, it is crucial to take into consideration the value of R-Squared with the achieved error values.

4.4. Experiment 4

In this experiment, we utilized our proposed model to perform COVID-19 forecasting in two categories of countries. Generally, the pandemic spread wave was followed by significant confirmed and death cases. This wave over-time is visualized as an exponential growth in the cumulative epidemic graph. In our study, we investigate how our model will perform with different spread forms. The first category A includes the counties with a single COVID-19 spread wave, which are Brazil, India, South Africa, and Saudi Arabia during the studying period, as shown in

Table 5. COVID-19 spread forecasting results in countries with single spread wave based statistical evaluation criteria.

Country	Cases	Period	MSE	MAE	RMSE	MASE	R_Squared
Brazil	Death Cases	Period 1	9.581 × 105	9.691 × 102	9.791 × 102	0.96	0.9581
		Period 2	1.271 × 106	1.111 × 103	1.131 × 103	3.47	0.9863
	Confirmed Cases	Period 1	4.341 × 107	5.881 × 103	6.591 × 103	0.71	0.9942
		Period 2	1.211 × 1010	1.101 × 105	1.101 × 105	7.93	0.9677
India	Death Cases	Period 1	4.911 × 103	6.281 × 101	7.011 × 101	3.8	0.9922
		Period 2	8.071 × 105	8.671 × 102	8.981 × 102	2.63	0.9879
	Confirmed Cases	Period 1	1.501 × 108	9.661 × 103	1.231 × 104	3.72	0.8662
		Period 2	1.841 × 1010	1.141 × 105	1.361 × 105	3.59	0.9551
South Africa	Death Cases	Period 1	1.211 × 103	2.331 × 101	3.471 × 101	4.15	0.9334
		Period 2	4.851 × 104	2.041 × 102	2.201 × 102	16.04	0.9878
	Confirmed Cases	Period 1	2.411 × 106	1.501 × 103	1.551 × 103	2.31	0.9237
		Period 2	1.711 × 108	1.201 × 104	1.311 × 104	5.56	0.9691
Saudi Arabia	Death Cases	Period 1	6.391 × 101	6.551 × 100	8.001 × 100	1.78	0.9989
		Period 2	1.561 × 101	3.401 × 100	3.961 × 100	0.46	0.9954
	Confirmed Cases	Period 1	9.371 × 105	8.031 × 102	9.681 × 102	2.24	0.9929
		Period 2	3.981 × 104	1.711 × 102	2.001 × 102	0.18	0.9981

.

The second category B includes the counties with double COVID-19 spread waves, which are Germany, Italy, Turkey, and Spain as shown in Table 2. In this experiment, we also utilize our proposed model to investigate both deaths and confirmed cases in two different periods, Period 1 starting from 22 January 2020 to 31 May 2020, and Period 2 starting from 22 January 2020 to 27 December 2020. Each period is split into 80% for training and 20% for testing.

As shown in Table 5, the forecasting results of Saudi Arabia’s death cases achieved the best performance. This country demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9989, 63.9, 6.55, and 8, respectively. Brazil achieved the second rank of the confirmed cases forecasting. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9942, 6.591 × 10³, 5.881 × 10³, and 4.341 × 10⁷, respectively. India achieved the lowest forecasting performance for confirmed cases during period 1. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.8662, 1.231 × 10⁴, 9.661 × 10³, and 1.501 × 10⁸, respectively.

As shown in Table 2, the forecasting results of Italy’s death cases achieved the best performance. During Period 2, it had R-Squared, RMSE, MAE, and MSE of 0.9987, 3.921 × 10², 3.021 × 10², and 1.541 × 10⁵, respectively. During Period 1, it had R-Squared, RMSE, MAE, and MSE of 0.9904, 60.39, 40.93, and 4.091 × 10³, respectively. Turkey achieved good forecasting results for death cases during period2. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9884, 3.101 × 10², 2.271 × 10², and 9.601 × 10⁴, respectively. Spain achieved the lowest forecasting performance for confirmed cases during period1. It demonstrated an R-Squared, RMSE, MAE, and MSE of 0.9141, 9.201 × 10², 8.381 × 10², and 8.461 × 10⁵, respectively.

Saudi Arabia achieved the lowest MASE error value of 0.18 for the forecasting of confirmed cases in period 2. South Africa achieved the highest MASE error value of 16.04 for the forecasting of death cases in period 2. Italy achieved the lowest MASE error value of 0.2 for the forecasting of death cases in period 1. Germany achieved the highest MASE error value of 83.27 for the forecasting of death cases in period 2.

We visualize the forecasting graphs for the second countries category, which are Turkey, Italy, Germany, and Spain. All forecasting results are shown in Figure 7a–d. Our proposed model achieved high performance to forecast the confirmed cases for along 14 days during the first spread wave for all four countries. On the other hand, it also achieved a high performance to forecast the Germany, Italy, and Spain confirmed cases during for 60 days during the second spread wave. However, the Turkey forecasting results achieved a low performance during the last 15 days.

From Experiment 4, we have proven the generalization capabilities of our proposed model and its superior performance. Our proposed model was able to forecast both deaths and confirmed cases for eight different countries, according to these counties’ variation of control policies. Some of them witnessed two spread waves during the same year, and some have witnessed a single spread wave. Our proposed model can forecast both confirmed and death cases for the two kinds of spread models. In addition, we noticed the superior performance of our proposed model to forecast the confirmed cases for 60 days during the first and second waves. We conclude that the countries from category A achieved lower MASE values than those from category B.

5. Conclusions

In this paper, an Encoder–Decoder BiLSTM deep learning model was utilized to forecast the COVID-19 confirmed cases, recovered cases, and death cases. First, we optimized the hyper-parameters for our proposed model based on Saudi Arabia’s reported cases. Primarily, our proposed model optimal values of the initial learning rate, hidden units, and dropout percentage were 0.0005, 128, and 0.2, respectively. We also investigated the effect of the features scaling technique, and we proved that the standardization technique was more suitable with COVID-19 time series forecasting performance.

We also investigated the COVID-19 spread through the year for two different periods. The forecasting evaluation was assessed based on 14 days and 60 days. We compared our proposed model and the fifteen previous methods in the literature, and our proposed model was the most stable. The Encoder–Decoder BiLSTM achieved higher performance than the traditional Encoder–Decoder BiLSTM model. Moreover, our proposed model achieved higher performance than traditional machine learning techniques, ARIMA models, and deep learning models applied to the same countries.

Our proposed model was applied to forecast several countries that witnessed a single-wave spread, Saudi Arabia, India, Brazil, South Africa, as well as double-wave spreads, Italy, Spain, Germany, and Turkey. Our proposed model achieved lower RMSE, MSE, and MAE values. Furthermore, it achieved the highest R-Squared values. For future work, it is recommended to study the capability of our proposed model to simulate the COVID-19 pandemic spread with different policies and lockdown regulations. Furthermore, there is a real need to investigate the recent hyper-parameter optimization techniques to speed up hyper-parameter selection. Our proposed model will be applied to other pandemic forecasting problems.