key: cord-0807132-dlltv6g0
authors: Alabdulrazzaq, Haneen; Alenezi, Mohammed N.; Rawajfih, Yasmeen; Alghannam, Bareeq A.; Al-Hassan, Abeer A.; Al-Anzi, Fawaz S.
title: On the accuracy of ARIMA based prediction of COVID-19 spread
date: 2021-07-15
journal: Results Phys
DOI: 10.1016/j.rinp.2021.104509
sha: a5218e10debded04ccb91388710a021cdcc28f32
doc_id: 807132
cord_uid: dlltv6g0

COVID-19 was declared a global pandemic by the World Health Organization in March 2020, and has infected more than 4 million people worldwide with over 300,000 deaths by early May 2020. Many researchers around the world incorporated various prediction techniques such as Susceptible–Infected–Recovered model, Susceptible–Exposed–Infected–Recovered model, and Auto Regressive Integrated Moving Average model (ARIMA) to forecast the spread of this pandemic. The ARIMA technique was not heavily used in forecasting COVID-19 by researchers due to the claim that it is not suitable for use in complex and dynamic contexts. The aim of this study is to test how accurate the ARIMA best-fit model predictions were with the actual values reported after the entire time of the prediction had elapsed. We investigate and validate the accuracy of an ARIMA model over a relatively long period of time using Kuwait as a case study. We started by optimizing the parameters of our model to find a best-fit through examining auto-correlation function and partial auto correlation function charts, as well as different accuracy measures. We then used the best-fit model to forecast confirmed and recovered cases of COVID-19 throughout the different phases of Kuwait’s gradual preventive plan. The results show that despite the dynamic nature of the disease and constant revisions made by the Kuwaiti government, the actual values for most of the time period observed were well within bounds of our selected ARIMA model prediction at 95% confidence interval. Pearson’s correlation coefficient for the forecast points with the actual recorded data was found to be 0.996. This indicates that the two sets are highly correlated. The accuracy of the prediction provided by our ARIMA model is both appropriate and satisfactory.

The first case of the novel coronavirus known as SARS-CoV-2, and commonly referred to as COVID-19, was reported in December 2019, in the city of Wuhan in China [1] . The virus was spreading rapidly and within two months of its emergence, the World Health Organization (WHO) declared it an international pandemic. From the time that COVID-19 was first identified, until early May 2020, it had spread to 216 countries and infected over 4 million people worldwide, with more than 300,000 deaths reported [2] . Due to the rapid rate of infection, many countries have adopted several precautionary measures to slow down the infection rate and not overwhelm its healthcare capacity. Some of these measures included the lock down of cities, imposing curfews, temporarily closing down non-essential businesses and issuing stay-at-home orders. The efficacy of these measures has yet to be determined, especially that in many cases, they have to be * Correspondence to: P.O. box 5969, Safat, 13060, Kuwait. weighed against the detrimental effect they impose on local and global economies [3] .

The prediction of the rate of infection of COVID-19 has become vital for decision and policy makers worldwide. It is important to estimate the rate as accurately as possible using reliable scientific techniques. A prediction of the number of infections would assist policy makers in a specific region to assess their current healthcare capacity and decide which measures need to be taken to curb and control the spread of COVID-19.

On February 24, 2020, the state of Kuwait reported the first confirmed case of COVID-19 in a citizen that had returned to the country from traveling abroad. Kuwait has an estimated total population of 4.26 million people [4] . Despite taking precautionary measures early on, such as suspending schools and non-essential businesses as well as imposing a partial curfew on March cases had grown to 5001 as of May 9, 2020 [5, 6] . The government had also imposed a 20 day total lockdown in early May [7] , however, on May 30, 2020, the total number of cases undergoing treatment grew to 16, 036, 194 of which were in critical condition [8] . In an effort to counter the economic effect the pandemic had on local businesses, Kuwait's government adopted a 5-phase plan of gradual preventive measures. Each phase would span a period of 3 weeks [9] . Every phase introduced a gradual recommencement of public and private services that had previously been halted. Table 1 shows the details pertaining to each phase of the preventive plan. At the end of each phase, the government would assess the efficacy of the preventive measures taken and re-evaluate the plan based on several criteria. These criteria included:

• Rate of COVID-19 transmission.

• Stability in number of infections observed.

• Decrease in number of patients in intensive care units.

• Decrease in number of patients hospitalized.

• Decrease in the number of positive results from COVID-19 swabs.

In epidemiology, it is common to use modeling techniques for predicting the spread of infectious diseases. There are many types of modeling techniques that researchers have applied in predicting epidemics. Modeling of infectious diseases is classified into 3 types: statistical methods for epidemic monitoring, mathematical and/or mechanistic state-space models, and empirical and/or machine learning models. Each type is then further classified into several categories. For example, statistical based methods include regression techniques, time-series auto-regressive techniques, spatial models, hidden Markov models, and statistical process control methods [10] . Mathematical modeling, that is formulated in time, can provide policy makers with crucial information that can be used to decide which precautionary measures should be enforced. It can also indicate the infection rate of the virus, the point at which it is most infectious, the severity of the infection, and the precautionary measures that need to be put in place to control the infections [11] . Statistical based mathematical modeling techniques are used by epidemiologists to predict and monitor disease outbreaks. Regression methods aim to identify an outbreak using time-series of periods where epidemics are non-existing. In these methods, the numbers of infected cases are monitored and an epidemic alert is be raised once a certain threshold is reached. However, their use is not recommended for the prediction of outbreaks in newly emerging epidemics [10] .

Auto-Regressive Integrated Moving Average (ARIMA), is a timeseries auto-regressive technique that calculates future short-term predictions from analyzing time-series of historical data. The ARIMA model was created by Box and Jenkins in the 1970s to describe changes of a time series in a mathematical approach [12] . ARIMA has been used in the past to predict several disease outbreaks such as Hemorrhagic Fever with Renal Syndrome (HFRS) [13] , Hand-Foot-Mouth Disease (HFMD) [14] , Hepatitis-B [15] , as well as the recently emerged COVID-19 virus [16] .

Since the emergence of COVID-19 in late 2019, researchers have been studying the pattern and rate of its infection using different mathematical modeling methodologies popularized in epidemiology research. In [17] , a prediction model was developed using Artificial Neural Networks (ANN) to estimate the growth of COVID-19 cases in the world based on geo-location and 2-week past data. This research compared the predicted numbers produced by their model with the actual values and found them to be closely matched. In [18] , a model called FPASSA-ANFIS is proposed. It improves the adaptive neuro-fuzzy inference system (ANFIS) using an enhanced flower pollination algorithm (FPA) by using the salp swarm algorithm (SSA). Their model's results showed good performance when comparing different accuracy measures against several other existing models. In [19] , a Neural Network model for COVID-19 spread prediction is proposed. The prediction model learns using NAdam training model and produced predictions for different countries and regions around the world. The proposed model achieved highly accurate results averaging 87.7% for most regions.

Some researchers focused on using compartmental models such as the classic susceptible, infected, recovered model (SIR) and the susceptible, exposed, infected, recovered model (SEIR) to model the spread of diseases such as Human Immunodeficiency Virus (HIV), Middle East Respiratory Syndrome (MERS), and Severe Acute Respiratory Syndrome (SARS) [20, 21] . In these models, the population would be divided and placed in one of the 3 or 4 compartments. An individual can be placed in one compartment at a certain time [22] . They can also be shifted between compartments based on their infectious status [23] . Several researchers have used these methods to model the spread of COVID-19 in different countries. In [11] , a dynamic SEIR model is proposed to estimate a daily reproduction number in China using time-dependent contact and diagnosis rates. In [24] , a stochastic transmission dynamic SEIR model was designed using multiple data sets for Wuhan as well as internationally exported cases from Wuhan. A dynamic extended SIR model was used in [25] to predict the trend of COVID-19 prevalence in Italy in comparison with Hunan, China. Another study focused on using a revised SIR model to estimate the number of infected individuals and mortality rate in Italy [26] . In [27] , SIR, SEIR, and HawkesN were applied using Poisson regression for the states of Indiana, New York and California. In [28] , an SIR model with delay in the context of the Atangana-Baleanu fractional derivative with Mittag-Leffler kernel was used. The author investigated and presented the asymptotic stability for the trivial and endemic equilibrium as well as the endemic equilibrium. In [29] , the author uses a stochastic model to predict COVID-19 by adding the terms of the stochastic perturbations to a deterministic model using Senegal as a case study. In [30] , a fractional SIRI model with time delay is considered. The research presents the global stability of the free fundamental point as well as the endemic point using the Lyapunov direct method. In [31] , an SIR model prediction was used to predict dates of peak infection rates for Kuwait from late February until 28 May 2020.

Several researchers have utilized the ARIMA model to forecast the spread of COVID-19 in many countries. In [32] , the researchers acquired Johns Hopkins epidemiological data to forecast COVID-19 incidence and prevalence. In [16] , a model was developed to estimate the prevalence point of COVID-19 in Spain, France, and Italy. Similarly, in [33] , a model using ARIMA forecasts the number of infections and recoveries from COVID-19 in Italy for May 2020. ARIMA model was also utilized in forecasting COVID-19 for India [34] and Iran [35] . In [36] , SutteARIMA, which is a method that combines ARIMA with -Sutte indicator, was used to forecast the date Spain would reach 200,000 cases. In [37] , ARIMA was used to forecast COVID-19 duration, infections, and deaths, as well as the effect on the Chinese economy in comparison with SARS virus. The researchers in [38] used ARIMA to predict prevalence, death, and recovery rates for 25 different countries.

Although the ARIMA model was utilized for predicting future COVID-19 infections; some researchers believe that it is not suitable for non-linear relationships especially in a complex and dynamic problem [39] . One limitation of the ARIMA model is its incapability of capturing non-linear patterns hidden within a time series [40] . However, the accuracy of an ARIMA based prediction should be investigated further, especially that it depends on the training set data and the parameters that are used in the model. Almost all ARIMA based model predictions use accuracy measures for selecting a best-fit model, however, the forecast values will not necessarily equate the actual values observed for the same time period. This can be due to several factors such as the different restrictions that are mandated by authorities to curb the spread and the degree to which the public adheres to these restrictions. The aim of this study is to evaluate how accurate the ARIMA model was in predicting COVID-19 infections and recoveries. The idea is to find a best-fit model early on and make the prediction for an extended period of time. Once the time period elapses, we then record the actual reported infections and recoveries, and compare them to the results produced by our best-fit prediction model. To the best of our knowledge, no ARIMA studies were found that compared an ARIMA model's prediction with actual values recorded for a significant time period. The novelty and contribution of this research addresses the aforementioned point by using Kuwait as a case study. Kuwait government's announcement and adoption of a 5-phase preventive plan provided us with the opportunity to use our ARIMA prediction model to forecast infection and recovery cases before the start of Phase I. Once the plan proceeded, we began collecting actual recorded data to be used for our comparison with the prediction.

The remainder of the paper is organized as follows. In section ARIMA based prediction, we describe, in general, how an ARIMA model works and then we provide justification for a best-fit ARIMA model which we used for our prediction. Section Results and Discussion describes the forecast produced by ARIMA for each of the 5 phases of Kuwait's preventive plan and presents a comparison between the forecast and the actual cases observed for the same time period. Finally, we draw our conclusion concerning the accuracy of the ARIMA based prediction in Conclusion.

We have acquired the training data set for the ARIMA model from the Kuwait government's official COVID-19 site [8] , as well as the COVID-19 site of WHO [41] . The data was pulled from the official sources using a Web-Application [42] and is updated daily. We used data from the date COVID-19 was first discovered in Kuwait on February 24, 2020, until May 30, 2020, the last day in the country's lock down.

Auto-Regressive Integrated Moving Average model, or ARIMA, is a time-series forecasting approach that is used in predicting the future value of a variable from its own past values. It uses auto-regression and moving average, and incorporates a differencing order to remove trend and/or seasonality. The model is expressed with the following equation:

where ′ represents the differenced series and is calculated by taking into account the lagged values of as well as the lagged errors. The ARIMA model contains 3 parameters ( , , ). Parameter in the ARIMA model represents the periods to lag for. For instance, a value of = 2 indicates that 2 previous time periods of the timeseries are used in the auto-regression part of the equation. Parameter represents the number of differencing transformations done to remove trend and/or seasonality therefore turning the time-series into a stationary one, i.e, making the mean and variance constant over time. This is an essential step to prepare the data for use in an ARIMA model. Parameter represents the lag of the error component of the ARIMA model. The error component is the part of the time-series that cannot be explained by trend or seasonality. Table 2 formally defines each of Table 2 Description of ARIMA parameters.

Parameter Description is the order of the auto regressive part of the model. is the order of differencing to remove trend/seasonality. = 0 means no differencing, = 1 means the series is differenced once and so on. is the order of the moving average part of the model; to smooth out random effects if any. these parameters. Incorporating the 3 parameters into Eq. (1) yields the following equation Eq. (2) [43] :

Selecting appropriate values for parameters and requires testing and optimization. To choose the values, one must inspect visual observations of the data to determine trend and/or seasonality. Visualizations in the form of auto-correlation function (ACF) and partial auto-correlation function (PACF) charts have proven to be useful in determining the values for parameters and .

Some rules may be followed in determining the values of and from ACF and PACF charts. For instance, if an ACF chart shows that the auto-correlation has drastically decreased at lag while the PACF chart shows a more sinusoidal decrease beyond lag , then should be zero and should be set to (the largest lag value). This is referred to as a MA( ), or a Moving Average model of order . On the other hand, if a PACF chart shows the partial auto-correlation to have drastically decreased at lag while the ACF chart shows a more sinusoidal decrease, then should be zero and should be set to . This is referred to as an AR( ), or an Auto Regression model of order [43, 44] . Relying on only manual visual inspections and following recommended guidelines may work in forecasting a small number of univariate time-series. However, in case large numbers need to be forecasted, manual visual inspection will not be possible. One solution for this problem is to incorporate the use of software packages that facilitate automated selection of ARIMA's parameters and chooses the best suited one. It should be noted that different software packages might yield different results because their performance relies on the implementation of the algorithm used.

In order to assess the prediction of an ARIMA model, we chose to implement a simulation in [45] . The language is an open source, platform-independent, programming language, specifically designated for statistical computing purposes and graphics. We have used the forecast package [43, 46] for predicting and ggplot2 [47] for graphing the results. The forecast package employs different automatic forecasting methods, such as exponential smoothing, Theta approach, cubic splines and it includes the ARIMA modeling approach. The package applies all the suitable models to the time-series and it features parameter optimization. The forecast package includes a function called auto.arima() which considers both seasonal and non-seasonal ARIMA models. The function returns the best fitting result based on values from Akaike Information Criterion (AIC), corrected AIC (AICc), or Bayesian Information Criterion (BIC). These values are used to compare the models in terms of quality, which is achieved by striking a balance between a model's fit and its complexity. For selecting parameter in non seasonal ARIMA ( , , ) model, the forecast package applies successive KPSS unit tests to the data to test for unit root, and when a significant result is found, the same method is applied to the differenced data. KPSS unit tests halt when encountering the first insignificant result [46] . Parameters p, q and c (if a constant c needs to be included), are chosen by minimizing AICc [48] .

The results produced by any ARIMA model are measured in terms of forecast accuracy. Observed data is divided into two sets; a training set to be used in estimating the parameters and a test set to be used in measuring the accuracy of the forecast. We can assume the observed data set is of size and the test set is of size , where < . Ideally, should be as large as the forecasted time-period. A forecast error is the difference between the observed value ( ) and its forecast (̂) [43] . The most commonly used scale-dependent measures of forecast accuracy are Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) defined in Eqs. (3) and (4), respectively. As for scale-independent measures, the most commonly used one is Mean Absolute Percentage Error (MAPE), defined in Eq. (5) [49] .

The process of choosing and optimizing ARIMA parameters starts out with plotting the data and looking for patterns and/or anomalies. If needed, transformations to stabilize the variances must be made. There are three steps that need to be taken when manually selecting a model. First, the data has to be differenced until it becomes stationary and this can be tested using unit-root tests. Next, ACF/PACF charts are produced for the differenced data, from which possible candidate models are selected. Once this is achieved, the selected models are compared in terms of their AICc values. The best-fit model would then be the model with the lowest AICc value. A comparison of the ACF plots of the residuals for close candidate models can be examined to determine which one resembles white noise [43] . Therefore, in our simulation we are considering RMSE, MAE, MAPE, AIC, AICc, and BIC values.

The auto.arima() function performs the three steps of the manual inspection process described above. When the function produces a result, the residuals of the resulting model are inspected. If they resemble white noise, the forecasts can then be calculated with optimal values for parameters p, d, and q. Furthermore, the auto.arima() function can be tweaked in order to find the best fitting model. The function, in its default form, might not take into account every possible model. Therefore, to ensure auto.arima() returns the optimal result, we can customize it to fit the data by setting function arguments such as approximation and step-wise to FALSE. Also, when the data is non-seasonal, the argument seasonal must also be set to FALSE [43] .

The auto.arima() function, when customized, seems to return a reasonable best-fit model. As a starting point, we chose to select the result model returned by auto.arima() on our confirmed cases data set. The result model returned by auto.arima() has the values = 1, = 2, and = 3, respectively. In order to verify whether or not this model is the best-fit, we applied the manual ARIMA function with various parameters for , , and following the steps described earlier.

It was noticed that some manual ARIMA models returned values for accuracy measures RMSE, MAE, and MAPE, that were slightly less (better) than those for the model selected by auto.arima() function. However, the values of manual ARIMA models for AICc were larger than the model returned by auto.arima() function. This suggests that some of the manual models over-fit the data. Furthermore, upon closer inspection of the ACF charts of residuals returned by these models, we have concluded that the model selected by auto.arima(1,2,3) was a best-fit model. The same process was performed for the recovered cases data set. The best-fit model returned by auto.arima() has the values = 2, = 2, and = 0. Similarly, results from our comparisons confirm that auto.arima(2,2,0) did provide a best-fit model while achieving a balance between accuracy and complexity of a model. In order to justify our chosen models of both confirmed and recovered, we show plots of the residuals in Fig. 1(a) and (b) . As observed, the sample auto-correlations for lags 1-20 do not exceed the significance bound. Moreover, the histogram of the residuals demonstrates a normal distribution.

To validate the selection of auto.arima(1,2,3) for confirmed cases, we examined its accuracy measures against ARIMA(2,2,1), which is one of the manual candidate models. ARIMA(2,2,1) returned less accurate results in all measures except for MAPE, where the value was 3.92. Fig. 2 shows a comparison of accuracy results between the two models. However, upon inspecting the ACF and PACF charts of the two models shown in Fig. 3(a) to (d), we observed that lags in ARIMA(2,2,1) were clearly exceeding the bounds. Furthermore, the AIC and AICc values of ARIMA(2,2,1) were larger than those of auto.arima (1,2,3) . Therefore, to avoid selecting a model that would over-fit our data set, we confirmed that auto.arima(1,2,3) was a more appropriate fit for our data. The same process was carried out for recovered cases. The accuracy results of auto.arima(2,2,0) are compared against ARIMA(1,2,1). Fig. 4 lists the accuracy measures being compared, and Fig. 5 (a) to (d) shows the comparison of ACF and PACF charts between the two models.

We measured the accuracy and quality of 7 candidate ARIMA models with different values for parameters , , and for confirmed cases as well as 7 others for recovered cases. The accuracy measures that were calculated are RMSE, MAE, MAPE, AIC, AICc, and BIC. Fig. 6 shows the values of the accuracy measures of all candidate models that were considered for confirmed cases. As for the candidate models for recovered cases; the accuracy measures are shown in Fig. 7 .

The accuracy measures for some manual ARIMA models returned values for RMSE, MAE, and MAPE, that are slightly more accurate than those for the model selected by the auto.arima() function, however, their AICc values were larger than the model returned by the auto.arima() function. This implies that some of the manual models might have over-fit the data. Furthermore, upon closer inspection of the ACF charts of residuals returned by manual ARIMA models, we have verified that auto.arima(1,2,3) for confirmed cases, was a bestfit model and is well suited for our prediction purposes. The same process was applied for recovered cases. Here, the best-fit model was auto.arima(2,2,0).

For the observed time period of February 24-May 30, 2020, the daily new, recovered, deceased and active cases as reported by the official Kuwaiti government, are shown in Fig. 8(a) to (d), respectively. We observed that there was a rapid and steady increase in the number of new cases, while the number of recovered cases showed a significant increase in the final week of May 2020. Active cases are the cases currently under treatment. The total number of active cases is obtained by deducting the total numbers of both deceased and recovered from the total number of new cases. The cumulative numbers of deceased cases, new cases, and recoveries are presented in Fig. 9 (a) to (c), respectively. The following sections compare the prediction results with the actual values reported for each phase of the government's preventive plan. Results of the forecast of auto.arima(1,2,3) model for confirmed cases in comparison with actual cases for the same time period, as well as the forecast of auto.arima(2,2,0) for recovered cases, are interpreted and discussed. Phase I Phase I of the government's preventive plan spans the period of May 31-June 20, 2020. Fig. 10(a) shows the prediction of cumulative confirmed cases for Phase I at 85% and 95% confidence intervals. In the best case scenario, by the end of Phase I, the forecast number of confirmed cases in the lower bound of 95% and 85% confidence intervals, would be between 39,818 and 41,663 cases. The actual number of cumulative confirmed cases for Phase I was reported to be 39,145. We observed a difference of 673 cases between the actual reported cases and the lower bound prediction of 39,818 at 95% confidence. The difference between the forecast point value by the end of Phase I and H. Alabdulrazzaq et al. the actual value reported is 7621 cases. Therefore, for the confirmed cases of Phase I, the prediction was higher than the value reported as shown in Fig. 11(a) . As for the forecast numbers of recovered cases, the result demonstrates recovery cases to be between 21,424 and 32,347 cases at 95% confidence by June 20, 2020. The actual value reported for recovery was 30,726 on the last day of Phase I. Fig. 11(b) indicates that auto.arima(2,2,0) had an accurate prediction of the number of recovery cases at the end of this phase.

Phase II spans the period from June 21-July 11, 2020. Fig. 12(a) shows the prediction of cumulative confirmed cases for Phase II at 85% and 95% confidence intervals. The forecast values for confirmed cases show that the higher bound predicted numbers at confidence intervals 85% and 95% are between 81,195 and 86,218. The forecast point value by the end of Phase II is estimated to be 67,300 confirmed cases. The lower bound number at 95% confidence intervals would be 48,382 cases, constituting an optimum situation. The actual number of cases reported at the end of this phase was 54,058. This value falls within the prediction bounds given by auto. arima(1,2,3) with a difference of only 5674 cases as shown in Fig. 13(a) . With respect to the forecast numbers of recovered cases, the result shows the best case scenario could return a recovery rate of 58,873 or 54,837 cases at 95% and 85% confidence respectively, by the end of Phase II. The forecast point for recoveries by July 11, 2020, is 43,673 cases and the actual number of recovery cases on that day is 43,960. In this phase, auto.arima(2,2,0) also produced an accurate prediction for the number of recovery cases. Fig. 13(b) shows the actual number of recovered cases in Phase II with respect to the prediction bounds.

Phase III spans the period of July 12-August 1, 2020. Fig. 14(a) shows the prediction of cumulative confirmed cases for Phase III. The forecast values for confirmed cases show that the higher bound numbers for 85% and 95% confidence intervals are 113,031 and 122,141 whereas the lower bound numbers are 62,638 and 53,529 cases. The forecast point value by the end of Phase III is estimated to be 87,835 confirmed cases. The actual value reported for confirmed cases in Phase III is 67,448. From Fig. 15(a) , we observe that the actual value is within the upper and lower bounds of auto.arima (1, 2, 3) prediction. The forecast values for recovered cases show that the higher bound estimates of recovered cases at 85% and 95% confidence are 80,859 and 88,233 cases, respectively. The worst possible outcome for recovery cases would be 32,687 at the lower bound of 95% confidence. The forecast point for recoveries by August 1, 2020 is 60,460. The actual reported number of recovered cases on the last day of Phase III is 58,524 cases. Again, the actual value reported is between auto.arima (2, 2, 0) prediction bounds at 95% confidence interval as depicted in Fig. 15(b) . 

Phase IV spans the period from August 2-August 22, 2020. Fig. 16(a) shows the prediction of cumulative confirmed cases for Phase IV at 85% and 95% confidence intervals. In the best case scenario, by the end of Phase IV, the forecast number of confirmed cases in the lower bound of 95% and 85% confidence intervals, would be between 55,896 to 69,829 cases. The highest forecast number could reach 160,842 confirmed cases. The forecast point value by the end of Phase IV is estimated to be 108,369 confirmed cases. Fig. 17(a) H. Alabdulrazzaq et al. shows the actual value reported on August 22, 2020 for confirmed cases is 79,957. This means that the actual value falls well within the lower and upper prediction bounds for both 95% and 85% confidence intervals. By the end of Phase IV, the forecast number of recovered cases in the lower bound of 95% and 85% confidence intervals, is between 34,604 and 45,927 cases. The highest forecast number could reach 119,891 recovered cases. The forecast point value by the end of Phase IV is estimated at 77,247 confirmed cases. The actual number of recovered cases is 71,769, as is illustrated in Fig. 17(b) . This shows that auto.arima(2,2,0) has produced a relatively accurate prediction.

The final phase of the preventive plan, Phase V, covers the period from August 23 September 12, 2020. At the time of writing this research, Kuwait's government opted to remain in Phase V. It should be noted that Phase V remains currently in effect. Fig. 18(a) shows the prediction of cumulative confirmed cases for the final phase at 85% and 95% confidence intervals. In a less than ideal circumstance, the higher bound number at 95% confidence intervals would be 201,949

cases. The lower bound numbers at confidence intervals 95% and 85% are between 55,859 and 75,255. The forecast point value by the end of Phase V is estimated to be 128,904 confirmed cases. The actual number of confirmed cases on September 12, 2020 is 94,211 cases as is shown in Fig. 19(a) . With respect to the forecast numbers of recovered cases, the result shows the best case scenario could return a recovery rate of 153,533 or 137,734 cases at 95% and 85% confidence respectively, by the end of Phase V. The worst possible outcome for recovery cases would be 34,537 at the lower bound of 95% confidence. The forecast point for recoveries by September 12, 2020 is 94,035 cases. The actual number of recovered cases on September 12, 2020 is 84,403, which is within the boundaries predicted as is depicted in Fig. 19(b) . Once more, the actual value recorded for recoveries is within the upper and lower bounds given by auto. arima(2,2,0) . At the time of writing this research, Kuwait government was still in Phase V of the plan. As such, we used our ARIMA models to forecast confirmed and recovered cases beyond September 12, 2020 and until December 4, 2020. Figs. 20 (a) and (b) show the confirmed and recovered predictions as well as the actual reported cases for this time period. Figs. 21 (a) and (b) illustrate the comparison between the actual values reported, the forecast point value, and the boundaries at 95% confidence interval for the entire duration of the prediction up to December 4, 2020. In comparison with the actual values reported for this period, the forecast values of our selected ARIMA models were within the upper and lower bounds at 95% confidence intervals for most phases. The Pearson's correlation coefficient was equal to 0.996 for the forecast data in correlation with the actual data recorded. This signifies a very high correlation between the two sets.

In Phase I, the actual values reported fell under the lower bound of the prediction. This can be attributed to the fluctuating number of COVID-19 tests conducted on the population during that time period. The number of administered COVID-19 tests was not consistent throughout each phase and was particularly scarce in the beginning. If the number of tests conducted was higher, then the actual values reported could have been closer to the forecast point.

For the most part, the population has adhered to the preventive measures mandated by the government. Individuals who violated these preventative measures paid a high monetary penalty. The government also enforced a mandatory self-quarantine period for individuals entering the country and was surveilling their location and progress through a government sponsored smartphone application. The population's adherence to the measures taken by the government had a significant effect in reducing the actual number of confirmed cases as was predicted by our ARIMA model. Nonetheless, our ARIMA model's prediction bounds at 95% confidence interval were accurate; the actual values recorded for the time period observed fell between the upper and lower bounds.

Pandemics are unpredictable by nature. This imposes a limitation on the use of ARIMA which belongs to a class of linear models. ARIMA cannot capture non-linear patterns hidden within a time-series [40] . One major threat to the validity of our results is the lack of consistency in the number of COVID-19 tests conducted daily by the health authorities in Kuwait. This number was not made publicly available at the beginning. If the number of daily tests was consistent, the actual number of confirmed cases could be significantly higher than what was announced. This, in turn, would affect the validity of our results obtained from the observed time period data set. Another factor that may have impacted the accuracy when comparing the predictions with the actual reported cases, is that initially, the COVID-19 virus was misdiagnosed and there were many uncertainties about the different symptoms associated with it. This could have caused an unknown number of false positives or false negatives in diagnosed cases.

Although ARIMA was used in epidemiology predictions for past diseases as well as the current pandemic of COVID-19, it has not been relied on heavily because it is regarded as unsuitable to use in complex and dynamic situations. This research aimed to test ARIMA model's accuracy in predicting COVID-19 confirmed and recovered cases using Kuwait as a case study. Kuwait government announced a 5-phase gradual preventive plan in late May 2020. After this announcement, we used the ARIMA model to predict confirmed and recovered cases for each one of the 5 phases. With careful considerations for the model's parameters, both manual and automated, we have found a best-fit model for our training data set obtained from official Kuwait government sources. The model we selected was examined in terms of accuracy measures such as RMSE, MAE, and MAPE, AIC, AICc, and BIC, as well as ACF and PACF charts. We provided validation for our selected model and used it to predict confirmed and recovered cases for Kuwait's 5-phase preventive plan. We then observed actual daily reported confirmed and recovered cases for the time period spanning the 5-phase plan to compare it with our prediction. Even with constant revisions and changes to the plan, our ARIMA model was reasonably accurate in its prediction. The actual values of confirmed and recovered cases were well within the upper and lower bounds of the prediction at 95% confidence interval for all phases except Phase I. We attribute this disparity to the inconsistent COVID-19 testing undertaken during the first phase of Kuwait's preventive plan. Pearson's correlation for the forecast and actual data sets was calculated to be 0.996, which indicates a very high correlation between the two data sets. Despite the dynamic conditions governing this disease and the continuous revisions of measures taken by Kuwait's government, the accuracy of the ARIMA model prediction is appropriate and satisfactory.

An interactive web-based dashboard to track COVID-19 in real time

Coronavirus disease (COVID-19) pandemic

Economic effects of coronavirus outbreak (COVID-19) on the world economy

World population by country

Schools suspended from 1st march to 12th march. 2020

Kuwait imposes partial curfew amid 12 new coronavirus cases

Kuwait imposes 20 day total curfew from may 10 to curb coronavirus

COVID-19 updates

Kuwait to move to 12 hour curfew from sunday

Mathematical modeling of infectious disease dynamics

An updated estimation of the risk of transmission of the novel coronavirus (2019-ncov)

Time series analysis: Forecasting and control

A new seasonal difference space-time autoregressive integrated moving average (SD-STARIMA) model and spatiotemporal trend prediction analysis for hemorrhagic fever with renal syndrome (HFRS)

Predicting the incidence of hand, foot and mouth disease in Sichuan province, China using the ARIMA model

Comparison of ARIMA and GM (1, 1) models for prediction of hepatitis B in China

Estimation of COVID-19 prevalence in Italy

Real-time neural network based predictor for cov19 virus spread

Optimization method for forecasting confirmed cases of COVID-19 in China

Neural network powered COVID-19 spread forecasting model

Epidemic models of contact tracing: Systematic review of transmission studies of severe acute respiratory syndrome and middle east respiratory syndrome

A study on the efficiency of the estimation models of COVID-19

Mathematical models in epidemiology

An introduction to compartmental modeling for the budding infectious disease modeler

Early dynamics of transmission and control of COVID-19: A mathematical modelling study

Extended sir prediction of the epidemics trend of COVID-19 in Italy and compared with hunan

A modified sir model for the COVID-19 contagion in Italy

The challenges of modeling and forecasting the spread of COVID-19

SIR epidemic model with Mittag-Leffler fractional derivative

Analysis of the stochastic model for predicting the novel coronavirus disease

Fractional SIRI model with delay in context of the generalized Liouville-Caputo fractional derivative

Building a sensible sir estimation model for COVID-19 outspread in Kuwait

Application of the ARIMA model on the covid-2019 epidemic dataset

COVID-19 disease outbreak forecasting of registered and recovered cases after sixty day lockdown in Italy: A data driven model approach

COVID-19): ARIMA based time-series analysis to forecast near future

Exponentially increasing trend of infected patients with COVID-19 in Iran: A comparison of neural network and ARIMA forecasting models

The date predicted 200.000 cases of COVID-19 in Spain using SutteARIMA

Risk prediction and assessment: Duration, infections, and death toll of the COVID-19 and its impact on China's economy

Modeling and short-term forecasts of indicators for COVID-19 outbreak in 25 countries at the end of march

Covid-19: A comparison of time series methods to forecast percentage of active cases per population

Carbon price forecasting with a novel hybrid ARIMA and least squares support vector machines methodology

Covid-19 dashboard by country

Kuwait corona virus platform:covid19

Forecasting: Principles and practice. OTexts

Notes on non-seasonal models

R: A language and environment for statistical computing

Automatic time series forecasting: The forecast package for r

Elegant graphics for data analysis

Automatic forecasting algorithms

Performance of state space and ARIMA models for consumer retail sales forecasting