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:
- 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}$;
- 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.
- 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$?
- 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?
