The novel coronavirus (nCOV-19) and the associated disease, coronavirus disease 2019 (COVID-19), presumably originated from Wuhan, China in December 2019 1 . The rapid spread of COVID-19 worldwide caused the World Health Organization (WHO) to declare it as a pandemic in March 2020 1 . As of May 25, 2020, more than 1.5 million people in the United States and nearly 6 million people worldwide have been confirmed with COVID-19 2 . Compared to other highly-contagious previously-identified coronavirusrelated diseases, such as Severe Acute Respiratory Syndrome (SARS) and Middle East Respiratory Syndrome (MERS), COVID-19 appears to be more contagious according to preliminary studies 3 .
In light of the present circumstances, it is of paramount importance for policymakers and healthcare professionals to plan for the future and be adequately equipped for circumstances that may arise due to the rapid spread of COVID-19. A crucial part of planning ahead in this scenario is to forecast the number of cases in the future, both in the short term and in the long term. While long term predictions and forecasting can be made with epidemiological models, short term predictions can be performed by time series analysis. In this study, we evaluate the performance of several time series analysis and forecasting models for predicting country-wise COVID-19 cases in the short term.
For forecasting using time series analysis, several models, such as Long Short-Term Memory (LSTM) networks 4 and Auto-regressive Integrated Moving Average (ARIMA) model 5 , have been widely used. Some preliminary studies for COVID-19 time series forecasting using ARIMA have also been done 6, 7, 8 . However, two avenues of research -forecasting for multiple days in the future and the comparison of different LSTM models with ARIMA model -are yet to be widely explored.
This study explores the performance of several LSTM models -vanilla LSTM 4 , stacked LSTM 9 , bidirectional LSTM 10 , convolutional neural network LSTM (CNN LSTM) 11 , and ConvLSTM 12 -along with ARIMA model. Moreover, since it is worthwhile to evaluate the performance of these time series forecasting methods for predicting multiple values in the future rather than simply single forecasts, we introduce the concept of k-period performance metrics. These k-period performance metrics extend existing performance metrics used for evaluating single forecasts to the general case of evaluating multiple forecasts in the future. While the proposed concept of k-period performance metrics can be used to extend any performance metric, we specifically choose Mean Absolute Percentage Error (MAPE) 13 and Median Symmetric Accuracy (MdSA) 14 for deriving the k-period performance metrics used in this study.
The main contributions of this study are:
 We evaluate and compare the performance of several LSTM architectures vis-à-vis the ARIMA model in time series analysis and forecasting of country-wise COVID-19 cases for 4 countries.
 We propose and introduce k-period performance metrics for evaluating and comparing the performance of time series forecasting algorithms where forecasts are made for multiple time periods in the future. While this definition can be used for any performance metric in time series forecasting, we extend the definition of MAPE and MdSA metrics to define k-period kMAPE (kMAPE) and k-period MdSA (kMdSA), respectively, in particular.
Time series data including the number of confirmed cases, recovered, and deaths due to COVID-19 for several countries have been made available publicly by researchers at Johns Hopkins University 2 (https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data). The dataset consists of time series information about the total daily cumulative number of confirmed cases, recovered, and deaths for almost all countries in the world as well as for individual US (United States of America) states. The data is collected from different sources, such as WHO, the Centers for Disease Control and Prevention (CDC) in the US, and the Chinese Center for Disease Control and Prevention (CCDC) to name a few. In our study, we have used the time series for daily cumulative confirmed COVID-19 cases from 4 countries -US, Italy, Spain, and Germany -from this dataset.
The dataset is updated daily with new information. For our experiments, we have used all data until May 25, 2020. Further, we have truncated the data in the beginning so that all time series points with number of cases less than 100 have been removed. In other words, the first time point in our time series corresponds to the first date when the number of confirmed cases in that particular country was 100 or above. This step has been done to ensure that the noise incorporated in the time series data due to lack of testing in several countries in the initial stages of the disease does not affect the time series data in general.
In this study, we have used several LSTM architectures and ARIMA model for the forecasting of confirmed COVID-19 cases in different countries. These models are briefly described here.
Long short-term memory (LSTM) networks 4 are recurrent neural networks (RNN) widely used in deep learning. LSTMs were designed to process sequences of data and improved upon traditional RNN by using memory cells that can store information in memory for long sequences, and a set of gates to control the flow of this memory information. These innovations allow LSTM to learn longer term dependencies in sequential data.
The following LSTM architectures were used in our experiments:
 Vanilla LSTM: A vanilla LSTM is the simplest LSTM model with a single LSTM unit having multiple neurons. In our experiments, the vanilla LSTM had = 100 neurons.  Stacked LSTM: A stacked LSTM 9 utilizes a stacked architecture of two or more vanilla LSTMs. In our experiments, two vanilla LSTMs with = 50 were stacked together.  Bidirectional LSTM: A bidirectional LSTM 10 architecture splits the neurons of a regular LSTM into two directions -positive time direction for forward states and negative time direction for backward states. As a result, input information from both the past and the future of the current time point is used. In our experiments, bidirectional LSTMs with = 50 neurons were used.
 CNN LSTM: Also known as Long-term Recurrent Convolutional Network (LRCN), CNN LSTM 11 architecture incorporates 1-D spatial structure along with temporal structure to learn patterns.
Here, 1-D spatial structures are created by grouping together a few time points and considering them as 1-D spatial structures. This unique architecture allows the network to identify patterns in small groups of time points along with the actual temporal patterns in the time series. In our experiments, the CNN LSTM consisted of = 50 neurons in the LSTM and = 64 convolutional filters along with a max pooling layer in the CNN.  ConvLSTM: ConvLSTM 12 (Convolutional LSTM) architecture allows convolutions at the gates of the LSTM to capture spatio-temporal patterns. Although it learns spatio-temporal patterns similar to CNN LSTM, the manner in which it does so is different (by using convolutions at the gates). In our experiments, the ConvLSTM used = 50 neurons and = 64 convolutional filters.
LSTMs are notoriously sensitive to initialization. To tackle the initialization problem, the models were trained on 100 different random initializations and the one with the best (least) validation kMAPE (k-day Mean Absolute Percentage Error) was selected for obtaining the performance on the test set. The data were divided into training, validation, and test sets of size − 20, 10, and 10 respectively (where is the number of time series sequences obtained from the data).
Further, we trained the LSTMs on varying sizes of input sequences. We used sequence sizes of = 3, 6, 10, 15 in our experiments. Adam optimizer with a learning rate of 0.1 was used for training the LSTM models and mean squared error was used as the loss function.
The Auto-regressive Integrated Moving Average (ARIMA) model is a time series analysis and forecasting method. The ARIMA model has three components: An ARIMA model has three parameters:
1. , the number of lag observations in the model (or the lag order) 2. , the number of times the raw observations are differenced (or the degree of differencing) 3. , the size of the window for moving average (or the order of moving average) In our experiments, these parameters were empirically chosen to be the following values: = 1, = 2, = 2. Values of the parameters were evaluated using a grid search on each parameter between values 1 and 3, and subsequently chosen empirically to reflect the best performance of the ARIMA model. For all experiments with the ARIMA model, the sequence length was kept as = 15.
Besides the commonly-used 1-day forecast into the future, we also performed experiments for predicting the time series values for periods (i.e., days in this case) in the future. We used = 1, 3, 5 in our
experiments. For both the LSTM models and the ARIMA model, the prediction of the 1st day in the future was added to the input sequence for predicting the 2nd day and so on. In general, if the input sequence is [ 1 , 2 , … , −1 , ] , where the size of the sequence is , the forecast for the next day is ̂+ 1 . For predicting the subsequent day, the input sequence is modified to [ 2 , 3 , … , −1 , ,̂+ 1 ] and the forecast for the 2nd day in the future is ̂+ 2 . In this way, forecasting for = days in the future gives us the k-day predictions as [̂+ 1 ,̂+ 2 , … ,̂+ ].
For evaluating the performance for time series forecasting models, several performance metrics are commonly used 15 . In this study, we have chosen the Mean Absolute Percentage Error (MAPE) and Median Symmetric Accuracy (MdSA) 14 , both of which are measures of error in prediction and expressed as a percentage. MAPE is given by,
where denotes the true value at time point , ̂ denotes the predicted value at time point , and denotes the number of 1-day predictions made.
MdSA is based on the log accuracy ratio and is given by,
where =1 ( ) denotes the median of the numbers , = 1, 2, … , . MdSA possesses several symmetry and robustness properties 14 , thus making it a good choice as a metric in our study. Further, the use of the median operation instead of the most commonly-used mean operation renders this metric immune to the effect of outliers, making it more robust compared to MAPE and other similar metrics.
Both the described metrics are designed for evaluating multiple 1-period forecasts. In our study, we also seek to predict and evaluate multiple k-period forecasts. Towards this goal, we extend the definition of any performance evaluation metric to incorporate the evaluation of k-period forecasts. Consider sequences [ 1 , 2 , … , ], = 1, 2, … , , each of length . For each of these sequences, we can obtain k-period predictions as described before. These predictions are given by [̂+ 1 ,̂+ 2 , … ,̂+ ], = 1, 2, … , , each of length .
We define the general k-period metric as, Specifically, we define the k-period MAPE (kMAPE) as, Similarly, k-period MdSA (kMdSA) is given by,
where the notations are the same as in Eqn. 2 and Eqn. 4.
It is to be noted that in the case of 1-day forecast in the future, the k-day performance metrics are identical to the original definition of the performance metrics. Figure 1 shows how performance metrics are calculated and how they can be extended to define the k-period performance metrics. 
The performance of the different LSTM models and ARIMA models were compared using kMAPE and kMdSA as described. For LSTM models, input sequence length = 3, 6, 10, 15 were used in our experiments while for the ARIMA model, only = 15 was used. The performance of the LSTM models with = 6 and the ARIMA model with = 15 are compared using the kMAPE metric in Table 3 . The values for ARIMA here are identical to the ones shown in Table 1 as they have the same = 15 but are repeated here for ease of comparison. Table 4 shows the comparison of the performance of the LSTM models with = 6 and the ARIMA model with = 15 using the kMdSA metric. The values for ARIMA are repeated in Table 4 from Table 2 
Short-and long-term forecasts and predictions of the number of cases in any epidemic or pandemic can be a useful tool for several government agencies and policymakers in different countries to prepare for the progression and spread of a disease. Given the rapid growth of COVID-19 cases worldwide, it is imperative that these tools be widely available to relevant agencies and personnel to aid in preparation (e.g. number of hospital beds and ventilators, requirement and availability of healthcare workers, and the demand for personal protective equipment (PPE)) which would be directly reflected in better outcomes for patients as well as the number of lives saved due to a deadly or fatal disease.
In this work, we have compared the ARIMA model with variants of a supervised model (LSTM) for forecasting country-wise COVID-19 cases for 4 different countries -United States of America, Italy, Spain, and Germany. Similar analyses can be performed with the data from any country and even at more granular levels, such as states, counties, cities, and districts to obtain an insight into the preparation required to combat the spread of the disease.
Our results indicate ARIMA performs reasonably well and is often superior to the variants of the supervised LSTM models used in this study. This makes the ARIMA model suitable for baseline performance comparisons and can be used for developing better models for time series analysis and forecasting.
As expected, the performance of all the models deteriorates slightly as the number of forecast periods (i.e., days in our case) in the future is increased. This drop in performance is attributed to the fact that the prediction for the 2 nd day uses the prediction for the 1 st day in its input sequence, and so on. As a result, the errors in prediction are propagated as the forecasts are made for more and more days into the future. This is the reason why time series forecasting models presented in this study are unsuitable for long-term forecasting.
It is also interesting to note that the performance of the vanilla LSTM is comparable and often better than the performance of the more complex LSTM models discussed here. This pattern indicates that more complex models are likely not necessary for this forecasting problem.
Since we have used only data starting with the day when 100 cases were reported for a country, the number of training sequences for each country has been different. Given the nature of deep learning models such as LSTM, improvements in performance often correlate with greater number of training data. As a result, it can be seen how the forecast performance of the LSTM models for US, Italy, and Germany are better than those for Spain since the number of training samples for Spain were fewer compared to the other countries. It is expected that the performance of all the LSTM models would improve over time as more and more data becomes available for all countries worldwide, thus making more training data available. The ARIMA model is, however, immune to these problems.
In light of the training set size, the choice of models for forecasting becomes obvious. If there are only a small number of training samples available, the ARIMA model can be chosen for forecasting. However, once more and more data is available for training the LSTM models, the best among these models can be chosen if its performance supersedes that of the ARIMA model.
On the other hand, whereas LSTM models in our experiment were fairly stable even for very short input sequences of length, = 3 and = 6, the sequence length for the ARIMA model is kept at = 15, since a short sequence length is deemed to be insufficient for exploring the correlation in the sequence by ARIMA. This makes LSTM models more attractive if we are dealing with sequences of shorter length.
From Figure 2 , it is also apparent that ARIMA models slightly overestimate the numbers while LSTM models are more likely to underestimate in their forecasts when compared to the actual ground truth data.
There are some limitations to this work. Out of several possible LSTM models, only five were evaluated in this study. Moreover, there is a plethora of time series analysis and forecasting algorithms that can be evaluated for this particular problem. The time series analysis and forecasting presented in this study are restricted to short term forecasts only. Forecasting over longer terms would require more parameters to be considered in the modeling that is beyond the scope of this work.
The data that has been made available by the public agencies of different countries might be prone to errors and the limited testing capabilities for COVID-19 in certain countries are likely to affect the results of the algorithms, thus yielding lower estimates. Moreover, the availability of testing and rate of infection within a country might vary between different demographics due to economic conditions, the cost of testing, and the availability of insurance among the different demographics 16 . The present study also does not model the changes in the public health policy of different countries which would possibly affect the future estimates for the numbers of patients suffering from COVID-19. However, it might be noted that these changes in public health policy would most likely affect long term forecasts rather than short term trends.
COVID-19 is a global public health crisis. Forecasts and predictions of cases for the future are essential for government agencies and policymakers to combat the situation by planning for the needs of healthcare personnel and the availability of healthcare equipment and PPEs. In this study, we have presented the time series analysis and forecasting of the number of COVID-19 cases in 4 countries -United States, Italy, Spain, and Germany -using the ARIMA model and several LSTM architectures. We have also proposed the concept of k-period performance metrics designed to evaluate the performance of a time series analysis model that forecasts for multiple periods in the future. In particular, we have extended the definition of MAPE and MdSA to define k-period MAPE (kMAPE) and k-period MdSA (kMdSA). Our results indicate that LSTM models perform comparably with the ARIMA model. However, both models have their advantages and disadvantages, and the choice of models is dictated by the availability of training samples as well as the sequence length. Further research would include the development and evaluation of other time series forecasting models as well as the evaluation of these models as more data becomes available.
Tables S1 and S2 show the performance of the different LSTM models with = 10 and the ARIMA model with = 15 using the kMAPE and the kMdSA metrics respectively. 