The coronavirus disease 2019 (COVID-19), i.e., a disease outbreak of atypical pneumonia that originated from Wuhan, China [1] , has caused globally at least 6 million confirmed cases, including an estimated 370,000 deaths in approximately 200 countries by the end of May 2020 [2] . The World Health Organization declared the COVID-19 crisis a pandemic on 11 March 2020 [3] .
In modelling the disease's transmission as well as to inform and evaluate control policies, it is particularly important to estimate its reproduction number. Early estimates for COVID-19 basic reproduction number R 0 , that denotes the transmission potential of infectious disease when introduced to a completely susceptible population, ranged from 1.4 to 6.49 [4] . The effective (or instantaneous) reproduction number R t , on the other hand, reflects the extent of transmission in the presence of population immunity or intervention. Thus, the estimation of R t is important for evaluating public measure success.
However, estimation of R t is sensitive to the model structure and parameter assumptions [5] . As a case in point, due to incorporation of more individual case information and travel data, the estimate for R 0 in Wuhan was revised upward from 2.2-2.7 to 5.7 [6] . On the other hand, data inavailability or poor quality often hinders the use of certain estimation methods, such as serial interval data that are usually needed to estimate R t (e.g., Fraser [7] , Wallinga and Teunis [8] , Cauchemez et al. [9] , White and Pagano [10] ).
In the course of calculating the exact value of R t , especially when the data has not yet reached its peak, precise assumptions and data estimates are needed. Nishiura et al. [11] discussed a likelihood-based approach to estimate R t from early epidemic growth data. Using the compartmental Susceptible-Infectious-Recovered (SIR) model, Bettencourt and Ribeiro [12] use the incidence data to estimate R 0 and R t . In this paper, based on the Susceptible-Infectious-Recovered-Dead (SIRD) model as a reference, we develop a novel approach to estimate R t of COVID-19. It uses information on the number of infected or active (I), recovered (R), and death (D) cases, which are readily available for all affected countries [13] , so that they can be accessed rather easily.
The reproduction number is estimated from reported cases under uncertainties using a two-stage estimation method based on the Extended Kalman Filter (EKF) and a low-pass filter. The method not only considers the nominal number of reported cases, but also its daily pattern. To show our method's practical ability, we apply it to COVID-19 cases in the Scandinavian countries, i.e., Denmark, Sweden, and Norway, and compare the results with two commonly used Bayesian methods due to Bettencourt and Ribeiro [12] and Cori et al. [7, 14] . We show that the results are indeed comparable.
Our paper is structured as follows. In Section 2, we discuss the mathematical model that will base the method. We then discuss the two-stage estimation method in Section 3. In Section 4 we apply the method to estimate the effective reproduction number of COVID-19 in Denmark, Sweden and Norway. We also compare the results with R 0 and R t calculated using the methods of [12] and [14] , respectively. We conclude our work in Section 5.
Our estimation method is based on the compartmental SIRD model that can be written as the following first-order nonlinear differential equationṡ
where S denotes the number of susceptible individuals, N is the total population, β is the average number of contacts per person per time, γ and κ are the recovery and death rate, respectively. Remark that the value of β is time-varying due to intervention, i.e., β=β(t). To use the model, we require information on the average infectious time T i and the case fatality rate CFR, so that
For COVID-19, we take T i = 9 as the infectious period on average lasts for 9 days (i.e., 7-11 days with 95% CI) [16] . CFR are taken based on the information from [13] . The time-varying effective reproduction number is then given by
The approximation is under the assumption that government intervention is taken at an early stage so that the susceptible is relatively the same over time as the total population. This is the case especially for emerging diseases. We modify the SIRD model by augmenting the following two equations into the system:Ḋ
The former equation takes into account the daily number of reported new cases DR [17] , while the latter one says that the effective reproduction number R t is assumed to be a piece-wise constant function with jump every one day time interval.
Discretizing the model using the forward Euler method, we obtain the fol-lowing discrete-time augmented SIRD model
Our method computes a new estimate of R t based on new reported cases. Since their frequency is low (could be once a day), the reported data can be interpolated using, e.g., a modified Akima cubic Hermite interpolation, such that it fits with the time step ∆t. In our simulation, the time step ∆t is chosen as 0.01, i.e., 100 time discretization within one day interval. The confidence interval of our estimated R t is determined by computing the reproduction number for different values of the infectious period T i within a certain interval.
To simplify the presentation, we define the augmented state vector
and as such, the discrete-time augmented SIRD model (8)- (13) can be written as follows
where f is the nonlinear term written in the right hand side of (8)-(13) and w is introduced as an uncertainty to model the inaccuracies due to simplification in the modelling. The uncertainty is assumed to be a zero mean Gaussian white noise with known covariance Q F . In practice, Q F can be considered as a tuning parameter for the EKF. Thus, the transmission model becomes a discrete-time stochastic augmented SIRD model. Reported cases, such as the number of active cases and the cumulative numbers of recovered and death, can be incorporated into the model using the following output vector
Here, v denotes uncertainties due to false testing results. We also assume the uncertainty to be a zero mean Gaussian white noise with known covariance R F . As well as Q F , R F can also be considered as a tuning parameter. Following the available data that include I, R, D, and DR, the data/measurement matrix C is taken to be
A two-stage filtering method is used to estimate the daily reproduction number R t . The method consists of the EKF and a low-pass filter. In the first stage of estimation, the EKF is used to estimate the state variables and the value of R t under uncertainties in the number of reported cases. Afterwards, the low pass filter is used to remove short term fluctuations of the reported cases that can be caused by delays in the reporting. For example, suddenly in Denmark there were 893 recovered patients reported on 1 April 2020, in contrast to the previous days from 16 February 2020 onwards when there was no recovery reported at all. Such an accumulated delay can cause a falsely decreasing value of R t .
The EKF is an extension of Kalman filter for nonlinear systems. The Kalman filter itself is based on a recursive Bayesian estimation and is an optimal linear filter. The idea of EKF is to linearize the non-linearity around its estimate. Due to that linearization, the optimality and stability of the EKF cannot be guaranteed. However, if the non-linearity is not severe, the EKF can give a reasonably good estimate.
Let us denotex(k) as an estimated vector state from the EKF. Applying first-order Taylor series expansion to f atx(k), we obtain
where J f (x(k)) is the Jacobian matrix of f , given by
where
The EKF consists of two steps: predict and update. The discrete-time stochastic augmented SIRD model is used to predict the next state and covariance and update them after obtaining new data/measurement. The EKF can be considered as one of the simplest dynamic Bayesian networks. While the EKF calculates estimates of the true values of states recursively over time using incoming measurements and a mathematical process model, recursive Bayesian estimation calculates estimates of an unknown probability density function recursively over time using incoming measurements and a mathematical process model [18] . Letx(n|m) denotes the estimate of x at time n given observations up to and including at time m ≤ n. The Kalman filter algorithm is given as follows [19] Predictx
Updateỹ (k + 1) = y(k + 1) − Cx(k + 1|k) (28)
x(k + 1|k + 1) =x(k + 1|k) + K(k + 1)ỹ(k + 1) (30)
Here P (k|k) denotes a posteriori estimate covariance matrix. In the second stage, a low pass filter based on a rational transfer function is used to remove short term fluctuation at time step k, and is given bŷ
where y n is a window length along the data. In our case, we choose y n = 3 ∆t .
To evaluate the quality of the estimate, we calculate Relative Root Mean Square Error (RRMSE) between the estimated and reported cases. The RRMSE is defined as
where N d is the number of observed days and X ∈ {I, D, R, DR}.
In this section, we apply our method to study viral transmission of COVID-19 in Denmark, Sweden, and Norway. All datasets and MATLAB code are available on GitHub (https://github.com/agusisma/covid19). As the country of residence of the corresponding author, Denmark has a high estimated percentage of reported symptomatic COVID-19 cases, i.e., 61% (42%-87% of 95% CI) [20] . Norway is at 73% with 43%-98% of 95% CI [20] . On the contrary, Sweden only has 15% (11%-19% of 95% CI) [20] . The countries also have a different approach in their public measures in responding to COVID-19, i.e., Sweden did not implement a strict lockdown, unlike its Nordic neighbouring countries.
We plot the observed incidence of COVID-19 in Denmark, Sweden, and Norway in Fig. 1 . We also plot in the same figure estimated numbers computed using our method, where good agreement is obtained. For all estimation, the process and observation covariance matrices are considered as tuning parameters and are chosen as Q F = diag(10 10 10 10 5 0.2) and R F = diag(100 10 10 5 1), respectively. These parameters are chosen such that the RRMSE between the estimated and reported data are sufficiently small. In our case study, the RRMSE are shown in Table 1 In applying our method, we also compare it with two commonly used methods to estimate transmission parameters, namely the sequential Bayesian method of Bettencourt and Ribeiro [12] that provides an approximation of the basic reproduction number, and the instantaneous method by Fraser [7] that is implemented with a Bayesian analysis [15, 14] . The former method exploits the new reported incidence, while the latter one uses the distribution of the serial interval.
First, we compare our method with that of Bettencourt and Ribeiro [12] , that allows sequential estimation of the basic reproduction number at the initial stage, i.e., when the growth is still exponential. While the two methods are based on the SIR model, Bettencourt and Ribeiro [12] rather use new incidence data, which in here need to be smoothed out using a five-day moving average filter. In Fig. 2 , we plot the comparison and summarise the basic reproduction numbers that are taken to be the maximum of the curves in Table 2 . It is interesting to note how the methods give rather similar estimations. This indicates that our method gives comparable results to those of [12] . Finally, we plot the time-varying effective reproduction number in Fig. 3 . Here, we compare our results with those using Cori et al. [14] . The method of [14] utilises the disease serial interval, which we approximate using a shifted Gamma distribution [14] with mean 4.7 and standard deviation 2.9 [21] . The prior belief for the value of R t is taken to be Gamma function with mean and standard deviation 5. We do not average out the data of daily new cases, but instead take the likelihood estimation of a new case at one day to depend also on the estimation of the previous three days. In Fig. 3 we obtain that the two methods give the plot of R t with the same trend, indicating that our method is also comparable with [14] . There is a delay of about four days in the trend, especially with the time when the reproduction number curve crossed the horizontal axis. The delay is caused by the peaks of new daily cases and active ones that also differ by about the same days.
A different trend especially at later times between the methods appears in Fig. 3(c) for Norway. The curve from our method is quite smooth, while it is rather fluctuating in that using Cori et al. [14] . The discrepancy is caused by the active and recovered cases that apparently were not updated regularly, in contrast to the new positive cases needed by the method of [14] . The unreported recovery cases were all released at once on May 22, 2020, see Fig. 1 (c).
Many mathematical models have been developed to estimate several types of reproduction numbers during epidemic outbreaks. We here provide a novel method exploiting reported active, recovered and death cases using the SIR model as an underlying approach. This new method offers several advantages, from the modeling point of view, the resulting R t value can follow the dynamics of the model suggested, so it is possible to develop it further if the model chosen has a higher complexity, besides that the estimation approach used can still be expanded in terms of statistical view. In the case that the data provided in time series do not change much or instead have drastic changes, such as accumulating at a certain time, the resulting R t value will show the same spikes and serrations. As a result, the latest information from data dynamics can be more elaborated. By applying the method to COVID-19 cases in the Scandinavian countries and comparing the results to commonly used methods due to [12] and [14] , we showed that our model is comparable, which expectedly will allow for fast assessment of the reproduction number in new outbreaks.
Using the method to forecast and critically assess incidence data in countries with high under-reporting, such as Indonesia, is addressed for future work.