key: cord-0509217-b2th82q3
authors: Pasquali, Sara; Pievatolo, Antonio; Bodini, Antonella; Ruggeri, Fabrizio
title: A stochastic SIR model for the analysis of the COVID-19 Italian epidemic
date: 2021-02-15
journal: nan
DOI: nan
sha: 81aa52a7013c28383d0911bba134fe43af7bc27f
doc_id: 509217
cord_uid: b2th82q3

We propose a stochastic SIR model, specified as a system of stochastic differential equations, to analyse the data of the Italian COVID-19 epidemic, taking also into account the under-detection of infected and recovered individuals in the population. We find that a correct assessment of the amount of under-detection is important to obtain reliable estimates of the critical model parameters. Moreover, a single SIR model over the whole epidemic period is unable to correctly describe the behaviour of the pandemic. Then, the adaptation of the model in every time-interval between relevant government decrees that implement contagion mitigation measures, provides short-term predictions and a continuously updated assessment of the basic reproduction number.

In this work we apply a stochastic version of the well-known SIR (with the three Susceptible, Infective and Removed compartments) epidemic model to the Italian COVID-19 data. The SIR model has been introduced about ten years after the 1918 influenza pandemic, also known as Spanish flu ( [24] ), and is still popular as a simple tool to approximate disease behaviour ( [34] ). Numerous extensions such as the simpler SEIR and SIRD models or the more complex SIDARTHE ( [19] ) have been subsequently developed to make more reliable assumptions on the epidemic dynamic. To recognize the role of heterogenous contact networks in the transmission dynamics of infectious diseases, an extended SIR model with seven compartments has been developed on the nodes of a network where each node accounts for the mean number of contacts in a typical day, [41] . For the outbreak of COVID-19, alternative models have been introduced ranging from elementary models ( [36] ) to very complex ones including spatial and temporal variations (e.g. [1] , [2] , [38] , [37] ). However, as a parsimonious model able to allowing the measurement and forecast of the impact of non-pharmacological interventions like social distancing, the SIR model still preserves a primary role in the analysis of the early phase of COVID-19 outbreak ( [3] , [8] ). In [12] for instance, the SIR model is combined with Bayesian parameter inference and the effects of governmental interventions in Germany are modelled as flexible change points in the spreading rate. A SIR model with time-dependent infectivity and recovery rate which are estimated by solving an inverse problem is considered in [27] . A hierarchical epidemiological model in which two observed time series of daily proportions of infected and removed cases are emitted from the underlying infection dynamics governed by a Markov SIR infectious disease process is proposed in [39] .

By introducing a time-varying transmission rate, the Authors cover the effects of different intervention measures in dissimilar periods in Italy. In [10] it is assumed that the susceptible population is a variable that can be adjusted to account for new infected individuals spreading throughout a community. With a similar aim, [7] include among the parameters of a SIRD model the initial number of susceptible individuals as well as the proportionality factor relating the detected number of positives with the actual (and unknown) number of infected individuals. An alternative way to account for possible random errors in reporting is proposed in [20] as well.

As mentioned before, several compartments may be added to the SIR model but, in previous work ([4] ) it was shown that while the parameters of the SIR model can be uniquely determined from the time evolution of the normalized curve of removed individuals, the same does not hold for more complex models. Thus, the SEIR and other models should not be used in the absence of additional information that might be obtained from clinical studies. In the present work, since we assume no clinical information, we will use the SIR model. In particular, we introduce a stochastic SIR model to describe the dynamics of the infection, coupled with an observation equation that relates the actual numbers of infected and removed to the observed ones. There are two sources of randomness: the first in the SIR model, introduced through uncertainty in the infection and recovery rate parameters, such that they depend on time; the second in the observation equation, where a random under-reporting error is assumed. This model is amenable to sequential updating of parameters and forecasting, a useful feature when data come in on a daily basis. The updating method is a particle filter with parameter learning [14] for the fraction of individuals in each compartment.

The article has two objectives: on the one side, the exploration of the suitability of a SIR model, while, on the other, it provides some results on the epidemic in Italy. The first finding is that a single SIR model is unable to capture the dynamics of the entire development of the COVID-19 epidemic, because of the numerous policy adjustments by health authorities and different government interventions that took place during the emergency. Hence we fitted one SIR model to each phase between selected government interventions, obtaining a good fit to the data. Notwithstanding, the predictive capability of this model remains very limited. The second finding is that it is very important to use the correct probability distribution for the observation error: failure to do so may produce parameter estimates that seemingly provide a good fit to the data but do not correctly describe the true underlying dynamics.

The article is organised as follows. In Section 2 we introduce the stochastic SIR model and in Section 3 we present a discretised state-space version and the observation equation, where the state is the fraction of susceptible, infected and removed individuals. A Rao-Blackwellised particle filter (RBPF) algorithm for state filtering, state prediction and parameter learning is illustrated in Appendix A, [14] . In Section 4 we study the filtering and prediction problem on data simulated from our model with the help of graphical displays addressing the following: the accuracy and the precision of both parameter estimation and filtering and the accuracy of the forecast (Sections 4.1 and 4.2); the problem of the simultaneous identifiability of the parameters in both the SIR model and the observation error distribution (Section 4.3). In Section 5 we apply our method to the Italian data of both the first and the second infection wave, obtaining a good fit, as stated above, but a forecast that may be valid only in the short term. In the same section we briefly compare the performance of a deterministic SIR model to our stochastic SIR model, showing that the latter is superior. We also consider the problem of assigning the correct observation error distribution using available information on the Infection Fatality Rate and compare the level of the filtered state to the result of a sample serological survey carried out by Istat (Italian national statistical office) and the Italian Health Ministry between May and July 2020. This comparison indicates that our model, if correctly calibrated, provides a realistic assessment of the state of the epidemics. A section with some final remarks concludes the article.

Consider a population and denote by S t the fraction of susceptible individuals at time t, by I t the fraction of infected individuals and by R t the fraction of recovered individuals (survivors and dead).

We suppose that the population is closed, then for every time t S t + I t + R t = 1. The deterministic SIR model can be written

where β is the disease transmission rate, that is, the fraction of all contacts, between infected and susceptible people, that become infectious per unit of time and per individual in the population, and γ is the removal rate. The reciprocal of γ is the duration of infection. The parameters β and γ allow to approximate the basic reproduction number (or ratio, also called basic reproductive number or ratio) that can be thought of as the expected number of infected people generated by an infected individual in a population where all individuals are susceptible to infection. Despite its conceptual simplicity, R 0 is usually estimated with complex mathematical models developed using various sets of assumptions ([13] ). In the above SIR model it holds

where parameters β and γ are unknown and have to be estimated. We suppose that they are subject to uncertainty and change in time as follows:

with w (1) t and w (2) t independent Wiener noises, that is, β t is supposed normally distributed with mean β 0 and variance σ 2 t and γ t is normally distributed with mean γ 0 and variance η 2 t. For alternative ways to introduce stochasticity, see [17] . The parameters σ and η measuring the noise intensity are assumed known and sufficiently small to obtain positive β t and γ t with probability approximately equal to one.

Substituting the expression (2) for β and γ in system (1), we obtain the following stochastic SIR model:

Unfortunately, the introduction of the noise in the parameters β and γ no longer grants the condition S t + I t + R t = 1. We can enforce it by substituting S t by 1 − I t − R t in the second and third equations and removing the first equation to obtain the reduced system

Denoting

t T the vector of independent Wiener processes and by θ 0 = (β 0 , γ 0 ) T the parameter vector, we can rewrite system (4) in vectorial form:

where

We called X t the state of the system, which for COVID-19 is unobservable, and introduce Y t to denote what can be actually observed, in accordance with the terminology derived from state-space modelling.

The vector Y t is characterised in the following section.

3 Parameter estimation, filtering, forecasting, and goodness-of-fit

To estimate parameter θ 0 we propose a Rao-Blackwellized particle filter (RBPF) algorithm based on the Euler discretization of the stochastic system (5):

where we also use t for discrete time to save notation. The method is described in Appendix A. This algorithm allows to jointly calculate, at each time step, the estimated parameter and the state of the system using a noisy observation of the state as input. It is widely recognised that in this pandemic collected data on infected and removed suffer of under-diagnosis and under-detection, ( [40] ) that is, the observations are subject to an observation error. We suppose that each component of the observation vector Y t+1 is given by the product of the corresponding component of X t+1 and a random variable:

where U t,1 and U t,2 are independent beta distributed random variables with shape parameters a and b. (In the following, by U , Y and X with no subscript we mean scalar random variables distributed as U t,i , Y t,i , and X t,i , respectively, i = 1, 2). The observation error in SIR models has been considered by other authors using different formulations (see, for example, [33] ). Finally, we assume that the initial distribution of θ 0 is Gaussian with mean µ 0 and covariance matrix Σ 0 .

The model can also be used to predict the future behaviour of the epidemic. Let y 1:t be the time series of observations up to time t; for a fixed initial state x 0 , the RBPF algorithm provides a sample x (i) 0:t , i = 1, . . . , M , to approximate the posterior distribution of the state p(x 0:t |y 1:t ). Furthermore, the conditional distributions of θ 0 given x 0:t , p(θ 0 |x 0:t ), is Gaussian with mean µ t = E(θ 0 |x 0:t ) and covariance matrix Σ t = Cov(θ 0 |x 0:t ). The RBPF algorithm produces a sample (µ To forecast X t+k given y 1:t , we aim at computing E(X t+k |y 1:t ). If we fix θ 0 and x 0:t , and run model (7) for k time steps, we obtain a value for X t+k as a function f k (x 0:t , θ 0 , ξ), in which ξ indicates the sequence of increments ∆W s , s = t + 1, . . . , t + k.

Using f k (x 0:t , θ 0 , ξ), the conditional expectation is E(X t+k |y 1:t ) = f k (x 0:t , θ 0 , ξ)p(ξ)p(θ 0 |x 0:t )p(x 0:t |y 1:t ) dξdθ 0 dx 0:t (9) where conditional independence of θ 0 on y 1:t given x 0:t allows for substitution of p(θ 0 |x 0:t , y 1:t ) by p(θ 0 |x 0:t ). Then, if for each i we draw θ (i) 0 from p(θ 0 |x (i) 0:t ) and ξ (i) from the distribution of the Wiener process increment, the predictive expectation of X t+k is approximated by

To visualize how well the SIR model fits the observations, we need a way to compare the simulated trajectories with the observed data. We have supposed that the observations are smaller than the true value of the state X, therefore we have to scale them by a factor that makes them comparable to the filtered state. A scaling factor is suggested by building a prediction interval of the state X at each observation time. Note that from (8) the random variable Y /X is a pivotal quantity with beta distribution and we may state that

where u q 2 and u 1− q 2 are the q 2 and the 1 − q 2 percentiles of the beta distribution of U . Then, the corresponding prediction interval for X, after observing y, is

and a natural scaling factor for a point prediction of X is the median of U , producing y/u 0.5 . The feature of (12) is that it does not depend on the SIR modelling assumption, but only on the observation error assumption, therefore it offers a way to see how well the SIR dynamic follows the (transformed) data.

To check the convergence of the method we fix a parameter value and simulate the observations (or data).

We start from an initial condition of 1% infected and 0.1% removed. We simulate data for the parameters β 0 = 0.3 and γ 0 = 0.1. The parameters σ and η in (2) are 0.03 and 0.01, respectively.

We run model (4) with an underlying time step of 1/24 day to generate 67 daily step states. The first 60 states will be used for the estimation procedure and the other 7 to check the goodness of the forecast. Then, we use equation (8) with a = 10 and b = 40 for the beta distribution of the observation error. This distribution is reported in the left panel of Figure 1 and the simulated observations are in the right panel of Figure 1 .

We apply the RBPF algorithm described in Appendix A with 200,000 particles and with a time step of 1/24 day. Since we have a single observation for each day, to run the algorithm we have to impute new hourly observations by means of a linear interpolation between two consecutive daily observations. The imputed observations are no longer indipendently distributed conditionally on the states, however this approximate procedure keeps the effective sample size of the RBPF algorithm at large values with no appreciable difference on the results. We denote the estimates of β 0 and γ 0 with information up to time t byβ t andγ t , see (A-5). Their time behaviour is shown in the left panel of Figure 3 . The behaviour ofR 0 (t), the estimated basic reproduction number (A-6), is displayed in the right panel of Figure 3 .

We denote the filtered or forecasted state asÎ t andR t , where the value of t determines whether we are filtering or forecasting, that is, if our observation period ends at time s, then when t > s,Î t andR t are forecasts; otherwise they are filtered states. From the RBPF we get the filtered states and, using (10), we get a forecast of the dynamics.Î t andR t are compared to the true state in the left panel of Figure 2 , where we see that both the fit up to day 60 and the forecast on days 61-67 are satisfactory. In particular, the forecast well represents the trend of the state.

In reality the true states are unobserved, and we can only compare the filtered states to the observations, by taking under-detection into account. Therefore, we compute daily prediction intervals for infected I t and removed R t , as in (12) 

In this section we analyze the variability of the filtered states and of the estimated parameters due to the variability of the data generated from the system (4), in order to get an impression of how far they can get from the true values, even if the true random under-reporting error distribution is used.

We use the same parameters of Section 4.1, that is, (β 0 , γ 0 ) = (0.3, 0.1), σ = 0.01, η = 0.03, initial values µ 0 = (0.5, 0.5) T and Σ 0 = diag(0.05, 0.02) and initial state X 0 = (1%, 0.1%).

We generate 500 dynamics from the system (4), from which we obtain 500 trajectories of observed infected and removed people. Then, we run the RBPF algorithm on every simulated dataset, with a time step of 1/24 day and for 200,000 particles. For every simulations we compute the trajectorieŝ

where we recall thatÎ t andR t are the estimated infected and removed individuals filtered by the RBPF algorithm, while I t and R t are the true states for the corresponding simulation.

For the sake of a clear representation, in Figure 5 we show one every five trajectories ofÎ t /I t (left In both panels of Figure 5 , after a transient phase with larger dispersion,Î t /I t andR t /R t end up fluctuating around 1, with a stable dispersion in the left panel and a decreasing dispersion in the right panel. Now consider the sum of the root mean square error (RMSE) betweenÎ t and I t and of the RMSE betweenR t and R t for all the trajectories, as a measure of distance between the estimate and the truth. Among all the trajectories we represented the one with the smallest and the one with the largest distance in the left and in the right panels of Figure 6 , respectively. The latter picture shows that the fit may be very unsatisfactory.

Finally, we report the scatter plot of all the pairs (β t ,γ t ) obtained in the 500 simulations ( Figure 7) at t = 60. We observe that they are dispersed around the true value of the parameter (0.3, 0.1). The pair corresponding to the trajectory of minimum distance (green dot) is closer to the true parameter (red dot) than the pair estimated from the trajectory of maximum distance (yellow dot). Then a good fit to the state trajectory is associated with a better estimate of unknown parameters, provided the remaining parameters, which have been assumed as known, are correctly assigned, as we will see in the next section.

An interesting feature of Figure 7 is that the ratioβ t /γ t shows a smaller variability thanβ t and γ t , around a straight line with slope close to three, the true value of R 0 . Therefore we expect to be able to estimate the basic reproduction number with better accuracy and precision than the infection and removal rate parameters.

A statistical model, belonging to a parametric family, is called identifiable if, for any two different values of parameters, there exists at least a measurable set in the sample space that is not assigned Figure 6 ).

The yellow dot represents the pair of parameters corresponding to the filtered state trajectory with maximum distance from the true states (right panel of Figure 6 ). The line is the least squares fit of β 60 againstγ 60 (slope 2.78).

the same probability by the two members of the family, that is, given θ 01 = θ 02 , there exists at least one set B such that

where Y 1:t denotes a finite-length trajectory of observations from (8) . For a deterministic model, this

property is called structural identifiability, which holds if there exists a map from the parameter to the output θ 0 → y 1:t (θ 0 ) which is injective, that is, given θ 01 = θ 02 , the two models y(θ 01 ) and y(θ 02 ) describe different output trajectories. By a differential algebra approach, [30] show that the following deterministic SIR model with its output

defined for a non-normalised population of size N , is structurally identifiable with respect to the unknown parameters β 0 , γ 0 and K. The parameter K > 1 accounts for under-reporting of infected and recovered and has the same function of the U random variables in (8) . Then, after adding noise to the output they go on to show that, despite structural identifiability, the parameters may not be practically identifiable, that is, a good or acceptable agreement between observations and fit is displayed for different values of the parameters when observations end before reaching the peak.

The way randomness has been included into this problem via model (7)-(8) is different from [30] , but we also observe practical identifiability. The problem of identifiability is also discussed in [17] for stochastic SIR models. We generate state trajectories composed by 30 daily values using model (4) with parameters β 0 = 0.1 and γ 0 = 0.03. Then, to obtain the actual observations we multiply each value by a number drawn from a beta distribution with parameters a = 10 and b = 40 (blue line in Figure 8 ). After running the RBPF algorithm using the same beta distribution we obtain results analogous to those in the previous sections, that is, a satisfactory fit of the observed dynamics and a good estimate of the parameters β 0 and γ 0 . Now we run the RBPF algorithm assuming an observation error with beta distribution with a mean larger than the truth, with parameters a = 10 and b = 30 (red line in Figure 8 ). We compare the filtered states with both the true ones and the observed data. First, we consider 500 simulations and look at the ratio between the filtered state and the state. For the sake of a clear representation, in Figure 9 we show one every five trajectories for both the infected and the removed individuals. These ratios are, generally, less than 1 denoting an underestimation of both infected and removed.

Then, we consider the ratio between the filtered states and the adjusted observations. The results of one every five trajectories are reported in Figure 10 for both infected and removed individuals. These ratios fluctuate around 1, denoting a satisfactory fit to the observed data. Figure 11 shows 

In this section we use data of the epidemic in Italy, collected by Protezione Civile (Civil Protection Department) from 24th February. As a first step we consider the data for the whole Italy from 1st

March up to 26th November 2020. Available data are the number of infected, dead, and recovered individuals. Removed people can be obtained by summing dead and recovered people. In Italy, all deaths of people infected with SARS-CoV-2 were classified as COVID-19 ( [31] ). The infected and removed individuals in Italy from 24th February to 26th November are represented in Figure 13 . The total residing population as of 31st December 2019 is 60,244,639 people, as certified by Istat. Figure   11 ). The yellow dot is the case corresponding to the maximum distance case (right panel of Figure   11 ). To deal with underdetection we consider observations as generated by (8), where we recall that U t,1

and U t,2 are independently beta distributed with common shape parameters a and b. In particular, we have I obs = Y 1 = U 1 X 1 and R obs = Y 2 = U 2 X 2 , where the additional notations I obs and R obs are introduced as meaningful names for what follows.

As we can see from Figure 13 the epidemic can be divided in two waves: the first officially began on 24th February and lasted until mid-summer, when the number of infected people began to rise again, as in the rest of Europe ( [5] ). The second wave is distinguished from the first also by the increased test capacity. Hence we consider the two waves as different models, with respect to both the SIR parameters and the observation error distribution parameters and we conventionally set the start of the second wave on 1st August.

To fix a and b we refer to the Infection Fatality Ratio (IFR), a fundamental quantity to estimate the seriousness of the epidemic, and to its crude estimate, the Case Fatality Ratio (CFR). The IFR is the ratio of COVID-19 deaths to total infections of SARS-CoV-2, asymptomatic and undiagnosed infections included, while the CFR is the ratio of COVID-19 deaths to confirmed cases. By the definition, the CFR is greater than the IFR. At any time t an estimate CF R t of the CFR and an estimate IF R t of the IFR are related by the simple relationship

where D t denotes all deaths by time t ( [18] ).

Since by (15)

the ratio IF R t /CF R t can be regarded as the under-reporting factor that we modelled as the U beta random variable introduced earlier.

Given estimates of IF R t as t = 1, . . . , T and the corresponding observed sequence of CF R t we would obtain a sample u 1 = IF R 1 /CF R 1 , . . . , u T = IF R T /CF R T and an estimate of a and b by any established method. Using the method of moments, for example, and considering an estimate of the IFR to substitute IF R t , we would get

whereū and s 2 u are the sample mean and variance of (u 1 , . . . , u T ) andū 1/CF R and s 2 1/CF R are the sample mean and variance of 1/CF R 1 , . . . , 1/CF R T . These equations show that the IFR affects both parameters.

The fatality ratio approach has the advantage that the IFR is a pure number and information on its value can be gathered from different populations. Then, in practice, we may estimate the sample mean and variance of 1/CF R 1 , . . . , 1/CF R T from the observed fatality and removal data, and for a selected IF R assignâ andb. If a range of values is available for the IFR from another source, such as a confidence or a credibility interval, we may repeat the analysis for IF R varying within the interval and evaluate the sensitivity of the results.

For Italy, we may use an indirect method to point at a plausible value of the IFR within this interval, taking advantage of a seroprevalence survey targeting IgG antibody conducted in Italy from May to July by Istat, the Italian national statistical office, and the Italian Health Ministry. Preliminary results obtained from 64,660 people were presented in early August ( [21] ). According to them, almost 1.5 million people in Italy or 2.5% of the population had developed coronavirus antibodies, a figure six times larger than official numbers reported. In short, the idea is to compare the 2.5% figure of people who developed antibodies to the healed people (who have antibodies) estimated from the filtered stateR t in an appropriate time interval. The infected compartment may also contain seropositive individuals, still, the fraction of people in this compartment had become small when Istat's survey started, so we consider only the recovered compartment. The reasoning behind this comparison is that if the assumed IFR is correct, then the observation error distribution derived from (17) is correct and the filtered states are realistic and they should be in agreement with the Istat survey result.

To be more specific, let R t = H t + D t , where H t and D t are the fractions (over the population) of healed and dead people by time t, respectively. Healed people can be seronegative if IgG antibodies are no longer in their system, but we can safely assume that a person enters the healed record soon so they s/he can be considered as seropositive when they do. Now, H t includes all healed since the start of the epidemic, hence a fraction of H t can be seronegative, depending on the duration d of seropositivity. Hence we should compare 2.5% to

values of H t are unknown. We may recover them fromR t and the available data on the fraction of deaths asĤ t =R t − D t /u 0.5 , where u 0.5 is the median of the distribution of the observation error.

Since Istat's survey was carried out between 25 May and 15 July 2020, we compare 2.5% tō

A plausible value of d is three months ( [15] ). This procedure rests on several assumptions and we only regard it as a way to check for gross deviations of our model from reality. indicates a point estimate of IFR of 0.68% (0.53%-0.82%) with high heterogeneity, and suggests that in many populations the IFR would be > 1% if excess mortality was taken into account.

For the first wave, we computed a and b from (17) for a range of IFR values from 0.1%-6%,

where the minimum value is the lowest we found in the relevant literature and has been suggested as lower bound of IFR in Europe by CEBM. The maximum value is still inspired by CEBM and by the considerations in [28] . Indeed, in Italy an estimated initial CFR of about 11-19% has been reported ( [11] , [35] , [31] , [9] ). This suggested to consider as a possible maximum initial value an IFR=6%, accounting for the lack of knowledge at the beginning of the first wave. In particular, we considered 0.1%, 0.35%, 1.3%, and 6%. Moving too far from the highest value gives a large discrepancy between Istat's 2.5% estimated seropositivity in the population and (18). Then we present results for IF R = 5%, for whichH = 2.55%. Then a = 11.72 and b = 81.39 and the corresponding beta density is shown in the right panel of Figure 14 . The initial condition for I t (R t ), for each trajectory, is given Figure 14 where also the prediction intervals for both infected and removed are reported. The prediction intervals are computed from (12) .

The plots ofβ t andγ t from (A-5) are in the top panel of Figure 15 and the plot ofR 0 (t) from (A-6) is in the bottom panel. The gap between the thick and the thin red lines in Figure 14 shows that a single SIR is not able to correctly describe the behaviour of the true dynamics. Therefore, we split the first wave into subintervals. For the partition we consider the DPCMs 1 with the greatest impact on social organization allowing for 10 days for the DPCM to have an effect on the epidemics (that is, change-points are the DPCM dates plus 10 days). In particular, we consider the following DPCM dates: 11th March, 22nd

March, 26th April and 3rd June, so the change-points are on 21st March, 1st April, 3rd May and 13rd

June.

For each time interval, except for the first, we use the filtered statex t (A-4) at the end of the previous interval as initial state and the values ofβ t andγ t at the end of the previous interval as starting parameters. Then the discontinuity in the update is determined only by the initial covariance matrix.

Since in the case of a single time interval we observed that the entries of Σ 0 are updated to very small values, then for the case of several intervals we do not restart each interval with Σ 0 = diag(0.05, 0.02), but compromise between the updated Σ 0 at the end of the previous interval and the initial Σ 0 , taking Σ 0 = diag(0.002, 0.001), for all the time intervals. This choice also allows us to avoid big jumps in the trajectories of the parameters at the change-points. The dynamics of the five different SIR models are represented as a whole dynamics in Figure 16 .

The trajectoriesβ t ,γ t andR 0 (t) show jumps at the change-points (Figure 17) , not very pronounced due to the choice of a small Σ 0 . After a few steps from each jump, the trajectories stabilize following a regular trend. The dynamics of infected individuals fits very well the observed infected divided by u 0.5 . After day 66 (corresponding to 26 April),β t is smaller thanγ t , soR 0 (t) < 1. This value is more realistic thanR 0 (t) in the single time interval case, which is always greater than 1 ( Figure 15 ).

This result is in agreement with the effective reproduction number published for the first time by

Istituto Superiore di Sanità (Italian National Istitute of Health, ISS) on 30 April ( [23] ): the effective reproduction numbers were reported for every Italian region (except for two because of bad quality data) and they were all smaller than one. Now, we analyze the forecast of the infected and removed dynamics for the first wave, computed as in (10) . We consider different cases: we suppose to have observations up to 10 days or 20 days after the second or third change-points, and we try to forecast the dynamics for 7 future days ( Figure 18 ).

It can be observed from Figure 18 that the forecast is satisfactory when it starts 10 days after beta distribution which is increasing with the IFR. Also the width of the prediction intervals decreases as IFR increases. The filtered dynamics well reconstruct the behaviour of the adjusted data in all the cases. The estimated parametersβ t andγ t are very similar for the different cases and only the parameters obtained using IFR=1.3% are reported in Figure 21 with the correspondingR 0 (t).

We observe a big jump inR 0 (t) on the date of the second change point, 10 days after the curfew in Lombardy, followed by a slow decrease, denoting that this measure did not produce the desired effect.

Then it was followed by the measure of 3rd November that allowsR 0 (t) to accelerate its decrease, approaching one, in agreement with what the ISS reported in its 25th November bulletin [22] .

In this section we consider for comparison model (1) combined with the observation equations (8) , that is, the state dynamics is completely deterministic. We repeat part of the analysis done with the stochastic SIR model on the first wave data, using the same observation error distribution, and estimate (β, γ) via maximum likelihood. In Figure 22 we show the single-phase deterministic SIR simulation along with the adjusted observations, in two situations: when the number of infected is descending after the peak and when the descent is almost complete. In both cases the deterministic SIR is unable to capture the dynamics. The single-phase stochastic SIR model (see Figure 14 ), although unsuitable, performs better thanks to the learning mechanism that makes adaptations both to the parameter and the filtered state. The piecewise deterministic SIR model, see Figure 23 , follows the scaled observations more closely, still it is rather less flexible than its stochastic counterpart (see Figure   16 ). This is a further demonstration that a single SIR is unsuitable to describe the true behaviour of the pandemic. updated estimates of key quantities such as the basic reproduction number, which can be acted upon by decision makers. We obtained a rather large basic reproduction number in the initial phase of the first wave, decreasing progressively in the following phases. This value may seem surprising, but other studies, such as [25] , confirm this behaviour which is common in European countries.

The stochastic SIR model might be used to evaluate the effect of mitigation measures by extending the predictions to a horizon of several weeks beyond the date of the next government decree, as an answer to the question: "what would the mid-term evolution of the epidemics have been if this specific measure had not been taken"? But, given the complexity of the phenomenon, which is only partially captured by the model, a great deal of caution is required in doing so. A possible way out is to enrich the model by adding compartments, but, as our model identifiability study has shown, solid prior information or data relevant to the required additional parameters are needed to obtain meaningful results.

With respect to prior information, our approach to the selection of the observation error distribution depends on an estimate of the IFR. For the first wave we have chosen larger values than for the second wave (up to 6% against up to 1.75%), which is in agreement with a report of the ISS [16] , published after we finished our analyis, where the monthly standardised CFR has been calculated.

When standardised with respect to the age and sex structure of the Italian population, the CFR is close to 9% in February-March 2020 and close to 4.5% in April. Then it falls to around 1% in June and July and increases to above 2% in October.

The authors acknowledge the many fruitful discussions with several colleagues, and in particular the colleagues at CNR-IMATI that participated in the COVID-19 modelling study group.

To estimate parameter θ 0 = (β 0 , γ 0 ) T and state X t we propose to apply a Rao-Blackwellized particle filter (RBPF) algorithm. We consider the Euler discretization of the stochastic system (5) reported in equation (7). Since the system is linear in θ 0 , we can apply the Kalman filter. Suppose that θ 0 = (β 0 , γ 0 ) T has a normal prior distribution with mean µ 0 and covariance matrix Σ 0 , then the distribution of θ 0 , given x 0:t+1 = (x 0 , x 1 , ..., x t+1 ) after t + 1 time steps, as t = 0, 1, 2, . . ., is normal with mean µ t+1 and covariance matrix Σ t+1 given by

The distribution of X t+1 given x 0:t is Gaussian with mean B t+1 and covariance matrix G t+1 given by

Recalling that the observations are obtained multiplying the state X t for the beta-distributed observation error term, as defined in equation (8), the RBPF algorithm can be summarized as follows:

STEP 1

• At time t = 0, we draw M initial values of X 0 from its prior distribution π (x 0 ) and obtain M values x (i) 0 , i = 1, 2, ..., M or, alternatively, we put x 0 equal to the initial observation.

• We consider a prior distribution for the parameter θ 0 , given by a normal distribution N (µ 0 , Σ 0 ), where µ 0 is a vector of initial parameters, and Σ 0 is a diagonal covariance matrix.

• To obtain candidate values of the state at importance sampling steps, we will use the distribution implied by the state-transition equation (7) after marginalising it with respect to θ 0 . At step one, a value forX • Denoting by y 1 the observation at time k = 1, we compute weights for each particle from the likelihood atx In order to resample the particles, we need to normalize the weights: R 0 (t) can be regarded as our best estimate of the basic reproduction number in the light of the observed data, not to be confused with the net or effective reproduction number. A very useful discussion on the meaning of R 0 is provided by [13] .

Modeling the Spatiotemporal Epidemic Spreading of COVID-19 and the Impact of Mobility and Social Distancing Interventions

Non-stationary Spatio-Temporal Modeling of COVID-19 Progression in The US. medRxiv

The challenges of modeling and forecasting the spread of COVID-19

On the uniqueness of epidemic models fitting a normalized curve of removed individuals

The Europe second wave of COVID-19 infection and the Italy "strange" situation

COVID-19 infection fatality ratio: estimates from seroprevalence

A Modified SIR Model for the COVID-19 Contagion in Italy

COVID-19: The unreasonable effectiveness of simple models

Interpreting SARS-CoV-2 seroprevalence, deaths, and fatality rate -Making a case for standardized reporting to improve communication

A SIR model assumption for the spread of COVID-19 in different communities

The COVID-19 Infection in Italy: a Statistical Study of an Abnormally Severe Disease

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Complexity of the Basic Reproduction Number (R0)

On sequential Monte Carlo sampling methods for Bayesian filtering

Persistence of IgG response to SARS-CoV-2. The Lancet Infectious Diseases

Il case fatality rate dell'infezione SARS-CoV-2 a livello regionale e attraverso le differenti fasi dell'epidemia in Italia

Simulation and Analysis Methods for Stochastic Compartmental Epidemic Models

Methods for estimating the case fatality ratio for a novel, emerging infectious disease

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Estimation of time-varying reproduction numbers underlying epidemiological processes: A new statistical tool for the COVID-19 pandemic

Primi risultati dell'indagine di sieroprevalenza sul SARS-CoV-2

Epidemia COVID-19. Aggiornamento nazionale25 novembre 2020 -ore 16:00

Bollettino-sorveglianza-integrata-COVID-19_25-novembre-2020.pdf, 2020

Epidemia COVID-19

Bollettino-sorveglianza-integrata-COVID-19_28-aprile-2020.pdf, 2020

A contribution to the mathematical theory of epidemics

The reproduction number of COVID-19 and its correlation with public health interventions

How deadly is the coronavirus? Scientists are close to an answer

Dynamics of COVID-19 using inverse problem for coefficient identification in SIR epidemic models

A systematic review and meta-analysis of published research data on COVID-19 infection-fatality rates

Global Covid-19 Case Fatality Rates

A note on tools for prediction under uncertainty and identifiability of SIR-like dynamical systems for epidemiology

Epidemiological characteristics of COVID-19 cases and estimates of the reproductive numbers 1 month into the epidemic

Estimating the infection and case fatality ratio for coronavirus disease (COVID-19) using age-adjusted data from the outbreak on the Diamond Princess cruise ship

Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in Germany

Modeling Epidemics With Compartmental Models

Clarification of Misleading Perceptions of COVID-19 Fatality and Testing Rates in Italy: Data Analysis

An evaluation of mathematical models for the outbreak of COVID-19

Correction to: System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19

System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19

Extended SIR Prediction of the Epidemics Trend of COVID-19 in Italy and Compared With Hunan

Estimating mortality from COVID-19: scientific brief

A data-driven network model for the emerging COVID-19 epidemics in Wuhan

• We resample M particles from a discrete distribution with support is available.STEP t + 1• For i = 1, . . . , M , sample candidate particlesx • Compute the weightsṽ (i) t+1 for each particle as the product of two distributions with density (A-3). Normalize the weights:• Update the posterior distribution of θ 0 given x are a sample from the distribution of interest. In detail, thet 's are a sample from p(x t |y 1:t ) and, by keeping track of the resampling history, the entire sample x (i) 0:t , i = 1, . . . , M is potentially available, hence a sample from p(x 0:t |y 1:t ). The mean of x (i) t over the particles approximates E(x t |y 1:t ) and we call it the filtered state:The µ The basic reproduction number is defined as R 0 = β 0 /γ 0 , therefore an estimate based on y 1:t is E(β 0 /γ 0 |y 1:t ), which is computed asR