Analyzing a SEIR-Type mathematical model of SARS-COVID-19 using piecewise fractional order operators

1.   Introduction

Modelling real data, particularly in fields like finance and biology, often involves latent randomness. Consequently, using deterministic models, such as those based on ordinary or partial differential equations, to represent this data may lead to inaccuracies. Alternatively, stochastic diffusion processes (SDPs) can be used, as they are designed to account for the random behavior inherent in these data sets. Typically, a diffusion process $X(t)$ is a solution of the stochastic differential equation (SDE)

where $W$ is a Wiener process, $a(t, X(t), \theta),$ $b(t, X(t), \theta)$ are called the drift and diffusion functions respectively, are known functions, and $\theta$ are unknown indexed parameters. Recently, (1.1) has been employed to develop various new homogeneous/non-homogeneous SDPs, including the Log-Logistic ^[1], Schumer ^[2], Lomax ^[3], Weibull ^[4], Square-Brennan-Schwartz ^[5], Brennan-Schwartz ^[6], Gamma ^[7], Gompertz ^[8], and Rayleigh ^[9] diffusion processes. Most of these studies have focused on modeling growth data, such as microorganism culture growth ^[1], the evolution of electricity net consumption ^[2], population growth ^[5], and the growth of the total stock of private car-petrol ^[7]. On the other hand, modeling data with a declining trend, such as mortality rates, unemployment rates, and infectious diseases, remains of great interest. However, there are few SDPs available for such data. For example, in ^[3] proposed the Lomax SDP to model the adolescent fertility rate in Morocco, and in ^[4] introduced the Weibull SDP to model the age dependency ratio in Morocco.

It is worth mentioning that homogeneous SDPs assume a constant rate for the phenomenon under consideration, whereas non-homogeneous SDPs allow for variability in the rate. The latter is more reasonable in many situations. For example, interest rates are typically a function of time, making it impractical to assume a constant rate. Similarly, while mortality rates can sometimes be constant over a short period, they often change over time due to various factors, notably education and health care improvements. Therefore, modeling such phenomena requires non-homogeneous SDPs to provide insightful analysis and more accurate trend forecasting.

The SDPs introduced in the above-cited references are solved using stochastic calculus methods. Moreover, statistical inference techniques, particularly the maximum likelihood (ML) estimates, are employed to estimate the parameters involved in these models. Consequently, the trend function of the process which is a function in these parameters, can be easily estimated immediately due to the invariance property. While the ML estimates are satisfactory, the derivation of these estimates requires the functional form of the processes, which can be summarized through the transition probability density function. Unfortunately, the transition probability density function is quite difficult to obtain in many processes, which makes obtaining the estimates of these parameters impossible, and approximation methods are required in these cases, for example, ^[10,11,12], among others.

In this article, we introduce a novel stochastic model related to a particular type of the generalized Gamma-like stochastic diffusion model. The generalized Gamma probability distribution, proposed in ^[13], is particularly useful for modeling diverse types of data, especially lifetime data. The probability density function of the generalized Gamma distribution is

where $a, \; \alpha, \; p > 0.$ Note that $a$ is a scale parameter, while the other two parameters are shape parameters. Additionally, this distribution includes several other common distributions, such as the exponential, gamma, and Weibull distributions, as special cases. However, we shall consider the following version of the generalized Gamma distribution with one parameter:

where $b_{1}(\alpha), \; b_{2}(\alpha) > 0.$

We are motivated to introduce a novel SDP that is designed to model various types of trend data, particularly decreasing and increasing trends. This contrasts with many existing SDPs in the literature, which typically model only a single type of trend. Interestingly, the proposed SDP, based primarily on a generalized gamma-like distribution with a single parameter, effectively accommodates both decreasing and increasing trend data. It is important to note that adding more parameters does not necessarily improve the modeling of such data. For instance, the two-parameter Weibull stochastic diffusion, as discussed in ^[4], does not significantly enhance the analysis of the data described in the data section. Moreover, introducing additional parameters can complicate the process of obtaining maximum likelihood estimates.

Our main contributions in this paper are as follows: Firstly, we define a novel SDP capable of modeling decreasing and increasing trend data through careful selection of its drift function, which involves a single parameter. Secondly, we thoroughly study the main characteristics of the proposed process by demonstrating its existence and uniqueness, determining the transition probability density function, calculating the moments of the process, and analyzing both the trend and conditional trends. Notably, the trend function is proportional to the probability density function given in (1.1). Thirdly, the SDP parameters are estimated using the maximum likelihood method. While the likelihood function appears intractable, which is common in various SDPs, it becomes manageable due to the kernel function being a log-normal distribution. By applying the log-likelihood function, the process simplifies after some elementary algebra. Fourthly, we conduct simulation experiments to demonstrate the consistency of the maximum likelihood estimates. Finally, two real-world applications exhibiting decreasing and increasing trends are analyzed, and the SDP outperforms several existing processes.

The paper is organized as follows: Section 2 provides a description of the novel SDP and its characteristics, including the solution of the process, the transition probability density function (TPDF), and the moments. In Section 3, the process parameters are obtained using the maximum likelihood estimation method with discrete sampling of the process. A small Monte Carlo experiment is conducted in Section 4. The performance of the proposed process is applied to a real-world data set in Section 5. Finally, some concluding remarks are given in Section 6.

2.   The model and its probabilistic characteristics

In this section, we introduce a new one-dimension stochastic diffusion process. Some features of the process, such as the existence of the solution, transition probability distributions, and moments, are explained and derived.

2.1.   The generalized Gamma-like curve (GGC) stochastic diffusion model

The GGC process is defined as the non-homogeneous diffusion process depending on $\{X(t), t \in [t_1, T], t_1 > 0\}$ assuming values in $(0, \infty)$ with infinitesimal moments given by

where $t_{1}$ is the initial time and $T$ is the last time. The above-described process can be formally viewed as a solution to the following SDP:

where $\sigma > 0$, and $\alpha\in \mathbb{R}\setminus\{0\}$ is a non-zero real constant. Clearly, the solution of (2.2) can be obtained via the one-dimensional Itô's integral as follows:

2.2.   Existence and uniqueness

Here, we show the existence and uniqueness of the solution for the GGC process given via in (2.2). Toward this goal, it is enough to verify that the infinitesimal moments satisfy uniform Lipschitz and linear growth conditions; see ^[11].

Theorem 2.1. The SDE in (2.2) possesses a unique solution.

Proof. First, we show that the GGC process satisfies a uniform Lipschitz. To do so, consider $x, y\in\mathbb{R}^+$ and $t\in[t_1, T].$ It then follows that

On the other hand, the process satisfies linear growth since for $y = 0$, we have that

Thus, there exists an almost surely (a.s.) continuous process $\{x(t), t\in[t_1, T]; t_1 > 0\}$ that is the unique solution of the SDE (2.2) with probability 1. □

2.3.   Probability distribution of GGC process

The determination of the probability distribution for the solution of the GGC process plays a central role in studying the fundamental characteristics of the proposed process, especially its mean function, which serves as a basis for trend analysis. Moreover, since the process parameters are unknown, estimating these quantities using methods like maximum likelihood requires knowledge of the probability distribution of the sample path. The explicit solution of the SDE (2.2) can be obtained by considering the transformation $Y(t) = \log(X(t)).$ Upon applying Itô's formula to $y,$ we have the following:

where $Y(t_1) = \log(X_1)$. On integrating the above equation, we obtain,

and hence

Therefore, the solution in terms of the original GGC process.

Observe that the process $Y(t)$ is a Gaussian process if and only if the initial condition $Y_{1}$ is a constant or is normally distributed. Since the initial condition is constant a.s, i.e., $P(Y_{1} = y_{1}) = 1,$ and $Y(t)$ is a Markovian process, it then follows that the finite-dimensional distribution of $Y(t)$ is normal. Consequently, the finite dimensional distribution of $X(t)$ is log-normal distribution. Additionally, the TPDF of $X(t)$ given $X(s)$ where $s < t$ follows a log-normal distribution denoted by $\Lambda_1(\mu(s, t, x_{ts}), \sigma^2 (t-s))$, where $\mu(s, t, x_s)$ is given by

Therefore, the TPDF of the process considered has the following form:

2.4.   Mean function of the process

Since $X(t)$ is distributed according to $\Lambda_1(\mu(s, t, x_{s}), \sigma^2 (t-s)),$ it then follows from the properties of the lognormal distribution that the $n^{th}$ conditional moment of $X(t)$ given $X(s)$ is

The conditional mean which considers as the trend function can be obtained using $(n = 1).$ That is

On the other hand, the mean function or the unconditional trend of the process is given by

where the above equation is obtained under the assumption that $P(X(t_1) = x_1) = 1.$ Notice that when $\alpha > 0,$ it then follows that the trend function is proportional to one-parameter generalized Gamma density. Similarly, other statistical measures such as the variance, Skewness, and Kurtois can be obtained. In particular, the variance of the process is

3.   Inference on the process

Once the process, along with its basic properties, is introduced and discussed, it becomes important to examine its significance in simulation and application. However, the presence of unknown parameters makes trend analysis in practice difficult. Consequently, these parameters need to be estimated, and the maximum likelihood estimate emerges as a suitable choice, given the known functional form of the TPDF. Furthermore, ML estimates possess advantageous properties such as invariance, efficiency, and asymptotic normality.

3.1.   ML estimates

Here, the two parameters involved in the drift and diffusion functions are estimated using the ML method. We consider discrete sampling observations of the process $x(t_1), x(t_2), \dots, x(t_n)$ at times $t_1, t_2, \dots, t_n = T$. For simplicity, put $t_{j+1}-t_{j} = h$ and use $x_{i}$ to refer to $x(t_i) = x_i$. The likelihood function of $\boldsymbol\theta = (\alpha, \sigma^{2})^{T}$ can be obtained from Eq (2.4), taking into account that the initial condition taking $P(X(t_1) = x_1) = 1,$ is

The log-likelihood equation is

The $\boldsymbol\ell(\alpha, \sigma^{2})$ can be maximized by solving the nonlinear likelihood equation obtained by differentiating with respect to $\boldsymbol\theta = (\alpha, \sigma^{2})^{T}$. The first partial derivatives of $\boldsymbol\ell(\alpha, \sigma^{2})$ are given by

where

Let $S(\boldsymbol\theta) = (\partial\boldsymbol\ell(\boldsymbol\theta)/\partial\alpha, \partial\boldsymbol\ell(\boldsymbol\theta)/\partial\sigma^{2})^{T}$ be the score function. Then, the MLE of $\widehat{\boldsymbol\theta} = (\widehat{\alpha}, \widehat{\sigma}^{2})$ can be obtained by solving $S(\boldsymbol\theta) = 0$. Unfortunately, the MLEs cannot be obtained in closed form, and numerical methods are required to obtain these estimates. From Eq (3.3), we obtain

On substituting (3.4) in Eq (3.2) the following nonlinear equation is obtained for the estimator $\hat{\alpha}$:

Let $g(\alpha)$ be the left-hand side equation (3.6). Therefore, the ML estimate of $\alpha$ can be achieved by solving the non-linear equation $g(\alpha) = 0.$

3.2.   Estimated trend functions and confidence bounds

Once the ML estimates of $\alpha$ and $\sigma^{2}$ are obtained, we proceed to provide estimates for the conditional mean and the mean of the process. Due to the invariance property of the ML estimates, see for example Theorem [5.28,308] in ^[14], the ML estimates for the conditional mean and the unconditional mean can be obtained. Let $\widehat{\alpha}$ and $\widehat{\sigma}^{2}$ be the ML estimates of $\alpha$ and $\sigma^{2},$ respectively, then the estimated conditional mean of the process (ECMF) is

Similarly, the estimated mean function (EMF) of the process is

where in the above equation we used the assumption that $P(X(t_{1}) = x(t_{1})) = 1.$

At time, we would like to obtain a confidence band for the CMF and MF of the process. From Eq (2.4), we have that for $t > s$, $X(t)|X(s)$ follows $\Lambda_1(\mu(s, t, x_{xs}), \sigma^2 (t-s).)$ Therefore, we have that

Consequently, a $(1-\alpha)100\%$ confidence band for $z$ is determined by $P(-Z_{\gamma/2}\leq Z\leq Z_{\gamma/2}) = 1-\gamma,$ for $\gamma\in(0, 1).$

Therefore, the $(1-\gamma)100\%$ confidence bound (CB) for $X(t)$ is

where,

4.   Simulation experment

Here, our primary goal is to illustrate the pattern of the GGC process by simulating sample trajectories. Additionally, we use these sample paths to evaluate the frequentist performance of the ML estimates for the trend function of the GGC process. The following algorithm describes the procedure for simulating a trajectory of the GGC process computed at $N$ time points.

Algorithm

1: Initialize, $t_{1},$ $T,$ $N,$ $x_{1}$, $\alpha,$ $\sigma$,

2: Compute $\Delta: = \frac{T-t_{1}}{N},$

3: For $j = 2, \dots, N,$ do

(a) Discretize the time interval $[t_{1}, T],$ into $N$ time points with each point computed as $t_{i} = t_{i-1}+(i-1)\Delta,$

(b) Generate a standard normal observation, i.e., $Z:N(0, 1),$ and the increments are computed as $W_{j} = W_{j-1}+\sqrt{\Delta}\; Z_{i},$

(c) The sample path of the GGC process can be obtained upon using (2.3), i.e., $X_{j} = X_{j-1}\left(\frac{t_{j}}{t_{j-1}}\right)^{ \alpha }\exp\left\{-\frac{10^3}{-10^2+\alpha} (t_{j}^{\frac{-10^2}{\alpha}+1}-t_{j-1}^{\frac{-10^2}{\alpha}+1})- \frac{\sigma ^2 }{2}(t_{j}-t_{j-1})\right\}\exp\left(\sigma (W_{j}-W_{j-1})\right).$

4: End do

5: Repeat step 3 $m$ times to obtain $m$ sample paths.

6: For $i = 1, \dots, m,$ use the sample path $X_{i}(t_{j})$ in step 5 to compute the ML estimates, using Eqs (3.2) and (3.3), $\widehat{\alpha}_{i}$ and $\widehat{\sigma}_{i},$ respectively.

7: The ML estimates are then given by: $\widehat{\alpha} = m^{-1}\sum_{i = 1}^{m}\widehat{\alpha}_{i}$ and $\widehat{\sigma} = m^{-1}\sum_{i = 1}^{m}\widehat{\sigma}_{i}$ with mean square errors respectively, $MSE(\widehat{\alpha}) = (m-1)^{-1}\sum_{i = 1}^{m}(\widehat{\alpha}_{i}-\alpha)^{2}$ and $MSE(\widehat{\sigma^{2}}) = (m-1)^{-1}\sum_{i = 1}^{m}(\widehat{\sigma^{2}}_{i}-\sigma^{2})^{2}$.

We generate $m = 50$ training sample paths using the algorithm described above, varying the parameters in the graphs. All computations in this paper were performed using R software ^[15].

(1) In the first scenario, we set the time interval $[t_1, T] = [1000, 1010],$ $N = 100,$ and $X(t_{1}) = 9000.$ Figure 1 shows that the trajectories vary in direction based on the sign of the parameter $\alpha$. The dispersion of the process is influenced by the value of $\sigma$; smaller values result in trajectories that closely follow a single curve, while larger $\sigma$ values lead to more dispersed yet similarly trending curves. Using the 50 sample paths, we obtain the maximum likelihood (ML) estimates of the process parameters and subsequently the ML estimates of the trend function. Figure 1 demonstrates that the estimated trend closely matches the actual trend function, indicating the reliability of our estimates.

(2) In the second scenario, we consider $[t_1, T] = [500, 510],$ $N = 100,$ and $X(t_{1}) = 3000.$ Similarly, Figure 2 reveals that the process exhibits either a decreasing or increasing trend. The sample paths tend to either closely follow a single curve or scatter depending on whether the dispersion parameter $\sigma$ is small or large, respectively. Additionally, using 50 sample trajectories, we obtain the ML estimates of $\alpha$ and $\sigma$, and consequently the estimate of the trend function. A closer examination shows that the ML estimates behave consistently across different scenarios.

5.   Applications

This section presents the application of our model to two real-world datasets: the number of infant deaths in the United Kingdom from 1977 to 2020 and $CO_2$ emissions (in kilotons) in Morocco from 1990 to 2020. These annualized data are sourced from the World Bank database. The first dataset is fitted using the GGC process and compared with two existing models: the two-parameter Weibull SDP ^[4] and the Lomax SDP ^[3], both of which are suitable for modeling data with decreasing trends. Similarly, the second dataset is fitted using the GGC process and compared with the log-logistic SDP ^[1], which is appropriate for modeling data with increasing trends.

5.1.   Application to number of infant deaths in the United Kingdom

According to the Royal College of Paediatrics and Child Health (RCPCH), infant death rates in all countries of the UK have significantly fallen over the past 40 years. Most childhood deaths occur during the first year of life, particularly in the first month (neonatal period). Newborn deaths account for 70% to 80% of infant deaths. The vast majority of newborn deaths are due to perinatal causes, especially preterm birth, and are strongly linked to maternal health and congenital anomalies. The remaining infant deaths occur due to a wide variety of causes, including sudden unexplained death in infancy (SUDI). The overall decline in infant death rates since 1980 likely reflects general improvements in health care, specifically prenatal and newborn care. Breastfeeding and safe sleeping positions are protective factors for infant survival, especially for premature babies. This study aims to use a stochastic diffusion process to model the number of infant deaths and to forecast the number of infant deaths in the United Kingdom. We use the data from 1977 to 2018, available at https://databank.worldbank.org/source/world-development-indicators, as training data to estimate the model parameters: $\hat{\alpha} = -1779.057$ and $\hat{\sigma} = 0.02208178.$ Consequently, the estimated trend functions (unconditional EMF and conditional ECMF) are obtained immediately due to the invariance of ML estimates. Additionally, we use the years 2019 and 2020 to forecast the number of infant deaths in the UK. The results for the forecasted values for these two years are presented in Table 1 and shown in Figure 3. These results demonstrate that the GGC process is quite effective in predicting the values for these years, particularly when employing ECMF. Figure 3 (left panel) displays the observed data, estimated trend function (EMF), and estimated confidence bands, revealing that the GGC process fits the current data well and provides accurate predictions. Figure 3 (right panel) further illustrates the conditional trend function estimation along with forecasts. Moreover, we compare the proposed diffusion process with the Weibull ^[4] and Lomax ^[3] diffusion processes using the current data. Table 2 lists the estimated parameters and the Akaike information criterion (AIC), showing that the GGC process has the lowest AIC value and therefore outperforms the two-parameter diffusion processes. This is also well demonstrated in Figure 3. Moreover, Figure 4 shows the fits made using the methods mentioned in Table 2, revealing that the GGC process performed much better than the Weibull ^[4] and Lomax ^[3] diffusion processes in modelling the current data.

The accuracy of the forecast can be quantified using measures such as the mean absolute error (MAE), the root mean square error (RMSE), and the mean absolute percentage error (MAPE), defined as follows: $\text{MAE} = \frac{1}{42}\sum_{i = 1}^{42}{ | x_(t_{i})-\widehat{x}(t_{i}) |} = 257.9876,$ $\text{RMSE} = \sqrt{\frac{1}{42}\sum_{i = 1}^{42}{ | x(t_{i})-\widehat{x}(t_{i}) |^2}} = 330.2669,$ and $\text{MAPE} = \frac{1}{42}\sum_{i = }^{42}{\frac{ | x(t_{i})-\widehat{x}(t_{i}) |}{ x(t_{i})}}\times 100 = 5.1237.$ Table 3 shows the observed data along with their corresponding forecast. The value obtained for MAPE is less than 10, and this indicates that we obtained a high-accuracy prediction.

5.2.   Application to $CO_2$ emissions (kt) in Morocco

Another application of our model is to analyze Carbon Dioxide ($CO_2$) emissions in Morocco, measured in kilotons (kt). $CO_2$ emissions primarily result from burning fuels and cement manufacturing. Since the Industrial Revolution, the increasing combustion of carbon-based fuels has significantly raised atmospheric $CO_2$ concentrations, accelerating global warming and contributing to man-made climate change. Additionally, $CO_2$ dissolves in water to form carbonic acid, which leads to ocean acidification. This subsection examines $CO_2$ emissions (kt) in Morocco from 1990 to 2020, with data sourced from https://databank.worldbank.org/source/world-development-indicators. We use the data from 1990-2018 as training data to estimate the model parameters: $\hat{\alpha} = 81.55085457$ and $\hat{\sigma} = 0.02977168$. Therefore, the estimated trend functions both (the unconditional EMF and the conditional ECMF) are derived directly due to the invariance of the ML estimates. Similarly, we use the data from 2019 and 2020 to project $CO_2$ emissions in Morocco. The predicted values for these two years are reported in Table 4 and illustrated in Figure 5. Once again, the results demonstrate that the GGC process is capable of accurately predicting the values for these years using both EMF and ECMF. Figure 5 (left panel) depicts the observed data, the estimated trend function (EMF), and the estimated confidence intervals, indicating that the GGC process accurately fits the data and yields precise forecasts. Figure 5 (right panel) further demonstrates the conditional trend function estimation and projections. Furthermore, we compare the proposed diffusion process with the log-logistic diffusion process ^[1] using the current dataset. Table 5 presents the estimated parameters and the AIC criterion, showing that the GGC process has a lower AIC value and thus surpasses the two-parameter diffusion process. This is also clearly illustrated in Figure 5. Furthermore, Figure 6 displays the fits performed using the methods outlined in Table 5, demonstrating that the GGC process outperformed the Log-Logistic diffusion process in modeling this data.

Similarly, the accuracy of the forecast can be assessed using the metrics: $\text{MAE} = 2567.19,$ $\text{RMSE} = 3115.123,$ and $\text{MAPE} = 6.032141.$ Table 6 shows the observed data along with their corresponding forecast. The MAPE value obtained is less than 10, indicating a high-accuracy prediction.

6.   Conclusions

We introduce a new stochastic diffusion process based on generalized Gamma-like curves, referred to as the GGC process. The GGC process is capable of modeling both increasing and decreasing trend data. We investigate and derive several structural properties, including the explicit solution of the process, the transition probability density function, and the moments of the process. The maximum likelihood method is employed to estimate the GGC process parameters. Simulation studies demonstrate that the ML estimates are consistent even in small samples. The potential of the GGC process for analyzing data with a declining trend is illustrated by fitting it to infant mortality data in the UK. The proposed diffusion process outperforms two existing popular stochastic diffusion processes that are designed to fit decreasing trend data. Additionally, it is applied to $CO_2$ emissions in Morocco, which exhibit an increasing trend. The GGC process again showed satisfactory performance in modeling this data and outperformed the log-logistic SDP, which is suitable for modeling growth data.

As with any stochastic diffusion process or model, the GGC process may be inappropriate for modeling unimodal or bathtub-shaped data. This can be addressed by selecting different infinitesimal moments or by incorporating additional suitable parameters.

For future work, we plan to model the first passage time using the GGC process with more applications. Additionally, we will explore other estimation methods, particularly the Bayesian approach, which may be more powerful in handling other application scenarios.

Author contributions

Conceptualization, M.K.S.; methodology, M.K.S.; writing-original draft preparation, M.K.S.; software, M.K.S and S.A.; validation, S.A.; formal analysis, M.K.S. and S.A. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

We would like to thank the Editor, Associate Editor, and the three anonymous referees for their insightful comments and valuable suggestions, which greatly improved our paper.

Conflict of interest

No potential conflict of interest was reported by the authors.