Case definitions were brought forth by WHO and different national health authorities and a confirmed case was defined as a case with respiratory specimens (commonly broncho-alveolar lavage) which tested positive for the 2019-nCoV in at least one of the following three methods: (a) isolation of 2019-nCoV, (b) at least two positive results by real-time, reverse-transcriptionpolymerasechain-reaction (RT-PCR) assay for 2019-nCoV, or (c) a genetic sequence that matches 2019-nCoV.
2019-nCoV closely resembles the SARS-CoV (about 80% sequence identity, and uses the same cell entry receptor -ACE-II), MERS-CoV and a bat coronavirus (with which it has 96% sequence identity) which are enveloped nonsegmented positive sense RNA viruses from the family Coronaviridae broadly distributed in mammals, and is likely to behave clinically and epidemiologically akin to SARS. Initial transmission of the disease involved seafood to human (zoonotic) transmission, which rapidly evolved into human to human transmission in secondary and tertiary cases. Human to human transmission of 2019-nCoV occurs through large droplets and contact and less so by means of fomites and aerosols. It has also been seen to survive for days on surfaces like cardboard, steel and plastic. described the early transmission dynamics of the disease among initially affected Chinese population. Their preliminary estimate of the incubation period distribution indicates a 14-day medical observation period or quarantine for exposed persons. The rapidity with which the virus has been noted to mutate, as evidenced by ability of human to human transmission, it is unlikely that a specific effective vaccine can be developed soon. Thus epidemiological interventions like lock-down or mass screening and quarantine might be more effective to contain the spread.
There have already been a lot of studies on this recent epidemic. They cover a wide range of aspects. Huang et al. (2020) described the spectrum of clinical presentation of 41 Chinese patients from Wuhan. The commonest symptoms were fever (98%), cough (76%) and myalgia (44%), while sputum production, headache and diarrhoea were rare presentations. Acute respiratory distress syndrome (ARDS), anemia, secondary infection, and acute cardiac injury were common complications noted more in patients needing intensive care. described clinical characteristics of 138 cases of NCIP (novel Coronavirus infected pneumonia) admitted in Zhongnan Hospital of Wuhan. They reiterated that fever, fatigue and dry cough were the commonest symptoms. Majority of the patients were in their fifth decade. 26% patients needed intensive care, with ARDS, shock and arrhythmia being the major complications.
The epidemic has been second most severe in Italy so far. Albarello et al. (2020) has discussed two of the cases from Italy which helped understanding some diagnostic aspects of the disease. Holshue et al. (2020) reported the first case of 2019-nCoV in the United States. Li and De Clercq (2020) explored therapeutic options with repurposing of older approved (favipiravir and ribavirin) or experimental antivirals (remdesivir and galidesivir) in view of highly conserved genetic sequence of the cat6alytic sites of the four 2019-nCoV enzymes which have high sequence similarity with those of SARS and MERS Coronaviruses. Baloxavir marboxil (influenza inhibitor of the cap-dependent endonuclease), interferons, ribavirin, Griffithsin (a red-alga-derived viral glycoprotein binding lectin), Lopinavir and Ritonavir (approved anti HIV drug), Nitazoxanide, Chloroquine are few of the drugs under scrutiny for use in COVID-19. Zhu et al. (2020) very recently published a systematic review of all ongoing clinical drug trials relating to this. Among many drugs studied so far, hydroxychloroquine, though pending FDA approval, appears promising and its use is supported by multiple studies like and Cortegiani et al. (2020) . Though development of targeted vaccine in this mutating virus may be effective in halting future epidemics, their use in the current emergent scenario seems to be limited, as it may take months or years to develop an effective, widely tested vaccine.
The epidemic spreading wildly, there have been a surge of statistical studies to predict the number of cases affected by the virus. Anastassopoulou et al. (2020) , Gamero et al. (2020) and are a few notable examples in this regard. Elsewhere, researched on the effect of individual reaction and governmental action (holiday extension, city lockdown, hospitalisation and quarantine) on the spread of the infection. Chinazzi et al. (2020) studied the effect of travel restrictions on the same.
We, in this study, work with an interesting aspect of the COVID-19 outbreak. We develop a novel mathematical model to study the trend pattern of the incidence, defined by the number of new cases identified everyday, for different provinces in China and for six different countries. The implication of our work is three-fold. First, it provides insight on the rate of growth of the number of affected people. In particular, we discovered that a quadratic trend was able to explain the data best. Second, we estimate the time when there is a significant change in the growth pattern of the disease. For China, both on country-level and on province-level, it was found that the growth ceased towards the end of February. In contrast, the other countries are seeing a surge in the growth rate from around the same time. The results from this analysis can help one identify possible factors which directly or indirectly dictate the change in the growth pattern of COVID-19 incidence, for better (enabling similar strategy initiation in different population or early in the same population), or for the worse (which might indicate challenges that are to be prioritised). Finally, we estimate the basic reproduction number for different countries using the fitted values of our model. It was found that the estimates are consistent with the range provided by WHO.
Data source:
All data related to COVID-19 were collected from JHU-CSSE (2020), the GitHub repository provided and maintained continuously by the Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE). We also use population data, both on province-level for China and on country-level for the world. These were collected respectively from Statista (2020) and Worldometer (2020) .
Variables:
The data consist of three main variables -number of confirmed cases, number of deaths and number of people recovered in different regions. Throughout the paper, they will be denoted by y k,t , d k,t and r k,t , respectively. In a similar fashion, population is denoted by P k,t . Here, k is used for a region and t is a time-point. The data is provided on daily level, starting from 22nd January, 2020 to 21st March, 2020.
In our analysis, we also study the effect of lockdowns on the growth rate of the number of affected cases. The Chinese province of Hubei went into complete lockdown on 23rd January, 2020. A few areas in other provinces also went to partial lockdown, and observed something called "closed management". This means that those areas would only keep one entrance and exit point open, and each household is allowed limited numbers of entrances and exits per day. In addition, night-time access was prohibited in some areas as well. Our method is able to quantify the effect of lockdowns (full or partial) in respective provinces, once we account for the general time trend. For that, we create a binary variable L k,t to denote pre (corresponds to L k,t = 0) and post (corresponds to L k,t = 1) lockdown periods.
Any other variable used in the analysis will be defined later, as needed.
Throughout this paper, we are primarily interested in the prevalence and incidence of the disease. Prevalence, in epidemiology, is defined as the proportion of people in a population affected by the disease. In this study, we calculate it as the number of cases per 1,000,000 people. Thus, if y k,t and P k,t denote respectively the number of affected cases and population (in millions) of region k at time t, then prevalence (π k,t ) is computed as π k,t = y k,t P k,t .
(2.1)
On the other hand, incidence is defined as the proportion of new cases per 1,000,000 people. With the same variables as above, incidence (θ k,t ) is computed as
Another contribution of this paper is to provide country-wise estimates for the basic reproduction number, hereafter denoted by R 0 . In epidemiology, R 0 of an infection is the expected number of cases directly caused by one affected case when all individuals in the population are equally susceptible to infection.
In order to estimate R 0 , we use the serial interval (I) as well. I, in epidemiology of infectious diseases, refers to the time between successive cases in a chain of transmission. One can estimate it by the time interval between infection and subsequent transmission. Thus, if the timing of infection transmission in deference to a person's clinical onset is T 1 , and the incubation period of a subsequent case is T 2 , then
3 Methods
For most parts of the statistical analysis in this paper, we use auto-regressive integrated moving average (ARIMA) methods with time-dependent parameters. The goal of this study is to identify if the coefficient estimate for the time trend changes after a certain point, and if so, we examine it further to find out the nature of change and potential causes. Our model is defined below using the notations from last section, but we drop the subscript k for simplicity. Formally, let θ t denote the incidence rate (see eq. (2.2)) of COVID-19 at time t ∈ {1, 2, . . . , T }. P t denote the overall population (in millions). L t is the dummy variable to signify lockdown. As mentioned earlier, every time point t denotes a single day. We shall use τ to denote the day when the first confirmed case is observed, i.e. y t = 0 for t < τ and y t > 0 for t ≥ τ . In the data, it is observed that the number of affected cases (denoted by y t ) grows with time, more than linearly at first and then stabilizing after a certain time. Keeping that in mind, we consider the following structure for the logarithm of θ t .
(3.1)
Here, f (t, τ ) is the trend function, γ captures the effect of lockdown, and u t is the error process. We want to point out that one can include more additive terms in the mean structure to analyze effects of other covariates (for example, policy change, healthcare system etc), if relevant information is available. In eq. (3.1), we assume an ARMA structure for the error process while f is assumed to be polynomial in nature. The maximum degree for the polynomial function was considered to be 3 at first. But, experimentation with the data suggested that a quadratic trend in (t − τ ) is best suited for the incidence model, across different regions. This quadratic function is defined in a time-dependent way. In particular, the coefficients for the linear and quadratic terms in f are considered to be different for (t − τ ) < η and (t − τ ) ≥ η. The value of η is crucial in this study. It is estimated from the data and it tells us when the trend of the growth changes its pattern. The corresponding coefficient estimates then provide an idea about how the growth is changing, for better or for worse. Combining everything, following is the equation of the model we use:
Here, p and q denote the auto-regressive (AR) and moving-average (MA) orders of the error process u t . t denotes a standard normal white noise process and for m = 1, 2,
Throughout the paper, we estimate η, p, q from the data using Akaike information criteria (AIC).
After fitting the model, we run a few residual diagnostics. First, a Box-Ljung test is performed on the residuals to test if any of the autocorrelation is significantly different from zero. These tests were done at 5% level of significance. Fitted vs residual plot, QQ plot of the residuals, partial autocorrelation function (PACF) plot of the residuals are also explored to confirm that the model is appropriate for the data. Next, predictive ability of the model is evaluated. For every site, we use all but last three days' data to train the model, and use the estimated coefficients to predict the log-incidence of those three days. The prediction, hereafter denoted as log θ t , is then transformed to incidence and subsequently the number of new cases for that day is computed. Root mean squared error (RMSE), as defined by eq. (3.4), is finally calculated to understand how well the model can forecast the future.
For every country, we use the number of new cases per day from our model and use that to estimate R 0 using maximum likelihood method. Note that we need to use the distribution of the serial interval (I) in this process. Unfortunately, that is still not known for 2019-nCoV. So, similar to Zhao et al. (2020), we use different distributions of I to estimate the reproduction number. There are three candidate distributions we have used in this study: (a) same as SARS (gamma distribution, with mean 8.4 and standard deviation 3.8), (b) same as MERS (gamma distribution, with mean 7.6 and standard deviation 3.4), and (c) average of SARS and MERS (gamma distribution, with mean 8.0 and standard deviation 3.6).
All of the above methods are implemented using RStudio version 1.2.5033, along with R version 3.6.2. All ARIMA related computations are done using the forecast package. The estimation of R 0 was performed with the help of R0 package.
We apply the above methods in two parts. First, we focus on the province-level daily data from China and observe how the trend of the incidence changes in different parts of the country. These results are described in Section 4.1. Next, we consider the country-level daily data from China, India, Iran, Italy, South Korea and the United States of America; and discuss the findings in Section 4.2.
There are 31 Chinese provinces affected by COVID-19. Maximumn number of affected cases, as expected, is in Hubei. Till date, there have been 67800 reported cases in Hubei, among which 3139 (4.63%) died. Next highest number of reported cases is in Guangdong (1400), although the number of deaths is only 8 (0.57%) there. Tibet and Qinghai have seen the lowest number of confirmed case, with 1 and 18 respectively. Their population are also the lowest of the country. On that note, prevalence, defined by eq. (2.1), is a better measure to understand the extent of outbreak. Hubei has observed a prevalence of staggering 1145.85 cases per million people. Next highest on the list is Beijing (23.40), followed by Zhejiang (21.54) and Jianxi (20.12). So far as the fatality is concerned, other than Hubei, Xinjiang (3.95%) and Hainan (3.57%) are the other two provinces to have a more than 3% rate. Recovery rate is more than 75% for all provinces.
Overall, there have been 81305 confirmed cases in China, among which 4.01% have died and 88.38% have recovered. A heatmap of the whole country based on the prevalence (changed to logarithmic scale) is shown in Figure 1 . Hubei is a clear outlier. It is also evident from the map that the provinces closer to Hubei were usually affected more than the others.
Model (3.2) is now applied to the log-incidence data of every province separately. Note that the length of every time series is 60, as we are dealing with daily data from 22nd January, 2020 to 21st March, 2020. As mentioned in Section 3, AIC is used to identify the date when the trend pattern of the log-incidence changes. Most suitable AR and MA orders were also chosen accordingly. In Table 1 , we present the most appropriate ARMA order, effect of lockdowns (partial or full) wherever appropriate, p-value from the Box-Ljung test on residuals and RMSE (on the number of new-cases) for three-day-prediction for the best models corresponding to different Chinese provinces. The results show that differencing was necessary for all provinces, thereby suggesting that all of the time series were nonstationary. For Qinghai, which has observed one of the lowest prevalence rates, the most appropriate ARMA order was found to be (6, 1, 0). This describes a significant dependence with up to lag 6, that is the previous 6 days. All other provinces showed dependence with at most prior 3 days. The last two columns of Table 1 report diagnostics results. The p-value of the Box-Ljung test on the residuals was more than 0.05 for all provinces. It establishes that the residuals are not autocorrelated, thereby showing the adequacy of the model to capture the time dependence pattern appropriately. On the other hand, RMSE of the predictions of new cases in a three-day-ahead time window was at most 7.35. Hence, the predictive ability of the model is good, so far as the Chinese provinces are concerned.
In the above table, we also present the effect of lockdown (partial or full) in provinces who imposed such rules. It was found to be significant only in a few provinces, and some of them still showed positive effect. This might be attributed to the theory that lockdown only ensured no exposure, but a large proportion of population were already affected and got tested thereafter.
The date on which the trend pattern of the log-incidence series changes significantly is one of the most interesting contributions of this study. The estimates and the standard errors of the linear and quadratic trend before and after the respective dates for all provinces are displayed in Table 2 .
The corresponding dates vary between 24th January, 2020 to 5th March, 2020. But, most of them are concentrated in the week of 17th to 24th February. The obvious next question is whether the trend changes to halt the growth or to increase it. That aspect was explored using the estimated coefficients and fitted values for the trend (0) Figure 2 , we show the log-prevalence of COVID-19 for different provinces. On each line in this graph, we have marked the date when the trend function for the log-incidence rate changes significantly. The direction of the change is also denoted in this graph, using × for decreasing and • for increasing. It is clear that all provinces but one (Anhui) observed a decreasing growth pattern after the respective date.
For this paper, we work with the COVID-19 cases from six countries -China, India, Iran, Italy, South Korea and the United States of America (US). These nations are considered for different reasons. China, being the epicenter of the outbreak, is definitely the most crucial country of the lot. While we have already performed a province-level analysis for China, we want to explore the efficiency of the modeling scheme on the country-level data as well. This would help us in confirming the conclusions from the previous section, thus establishing robustness of the model. Italy, Iran and South Korea were hit early with the outbreak. While they have continued to be three of the top ten countries with highest number of confirmed cases, they showed different patterns when it comes to fatality rates. Italy and Iran have 9.01% and 7.55% fatality rates respectively, both being much higher than the overall average. South Korea, on the contrary, has observed only 1.16% fatality rate. Moreover, Korea have already adopted several screening and healthcare policies, and that may have helped them with the situation. It is also worth mention that Iran has a recovery rate which is on a higher side (37.05%). This is second among all countries in consideration in our paper, bettered by only China who are reporting a 88.38% recovery rate as of 21st March, 2020.
Last but not the least, US and India got affected by the virus a bit later than the aforementioned countries and the growth rate is still believed to be at an early stage. Alarmingly, US is already at the third spot in the global list. India, meanwhile, is the only country in our study from outside the top ten. Despite being a nation of more than a billion people, it has reported only 330 confirmed cases, 5 having passed away till date. Refer to Table 3 for a detailed comparison of these numbers. Our main aim is to explore whether our model can provide valuable insight about the trend of the log-incidence of all of these countries. In addition, we also analyze the effect of lockdowns in Italy and China.
Looking at Table 4 , we can say that the complete lockdown helped Italy reduce their number of affected cases significantly. The lockdown in China, which was considered to be 2nd February, 2020 (as major provinces implemented their restrictions by then), was not significant in the model. This is in line with what we saw in Section 4.1. We also arrive at a similar conclusion about the nonstationarity and the most appropriate ARMA structure of the time series. South Korea showed a significant autocorrelation with prior three days. But, no other country has displayed the same beyond lag 1. It was further confirmed by the diagnostics we performed. The Box-Ljung test on the residuals turned out to be nonsignificant (p-values were greater than 0.05) for all countries. The prediction RMSE was good for most countries, except for Italy and US. Next, we look at the dates when the trend function changed its curve significantly for these countries. The estimates and corresponding standard errors for the coefficients of the trend function of the models are presented in Table 5 .
For China, it was 19th February after which the log-incidence started to decrease (refer to Figure 3 ). This is in accordance with what we observed in the previous section.
Interestingly, contrary to what we noticed for China, the other countries showed a different pattern. For all five of these countries, the growth has not yet started to cease. Rather, it has been increasing more. Situations have been getting worse in Iran and Italy since the last week of February. US have seen the rise in confirmed cases after 10th March and it continues to be disastrous. India, albeit have managed to contain the pandemic so far, are seeing an increasing trend since 16th March and the number may rise drastically in coming days. For South Korea, on the other hand, the change in trend pattern came on 20th February. While the graph depicts that the situation has more or less stabilized over the last few days, statistically we cannot confirm it yet. (0) 
R 0 was estimated for all of the above countries, using different serial interval distributions, as mentioned in Section 3.3. These results are reported in Table 6 . Earlier, Chen et al. (2020) estimated R 0 to be 2.30 from reservoir to person and 3.58 from person to person. reported that the reproduction number varies between 2.24 to 3.58. Our estimates for all countries but the US fall in this range. For India, it is the lowest. But, that is also in accordance with World Health Organization (2020) who reported that R 0 for the human-to-human transmission ranges from 1.4 to 2.5. R 0 for US is surprisingly high and might be a reason to worry. As Ridenhour et al. (2018) pointed out, R 0 may vary geographically "owing to changes in environment, population structure, viral evolution, and immunity". Different healthcare and immigration policies might also be responsible in this regard. 
It is important to note that the time passed since the outbreak has been more than 2 months and we are yet to find a vaccine or an approved drug for treatment or pre-exposure prophylaxis for COVID-19. Thus, at the moment, epidemiological rather than pharmaceutical solutions need to be implemented, meanwhile promoting active research in finding a curative medicine. These measures need precise depiction of the ongoing disease and realistic growth predictions for proactive measures and continuous re-assessment of implemented regulations. This study endeavoured to sketch this very scenario using a statistical model taking in account the most recent data. Using appropriate diagnostics, we confirmed that the model is able to capture the time-dependence pattern and the time trend of the disease-incidence. Combining all of the above analysis, it can be inferred that by the beginning of March, all provinces in China managed to contain the virus and consequently, the growth in the number of confirmed cases has been minimal afterwards. On the global scenario, all other countries are in a bad state. The growth rates of the incidence for these countries are increasing at the moment. As we analyzed the effect of lockdown, we found that for China, both on country-level and on province-level, it was not significant. Italy, however, was benefitted from going into complete lockdown.
We further want to mention that the spread of the COVID-19 across the globe has been described to occur in four stages. Stage 1: Cases are imported from affected countries. At this stage there is no spread of the disease locally. Stage 2: There is local transmission from infected persons, relatives or acquaintances of those who travelled abroad. During this time, source is known and it is easier to do contact tracing and self-quarantine. Countries like India are in stage 2. Stage 3: Community transmission has already begun. People with no prior history of travel or contact with affected persons test positive, meaning the source cannot usually be traced. We think US is one of the countries currently in Stage 3, as the number of cases has been growing substantially since 10th March, 2020. Stage 4: Phase of the infection where it reaches epidemic proportions, such that containing the spread is extremely difficult. According to our analysis, Italy, where the growth has an increasing trend since 22nd February and is yet to show a sign of decrease, is in stage 4.
Clearly, the modeling approach gives an insight on what epidemiological stage a region is in. This has the potential to help in prompting implementation or change in approach to address COVID-19 in different countries.
The main focus of this study has been on the trend pattern of the incidence of COVID-19. While we incorporated information on lockdown as best as we could, more detailed information on screening and testing measures, policies passed on different countries, effective preexposure drug profile access and robustness of healthcare services can be more effective in analyzing the data. Should that be available, one can easily include that information in the mean function of the model (cf. eq. (3.2)), given the flexible structure we consider.
On the other hand, in order to estimate R 0 in the most appropriate way, a detailed account of case-level data and case-to-case transmission data is needed. The same can be said about the estimation of serial interval as well. A future direction of this paper will be to use case-level data, once it is available, to estimate these two variables, which are crucial in an epidemiological study of an infectious disease. Figure 2: log-prevalence of COVID-19, along with the direction of change in the trend function for log-incidence rate (plotted against the date when that change occurs), for different provinces in China. Anhui is the only province to show increasing behavior.