We consider state and parameter estimation for compartmental models having both time-varying and time-invariant parameters. Though the described Bayesian computational framework is general, we look at a specific application to the susceptible-infectious-removed (SIR) model which describes a basic mechanism for the spread of infectious diseases through a system of coupled nonlinear differential equations. The SIR model consists of three states, namely, the three compartments, and two parameters which control the coupling among the states. The deterministic SIR model with time-invariant parameters has shown to be overly simplistic for modelling the complex long-term dynamics of diseases transmission. Recognizing that certain model parameters will naturally vary in time due to seasonal trends, non-pharmaceutical interventions, and other random effects, the estimation procedure must systematically permit these time-varying effects to be captured, without unduly introducing artificial dynamics into the system. To this end, we leverage the robustness of the Markov Chain Monte Carlo (MCMC) algorithm for the estimation of time-invariant parameters alongside nonlinear filters for the joint estimation of the system state and time-varying parameters. We demonstrate performance of the framework by first considering a series of examples using synthetic data, followed by an exposition on public health data collected in the province of Ontario.

§ INTRODUCTION

This paper presents an algorithm for the combined state and parameter estimation of stochastic compartmental models having both time-varying and time-invariant parameters. The methodology described herein was originally motivated by system identification in nonlinear dynamical systems with a specific interest in mechanical systems [1, 2, 3].  Bisaillon et al. [4] describe a fully Bayesian model selection framework which leverages the use of nonlinear filters and Markov Chain Monte Carlo (MCMC) methods for the combined state and parameter estimation of dynamical systems with both time-varying and time-invariant parameters. Even though the methodology is inspired by mechanical systems, it operates on general nonlinear dynamical systems, which naturally lends itself to infectious disease modelling using compartmental models. It is particularly relevant due to the fact that many system parameters in a compartmental model may change over the course of a given epidemic or pandemic, whereas others will tend to remain constant (time-invariant). The complimentary use of nonlinear filters and MCMC provides a robust computational framework which systematically alleviates many limitations of traditional methods for time-varying parameter estimation using just nonlinear filters.

There are many examples of the Kalman filter and other nonlinear filters being used for tracking case counts
[5, 6, 7, 8, 9]. When accounting for model discrepancy using stochastic compartmental models, Mamis & Farazmand [10] suggested that the use of correlated noise should be favoured over white noise. Notably, the mechanism by which we estimate the time-varying system parameters are mathematically analogous to a stochastic compartmental model with a coloured noise process (as explained later). Traditional methods for combined state and parameter estimation using nonlinear filters operate by augmenting the state vector to include the system parameters. By appending the time-varying parameters to the state vector, these parameters can be jointly estimated alongside the system states [11, 12]. In order to ensure filter convergence, the parameters are perturbed by an artificial random noise, ensuring that some uncertainty is maintained in the parameter estimates as time progresses. While this  artificial noise resolves the issue of filter divergence, it introduces three new points of concern. First, by treating the system parameters as states in the augmented system, the degree of nonlinearity in the state space model is increased, which in practice may limit the modeller's ability to use certain filters [2]. Second, by permitting both time-varying and time-invariant parameters to vary in time for the purpose of estimation, artificial (or non-physical) dynamics are introduced into the system. Third, by perturbing the parameters by an artificial random noise process, the covariance of this system noise contains additional parameters that must be tuned or estimated in a robust fashion.

The Bayesian computational framework described here alleviates each of these concerns through the combination of sampling methods and nonlinear filters. Relying on MCMC for the time-invariant parameter estimation means that the state vector need only be augmented to include the original system states and time-varying parameters. By only appending the subset of parameters that are known to (or assumed to) vary in time to the state vector, the degree of nonlinearity in the augmented state space model will typically be reduced vis-a-vis the case where all parameters are modelled in this fashion. Furthermore, the parameters that are known to be time-invariant behave as such within the state estimation procedure. Hence, the strength of the artificial dynamics introduced by the estimation procedure is reduced. Finally, the question of estimating the process noise covariance matrix is critically addressed by this setup. By treating the parameters of the covariance matrix as time-invariant parameters to be jointly estimated alongside the set of time-invariant system parameters by MCMC, we provide a robust, fully Bayesian approach to estimate the artificial noise strength as dictated by the system dynamics and data. The caveat of this approach is the increased computational overhead associated with MCMC. The time-scale of infectious disease spread and the frequency of data collection is on the order of days compared to mechanical systems which may be on the order of seconds. Hence, for infectious disease modelling, where uncertainties in the states and parameters are abundant, the increased time required for computation is justified.

In addition to addressing some shortcomings of existing methods for combined state and parameter estimation using filtering methods, this framework critically allows us to directly address some major criticisms of the use of SIR models for long-term predictions of infectious disease spread. Some researchers have discussed the inadequacy of the SIR model to capture the complexity of a global phenomenon like COVID-19 [13, 14], whereas others have discussed the pitfalls of more complex models regarding the identifiability of parameters and the observability of additional compartments [15, 16]. These studies essentially point out that the complexity of the model and availability of data should be considered hand-in-hand to select optimal models capable of meaningful predictions. For aggregated population-level modelling, the inclusion of a time-varying infection rate parameter, estimated within a robust fully Bayesian computational framework addresses these sources of concern concurrently. Note that the definition of a meaningful prediction itself is predicated on the desired spatial and temporal scales, sources of uncertainty and other modelling assumptions [17]. Where the simplicity of a deterministic SIR model with constant parameters is mathematically incapable of generating multiple waves of infection, an SIR model augmented by a time-varying infection rate permits such dynamics to be captured. Furthermore, in addition to the time-varying nature of the disease transmission, other sources of model discrepancy are also captured through this parameter, which is akin to a model error term [12]. Another prominent issue impacting our understanding of the course of the pandemic and our ability to forecast the future trajectory of case counts is the unknown portion of cases that are undetected despite widespread testing efforts [18, 19, 20]. Such predictions are critical to informing policy decisions regarding interventions and healthcare resource allocations based on predicted case counts and hospitalizations [21, 22, 23]. To address these concerns in modelling efforts, some researchers have attempted to systematically address the under reporting of cases by introducing explicit compartments that differentiate between detected and undetected cases [24, 25, 26]. In our case, we have chosen to parameterize the measurement operator, introducing a multiplication factor to be estimated which quantifies the proportion of active cases are actually reported. While this multiplication factor is treated as a time-invariant parameter in the current investigation, it could be treated as a time-varying parameter within the framework to capture potential fluctuations in testing rates.

Section <ref> outlines the mathematical methods employed in this study. The SIR model and its discrete state space form are described in Section <ref>. The Bayesian computational framework and its implementation are described in Section <ref>. The results are presented in Section <ref>. Section <ref> considers model calibration using a series of synthetic datasets that mimic observed trends in public health data collected and reported at the population level during the COVID-19 pandemic. Using known time-invariant parameters and defining three unique time-varying infection rate parameters representing i) random variations in the infection rate, ii) gradual seasonal variations in the infection rate, modelled by a sinusoidal function, and iii) a sudden and significant change in the infection rate to simulate a lockdown, modelled by a sigmoid function. These cases permits the assessment of the algorithmic performance vis-a-vis its ability to recover the ground truth parameter values from which the data were generated. We then apply the same SIR model and inference setup to publicly reported active case data [27] from the province of Ontario as a case study in Section <ref>.

§ MATERIAL AND METHODS

This section contains a description introducing the compartmental models used in this study, and the available public health data that will be essential to the estimation of time-varying and time-invariant parameters. Subsequently, the Bayesian algorithm for parameter estimation is described for a general state-space model.

§.§ SIR Model

The SIR model consists of three compartments, namely, susceptible ($S$), infectious ($I$), and removed ($R$), also commonly referred to as recovered or resolved. The three compartments classify the population based on their current disease state. The relationship between the three compartments, and dynamics of how the population interacts with the disease in question can be represented by the flowchart in Figure <ref>. Mathematically, the dynamics are described by the following system of first-order coupled nonlinear ordinary differential equations (ODEs)

\begin{align*}
\frac{dS}{dt} &= -\beta S(t) I(t) \stepcounter{equation}\tag{{\theequation}a} \label{eq:3compartment1} ,\\
\frac{dI}{dt} &= \beta S(t) I(t) - \gamma I(t) \tag{{\theequation}b}\label{eq:3compartment2}, \\
\frac{dR}{dt} &=  \gamma I(t) \tag{{\theequation}c}\label{eq:3compartment3}.
\end{align*}

where the compartments $S$, $I$, and $R$ are model states, and $\beta$ and $\gamma$ are the system parameters. This system of coupled ODEs can be rewritten as a discrete state space model, which can be written as the following system of algebraic equations when using the forward Euler method,

\begin{align}
S_{k+1} &= S_{k} + \Delta t \left( -\beta S_{k}I_{k} \right),
\\
I_{k+1} &= I_{k} + \Delta t \left( \beta S_{k}I_{k} - \gamma I_{k} \right),
\\
R_{k+1} &= R_{k} + \Delta t \left( \gamma I_{k} \right),
\end{align}

the system states at discrete time instance $t_k$ along the computational grid are $S_k$, $I_k$, and $R_k$, $\Delta t$ is the timestep.

The susceptible compartment consists of all individuals that are vulnerable to becoming infected with the disease. The infectious compartment accounts for all the active infections. An individual will move from the susceptible compartment to the infectious compartment if and only if they engage in a disease transmitting interaction with an individual who is infectious. This is represented by the quadratic nonlinear term $\beta S I$, where the parameter $\beta$ reflects the average number of daily interactions an individual is expected to have, and the probability that any one of those interactions between a susceptible and an infectious individual will result in the disease being transmitted. The parameter $\beta$, may vary in time as a result of the implementation of non-pharmaceutical interventions, such as lockdown measures which reduces the number of daily interactions, as well as a result of various public health measures such as mask mandates or vaccination, which reduce the probability that a given interaction will result in the transmission of the disease. The removed compartment in this model is a terminal compartment, as we do not consider reinfection. Furthermore, we account for both recovered and dead individuals within the removed compartment. These two aspects can are generalized in other basic compartmental such as the SIRS model
(e.g., [28, 29]) which allows for reinfection after a period of temporary immunity, and the SIRD model (e.g., [30, 31]) that has two separate terminal compartments differentiating recovered and deceased individuals. In any event, the removal of individuals from the infectious compartment is controlled by the term $\gamma I$, wherein the parameter $\gamma$ represents a rate equal to the inverse of the average duration of infection. This parameter, $\gamma$, will tend to be constant in time. However, for certain diseases, it may still vary due to new strains or variants of a virus emerging, or as treatments become available. For the sufficiently short time-horizon considered in this study, the parameter $\gamma$ can reasonably be taken to be time-invariant.

Three-compartment SIR model representing Eqs. (<ref> - <ref>).

The SIR model may become stochastic by introducing of a model error term. The model error is implemented systematically to ensure conservation [32, 33, 34]. This model error permits unmodelled disease dynamics to be accounted for by the stochastic volatility of the parameters. Alternatively, we would like to model the parameter $\beta$ as a time-varying quantity allowing the volatility in the population-wide disease dynamics to be reflected by recursively estimating it. To permit the estimation of the time-varying parameter $\beta$, we append it to the state vector, augmenting the dimension of the state space and introducing an additional coupled differential equation to the system in Eq. (<ref> - <ref>),

\begin{equation}
\frac{d \beta}{dt} = q_\xi \xi(t)
\end{equation}

where this rate of change of the time-varying parameter $\beta$, is described by a white noise Gaussian process, thereby resulting.in the following discretized system of equations,

\begin{align}
S_{k+1} &= S_{k} + \Delta t \left( -\beta_{k} S_{k}I_{k} \right),
\\
I_{k+1} &= I_{k} + \Delta t \left( \beta_{k} S_{k}I_{k} - \gamma I_{k} \right),
\\
R_{k+1} &= R_{k} + \Delta t \left( \gamma I_{k} \right),
\\
\beta_{k+1} &= \beta_{k} + q_\xi \sqrt{\Delta t}  \varepsilon_{k}.
\end{align}

where $\varepsilon_k$ is a zero mean Gaussian random variable with unit variance. The parameter $q_\xi$ controls the strength of the artificial noise and is an additional static parameter to be estimated. In the form of a general nonlinear state space model, we can summarize Eqs. (<ref>) as

\begin{equation}
\mathbf{x}_{k+1} = \mathbf{g}_k(\mathbf{x}_k, \bm{\Phi}_s,\mathbf{q}_k).\label{eq:model_op}
\end{equation}

where $\mathbf{g}_k$ is the nonlinear model operator, $\mathbf{x}_k = \{x_{1,k},x_{2,k},x_{3,k},x_{4,k}\}^T = \{S_k,I_k,R_k,\beta_k\}^T$ is the augmented state vector, $\bm{\Phi}_s$ is the set of time-invariant parameters, and $\mathbf{q}_k = \{0,0,0,q_\xi\}^T$ is the stochastic forcing.

All of the parameters are to be estimated using aggregated public health data, reporting daily public testing results for COVID-19. The Bayesian framework used herein provides a robust approach for blending measurements with model predictions. The public health data is characterized by the following general measurement equation  [12, 35]

\begin{equation}
\mathbf{d}_{j} = \mathbf{h}_{j} (\mathbf{x}_{d(j)},\bm{\Phi}_s, \mathbf{\epsilon}_{j}) \label{eq:meas}
\end{equation}

where $\mathbf{d}_{j}$ is the noisy observation of the state vector at time step $t_{d(j)}$ and $\mathbf{\epsilon}_j$ is the measurement noise modelled as independent zero-mean Gaussian random vector with a known covariance matrix. For the current investigation, the measurement operator defines a multiplicative Gaussian noise, with a parameter $\rho$, which characterizes the observed data as a partial account of the true number of infected individuals. The noise in the data accounts for false positives and negatives, while the parameter $\rho$ accounts for the fact that the number of recorded positive tests do not reflect a complete observation of the entire population at all time instance $t_{d(j)}$.

The computational grid in Eq. (<ref>) is indexed by $k$, depicted by filled circles ($\textcolor{tabblue}{\bullet}$). The data in Eq. (<ref>) are indexed by $j$, depicted by empty circles ($\circ$).

The simplicity of the SIR model permits the convenient computation of the effective reproduction number. The effective reproduction number is a critical metric which may be interpreted as the average number of new infections that will result from each active infection. This is conveniently computed as a function of the two system parameters [36],

\begin{equation}
R_e(t) = \frac{S(t)\beta(t)}{\gamma}. \label{eq:Re}
\end{equation}

where $S(t) \beta (t)$ indicates the daily number of disease transmitting interactions a given infectious individual will have with susceptible individuals, and $1/ \gamma$ is the average number of days that an individual will be infectious. This metric is a generalized form of the basic reproduction number, where a value of $R_e(t) \geq 1$ indicates an ongoing outbreak [36].

§.§ Framework for estimating time-varying and time-invariant parameters

If we consider an augmented state $\mathbf{x}_k = \{\mathbf{u}_k, \bm{\Phi}_{t(k)}\}$, the joint posterior of the system state, time-invariant and time-varying parameters of model $\mathcal{M}$, conditioned on data $\mathbf{D} = \{\mathbf{d}_0, \hdots,\mathbf{d}_J\}$ can be written as [4],

\begin{equation}
\text{p}(\mathbf{u}_1,...,\mathbf{u}_k,\bm{\Phi}_{t(0)},\hdots,\bm{\Phi}_{t(k)}, \bm{\Phi}_s | \mathbf{D}, \mathcal{M})  = \text{p}(\mathbf{x}_1,...,\mathbf{x}_k, \bm{\Phi}_s | \mathbf{D}, \mathcal{M}).
\end{equation}

This joint posterior pdf can then be decomposed by conditioning the augmented state on the static parameters as

\begin{equation}
\text{p}(\mathbf{x}_1,...,\mathbf{x}_k, \bm{\Phi}_s | \mathbf{D}, \mathcal{M}) = \text{p}(\mathbf{x}_1,...,\mathbf{x}_k |\bm{\Phi}_s,  \mathbf{D}, \mathcal{M}) \text{p}(\bm{\Phi}_s | \mathbf{D},\mathcal{M}).
\end{equation}

where $\text{p}(\mathbf{x}_1,...,\mathbf{x}_k |\bm{\Phi}_s,  \mathbf{D}, \mathcal{M})$ is the conditional posterior pdf of the augmented state and $\text{p}(\bm{\Phi}_s | \mathbf{D},\mathcal{M})$ is the posterior pdf of the time-invariant parameters. The computation of the former involves the use of state estimation, while the latter is obtained using MCMC. For a given joint sample of the static parameters, the state estimation procedure using nonlinear filters permits the forecasting of uncertainty in the augmented state in time between measurements and the updating of the state estimates where the computational grid $k$ and measurement $d(j)$ coincide. The state estimation procedure is briefly explained in Appendix <ref> and <ref>, with details specific to the current implementation using the Extended Kalman Filter.

The time-invariant parameter posterior pdfs are obtained using MCMC. In the current implementation, we use the Transitional MCMC (TMCMC) algorithm [37] to generate samples from the unnormalized time-invariant parameter posterior,

\begin{equation}
\text{p}(\bm{\Phi}_s | \mathbf{D}, \mathcal{M}) = \frac{\text{p}(\mathbf{D} | \bm{\Phi}_s, \mathcal{M}) \text{p}( \bm{\Phi}_s | \mathcal{M})}{\text{p}(\mathbf{D} | \mathcal{M})} \propto \text{p}(\mathbf{D} | \bm{\Phi}_s, \mathcal{M}) \text{p}( \bm{\Phi}_s | \mathcal{M}) ,\label{eq:bayesmodel}
\end{equation}

where $ \text{p}(\bm{\Phi}_s|\mathcal{M})$ is the parameter prior distribution, and $\text{p}(\mathbf{D} |\bm{\Phi}_s,\mathcal{M})$ is the likelihood function. For the purposes of the current investigations, the model evidence $\text{p}(\mathbf{D} |\mathcal{M})$ is simply a normalization factor. Priors may be assigned to the time-invariant parameters in accordance with some known information from clinical studies, or may simply be used to enforce the parameters in the SIR model to be non-negative. For independent measurement noise, the likelihood function in Eq. (<ref>) is evaluated as the product of the likelihood computed at each data point,

\begin{equation}
\text{p}(\mathbf{D} | \bm{\Phi}_s, \mathcal{M}) = \prod_{j=1}^J \text{p}(\mathbf{d}_j | \bm{\Phi}_s, \mathcal{M}), \label{eq:bayesmodel}
\end{equation}

\begin{align}
\text{p}(\mathbf{d}_j | \bm{\Phi}_s, \mathcal{M}) = \int_{-\infty}^{\ \infty} \text{p}( \mathbf{x}_{d(j)} | \mathbf{x}_{d(j)-1}, \bm{\Phi}_s) \text{p}(\mathbf{d}_j | \mathbf{x}_{d(j)} , \bm{\Phi}_s) \text{d} \mathbf{x}_{d(j)}.\label{eq:lik}
\end{align}

which necessitates the use of state estimation [2, 3]. The state estimation can be performed using Bayesian recursive filters, and expressions for the likelihood function in Eq. (<ref>) are available in [3] for various choices of filters.

§ RESULTS AND DISCUSSION

To demonstrate the performance of this Bayesian computational framework, we consider examples using both synthetic data in Section <ref> followed by public health data in Section <ref>.

For all cases, the measurement operator will be derived from daily testing data, which are partial/incomplete measurements of the infectious compartment, $I$, at time index $k = d(j)$ (i.e. $i_k$ or $x_{2,k}$), corrupted by multiplicative noise. This ensures that the noise in the measurements is higher when case counts are higher. Proportioning the noise strength to the number of reported cases is a more intuitive model compared to additive noise, for instance. Here the noisy and incomplete measurements of the infectious compartment are given by

\begin{equation}
\mathbf{d}_j = \rho i_{d(j)} (1 + \epsilon_j), \qquad \epsilon_j \sim \mathcal{N}(0,{q_\epsilon}^2), \qquad j = 1, \hdots, J
\end{equation}

Using synthetic data representing a sparse, noisy, and partial observation of the state allows us to demonstrate the algorithm's ability to concurrently estimate the model parameters and model states (observed and unobserved). The synthetic data is generated using known parameter values for the time-invariant parameters, and using known functions of time for time-varying parameters. Hence, we are able to test the ability of the algorithm to recover the underlying true parameters and states. We generate three sets of synthetic data, which result from i) a randomly varying infection rate in Section <ref> modelled by a realization of a Wiener process (the integral of a white noise Gaussian process)
[38], ii) a seasonally varying infection rate in Section <ref> modelled by a sinusoidal function with a period of one year, and iii) a simulated lockdown in Section <ref> modelled by the sum of two sigmoid functions representing the effect of the sudden enforcement of a lockdown followed by the gradual relaxation of restrictions as observed through monitoring of human mobility patterns during the COVID-19 pandemic [39, 40, 41]. Note that in the case where the infection rate parameter happens to be time-invariant, the generality of the framework permits the accurate estimation of this parameter through the state estimation procedure, though these results are omitted for brevity.

The public health data consists of COVID-19 testing data recorded over a period of one year in the province of Ontario [27]. This real data will inherently embody a combination of the three aforementioned effects that were tested using the synthetic data. It is crucial that the joint estimation of the state, time-varying parameter, and time-invariant parameters be successful when applied to real data, with the caveat that the results cannot be validated as easily as with synthetic data, as the true values of the states and parameters are ultimately unknown.

For the estimation of the time-invariant parameters, we assign parameter prior distributions within the Bayesian framework as listed in Table <ref>. Each of the parameters have a direct physical interpretation, and this prior pdf allows us to impose restrictions on the parameter values. We do not favour any specific value for the ratio of observed infections to active infections, $\rho$ and instead use a uniform distribution to simply enforce a lower bound of a 0% detection rate ($\rho =0$), and an upper bound of 100% ($\rho =1$) with all values within that range being equally likely a priori. The artificial noise strength $q_\xi$ that perturbs the time-varying infection rate parameter $\beta$, and the estimate of the initial mean of $\beta$ must both be non-negative. Hence, we enforce a lower bound of zero through a uniform distribution for these parameters. The upper bound of the uniform distribution is set to 1 for both parameters as it provides broad support for the parameters given their respective scales. Finally, the initial mean of the infectious compartment is assigned a uniform distribution enforcing that the proportion of individuals initially in the $I$ compartment must be no fewer than 0%, and no more than 100%. The average rate of recovery parameter $\gamma$ is taken as a known value of 1/14, and the multiplicative noise strength $q_\epsilon$ is taken as 5%. The noise strength is estimated based on the deviation of the reported daily cases from the weekly averages in public health data from Ontario [27]. Considering a trailing weekly average, the noise strength was as high as 9%, whereas considering a midpoint estimate results in a noise strength closer to 2%, hence we have selected a value intermediate to these two values. Table <ref> also indicates the rate of recovery parameter and the measurement noise strength, which are taken as known and are not estimated here. This choice was based on the reporting of active cases in the real data, which are predicated on an assumed infectious period of 14 days. An active case that does not result in active hospitalization or death is considered to be resolved 14 days after the reported date of symptom onset date. This assumption is maintained in the synthetic noise case for consistency.

Parameter     Symbol      Prior

Ratio of observed infections to active infections     $\rho$     $\mathcal{U}(0,1)$

Perturbation strength in parameter $\beta$      $q_\xi$      $\mathcal{U}(0,1)$

Initial mean estimate of parameter $\beta$     $\mathbb{E}[\beta_0]$      $\mathcal{U}(0,1)$

Initial mean estimate of state $I$     $\mathbb{E}[I_0]$     $\mathcal{U}(0,1)$

Rate of recovery     $\gamma$    1/14

Measurement noise strength      $q_\epsilon$      0.05

Parameter prior pdfs of uncertain parameters.

§.§ Synthetic data

Here we demonstrate the capability of the computational framework in jointly estimating the above time-invariant parameters, system states and time-varying parameter. The three cases considered here seek to capture multiple waves of infection, which a deterministic SIR model with time-invariant parameters is unable to adequately model. In addition to representing three unique characteristics of disease transmission that were prevalent during the COVID-19 pandemic, these distinct forms of a time-varying infection rate are formulated to test different abilities of the Bayesian computational framework.  The first case in Figure <ref>, we see an infection rate whose evolution in time is generated from a Wiener process. The random fluctuations in the infection rate at each individual timestep are relatively small but compounded over time leading to significant variability over the course of the simulated one-year period. The resulting data generated using this infection rate parameter, shown in Figure <ref>, is by comparison remains relatively smooth over time. This will test the ability of the algorithm to track the highly oscillatory nature of the unobserved time-varying parameter through indirect partial observations of the infectious compartment. The two subsequent cases, the infection rate parameter in Figures <ref> and <ref> vary according to prescribed continuous functions in time. Through partial observations in Figure <ref>, the algorithm will be tasked with tracking a large, but gradual variation in the value of the time-varying parameter value. Similarly, from the observations in Figure <ref> a similar change in parameter value will have to be tracked but the simulated lockdown results in a much more sudden change. The ability to track the time-varying changes in the infection rate parameter value hinges critically on the estimation of the optimal strength of the artificial noise, enabled by the Bayesian computational framework.

The prescribed functions of the infection rate $\beta(t)$ (left column), the compartments $S(t)$, $I(t)$, and $R(t)$ (middle column) and the synthetic data (right column). The first row, panels (a), (b), and (c) show synthetic data case 1 with a randomly varying infection rate. The second row, panels (d), (e), and (f) show synthetic data case 2 with a seasonally varying infection rate. The third row in panels (g), (h), and (i) show synthetic data case 3 with a simulated abrupt lockdown with the subsequent gradual removal of restrictions. All synthetic data represents detected active cases, which here accounts for 25% of the infectious compartment, sampled once per day, and corrupted by 5% multiplicative noise.

§.§.§ Synthetic data case 1: random variations

The case of a noisy infection rate illustrated in Figures <ref> and <ref> is presented first. Given that the fluctuations in the time-varying infection rate being given by a Wiener process, in this case alone, the data-generating model is identical to the model being used for inference in Eq. (<ref>). Hence, in this case, there is an underlying known true value for each parameter. In the subsequent cases, where the infection rate is defined by a time-varying function, there will be no such ground truth for the artificial noise strength, $q_\xi$.

Figure <ref>, which presents a matrix of the marginal and pairwise joint posterior pdfs of the time-invariant parameters. The MAP estimates are indicated by a dotted line, which can be visually compared against the true parameter value, indicated by a dashed line. All parameter values are captured within the support of the posterior pdfs, with the MAP estimate being closely aligned with the true parameter value. Furthermore, note that the posterior pdfs appear to be nearly Gaussian, with the ratio of detected-to-undetected cases parameter, $\rho$ exhibiting high correlation with the estimate of the initial mean number of infectious individuals.

All true parameter values are captured within the support of the posterior pdfs, and the MAP values all appear to be all closely aligned with the true parameter values. To provide a numerical comparison against the ground truth value of the parameters, the MAP estimate, the mean, and standard deviation are summarized in Table <ref>.

Marginal parameter posterior pdfs and scatterplots of pairwise joint parameter posterior pdfs. True parameter values are indicated by dashed lines ($\textcolor{tabred}{\blacksquare}$) and MAP estimates are indicated by dotted lines ($\textcolor{tabgreen}{\blacksquare}$).

Parameter .      $\rho$     $q_\xi$      $\mathbb{E}[\beta_0]$      $\mathbb{E}[I_0]$

Ground truth     0.2500     5.000  $\times 10^{-3}$     0.1200     1.000  $\times 10^{-3}$

MAP      0.2519     4.995 $\times 10^{-3}$    	0.1208    0.968 $\times 10^{-3}$

Mean     0.2627      4.957  $\times 10^{-3}$    0.1336     0.945 $\times 10^{-3}$

Standard deviation     0.0130      0.516 $\times 10^{-3}$     0.0289     0.060 $\times 10^{-3}$

Pertinent statistics of the parameter posterior pdfs of time-invariant parameters in Figure <ref>.

In Figure <ref>, we provide the mean state estimate and a three standard deviation uncertainty interval for the infectious compartment $I(t)$ (panel <ref>) and the time-varying infection rate parameter $\beta(t)$ (panel <ref>), plotted against the known ground truth from the data-generating model. We also plot the effective reproduction number, a quantity that is derived using the joint state estimates of $S(t)$ and $\beta(t)$ (panel <ref>). We therefore omit the state estimation results for the susceptible and recovered compartment for brevity. The results shown for the augmented state are conditioned on the maximum a posteriori (MAP) estimate of the time-invariant parameters.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for (a) the infectious compartment, $I$, (b) the time-varying infection rate parameter, $\beta$, and (c) the effective reproduction number $R_e$. The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations.

Critically, the joint estimation of the time-invariant, the time-varying parameters and the system states permits the forecasting of future cases under uncertainty. Figure <ref> illustrates a three-week forecasts of the $I$ compartment, representing the true number of active cases in panel <ref>, and the time-varying infection rate in panel <ref>.

Note that $\beta$ is constant as the mean forecast is obtained by setting the process noise to zero, which can be understood by referring to Eq. (<ref>). The uncertainty (quantified here by the standard deviation bounds), nevertheless grows with time according to Eq. (<ref>), which provides an expression for forecasting the covariance matrix. The resulting uncertainty bounds grow in time, such that they envelope the possible realizations of the randomly varying infection rate. Since the mean estimate of the time-varying parameter is constant throughout the forecasts, the mean estimate of $I$ is given according to the solution of the SIR model with a constant parameter $\beta$. While the agreement of the predictions shown in panel <ref>, hinges on the successful estimation of the $I$ compartment leading up to the data cut-off, it is also continent upon the specific realization of $\beta$. In this case, the future values of $\beta$ happened to be such that the mean estimate exhibits reasonable agreement. The agreement of the mean estimate of $I$ will not be as close for all realizations of $\beta$, however, the uncertainty bounds grow in time, enveloping the number of infections resulting from the possible realizations.

Three week forecast result from day 365 to day 385 for (a) the infectious compartment, $I$ and (b) the time-varying infection rate parameter, $\beta$. The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations. The dashed vertical line indicates the date of the last available observation.

Continuing with this synthetic data case, we intend to highlight the critical importance of the framework's ability to estimate all the model parameters. Through a series of examples, we illustrate the negative impact that imposing incorrect assumptions regarding parameter values has on the joint estimates of the system state, time-varying parameter, and time-invariant parameters. This will be shown through the inability to accurately capture the ground truth in the presence of such assumptions.

The time-varying infection rate parameter $\beta(t)$ was generated from a Wiener process having a noise strength with a magnitude of 0.005. Here we consider the results in the case where rather than estimating this critical parameter through the inference procedure, we relied on manual tuning, selecting sub-optimal values. We illustrate how a low value of 0.001 and a high value of 0.01 affect the results in Figure <ref>. For this illustrative example, we have taken for granted that we know the initial conditions, thus the only time-invariant parameter that remains to be estimated is the ratio of detected to undetected infectious cases, $\rho$. Panels <ref> and <ref> depict the consequences of insufficient perturbation strength for time-varying parameter estimation. The  the state estimates of $\beta$ clearly illustrate how the low artificial noise strength severely limits the ability of the filter to track the true value of $\beta$, with the estimate of $\beta$ being much more gradually evolving than the truth (panel <ref>). Through the joint estimation of the ratio of detected to undetected infectious cases, the model still manages to exhibit reasonable data-fit. However, the resulting effect leads to a severe underestimate of the infectious compartment (panel <ref>). The case where the artificial noise strength is too high is illustrated in panels <ref> and <ref>. Here, the true infectious compartment is tracked well by the state estimation procedure (though being noisy). This noise is even more apparent in the estimates of the infection rate parameter (panel <ref>), which tracks the trend well. However, there is significant noise in the mean estimate, and the plotted three standard deviation uncertainty intervals indicate significant uncertainty in the estimate.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for the infectious compartment, $I$ (left column) and the time-varying infection rate parameter, $\beta$ (right column). The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations. Recall the generating value of the artificial noise strength, $q_\xi$ is 0.005. Panels (a) and (b) have $q_\xi = 0.001$. Panels (c) and (d) have $q_\xi = 0.01$.

To complement this exercise, we now propose a model where we assume we know $\rho$ (the ratio of detected to undetected cases), but decide to estimate the artificial noise strength needed to track the evolution of the infection rate, $\beta$ in time in Figure <ref>. Recall, the data was generated using a 25% detection rate. Hence, we will see how the estimation results change if we were to assume a 100% detection rate, or 10% detection rate. Intuitively, we will expect the assumed 100% detection rate will result in an underestimate of the true number of infections, and conversely, the 10% detection rate will result in an overestimate. This conjecture is confirmed in panels <ref> and <ref>. What is of interest here, is understanding how these assumptions influence the estimation of the infection rate parameter, from which the effective reproduction number is derived. In panel <ref>, we can see that overestimating the detection ratio has a minimal influence on the ability to track the infection rate, whereas in panel <ref>, underestimating the detection ratio results in poor tracking. This occurs because of the estimation of the artificial noise strength. Recall, the generating value of the artificial noise strength $q_\xi$ was $5.000 \times 10^{-3}$. The results in panel <ref>, is obtained from a MAP estimate of $4.351\times 10^{-3}$, whereas the MAP estimate corresponding to the results in panel <ref> was $6.585\times 10^{-2}$ , an order of magnitude larger. The reason for the large value is not because the framework was unable to estimate the parameter, but rather because in this situation, a larger magnitude process noise is required in order to achieve satisfactory data-fit.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for the infectious compartment, $I$ (left column) and the time-varying infection rate parameter, $\beta$ (right column). The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations. Recall the generating value of the ratio of detected to undetected cases, $\rho$ is 0.25. Panels (a) and (b) have $\rho = 1.0$. Panels (c) and (d) have $\rho = 0.1$.

Finally, we consider a series of combinations of incorrect initial conditions for the augmented state vector in Figure <ref>. In panels <ref> and <ref> we have the correct initial value of $\beta_0 = 0.12$, but an overestimate of $I_0 = 0.01$. In panels <ref> and <ref>, once again, we have the correct initial value of $\beta_0 =0.12$, but an underestimate of $I_0=0.0001$. These four panels exhibit similar behaviour to what was depicted in Figure <ref>. The incorrect initial estimates of the infectious compartment adversely impacts the joint estimates of the state $I(t)$ and the ratio of observed to unobserved cases, $\rho$. Conversely, for a correct initial values of  $I_0 = 0.01$, for an overestimated initial value of $\beta_0 = 0.24$ in panels <ref> and <ref>, or an underestimated initial value of $\beta_0 = 0.06$ in panels <ref> and <ref>, the time-varying parameter is capable of quickly correcting for an erroneous initial condition.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for the infectious compartment, $I$ (left column) and the time-varying infection rate parameter, $\beta$ (right column). The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations. Recall the initial conditions of the generating model are $I_0 = 0.001$ and $\beta_0 = 0.12$. Panels (a) and (b) have $\mathbb{E}[I_0] = 0.005$ and $\mathbb{E}[\beta_0] = 0.12$. Panels (c) and (d) have $\mathbb{E}[I_0] = 0.0002$ and $\mathbb{E}[\beta_0] = 0.12$. Panels (e) and (f) have $\mathbb{E}[I_0] = 0.001$ and $\mathbb{E}[\beta_0] = 0.24$. Panels (g) and (h) have $\mathbb{E}[I_0] = 0.001$ and $\mathbb{E}[\beta_0] = 0.06$.

§.§.§ Synthetic data case 2: seasonal variations

The infection rate exhibiting seasonal variation from Figures <ref> and <ref> is now considered. Because the infection rate parameter is now defined by a sinusoidal function, there is no longer an underlying true value of the artificial noise strength, $q_\xi$, perturbing $\beta$. Instead, the value of $q_\xi$is estimated such that it introduces minimal artificial dynamics, while being sufficiently large to permit the nonlinear filters to track the value of $\beta$ as it varies in time. This is implicitly imposed by the computation of the likelihood using (see Eq. (<ref>)).

In contrast to the previous case, the posterior pdfs of the time-invariant parameters shown in Figure <ref> exhibit clear non-Gaussian behaviour. The pdfs of the ratio of detected cases and the initial number of infectious individuals exhibit evident skewness, and increased uncertainty compared to their counterparts in Figure <ref> for the previous case. Nevertheless, the true parameter values are all captured by the posterior pdfs. Note again that there is no ground truth for the artificial noise strength $q_\xi$.
The mean estimates in Table <ref> exhibit slightly better agreement with the ground truth values than the MAP estimates plotted in Figure <ref>.

Marginal parameter posterior pdfs and scatterplots of pairwise joint parameter posterior pdfs. True parameter values are indicated by dashed lines ($\textcolor{tabred}{\blacksquare}$) and MAP estimates are indicated by dotted lines ($\textcolor{tabgreen}{\blacksquare}$).

Parameter .      $\rho$     $q_\xi$      $\mathbb{E}[\beta_0]$      $\mathbb{E}[I_0]$

Ground truth     0.2500     -     0.1551     1.000 $\times 10^{-3}$

MAP      0.2414     3.307 $\times 10^{-3}$    	0.1701    1.186 $\times 10^{-3}$

Mean     0.2554      3.387  $\times 10^{-3}$    0.1610    1.055 $\times 10^{-3}$

Standard deviation     0.0225      0.437 $\times 10^{-3}$     0.0264     0.118$\times 10^{-3}$

Pertinent statistics of the parameter posterior pdfs of time-invariant parameters in Figure <ref>.

In Figure <ref>, the resulting state estimation results performed using the MAP estimate of the time-invariant parameter tracks the gradual change in the infection rate reasonably well. While the tracking of $\beta(t)$ (in panel <ref>) appears to degrade near the end of the simulation, the infectious compartment, $I(t)$, and the the effective reproduction number, $R_e(t)$, are estimated accurately (having low uncertainty in their estimates) as shown in panels <ref> and <ref>, respectively.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for (a) the infectious compartment, $I$, (b) the time-varying infection rate parameter, $\beta$, and (c) the effective reproduction number $R_e$. The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations.

§.§.§ Synthetic data case 3: abrupt lockdown measures

The following data represents a sudden drop in the infection rate parameter caused by introduction of public health measures such as social distancing or lockdown policies. This sudden drop is followed by a gradual return to the baseline parameter value due to the gradual relaxation of the measures. In a similar fashion to the previous example, the infection rate here in the data-generating model is given by an explicit function of time. Due to the significant and sudden change in parameter value, an optimal artificial noise strength value would permit the signal to quickly adjust to the sudden change in the time-varying parameter while minimally inflating the process noise.
The parameter posterior pdfs are presented in Figure <ref>, and summarized in Table <ref>. Note again in this case, due to the skewness of the posterior pdf, the mean exhibits slightly closer agreement with the ground truth of the time-invariant parameters than the MAP estimates. As one would expect, the artificial noise strength $q_\xi$ has a larger value in this example, than was reported for the more gradual seasonal variations of the infection rate.

Marginal parameter posterior pdfs and scatterplots of pairwise joint parameter posterior pdfs. True parameter values are indicated by dashed lines ($\textcolor{tabred}{\blacksquare}$) and MAP estimates are indicated by dotted lines ($\textcolor{tabgreen}{\blacksquare}$).

Parameter .      $\rho$     $q_\xi$      $\mathbb{E}[\beta_0]$      $\mathbb{E}[I_0]$

Ground truth     0.2500     -     0.1400     1.000 $\times 10^{-3}$

MAP      0.2392     5.201 $\times 10^{-3}$    	0.1334    1.034 $\times 10^{-3}$

Mean     0.2487      5.541  $\times 10^{-3}$    0.1419     1.009 $\times 10^{-3}$

Standard deviation     0.0171      0.889 $\times 10^{-3}$     0.0319     0.085$\times 10^{-3}$

Pertinent statistics of the parameter posterior pdfs of time-invariant parameters in Figure <ref>.

For the same reasons explained in the case of the seasonally varying infection rate parameter, the estimates of $\beta(t)$ appear to degrade as time progresses in panel <ref>. The partially observed infectious compartment, and estimated effective reproduction number are both tracked well in time in panels <ref> and <ref>.

State estimation results $\text{p}(\mathbf{x}_1,\hdots,\mathbf{x}_K \vert \mathbf{D}, \bm{\Phi}_s^{\text{MAP}},\mathcal{M})$ for (a) the infectious compartment, $I$, (b) the time-varying infection rate parameter, $\beta$, and (c) the effective reproduction number $R_e$. The ground truth is indicated by a dashed line ($\textcolor{tabred}{\blacksquare}$), the mean estimates are indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the shaded area reflects the mean $\pm$ 3 standard deviations.

§.§ Public health data: Ontario (March 1, 2020 - February 28, 2021)

The province of Ontario publishes a number of COVID-19 metric in their public Data Catalogue. The data used in this example is the reported number of active cases, from the column labelled confirmed positive in the data set [27]. Cases that are classified as active, is based on an assumed infectious period of 14 days. Thus, a case is considered positive for 14 days following the reported date of symptom onset or estimated episode date. We have selected one-year worth of testing data, collected between April 1, 2020 and March 31, 2021, such that the start date of our inference procedure aligns with the date this definition of an active case was adopted. Similarly, the selected end date aligns with the more widespread availability of vaccines in the province, thus pharmaceutical interventions have a minimal effect on the observed trends. Though the province of Ontario administered its first vaccines on December 14, 2020, the widespread distribution of vaccines to the general population did not begin until after the one-year period we consider in this case study. Hence, we intentionally omit the explicit modelling of the vaccination, and include any effect the early doses had on the population-level metrics to be absorbed by the time-varying infection rate parameter.

Each new case will be considered active for the average duration of infection, given by the inverse of the recovery rate parameter $\gamma$. As noted in Table <ref>, the rate of recovery corresponds to an infectious period of 14 days.

Ontario public health data reporting active cases for the 365 days from April 1, 2020 to March 31, 2021.

Marginal parameter posterior pdfs and scatterplots of pairwise joint parameter posterior pdfs. MAP estimates are indicated by dotted lines ($\textcolor{tabgreen}{\blacksquare}$).

Through the examples with synthetic data, we demonstrated the benefits of jointly estimating the state, time-varying parameters using nonlinear filters, and the time-invariant parameters using MCMC. Figure <ref> shows the marginal and pairwise joint posterior pdfs of the time-invariant parameters. As in the case of synthetic data (Figures <ref>, <ref>, and <ref>), the mean initial estimate $\mathbb{E}[I_0]$ and the fraction of observed infections show substantial correlation.

For the case with public health data, there is no known truth to provide a basis for comparison.  Hence, instead, we provide the estimated effective reproduction number $R_e$ in Figure <ref> based on the joint estimates of the susceptible compartment, $S$, and time-varying infection rate parameter, $\beta$ as in Eq. (<ref>). We have superimposed the the estimates published by Ontario Public Health [42]. At first glance the two sets of results appear to be quite disparate. However, comparing the two sets of results relative to an effective reproduction number of $R_e = 1$ (indicated by a dotted in Figure <ref>) the observed trends between the means of the two results are actually qualitatively quite similar in shape. The mean prediction of our results suggest that the reproduction number varies much more dramatically in time, compared to the published results which instead suggest the reproduction number hovered around 1 over the course of the one-year period studied here. Critically, note the difference in the reported confidence interval for the published results, compared to the uncertainty bounds obtained by in the current study. The three-standard deviation uncertainty bounds which envelopes the reported mean prediction almost fully captures the effective reproduction number published by Ontario Public Health . Given the modelling simplifications that are imposed when modelling a global phenomenon using compartmental models, and uncertainties in the system states and parameters, it is critical that modelling approaches adequately reflect the degree of uncertainty in the predictions, particularly when they form a basis for decision-making.

Our effective reproduction number estimate is indicated by a solid line ($\textcolor{tabblue}{\blacksquare}$), and the Ontario Public Health estimates are shown in a dashed line ($\textcolor{tabred}{\blacksquare}$) for 365 days from April 1, 2020 to March 31, 2021. The shaded area for our estimate reflects the mean $\pm$ 3 standard deviations, and the shaded area from Ontario Public Health in orange reflects the reported upper and lower bound confidence intervals. A dotted line ($\blacksquare$) is included included at $R_e = 1$, to assist in comparing the estimated timing of outbreaks.

§ CONCLUSION

Though we demonstrate the performance of the Bayesian computational framework for COVID-19 inspired problems, the framework itself is general. It is portable to other compartmental model structures, and more broadly to other models of dynamical systems. The generality of the framework makes it useful for both retrospective analyses and studies for post-pandemic recovery, as well as a versatile for future pandemic preparedness. Through three illustrative examples using synthetic data, the methodology proved to be adept at recovering the time-invariant and time-varying data-generating model parameters. For the case where the infection rate parameter was modelled by a Wiener process, the algorithm successfully estimated the strength of the artificial noise jointly with the parameter controlling the proportion of detected cases, as well as the mean initial conditions for the augmented state vector. For the cases where the infection rate instead had a known functional form, there is no ground truth value for the artificial noise strength. The estimation procedure bypasses the need for manual tuning of the artificial noise strength, by including this parameter among the time-invariant parameters and the Bayesian framework automatically favours the lowest noise strength that permits the filter to track the observed signal optimally. Even for a simple SIR model, the generality of this method relaxes many modelling assumptions, particularly those related to the initial conditions, which are of critical importance. Applying this validated methodology to public health data, the algorithm was shown to produce qualitatively similar results to those published by public health agencies. For instance, estimating the ratio of observed infections to active infections being at 11.2% with a standard deviation of 0.6%. The estimating this quantity as a time-invariant parameter through the framework leads to the estimates of the basic reproduction number differing in variability relative to the values reported by the agency.

§ STATE ESTIMATION

In this application, the state estimation procedure serves two purposes; to estimate the augmented state of the system (compartments and time-varying parameter) given the proposed model, the set of static parameters, and the incoming noisy sensor data, and to compute the likelihood function for the joint estimation of the time-invariant parameters. As we augment the state to include the time varying parameters, we have now introduced additional nonlinearity into the state space model, and hence, nonlinear filters are required for state estimation. Here we use the extended Kalman filter (EKF) for all cases, noting that this imposes the assumption that the state pdfs are approximately Gaussian, which generally holds for weakly nonlinear systems, complemented by dense/frequent observations [43].

The EKF imposes a Gaussian assumption on the state such that $\mathbf{x}_k \sim \mathcal{N}(\mathbf{x}_k^a,\mathbf{P}_k^a)$, where $\mathbf{x}_k^a$ and $\mathbf{P}_k^a$ are the update (or analysis) mean vector and covariance matrix, respectively. Similarly, the process noise is assumed to be normally distributed with $\mathbf{q}_k \sim \mathcal{N}(\mathbf{0},\mathbf{Q}_k)$. For points along the computational grid, indexed by $k$, that do not coincide with a data point, indexed by $d(j)$, the current mean state estimate and its uncertainty are forecasted according to [35]

\begin{align}
\mathbf{x}_{k+1}^f&=\mathbf{g}_k(\mathbf{x}_k,\bm{\Phi}_s, \mathbf{0}), \label{eq:forecast_x} \\
\mathbf{P}_{k+1}^f&= \mathbf{A}_k\mathbf{P}_k^a\mathbf{A}_k^T + \mathbf{B}_k\mathbf{Q}_k\mathbf{B}_k^T. \label{eq:forecast_P}
\end{align}

where the Jacobian matrices

\begin{align}
\mathbf{A}_k &= \left. \frac{\partial \mathbf{g}_{k}(\mathbf{x}_k^a,\bm{\Phi}_s, \mathbf{q}_k)}{\partial \mathbf{x}_{k}} \right\vert_{\mathbf{x}_k = \mathbf{x}_{k}^a, \mathbf{q}_k = \mathbf{0}} \label{eq:ekf_update1}, \\
\mathbf{B}_k &=  \left. \frac{\partial \mathbf{g}_{k}(\mathbf{x}_{k},\bm{\Phi}_s, \mathbf{q}_k)}{\partial \mathbf{q}_k} \right\vert_{\mathbf{x}_k = \mathbf{x}_k^a, \mathbf{q}_k = \mathbf{0}}  \label{eq:ekf_update2}.
\end{align}

When the computational grid coincides with a measurement, the analysis step in Eqs. (<ref> - <ref>), is performed, allowing for the mean and covariance of the state estimate to be updated. The measurement operator is given in Eq. (<ref>), where the measurement noise is assumed to be Gaussian $\varepsilon_j \sim \mathcal{N}(\mathbf{0},\bm{\Gamma}_j)$, where $\Gamma_j$ is the covariance matrix, resulting in the following equations [35]

\begin{align}
\mathbf{x}_{d(j)}^{a}&=\mathbf{x}_{d(j)}^f + \mathbf{K}_{d(j)}(\mathbf{d}_j - \mathbf{h}_{d(j)}(\mathbf{x}_{d(j)}^f,\mathbf{0})), \label{eq:a5} \\
\mathbf{P}_{d(j)}^a&= (\mathbf{I} - \mathbf{K}_j\mathbf{C}_j)\mathbf{P}_{d(j)}^f, \label{eq:a6} \\
\mathbf{K}_{d(j)}&= \mathbf{P}_{d(j)}^f \mathbf{C}_j^T \left[\mathbf{C}_j \mathbf{P}_{d(j)}^f \mathbf{C}_j^T + \mathbf{D}_j \bm{\Gamma}_{j} \mathbf{D}_j^T  \right]^{-1},  \label{eq:a7}
\end{align}

where the Jacobian matrices $\mathbf{C}_j $ and $\mathbf{D}_j$ are evaluated as

\begin{align}
\mathbf{C}_j &= \left. \frac{\partial \mathbf{h}_{d(j)}(\mathbf{x}_{d(j)}, \bm{\varepsilon}_j)}{\partial \mathbf{x}_{d(j)}} \right\vert_{\mathbf{x}_{d(j)} = \mathbf{x}_{d(j)}^f, \mathbf{\varepsilon}_{j} = \mathbf{0}} \label{eq:ekf_update1}, \\
\mathbf{D}_j &=  \left. \frac{\partial \mathbf{h}_{d(j)}(\mathbf{x}_{d(j)}, \bm{\varepsilon}_j)}{\partial \bm{\varepsilon}_j} \right\vert_{\mathbf{x}_{d(j)} = \mathbf{x}_{d(j)}^f, \mathbf{\varepsilon}_j = \mathbf{0}}  \label{eq:ekf_update2}.
\end{align}

Critically, when data are available, we can also compute the likelihood function in Eq. (<ref>) according to [3, 43, 2, 1]

\begin{align}
\text{p}(\mathbf{d}_j\vert \bm{\Phi}_s,\mathcal{M}_i) =  \mathcal{N}\left(\mathbf{d}_j | \mathbf{h}_{j}\left(\mathbf{x}_{d(j)}^f, \mathbf{0}\right), \bm{\Sigma '}\right)  \label{eq:int_ekf}
\end{align}

\begin{equation}
\bm{\Sigma '} = \mathbf{C}_j \mathbf{P}_{d(j)}^f \mathbf{C}_j^T + \mathbf{D}_j \bm{\Gamma}_j \mathbf{D}_j^T
\end{equation}

§ JACOBIAN MATRICES

The model operator in Eq. (<ref>) can be written in terms of the elements of the state vector $\mathbf{x}_k$,

\begin{align}
x_{1,k+1} &= x_{1,k} + \Delta t \left( -x_{1,k} x_{2,k} x_{4,k}  \right),
\\
x_{2,k+1} &= x_{2,k} + \Delta t \left( x_{1,k} x_{2,k} x_{4,k} - \gamma x_{2,k} \right),
\\
x_{3,k+1} &= x_{3,k} + \Delta t \left( \gamma x_{2,k} \right),
\\
x_{4,k+1} &= x_{4,k} + q_\xi \sqrt{\Delta t}  \xi_{k}.
\end{align}

The above equations lead to the following expressions for the Jacobian matrices,

\begin{equation}
\mathbf{A}_k = { \left[ \begin{array}{cccc}
1 + \Delta t \left( - x_{2,k} x_{4,k}  \right)
\Delta t \left( -x_{1,k}x_{4,k}  \right)
\Delta t \left( -x_{1,k} x_{2,k}   \right)
\\
\Delta t \left( x_{2,k} x_{4,k}  \right)
1 + \Delta t \left( x_{1,k} x_{4,k} - \gamma \right)
\Delta t \left( x_{1,k} x_{2,k} \right)
\\
\Delta t \left( \gamma \right)
\\
\end{array}\right]},
\qquad
\mathbf{B}_k =  { \left[ \begin{array}{ccccc}
0 \\ 0 \\ 0 \\ q_\xi \sqrt{\Delta t}
\end{array}\right]}
\end{equation}

The measurement operator represents a partial measurement of the infectious compartment. The ratio of detected to undetected infections is reflected in the parameter $\rho$,

\begin{equation}
d_j = \rho x_{2,j}\left(1 + \epsilon_j \right)
\end{equation}

The Jacobian matrices of the measurement operator are

\begin{equation}
\mathbf{C}_j =  { \left[ \begin{array}{ccccc}
0&1 & 0 &  0
\end{array}\right]},
\qquad
\mathbf{D}_j = \left[ \begin{array}{c}
\end{array}\right]
\end{equation}

[1]
Mohammad Khalil, Dominique Poirel, and Abhijit Sarkar.
Probabilistic parameter estimation of a fluttering aeroelastic system
in the transitional Reynolds number regime.
Journal of Sound and Vibration, 332(15):3670–3691, 2013.

[2]
Mohammad Khalil, Abhijit Sarkar, Sondipon Adhikari, and Dominique Poirel.
The estimation of time-invariant parameters of noisy nonlinear
oscillatory systems.
Journal of Sound and Vibration, 344:81–100, 2015.

[3]
Philippe Bisaillon, Rimple Sandhu, Mohammad Khalil, Chris Pettit, Dominique
Poirel, and Abhijit Sarkar.
Bayesian parameter estimation and model selection for strongly
nonlinear dynamical systems.
Nonlinear Dynamics, 82(3):1061–1080, 2015.

[4]
Philippe Bisaillon, Brandon Robinson, Mohammad Khalil, Chris L Pettit,
Dominique Poirel, and Abhijit Sarkar.
Robust bayesian state and parameter estimation framework for
stochastic dynamical systems with combined time-varying and time-invariant
arXiv preprint arXiv:2210.08730, 2022.

[5]
Xiaoshu Zeng and Roger Ghanem.
Dynamics identification and forecasting of COVID-19 by switching
Kalman filters.
Computational Mechanics, 66:1179–1193, 2020.

[6]
Francisco Arroyo-Marioli, Francisco Bullano, Simas Kucinskas, and Carlos
Tracking R of COVID-19: A new real-time estimation using the Kalman
PloS one, 16(1):e0244474, 2021.

[7]
Koushlendra Kumar Singh, Suraj Kumar, Prachi Dixit, and Manish Kumar Bajpai.
Kalman filter based short term prediction model for COVID-19
Applied Intelligence, 51(5):2714–2726, 2021.

[8]
Qiwen Sun, T Miyoshi, and Serge Richard.
Analysis of COVID-19 in Japan with extended SEIR model and ensemble
Kalman filter.
Journal of Computational and Applied Mathematics, 419:114772,

[9]
Geir Evensen, Javier Amezcua, Marc Bocquet, Alberto Carrassi, Alban Farchi,
Alison Fowler, Pieter L Houtekamer, Christopher K Jones, Rafael J de Moraes,
Manuel Pulido, et al.
An international initiative of predicting the SARS-CoV-2 pandemic
using ensemble data assimilation.
Foundations of Data Science, 3(3), 2021.

[10]
Konstantinos Mamis and Mohammad Farazmand.
Stochastic compartmental models of the COVID-19 pandemic must have
temporally correlated uncertainties.
Proceedings of the Royal Society A, 479(2269):20220568, 2023.

[11]
Robert H Shumway and David S Stoffer.
Time series analysis and its applications, volume 3.
Springer, New York, 2000.

[12]
Geir Evensen.
Data assimilation: the ensemble Kalman filter.
Springer Science & Business Media, 2009.

[13]
Shiva Moein, Niloofar Nickaeen, Amir Roointan, Niloofar Borhani, Zarifeh
Heidary, Shaghayegh Haghjooy Javanmard, Jafar Ghaisari, and Yousof Gheisari.
Inefficiency of SIR models in forecasting COVID-19 epidemic: a case
study of Isfahan.
Scientific reports, 11(1):4725, 2021.

[14]
Glenn Ellison.
Implications of heterogeneous SIR models for analyses of COVID-19.
Technical report, National Bureau of Economic Research, 2020.

[15]
Gemma Massonis, Julio R Banga, and Alejandro F Villaverde.
Structural identifiability and observability of compartmental models
of the COVID-19 pandemic.
Annual reviews in control, 51:441–459, 2021.

[16]
Luca Gallo, Mattia Frasca, Vito Latora, and Giovanni Russo.
Lack of practical identifiability may hamper reliable predictions in
COVID-19 epidemic models.
Science advances, 8(3):eabg5234, 2022.

[17]
Inga Holmdahl and Caroline Buckee.
Wrong but useful—what covid-19 epidemiologic models can and cannot
tell us.
New England Journal of Medicine, 383(4):303–305, 2020.

[18]
Siuli Mukhopadhyay and Debraj Chakraborty.
Estimation of undetected COVID-19 infections in India.
MedRxiv, pages 2020–04, 2020.

[19]
Jeffrey Shaman.
An estimation of undetected COVID cases in France.
<https://www.nature.com/articles/d41586-020-03513-9>, 2021.

[20]
Claire Rippinger, Martin Bicher, Christoph Urach, Dominik Brunmeir,
N Weibrecht, G Zauner, G Sroczynski, B Jahn, N Mühlberger, U Siebert,
et al.
Evaluation of undetected cases during the COVID-19 epidemic in
BMC Infectious Diseases, 21(1):1–11, 2021.

[21]
Jianfei Zhang, Harini Sanjay Pathak, Anne Snowdon, and Russell Greiner.
Learning models for forecasting hospital resource utilization for
COVID-19 patients in Canada.
Scientific Reports, 12(1):8751, 2022.

[22]
Marietta M Squire, Megashnee Munsamy, Gary Lin, Arnesh Telukdarie, and Takeru
Modeling hospital energy and economic costs for COVID-19 infection
control interventions.
Energy and Buildings, 242:110948, 2021.

[23]
Tetyana Kendzerska, David T Zhu, Michael Pugliese, Douglas Manuel, Mohsen
Sadatsafavi, Marcus Povitz, Therese A Stukel, Teresa To, Shawn D Aaron,
Sunita Mulpuru, et al.
Trends in all-cause mortality and inpatient and outpatient visits for
ambulatory care sensitive conditions during the first year of the COVID-19
pandemic: A population-based study.
Journal of Hospital Medicine, 17(9):726–737, 2022.

[24]
Giuseppe Gaeta.
A simple SIR model with a large set of asymptomatic infectives.
arXiv preprint arXiv:2003.08720, 2020.

[25]
Qing Liu and Longbing Cao.
Modeling time evolving COVID-19 uncertainties with density dependent
asymptomatic infections and social reinforcement.
Scientific Reports, 12(1):5891, 2022.

[26]
KD Olumoyin, AQM Khaliq, and KM Furati.
Data-driven deep-learning algorithm for asymptomatic COVID-19 model
with varying mitigation measures and transmission rate.
Epidemiologia, 2(4):471–489, 2021.

[27]
Government of Ontario.
Status of COVID-19 cases in Ontario.
Accessed: 2023-08-24.

[28]
Guanyu Hu and Junxian Geng.
Heterogeneity learning for SIRS model: an application to the
arXiv preprint arXiv:2007.08047, 2020.

[29]
Mevin B Hooten, Jessica Anderson, and Lance A Waller.
Assessing North American influenza dynamics with a statistical
SIRS model.
Spatial and spatio-temporal epidemiology, 1(2-3):177–185,

[30]
Giuseppe C Calafiore, Carlo Novara, and Corrado Possieri.
A time-varying SIRD model for the COVID-19 contagion in italy.
Annual reviews in control, 50:361–372, 2020.

[31]
Jesús Fernández-Villaverde and Charles I Jones.
Estimating and simulating a SIRD model of COVID-19 for many
countries, states, and cities.
Journal of Economic Dynamics and Control, 140:104318, 2022.

[32]
Elisabetta Tornatore, Stefania Maria Buccellato, and Pasquale Vetro.
Stability of a stochastic SIR system.
Physica A: Statistical Mechanics and its Applications,
354:111–126, 2005.

[33]
Xianghua Zhang and Ke Wang.
Stochastic SIR model with jumps.
Applied Mathematics Letters, 26(8):867–874, 2013.

[34]
Chunyan Ji and Daqing Jiang.
Threshold behaviour of a stochastic SIR model.
Applied Mathematical Modelling, 38(21-22):5067–5079, 2014.

[35]
Andrew H Jazwinski.
Stochastic processes and filtering theory.
Courier Corporation, 2007.

[36]
Fred Brauer, Pauline Van den Driessche, Jianhong Wu, and Linda JS Allen.
Mathematical epidemiology, volume 1945.
Springer, 2008.

[37]
Jianye Ching and Yi-Chu Chen.
Transitional Markov chain Monte Carlo method for Bayesian model
updating, model class selection, and model averaging.
Journal of engineering mechanics, 133(7):816–832, 2007.

[38]
Peter E Kloeden and Eckhard Platen.
Stochastic differential equations.
In Numerical solution of stochastic differential equations,
pages 103–160. Springer, Berlin, 1992.

[39]
Kevin Linka, Mathias Peirlinck, Francisco Sahli Costabal, and Ellen Kuhl.
Outbreak dynamics of COVID-19 in Europe and the effect of travel
Computer methods in biomechanics and biomedical engineering,
23(11):710–717, 2020.

[40]
Kyra H Grantz, Hannah R Meredith, Derek AT Cummings, C Jessica E Metcalf,
Bryan T Grenfell, John R Giles, Shruti Mehta, Sunil Solomon, Alain Labrique,
Nishant Kishore, et al.
The use of mobile phone data to inform analysis of COVID-19
pandemic epidemiology.
Nature communications, 11(1):4961, 2020.

[41]
Jed A Long and Chang Ren.
Associations between mobility and socio-economic indicators vary
across the timeline of the covid-19 pandemic.
Computers, environment and urban systems, 91:101710, 2022.

[42]
Government of Ontario.
Effective reproduction number (Re) for COVID-19 in Ontario.
Accessed: 2022-01-14.

[43]
M Khalil, Abhijit Sarkar, and S Adhikari.
Nonlinear filters for chaotic oscillatory systems.
Nonlinear Dynamics, 55(1):113–137, 2009.