# Parameter Estimation for the SIRD model via Kalman Filter (Part I)

## 0 Introduction

This is my third attempt in tuning a good prediction SIRD model for the COVID-19 outbreak. The model in question is the following: $$\tag{12}\begin{cases} S_{t}&=S_{t-1}-\alpha\frac{S_{t-1}I_{t-1}}{N} \\ I_{t}&=I_{t-1}+\alpha\frac{S_{t-1}I_{t-1}}{N}-\beta I_{t-1}-\gamma I_{t-1} \\ R_{t}&=R_{t-1}+\beta I_{t-1} \\ D_{t}&=D_{t-1}+\gamma I_{t-1} \\ \end{cases} \qquad \text{for} \,\, t=1,2,\dots$$

0.1 first attempt

In my first try I've used for the estimation of the parameters $$\alpha, \beta, \gamma$$ a very simple strategy by searching the least squares parameters that minimize the quadratic cost of the 1-step prediction along all the observation horizon. The complete dissertation of this problem is here (note: in this thread I will use the same notations that I've used in the cited thread). The estimation problem is very easy to resolve, but leads to a result with very poor performances.

0.2 second attempt

As Sextus Empiricus suggests to me, using the 1-step prediction, i.e. $$\tag{13}\hat{y}_t(\theta)=\varphi_t \theta + y_{t-1} \qquad \text{for} \,\, t=0, 1, 2,\dots$$ is not a good idea because, once one has found the parameters, during the simulation phase the model haves to feedback the proper previous predictions. This means, in principle, that the right way to tune the parameters $$\theta$$ is to minimize some suitable cost function of the prediction errors generated by a "batch" predictor of the form
$$\tag{14}\hat{y}_t(\theta)=\hat\varphi(\hat{y}_{t-1}) \theta + \hat{y}_{t-1} \qquad \text{for} \,\, t=0, 1, 2,\dots$$ where $$\tag{15}\hat\varphi(\hat{y}_{t-1}) \triangleq \begin{bmatrix} -\frac{\hat{S}_{t-1}\hat{I}_{t-1}}{N} & 0 & 0 \\ \phantom{-}\frac{\hat{S}_{t-1}\hat{I}_{t-1}}{N} & -\hat{I}_{t-1} & -\hat{I}_{t-1} \\ 0 & \phantom{-}\hat{I}_{t-1} & 0\\ 0 & 0 & \phantom{-}\hat{I}_{t-1} \end{bmatrix} \qquad \text{for} \,\, t=1,2,\dots$$ and $$\hat{S}_{t-1}, \hat{I}_{t-1}$$ are predictions (and not observed values) generated by the model. As a consequence, the new estimation problem $$\tag{16}\theta_\text{LS}\triangleq \arg\min_{\theta \in \mathbb{R^3}} V_T (\theta)$$ with the same quadratic cost $$V_T(\theta)\triangleq \frac{1}{2}\sum _{t=0}^T \|y_t-\hat{y}_t(\theta) \|^2$$ is not easy anymore because now the considered model $$(14)$$ is not linear in his parameter $$\theta$$.

To resolve $$(16)$$, under the suggestion Sextus Empiricus, I've chosen a numerical approach. More precisely, I've used a simple gradient descent equipped with Armijo's adaptive learning rate / step update, with a simple finite-difference approximation of the gradient. The optimization works pretty well, but from the starting guess $$\theta_{\text{LS},0}=0$$ converges to a parameter $$\theta_{\text{LS},\infty}$$ with a non-zero cost that doesn't satisfies me in the final simulation phase.

I've tried several different optimizations by using randomly different initializations $$\theta_{\text{LS},0}$$ (like a MultiStart global optimization policy) with even worse results.

0.3 third attempt

I'm quite sure that the optimizer works well, so I believe that to improve the performance of the final simulator it is necessary to increase the expressivity of the model $$(12)$$. To make a step ahead towards an acceptable solution, under the observations of Sextus Empiricus, now I want make an attempt where I increase the difficult of the estimation problem by relaxing the time-invariance hypothesis of the parameters. Thus, now I consider this more complicated model $$\tag{17}\begin{cases} S_{t}&=S_{t-1}-\alpha_t\frac{S_{t-1}I_{t-1}}{N} \\ I_{t}&=I_{t-1}+\alpha_t\frac{S_{t-1}I_{t-1}}{N}-\beta_t I_{t-1}-\gamma_t I_{t-1} \\ R_{t}&=R_{t-1}+\beta_t I_{t-1} \\ D_{t}&=D_{t-1}+\gamma_t I_{t-1} \\ \end{cases} \qquad \text{for} \,\, t=1,2,\dots$$ The actual "dynamic" estimation problem can be solved by using a Kalman filter.

note: as mentioned by Sextus Empiricus, another possible approach to the problem is to tune the parameters by using the explicit solution of the SIRD model. I don't want follow this way because, as a student in control system engineering, I'm interested in learn how to identify generic dynamical systems in their implicit state-space representations. The SIRD problem is a "toy example" that I'm trying to resolve better as possible in my quarantine free-time.

## 1 Kalman Filtering

1.1 dynamic system model

Consider the following linear, time-variant and autonomous state-space model for a dynamic system $$\tag{18}\begin{cases} x_{t+1} &= A_t x_t+w_t \\ y_t &= C_t x_t+v_t \end{cases}\qquad \text{for} \,\, t=0,1,2,\dots$$ where $$x \in \mathbb{R}^n$$ is the state vector, $$y \in \mathbb{R}^p$$ is the output vector of the system and the parameters $$A\in \mathbb{R}^{n \times n}$$ and $$C \in \mathbb{R}^{p \times n}$$ are matrices. The signals $$w\in \mathbb{R}^n$$ and $$v\in \mathbb{R}^p$$ are white noises characterized respectively by null mean values and covariances $$Q\in \mathbb{R}^{n \times n}$$ and $$R\in \mathbb{R}^{p \times p}$$, i.e. \begin{align} w_t &= \text{wn}(0,Q_t) \\ v_t &= \text{wn}(0,R_t) \end{align}\qquad \text{for} \,\, t=0,1,2,\dots \tag{19}

1.2 Kalman filter definition

The Kalman filter is another dynamic system that estimates dynamically the state of the system $$(18)$$ by watching the model $$(18)$$ and a dataset $$D\triangleq\{y_0, y_1, y_2, \dots\}$$ of observed outputs. The estimation process is a recursion starting from the initial state $$x_0$$, that is unknown and modelled as a random variable with mean value $$\hat{x_0}\in \mathbb{R}^n$$ and covariance $$P_0 \in \mathbb{R}^{n \times n}$$, i.e. $$\tag{20} x_0 \sim (\hat{x}_0, P_0)$$
At every time instant $$t=1,2,3,\dots$$ the recursion is splitted in two different steps:

1. correction: the filter computes the estimate $$\hat{x}_{t|t}$$ of the state $$x_t$$ by considering all the observed past values $$y_0,y_1, \dots, y_{t-1}$$, summarized in a predicted estimate $$\hat{x}_{t|t-1}$$, and the actual observation $$y_t$$. Such estimate $$\hat{x}_{t|t}$$ is characterized by the covariance $$P_{t|t}$$;
2. prediction: the filter computes the estimate $$\hat{x}_{t+1|t}$$ of the state $$x_{t+1}$$ by considering the first equation of the model $$(18)$$. Such estimate $$\hat{x}_{t+1|t}$$ is characterized by the covariance $$P_{t+1|t}$$.

the Kalman filter algorithm is the following \boxed{\begin{align} \text{FOR}\quad &t=1,2,3,\dots :\\ &e_t \triangleq y_t - C_t \hat{x}_{t|t-1} \\ &S_t \triangleq C_t P_{t|t-1}C_t'+R_t \\ &L_t \triangleq P_{t|t-1}C_t' S_t ^{-1} \\ &\text{1) correction}\\ &\hat{x}_{t|t} \triangleq \hat{x}_{t|t-1}+L_t e_t \\ &P_{t|t} \triangleq (I-L_t C_t)P_{t|t-1}(I-L_t C_t)'+L_t R_t L_t'\\ &\text{2) prediction}\\ &\hat{x}_{t+1|t} \triangleq A_t \hat{x}_{t|t}\\ &P_{t+1|t} \triangleq A_t P_{t|t}A_t'+ Q_t \end{align}}\tag{21} note: if $$x_0$$, $$w_t$$, $$v_\tau$$, are gaussians and incorrelated for all the possible choices of $$t, \tau$$ then at every time istant $$t$$ the estimate $$\hat{x}_{t|t}$$ is by definition the MMSE estimate of the state $$x_t$$, so the Kalman Filter finds the optimal estimation of the state $$x_t$$.

note: the signal $$e_t$$ is called innovation and it is the prediction error of the output of the system $$(18)$$ (since the predicted output generated by the Kalman filter is $$\hat{y}_{t|t-1}=C_t \hat{x}_{t|t-1}$$). Here the matrix $$S_t$$ is the covariance of the innovation (and not the number of susceptibles at time $$t$$, I hope that it is clear from the context when I'm talking about the innovation covariance or the number of susceptibles). The matrix $$L_t$$ is called correction gain since it is used to define the corrected estimate $$\hat{x}_{t|t}$$. The matrix $$I$$ is an identity matrix $$n \times n$$.

note: the initialization is $$\hat{x}_{1|0}\triangleq \hat{x}_0$$, $$P_{1|0} \triangleq P_0$$.

1.3 application to the SIRD estimation problem: static estimation

Before to jump in to the more general problem that involves the model $$(17)$$, I want to use the Kalman filter to make a more simpler "static" estimation of the parameters.
In order to use the Kalman filter, the idea is to define the following non-linear, time-variant, static system $$\begin{cases}x_{t+1} &=x_t+w_t \\ y_t &=h_t(x_t)+v_t\end{cases}\qquad \text{for} \,\, t=0,1,2,\dots \tag{22}$$ where $$x_t \triangleq \begin{bmatrix} \alpha_t \\ \beta_t \\ \gamma_t \\ \end{bmatrix} \qquad y_t \triangleq \begin{bmatrix} S_t \\ I_t \\ R_t \\ D_t \end{bmatrix} \qquad h_t(x_t) \triangleq \begin{bmatrix} {S}_{t-1}-\frac{{S}_{t-1}{I}_{t-1}}{N}\alpha_t \\ {I}_{t-1}-{I}_{t-1}\beta_t-{I}_t\gamma_t+\frac{{S}_{t-1}{I}_{t-1}}{N}\alpha_t \\ {R}_{t-1}+{I}_{t-1}\beta_t \\ {D}_{t-1}+{I}_{t-1}\gamma_t \end{bmatrix} \\ \text{for} \,\, t=0,1,2,\dots \tag{23}$$ the system $$(22)$$ is not linear, so the idea to deal with this problem is to consider a simple linearization around the estimate $$\hat{x}_{t|t-1}$$of the output equation (like in an extended Kalman filter). $$y_t \approx h_t(\hat{x}_{t|t-1})+\frac{\partial h_t}{\partial x_t}\Bigg|_{x_t=\hat{x}_{t|t-1}}(x_t-\hat{x}_{t|t-1})+v_t \tag{24}$$ by defining the new coordinates $$\tilde{y}_t \triangleq y_t-h_t(\hat{x}_{t|t-1}) \qquad \tilde{x}_t \triangleq x_t -\hat{x}_{t|t-1} \tag{25}$$ and the time-variant matrix $$C_t \triangleq \frac{\partial h_t}{\partial x_t}\Bigg|_{x_t=\hat{x}_{t|t-1}}=\begin{bmatrix} -\frac{\hat{S}_{t-1}\hat{I}_{t-1}}{N} & 0 & 0 \\ \phantom{-}\frac{\hat{S}_{t-1}\hat{I}_{t-1}}{N} & -\hat{I}_{t-1} & -\hat{I}_{t-1} \\ 0 & \phantom{-}\hat{I}_{t-1} & 0\\ 0 & 0 & \phantom{-}\hat{I}_{t-1} \end{bmatrix} \tag{26}$$ the system $$(22)$$ is approximated by the linear system $$\begin{cases}\tilde{x}_{t+1} &=\tilde{x}_t+w_t \\ \tilde{y}_t &=C_t \tilde{x}_t+v_t\end{cases}\qquad \text{for} \,\, t=0,1,2,\dots \tag{27}$$ at this point it is possible to apply the Kalman filter $$(21)$$ to get an estimation for the parameters $$\alpha, \beta, \gamma$$.

1.4 application to the SIRD estimation problem: dynamic estimation

In order to take into account the time dependency of the parameters it is sufficient to modify the first equation of the system $$(22)$$ in a suitable way that reflects approximately the real dynamic of the parameters
$$\begin{cases}x_{t+1} &=f_t(x_t)+w_t \\ y_t &=h_t(x_t)+v_t\end{cases}\qquad \text{for} \,\, t=0,1,2,\dots \tag{28}$$ in other words, the dynamic of the parameters $$\alpha_t, \beta_t,\gamma_t$$ is described by the motion model $$f_t(x_t)$$. To get an estimate of theese parameters it is possible to use an extended Kalman filter like in the previous section 1.3.

## 2 My Questions

Before to writing down some code to get experimental responses, I'd like to have some suggestions about the previous dissertation. I have also two questions.

1. I'm not sure that $$C_t$$ is well defined. Let's assume that his definition $$(26)$$ is correct. Since the prediction model is $$(14)$$, that is a dynamic model, the quantities $$\hat{S}_{t-1}$$ and $$\hat{I}_{t-1}$$ are functions of all the previous predictions $$\hat{y}_0 (\theta=\hat{x}_{t|t-1}),\hat{y}_1 (\theta=\hat{x}_{t|t-1}), \dots, \hat{y}_{t-2} (\theta=\hat{x}_{t|t-1})$$. This means that at every time step the computational burden to get $$C_t$$ increase because it requires the simulation with the actual parameter $$\theta=\hat{x}_{t|t-1}$$ of the system $$(12)$$ up to time $$t-2$$ (that is increasing). It is possible to define a more efficient recursive formula to define $$C_t$$?
2. How can we define a reasonable motion model $$f_t(x_t)$$ that can take into account the observations made by Sextus Empiricus in the reply in my first thread?

Nice write-up. I think eq#23 is not quite correct though. The observation $$y_t$$ should be a function of the state $$x_t$$ and noise $$v_t$$, like you said in eq#22 $$y_t = h(x_t) + v_t$$.

The issue with eq#23 is that $$y_t$$ depends on $$y_{t-1}$$. This would break the Markov assumption that $$x_t$$ is a sufficient state, which is foundational to the Kalman filter formulation.

I think the correct approach is to use $$x_t := \begin{bmatrix} S_t \\ I_t \\ R_t \\ D_t \\ \alpha_t \\ \beta_t \\ \gamma_t \end{bmatrix}$$

This would be considered "augmenting the state" (SIRD) with the parameters ($$\alpha$$,$$\beta$$,$$\gamma$$) to do a "joint-state-parameter" estimation.

Your dynamic would be $$x_{t+1} =\begin{bmatrix} {S}_{t}-\frac{{S}_{t}{I}_{t}}{N}\alpha_t \\ {I}_{t}-{I}_{t}\beta_t-{I}_t\gamma_t+\frac{{S}_{t}{I}_{t}}{N}\alpha_t \\ {R}_{t}+{I}_{t}\beta_t \\ {D}_{t}+{I}_{t}\gamma_t \\ f_\alpha(x_t)\\f_\beta(x_t) \\ f_\gamma(x_t) \end{bmatrix} + w_t$$

You can treat the SIRD part of the dynamic as deterministic if you so desire by setting its entries of the covariance $$Q_t$$of $$w_t$$ to $$0$$, though that places all the uncertainty on the initial condition and parameters, which I don't think accurately reflects reality - the SIRD model isn't perfect even with best-fit parameters.

Your observation model would be the actual data, which we'd consider to be measurements of the SIRD part of the state, corrupted by noise (which you can remove by making its covariance $$R_t \approx 0$$ if you think the SIRD records are perfect). $$y_t = \begin{bmatrix} S_t \\ I_t \\ R_t \\ D_t \end{bmatrix} + v_t$$

As for what you should design for $$f_\alpha(x_t)$$ and the other parameter dynamics, I think a random walk is fine (i.e. $$f_\alpha(x_t) = \alpha_t$$ so all the movement is due to the $$\alpha$$-component of $$w_t$$). But really I think designing a more sophisticated structure for that is a matter of studying epidemiology more. By creating a structured dynamic for the parameters, you are essentially making a more detailed model than what SIRD alone captures. Such an advanced design should probably be rooted in hard science, else it just becomes excess model capacity ripe for overfitting / poor generalization.

In other words, if you want a model for how the parameters would move, you should come up with an explanation / mechanism for why they would move. What do these parameters mean - what dynamical factors would influence them?

By the way, the selection of the noise covariances is pretty critical to the usefulness of the Kalman filter. If you are using them as tuning knobs rather than estimating them from data, you might as well just go back to your second approach (nonlinear regression of a multi-step predictor) and just use Gauss-Newton's method instead of gradient descent for optimization. The extended Kalman filter is essentially a covariance-weighted recursive application of Gauss-Newton least-squares.