Data Forecasts of the Epidemic COVID-19 by Deterministic and Stochastic Time-Dependent Models

According to WHO [1] report, as of 7:06pm CEST, 15 September 2022, there have been 607,745,726 confirmed cases of COVID-19, including 6,498,747 deaths, globally. The quick spread of this epidemic makes a worldwide impact on not only the human health but the economics and developments. Even worse, the uncontrollability of asymptomatic infections causes a problem that we may underestimate the transmission of the disease. [2] estimates that the pooled percentage of asymptomatic infections among the confirmed cases was $40.50\%$ globally. Nevertheless, tough tactics, such as lockdown and traffic halt, to block the spread of the disease are not long plans because these require large social cost. Fortunately, vaccinations help the immunity and lower the probability of getting infected or dead from this disease [3]. As of 12 September 2022, a total of 12,613,484,608 vaccine doses, reported to WHO [1], have been administered. Hence, it is important to consider the contribution of vaccinations to make more precise policy decision.

In [4] and [5], traditional epidemic models, such as SIS, SIR and SEIR models, and mathematical models of some well-known pandemics, such as smallpox, HIV and malaria, have been discussed in detail. In [6] and [7], the authors considered asymptomatic infections and the quarantine policy when studying the transmission of SARS in 2003. [8], [9] and [10] also investigate mathematical models of COVID-19 considering also the asymptomatic infections. Especially, [9] and [10] discretize the differential equations to difference equations and assume that some parameters are time-dependent. They use the data from different countries, such as China, the USA and Brazil, and apply a popular machine learning method: the Finite Impulse Response (FIR) filters to track the time series of these parameters and then predict the desired values, such as the number of infections and recoveries. Furthermore, [11] considers vaccinations against COVID-19 and uses an Ensemble Kalman Filter to conduct data forecast.

In contrast of deterministic models, stochastic models are more realistic in epidemiology. Many stochastic models of epidemics are also proposed in recent years. The most common two ways of stochastic modeling of epidemics are to conduct a Markov branching process [12], [13] and [14] and to conduct stochastic differential equations [15], [16] and [17]. The quantities in interest shall be probability distributions instead of exact values. For example, [18] provides formulae of the final size in its deterministic model, while [19], [20] and [21] give algorithms calculating the final size probability distribution.

The probability distributions of the final size and the maximum size are important topics in stochastic epidemic models. [12], [19] and [22] are good materials for a review of calculations of these two distributions. Many optimized algorithms calculating these distributions have been proposed, such as in [20] and in [21]. It was calculated that these probability distributions are bimodal.

In this stage of the epidemic COVID-19, the quarantine policy is loosen gradually in many countries and so the data of quarantine is more unavailable. Hence, we will only consider an epidemic prevention policy: vaccination in our models in this paper. In this paper, we will conduct data forecasts to see whether the epidemic will ease or become more severe in the short term and how large the epidemic will be in the long run.

For the first question, we conduct a deterministic SAIVRD model using ordinary differential equations, which involves six compartments: susceptible people (S), (unconfirmed) asymptomatically infectious cases (A), (confirmed) isolated infectious cases (I), recovered people (R), vaccinated-immune people (V) and deaths (D). In this model, we will analyze the existence and the asymptotic stability of the disease-free and endemic equilibria related to the basic reproduction number $R_{0}$. It is shown that if $R_{0}>1,$ then there is a unique endemic equilibrium which is locally asymptotically stable, and if $R_{0}<1,$ then the disease-free equilibrium is locally asymptotically stable. We also assume the time-dependence of some parameters and track their time series to predict the time-varying parameters and values of confirmed cases $I$ using FIR filters. Using the data set [23], [24] and [25] from the USA for $40$ days from June 30-th 2022 to August 8-th, 2022 to predict the next $10$ days, we get the relative prediction errors of $I$ within $3\%$ in $5$ days and $7\%$ in $10$ days, and the trend of change is closed to the real data of $I$. This shows that our proposed model meets the real situation.

For the second question, we conduct a time-dependent Markov SARV model with no demography assuming that there will not be another outbreak and considering the administration of vaccinations. The four compartments are susceptible people (S), asmytomatic infections (A), isolations/recoveries/deaths (R) and vaccine-immune (V). Under this model, we approximate the probability distributions of the final size and the maximum size with one initial infectious person using the time-varying parameters obtained from the first forecast.

The desired distributions are also bimodal in our model with time-independent parameters in our simulations of smaller populations. For the model with time-dependent model, what are different are that there is another peak in the final size distribution, that the probability of ‘‘minor outbreak’’ [15] is higher with time-dependent parameters and that the maximum size distribution is oscillating. The distributions are similar between different populations and so we can estimate the “ratio” of the final size and the maximum size distributions by observing those in small populations. Under recent transmissibility of this disease in the USA, when an initial infection is introduced into all-susceptible large population, ‘‘major outbreak’’ [15] occurs with around $95\%$ of the population; with high probability the epidemic is maximized to around $30\%$ of the population. Moreover, to estimate the extinction probability of this epidemic in a large population, we assume the infinite population. As shown in [21] and [26], the plateau value of cumulative probability of the final size in a large population is approximately the extinction probability. It is also interesting to note that this probability is approximately the reciprocal of the basic reproduction number $R_{0}^{\mathrm{SARV}}$ (defined in the SARV model) whenever $R_{0}^{\mathrm{SARV}}>1$, which demonstrates numerically the result in [15]. These data forecasts can serve as reference to public health policy.

We organize the article as follows. We derive our deterministic SAIVRD model, define the basic reproduction number $R_{0}$ and investigate the existence and asymptotic stability of the disease-free and endemic equilibria in section 2. Section 3 uses the FIR filter to train the time-varying rates to predict the trend of the epidemic in the short term. The probability distributions of the final size and the maximum size in the Markov SARV model with both time-independent and time-dependent parameters will be calculated in section 4. Section 5 provides numerical results. Section 6 states our concluding remarks.

We firstly introduce our deterministic SAIVRD model in this section. The six compartments are defined in section 1. The model is under the following assumptions: there must be a period that an infected case is not confirmed yet which is asymptomatically infectious; asymptomatically infectious cases will not enter $R$ or $D$; a dead person cannot infect others; the system is closed, that is, there is no migration (the border is well-controlled); confirmed cases are well isolated so that they have very low adequate contact with others; nosocomial infections and re-infections are omitted.

The parameters are set as following and to be non-negative:

Let $N_{0}$ be the initial total population. The flow chart is shown as following.

Then the model equations are

with the initial conditions $S(0)=S_{0}>0,A(0)=A_{0}>0,I(0)=V(0)=R(0)=D(0)=0.$ Clearly, we have the result in theorem 2.1.

Since the number of deaths can be calculated by $D=N-S-A-I-V-R$, the model can be reduced to have five variables $S,A,I,V$ and $R$.

We now derive the basic reproduction number $R_{0}$ in our SAIVRD model. The Jacobian matrix of the reduced model from eq. 1 at an equilibrium $(S^{*},A^{*},I^{*},V^{*},R^{*})$ is given by

whose characteristic polynomial is

At the disease-free equilibrium (DFE) $(\displaystyle\frac{B}{v+\mu},0,0,\frac{v}{\mu}\frac{B}{v+\mu},0),$ the Jacobian of eq. 1 is

whose eigenvalues are $-(v+\mu),-(\rho+\gamma+\mu),-\mu,-\mu,\displaystyle\frac{B\beta}{N_{0}(v+\mu)}-(w+\mu).$ Hence, we have $\dfrac{B\beta}{N_{0}(v+\mu)}-(w+\mu)<0$ if and only if the DFE is asymptotically stable. Hence, if we put

then we have theorem 2.2.

We then define the basic reproduction number by eq. 3 and explain why this definition of $R_{0}$ is reasonable in epidemiology as follows. Firstly, the value of $R_{0}$ has positive correlations with $B$ and $\beta$, which means that there are more people will be infected when more people are born or when more people are contacted by infectious people. Secondly, $R_{0}$ has negative correlations with $w,v$ and $\mu$, which means that if we can improve the efficiency of disease screening (to shorten the period of infecting others), improve the efficiency of vaccinations or the infected people die fast, then fewer people would be infected.

It is also worthwhile to note that the rates of recovery $\gamma$ and death $\rho$ of the confirmed cases are not involved in $R_{0}$. This is because the class $I$ is well-controlled such that people in the class $I$ cannot cause further infections.

There may be an equilibrium in which $A^{*}$ and $I^{*}$ are positive in eq. 1, whose stability means that the human will coexist the virus finally and which is called an endemic equilibrium. The existence and asymptotic stability of the endemic equilibria are stated in theorem 2.3.

In this section, we assume the time-dependence of some parameters in our SAIVRD model derived in section 2 to conduct the predictions of infections and time-varying rates in the short term. Since the data is updated in days, it is reasonable to consider the discrete-time model instead of the ordinary differential equations.

Suppose that the parameters $B$ and $\mu$ are constants and that the other parameters $\rho(t),\gamma(t),w(t),v(t)$ and $\beta(t)$ are time-dependent. Then we have the model difference equations in eq. 4.

We will use the known data $S(t),A(t),I(t),V(t),R(t)$ and $D(t)$ for $t=0,1,\cdots,\\ T-1$ to track these five time series: $\rho(t),\gamma(t),w(t),v(t)$ and $\beta(t).$ Then we predict the values of $S(T),A(T),I(T),V(T),R(T),D(T),\rho(T-1),\gamma(T-1),w(T-1),v(T-1)$ and $\beta(T-1).$

For $t=0,1,\cdots,T-2,$ we solve $\rho(t),\gamma(t),w(t),v(t)$ and $\beta(t)$ in eq. 5.

In our numerical experiment in section 5, some $\beta(t)$ are negative, which is epidemiologically unreasonable. In this case, we reset $\beta(t)$ to be $0$.

We will apply a data set in the USA from [23], [24] and [25] that is of large number of infections and so the Finite Impulse Response (FIR) filters is an appropriate method to predict the time-varying rates [9]:

where $J_{i},i=1,2,3,4,5$ are orders of the five FIR filters and $x_{j},y_{j},z_{j},u_{j},b_{j}$ are coefficients of the impulse responses of these five FIR filters. The coefficients can be determined by considering the following minimization problems.

using rigid regression, where $\alpha_{i},i=1,2,3,4,5$ are regulation parameters, which avoid overfitting. After predicting the values of $\hat{\rho}(T-1),\hat{\gamma}(T-1),\hat{w}(T-1),\hat{v}(T-1)$ and $\hat{\beta}(T-1),$ we plug these data back into eq. 4 to get the desired values $\hat{S}(T),\hat{A}(T),$ $\hat{I}(T),\hat{V}(T),\hat{R}(T)$ and $\hat{D}(T)$ by replacing $t$ with $T-1$.

After doing prediction for one day, we can treat this predicted day as known data to predict more one day, and step by step we calculate more days. More precisely, we will predict the values of $\hat{\rho}(t),\hat{\gamma}(t),\hat{w}(t),\hat{v}(t),\hat{\beta}(t),\hat{S}(t+1),\hat{A}(t+1),\hat{I}(t+1),\hat{V}(t+1),\hat{R}(t+1)$ and $\hat{D}(t+1)$ for $t=T-1,T,\cdots,T+d-1.$

In this process, we still need to determine the orders $J_{i},i=1,2,3,4,5$ and find the formulae to solve the optimization problems in eq. 7. Given the orders, we can use a result from [10] named ‘‘Normal Gradient Equation’’ and described as follows, to compute the desired coefficients. Let $f(0),f(1),\cdots,f(T-2)$ be the known data, and the FIR filter $\hat{f}(t)=\sum_{j=1}^{J}x_{j}f(t-j)+x_{0}$ be the predicted data, where $t=J,\cdots,T-2$. Define the function $F$ by

where $\alpha$ is the regulation parameter. If $\det(\alpha I+A)\neq 0,$ then

where $I$ is $(J+1)\times(J+1)$ identity matrix,
$A=\begin{bmatrix}(T-J-1)&\sum_{t=J}^{T-2}f(t-1)&\cdots&\sum_{t=J}^{T-2}f(t-J)\\ \sum_{t=J}^{T-2}f(t-1)&\sum_{t=J}^{T-2}(f(t-1))^{2}&\cdots&\sum_{t=J}^{T-2}f(t-1)f(t-J)\\ \vdots&\vdots&\ddots&\vdots\\ \sum_{t=J}^{T-2}f(t-J)&\sum_{t=J}^{T-2}f(t-J)f(t-1)&\cdots&\sum_{t=J}^{T-2}(f(t-J))^{2}\end{bmatrix}$ and
${\bf{b}}=\begin{bmatrix}\sum_{t=J}^{T-2}f(t)\\ \sum_{t=J}^{T-2}f(t)f(t-1)\\ \vdots\\ \sum_{t=J}^{T-2}f(t)f(t-J)\end{bmatrix}$.

We give the process of calculation in algorithm 1 to algorithm 3 [10]. Firstly, we use algorithm 1 to find the best choice of orders. Then after some pre-processing of data, we use algorithm 2 to predict the next-day data of $S,A,I,V,R,D$ and the time-varying rates. Finally, we predict the data for several days in algorithm 3.

The procedure of algorithm 1 is described as follows. Firstly, we input the known data (the time-dependent parameters), and decide how large $\ell_{T}$ the training set we want to extract from is, bounds $U_{J}$ and $L_{J}$ and regulation parameters. Then we use different $J$’s to find the best $J$ which gives the least residual.

We next conduct short-term prediction in algorithm 2. Due to the limitation of available data from [25], we need to make more modifications and assumptions as follows. Firstly, the cumulative confirmed cases $\tilde{I}$ and cumulative vaccinations $\tilde{V}$ provided in our data set [25] is different from the definition of $I$ and $V$ in our model. Our definition of $I$ is the simultaneous confirmed cases, and $V$ is the number of valid immunity by vaccination. To adapt our definition, we update $I$ to be $(\tilde{I}-R-D)(1-\mu)$ and update $V$ to be $\tilde{V}\sigma(1-\mu),$ where $\sigma$ is the efficiency of vaccination. Secondly, the asymptotic infections $A$ is updated by $A=pI,$ where $p$ is the estimated ratio between asymptomatic and symptomatic infections. Thirdly, $S(t)$ is updated by $S(t)=(N_{0}-A(t)-I(t)-V(t)-R(t)-D(t)+tB)(1-\mu)$ for $t=0,1,\cdots,T-1.$

After these settings, we can predict the data for one day by algorithm 2. Firstly, we calculate the time-varying rate by known data and reset $\beta(t)$ if it is negative. Then use algorithm 1 to find the best FIR orders for prediction and calculate $\operatornamewithlimits{argmin}\limits_{{\bf{x}}\in\mathbb{R}^{J+1}}F({\bf{x}},\alpha)$ to predict the rates. Finally, plug these predicted rates into eq. 4 to predict the data of the six compartments.

When we get the prediction for one day, we can treat them as the known data to predict more days by algorithm 3. Firstly, apply algorithm 2 to predict the data for one day. Then iteratively treat the predicted data as known data to conduct predictions for several days and our goal of this section is accomplished so far.

To estimate how large the epidemic caused by a single infected person is in the long run, we calculate the probability distributions of the final size and the maximum size. The final size is defined to be the cumulative number of the infections until the end of epidemic and the maximum size is defined to be the maximum number of the simultaneous infections during the prevalence of the epidemic. We propose a simplified stochastic SARV model whose four compartments are defined in section 1. We do not consider the demography and assume that there will not be another outbreak, that each asymptomatic infectious case transmits the disease to others according to a Poisson process with parameter $\lambda$, that when a susceptible is infected, it enters $A$, that $A$ enters $R$ (isolated, dead or recovered) after the period $T$ exponentially distributed with parameter $\alpha$ and that $\theta$ is the average ratio of the number of vaccinations to the number of jumps [12].