COVID-19 is the latest infectious virus that has become a global pandemic and brought the 1 global economies to their knees. Precise analysis and forecast of the disease spread can help with 2 resource planning and create strategies to slow down the progress of this deadly virus. This paper 3 explores a variety of machine learning models, from heuristic statistical techniques to advanced deep 4 learning methods, to forecast the COVID-19 dynamic. To measure the daily spread of COVID-19, 5 we opt for two target variables: the number of daily positive cases and the number of daily deaths. 6 Although the chance of irregularities and reporting lags is high, it is more sensitive to short-term 7 time series forecasting. These two variables look for stable and reliable estimates for COVID-19 8 spread. The peculiarity of the data is that it is time series but without one complete period, thereby 9 preventing us from directly using established forecasting methods. Thus, our analysis uses some 10 non-time series methods by including time factors and a few time series methods with the inclusion 11 of exogenous variables by tailoring the data into the appropriate format. We aim to ﬁnd an optimal 12 model for each family of models where possible. To illustrate the results, India has been chosen for 13 the case study, as this country presently recorded the fastest pace of COVID-19 spread in the ﬁrst six 14 months of the pandemic. A comparative study has been included with different evaluation metrics. 15 The metrics such as Mean absolute error (MAE), Mean squared error (MSE), Median squared error 16 (MEME), and Mean squared log error (MSLE) has been used for evaluating the spread of COVID-19. 17 We have compared methods such as Liner Regression, Elastic net regularization, Random-forest 18 regressor, XGBoost regressor, Simple exponential smoothing, and so on. Among these methods, the 19 Random-forest regressor shows the highest MAE (11351.8833), MSE(11827.2160), MEME(9998.6333), 20 and MSLE(0.0220) values than the other state-of-the-art methods. Our study indicates that more 21 complex models may not be more reliable compared to simpler ones for forecast COVID-19 spread. 22 We have used python to analyze our results. 23 Keywords: COVID-19; Machine Learning; Linear Regression; Random Forest; Time Series; XGBoost 24 1. Introduction 25 It has been nearly 15 months [ 5] since the world has noticed the most devastating pandemic of 26 the21stcentury – COVID-19. COVID-19 has cost more deaths and misery all over the world than 27 anything else in the past century [ 2]. Major lifestyle changes have been observed as most of the world 28 was shut down in the ﬁrst year of the pandemic. For instance, schools were closed and people were 29 stuck in their homes. It led to devastating changes in human civilization in which social interaction 30 was forbidden, the global economy was on its knees, and people lost jobs in nearly every sector. In 31 addition, ﬁrst responders, medical professionals, and critical workers have been on their toes and 32 Submitted to Mach. Learn. Knowl. Extr. , pages 1 – 20 www.mdpi.com/journal/make Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 2 of 20 continuously active in the war against COVID-19. The latter is true even in the present day when 33 vaccines are readily available in many parts of the world, most economies have reopened, and many 34 lives have returned mostly to normal. Despite currently being in a less intense phase of the pandemic, 35 however, COVID-19 is still adversely affecting the world, through supply chain shortages and delays, 36 worker shortages, and repeated stresses on healthcare systems due to waves caused by new variants 37 of the virus. In addition, many individuals infected by COVID-19 suffer long-term effects of the virus 38 long after the infection is over. For instance, this virus has caused many cardiovascular diseases [7]. 39 The motivation of the other side of this research is the speed prediction in the growth of 40 the COVID-19 virus. Prediction for the daily spread of COVID-19 in future days can be helpful 41 for government and medical staff to be prepared for current and future waves of the pandemic. 42 Furthermore, predictions of new conﬁrmed cases and new deaths can help predict the daily spread of 43 COVID-19. Because these two variables each have their own pros and cons, no one variable can be 44 chosen over the other. The number of reported positive cases can be biased, as it is highly dependent 45 on the number of tests done, which is highly dependent on the number of test centers available in the 46 geographical unit. The cause of death, however, can be something other than COVID-19, i.e. dying of 47 COVID versus DYING with COVID. In addition, the mortality rate of COVID-19 is not static, as there 48 is often a lag in reporting both numbers. Despite all these pitfalls, these two variables are still the most 49 ideal for our analysis, as there is a no better indicator to capture COVID-19 exposure [ 10–15,31]. The 50 combined use of both metrics is a novel focus in this paper, as the prediction of COVID deaths has not 51 been addressed in prior works. 52 Section 2 provides a literature review of similar work and their results. In Section 3, we discuss the 53 data used to carry out our analysis. Sections 4, 5, and 6 are devoted to predication for daily conﬁrmed 54 cases. In Sections 4 and 5, we explore some predictive and time-series forecasting models respectively. 55 We compare all models built in these sections in Section 6.1. Brief details on the results for the number 56 of daily deaths are given in Section 6.2, and the ﬁnal conclusions are made in Section 7. 57 2. Related work 58 A prediction model is used to analyze future conditions based on the data available. Many 59 predictive modeling methods use statistics to predict events [ 17,18]. Forecasting always plays an 60 important role in assisting the predictive outcomes of many models to analyze the accuracy of the 61 prediction framework. This is estimated across different study populations, ecosystems, and locations 62 for further improvement of the model [ 19]. Yang et al. [ 20] proposed a new method to identify the 63 forecast of the COVID-19 virus using the SEIR and AI model and showed a good quality assessment of 64 95%. Liang et al., [ 21] used the Statistical software: LASSO, a logistic regression model to forecast the 65 risk of critical illness of the patient who is affected by COVID-19. An accuracy of 88% was achieved by 66 this method. Yan et al., [ 22] used the Machine learning tool: XG Boost to relieve the clinical burden and 67 reduce the mortality rate of the people who are affected by COVID-19. Another interesting method 68 proposed by Gong et al., [ 23] used statistical analysis for predicting the forecast of COVID-19. However, 69 the accuracy achieved was not higher than the other methods. Chatterjee et al. [ 24] proposed a new 70 method namely SEIR to predict the presence of COVID-19 in the people. Tomar and Gupta [ 25] used 71 the LSTM method for prediction purposes. Another method that used the LSTM was proposed by 72 Chimmula & Zhang [ 26]. IHME COVID-19 Health Service Utilization Forecasting Team & Murray [ 27] 73 analyzed the presence of COVID -19 using the statistical model. 74 Many machine learning methods were used to predict the forecast of the spread of the COVID-19 75 virus. Pandey et al. [ 28] used SEIR and regression models for COVID-19 outbreak predictions. A 76 machine learning forecasting model achieved high accuracy in predicting the outbreak [ 29]. Deep 77 learning models for the prediction and analysis of COVID-19-positive cases were proposed by Ghosal 78 et al [ 30]. Another yet interesting method that used LSTM and RNN for predicting and analyzing the 79 COVID-19 positive cases proposed by Arora et al., [31] showed a better performance. 80 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 3 of 20 The prediction for the number of conﬁrmed cases has been carried out by many researchers. [ 1] 81 have proposed a mathematical model to predict the dynamic of COVID-19 for India. In the initial 82 days of COVID-19 spreading, [ 3] used ARIMA, a wavelet-based forecasting model, and a hybrid 83 implementation of both models. A deep learning model, LSTM, has been explored by [ 5,6,25] to 84 predict the number of conﬁrmed cases. To the best of our knowledge, no one has tried to predict 85 the number of COVID-19 deaths so far. For predicting the number of daily deaths, number of daily 86 positive cases, number of daily recovered cases, and cumulative number of conﬁrmed cases, [ 10] 87 used a support vector machine model. Similarly, [ 11] was able to predict cumulative daily counts of 88 conﬁrmed cases, deaths, and recoveries. A few other surveys have newly been divulged, but they did 89 not cover much observation of many machine learning and deep learning uses. 90 3. Exploring the data and feature engineering 91 The main objective of this paper is to compare different models for forecasting COVID-19 spread. 92 Thus, we require data from a geographical unit for the case study. We chose the following source for 93 data - https://github.com/GoogleCloudPlatform/covid-19-open-data, as it is available for various 94 countries and at different geographic levels. This source has multiple datasets such as epidemiology 95 (COVID-19 data statistics), demographics, economy, weather, health, mobility, government response 96 data, etc. We use a compiled version of all these datasets in our analysis. 97 First, we chose the top three most infected countries over the ﬁrst six months of the pandemic - 98 the US, India, and Brazil - for the analysis. More than 40 million positive cases were reported in each of 99 these countries over such a time period. We have multiple choices to describe the COVID spread such 100 as the number of reported positive cases and the number of deaths. Both numbers have drawbacks, as 101 detailed in the previous section, so to avoid misleading results, we consider both features together. 102 The raw data is available from January 1, 2020, to the present day. For most of the countries, 103 COVID-19 data was not updated for the initial days. Therefore, we consider data from February 15, 104 2020, as the starting point for our analysis. While there are presently over two years of data with which 105 to work thereafter, we will only focus on data going up to September 1, 2020. This is due to the focus 106 of this paper being on comparing machine learning algorithms in overall effectiveness in predicting 107 the spread of COVID-19 rather than predicting present spread levels. Since we are trying to compare 108 algorithms with the use of ground truth data, it is ideal to narrow the focus of the overall timeline to 109 the ﬁrst few months of the pandemic rather than every stage and wave encountered thus far. 110 Mar 2020Apr May Jun Jul Aug Sep date020000400006000080000100000Number of daily new casesNumber of daily new cases Brazil India US (a)positive cases Mar 2020Apr May Jun Jul Aug Sep date0100020003000400050006000Number of daily deathsNumber of daily deaths Brazil India US (b)deaths Figure 1. Time series plots for COVID-19 spread data Figure ??depicts plots of the number of positive cases and the number of deaths reported daily for 111 all three countries. Over the ﬁrst few months of the pandemic, the US was the most infected country, 112 although the trend eventually gravitated towards Brazil. Towards the end of the timeline, there was 113 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 4 of 20 exponential growth in India. Since India seemed to be the worst country in terms of daily COVID-19 114 cases at the end of the timeline, we have chosen India for the case study to compare various prediction 115 techniques. 116 After ﬁltering the data for India at the country level, we pre-processed the data for modeling. 117 For the remainder of this paper, the word data stands for Indian data at the country level. There are 118 a few drawbacks and limitations to the raw data. For instance, the population and related variables 119 have static data over time across all the rows. Since we have exponential and substantial growth in 120 the death poll, we should not consider the given demographic data as is. Thus, we distribute total 121 deceased counts uniformly across gender and ten buckets of age data. This helps us to update the 122 daily population and related variables in a meaningful way. 123 One of the most signiﬁcant features we have available is mobility data. This tells us about the 124 change in footfall and visitation patterns of consumers at different locations such as stores, parks, 125 restaurants, cafes, workplaces, and homes. 126 Since we should evaluate the model performance on unseen data, we implement a supervised 127 learning process by dividing our data into training and test data. For each class of models, we adopted 128 the same train and test data for the sake of comparison of forecast values. We use the time frame 129 leading up to September 1, 2020, as the training partition and the time frame thereafter as the test 130 partition. 131 In the following sections, we showcase the usage of selected statistical and machine-learning 132 models to predict the spread of COVID-19. 133 4. Non-time-series predictive models to predict number of daily positive cases 134 In this section, we explore and implement some classes of predictive models to forecast the 135 number of daily positive cases. To impose the time factor, we construct a new variable called delay , 136 which is the difference in days from the oldest date in the data. This variable is included in the list of 137 predictors for all the models covered in this section. 138 In the following subsections, we aim to get the optimal model from each class. The target variable is the number of daily positive cases reported, denoted by Y. The value of Ymust be non-negative, so in order to avoid predictions by models from being negative, we implemented the transformation Y !log(1+Y) (1) for the target variable. 139 Most models in the section have one or more hyperparameters, which when properly tuned can 140 provide us with an optimal model. Thus, we use a model-tuning approach to ﬁnd the best values 141 of hyper-parameters. We deﬁne the search space for hyperparameters with scoring criteria as mean 142 squared error. Once the model and tuning parameter values have been deﬁned, we need to specify 143 the type of resampling. We opt for repeated k-fold cross-validation with 5 folds, repeated 10 times to 144 get the best values of hyper-parameters. The model corresponding to these hyper-parameters is the 145 optimal model, due to having the smallest amount of mean squared error. Each of these models is 146 implemented in Python using various libraries detailed in the following subsections. 147 4.1. Linear regression 148 Linear regression [ 32] can be used to ﬁnd the linear relationship between a target variable and 149 one or more independent variables. This model is a basic regression model for comparison and can be 150 treated as a baseline model. This model is created using the OLS (ordinary least squares) library in the 151 statsmodels Python library. 152 The standard regression model is represented in the equation.4.1 : yt=x0 tbut(t=1, 2, .... T) (2) Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 5 of 20 Where ytrepresents the t0thobservation of the dependent and response variable. X1 is the column 153 vector of the observation K which is the independent and regression variable. The index t is the time 154 series data. bis the Kx1 vector to be estimated and utis the stochastic term. 155 The ﬁrst regression model is built by using all predictors. The importance of predictors is given in 156 Figure 2. 157 rainfall dew_point relative_humidity school_closing workplace_closing cancel_public_events restrictions_on_gatherings public_transport_closing stay_at_home_requirements restrictions_on_internal_movement international_travel_controls income_support fiscal_measures international_support public_information_campaigns testing_policy emergency_investment_in_healthcare investment_in_vaccines stringency_index mobility_retail_and_recreation mobility_grocery_and_pharmacy mobility_parks mobility_transit_stations mobility_workplaces mobility_residential delay10000 5000 050001000015000 Figure 2. Coefﬁcients of regression equations with 95% conﬁdence interval Some predictors are found to have large p-values, and their corresponding correlation coefﬁcients 158 are nearly zero. Such predictors are not signiﬁcant. We choose the level of signiﬁcance a=0.05 and 159 skip the predictors with p-values greater than a. Table 1 shows the values of R2and adjusted R2160 for both regression models: one with all predictors and one with only signiﬁcant predictors. Both 161 models have fairly high values for R2and adjusted R2, but both values seemed to worsen when we 162 skip insigniﬁcant predictors. 163 R2Adjusted R2 Model with all predictors 0.989 0.987 Model with signiﬁcant predictors 0.986 0.985 Table 1. R2and adjusted R2values for different linear regression models 164 Figure 3 compares the results of both models against the actual values. To our surprise, the model with 165 all predictors outperformed the one with only signiﬁcant predictors from every angle, since the red 166 line is closer to the black one (actual values) than the blue for all given date ranges. Thus, to compare 167 the linear regression model with other classes of models, we use only the model with all predictors 168 onward. 169 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 6 of 20 01 Sep 202006 02 03 04 05 date740007600078000800008200084000860008800090000Number of daily confirmed casesLinear Regression Actual Predicted with all predictors Predicted with significant predictors Figure 3. Comparison of actual daily case counts with predicted counts from two regression models: one with all predictors and one with signiﬁcant predictors 4.2. Elastic net regularization 170 To overcome model complexity and overﬁtting that can occur in simple linear regression, two 171 other penalized regression models - Ridge ( L2regularization) and Lasso regression ( L1regularization) 172 - have been widely used. The overﬁtting occurs due to the large model parameters. The elastic net 173 regularization is used as same as the ridge or Lasso. If the mixing parameter is zero, then we can use 174 ridge regression. If the mixing parameter is one, then we can use the lasso regression [33]. 175 In the section, using the linear model package of Python’s scikitlearn library, we ﬁt a model 176 known as elastic net regularization, which is the generalization of the two penalized regression models. 177 This class of models has two hyper-parameters: 178 a: mixing parameter, which controls the type of regression 179 l: shrinkage parameter which is the amount of the shrinkage. 180 The search space is chosen as 181 a2 f 0.1, 0.2, . . . , 1g, l2 f 105, 104, . . . , 101, 1, 101, 102g. After hyperparameter tuning, the optimal values turned out to be a=0.2and l=0.1. Thus, we 182 consider this model for this class of models to compare in the next section. 183 4.3. Random forest regressor 184 Random forest [ 34] is a supervised machine learning algorithm used for classiﬁcation and 185 regression. This is a bagging ( bootstrap aggregat ing) ensemble learning method that combines 186 (i.e., aggregates) the predictions from multiple decision tree algorithms with varying bootstrapped 187 subsets of data to make more accurate predictions than any individual one. To ensure that the model 188 does not rely on any individual predictor, the number of predictors used for a split is controlled by 189 hyperparameters speciﬁc to the random forest, including: 190 n_estimators = number of trees in the forest, 191 max_features = number of maximum features to consider at every split, 192 max_depth = maximum number of levels in the tree, 193 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 7 of 20 min_samples_split = minimum number of samples required to split a node, 194 min_samples_leaf = minimum number of samples required at each leaf node, and 195 bootstrap = method of selecting samples for training each tree. 196 To ﬁnd the best hyperparameter value, we choose the following parameter space: 197 n_estimators2 f 50, 100, 200, 500, 1000 g max_features 2 f0auto0,0sqrt0g max_depth2 f 5, 20, 50, 100g min_samples_split 2 f 2, 5, 10g min_samples_leaf 2 f 1, 2, 4g. After tuning, the optimal random forest regressor uses the following optimal values: 198 n_estimators =200 max_features =0auto0 max_depth =50 min_samples_split =2 min_samples_leaf =5. We consider this model from this class of models for comparison in Section 6. 199 4.4. XGBoost regressor 200 The XGBoost [ 35] is a widely used supervised machine learning model that is an implementation 201 of the gradient boosting decision tree algorithm. The validity of this statement can be inferred by 202 knowing about its (XGBoost) objective function and base learners. The objective function contains a loss 203 function and a regularization term. It tells about the difference between actual values and predicted 204 values, i.e how far the model results are from the real values. The most common loss function 205 in XGBoost for regression problems is reg:linear, and that for binary classiﬁcation is reg:logistics. 206 Ensemble learning involves training and combining individual models (known as base learners) to 207 get a single prediction, and XGBoost is one of the ensemble learning methods. XGBoost expects to 208 have the base learners which are uniformly bad at the remainder so that when all the predictions are 209 combined, bad predictions cancels out and better one sums up to form ﬁnal good predictions. This 210 algorithm has the following hyperparameters: 211 n_estimators = number of gradients boosted trees, 212 objective = a learning objective function corresponding to the learning task, 213 learning_rate = step size shrinkage for tree booster, 214 max_depth = maximum tree depth for base learners, 215 min_child_weight = minimum sum of instance weight (hessian) needed in a child, 216 min_samples_leaf = minimum number of samples required at each leaf node, and 217 bootstrap = method of selecting samples for training each tree. 218 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 8 of 20 To ﬁnd the best value of hyper-parameters, we choose the following search space: 219 n_estimators2 f 50, 100, 200, 500, 1000 g objective2 f0reg:squarederror0,0reg:squaredlogerror0g learning_rate 2 f 0.2, 0.5, 0.8g max_depth2 f 5, 20, 50, 100g min_child_weight 2 f 3, 4, 5g silent2 f 0, 1g subsample2 f 0.2, 0.7g colsample _bytree 2 f 0.2, 0.7g. The optimal XGBoost regressor corresponds to the values of following hyper-parameters: 220 n_estimators =50 objective =0reg:squarederror0 learning_rate =0.5 max_depth =5 min_child_weight =5 silent =0 subsample =0.7 colsample _bytree =0.7. We consider this model for comparison in Section 6 using the xgboost Python library. 221 4.5. Recurrent neural network (RNN) 222 A neural network is a predictive model that uses layers of neurons to map inputs to outputs 223 using the multiplication of weights and neuron values followed in some cases by activation functions. 224 The weights are optimized using backpropagation. The latter is used to add non-linearity to a model, 225 thereby serving as a stark contrast to linear regression, in which inputs and outputs can only correlate 226 linearly. 227 A typical neural network has input, output, and hidden layers. The former two are 228 self-explanatory, while hidden layers connect the two. A recurrent neural network is a variation 229 of this that involves time. While input, hidden, and output layers can connect to one another like 230 before, an RNN can also connect between hidden layers of adjacent time steps, thereby allowing neural 231 network modeling of simple time-series problems. However, in our study, RNNs [ 36] are fairly limited 232 in that a particular point in time only has a connection to adjacent time steps, and thus the information 233 for one particular data can only be directly inﬂuenced by the most immediate previous day. 234 We implement RNN, as well as the following two methods, using the keras API of the Tenspr f low 235 deep learning framework. 236 4.6. Long short-term memory network (LSTM) 237 The long short-term memory (LSTM) [37] network is an advanced deep learning method based 238 on RNN to forecast time-series data. Instead of neurons, LSTM networks have memory blocks that are 239 connected through layers. A block has components that make it smarter than a classical neuron and 240 a memory for recent sequences. A block contains gates that manage the block’s state and output. A 241 block operates upon an input sequence and each gate within a block uses sigmoid activation units 242 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 9 of 20 to control whether they are triggered or not, making the change of state and addition of information 243 ﬂowing through the block conditional. 244 Using LSTM, we can frame this problem as the following regression problem: what will be the 245 number of positive cases tomorrow given the number of positive cases today and previous k1days? 246 The parameter kis known as look-back, which decides how many previous time steps we want to 247 include. For simplicity, we choose k=1. Therefore, we must convert our univariate data into bivariate, 248 where the ﬁrst variable indicates the number of the present day’s positive cases and the second variable 249 stands is the number of positive cases predicted on the next day. Since this method is sensitive to the 250 scale of data, we, therefore, normalize the data to lie between 0 to 1. To build this model, we use the 251 default settings. 252 4.7. Gated recurrent unit (GRU) 253 A 254 5. Time-series forecasting method to forecast number of daily positive cases 255 In this section, we explore some time series methods to predict daily cases. These models are 256 forecasting methods that are completely based on the demand history of the item which has been 257 forecasted. These methods work by capturing the patterns in the historical data and extending the 258 application into the future. They are appropriate when you can assume a reasonable amount of 259 continuity between the past and the future. A common approach to model time series is to treat the 260 current time step Ytas a variable dependent on previous time steps Ytk. 261 5.1. Exponential smoothing 262 Exponential smoothing [ 38] is a powerful time series forecasting method for univariate data. 263 There are many different kinds of exponential smoothing methods, such as: 264 Simple exponential smoothing, 265 Double exponential smoothing (Holt method), 266 Triple exponential Smoothing (Holt-Winters method). 267 These methods are implemented using the tsa(Time Series Analysis) packages of the statsmodels 268 Python library. Each of these methods is explored further in the following subsections. 269 5.1.1. Simple exponential smoothing 270 As the name suggests, simple exponential smoothing is the simplest method. It is widely used 271 when our univariate time series data has no clear trend or no seasonal pattern. This method forecasts 272 using weighted averages with the largest weights associated with the most recent observations and 273 the smallest weights to the oldest observations. The weights decrease rate is controlled by a parameter 274 known as a smoothing parameter, denoted by a. The value of alies between 0 to 1, where a larger 275 value requires the model to pay close attention to the most recent past observations. The extreme cases 276 are: 277 a=0:Becomes an average since all weights are equal and the next predicted value is equal to 278 the average of historical data, 279 a=1:Becomes a naive method since a weight’s most recent observation is one and all others 280 are zero. Thus, the next predicted value is the same as the recent observation. 281 5.1.2. Double exponential smoothing (Holt method) 282 This is an extension of simple exponential smoothing. Double exponential smoothing was 283 proposed by Holt in 1957. We use simple exponential smoothing when there is no clear trend or 284 seasonality, but if we know the trend of data, we can use this extended method. Holt’s method 285 involves the following two parameters: 286 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 10 of 20 a=smoothing parameter, 287 b=trend smoothing parameter. 288 Both parameters take values between 0 to 1. There is also an option to choose a trend type. It can be 289 either additive or multiplicative, indicating a linear trend or exponential trend, respectively. In Section 290 5, we found the admissible value for smoothing parameter a. Thus, we consider the ﬁxed value of 291 a=0.8 and then determine the optimal trend type with ﬁxed values of aand b. 292 5.1.3. Triple exponential Smoothing (Holt-Winters method) 293 This is the most advanced exponential smoothing method, as it is ideal for data with clear trends 294 and seasonality. It has the power to add support for seasonality in a model. There are four important 295 aspects of time series namely level, trend, seasonality, and noise. The level will always be up and 296 down whereas the trend changes in level in some sort of pattern. The commonly observed trends are 297 linear, square, exponential, logarithmic, square root, inverse, and 3rd-degree or higher polynomials. 298 Like the trend in double exponential smoothing, we have two variations for seasonality: 299 Additive method: the seasonal variations are constant, 300 Multiplicative method: the seasonal variations changes with time. 301 5.2. Auto Regressive Integrated Moving Average (ARIMA) 302 Auto-Regressive Integrated Moving Average (ARIMA) model [ 39] is one of the most widely 303 used families of models for time series. These models are a generalization of two processes: An 304 auto-Regressive (AR) process and a Moving Average (MA) process. Some people consider this as a 305 combination of three models by counting differencing as a model. In ARIMA, we initially assume 306 that the time series is stationary; if it is not, then we take the differences between two consecutive 307 observations until the time series becomes stationary. An ARIMA model is classiﬁed by three following 308 parameters: 309 p: number of autoregressive terms, 310 d: number of nonseasonal differences needed to make time series stationary, 311 q: number of lagged forecast errors in the prediction equation. 312 This model considers the independent variable that can inﬂuence our time-series data. In the following 313 subsections, we consider two versions of ARIMA, based on the inclusion of exogenous variables. Both 314 versions are implemented using the pmdarima package in Python. 315 5.2.1. ARIMA without exogenous variables 316 Here, we build an ARIMA model with the count of daily positive cases as the only training data. 317 To optimize the parameters p,d, and q, we use a built-in function known as autoarima rather than 318 deﬁning the explicit values for p,d, and q. The autoarima is mainly used for identifying the most 319 optimal parameters for the ARIMA model. It settles on a single-ﬁtted ARIMA model. This method is 320 completely based on the commonly used R function. 321 5.2.2. ARIMA with exogenous variables 322 As exogenous variables, we use all the independent variables used in Section 4 except for delay 323 variables. The reason to skip this variable is that we created this variable to impose a time factor, which 324 is not required for ARIMA. Autoarima is used here as well. 325 5.2.3. Seasonal ARIMA 326 Seasonal ARIMA (SARIMA) is an ARIMA model in which 327 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 11 of 20 01 Sep 202006 02 03 04 05 date7000075000800008500090000Number of daily cases Actual Predicted without exogenous variables With exogenous variables Figure 4. Comparison of SARIMA models 01 Sep 202006 02 03 04 05 date7250075000775008000082500850008750090000Number of daily confirmed casesSimple Exponential Smoothing Actual Predicted with = 0.2 Predicted with = 0.5 Predicted with = 0.8 Figure 5. Comparison of predicted values with different smoothing parameters a 6. Results and analysis 328 In this section, we review the models with the following metrics for evaluating predictions and 329 also the analysis for each method 330 Mean absolute error (MAE): average of the absolute differences between predicted and actual 331 values. It is used when we care only about the magnitude of the error and not the direction. 332 Mean squared error (MSE): also gives the idea of the magnitude of error, like MAE. It is the 333 average of squared differences between predictions and actual values. 334 Median squared error (MEME): median of squared differences between predicted and actual 335 values. Since the mean is not robust. The mean is much more sensitive to extreme values than 336 the median. Therefore we consider MEME as an alternative evaluation metric. 337 Mean squared log error (MSLE): squared differences between the log-transformed actual and 338 predicted values. It provides the idea of the relative difference between the true and predicted 339 values. 340 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 12 of 20 We compare the different simple exponential smoothing models and we choose a variety of values 341 fora. The resultant predicted values are given in Figure 5. 342 For most of the dates, predicted values from the model with a=0.8are the closest to actual 343 values. Therefore from this family, we choose the simple exponential smoothing model with a=0.8to 344 compare it with other classes of models. 345 The double exponential smoothing method is implemented as shown in Figure 6. As we can see, 346 there is no substantial difference when changing the trend type. So, we select additive trend type and 347 plot for different values for bin Figure 7. 348 01 Sep 202006 02 03 04 05 date65000700007500080000850009000095000Number of daily confirmed casesHolt Method Actual Predicted with = 0.2 and multiplicative trend Predicted with = 0.2 and additive trend Predicted with = 0.8 and multiplicative trend Predicted with = 0.8 and additive trend Figure 6. Comparison of predicted values with different trend smoothing parameters band trend type 01 Sep 202006 02 03 04 05 date65000700007500080000850009000095000Number of daily confirmed casesHolt Method with = 0.8 and additive trend type Figure 7. Comparison of predicted values with different trend smoothing parameters b As indicated in the ﬁgure, there is no admissible choice for b. Therefore, we will consider all three 349 methods with b=0.2, 0.5, and 0.8 in Section 6. 350 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 13 of 20 The predicted values of the triple exponential smoothing method is plotted in Figure 8 for a 351 different type of trend and seasonality. 352 01 Sep 202006 02 03 04 05 date7800080000820008400086000880009000092000Number of daily confirmed casesHolt-Winter Method Actual Predicted with additive seasonal and additive trend Predicted with multiplicative seasonal and additive trend Predicted with additive seasonal and multiplicative trend Predicted with multiplicative seasonal and multiplicative trend Figure 8. Comparison of predicted values with a different type of trend and seasonality As the ﬁgure indicates, the Holt-Winters method with the additive trend and additive seasonality 353 is found to be the best. 354 In Figure 9, we compare both ARIMA models, one without exogenous and one with, against 355 ground truth values. 356 01 Sep 202006 02 03 04 05 date7000075000800008500090000Number of daily cases Figure 9. Comparison of ARIMA models Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 14 of 20 Model Variables used for trainingConﬁguration Evaluation Linear regressionModels with predictors and without predictorsa0.05 Actual values with the predicted values Elastic net regulationRidge Lasso a= 0.2 g=0.1 Mean absolute error (MAE), Mean squared error (MSE), Median squared error (MEME), Mean squared log error (MSLE) Random forest regressornestimators , max features , max depth, min sample split, min samples lea fnestimators =200, max features =auto, max depth = 50, min samples split = 2, min samples lea f= 5Mean absolute error (MAE ), Mean squared error (MSE ), Median squared error (MEME ), Mean squared log error (MSLE ) Recurrent neural network (RNN)uses layers of neurons to map inputs to outputsKeras API Mean absolute error (MAE), Mean squared error (MSE), Median squared error (MEME), Mean squared log error (MSLE) Long short-term memory network (LSTM)RNN k = 1 and Normalize data = 0 and 1Mean absolute error (MAE), Mean squared error (MSE), Median squared error (MEME), Mean squared log error (MSLE) Gated recurrent unit (GRU) Table 2. Analysis of non-time-series predictive models to predict the number of daily positive cases As indicated in the ﬁgure, there is no admissible choice between these two ARIMAs. For some 357 dates, ARIMA without exogenous variables outperforms the one with exogenous variables. Therefore 358 we will consider both models for comparison in Section 6. 359 6.1. A comparative study of models to predict the number of daily positive cases 360 In Sections 4 and 5, we have explored many methods to predict the number of daily positive cases. 361 For many classes of models, we have succeeded in obtaining an optimal model. In this section, we 362 compare all models together with multiple evaluation methods. 363 First, we compared two linear regression models and opted for the model with all predictors. 364 In addition, we calculated the best hyper-parameters within the deﬁned search spaces for elastic net 365 regularization, random forest regressor, and XGBoost regressor families. For each family, we have 366 an optimal model corresponding to the best hyper-parameters. We have also built an LSTM model, 367 forming a total of ﬁve models from Section 4. However, the main disadvantage of the linear regression 368 model is over-ﬁtting. The elastic net regularization can cause a small bias in the model where the 369 prediction is too dependent upon a particular variable. In fact, the random forest algorithm may 370 change considerably by a small change in the data. 371 In Section 5, we explored some time-series forecasting methods. For the simple exponential 372 smoothing method, we have chosen the model with smoothing parameter a=0.8. For the Holt 373 method, we did not obtain anyone’s admissible method. Thus, we decided to have three models 374 with smoothing parameter a=0.8, additive trend type, and corresponding to the trend’s smoothing 375 parameters b=0.2, 0.5 , and 0.8. For Holt-winter’s method, we have selected the one with the additive 376 trend and additive seasonality. For ARIMA family, we have two models with and without exogenous 377 variables. Thus, we have seven models from Section 5. 378 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 15 of 20 Model Variables used for training Conﬁguration Evaluation Exponential smoothingSimple exponential smoothing, Double exponential smoothing (Holt method), Triple exponential Smoothing (Holt-Winters method)TSA (Time Series Analysis)MAE, MSE, MESE, MSLE Auto-Regressive Integrated Moving Averagep : number of autoregressive terms, d : number of nonseasonal differences needed to make time series stationary, q : number of lagged forecast errors in the prediction equationARIMA with and without exogenous variables. Seasonal ARIMAMAE, MSE, MESE, MSLE Table 3. Analysis of time-series predictive models to predict the number of daily positive cases Model MAE MSE MESE MSLE Linear regression 4804.8860 6172.9314 2723.2462 0.0054 Elastic net regularization 7265.5959 8245.1422 5342.1506 0.0100 Random forest regressor 11351.8833 11827.2160 9998.6333 0.0220 XGBoost regressor 10130.6125 10566.9168 9346.6719 0.0173 Simple exponential smoothing 4507.6726 5045.6480 4851.4896 0.0040 Holt with b=0.2 3552.8030 4266.5536 2670.3701 0.0030 Holt with b=0.5 4168.4262 5401.2516 2615.0862 0.0050 Holt with b=0.8 5120.1373 6533.2962 5305.5930 0.0076 Holt-winters 1629.8258 2253.0399 506.1216 0.0007 Arima 4918.0511 5459.2333 4427.1078 0.0048 Arima with exogenous variables 4061.0362 4766.3267 3037.7827 0.0033 Sarima 4918.0511 5459.2333 4427.1078 0.0048 RNN 7604.9391 7895.1482 8395.9531 0.0098 GRU 4490.1203 5020.4703 5372.7188 0.0039 LSTM 6238.7969 6588.8430 7022.9141 0.0067 Table 4. Comparison of models from different classes with different evaluation metrics 379 01 Sep 202006 02 03 04 05 date78000800008200084000860008800090000Number of daily deathsPrediction from final model Actual Predicted Figure 10. Prediction of daily conﬁrmed cases for ﬁrst six days of September 2020 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 16 of 20 500010000 50001000015000 010000 Linear regression Elastic net regularization Random forest regressor XGBoost regressor Simple exponential smoothing Holt with =0.2 Holt with =0.5 Holt with =0.8 Holt-winters Arima Arima with exogenous variables Sarima Sarima with exogenous variables RNN GRU LSTM0.0000.025Sharing both axes MAE MSE MEME MSLE Figure 11. Comparison of models from different classes with different evaluation metrics 6.2. Predicted number of daily deaths 380 In this section, we predict daily deaths on the same line using the methods from previous sections. 381 We provide the ﬁnal results in the following table and graphs. There are different methods to handle 382 the computational cost and missing data. In these models such as XGBoost, and Random-forest, the 383 missing values are interpreted as data that contain information (ie, data that are missing for a reason) 384 instead of data that are missing at random. 385 Model MAE MSE MESE MSLE predicted_lm1 75.9125 84.5736 87.1517 0.0067 predicted_el1 47.9260 55.1785 51.6010 0.0028 predicted_rf1 136.8500 152.6284 167.7500 0.0244 predicted_xgb1 196.8140 230.1556 222.5404 0.0661 predicted_ses_0.8 64.3146 93.7673 37.1779 0.0096 predicted_holt_0.2 73.1619 104.6138 44.9703 0.0123 predicted_holt_0.5 81.5152 125.5030 53.2248 0.0187 predicted_holt_0.8 94.6749 138.4354 49.2252 0.0234 predicted_hw1 29.8014 39.5107 18.9877 0.0014 predicted_autoarima 51.8795 55.8371 57.1838 0.0029 predicted_autoarima_ex 69.7804 75.0977 70.4612 0.0053 predicted_lstm 189.0365 191.7943 184.3370 0.0397 Table 5. Comparison of models from different classes with different evaluation metrics 386 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 17 of 20 Figure 12. Comparison of models from different classes with different evaluation metrics Figure 13. Prediction of daily deaths for ﬁrst six days of September 2020 7. Conclusions 387 In this paper, we aimed to forecast COVID-19 spread by predicting the number of daily positive 388 cases and daily deaths. For the case study, we considered the data of India at the national level. Any 389 country with any geographic level can be analyzed in the same manner according to the availability of 390 data. We used many different models to predict the number of daily positive cases and the number 391 Version December 16, 2022 submitted to Mach. Learn. Knowl. Extr. 18 of 20 of daily deaths. From there, we compared with respect to different evaluation metrics. Our study 392 indicates that more complex models do not outperform simpler ones for COVID-19 spread. Regardless, 393 we have observed that the Holt-winters model is an optimal model for both the number of daily 394 positive cases and the number of daily deaths. 395