After a period of confinement due to the presence of COVID-19 and facing economic and social pressures, societies start to open up, seeking to return to productive, sport, and recreational activities. As the interactions between people increase, the impending menace of outbreaks materializes. Naturally, there is a tendency to apply once again lockdowns, in what has been called the hammer and the dance [22] . In this context, it is essential to have easy to apply indicators for people to use as a reference. The basic reproduction number, R t , fulfill such a role [15] . When R t is higher than one, the number of infected people grows exponentially, i.e., their number will double in a short period. When R t is less than one, the epidemic will tend to disappear. However, estimating R t accurately at the required level of geospatial resolution is a complex problem.
Although applicable to any country, let us take the case of Mexico, as an example. The records generated by the epidemiological surveillance system contain information that includes, among other predictors, the number of confirmed positives, deaths, and suspects. Daily, the Ministry of Health informs the public about the status of its records [26] . However, the data it discloses updates records of events that occurred in the past, sometimes as far as 50 or 60 days ago. At other times, with a significant frequency, the records that were previously released are discarded. Although publishers often drop these erroneous entries overnight, there have been cases of records eliminated after more than 50 days.
Besides the integrity of the information, there are other difficulties in tracking the epidemy inherent to the pandemic and interesting for researchers, decision-makers, and the general public. SARS-CoV-2 is an airborne virus [4] , which infects some people without causing symptoms [19] . On a significant number of occasions, people begin to spread COVID-19 before they start to feel sick [13] . Also, each infected person reacts differently and will have, if anything, a different latency and incubation period [12] . People will have a different contagious period, manifested with inequal intensity during that time [7] . Although the symptoms are known, one may reveal them differently. People will require different types of medical attention, which may or may not require hospitalization [11] . In some cases, someone ill may need or not a ventilator [18] . Eventually, a given person may recover, possibly with sequels, or will pass away [24] . About the whole process, we begin to have some statistical knowledge on which we can develop models. In this paper, propose a data-driven approach that leverage experience to create a simple, yet effective nowcasting method for R t that can be used by policy-makers as well by the general public. Our main contribution is an approach to use past observations to generate plausible sequences of estimates for the number of confirmed positive cases that could have possibly occurred in the recent days to compute variations of the basic reproduction number.
We base our method on the statistical behavior of previous observations. As more information arrives, the estimation naturally becomes more precise about the final count of confirmed positives. In the next section, we review the literature about related methods. Then, in §3, we discuss the intrinsic delay in information flow that exists in the process of detecting a COVID-19 confirmed positive and detail our approach to estimate plausible sequences for the number of infected people. In §4, we review the underlying method we employ to estimate the basic reproduction number out of the possible sequences. We show our implementation of the nowcasting method for R t in §5. Finally, we conclude our study summarizing our findings and delineating some potentials lines of research.
Though recent, COVID-19 has kickstarted some novel ideas to track it reliably. The research effort to nowcast the basic reproduction number can be classified in either mechanistic approaches, Bayesian approaches, or a hybrid combination of both.
Wang et al. [27] developed a hybrid model to complement the dynamics of the SIR (Suspected, Infectious, Recovered) model with spatiotemporal analysis. The space-time component is modeled, at the start, with a Poisson distribution to describe rare events. Then, they complemented it with a negative binomial random model during over-dispersion. Balabdaoui and Mohr [5] propose an age-stratified discrete compartment model as an alternative to SIR type models. Their approach follows the trajectory of individuals that includes the exposed, the asymptomatic, the symptomatic infectious, the symptomatic in self-isolation, the patients in the intermediate care unit, and the patients in the intensive care unit. Masjedi et al. [16] compares phenomenologic and mechanistic models. The former based on generalized Richards models [23] (an extension of sigmoid functions) and the latter on a modified SEIR (Susceptible, Exposed, Infectious, and Recovered) model. They fit the models with observed data to forecast the next month. They observe that although phenomenologic models fit the data, they are not reliable for decision-making. In contrast, SEIR models predicted the phenomena better. Contaldi [9] presents SIRFH, an extension of the SIR model that tracks hospitalizations and hospital-based fatalities introducing additional differential equations. The estimation for the basic reproduction number derives from the solution to this extended model. Finally, Annan and Hargreaves [3] produce a nowcasting method based on the SEIR model. To calibrate the parameters, they use observational data and a Bayesian approach. Annan and Hargreaves' analysis includes the uncertainties associated with deaths' stochastic nature, the reporting errors, and the model itself. 
Altmejd et al. [2] present a model based on the removal method [21] , where one extracts batches of a fixed population. Their models deal with lags arising from the calendar patterns, where events reported during the weekends are less. Their Bayesian approach uses a likelihood that considers the number of reports by day of the week, and priors with improper uniform distribution. Their model provides better estimates than seven days averages. Schneble et al. [25] present a nowcasting model based on the number of deaths, as quantifying their correct number is more reliable than for infected people. Their epidemic spread model considers region and age-specific Poisson distributions, where they consider lag to report. They model the effect of age, gender, weekday, and location as a quasi Poisson distribution. Then, they infer a posterior using a Gaussian prior. For nowcasting, they model the delay as a random variable which will provide death counts. They distribute these death counts as a quasi-binomial distribution. Chitwood et al. [8] propose to use a Bayesian framework for nowcasting. They take into account delayed and incomplete reporting. They assume that one can understand the COVID-19 complex spread system by examining the individual components. In that model, they consider the uncertainty that results from available diagnosis and delays in the estimation of disease progression and reporting systems. Lastly, Abbot et al. [1] employ a quasipoisson regression model to estimate the spread rate. Interestingly, they base their analysis on the reported dates for the confirmed positives and infer the symptom onset through statistical modeling.
In our approach, we characterize the frequency at which the counting updates of COVID-19 confirmed positives occur. In this section, we analyze the origin of such delays and describe the form we model them.
Declaring a person confirmed positive involves a complex process that may take days, even nowadays, when it is of paramount importance to achieve certainty for decision-making. Just consider the case of a person showing symptoms related to COVID-19 [17] that decides to visit the physician. After an interview to collect some necessary clinic information, the physician chooses to take either a sample from the nasopharynx using a long swab [20] or a CT (Computer Tomography) [14] . In some places, the sample can be analyzed via the RT-PCR(reverse-transcription polymerase chain reaction) [28] in situ with results on the same day but frequently it may take a week or longer to be processed. Afterward, the results will be uploaded in computer systems and summarized for analysis.
In Figure 1 , we illustrate the effect of delays in reporting using the data set made public by the Mexican Health Ministery [26] . The horizontal and vertical axes show the day of onset and the number of confirmed positive cases. Each layer corresponds to the number of cases added to a prior date. Although the number of updates may be significant for a given day, they eventually converge to the total number of confirmed positives for that day, C t (D), for D large, and where t expresses the day of interest (see Figure 2 (a)). If we divide the daily accumulated of confirmed positives C t (δ), for a given day δ, by C t (D), the cumulated distribution will tend to one. We illustrate this in Figure 2 (b)-(c), where we show both, the rate of daily change and the cumulative change. Our approach aims to characterize the variations we observe in these distributions to develop a model for nowcasting.
We aim to estimate the number of confirmed positive cases C t (D) for the day t using the following δ days of reports available. In principle, we would learn about C t (D) when D is cosiderably large. But in practice, D can be as short as one month and a half of daily updates. Given the number of confirmed positives δ days after day t, C t (δ), the number of confirmed positives on day C t (δ + 1) can be expressed as
where ρ t (δ) is the rate of change from one day δ to the next δ + 1, for reference day t. If we solve the recursion, we will have the expression
where one assumes that the daily rate changes over time. In the cases we are studying, the curves expressing the rate of change of the number of confirmed positive relative to the day before, for a different starting day, seem to be somewhat consistent over the samples. In our case, we model ρ t as a random variable, for which we may be able to fit some standard distribution to the experimental samples. Then, on the day t + δ , the best-guess prediction for the number of confirmed positive,
where our newly defined random variable ρ D t (δ) expresses the rate of change from day δ + t to day D. In our approach, we model ρ D t (δ) as a random variable with different parameters for each day δ, for more fine-grained or longer-term prediction. One may find the relationship between ρ D t (δ) and ρ t (δ) by noting that (2 and (3) 
Expanding C t (δ) using the recurrence relationship in (1), we have
from where, after eliminating for the common factors, solving for ρ D t (δ) results in
Given a particular sequence of the observed number of infected people {C 0 (δ), C 1 (δ − 1) . . . , C t (0)}, and the argument of the number of days the report has been updated, we aim to nowcast the basic reproduction number R t , i.e., given the distribution of the rate of change ρ D t (δ), we generate ensembles of sequences aiming to estimate {C 0 (D), C 1 (D) . . . , C t (D)} before proceeding to calculate R t . We first review EpiEstim, a method proposed by Cori et al. [10] , to estimate R t from the observed number of cases.
Cori et al. [10] proposed a Bayesian framework to compute R t , where the number of infected people observed at day t, C t , follows a Poisson process. In a simplification, they assume that the daily observations of infected people are independent. Thus, one may express the likelihood of observing a sequence of infected people between day t − τ − 1 and day t as [10] 
where the transmisibility R t,τ is assumed to be constant over the period [t − τ + 1, t], Λ t = t s=1 C t−s w s is the total infectiousness of infected people at time t, and w = (w 1 , . . . , w t ) T is a mass density probability profile of infectivity profile for an individual. Cori et al. [10] assume that the basic reproduction number R t,τ is a random variable which probability follows a Gamma distribution as [10] 
where a and b are the parameters of shape and scale. Since the Poisson and Gamma probability distributions are conjugate, one can express the posterior in closed form, again as a Gamma distribution, as [10] (9)
from where the mean α and standard deviation β are given by
Given C t (δ), the information about the number of infected people δ days after the day of interest t, and the model for the probability function for ρ D t (δ), we produce N random samples which will correspond to the number of people infected that day. We then compute R t for each of the sequences using the model proposed by Cori et al. [10] . Finally, we calculate the mean and standard deviation for R t to provide the most likely value and uncertainty at one standard deviation. To take into account the difference between the accepted values for the average incubation (five days) [12] and latency periods (three days) [13] , we represent them two days before t.
We took the data set for COVID-19 cases provided by the Mexican Health Ministery corresponding to July 11, 2020. The data set contains 723,668 records, out of which 295,268 correspond to confirmed positives. As time passes by, the number of confirmed positives for a given day t is updated. In Figure 1 , we illustrate how each day the updates stack up a layer of updated registers toward the past. As we accumulate the number of confirmed positive updates, we observe that the total quantity levels off and reaches a maximum at C t (D) (see Figure 2 ). About 98% of reports are filled out by day 33. When we divide the daily updates for the day t by C t (D), we obtain the normalized updates by day and accumulated registers illustrated in Figure 2 .
We then proceed to fit Gamma distributions to the variation of ρ D t (δ). We show illustrations of this fit for δ = 1, 3, 7, 15, 25 and 35 in Figure 3 . Note that δ = 0 is not present as generally the number of reported confirmed positive for C t (0) = 0 causing ρ t (0) to be undefined. Once we have the models for ρ t (δ), we may proceed to generate estimates for the number of confirmed positives for C t (D) using (3). The mean and standard deviation statistics will provide us with the most likely value and an estimate for the uncertainty. We use the same set of randomly generated values to obtain sequences, which we evaluate using the method proposed by Cori et al. [10] to obtain the instantaneous R t . Our implementation considers the pre-symptomatic transmission, i.e., the incubation period, or the time it takes for an infected person to start showing symptoms, is greater than the latent period, or the time from which an infected person can spread to others. Following Bar-On et al. [6] , we assume that the latent period lasts for three days and the incubation period for five days.
We compare the performance of our nowcasting with the proposed by Abbott et al. [1] (see Figure 4 ). In their case, the nowcasting tends to closely follow the number of reported confirmed positives, which gives the undesirable effect of resulting in a descending R t , when it is not. Our proposal, on the other hand, increases its certainty naturally as more information is available. As we increase δ in our nowcasting exercise, there is a tendency to observe fewer cases, because there has not been enough time for the information to arrive. We set a dynamic threshold to stop the nowcasting estimation when for a particular day of analysis, δ, the number of confirmed positive cases is less than 30. Also, we have observed that as the number of confirmed positive is less, the normalized cumulative curves tend to be noisier. In Figure 5 , we illustrate what happens for entities in Mexico where the number of confirmed positive cases is 59,667 (Mexico City), 44,114 (State of Mexico), 15,909 (Tabasco), and 2,667 (Querétaro). We believe that our method works best when the number of positive cases is beyond 2,600 for the observed interval of 90 days. Using this threshold, there are still currently 30 States (out of 32) and 32 Municipios (out of 2450) in Mexico subject to our analysis.
To foster further research, allowing other researchers to verify our results and serve as a stepping stone, we make our code publicly available at https://www.github.com/joaquinsalas/ nowcastingRt.
In this document, we have presented a nowcasting method to estimate the number of confirmed positives. We have shown that this may be the foundation to generate plausible sequences out of which one may determine useful epidemy tracking indicators, such as the basic reproduction number. Our method naturally expresses uncertainty due to the lack of information but eventually gains certainty as more data accumulates.
Our method's strength is that it is based on the self-reported onset of symptoms, in contrast to other methods that use the number of confirmed positives cases accumulated by the report's day to infer this quantity. A potential drawback of our approach is that it relies on a regularity of the update cycle. As researchers implement more sophisticated systems for testing and reporting, the statistics may change. To remedy this potential effect, one may eliminate old observations and We fed the data to the method provided by Abbott et al. [1] , which seems to follow the reported daily cases. Finally, we illustrate the output of our approach, including its area of uncertainty. Thanks to Dagoberto Pulido for implementing Abbott et al. [1] to generate (c).
update the distributions for ρ t regularly. Due to the difference between the incubation and latent periods, and delays in the detection and reporting cycle, our model estimates R t up to several days in the past. We decided to take no further assumptions about the progression of the epidemy. Although potentially some form of state estimation may be possible to implement to fill the gap. We believe that it is crucial to continue developing solutions to quickly, robustly, and reliably estimate indicators such as the basic reproduction number. A possible direction for future research may be to determine the disaggregation level to continue to generate a reliable indicator. The resulting nowcasting methods should compensate for the delays inherent in producing and processing information about this critical, global, and urgent problem. Also, we are planning to study the extend at which our model can be incorporated into dynamics-based models. This enhancement could offer improved nowcasting. 