Fractional Stochastic Differential Equation Approach for Spreading of Diseases

Abstract

The nonlinear fractional stochastic differential equation approach with Hurst parameter H within interval $H\in (0,1)$ to study the time evolution of the number of those infected by the coronavirus in countries where the number of cases is large as Brazil is studied. The rises and falls of novel cases daily or the fluctuations in the official data are treated as a random term in the stochastic differential equation for the fractional Brownian motion. The projection of novel cases in the future is treated as quadratic mean deviation in the official data of novel cases daily since the beginning of the pandemic up to the present. Moreover, the rescaled range analysis (RS) is employed to determine the Hurst index for the time series of novel cases and some statistical tests are performed with the aim to determine the shape of the probability density of novel cases in the future.

1. Introduction

The coronavirus disease emerged at the end of 2019 and since then it has become a global crisis. The very quick spreading of the disease made the study of its spreading of great importance for its control and prevention [1,2,3,4,5,6]. On the other hand, the mathematical modeling for the spreading of diseases has been performed since the works of Daniel Bernoulli in the 17th century for spreading of the small pox. Later, the Bernoulli model has been modified with the aim that unbounded growth does not occur and randomness effects are included in the model. At first, a factor was introduced that had the effect of making the variation rate of population growth N, negative, being, in this case that the model is given by $dN/dt=-\alpha N(1-\nu N)(1-\eta N)$, where the parameters $\alpha$ and $\eta$ are in the range $\alpha >0$ and $0<\eta <\nu$. Moreover, taking some randomness in these parameters, we obtain a more realistic model to handle real situations [7,8,9,10,11,12]. In a general way, whether the size of the population at time t is $N(t)$, being the relative rate of growth of the population at time t, we have the effective population in time t, $N_e(t)$, subject to some randomness due to the environmental effects: $N_e(t)=N(t)+\xi (t)$, where the first term $N(t)$ is deterministic, while the second term $\xi (t)$ is random (reflecting the environmental randomness effect).

The purpose of this paper is the epidemic modeling of new cases of coronavirus using the nonlinear stochastic differential equation approach with fractional white noise and Hurst parameter at range $H\in (0,1)$, thus showing that the fluctuations of officially released COVID data follow non-trivial self-similar and Gaussian statistics, as a fractional Gaussian driving noise would imply. We analyze the rises and falls of novel cases daily of COVID-19 reported by the public healthy agencies in countries as Brazil, where the low number tests performed in the population lead to a great uncertainty in the official data. Due to the possibility that the dynamics of novel cases might exhibit long-range correlations; thus, it may be described in terms of long-memory processes such as the fractional Brownian motion [13]. Since in many problems related to network traffic analysis, finances and many other fields, the processes under study seem empirically to display self-similar properties and the long-range dependent properties and since the fractional Brownian motions are the simplest processes of this kind, it is important to have a systematic study of these processes and to use them to construct other stochastic processes.

The use of standard stochastic differential equations has already been used in the literature to study the behavior of time series for price dynamics in refs. [14,15,16,17]. In addition, the use of standard Brownian motion that corresponds to the case of Hurt index $H=1/2$ was employed recently for studying the spreading of COVID-19 in ref. [3]. The stochastic dynamics smoking with non-Gaussian noise was studied in refs. [18,19,20]. The dynamics of a stochastic COVID-19 epidemic model with jump-diffusion was then analyzed in ref. [18]. The stochastic susceptible-infected-recovered (SIR) dynamics expressed using an ordinary Itô-stochastic differential equation has been proposed in refs. [21,22,23,24,25,26,27]. The statistical analysis described by the master equation and transition rates for the infection process was made in ref. [28]. Moreover, the statistical inference in the epidemic model for the Ebola virus was performed in ref. [29]. The discrete-time epidemic model with binomial distribution has been used to study the transmission rate of COVID-19 in [30]. A statistical model for the health care impact of COVID-19 in India has been studied in [31]. The impact of mobility restriction in COVID-19 superspreading was analyzed in [32]. The effect of immunization through vaccination on the susceptible–infected–susceptible (SIS) epidemic spreading model was analyzed in [33]. Here, we propose the analysis based on the fractional Brownian motion to study the dynamics of novel cases. Moreover, we aim to prove that noise has long memory. The plan of this paper is as follows: in Section 2, we describe fractional white noise analysis as a model for dynamics of the spreading of disease. In Section 3, we present the numerical results for the model and the long range memory by rescaled range analysis (RS). In Section 4, we present our conclusions and final remarks.

2. Fractional White Noise Analysis

The fractional Brownian motion (BMF) is defined as a family of Gaussian processes, being introduced by Kolmogorov in ref. [34] and indexed by the Hurst parameter H into the range $H\in (0,1)$, being the first application made by Hurst in 1951 to model the long-term storage capacity of reservoirs along the Nile River [35,36]. The fractional Brownian motion was defined in refs. [37,38,39], where, for each $H\in (0,1)$, a real-valued Gaussian process $B_H(t)$, $t\le 0$ is defined by $\langle B_H(t)\rangle =0$ and $\langle B_H(t)B_H(s)\rangle =\frac{1}{2}\left[t^{2H}+s^{2H}-{(t-s)}^{2H}\right]$ for all $s,t\in {\mathbb{R}}_{+}$. If $H=1/2$, the fractional Brownian motion is the standard Brownian motion or the Wiener process. We consider H restricted to range $H\in (0,1)$ and apply this equation as well as its modifier for the infected numbers by coronavirus $N(t)$ in each day. Hence, we use the nonlinear fractional stochastic equation with a drift term in the form of a logistic model with a random term to perform the analysis of the spreading of novel cases of coronavirus $N(t)$, given by

$$\begin{array}{c}\hfill dN(t)=\mathfrak{f}(t)dt-\alpha (t)N(t)\left(1-\nu (t)N(t)\right)\left(1-\eta (t)N(t)\right)dt\\ \hfill +\beta (t)\diamond d{B}_H(t),\end{array}$$

(1)

where the drift term, $\mathfrak{f}(t)$, is a better adjustment using a third order polynomial as $\mathfrak{f}(t)$ = −30,827 $+1177t-8.0186t^2+0.02t^3$ in the later considered cases obtained by least squares fit to the data. The add of the least square fit $\mathfrak{f}(t)$ to the logistic model with noise is to take into account the isolation measures and slow vaccination of the populations of each country, which smooths the curve of the growing of novel cases. Due to this, the data are better fitted with a polynomial than fitting directly the solution of the model without noise.

Whether $t\to B_H(t)$ is differentiable into the Schwartz space $S(\mathbb{R})$ of rapidly decreasing smooth functions, ${(S_H)}_H^{*}$ being the dual of ${(S)}_H$

$$\frac{d{B}_H(t)}{dt}={\mathcal{W}}_H(t),\mathrm{in}{(S)}_H^{*},$$

(2)

where $\mathcal{W}{(t)}_H$ is a fractional white noise defined by

$${\mathcal{W}}_H=\sum _{i=1}^{\infty }\left[{\int }_{\mathbb{R}}{e}_i(v)\varphi (t,v)dv\right]{\mathcal{H}}_{\epsilon (i)}(\omega )$$

(3)

where ${\left\{e_n\right\}}_{n=1}^{\infty }$ is the orthonormal basis of $L_{\varphi }^2(\mathbb{R})$ and

$${\mathcal{H}}_{\epsilon (i)}(\omega )={\int }_{\mathbb{R}}{e}_i(t)d{B}_H(t).$$

(4)

The $B_H$-integration theory of Equation (1) is based on using ordinary products path-wise and leads to an integral that, for integrands $g(t,\omega )$, is defined by

$$\begin{array}{c}\hfill {\int }_a^{b}g(t,\omega )\delta B_H(t)=\underset{|\Delta |\to 0}{lim}\sum _{k=0}^{n-1}g(t_k,\omega )\left[B_H(t_{k+1})-B_H(t_k)\right],\end{array}$$

(5)

where $\Delta$: $a=t_0<t_1<···<t_n=b$ is a partition of the interval $[a,b]$ and $|\Delta |={max}_{0\le k\le n-1}\left(t_{k+1}-t_k\right)$. These integrals do not have an expectation zero and are called fractional path-wise integrals. Here, we consider the $B_H$-integral considered in ref. [40], defined by

Definition 1.

Suppose $g:\mathbb{R}\to {(S)}_H^{*}$ is a given function such that $g(t)\diamond W_H(t)$ is integrable in ${(S)}_H^{*}$. Then, we define the fractional stochastic integral of Itô type, ${\int }_{\mathbb{R}}g(t)d{B}_H(t)$, by

$$\begin{array}{c}\hfill {\int }_{\mathbb{R}}g(t)d{B}_H(t):={\int }_{\mathbb{R}}g(t)\diamond d{W}_H(t).\end{array}$$

(6)

For instance, suppose

$$\begin{array}{c}\hfill g(t)=\sum _{k=1}^n{g}_k(\omega ){\chi }_{[t_k,t_{k+1})}(t),g_k\in {(S)}_H^{*},\end{array}$$

(7)

where

$$\begin{array}{c}\hfill {\chi }_{[0,t]}(\tau )=\{\begin{array}{c}1\mathit{if}0\le \tau \le 1,\\ -1\mathit{if}t\le \tau \le 0,\mathit{except}t=\tau =0\\ 0\mathit{otherwise}.\end{array}\end{array}$$

(8)

We concentrate on the $B_H$-integral given in [39]

$$\begin{array}{c}\hfill {\int }_{\mathbb{R}}g(t)d{B}_H(t)=\sum _{k=0}^n{g}_k(\omega )\diamond \left[B_H(t_{k+1})-B_H(t_k)\right],\end{array}$$

(9)

where the definition here is an extension of the fractional Itô integral introduced in ref. [37]. Hence, we concentrate on the $B_H$-integral defined by

$$\begin{array}{c}\hfill {\int }_a^{b}g(t,\omega )d{B}_H(t)=\underset{|\Delta |\to 0}{lim}\sum _{k=0}^{n-1}g(t_k,\omega )\diamond \left[B_H(t_{k+1})-B_H(t_k)\right],\end{array}$$

(10)

where ⋄ denotes the Wick’s product [37,38,39,40,41,42,43,44]. These integrals use the Wick’s products and are known as the fractional Itô integrals with reference to the corresponding situation for standard Brownian motion [38,40]. We consider the white noise calculus based on $B_H(t)$, $H\in (0,1)$. Due to the behavior of novel cases, we consider the deterministic part adjusting to the daily novel cases using a suave curve. The aim is simply to provide a smoothed version of the time series using standard methods as the lowest regression. Furthermore, we aim to calculate the ensemble of all possible trajectories using the fractional Brownian motion as is made in Figure 1.

The zigzag behavior in the range of large t values reflects in an increase of the uncertainty in the data and low-number of test performed in the population. For modeling of the zigzag behavior of novel cases, we add a random noise term in the logistic model for the spreading of an infectious number as given by Equation (1) in order to simulate the effect of rises and falls of novel cases [3]. The $\alpha$, $\nu$, and $\eta$ parameters comes from the logistic model with a threshold. We use the values obtained by Bernoulli for smallpox $\alpha =\nu =0.13$; however, we should have $\alpha (t)$ and $\nu (t)$. The range of $\eta (t)$ values is into $0<\eta (t)<\nu (t)$. $\mathfrak{f}(t)$ comes from the adjustment of least squares to the data into the period considered. The daily fluctuations with a weekly cycle of reported cases have been observed in many countries and are due to the diagnostic and data reporting practices [45]. We obtain the time series of the model Equation (1) using the values obtained by Bernoulli given by $\alpha =\nu =0.13$, $\eta =0.063$. Furthermore, we use a low intensity of fractional white noise $\beta$ as $\beta =1.0\times {10}^{-6}$.

3. Numerical Results

We write the BMF model that is a Gaussian process $W_H(t)$, $t>0$ with zero mean and stationary increments and variance given by $\langle W_H^2(t)\rangle =t^{2H}$. The sample path of the BMF model is a fractal curve with fractal dimension $d=1/H$, where $H=1/2$ corresponds to the standard Brownian motion whose increments $\Delta W_H(t)=W_H(t+1)-W_H(t)$ are statistically independent. For $H\ne 1/2$, the increments $\Delta W_H(t)$ are known as fractional white noise, showing long-range correlation as $\langle \Delta W_H(t+h)\Delta W_H(t)\rangle \simeq 2H(2H-1)h^{2H-2}$, for $h\to \infty$. For $H\in (1/2,1)$, the increments are positively correlated, thus exhibiting a persistent behavior. For $H\in (0,1/2)$, the increments are negatively correlated, showing anti-persistence. For $H=1/2$, we write the Wiener increment as $d{W}_H(t)\sim {\left(\sqrt{dt}{R}_G\right)}^{2H}$, where $R_G$ is an aleatory generator number with Gaussian distribution of mean zero and variance ${\sigma }_{W_H}^2=1$. In the following, we write Equation (1) in the discrete form with the time step $\Delta t$ used in the simulations as $\Delta t=0.001$, with t running along chosen units (say seconds, minutes, or days). In Figure 1, we plot the time series of the model Equation (1) for values of noise intensity as $\beta (t)=1.0\times {10}^{-6}{t}^3$. We have the time series of the variation of novel cases oscillating quickly as displayed in the figure. The results of simulations show that value of H of $0.3$ ($H\le 0.5$) is insufficient for good accuracy of the model and the value of $H=0.55$ looks better and $H=0.46789$ may be quite suitable for a best fit if the value of H does not depend on the time, which is probably not the case. We plot the time series of the model for values of Hurst parameter $H>1/2$: $H=0.55$. The results for $H=1/2$ correspond to random walk [3]. The intensity of the value of fractional white noise used in the calculations is such as $\beta (t)=1.0\times {10}^{-6}{t}^3$.

In Figure 2, we plot the numerical result of the half-width of the distribution of $N(t)$ as a function of time t, $\sigma (t)$. We calculate the variance of the distribution where the standard deviation gives an estimation of novel cases in each day t. The range of data considered here is from 14 March 2021 up to 21 July 2021. The difference to the real data is due to the approach used. We calculate the nth order moments ${\mu }_n=\langle {(x-m_1)}^n\rangle$ until the fourth order about the mean or central moments, where we have the following relations: $c_1={\mu }_1$, $c_2={\mu }_2$, $c_3={\mu }_3$, $c_4={\mu }_4-3{\mu }_2^2$. Normalized measures commonly used indicating a deviation from the Gaussian are the kurtosis ${\lambda }_4$, defined as ${\lambda }_4=({\mu }_4/{\sigma }^4)-3$ and the skewness, ${\lambda }_3$. In Figure 3, we display that the behavior of ${\lambda }_4(t)$. ${\lambda }_4(t)$ is numerically calculated by solving Equation (1)). Moreover, the kurtosis relates to the deviation of the tail of the distribution as compared to Gaussian $P(N(t),t)=\frac{1}{\sqrt{4\pi \beta (t)}}{e}^{-{\left[N(t)\right]}^2/4\beta (t)}$, whose solution would correspond to Equation (1) for $\nu (t)=\eta (t)=0$ and $\mathfrak{f}(t)\equiv 0$. The range of negative values obtained for the kurtosis indicates that the shape of the distribution is closest to Wigner’s distribution [11]. Furthermore, at range of small t considered, where the kurtosis is close to zero, we have that the distribution is closest to a Gaussian (${\lambda }_4=0$).

Rescaled Range Analysis

There are several estimators for the Hurst exponent, such as rescaled range (RS) analysis and multifractal and detrended fluctuation analysis [46,47] (DMFA) and (DFA) that differ basically by the choice of the fluctuation measure. The obtained exponent by DMFA and DFA is similar to the Hurst exponent, except that both DMFA and DFA may also be applied to signals who have underlying statistics. In Figure 4, we present the Log-Log graphic with the objective of determining the Hurst exponent. The window of data considered for analysis were the last 120 days before 9 January 2021. We use the RS analysis, where the Hurst exponent is used as a measure of long-term memory of time series. It relates to the autocorrelations of the time series and the rate at which these decrease as the lag between pairs of values increases. The RS analysis is the oldest and best-known method to estimate the Hurst indexes. In general, the RS analysis is an estimator of the Hurst exponent for random walks; however, in this nonlinear logistic dynamic, the rescaled range is a reliable exponent. The start point of the method is to determine the standard deviation using the mean field approach as $\sigma =\sqrt{\langle N^2\rangle -{\langle N\rangle }^2}$ and the range R as

$$\begin{array}{c}\hfill R={\left[\sum _{t=1}^n{N}_t-n\langle N\rangle \right]}_{max}-{\left[\sum _{t=1}^n{N}_t-n\langle N\rangle \right]}_{min}\\ \hfill t=1,2,...,n.\end{array}$$

(11)

The rescaled range series (RS) is given by

$$\begin{array}{c}\hfill (R/s)=\frac{R_t}{\sigma },log\left(\frac{R_t}{\sigma }\right)=Hlog(t),\end{array}$$

(12)

where H is the Hurst exponent. We obtain the value of the Hurst exponent using the RS method as $H=0.4678(9)$. The presented numerical results show the approximation of the Hurst parameter H close to $H=1/2$, which is related to the classical Brownian motion and common distribution shape parameter as moments, kurtosis, and skewness. Although the value obtained for the data set is very close to $0.5$, we should have H time dependent, $H(t)$, that is, the estimator is relevant in the current context and may be time dependent. The value found for the data set is within the interval $H\in (0,0.5)$, thus indicating an anti-persistence for the time series of case numbers $N(t)$ reported daily.

4. Conclusions

In brief, we analyzed the spreading of novel cases of coronavirus using the nonlinear stochastic differential equation with fractional white noise (nonlinear fractional Brownian motion approach) to model the rises and falls of novel cases of coronavirus in countries such as Brazil. The variability in the time series, for instance in Brazil, is mostly due to the case reports not being filed on weekends, leading to a large-amplitude oscillation within a period of seven days. Furthermore, the uncertainty in the results generated by the low number of tests performed in the populations and to underreporting generates a large uncertainly in the official results and, consequently, the addition of randomness in the differential equations becomes necessary for the growth of infected numbers. We find a small deviation in the Hurst indexes from $H=1/2$ as $H\approx 0.47$, indicating a small deviation from the standard Brownian motion approach. An investigation into whether fractional Gaussian noise is indeed the best (phenomenological) explanation for the fluctuations of COVID data can be made in a future work. The starting point would be more established schemes such as DFA, with as few free parameters as possible, and a much larger database, i.e., COVID data from many countries. However, we must be careful in using fractional, i.e., Gaussian noise with algebraic correlation, and further investigate in a data-driven manner what correlation structure the fluctuations display across different data sets.