The SIR model [KMW27] is a simple compartmental model that is widely used to model infectious diseases [Het00] . Letting S(t), I(t), and R(t) denote the number of susceptible, infectious and removed (or recovered) individuals at time t, and lettingṠ(t),İ(t), andṘ(t) denote their time derivatives, the SIR model consists in the following three-dimensional continuous-time autonomous dynamical systemṠ (t) = − β N S(t)I(t) (1a)
where N = S(t) + I(t) + R(t) is the constant total population and β and γ are parameters. The SIR model, and several (sometimes deep) variations thereof, have been applied in several works to model the COVID-19 dynamics (see, e.g., [ LGWR20, Atk20, Koz20, Nes20, CNP20, FP20, CFP20]) with known limitations (see [RVHL20, BFG + 20, WF20]). Sometimes, an SIR-like model is used to make long-term predictions (see [BD20] ). However, at the time of writing this paper, it appears that studies are still rare (see, e.g., [SM20] ) where the SIR model parameters and initial conditions are learned on a "train" part of the available data in order to predict a "test" part of the data, making it possible to assess the prediction accuracy of the model. In this paper, we adapt the SIR model to the situation where (i) S, I and R are hidden variables but I(t) is observed through a "proxy" H(t) = αI(t), where α is unknown but constant, and (ii) not only β and γ but also the total population N are unknown and have thus to be estimated. In the context of the COVID-19 application, H will stand for the total number of lab-confirmed hospitalized patients. The proposed adapted SIR model, which we term SH model, is given in (8). It has two state variables (S-a scaled version S-and H) and two parameters (β-which lumps together the parameters β, N , and α-and γ).
We leverage the proposed SH model as follows in order to make hospitalization predictions. Given observed values (H o (t)) t=ti,...,tc , we estimate the parametersβ, γ, and the initial condi-tionsS(t i ) and H(t i ) of the SH model. Then we simulate the SH model in order to predict (H(t)) t=tc+1,...,t f for a specified final prediction time t f . This approach thus combines the areas of parameter estimation (for obvious reasons), data assimilation (for the generation of the initial conditions) and machine learning (for the train-test approach).
In Section 6, we will use a COVID-19 dataset for Belgium 1 that provides us with the following data for t = t s , . . . , t e , where t s is 2020-03-15 and t e is 2020-07-15:
• H o (t): number of COVID-19 hospitalized patients on day t (TOTAL_IN column);
• E o (t): number of COVID-19 patients entering the hospital (number of lab-confirmed hospital intakes) on day t (NEW_IN column);
• L o (t): number of COVID-19 patients discharged from the hospital on day t (NEW_OUT column).
The subscript o stands for "observed". We will also mention results obtained with a dataset for France 2 where t s is 2020-03-18 and t e is 2020-07-17.
In the data for Belgium, there is a mismatch between H o (t) and
. This can be due to the patients who get infected at the hospital (they would be counted in H o without appearing in E o ) and to the patients who die at the hospital (they would be removed from H o without appearing in L o ). In order to remedy this mismatch, we redefine
For the French data, we sum the "rad" (daily number of new home returns) and "dc" (daily number of newly deceased persons) columns to get L o (t). Since there is no column for
Several other COVID-19-related data are available. In particular, the daily number of infected individuals, I o (t), is also reported by health authorities. However, a visual inspection reveals that the graph of I o is less smooth than the graph of H o . A possible reason is that I o is affected by two technical sources of variation: the fraction of tested persons and the accuracy of the tests. In contrast, the reported number of COVID-19 hospitalized individuals, H o , is expected to be much more accurate. Moreover, for the authorities, predicting H is more crucial than predicting I. Therefore, as in [Koz20] , we focus on H.
We assume that, for all t,
where α is unknown but constant over time. In other words, (2) posits that a constant fraction of the infected people is hospitalized. Equation 2 is reminiscent of [CFP20, (3) ], where the number of dead individuals plays the role of H and α is time dependent.
We assume the following observation models with additive noise:
Assuming that the noises are independent Gaussian centered random variables confers a maximum likelihood interpretation to some subsequent estimators, but this assumption is very simplistic.
Multiplying (1a) and (1b) by α, and multiplying the numerator and denominator of (1a) by α, we obtain
LettingS
and using (2), we obtain the simplified SIR model
which we term the SH model. (The "S" in this SH model can be interpreted as the number of individuals susceptible of being hospitalized.) The SH model has only two parameters (β and γ), one hidden state variable (S) and one observed state variable (H) with observation model (3a). Note that, in the SH model (8), the number of patients entering the hospital by unit of time is
and the number of patients leaving the hospital by unit of time is
The goal is now to leverage the SH model (8) in order to predict future values of H based on its past and current observations (H o (t)) t=ts,...,tc . To this end, we have to estimate (or "learn") four variables, which we term estimands: the two parametersβ and γ and the two initial valuesS(t i ) and H(t i ), where t i is the chosen initial time for the SH model (8). One possible approach is to minimize some error measure between the simulated values (H(t)) t=ti,...,tc and the observed values (H o (t)) t=ti,...,tc as a function of the four estimands. However, the error measure is not available as a closed-form expression of the four estimands, and this makes this four-variable optimization problem challenging. We now show that it is possible to estimate H(t i ) and γ separately. This leaves us with an optimization problem in the two remaining estimandsβ andS(t i ), making it possible to visualize the objective function by means of a contour plot. te] , and this data cannot be used to estimate the variables and simulate the SH model. The SH model is initialized at t i , and we refer to the data ( [ti,tc] as the train set, though it is legitimate to widen it to t ∈ [t s , t c ].
To recap, we have t s ≤ t i < t c < t e . The provided dataset goes from t s to t e . The test set is (H o (t), E o (t), L o (t)) t∈[tc+1,H o (t), E o (t), L o (t)) t∈
It is reasonable to believe that H in (3a) is small in practice. Hence we simply take
We have L(t) = γH(t), see (10). In view of the observation model (3), we can estimate γ by a ratio of means:γ
Several other estimators are possible, such as the least square estimator, or the total least squares estimator which is the maximum likelihood estimator of γ for the iid Gaussian noise model (3). Note that t i in the expression ofγ can legitimately be replaced by any time between t s and t c . Only data in the test set, i.e., occurring after t c , are unavailable in the variable estimation phase.
We now have to estimate the two remaining estimands, namelyβ andS(t i ). We choose the following sum-of-squared-errors objective function
where the c coefficients are parameters, all set to 1 in our experiments unless otherwise stated. In (11), H(t), E(t) as in (9), and L(t) as in (10), are given by the (approximate) solution of the SH model (8) in which (i) H(t i ) and γ take the values estimated as above, and (ii)β andS(t i ) take the values specified in the argument of φ. In order to compute the required (approximate) solution of the SH model (8), we use the explicit Euler integration with a time step of one day, yielding, for
Now that the objective function φ (also termed "cost function" or "loss function") is defined, we let the estimated (β,S(t i )) be the (approximate) minimizer of φ returned by some optimization solver.
Recall that the time range between t i and t c is the train period and the time range between t c + 1 and t e is termed the test period.
In order to predict the values of H over the test period, we apply the above procedure to estimate the four estimand variablesβ, γ,S(t i ), and H(t i ), and we compute the solution H(t) of (12) for t from t i to t e . The prediction is then (H(t)) t=tc+1,...,te . The discrepancy between (H(t)) t=tc+1,...,te and (H o (t)) t=tc+1,...,te reveals the accuracy of the prediction.
As an alternative to Section 4.4, we now present a method to estimateβ independently. We do not recommend this alternative, but it sheds light on the various forecast accuracies observed in Section 6.
From (8a) and (9), we obtain d dt
.
Hence a possible estimator forβ is
and, from (9), a possible simple estimator for the remaining estimand is S (t i ) = E(ti)
β H(ti)
.
We can now investigate how the error terms in the observation model (3) impact β . We assume throughout that the errors in H o (t + 1) − H o (t) and E o (t + 1) − E o (t) are comparable. Except at the very beginning of the outbreak, E o (t)H o (t) is considerably smaller than (H o (t)) 2 , and thus the second term of (13) drives the error.
Consequently, the estimation ofβ should be the most accurate when E o (t)H o (t) is the largest. This occurs slightly before the peak of H o (t). This means that the estimation ofβ should be the most accurate for a train period slightly before the peak. However, this does not mean that this position of the train period gives the most accurate forecasts, as we will see below.
Let us consider the situation where the train period is located before the peak. Then the estimation ofβ is less accurate, and this impacts S (t i ). At the initial time t i , this does not impact the right-hand term of (8b) in view of the definition of S (t i ). However, an overestimation ofβ will induce an underestimation ofS(t i ) and, in view of (8a), a subsequent even stronger underestimation ofS(t). Hence the first term of (8b) will be underestimated. As a consequence, the peak in H will appear sooner and lower. The case of an underestimation ofβ leads to the opposite conclusion, namely a peak in H that appears later and higher. In summary, the further before the peak the train period is located, the more inaccurate the position and height of the peak is expected to be.
Finally, let us consider the situation where the train period is located after the peak. Then we can make the same observations as in the previous paragraph, except that predicting the peak is now irrelevant. Moreover, we are in the decrease phase, where the first term of (8b) (which involvesβ andS(t)) is smaller than the second term (which does not involve these quantities). Consequently, the possibly large estimation errors onβ andS(t) will only slightly affect the forecast of H(t).
An alternative to Sections 4.2-4.4 is to reconsider (11) as a function of all four estimands:
In (14), H(t), E(t) as in (9), and L(t) as in (10), are given by the solution of the discrete-time SH model (12) where the parametersβ and γ and the initial conditionsS(t i ) and H(t i ) take the values specified in the argument ofφ. Minimizingφ is a more challenging problem than mimizing φ (11) in view of the larger number of optimization variables. It may be essential to give a good initial guess to the optimization solver, and a natural candidate for this is the values obtained by the procedure described in Sections 4.2-4.4.
In our preliminary experiments, we have found that this alternative does not present a clear advantage in terms of the prediction mean absolute percentage error (MAPE). The results reported in Section 6 are obtained with the sequential prediction approach of Section 4, unless otherwise specified.
We now apply the method of Section 4 (by default) or a method of Section 5 (when specified) to the data of Section 2 available for Belgium (by default) and France (when specified).
The methods are implemented in Python 3 and run with Anaconda 2019.10. The code to reproduce the results is available from https://sites.uclouvain.be/absil/2020.05.
We first check how well the SH model (8) can fit the available data for Belgium. For this experiment, we use the method of Section 5.2 with c E = c L = 0 in order to get the best possible fit (in the least squares sense) to the H o curve. The result is shown in Figure 1 .
The fitting error is remarkably small (MAPE below 15%). For the French data, the fit is even better in terms of MAPE (about 3%).
Note that the parameters of the SH model are constant with respect to time in our experiments. This contrasts with [Koz20] where there are two phases, and with [Nes20] where the infection rate is piecewise constant with several pieces.
We stress that Figure 1 tells us little about the prediction capability of the model. If the fit over some period is bad, then predictions (i.e., forecasts) over that period can only be bad. But if the fit is good (as it is the case here), the predictions can still be bad due to their sensitivity with respect to the data preceding the to-be-predicted period. For example, a better fit (in the RMSE sense) than in Figure 1 can be obtained with a polynomial of degree 8; however, its prediction capability is abysmal.
In order to assess the prediction capability of the model, we have to learn the estimand variables over a train period that we make available to the algorithm, use the learned estimand variables in order to predict H over a subsequent test period (whose data is not available to the algorithm), and finally compare the prediction with the data on the test period. This is what we proceed to do in the rest of this Section 6. Total_in H_train H_pred Figure 1 : Belgium, fitting the SH model to the H o (total hospitalized) curve. In this experiment, the train set is the whole dataset, hence there is no test (prediction) curve. Reproduce with SHR_16PA_py_BEL_1sttraintstart1_1sttraintstop123_c100.zip.
We start with a prediction experiment where the train period is around the peak. According to Section 5.1, this is a promising location. A contour plot of the objective function φ (11) is given in Figure 2 . In order to make the minimizer easier to visualize, the plot shows equispaced-level curves of log(φ − 0.99 φ * ), where φ * is an approximation of the minimal value of φ. Based on a visual inspection, we choose (1e-5,1e4) as the initial guess of the optimization solver. The optimization solver is scipy.optimize.fmin with its default parameters.
The middle plot of Figure 2 shows (H o (t)) t=ts,...,te (observed hospitalizations, gray solid line), (H(t)) t=ti,...,tc (hospitalizations given by the model over the train period, blue dashed line), and (H(t)) t=tc+1,...,te (hospitalizations predicted over the test period, in red dash-dot line). In order to give a sense of the sensitivity of the results, we superpose the curves obtained for three slightly different train periods. The test MAPE for the three curves are 27%, 7%, and 8%.
The right-hand plot of Figure 2 shows the evolution ofS(t).
In Figure 3 , we superpose the results obtained with various train periods of 14 days. The figure corroborates the comments of Section 5.1. In particular, if the train period is fully located before the peak, then the predictions are rather inaccurate. Placing the train period around the peak gives excellent prediction results. When the train period is fully located in the decrease phase, the estimation ofβ andS(t i ) is seen to be very sensitive, but this does not affect much the quality of the prediction of H(t). Figure 4 is the counterpart of Figure 3 for France. These experiments are also compatible with the comments of Section 5.1. We also considered some departments separately, with similar results.
A disconcerting aspect is the evolution of the estimated γ as a function of the location of the train period. In Figure 3 (Belgium), the estimation of γ is grouped around 0.08 for several train periods. However, in Figure 4 , the estimation of γ keeps decreasing, indicating that the daily number of patients leaving the hospital is an increasingly small fraction of the number of patients at the hospital.
The experiments in Section 6 have shown that the proposed method has a remarkably good fitting capability over the whole available data, and also a remarkably good predictive value over certain time ranges for the Belgian data. However, there are also time ranges where the prediction is very inaccurate, and the accuracy is also found to be lower for the French data. The predictions returned by the model should thus be taken with much caution. In keeping with this warning, we refrained from displaying predictions beyond the end time of the datasets. The Python code is freely available to make such predictions, but there is no warranty on their accuracy.
Another source of caution is that we cannot rule out the situation where the considered objective function would be multimodal. The optimization solver might thus get stuck in a local nonglobal minimum, yielding a suboptimal fit of the train data and possibly a poorer prediction than what an omniscient solver would achieve. Moreover, even if the objective function is unimodal, the stopping criterion of the solver may trigger before an accurate approximation of the minimum is reached.
If the proposed model is used to guide prevention policies, then further caveats are in order.
We have seen that the estimation ofβ is very sensitive. Hence the proposed model can hardly help assess the impact of prevention measures onβ. Without knowing sufficiently accurately the impact of prevention measures onβ, we may not aptly use the model to predict their impact on the evolution of the hospitalizations. Yet another caveat is that it may be tempting to deduce from the excellent fit with a constantparameter model (Figure 1 ) that the evolution of the prevention measures over the dataset period has had no impact onβ. But the deduction is flawed. Indeed, in view of the comments made in Section 5.1, the available data could also be very well explained with fairly large jumps inβ during the decrease phase.
In spite of all these caveats, the hospitalization forecasts returned by the method, and also the evolution ofS(t), might be of practical use in the context of various disease outbreaks, e.g., for resource planning. To this end, it will be important to understand which specific features of the COVID-19 outbreak in Belgium made it possible to forecast so accurately the hospitalization decrease several months ahead.