Abstract
In this article, we analyze the growth pattern of Covid-19 pandemic in India from March 4th to May 15th using regression analysis (exponential and polynomial), auto-regressive integrated moving averages (ARIMA) model as well as Exponential Smoothing and Holt-Winters models. We found that the growth of Covid-19 cases follows a power regime of (t2, t,…) after the exponential growth. We have found the optimal change points from where the covid-19 cases shifts their course of growth from exponential to quadratic and then quadratic to linear. We have also found the best fitted regression models using the various criteria like-significant p-values, coefficients of determination R2 values and ANOVA etc. Further, we have searched the best fitting ARIMA model for the data using the AIC (Akaike Information Criterion) and CAIC (Consistent Akaike Information Criterion) and forecasted the number of cases for future days. We have used the usual exponential smoothing and Holt-Winters models for the data. We further found that the ARIMA(2,2,0) model is the best-fitting model for Covid-19 cases in India.
1. Introduction
The Covid-19 pandemic has created a lot of havoc in the world. It is caused by a virus called SARS-CoV-2, which comes from the family of coronaviruses and is believed to be originated from the unhygienic wet seafood market in Wuhan, China but it has now infected more than 200 countries of the world. With around 5 million people affected around the world (As of 20th May, 2020, *to be updated*), it has forced people to stay in their homes and has caused huge devastation in the world economy. (Ref [25], [28]).
In India, the first case of Covid-19 was reported on 30th January, which is linked to the Wuhan city of China (as the patient has travel history to the city). On 4th March, India saw a sudden hike in the number of cases and since then, the numbers are increasing day by day. As of 19th May (*to be updated*), India has more than 106,000 cases with more than 3300 deaths. (Ref [27]).
Since the outbreak of the pandemic, scientists across the world have been indulged in the studies regarding the spread of the virus. Lin et al. [1] suggested the use of the SEIR (Susceptible-Exposed-Infectious-Removed) model for the spread in China and studied the importance of government-implemented restrictions on containing the infection. As the disease grew further, Ivorra et al. [2] suggested a θ-SEIHRD model that took into account various special features of the disease. It also included asymptomatic cases into account (around 51%) in order to forecast the total cases in China (around 168500). Giordano et al. [3] also suggested an extended SIR model called SIDHARTHE model for cases in Italy which was more customized for Covid-19 in order to effectively model the course of the pandemic to help plan a better control strategy. Apart from the epidemiological models, various time-series models were also suggested in order to model the cases and predict future cases. Ceylan [4] suggested the use of Auto-Regressive Integrated Moving Average (ARIMA) model to develop and predict the epidemiological trend of Covid-19 for better allocation of resources and proper containment of the virus in Italy, Spain and France. Chintalapudi et al. [5] suggested its use for predicting the number of cases and deaths post 60-days lockdown in Italy. Petropoulos and Makridakis [6] suggested the use of exponential smoothing method to model the trend of the virus, globally. Kumar et al. [29] gave a review on the various aspects of modern technology used to fight against COVID-19 crisis.
Since the spread of the virus started to grow in India, various measures were taken by the Indian Government in order to contain it. A nationwide lockdown was announced on March 25th to April 14th, which was later extended to May 3rd. The whole country was divided into containment zones (where large number of cases were observed from a relatively smaller region), red zones (districts where risk of transmission was high and had higher doubling rates), green zones (districts with no confirmed case from last 21 days) and orange zones (which didn’t fall into the above three zones). After the further extension of the lockdown till May 17th, various economic activities were allowed to start (with high surveillance) in areas of less transmission. Further, the lockdown is now extended to May 31st and some more economic activities have been allowed as per the transmission rates, which are the rates at which infectious cases cause new cases in the population, i.e. the rate of spread of the disease.
On the other hand, Indian scientists and researchers are also working on addressing the issues arising from the pandemic, including production of PPE kits and tests kits as well as studying the behaviour of spread of the disease and other aspects of management. Various mathematical and statistical methods have been used for predicting the possible spread of Covid-19. The classical SIR model by Ranjan [7], SEIR model by Ranjan [8] and SIQR model by Tiwari [9] suggested the increasing trend of the virus and predicted the peaks of the pandemic. Early researches showed the pandemic to reach its peak by mid-May. They also showed that the basic reproduction number R0 and the doubling rates are lower in India, with comparison to European nations and USA. A tree-based model was proposed by M.K. & Bhatnagar [10] and Bhatnagar [11] in order to study and predict the trends. They suggest that lockdown and social-distancing in India has played a significant role to control the infection rates. Kumar et al. [12] suggested the use of ARIMA and Richard models for the same and have used them to predict the cases and deaths in future days. Chatterjee et al. [13] suggests growth of the pandemic through power law and its saturation at the later stages. Due to the complexities in the epidemic models of Covid-19, various researchers have been focusing on the data in order to forecast the future cases. Chatterjee et al. [13,14] and Ziff & Ziff [15] suggest that after exponential growth, the total count follows a power regime of t3, t2, t and t0.5 before flattening out, where t refers to time. It can therefore be realized that there is an urgent need to model and forecast the growth of Covid-19 in India as the virus is in the growing stage here.
In India, the most affected states are Maharashtra with over 37,000 cases (as of 19th May,2020 *need to be updated*), Tamil Nadu (over 12,000 cases) and Gujarat (over 12,000 cases). The greatest number of cases per million have been seen in the national capital of Delhi (532.65 cases per million). (Refer [20] for population estimates). Various states such as Arunachal Pradesh, Goa, Mizoram and Manipur have been declared Covid-19 free states as they have treated all their cases since more than 14 days. States of Nagaland and Sikkim and Union Territories of Lakshwadeep Islands and Daman and Diu are yet to report a single case. These large variations suggest the effectiveness of lockdown and sealing of state borders in containing the virus. In the latest research. Singh & Jadaun [25] studied the significance of lockdown in India and suggested that the new Covid-19 cases would stop by the end of August in India with around 350,000 total cases. While some states may see an early stopping of new cases such as Telangana (mid-June), Uttar Pradesh and West Bengal (July-end) etc., the badly affected states of Maharashtra, Tamil Nadu and Gujarat will achieve this by August-end.
Since a proven vaccine and medication is yet to be developed by the researchers then in such a scenario, modelling the present situation and forecasting the future outcome becomes crucially important in order to utilize our resources in the most optimal way. Therefore, the article aims to study the growth curve of Covid-19 cases in India and forecast its future observations. Since the disease is still in its growing age and very dynamic in nature, no model remains perfectly valid for future. We need to develop the understanding of the present situation of the pandemic.
In this article, we first study the growth curve using regression methods (exponential, linear and polynomial etc.) and propose an optimal model for fitting the cases till May 15th. Further, we propose the use of time-series models for forecasting the future observations on Covid-119 cases. Here we reach the best-fitted ARIMA model for forecasting the Covid-19 cases. We also compare these results with Exponential Smoothing (Holt-Winters) model. This study will help us to understand the course of spread of SARS-CoV-2 in India better and help the government and the people to optimally use the resources available to them.
2. Statistical Methodologies
In this section, we briefly present the statistical techniques used for analyzing the Covid-19 cases in India. Here, we used usual regression (exponential, polynomial), times series (ARIMA and GARCH) and exponential smoothing models.
2.1 Exponential-Polynomial Regression
Regression is a statistical technique that attempts to estimate the strength and nature of relationship between a dependent variable and a series of independent variables. Regression analyses may be linear and non-Linear. A regression is called linear when it is linear in parameters e.g. y = β0 + β1t+∈ and y = β0 + β1x + β2x2 + β3x3 + ∈, ∈ ~N(0, σ2), where is response variable, t denotes the indepenet variable, β0 is the intercept and other βs are known as slop.
A non-linear regression is a regression when it is non-linear in its parameters e.g. . In the beginning of the spread of a disease, we see that the new cases are directly proportional to the existing infected cases and may be represented by
, where k is the proportionality constant. Solving this differential equation, we get that, at the beginning of a pandemic,
![Embedded Image](https://www.medrxiv.org/sites/default/files/highwire/medrxiv/early/2020/05/25/2020.05.20.20107540/embed/graphic-1.gif)
Thus, at the beginning of a disease, the growth curve of the cases grows exponentially.
As the disease spreads in a region, governments start to take action and people start becoming conscious about the disease. Thus, after some time, the disease starts to follow a polynomial growth rather than continuing to grow exponentially.
In order to fit an exponential regression to our data, we linearize the equation by taking the natural logarithm of the equation and convert it to a linear regression in first order.
We estimate the parameters of a linear regression of order n as following-
Let the model of linear regression of order p be: with ∈ ~N(0, σ2) and i = 1, 2,..,N. Let
represent the residual sum of square (RSS).
By minimizing the RSS, we get the best estimates of these coefficients by solving the following normal equations; ,
,…,
. This technique is referred to as the ordinary least squares (OLS). We will use this technique of the OLS in order to estimate the coefficients of our proposed model. (Refer Montgomery et al. [17])
Since we know that the growth curve of the disease changes after some time point, exponential to polynomial, we propose to use the following joint regression model with change point μ,
where we take
, f2 = β0 + β1t + β1t2 +… + βp tp + ∈, ∈ ~N(0, σ2) and p is the order of the polynomial regression model and stands for the time (an independent variable).
During the analysis, we found that a suitable choice of f2(t) is a quadratic or cubic model. One the order of the polynomial is kept fixed, as optimum values of the change point can be easily obtained. We can obtain the OLS of the parameters of the model (1) as given below.
The least square estimates (LSEs) of the parameters, Θ = {θ1, θ2, μ, β0, β1, β2, β3} can be obtained by minimizing the residual sum of squares as given by
where,
and
are the estimates value of yi from the exponential and polynomial regression models, respectively and N is the size of the data set.
The LSEs of Θ = {θ1, θ2, μ, β0, β1, β2, β3} can be obtained as the simultaneous solutin of the following normal equations, ,
,
,
,
,
, and
. Solution to these equations is difficult since the parameter μ fi is decenter time point. We suggest to use the following algorithm while μ is kept fixed.
Algorithm 1
Set μ = j; j = 1,2,…, N
For a given μ, obtain LSEs of θ1, θ2 using the data {(y1, t1), (y2, t2),…, (yj, tj)}
For a given μ, obtain LSEs of β0, β1, β2, β3 using the data {(yj+1, tj+1), (yj+2, tj+2),…, (yN, tN)}
Compute RSSj using the estimates of {θ1, θ2, β0, β1, β2, β3} and fixed μ
Repeat steps 2-4 for all j = 1,2,…, N and obtain the RSS for each iteration.
Search j* which is a j that corresponds to the minimum RSS.
Take μ = j* as an optimum value of the parameter μ
In order to find the optimal value of μ, i.e. the turning point between the exponential and polynomial growth, we will use the technique of minimizing the residual sum squares in section 3.
2.2 ARIMA Model
The Auto-Regressive Integrated Moving Averages method gauges the strength of one dependent variable relative to other changing variables. An ARIMA model can be understood by outlining each of its components as follows: Autoregression (AR) refers to a model that shows a changing variable that regresses on its own lagged, or prior values. Integrated (I) represents the differencing of raw observations to allow for the time series to become stationary, i.e., data values are replaced by the difference between the data values and the previous values. Moving average (MA) incorporates the dependency between an observation and a residual error from a moving average model applied to lagged observations. Here, we may refer the reader to follow Box et al. [18] and Box et al. [19] for more details on ARIMA model, estimation and its application
ARIMA is a joint model of two models, autoregressive ARp and moving average MAq and is integrated using the difference variable, d. The mathematical model for ARIMA can be expressed as follows; Let us consider a backshift operator B, which causes the back shifting of the observation by 1 interval. For any time-series {yt}, we have: Byt = yt−1, B2yt = yt−1,…, Bnyt = yt−n.
The general model for ARIMA/Seasonal ARIMA can be expressed as:
where, øP (Bs) = 1 − ø1 (Bs) − ø2 (B2s) −… − øP (BPs), φp (B) = 1 − φ1 (B) − φ2 (B2) −… − φp (Bp), θq (B) = 1 −θ1 (B) − φ2(B2) −… − φp (Bq), ϑQ (Bs) = 1 − ϑ1(Bs) − ϑ2(B2s) −… − ϑP (BPs).
Since we are using ARIMA only, the components of SARIMA, i.e. P=D=Q=s=0 for our model.
So, our model becomes:
where the parameters (p, d and q) represents the lag order, degree of differencing and the order of moving average, respectively.
We obtain the optimal values of p, d and q by using the AIC (Akaike Information Criterion) and CAIC (Consistent Akaike Information Criterion), for more details see https://en.wikipedia.org/wiki/Akaike_information_criterion. These information criteria may be used for selecting the best fitting models. Lower the values of criteria, higher will be its relative quality. The AIC, CAIC and BIC are given by
where K=number of model parameters, ℓ = maximized value of log − likelihood function and N=no. of data points.
2.3 Exponential Smoothing
Exponential smoothing is one of the simple techniques to model time-series data where the past observations are assigned weights that are exponentially decreasing over time. We propose the following models, for modelling of Covid-19 cases (see Holt [23] and Winters [24]).
For single exponential smoothing, let the raw observations be denoted by {yt} and {st} denote the best estimate of trend at time t. Then, s0 = y0, st = αyt + (1 − α) (st−1), where α ∈ (0, 1) denotes the data smoothing factor.
For double exponential (Holt-Winters) smoothing, let the raw observations be denoted by {yt}, smoothened values by {st} and {bt} denotes the best estimate of trend at time t. Then,
where α ∈ (0, 1) denotes the data smoothing factor and β ∈ (0, 1) denotes the trend smoothing factor. For the forecast beyond N days is represented by, we denote it by FN+m and calculated by
![Embedded Image](https://www.medrxiv.org/sites/default/files/highwire/medrxiv/early/2020/05/25/2020.05.20.20107540/embed/graphic-8.gif)
3. Analysis of Covid-19 cases in India
For this study, we have used the data available at GitHub, provided by John Hopkins University (see [16]). We have used the data from March 4th to May 15th (*to be updated later*) for continuity of the data.
3.1 Exponential-Polynomial Regressions
We know that at the beginning of the spread of the disease in India, the growth was exponential and after some time, it was shifted to polynomial. We first obtain optimum turning point of the growth, i.e. when did the growth rate of the disease shifted to polynomial regime from the exponential. We consider both quadratic and cubic regression model for second part of the data. We will also discuss the types of polynomial growth (with their equations) in India.
In order to find the turning point of the growth curve, we follow the Algorithm 1, given in the previous section. Using that, we evaluate the RSS for all the days (from March 4th) and find the date on which it is minimum. The change points of growth curve for cubic and quadratic regressions are presented in Figure 1 depending upon the size of the data set. From Figure 1, we can confirm that the growth rate of Covid-19 cases was exponential till April 5th and then after it follows the polynomial growth regime while we use the covid-19 cases till May 2nd.
Turning point of growth curve for cubic and quadratic regression beyond change point using the Covid-19 cases from 4th March to a given day.
Trend of RSS and optimum μ for exponential-quadratic regression model.
We call the region of exponential growth in India as Region I. The coefficients of the model are presented in Table 2.
Regression Table for Region I (Exponential Regression)
We see that after the exponential regime (till April 5th), the growth curve follows a polynomial growth till May 2nd. After this, we again see a change in the behavior of the growth curve. In Tables 3-4, we try to model these growth curves through regression analysis.
Regression models fitting for Region II (5 April – 2 May).
Regression models fitting for Region III (3rd May – 15 May).
Having evaluated the coefficients for various models (i.e. linear, quadratic and cubic) as well as the important statistics (i.e. R2 values, p-values of the models as well as individual coefficients and F-statistic), we will select the best fitting models. In order to select the best fitting models for Region II (April 6th to May 2nd) and III (May 3rd to May 15th), we have the following steps. We select that model which has high R2 values, significant p-value, high F-statistic and where the p-values of all the variables are significant.
We see for Region II, from Table 3 that the linear model is having a relatively lower F-statistic and R2 values in comparison to the Quadratic and Cubic models. So, we eliminate the possibility of linear fitting. Further, we see that the p-values, F-statistics and the R2 values are quite significant in both Quadratic as well as the Cubic models. But, if we look at the individual p-values of the coefficients, we see that the individual p-values are not significant for the Cubic model. On the other hand, the individual p-values are significant for the Quadratic model. Thus, we can conclude that the Quadratic model is the best fitting model for Region II (April 6th to May 2nd).
For Region III, from Table 4, that all the three models have high F-statistic values, high p-values and high R2 values. But we notice that the coefficient individual p-values are not significant in both Quadratic and Cubic models. Thus, we conclude that the Linear model is the best fitting model for Region III (May 3rd to May 15th).
ANOVA Table for Region II (Quadratic Regression)
Both the ANOVA tables for Region II and III suggest significant p-values for its coefficients and suggest that the models fit well the respective regions.
Thus, according to our study, the growth of the virus was exponentially increasing from March 4th to April 5th. Then after, the virus grew by following a quadratic rate from April 6th to May 2nd. Since May 3rd, we have been experiencing a linear growth, see Table 6 for best fitted regression models. Figure 2 shows the best fitted Regression models to the daily cumulative cases of Covid-19 in India till 15 May.
Course of Covid-19 growth in India
Fitted Regression models to the daily cumulative cases of Covid-19 in India till 15 May (Bar chart shows the daily confirmed cases).
3.2 Time series Models fitting
First, we obtain the optimal ARIMA parameters (p, d, q) by using the AIC and CAIC Criteria. We take various possible combinations of (p, d, q) and compute the AIC and CAIC Criteria. Then, select the best fitted ARIMA model that has the lowest AIC and CAIC among all considered models. According to the AIC and CAIC, the ARIMA (2,2,0) is the best fitted models for the Covid-19 cases, India. Estimates of ARIMA (2,2,0) parameters and MAPE are shown in Table 8.
AIC and CAIC for ARIMA models for Covid-19 cases, India
Estimates of ARIMA (2,2,0) parameters and MAPE.
Estimates of the Holt-Winters exponential smoothing and exponential smoothing models are given in Table 10. According to the MAPE and accuracy measures, the ARIMA (2,2,0) is a better model than the Holt-Winters exponential smoothing and exponential smoothing models. From this, we can conclude that the ARIMA model is the best fit for the cases of Covid-19, followed by Holt-Winters model. The Forecasting are shown in Table 12 and Figure 3. We observe that the ARIMA model captures the trend well but it underestimates the actual covid-19 cases. We therefore suggest to update the ARIMA model in future or to use some generalized versions of the ARIMA models in future studies.
Estimates and MAPE of exponential smoothing models.
Forecast using ARIMA and Holt-Winters models for 10 days
Fitted ARIMA (2, 2, 0) and exponential smoothing models and forecasting from ARIMA for Covid-19 cases in India (stars show the actual observations).
4. Conclusions
From the regression analysis, we conclude that the spread of convi-19 disease grew exponentially from March 3rd to April 5th. Further, from April 6th to May 2nd, the cases followed a quadratic regression. From May 3rd to May 15th, we see a linear growth of the virus with average daily cases of 3584.
Verma et al. [14] showed the four stages of the epidemic, S1: exponential, S2: power law, S3: linear and S4: flat. Therefore, covid-19 pandemic in India has entered in stage S3 of linear growth. In the days to come, it is highly likely that the total cases may start to follow a square root equation, i.e. . And this may lead to reduction in the daily number of cases
leading to flattening of the curve.
In time series analysis, we conclude that the ARIMA (2, 2, 0) is the best fitting model for the cases of Covid-19 with an accuracy of 96.13%. This is also suggested by the Exponential Smoothing and Holt-Winters models. The basic exponential smoothing is not very accurate for our case but we see that the Holt-Winters model is around 96.073% accurate. Both ARIMA (2, 2, 0) and Holt-Winters models suggest a rise in the number of cases in the coming days. We also observed that the ARIMA model underestimates the actual observations. Therefore we suggest updating the ARIMA model time to time or using some generalized ARIMA models.
We may also conclude that the cases of Covid-19 will rise in the coming days and but slowly, we may head towards the reduction in the daily number of cases. But this should be accompanied by following of proper safety measures and following the guidelines of the government of India. With the gradual relaxation of lockdown measures, if proper precautions are not taken, we may see an increase in the daily cases. We must learn to lead our lives by following all the precautions once the lockdown measures are relaxed.
Data Availability
The data is available at GitHub, provided by John Hopkins University. 15. https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
Acknowledgments
Dr. Vikas Kumar Sharma greatly acknowledges the financial support from Science and Engineering Research Board, Department of Science & Technology, Govt. of India, under the scheme Early Career Research Award (file no.: ECR/2017/002416).