Estimation of the “effective” SIRD model parameters
Here we note that the new cases of recovered and deaths at each time time t appear with a time delay with respect to the actual number of infected cases. This time delay is generally unknown but an estimate can be given by clinical studies. However, one could also attempt to provide a coarse estimation of these parameters based only on the reported data by considering the first period of the outbreak and in particular the period from the 11th of January to the 16th of January where the number of infected cases appear to be constant. Thus, based on Eqs (3) and (4), and the above assumption, the “effective” per day recovery rate β and the “effective” per day mortality rate γ were computed by solving the least squares problems (see Eqs (2) and (4): γ=[(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))T(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))]−1(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))TΔD(t),(14)  and β=[(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))T(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))]−1(CΔI(t−1)−CΔD(t−1)−CΔR(t−1))TΔR(t),(15) 
As noted, these values do not correspond to the actual per day mortality and recovery rates as these would demand the exact knowledge of the corresponding time delays. Having provided an estimation of the above “effective” approximate values of the parameters β and γ, an approximation of the “effective” infected rate α, that is not biased by the assumption of S = N, can be obtained by using the SIRD simulator. In particular, in the SIRD model, the values of the β and γ parameters were set equal to the ones found using the reported data solving the corresponding least squares problems given by Eqs (14) and (15). As initial conditions we have set one infected person on the 16th of November and ran the simulator until the last date for which there are available data (here up to the 10th of February). Then, the optimal value of the infection rate α that fits the reported data was found by “wrapping” around the SIRD simulator an optimization algorithm (such as a nonlinear least-squares solver) to solve the problem: argminα{∑t=1M(w1ft(α;β,γ)2+w2gt(α;β,γ)2+w3ht(α;β,γ)2)},(16)  where ft(α;β,γ)=CΔISIRD(t)−CΔI(t),gt(α;β,γ)=CΔRSIRD(t)−CΔR(t),ht(α;β,γ)=CΔDSIRD(t)−CΔD(t) where, CΔXSIRD(t), (X = I, R, D) are the cumulative cases resulting from the SIRD simulator at time t; w1, w2, w3 correspond to scalars serving in the general case as weights to the relevant functions. For the solution of the above optimization problem we used the function “lsqnonlin” of matlab [20] using the Levenberg-Marquard algorithm.