In December 2019, clusters of pneumonia cases caused by the novel Coronavirus were identified at the Wuhan, Hubei province in China [58, 48] after almost hundred years of the 1918 Spanish flu [118] . Soon after the emergence of the novel beta coronavirus, World Health Organization (WHO) characterized this contagious disease as a "global pandemic" due to its rapid spread worldwide [99] . Many scientists have attempted to make forecasts about its impact. However, despite involving many excellent modelers, best intentions, and highly sophisticated tools, forecasting COVID-19 pandemics is harder [64] , and this is primarily due to the following major factors:
-Very less amount of data is available; -Less understanding of the factors that contribute to it; -Model accuracy is constrained by our knowledge of the virus, however. With an emerging disease such as COVID-19, many transmission-related biologic features are hard to measure and remain unknown; -The most obvious source of uncertainty affecting all models is that we don't know how many people are or have been infected;
-Ongoing issues with virologic testing mean that we are certainly missing a substantial number of cases, so models fitted to confirmed cases are likely to be highly uncertain [55] ; -The problem of using confirmed cases to fit models is further complicated because the fraction of confirmed cases is spatially heterogeneous and time-varying [123] ; -Finally, many parameters associated with COVID-19 transmission are poorly understood.
Amid enormous uncertainty about the future of the COVID-19 pandemic, statistical, machine learning, and epidemiological models are critical forecasting tools for policymakers, clinicians, and public health practitioners [26, 76, 126, 36, 73, 130] . COVID-19 modeling studies generally follow one of two general approaches that we will refer to as forecasting models and mechanistic models. Although there are hybrid approaches, these two model types tend to address different questions on different time scales, and they deal differently with uncertainty [25] . Compartmental epidemiological models have been developed over nearly a century and are well tested on data from past epidemics. These models are based on modeling the actual infection process and are useful for predicting long-term trajectories of the epidemic curves [25] . Short-term Forecasting models are often statistical, fitting a line or curve to data and extrapolating from there -like seeing a pattern in a sequence of numbers and guessing the next number, without incorporating the process that produces the pattern [23, 24, 26] . Well constructed statistical frameworks can be used for short-term forecasts, using machine learning or regression. In statistical models, the uncertainty of the prediction is generally presented as statistically computed prediction intervals around an estimate [51, 65] . Given that what happens a month from now will depend on what happens in the interim, the estimated uncertainty should increase as you look further into the future. These models yield quantitative projections that policymakers may need to allocate resources or plan interventions in the short-term.
Forecasting time series datasets have been a traditional research topic for decades, and various models have been developed to improve forecasting accuracy [27, 10, 49] . There are numerous methods available to forecast time series, including traditional statistical models and machine learning algorithms, providing many options for modelers working on epidemiological forecasting [24, 25, 20, 26, 80, 22, 97] . Many research efforts have focused on developing a universal forecasting model but failed, which is also evident from the "No Free Lunch Theorem" [125] . This chapter focuses on assessing popularly used short-term forecasting (nowcasting) models for COVID-19 from an empirical perspective. The findings of this chapter will fill the gap in the literature of nowcasting of COVID-19 by comparing various forecasting methods, understanding global characteristics of pandemic data, and discovering real challenges for pandemic forecasters.
The upcoming sections present a collection of recent findings on COVID-19 forecasting. Additionally, twenty nowcasting (statistical, machine learning, and hybrid) models are assessed for five countries of the United States of America (USA), India, Brazil, Russia, and Peru. Finally, some recommendations for policy-making decisions and limitations of these forecasting tools have been discussed.
Researchers face unprecedented challenges during this global pandemic to forecast future real-time cases with traditional mathematical, statistical, forecasting, and machine learning tools [76, 126, 36, 73, 130] . Studies in March with simple yet powerful forecasting methods like the exponential smoothing model predicted cases ten days ahead that, despite the positive bias, had reasonable forecast error [91] . Early linear and exponential model forecasts for better preparation regarding hospital beds, ICU admission estimation, resource allocation, emergency funding, and proposing strong containment measures were conducted [46] that projected about 869 ICU and 14542 ICU admissions in Italy for March 20, 2020. Health-care workers had to go through the immense mental stress left with a formidable choice of prioritizing young and healthy adults over the elderly for allocation of life support, mostly unwanted ignoring of those who are extremely unlikely to survive [35, 100] . Real estimates of mortality with 14-days delay demonstrated underestimating of the COVID-19 outbreak and indicated a grave future with a global case fatality rate (CFR) of 5.7% in March [13] . The contact tracing, quarantine, and isolation efforts have a differential effect on the mortality due to COVID-19 among countries. Even though it seems that the CFR of COVID-19 is less compared to other deadly epidemics, there are concerns about it being eventually returning as the seasonal flu, causing a second wave or future pandemic [89, 95] .
Mechanistic models, like the Susceptible-Exposed-Infectious-Recovered (SEIR) frameworks, try to mimic the way COVID-19 spreads and are used to forecast or simulate future transmission scenarios under various assumptions about parameters governing the transmission, disease, and immunity [56, 52, 8, 29, 77] . Mechanistic modeling is one of the only ways to explore possible long-term epidemiologic outcomes [7] . For example, the model from Ferguson et al. [40] that has been used to guide policy responses in the United States and Britain examines how many COVID-19 deaths may occur over the next two years under various social distancing measures. Kissler et al. [70] ask whether we can expect seasonal, recurrent epidemics if immunity against novel coronavirus functions similarly to immunity against the milder coronaviruses that we transmit seasonally. In a detailed mechanistic model of Boston-area transmission, Aleta et al. [4] simulate various lockdown "exit strategies". These models are a way to formalize what we know about the viral transmission and explore possible futures of a system that involves nonlinear interactions, which is almost impossible to do using intuition alone [54, 85] . Although these epidemiological models are useful for estimating the dynamics of transmission, targeting resources, and evaluating the impact of intervention strategies, the models require parameters and depend on many assumptions.
Several statistical and machine learning methods for real-time forecasting of the new and cumulative confirmed cases of COVID-19 are developed to overcome limitations of the epidemiological model approaches and assist public health planning and policy-making [25, 91, 6, 26, 23] . Real-time forecasting with foretelling predictions is required to reach a statistically validated conjecture in this current health crisis. Some of the leading-edge research concerning real-time projections of COVID-19 confirmed cases, recovered cases, and mortality using statistical, machine learning, and mathematical time series modeling are given in Table 1 . 
A univariate time series is the simplest form of temporal data and is a sequence of real numbers collected regularly over time, where each number represents a value [28, 18] . There are broadly two major steps involved in univariate time series forecasting [60] :
-Studying the global characteristics of the time series data; -Analysis of data with the 'best-fitted' forecasting model.
Understanding the global characteristics of pandemic confirmed cases data can help forecasters determine what kind of forecasting method will be appropriate for the given situation [120] . As such, we aim to perform a meaningful data analysis, including the study of time series characteristics, to provide a suitable and comprehensive knowledge foundation for the future step of selecting an apt forecasting method. Thus, we take the path of using statistical measures to understand pandemic time series characteristics to assist method selection and data analysis. These characteristics will carry summarized information of the time series, capturing the 'global picture' of the datasets. Based on the recommendation of [32, 122, 75, 74] , we study several classical and advanced time series characteristics of COVID-19 data. This study considers eight global characteristics of the time series: periodicity, stationarity, serial correlation, skewness, kurtosis, nonlinearity, long-term dependence, and chaos. This collection of measures provides quantified descriptions and gives a rich portrait of the pandemic time-series' nature. A brief description of these statistical and advanced time-series measures are given below.
A seasonal pattern exists when a time series is influenced by seasonal factors, such as the month of the year or day of the week. The seasonality of a time series is defined as a pattern that repeats itself over fixed intervals of time [18] . In general, the seasonality can be found by identifying a large autocorrelation coefficient or a large partial autocorrelation coefficient at the seasonal lag. Since the periodicity is very important for determining the seasonality and examining the cyclic pattern of the time series, the periodicity feature extraction becomes a necessity. Unfortunately, many time series available from the dataset in different domains do not always have known frequency or regular periodicity. Seasonal time series are sometimes also called cyclic series, although there is a significant distinction between them. Cyclic data have varying frequency lengths, but seasonality is of a fixed length over each period. For time series with no seasonal pattern, the frequency is set to 1. The seasonality is tested using the 'stl' function within the "stats" package in R statistical software [60] .
Stationarity is the foremost fundamental statistical property tested for in time series analysis because most statistical models require that the underlying generating processes be stationary [27] . Stationarity means that a time series (or rather the process rendering it) do not change over time. In statistics, a unit root test tests whether a time series variable is non-stationary and possesses a unit root [93] . The null hypothesis is generally defined as the presence of a unit root, and the alternative hypothesis is either stationarity, trend stationarity, or explosive root depending on the test used. In econometrics, Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests are used for testing a null hypothesis that an observable time series is stationary around a deterministic trend (that is, trend-stationary) against the alternative of a unit root [108] . The KPSS test is done using the 'kpss.test' function within the "tseries" package in R statistical software [117] .
Serial correlation is the relationship between a variable and a lagged version of itself over various time intervals. Serial correlation occurs in time-series studies when the errors associated with a given time period carry over into future time periods [18] . We have used Box-Pierce statistics [19] in our approach to estimate the serial correlation measure and extract the measures from COVID-19 data. The Box-Pierce statistic was designed by Box and Pierce in 1970 for testing residuals from a forecast model [122] . It is a common portmanteau test for computing the measure. The mathematical formula of the Box-Pierce statistic is as follows:
where n is the length of the time series, h is the maximum lag being considered (usually h is chosen as 20) , and r k is the autocorrelation function. The portmanteau test is done using the 'Box.test' function within the "stats" package in R statistical software [63] .
Nonlinear time series models have been used extensively to model complex dynamics not adequately represented by linear models [67] . Nonlinearity is one important time series characteristic to determine the selection of an appropriate forecasting method. [115] There are many approaches to test the nonlinearity in time series models, including a nonparametric kernel test and a Neural Network test [119] . In the comparative studies between these two approaches, the Neural Network test has been reported with better reliability [122] . In this research, we used Teräsvirta's neural network test [112] for measuring time series data nonlinearity. It has been widely accepted and reported that it can correctly model the nonlinear structure of the data [113] . It is a test for neglected nonlinearity, likely to have power against a range of alternatives based on the NN model (augmented single-hidden-layer feedforward neural network model). This takes large values when the series is nonlinear and values near zero when the series is linear. The test is done using the 'nonlinearityTest' function within the "nonlinearTseries" package in R statistical software [42] .
Skewness is a measure of symmetry, or more precisely, the lack of symmetry. A distribution, or dataset, is symmetric if it looks the same to the left and the right of the center point [122] . A skewness measure is used to characterize the degree of asymmetry of values around the mean value [83] . For univariate data Y t , the skewness coefficient is
whereȲ is the mean, σ is the standard deviation, and n is the number of data points. The skewness for a normal distribution is zero, and any symmetric data should have the skewness near zero. Negative values for the skewness indicate data that are skewed left, and positive values for the skewness indicate data that are skewed right. In other words, left skewness means that the left tail is heavier than the right tail. Similarly, right skewness means the right tail is heavier than the left tail [69] . Skewness is calculated using the 'skewness' function within the "e1071" package in R statistical software [81] .
Kurtosis is a measure of whether the data are peaked or flat, relative to a normal distribution [83] . A dataset with high kurtosis tends to have a distinct peak near the mean, decline rather rapidly, and have heavy tails. Datasets with low kurtosis tend to have a flat top near the mean rather than a sharp peak. For a univariate time series Y t , the kurtosis coefficient is 1
The kurtosis for a standard normal distribution is 3. Therefore, the excess kurtosis is defined as
So, the standard normal distribution has an excess kurtosis of zero. Positive kurtosis indicates a 'peaked' distribution and negative kurtosis indicates a 'flat' distribution [47] . Kurtosis is calculated using the 'kurtosis' function within the "Performance-Analytics" package in R statistical software [90] .
Processes with long-range dependence have attracted a good deal of attention from a probabilistic perspective in time series analysis [98] . With such increasing importance of the 'self-similarity' or 'long-range dependence' as one of the time series characteristics, we study this feature into the group of pandemic data characteristics. The definition of self-similarity is most related to the self-similarity parameter, also called Hurst exponent (H) [15] . The class of autoregressive fractionally integrated moving average (ARFIMA) processes [44] is a good estimation method for computing H. In an ARIMA(p, d, q), p is the order of AR, d is the degree first differencing involved, and q is the order of MA. If the time series is suspected of exhibiting long-range dependency, parameter d may be replaced by certain non-integer values in the ARFIMA model [21] . We fit an ARFIMA(0, d, 0) to the maximum likelihood, which is approximated by using the fast and accurate method of Haslett and Raftery [50] . We then estimate the Hurst parameter using the relation H = d + 0.5. The self-similarity feature can only be detected in the RAW data of the time series.
The value of H can be obtained using the 'hurstexp' function within the "pracma" package in R statistical software [16] .
Many systems in nature that were previously considered random processes are now categorized as chaotic systems. For several years, Lyapunov Characteristic Exponents are of interest in the study of dynamical systems to characterize quantitatively their stochasticity properties, related essentially to the exponential divergence of nearby orbits [39] . Nonlinear dynamical systems often exhibit chaos, characterized by sensitive dependence on initial values, or more precisely by a positive Lyapunov Exponent (LE) [38] . Recognizing and quantifying chaos in time series are essential steps toward understanding the nature of random behavior and revealing the extent to which short-term forecasts may be improved [53] . LE, as a measure of the divergence of nearby trajectories, has been used to qualifying chaos by giving a quantitative value [14] . The algorithm of computing LE from time-series is applied to continuous dynamical systems in an n-dimensional phase space [101] . LE is calculated using the 'Lyapunov exponent' function within the "tseriesChaos" package in R statistical software [9] .
Time series forecasting models work by taking a series of historical observations and extrapolating future patterns. These are great when the data are accurate; the future is similar to the past. Forecasting tools are designed to predict possible future alternatives and help current planing and decision making [10] . There are essentially three general approaches to forecasting a time series [82] :
1. Generating forecasts from an individual model; 2. Combining forecasts from many models (forecast model averaging); 3. Hybrid experts for time series forecasting.
Single (individual) forecasting models are either traditional statistical methods or modern machine learning tools. We study ten popularly used single forecasting models from classical time series, advanced statistics, and machine learning literature. There has been a vast literature on the forecast combinations motivated by the seminal work of Bates & Granger [12] and followed by a plethora of empirical applications showing that combination forecasts are often superior to their counterparts (see, [17, 114] , for example). Combining forecasts using a weighted average is considered a successful way of hedging against the risk of selecting a misspecified model [31] . A significant challenge is in choosing an appropriate set of weights, and many attempts to do this have been worse than simply using equal weightssomething that has become known as the "forecast combination puzzle" (see, for example, [109] ). To overcome this, hybrid models became popular with the seminal work of Zhang [127] and further extended for epidemic forecasting in [24, 26, 23] . The forecasting methods can be briefly reviewed and organized in the architecture shown in Figure 1 . 
The autoregressive integrated moving average (ARIMA) is one of the well-known linear models in time-series forecasting, developed in the early 1970s [18] . It is widely used to track linear tendencies in stationary time-series data. It is denoted by ARIMA(p, d, q), where the three components have significant meanings. The parameters p and q represent the order of AR and MA models, respectively, and d denotes the level of differencing to convert nonstationary data into stationary time series [78] . ARIMA model can be mathematically expressed as follows:
where y t denotes the actual value of the variable at time t, ǫ t denotes the random error at time t, β i and α j are the coefficients of the model. Some necessary steps to be followed for any given time-series dataset to build an ARIMA model are as follows:
-Identification of the model (achieving stationarity).
-Use autocorrelation function (ACF) and partial ACF plots to select the AR and MA model parameters, respectively, and finally estimate model parameters for the ARIMA model. -The 'best-fitted' forecasting model can be found using the Akaike Information Criteria (AIC) or the Bayesian Information Criteria (BIC). Finally, one checks the model diagnostics to measure its performance.
An implementation in R statistical software is available using the 'auto.arima' function under the "forecast" package, which returns the 'best' ARIMA model according to either AIC or BIC values [61] .
Wavelet analysis is a mathematical tool that can reveal information within the signals in both the time and scale (frequency) domains. This property overcomes the primary drawback of Fourier analysis, and wavelet transforms the original signal data (especially in the time domain) into a different domain for data analysis and processing. Wavelet-based models are most suitable for nonstationary data, unlike standard ARIMA. Most epidemic time-series datasets are nonstationary; therefore, wavelet transforms are used as a forecasting model for these datasets [26] . When conducting wavelet analysis in the context of time series analysis [5] , the selection of the optimal number of decomposition levels is vital to determine the performance of the model in the wavelet domain. The following formula for the number of decomposition levels, W L = int[log(n)], is used to select the number of decomposition levels, where n is the time-series length.
The wavelet-based ARIMA (WARIMA) model transforms the time series data by using a hybrid maximal overlap discrete wavelet transform (MODWT) algorithm with a 'haar' filter [88] . Daubechies wavelets can produce identical events across the observed time series in so many fashions that most other time series prediction models cannot recognize. The necessary steps of a wavelet-based forecasting model, defined by [5] , are as follows. Firstly, the Daubechies wavelet transformation and a decomposition level are applied to the nonstationary time series data. Secondly, the series is reconstructed by removing the high-frequency component, using the wavelet denoising method. Lastly, an appropriate ARIMA model is applied to the reconstructed series to generate out-of-sample forecasts of the given time series data. Wavelets were first considered as a family of functions by Morlet [121] , constructed from the translations and dilation of a single function, which is called "Mother Wavelet". These wavelets are defined as follows:
where the parameter m ( = 0) is denoted as the scaling parameter or scale, and it measures the degree of compression. The parameter n is used to determine the time location of the wavelet, and it is called the translation parameter. If the value |m| < 1, then the wavelet in m is a compressed version (smaller support is the time domain) of the mother wavelet and primarily corresponds to higher frequencies, and when |m| > 1, then φ ( m, n)(t) has larger time width than φ(t) and corresponds to lower frequencies. Hence wavelets have time width adopted to their frequencies, which is the main reason behind the success of the Morlet wavelets in signal processing and time-frequency signal analysis [86] . An implementation of the WARIMA model is available using the 'WaveletFittingarma' function under the "WaveletArima" package in R statistical software [87] .
Fractionally autoregressive integrated moving average or autoregressive fractionally integrated moving average models are the generalized version ARIMA model in time series forecasting, which allow non-integer values of the differencing parameter [44] . It may sometimes happen that our time-series data is not stationary, but when we try differencing with parameter d taking the value to be an integer, it may over difference it. To overcome this problem, it is necessary to difference the time series data using a fractional value. These models are useful in modeling time series, which has deviations from the long-run mean decay more slowly than an exponential decay; these models can deal with time-series data having long memory [94] . ARFIMA models can be mathematically expressed as follows:
where B is is the backshift operator, p, q are ARIMA parameters, and d is the differencing term (allowed to take non-integer values). An R implementation of ARFIMA model can be done with 'arfima' function under the "forecast"package [61] . An ARFIMA(p, d, q) model is selected and estimated automatically using the Hyndman-Khandakar (2008) [59] algorithm to select p and q and the Haslett and Raftery (1989) [50] algorithm to estimate the parameters including d.
Exponential smoothing state space methods are very effective methods in case of time series forecasting. Exponential smoothing was proposed in the late 1950s [124] and has motivated some of the most successful forecasting methods. Forecasts produced using exponential smoothing methods are weighted averages of past observations, with the weights decaying exponentially as the observations get older. The ETS models belong to the family of state-space models, consisting of three-level components such as an error component (E), a trend component (T), and a seasonal component(S). This method is used to forecast univariate time series data. Each model consists of a measurement equation that describes the observed data, and some state equations that describe how the unobserved components or states (level, trend, seasonal) change over time [60] . Hence, these are referred to as state-space models. [59] . An R implementation of the model is available in the 'ets' function under "forecast" package [61] .
As an extension of autoregressive model, Self-exciting threshold autoregressive (SE-TAR) model is used to model time series data, in order to allow for higher degree of flexibility in model parameters through a regime switching behaviour [116] . Given a time-series data y t , the SETAR model is used to predict future values, assuming that the behavior of the time series changes once the series enters a different regime. This switch from one to another regime depends on the past values of the series. The model consists of k autoregressive (AR) parts, each for a different regime. The model is usually denoted as SETAR (k, p) where k is the number of threshold, there are k +1 number of regime in the model and p is the order of the autoregressive part. For example, suppose an AR(1) model is assumed in both regimes, then a 2-regime SETAR model is given by [41] :
where for the moment the ǫ t are assumed to be an i.i.d. white noise sequence conditional upon the history of the time series and c is the threshold value. The SETAR model assumes that the border between the two regimes is given by a specific value of the threshold variable y t−1 . The model can be implemented using 'setar' function under the "tsDyn" package in R [34] .
Bayesian Statistics has many applications in the field of statistical techniques such as regression, classification, clustering, and time series analysis. Scott and Varian [105] used structural time series models to show how Google search data can be used to improve short-term forecasts of economic time series. In the structural time series model, the observation in time t, y t is defined as follows:
where β t is the vector of latent variables, X t is the vector of model parameters, and ǫ t are assumed follow Normal distributions with zero mean and H t as the variance. In addition, β t is represented as follows:
where δ t are assumed to follow Normal distributions with zero mean and Q t as the variance. Gaussian distribution is selected as the prior of the BSTS model since we use the occurred frequency values ranging from 0 to ∞ [66] . An R implementation is available under the "bsts" package [103] , where one can add local linear trend and seasonal components as required. The state specification is passed as an argument to 'bsts' function, along with the data and the desired number of Markov chain Monte Carlo (MCMC) iterations, and the model is fit using an MCMC algorithm [104] .
The 'Theta method' or 'Theta model' is a univariate time series forecasting technique that performed particularly well in M3 forecasting competition and of interest to forecasters [11] . The method decomposes the original data into two or more lines, called theta lines, and extrapolates them using forecasting models. Finally, the predictions are combined to obtain the final forecasts. The theta lines can be estimated by simply modifying the 'curvatures' of the original time series [110] . This change is obtained from a coefficient, called θ coefficient, which is directly applied to the second differences of the time series:
where Y " data = Y t − 2Y t−1 + Y t−2 at time t for t = 3, 4, · · · , n and {Y 1 , Y 2 , · · · , Y n } denote the observed univariate time series. In practice, coefficient θ can be considered as a transformation parameter which creates a series of the same mean and slope with that of the original data but having different variances. Now, Eqn. (2) is a second-order difference equation and has solution of the following form [62] :
where a θ and b θ are constants and t = 1, 2, · · · , n. Thus, Y new (θ) is equivalent to a linear function of Y t with a linear trend added. The values of a θ and b θ are computed by minimizing the sum of squared differences:
Forecasts from the Theta model are obtained by a weighted average of forecasts of Y new (θ) for different values of θ. Also, the prediction intervals and likelihoodbased estimation of the parameters can be obtained based on a state-space model, demonstrated in [62] . An R implementation of the Theta model is possible with 'thetaf' function in "forecast" package [61] .
The main objective of TBATS model is to deal with complex seasonal patterns using exponential smoothing [33] . The name is acronyms for key features of the models: Trigonometric seasonality (T), Box-Cox Transformation (B), ARMA errors (A), Trend (T) and Seasonal (S) components. TBATS makes it easy for users to handle data with multiple seasonal patterns. This model is preferable when the seasonality changes over time [60] . TBATS models can be described as follows:
where y (µ) t is the time series at time point t (Box-Cox Transformed), s (i) t is the i-th seasonal component, l t is the local level, b t is the trend with damping, d t is the ARMA(p, q) process for residuals and e t as the Gaussian white noise. TBATS model can be implemented using 'tbats' function under the "forecast" package in R statistical software [61] .
Forecasting with artificial neural networks (ANN) has received increasing interest in various research and applied domains in the late 1990s. It has been given special attention in epidemiological forecasting [92] . Multi-layered feed-forward neural networks with back-propagation learning rules are the most widely used models with applications in classification and prediction problems [129] . There is a single hidden layer between the input and output layers in a simple feed-forward neural net, and where weights connect the layers. Denoting by ω ji the weights between the input layer and hidden layer and ν kj denotes the weights between the hidden and output layers. Based on the given inputs x i , the neuron's net input is calculated as the weighted sum of its inputs. The output layer of the neuron, y j , is based on a sigmoidal function indicating the magnitude of this net-input [43] . For the j th hidden neuron, the calculation for the net input and output are:
ω ji x i and y j = f (net h j ).
For the k th output neuron:
with λ ∈ (0, 1) is a parameter used to control the gradient of the function and J is the number of neurons in the hidden layer. The back-propagation [102] learning algorithm is the most commonly used technique in ANN. In the error back-propagation step, the weights in ANN are updated by minimizing
where, d pk is the desired output of neuron k and for input pattern p. The common formula for number of neurons in the hidden layer is h = (i+j) 2 + √ d, for selecting the number of hidden neurons, where i is the number of output y j and d denotes the number of i training patterns in the input x i [128] . The application of ANN for time series data is possible with 'mlp' function under "nnfor" package in R [72] .
Autoregressive neural network (ARNN) received attention in time series literature in late 1990s [37] . The architecture of a simple feedforward neural network can be described as a network of neurons arranged in input layer, hidden layer, and output layer in a prescribed order. Each layer passes the information to the next layer using weights that are obtained using a learning algorithm [128] . ARNN model is a modification to the simple ANN model especially designed for prediction problems of time series datasets [37] . ARNN model uses a pre-specified number of lagged values of the time series as inputs and number of hidden neurons in its architecture is also fixed [60] . ARNN(p, k) model uses p lagged inputs of the time series data in a one hidden layered feedforward neural network with k hidden units in the hidden layer. Let x denotes a p-lagged inputs and f is a neural network of the following architecture:
where c 0 , a j , w j are connecting weights, b j are p-dimensional weight vector and φ is a bounded nonlinear sigmoidal function (e.g., logistic squasher function or tangent hyperbolic activation function). These Weights are trained using a gradient descent backpropagation [102] . Standard ANN faces the dilemma to choose the number of hidden neurons in the hidden layer and optimal choice is unknown. But for ARNN model, we adopt the formula k = [(p+1)/2] for non-seasonal time series data where p is the number of lagged inputs in an autoregressive model [60] . ARNN model can be applied using the 'nnetar' function available in the R statistical package "forecast" [61] .
The idea of ensemble time series forecasts was given by Bates and Granger (1969) in their seminal work [12] . Forecasts generated from ARIMA, ETS, Theta, ARNN, WARIMA can be combined with equal weights, weights based on in-sample errors, or cross-validated weights. In the ensemble framework, cross-validation for time series data with user-supplied models and forecasting functions is also possible to evaluate model accuracy [106] . Combining several candidate models can hedge against an incorrect model specification. Bates and Granger(1969) [12] suggested such an approach and observed, somewhat surprisingly, that the combined forecast can even outperform the single best component forecast. While combination weights selected equally or proportionally to past model errors are possible approaches, many more sophisticated combination schemes, have been suggested. For example, rather than normalizing weights to sum to unity, unconstrained and even negative weights could be possible [45] . The simple equal-weights combination might appear woefully obsolete and probably non-competitive compared to the multitude of sophisticated combination approaches or advanced machine learning and neural network forecasting models, especially in the age of big data. However, such simple combinations can still be competitive, particularly for pandemic time series [106] . A flow diagram of the ensemble method is presented in Figure 2 . The ensemble method by [12] produces forecasts out to a horizon h by applying a weight w m to each m of the n model forecasts in the ensemble. The ensemble forecast f (i) for time horizon 1 ≤ i ≤ h and with individual component model forecasts f m (i) is then
The weights can be determined in several ways (for example, supplied by the user, set equally, determined by in-sample errors, or determined by cross-validation). The "forecastHybrid" package in R includes these component models in order to enhance the "forecast" package base models with easy ensembling (e.g., 'hybridModel' function in R statistical software) [107] .
The idea of hybridizing time series models and combining different forecasts was first introduced by Zhang [127] and further extended by [68, 24, 26, 23] . The hybrid forecasting models are based on an error re-modeling approach, and there are broadly two types of error calculations popular in the literature, which are given below [84, 30] :
In the additive error model, the forecaster treats the expert's estimate as a variable,Ŷ t , and thinks of it as the sum of two terms:
where Y t is the true value and e t be the additive error term.
In the multiplicative error model, the forecaster treats the expert's estimateŶ t as the product of two terms:
where Y t is the true value and e t be the multiplicative error term. Now, even if the relationship is of product type, in the log-log scale it becomes additive. Hence, without loss of generality, we may assume the relationship to be additive and expect errors (additive) of a forecasting model to be random shocks [23] . These hybrid models are useful for complex correlation structures where less amount of knowledge is available about the data generating process. A simple example is the daily confirmed cases of the COVID-19 cases for various countries where very little is known about the structural properties of the current pandemic. The mathematical formulation of the proposed hybrid model (Z t ) is as follows:
where L t is the linear part and N t is the nonlinear part of the hybrid model. We can estimate both L t and N t from the available time series data. LetL t be the forecast value of the linear model (e.g., ARIMA) at time t and ǫ t represent the error residuals at time t, obtained from the linear model. Then, we write
These left-out residuals are further modeled by a nonlinear model (e.g., ANN or ARNN) and can be represented as follows:
where f is a nonlinear function, and the modeling is done by the nonlinear ANN or ARNN model as defined in Eqn. (5) and ε t is supposed to be the random shocks. Therefore, the combined forecast can be obtained as follows:
whereN t is the forecasted value of the nonlinear time series model. An overall flow diagram of the proposed hybrid model is given in Figure 3 . In the hybrid model, a nonlinear model is applied in the second stage to re-model the left-over autocorrelations in the residuals, which the linear model could not model. Thus, this can be considered as an error re-modeling approach. This is important because due to model misspecification and disturbances in the pandemic rate time series, the linear models may fail to generate white noise behavior for the forecast residuals. Thus, hybrid approaches eventually can improve the predictions for the epidemiological forecasting problems, as shown in [24, 26, 23] . These hybrid models only assume that the linear and nonlinear components of the epidemic time series can be separated individually. The implementation of the hybrid models used in this study are available in [1] . 
Five time series COVID-19 datasets for the USA, India, Russia, Brazil, and Peru UK are considered for assessing twenty forecasting models (individual, ensemble, and hybrid). The datasets are mostly nonlinear, nonstationary, and non-gaussian in nature. We have used root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE), and symmetric MAPE (SMAPE) to evaluate the predictive performance of the models used in this study. Since the number of data points in both the datasets is limited, advanced deep learning techniques will over-fit the datasets [51] .
We use publicly available datasets to compare various forecasting frameworks. COVID-19 cases of five countries with the highest number of cases were collected [2, 3] . The datasets and their description is presented in Table 2 . 
Characteristics of these five time series were examined using Hurst exponent, KPSS test and Terasvirta test and other measures as described in Section 3. Hurst exponent (denoted by H), which ranges between zero to one, is calculated to measure the longrange dependency in a time series and provides a measure of long-term nonlinearity. For values of H near zero, the time series under consideration is mean-reverting. An increase in the value will be followed by a decrease in the series and vice versa. When H is close to 0.5, the series has no autocorrelation with past values. These types of series are often called Brownian motion. When H is near one, an increase or decrease in the value is most likely to be followed by a similar movement in the future. All the five COVID-19 datasets in this study possess the Hurst exponent value near one, which indicates that these time series datasets have a strong trend of increase followed by an increase or decrease followed by another decline.
KPSS tests are performed to examine the stationarity of a given time series. The null hypothesis for the KPSS test is that the time series is stationary. Thus, the series is nonstationary when the p-value less than a threshold. From Table 3 , all the five datasets can be characterized as non-stationary as the p-value < 0.01 in each instances. Terasvirta test examines the linearity of a time series against the alternative that a nonlinear process has generated the series. It is observed that the USA, Russia, and Peru COVID-19 datasets are likely to follow a nonlinear trend. On the other hand, India and Brazil datasets have some linear trends. Further, we examine serial correlation, skewness, kurtosis, and maximum Lyapunov exponent for the five COVID-19 datasets. The results are reported in Table 4 . The serial correlation of the datasets is computed using the Box-Pierce test statistic for the null hypothesis of independence in a given time series. The p-values related to each of the datasets were found to be below the significant level (see Table 4 ). This indicates that these COVID-19 datasets have no serial correlation when lag equals one. Skewness for Russia COVID-19 dataset is found to be negative, whereas the other four datasets are positively skewed. This means for the Russia dataset; the left tail is heavier than the right tail. For the other four datasets, the right tail is heavier than the left tail. The Kurtosis values for the India dataset are found positive while the other four datasets have negative kurtosis values. Therefore, the COVID-19 dataset of India tends to have a peaked distribution, and the other four datasets may have a flat distribution. We observe that each of the five datasets is non-chaotic in nature, i.e., the maximum Lyapunov exponents are less than unity. A summary of the implementation tools is presented in Table 5 . 
We used four popular accuracy metrics to evaluate the performance of different time series forecasting models. The expressions of these metrics are given below.
where y i are actual series values,ŷ i are the predictions by different models and n represent the number of data points of the time series. The models with least accuracy metrics is the best forecasting model.
This subsection is devoted to the experimental analysis of confirmed COVID-19 cases using different time series forecasting models. The test period is chosen to be 15 days and 30 days, whereas the rest of the data is used as training data (see Table 2 ). In first columns of Tables 6 and 7, we present training data and test data for USA, India, Brazil, Russia and Peru. The autocorrelation function (ACF) and partial autocorrelation function (PACF) plots are also depicted for the training period of each of the five countries in Tables 6 and 7 . ACF and PACF plots are generated after applying the required number of differencing of each training data using the r function 'diff'. The required order of differencing is obtained by using the R function 'ndiffs' which estimate the number of differences required to make a given time series stationary. The integer-valued order of differencing is then used as the value of 'd' in the ARIMA(p, d, q) model. Other two parameters 'p' and 'q' of the model are obtained from ACF and PACF plots respectively (see Tables 6  and 7) . However, we choose the 'best' fitted ARIMA model using AIC value for each training dataset. Table 6 presents the training data (black colored) and test data (red-colored) and corresponding ACF and PACF plots for the five time-series datasets. Further, we checked twenty different forecasting models as competitors for the short-term forecasting of COVID-19 confirmed cases in five countries. 15-days and 30-days ahead forecasts were generated for each model, and accuracy metrics were computed to determine the best predictive models. From the ten popular single models, we choose the best one based on the accuracy metrics. On the other hand, one hybrid/ensemble model is selected from the rest of the ten models. The bestfitted ARIMA parameters, ETS, ARNN, and ARFIMA models for each country are reported in the respective tables. Table 7 presents the training data (black colored) and test data (red-colored) and corresponding plots for the five datasets. Twenty forecasting models are implemented on these pandemic time-series datasets. Table 5 gives the essential details about the functions and packages required for implementation. Results for USA COVID-19 data: Among the single models, ARIMA (2, 1, 4) performs best in terms of accuracy metrics for 15-days ahead forecasts. TBATS and ARNN (16, 8) also have competitive accuracy metrics. Hybrid ARIMA-ARNN model improves the earlier ARIMA forecasts and has the best accuracy among all hybrid/ensemble models (see Table 8 ). Hybrid ARIMA-WARIMA also does a good job and improves ARIMA model forecasts. In-sample and out-of-sample forecasts obtained from ARIMA and hybrid ARIMA-ARNN models are depicted in Fig. 4(a) .
Out-of-sample forecasts are generated using the whole dataset as training data. ARFIMA(2,0,0) is found to have the best accuracy metrics for 30-days ahead forecasts among single forecasting models. BSTS and SETAR also have good agreement with the test data in terms of accuracy metrics. Hybrid ARIMA-WARIMA model and has the best accuracy among all hybrid/ensemble models (see Table 9 ). In-sample and out-of-sample forecasts obtained from ARFIMA and hybrid ARIMA-WARIMA models are depicted in Fig. 4(b) . Results for India COVID-19 data: Among the single models, ANN performs best in terms of accuracy metrics for 15-days ahead forecasts. ARIMA(1,2,5) also has competitive accuracy metrics in the test period. Hybrid ARIMA-ARNN model improves the ARIMA(1,2,5) forecasts and has the best accuracy among all hybrid/ensemble models (see Table 10 ). Hybrid ARIMA-ANN and hybrid ARIMA-WARIMA also do a good job and improves ARIMA model forecasts. In-sample and out-of-sample forecasts obtained from ANN and hybrid ARIMA-ARNN models are depicted in Fig. 5(a) . Out-of-sample forecasts are generated using the whole dataset as training data (see Fig. 5 ). ANN is found to have the best accuracy metrics for 30-days ahead forecasts among single forecasting models for India COVID-19 data. Ensemble ANN-ARNN-WARIMA model and has the best accuracy among all hybrid/ensemble models (see Table 11 ). In-sample and out-of-sample forecasts obtained from ANN and ensemble ANN-ARNN-WARIMA models are depicted in Fig. 5(b) . Results for Brazil COVID-19 data: Among the single models, SETAR performs best in terms of accuracy metrics for 15-days ahead forecasts. Ensemble ETS-Theta-ARNN (EFN) model has the best accuracy among all hybrid/ensemble models (see Table 12 ). In-sample and out-of-sample forecasts obtained from SETAR and ensemble EFN models are depicted in Fig. 6 (a). WARIMA is found to have the best accuracy metrics for 30-days ahead forecasts among single forecasting models for Brazil COVID-19 data. Hybrid WARIMA-ANN model has the best accuracy among all hybrid/ensemble models (see Table 13 ). Insample and out-of-sample forecasts obtained from WARIMA and hybrid WARIMA-ANN models are depicted in Fig. 6(b) . Results for Russia COVID-19 data: BSTS performs best in terms of accuracy metrics for a 15-days ahead forecast in the case of Russia COVID-19 data among single models. Theta and ARNN(3,2) also show competitive accuracy measures. Ensemble ETS-Theta-ARNN (EFN) model has the best accuracy among all hybrid/ensemble models (see Table 14 ). Ensemble ARIMA-ETS-ARNN and ensemble ARIMA-Theta-ARNN also performs well in the test period. In-sample and out-ofsample forecasts obtained from BSTS and ensemble EFN models are depicted in Fig. 7(a) . SETAR is found to have the best accuracy metrics for 30-days ahead forecasts among single forecasting models for Russia COVID-19 data. Ensemble ARIMA-Theta-ARNN (AFN) model has the best accuracy among all hybrid/ensemble models (see Table 15 ). All five ensemble models show promising results for this dataset.
In-sample and out-of-sample forecasts obtained from SETAR and ensemble AFN models are depicted in Fig. 7(b) . Results for Peru COVID-19 data: WARIMA and ARFIMA(2,0.09,1) perform better than other single models for 15-days ahead forecasts in Peru. Hybrid WARIMA-ARNN model improves the WARIMA forecasts and has the best accuracy among all hybrid/ensemble models (see Table 16 ). In-sample and out-of-sample forecasts obtained from WARIMA and hybrid WARIMA-ARNN models are depicted in Fig. 8(a) . ARFIMA(2,0,0) and ANN depict competitive accuracy metrics for 30-days ahead forecasts among single forecasting models for Peru COVID-19 data. Ensemble ANN-ARNN-WARIMA (AAW) model has the best accuracy among all hybrid/ensemble models (see Table 17 ). In-sample and out-of-sample forecasts obtained from ARFIMA(2,0,0) and ensemble AAW models are depicted in Fig. 8(b) . Results from all the five datasets reveal that none of the forecasting models performs uniformly, and therefore, one should be carefully select the appropriate forecasting model while dealing with COVID-19 datasets.
In this study, we assessed several statistical, machine learning, and composite models on the confirmed cases of COVID-19 data sets for the five countries with the highest number of cases. Thus, COVID-19 cases in the USA, followed by India, Brazil, Russia, and Peru, are considered. The datasets mostly exhibit nonlinear and nonstationary behavior. Twenty forecasting models were applied to five datasets, and an empirical study is presented here. The empirical findings suggest no universal method exists that can outperform every other model for all the datasets in COVID-19 nowcasting. Still, the future forecasts obtained from models with the best accuracy will be useful in decision and policy makings for government officials and policymakers to allocate adequate health care resources for the coming days in responding to the crisis. However, we recommend updating the datasets regularly and comparing the accuracy metrics to obtain the best model. As this is evident from this empirical study that no model can perform consistently as the best forecasting model, one must update the datasets regularly to generate useful forecasts. Time series of epidemics can oscillate heavily due to various epidemiological factors, and these fluctuations are challenging to be captured adequately for precise forecasting.
All five different countries, except Brazil and Peru, will face a diminishing trend in the number of new confirmed cases of COVID-19 pandemic. Followed by the both short-term out of sample forecasts reported in this study, the lockdown and shutdown periods can be adjusted accordingly to handle the uncertain and vulnerable situations of the COVID-19 pandemic. Authorities and health care can modify their planning in stockpile and hospital-beds, depending on these COVID-19 pandemic forecasts.
Models are constrained by what we know and what we assume but used appropriately, and with an understanding of these limitations, they can and should help guide us through this pandemic. Since purely statistical approaches do not account for how transmission occurs, they are generally not well suited for long-term predictions about epidemiological dynamics (such as when the peak will occur and whether resurgence will happen) or inference about intervention efficacy. Several forecasting models, therefore, limit their projections to two weeks or a month ahead.
In this research, we have focused on analyzing the nature of the COVID-19 time series data and understanding the data characteristics of the time series. This empirical work studied a wide range of statistical forecasting methods and machine learning algorithms. We have also presented more systematic representations of single, ensemble, and hybrid approaches available for epidemic forecasting. This quantitative study could be used to assess and forecast COVID-19 confirmed cases, which will benefit epidemiologists and modelers in their real-world applications.
Considering this scope of the study, we can present a list of challenges of pandemic forecasting (short-term) with the forecasting tools presented in this chapter:
-Collect more data on the factors that contribute to daily confirmed cases of COVID-19. -model the entire predictive distribution, with particular focus on accurately quantifying uncertainty [55] . -There is no universal model that can generate 'best' short-term forecasts of COVID-19 confirmed cases. -Continuously monitor the performance of any model against real data and either re-adjust or discard models based on accruing evidence. -Developing models in real-time for a novel virus, with poor quality data, is a formidable task and real challenge for epidemic forecasters. -Epidemiological estimates and compartmental models can be useful for longterm pandemic trajectory prediction, but they often assume some unrealistic assumptions [64] . -Future research is needed to collect, clean, and curate data and develop a coherent approach to evaluate the suitability of models with regard to COVID-19 predictions and forecast uncertainties.
For the sake of repeatability and reproducibility of this study, all codes and data sets are made available at https://github.com/indrajitg-r/Forecasting-COVID-19-cases.