Numerical Methods in Python for Solving Systems of Stochastic Differential Equations: An Application to the SIRD Model in Latin America

Abstract

In this chapter, we apply the exploration of the Euler–Maruyama, Milstein, and Runge–Kutta methods to solve systems of stochastic differential equations associated with the stochastic SIRD model. We simulate sample trajectories of each variable using Python and data collected in Peru during the years 2020 and 2021, marked by the onset of the COVID-19 pandemic. Our research involves comparing stochastic and deterministic systems obtained from the SIRD model in the Peruvian context. We use different values for the intensity of white noise to assess the impact of stochasticity on the dynamics of the SIRD model. Presenting random simulations alongside deterministic ones provides a comprehensive understanding of the effect of randomness in the context of infectious disease modeling.

1. Introduction

In our investigation into the dynamics of infectious diseases, we leverage advanced numerical methods to address the complexities posed by systems governed by stochastic differential equations (SDE). Our primary focus is on the widely recognized SIRD (Susceptible-Infectious-Recovered-Deceased) model, particularly relevant in the Latin American context.

Python serves as our computational tool of choice, ensuring the efficient implementation of numerical methods. Our study is grounded in data collected in Peru during 2020 and 2021, a period marked by the profound impact of the COVID-19 pandemic. This temporal context provides a unique opportunity for a meticulous examination of the influence of stochasticity on the SIRD model, as evidenced by our previous works [1, 2].

Our methodology involves simulating trajectories of SIRD variables in Python, considering varying intensities of white noise resembling Brownian motion. We employ the Euler–Maruyama, Milstein, and Runge–Kutta methods, each strategically chosen for its precision and relevance to infectious disease dynamics within the framework of stochastic epidemiological systems.

In the mathematical modeling of COVID-19 transmission, our study aligns with prior research [1], which utilized the deterministic SIR model for COVID-19 transmission in Peru. In parallel, our investigation [2] employs the stochastic SIRD model, utilizing the Milstein method to simulate and demonstrate results with Peruvian data. This approach promises a nuanced exploration of SIRD dynamics under diverse conditions, dissecting the interplay between deterministic trends and stochastic fluctuations. Our work contributes valuable insights to the field of infectious disease dynamics, shedding light on the complexities inherent in real-world scenarios.

2. SIRD model

The SIRD model, representing Susceptible-Infectious-Recovered-Deceased, offers a theoretical foundation for capturing the dynamics of infectious diseases in a population. This model categorizes individuals into four compartments based on their health status concerning the specific disease.

1. Susceptible (S): Individuals in this group are susceptible to the infectious agent but have not yet been exposed or infected. The number of susceptibles diminishes as individuals encounter the disease.

2. Infectious (I): This compartment includes individuals currently infected and capable of transmitting the disease to susceptible individuals. The duration of infectivity significantly influences disease spread.

3. Recovered (R): Individuals in this category have recovered from the infection and are assumed to have acquired immunity. Recovered individuals are no longer susceptible or infectious.

4. Deceased (D): Unfortunately, some individuals may succumb to the disease, transitioning to the deceased category. The mortality rate is a crucial parameter governing this transition.

The SIRD model’s dynamics are commonly expressed through a system of differential equations, where transition rates between compartments depend on parameters such as transmission, recovery, and mortality rates. This model provides valuable insights into the evolution of infectious diseases within a population, aiding researchers and policymakers in assessing intervention impact and disease control strategies.

In epidemiology, the SIRD model stands as a fundamental tool, offering a quantitative understanding of interactions among population compartments during infectious disease outbreaks. Its widespread use is attributed to its versatility and simplicity, making it a valuable framework for analyzing and predicting the spread of infectious diseases.

2.1 Deterministic model

In a given population N, if there are individuals affected by a virus, the population is categorized into the infected group (I), the susceptible group (S), the recovered group (R), and the deceased group (D). These groups undergo changes over time following the depicted scheme in Figure 1.

Figure 1.

The progression of a pandemic’s development over time.

Expressed as a modification of the mathematical model proposed in 1927 by Kermack and McKendrick [3], it is represented by the system of four equations:

$$\left\{\begin{array}{ll}{S}^{\prime }\left(t\right)& =-\mathit{\lambda S}\left(t\right)I\left(t\right),\\ I^{\prime }\left(t\right)& =\mathit{\lambda S}\left(t\right)I\left(t\right)-\mathit{\gamma I}\left(t\right)-\mathit{\mu I}\left(t\right),\\ R^{\prime }\left(t\right)& =\mathit{\gamma I}\left(t\right),\\ D^{\prime }\left(t\right)& =\mathit{\mu I}\left(t\right).\end{array}\right.$$E1

Referred to as the SIRD model, the variables and parameters are detailed in Table 1.

Variables and parameters	Meaning
S(t)	Susceptible population at time t,
I(t)	Infected population at time t,
R(t)	Recovered population at time t,
D(t)	Dead population at time t,
t	The time expressed in days,
$\lambda$	Transmission rate,
$\gamma$	Rate of recovery,
$\mu$	Mortality rate.

Table 1.

Variables and parameters associated with the system (1).

All variables are required to meet the compatibility condition $N=S\left(t\right)+I\left(t\right)+R\left(t\right)+D\left(t\right)$, ensuring the constancy of the total population.

2.2 Stochastic model

Given the significant impact of environmental factors on COVID-19 spread, coupled with the inherent uncertainty in data accuracy, the stochastic aspect becomes particularly intriguing. The complexity is heightened by scenarios where infected individuals, especially children, may not show symptoms but contribute to transmission. Additionally, asymptomatic individuals, unmonitored yet active in community transmission, add another layer of complexity.

Transmission through contact with contaminated surfaces further widens the scope of possibilities. Acknowledging literature on modeling environmental fluctuations’ influence on population dynamics [4, 5], two approaches emerge. One involves introducing random perturbations for each variable, delving into the less-understood connections among susceptible, infected, and recovered individuals through white noise. Following [6], a stochastic system is proposed as follows:

$$\left\{\begin{array}{ll}\mathit{dS}\left(t\right)& =-\mathit{\lambda S}\left(t\right)I\left(t\right)\mathit{dt}+{\sigma }_{1}S\left(t\right){\mathit{dW}}_1\left(t\right),S\left(0\right)=s_0,\\ \mathit{dI}\left(t\right)& =\left(\mathit{\lambda S}\left(t\right)I\left(t\right)\mathit{dt}-\mathit{\gamma I}\left(t\right)-\mathit{\mu I}\left(t\right)\right)\mathit{dt}+{\sigma }_{2}I\left(t\right){\mathit{dW}}_2\left(t\right),I\left(0\right)=i_0,\\ \mathit{dR}\left(t\right)& =\mathit{\gamma I}\left(t\right)\mathit{dt}+{\sigma }_{3}R\left(t\right){\mathit{dW}}_3\left(t\right),R\left(0\right)=r_0,\\ \mathit{dD}\left(t\right)& =\mathit{\mu I}\left(t\right)\mathit{dt}+{\sigma }_{4}D\left(t\right){\mathit{dW}}_4\left(t\right),D\left(0\right)=d_0.\end{array}\right.$$E2

Here, $s_0,i_0,r_0,d_0\ge 0$; $W_i\left(t\right)$ are independent standard Brownian motions, and ${\sigma }_i^2\ge 0$ represent the intensities of for i = 1, 2, 3, 4.

The entire system is established on a comprehensive probability space $\left(\mathrm{\Omega },\mathrm{ℱ},P\right)$ accompanied by a filtration ${\mathrm{ℱ}}_t\ge 0$ that adheres to standard conditions (i.e., it is right-continuous and increasing, with ${\mathrm{ℱ}}_0$ encompassing all null sets on P).

3. Numerical methods

Stochastic differential equations (SDEs) serve the purpose of describing phenomena influenced by both deterministic and random terms.

$$\left\{\begin{array}{l}\mathit{dX}\left(t\right)=f\left(X\left(t\right)\right)\mathit{dt}+g\left(X\left(t\right)\right)\mathit{dW}\left(t\right),t\ge 0,\\ X\left(0\right)=X_0.\end{array}\right.$$E3

Eq. (3) consists of a right-hand side depending on two components, the function $f:{\mathrm{ℝ}}^d\to {\mathrm{ℝ}}^d$, recognized as the drift, representing the deterministic aspect and the function $g:{\mathrm{ℝ}}^d\to {\mathrm{ℝ}}^{d\times m}$, referred to as the diffusion, governing the stochastic part.

The term W(t) in (3) denotes an m-dimensional standard Wiener process. Because of its lack of differentiability everywhere (with probability 1), the expression in Eq. (3) serves as a concise notation for its integral form:

$$X\left(t\right)=X\left(0\right)+{\int }_0^{t}f\left(X\left(s\right)\right)\mathit{ds}+{\int }_0^{t}g\left(X\left(s\right)\right)\mathit{dW}\left(s\right).$$E4

If the stochastic integral in (4) is an Itô integral, the equation is an Itô stochastic differential equation. If it is a Stratonovich integral, (4) represents a Stratonovich stochastic differential equation. Unless stated otherwise, our discussion primarily revolves around Itô SDEs. Notably, Itô SDEs can be equivalently formulated in Stratonovich form:

$$\mathit{dX}\left(t\right)=\left(f\left(t\right)-\frac{1}{2}{g}^{\prime }\left(t\right)\right)\mathit{dt}+g\left(t\right)\circ \mathit{dW}\left(t\right).$$E5

This interchangeability allows transitioning between Itô and Stratonovich SDEs, providing equivalent formulations in terms of the solution [7, 8].

3.1 Euler: Maruyama

The Euler–Maruyama (EM) method is recognized as one of the simplest computational techniques for approximating Stochastic Differential Equations (SDEs), serving as the stochastic counterpart to the Euler method employed in ordinary differential equations. Before delving into the stochastic scenario, let us provide a brief introduction to this method.

To approximate the solution of a Stochastic Differential Equation (SDE) expressed in the form (3) over the interval [0, T], the interval is initially discretized into L parts $\left(L\in \mathrm{\mathbb{N}}\right)$. Here, $\mathrm{\Delta }t=T/L$ represents the time step size, and ${\tau }_j=j\mathrm{\Delta }t,$ j = 0, …, L. Our numerical approximation is denoted as $X_j$. The Euler–Maruyama method, when applied to an SDE of the type (3), is precisely defined by the following recursive formula:

$$X_{j+1}=X_j+f\left(X_j\right)\mathrm{\Delta }t+g\left(X_j\right)\mathrm{\Delta }{W}_{j+1},j=0,1,\dots ,L-1.$$E6

In this context, we compute the discretized trajectories of the Brownian motion and use them to generate the increments $\mathrm{\Delta }{W}_{j+1}=W\left({\tau }_{j+1}\right)-W\left({\tau }_j\right)$ needed in the above recursive formula. For convenience, we always choose the time step $\mathrm{\Delta }t={\tau }_{j+1}-{\tau }_j$ as an integer multiple $R\ge 1$ of the increment $\delta =t_{j+1}-t_j$ for the trajectory that generates the Brownian sample, that is, $\mathrm{\Delta }t=\mathit{R\delta t}$. This ensures that the set of points $\left\{{\tau }_j\right\}$ in each discretized Brownian trajectory contains the points $\left\{{\tau }_j\right\}$ in each calculated solution of the Euler–Maruyama method. The increments $W\left({\tau }_{j+1}\right)-W\left({\tau }_j\right)$ are given by:

$$W\left({\tau }_{j+1}\right)-W\left({\tau }_j\right)=W\left(\left(j+1\right)\mathit{R\delta t}\right)-W\left(\mathit{\text{jRδt}}\right)=\sum _{k=\mathit{jR}+1}^{\left(j+1\right)R}\delta W_k.$$E7

Finally, we describe the Euler–Maruyama method as follows:

Euler–Maruyama Method (3.1):

$$\left\{\begin{array}{ll}{X}_0& =x_0,\mathrm{\Delta }t=\mathit{R\delta t}=R\left(t_{j+1}-t_j\right)={\tau }_{j+1}-{\tau }_j,\\ X_{j+1}& =X_j+f\left(X_j\right)\mathrm{\Delta }t+g\left(X_j\right)\left\{W\left({\tau }_{j+1}\right)-W\left({\tau }_j\right)\right\},j=0,1,\dots ,L-1.\end{array}\right.$$E8

The method mentioned above is of the explicit one-step type since each value $X_{j+1}$ is expressed in terms of a single previous value, $X_j$. The denomination Euler–Maruyama method is then the summa of the work of Euler (1707–1783) and Gisiro Maruyama (1916–1986). Subsequently, the Euler–Maruyama method emerges as a first-order truncation of the Itô-Taylor expansion of the precise solution to (3).

3.2 Milstein method

Additional formulae for approximating SDEs (3) may arise by including further terms from the Itô-Taylor expansion in the numerical method. An instance of this is the Milstein method, which, for scalar problems, can be obtained as follows. Applying the Itô formula to the expressions $f\left(X_s\right)$ and $g\left(X_s\right)$, where these represent the coefficients of the Stochastic Differential Equation (SDE), we can derive

$$\begin{array}{ll}{X}_{t_{j+1}}=& X_{t_j}+{\int }_{t_j}^{t_{j+1}}(f\left(X_{t_j}\right)+{\int }_{t_j}^s\left(f^{\prime }\left(X_u\right)f\left(X_u\right)+\frac{1}{2}{f^{\prime }}^{\prime }\left(X_u\right)g{\left(X_u\right)}^2\right)\mathit{du}\\ & +{\int }_{t_j}^s{f}^{\prime }\left(X_u\right)g\left(X_u\right){\mathit{dW}}_u)\mathit{ds}+{\int }_{t_j}^{t_{j+1}}(g\left(X_{t_j}\right)+{\int }_{t_j}^s\left(g^{\prime }\left(X_u\right)f\left(X_u\right)+\frac{1}{2}{g^{\prime }}^{\prime }\left(X_u\right)g{\left(X_u\right)}^2\right)\mathit{du}\\ & +{\int }_{t_j}^s{g}^{\prime }\left(X_u\right)g\left(X_u\right){\mathit{dW}}_u){\mathit{dW}}_s.\end{array}$$E9

We can ignore the double integrals of the form ${\mathit{dW}}_s\mathit{du}$ and ds du. Then, we obtain the following expression:

$$X_{t_{j+1}}\approx X_{t_j}+{\int }_{t_j}^{t_{j+1}}f\left(X_{t_j}\right)\mathit{ds}+{\int }_{t_j}^{t_{j+1}}\left(g\left(X_{t_j}\right)+{\int }_{t_j}^s{g}^{\prime }\left(X_u\right)g\left(X_u\right){\mathit{dW}}_u\right){\mathit{dW}}_s.$$E10

This can be expressed as:

$$X_{t_{j+1}}\approx X_{t_j}+f\left(X_{t_j}\right)\mathrm{\Delta }t+g\left(X_{t_j}\right)\mathrm{\Delta }{W}_{j+1}+{\int }_{t_j}^{t_{j+1}}{\int }_{t_j}^s{g}^{\prime }\left(X_u\right)g\left(X_u\right){\mathit{dW}}_u{\mathit{dW}}_s.$$

The first three terms are familiar from the Euler–Maruyama method. The fourth term can be approximated as:

$${\int }_{t_j}^{t_{j+1}}{\int }_{t_j}^s{g}^{\prime }\left(X_u\right)g\left(X_u\right){\mathit{dW}}_u{\mathit{dW}}_s\approx g^{\prime }\left(X_{t_j}\right)g\left(X_{t_j}\right){\int }_{t_j}^{t_{j+1}}{\int }_{t_j}^s{\mathit{dW}}_u{\mathit{dW}}_s.$$

The integral on the right-hand side is known to be:

$${\int }_{t_j}^{t_{j+1}}{\int }_{t_j}^s{\mathit{dW}}_u{\mathit{dW}}_s=\frac{1}{2}\left({\left(\mathrm{\Delta }{W}_{j+1}\right)}^2-\mathrm{\Delta }t\right)$$

Substituting into our previous approximation, we finally obtain the Milstein method Milstein Method (3.2):

$$\left\{\begin{array}{l}{X}_0=x_0\\ X_{j+1}=X_j+f\left(X_j\right)\mathrm{\Delta }t+g\left(X_j\right)\mathrm{\Delta }{W}_{j+1}+\frac{1}{2}{g}^{\prime }\left(X_j\right)g\left(X_j\right)\left({\left(\mathrm{\Delta }{W}_{j+1}\right)}^2-\mathrm{\Delta }t\right)\end{array}\right.$$E11

while for systems of SDEs,

$$\mathit{dX}\left(t\right)=g^0\left(X\left(t\right)\right)\mathrm{\Delta }t+\sum _{k=1}^m{g}^k\left(X\left(t\right)\right){\mathit{dW}}^k\left(t\right),$$

where $g^k:{\mathrm{ℝ}}^d\to {\mathrm{ℝ}}^d,$, k = 0, 1, …, m, the method is given by

$$\begin{array}{ll}{X}_{n+1}=& X_n+g^0\left(X_n\right)\mathrm{\Delta }t+\sum _{k=1}^m{g}^k\left(X_n\right)\mathrm{\Delta }{W}_{n+1}^k\\ & +\sum _{j_1,j_2=1}^m{L}^{j_1}{g}^{j_2}\left(X_n\right){\int }_{t_n}^{t_{n+1}}{\int }_{t_n}^t{\mathit{dW}}^{j_1}\left(s\right){\mathit{dW}}^{j_2}\left(t\right),\end{array}$$

with $L^k={\sum }_{i=1}^d{g}^{k,i}\frac{\partial }{\partial x_i}$, k = 1, 2, …, m.

The Milstein method, also a one-step explicit approach, achieves increased accuracy compared to the Euler–Maruyama method due to the inclusion of additional terms originating from the truncation of the Itô-Taylor expansion.

3.3 Runge–Kutta

A category of stochastic Runge–Kutta methods (SRK methods) has been derived by appropriately perturbing deterministic Runge–Kutta methods. The formulation we adopt is based on the framework outlined in [9], which characterizes SRK methods as follows:

$$X_{n+1}=X_n+\mathrm{\Delta }t\sum _{i=1}^s{b}_{i}f\left({\widehat{X}}_i\right)+\mathrm{\Delta }{W}_{n+1}\sum _{i=1}^s{q}_{i}g\left({\widehat{X}}_i\right),$$E12

Where the intermediate values are used to estimate the value of $X\left(t_n+c_i\mathrm{\Delta }t\right)$ for every index i = 1, 2, …, s, are provided by

$${\widehat{X}}_i=X_n+\mathrm{\Delta }t\sum _{j=1}^s{a}_{\mathit{ij}}f\left({\widehat{X}}_j\right)+\mathrm{\Delta }{W}_{n+1}\sum _{j=1}^s{\gamma }_{\mathit{ij}}g\left({\widehat{X}}_j\right),$$E13

A more compact effective representation of SRK methods (12)–(13) is obtained by using a Butcher tableau:

with $\mathbf{A}=\left({\mathbf{a}}_{\mathbf{i}\mathbf{,}\mathbf{j}}\right),\mathrm{\Gamma }=\left({\gamma }_{i,j}\right)\in {\mathrm{ℝ}}^{s\times s}$, and vector of weights $\mathbf{b}=\left({\mathbf{b}}_{\mathbf{j}}\right)$, $\mathbf{q}=\left({\mathbf{q}}_{\mathbf{j}}\right)\in {\mathrm{ℝ}}^s$ and vector of abscissae $\mathbf{c}=\left({\mathbf{c}}_{\mathbf{j}}\right)\in {\mathrm{ℝ}}^s$.

If $\mathrm{\Gamma }$ is the zero matrix and q the null vector, then SRK methods (12)–(13) recovering the equivalent family of deterministic Runge–Kutta methods is also applicable. In the stochastic scenario, the computational cost of the method is influenced by the composition of matrices in (13). Specifically, if A and $\mathrm{\Gamma }$ are strictly lower triangular, the corresponding method is explicit.

3.4 Convergence analysis

Convergence in the stochastic context differs due to the randomness of the error, denoted as $\Vert X_n-X\left(t_n\right)\Vert$, subject to random fluctuations. The error is measured using the expectation operator.

Assuming analysis within the interval [0, T], partition it into L subintervals, with $\mathrm{\Delta }t=\frac{T}{L}$. The partition intersects a set of grid points ${\mathrm{ℐ}}_{\mathrm{\Delta }t}=\left\{t_j=j\mathrm{\Delta }t|j=0,1,\dots ,L\right\}$. Typically, these points are a subset of Wiener points, $\mathrm{\Delta }t=\mathit{R\delta t},R\in \mathrm{\mathbb{N}}$.

Definition 1: A numerical method for SDEs (3) providing the numerical solution ${\left\{X_n\right\}}_{n=0}$ is strongly convergent if

$$\underset{\mathrm{\Delta }t\to 0}{lim}\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mathbb{E}\left[||X_n-X\left(t_n\right)||\right]=0.$$E14

Moreover, we say that its strong order of convergence is p if there exist two positive numbers C and $\mathrm{\Delta }{t}^{\ast }$ such that

$$\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mathbb{E}\left[||X_n-X\left(t_n\right)||\right]\le C\mathrm{\Delta }{t}^p,\text{for}\mathrm{any}\mathrm{\Delta }t\le \mathrm{\Delta }{t}^{\ast }.$$E15

It is noteworthy to mention that in numerous real-world scenarios, an integer p that satisfies (15) may not be present, particularly in cases involving low regularity problems, such as certain stochastic partial differential equations. An alternative concept of convergence is presented as follows:

Definition 2: A numerical method for SDEs (3) providing the numerical solution is weakly convergent if

$$\underset{\mathrm{\Delta }t\to 0}{lim}\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[\mathrm{\Phi }\left(X_n\right)\right]-\mathbb{E}\left[\mathrm{\Phi }\left(X\left(t_n\right)\right)\right]\mid =0$$E16

for any test function $\mathrm{\Phi }$ within a designated function space S. Furthermore, we assert that its weak order of convergence is p if there exist two positive numbers C and $\mathrm{\Delta }{t}^{\ast }$, satisfying

$$\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[\mathrm{\Phi }\left(X\right)\right]-\mathbb{E}\left[\mathrm{\Phi }\left(X\left(t_n\right)\right)\right]\mid \le C\mathrm{\Delta }{t}^p,\text{for}\mathrm{any}\mathrm{\Delta }t\le \mathrm{\Delta }{t}^{\ast }.$$E17

In numerous real-world scenarios, a suitable choice for the functional space S is given by the space of algebraic polynomials limited to a specified degree. It is crucial to underscore the connection between the two previously mentioned convergence concepts. Assuming that $\mathrm{\Phi }\left(x\right)=x$, a method demonstrating strong convergence is inherently weakly convergent.

$$\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[\mathrm{\Phi }\left(X_n\right)\right]-\mathbb{E}\left[\mathrm{\Phi }\left(X\left(t_n\right)\right)\right]\mid =\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[X_n\right]-\mathbb{E}\left[X\left(t_n\right)\right]\mid$$E18

Since

$$\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[X_n\right]-\mathbb{E}\left[X\left(t_n\right)\right]\mid =\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[X_n-X\left(t_n\right)\right]\mid \le \underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mathbb{E}\left[\Vert X_n-X\left(t_n\right)\Vert \right],$$E19

We obtain

$$\underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mid \mathbb{E}\left[\mathrm{\Phi }\left(X\right)\right]-\mathbb{E}\left[\mathrm{\Phi }\left(X\left(t_n\right)\right)\right]\mid \le \underset{t_n\in {\mathrm{ℐ}}_{\mathrm{\Delta }t}}{sup}\mathbb{E}\left[\Vert X_n-X\left(t_n\right)\Vert \right],$$E20

as we aimed to prove.

Below, we present some results on the convergence of the numerical methods discussed in this section. The proofs are omitted and can be found in [9].

Theorem 1: Euler–Maruyama method (8) is strongly convergent with a strong order of $\frac{1}{2}$.

Theorem 2: Euler–Maruyama method (8) is weakly convergent with weak order 1.

It is also possible to demonstrate the following, [9]:

• Both the strong and weak orders of convergence for the Milstein method (10) are equivalent and set at 1.

• Consequently, the Milstein method represents an improvement over the Euler–Maruyama method, particularly in terms of strong convergence.

• The convergence of explicit SRK methods (12)–(13) relies on the convergence of the corresponding deterministic RK methods, subject to an additional condition. This additional condition is necessary to ensure that

$$\sum _{i=1}^s{b}_i=\sum _{i=1}^s{q}_i=1.$$

4. An application to the case of Peru

4.1 Data and equations

Based on the national census conducted on October 22, 2017, the population of Peru is reported to be 31,237,385 individuals. Subsequent yearly estimates, as of June 30, according to [10], project the population to approximately N = 32,625,948, as indicated by this institution.

To apply the system of eqs. (1) for parameter calculation, we require the values of S(t) I(t), R(t), and D(t). These values have been derived from official reports of the Ministry of Health (MINSA) [11] spanning from March 5 to September 30, 2020. A subset of this data is presented in Table 2 for illustrative purposes.

Day	Date	Susceptible S(t)	Infected I(t)	Recovered R(t)	Deceased D(t)
0	Mar-05	32,625,948	0	0	0
1	Mar-06	32,625,947	1	0	0
12	Mar-17	32,625,831	116	0	1
26	Mar-31	32,624,883	630	394	41
27	Apr-01	32,624,625	829	447	47
33	Apr-07	32,622,994	1734	1100	120
40	Apr-14	32,615,645	7204	2869	230
47	Apr-21	32,608,111	10,371	6982	484
56	Apr-30	32,588,972	25,520	10,405	1051
57	May-01	32,585,489	28,206	11,129	1124
71	May-15	32,541,453	54,956	27,147	2392
87	Jun-01	32,455,909	96,898	68,507	4634
117	Jul-01	32,337,471	100,372	178,245	9860
148	Aug-01	32,203,765	111,940	290,835	19,408
179	Sept-01	31,962,511	154,001	480,177	29,259
208	Sept-30	31,807,651	95,234	690,528	32,535

Table 2.

Data related to COVID-19 in Peru (mar. 05 to sept. 30, 2020).

In order to calculate the parameters, we use the approximation of the derivative by central difference, implying that for each day of the sample data t = 1, 2, …, 208 = T, we calculate the parameters ${\lambda }_t$, ${\gamma }_t$, and ${\mu }_t$, which give us a sequence of data. Then, we calculate the simple average of such values to have the final parameters, as shown in Table 3.

t = 0	S(0)	I(0)	R(0)	D(0)	$\mathit{\lambda }$	$\mathit{\gamma }$	$\mathit{\mu }$
March 17	32,625,831	116	0	1	$4.244\times {10}^{-9}$	$4.576\times {10}^{-2}$	$1.990\times {10}^{-3}$

Table 3.

Parameters of the SIR model for COVID-19 in Peru (mar. 17 to sept. 30, 2020).

To obtain an approximate solution for the system (2), we consider $\mathrm{\Delta }{W}_l\sim {\zeta }_l\sqrt{\mathit{dt}}$, where ${\zeta }_l$ is a standard Gaussian random variable with mean 0 and variance 1. Subsequently, we discretize the system using the Euler–Maruyama, Milstein’s, and Runge–Kutta methods, as outlined in (8), (10), (12), and (13). The Python codes obtained for both the deterministic case and the stochastic system are presented in the following section.

4.2 Codes in python and simulations

Below, we present the complete codes for both deterministic and stochastic SIRD models.

4.2.1 Deterministic case

import numpy as np.

import matplotlib.pyplot as plt.

# Model parameters.

T0 = 0, T = 365, n = 36,500, dt = (T - T0) /n.

lambda_val = 0.00000000424404.

gamma = 0.0457600061.

mu = 0.0019902810.

initial_conditions = np.array ([32,625,831, 116, 0, 1]).

# Define the system of differential equations.

def sird_model (S, I, R, D):

dS = − lambda_val * S * I * dt.

dI = (lambda_val * S * I – gamma * I – mu * I) * dt.

dR = gamma * I * dt.

dD = mu * I * dt.

return np.array([dS, dI, dR, dD]).

# Initialize matrices to store results.

population_states = np.zeros((n, 4)).

population_states[0,:] = initial_conditions.

# Simulation of the model.

for i in range(1, n):

d_population = sird_model(*population_states[i – 1,:])

population_states[i,:] = population_states[i – 1,:] + d_population.

lot(np.arange(T0, T, dt), population_states[:, 0], label = ‘Susceptible’).

plt.plot(np.arange(T0, T, dt), population_states[:, 1], label = ‘Infected’).

plt.plot(np.arange(T0, T, dt), population_states[:, 2], label = ‘Recovered’).

plt.plot (np.arange(T0, T, dt), population_states [:, 3], label = ‘Deceased’).

plt.lege.

# Plotting results.

plt.pnd().

plt.grid(True).

plt.xlabel(‘Time’).

plt.ylabel(‘Population’).

plt.title(‘SIRD Model Simulation’).

plt.show().

We observe in Figure 2 that the simulation of the deterministic case using Python reproduces the same patterns obtained in MATLAB, which aligns with the literature [1, 2].

Figure 2.

Simulation results of the SIRD model in the deterministic case. The plot shows the evolution of the susceptible, infected, recovered, and deceased populations over time.

4.2.2 Euler–Maruyama code

import numpy as np.

import matplotlib.pyplot as plt.

# Model parameters.

T0 = 0, T = 365, n = 36,500, dt = (T - T0) / n.

lambda_val = 0.00000000424404.

gamma = 0.0457600061.

mu = 0.0019902810.

sigma1 = 0.05.

sigma2 = 0.05.

sigma3 = 0.05.

sigma4 = 0.05.

initial_conditions = np.array([32,625,831, 116, 0, 1]).

# Stochastic SIRD Model using Euler-Maruyama.

np.random.seed(100).

S = np.zeros(n).

I = np.zeros(n).

R = np.zeros(n).

D = np.zeros(n).

dW = np.sqrt(dt) * np.random.randn(n, 4).

S[0] = 32,625,831.

I[0] = 116.

R[0] = 0.

D[0] = 1.

def f(S, I, R, D, lambda_val, gamma, mu):

dS = −lambda_val * S * I.

dI = lambda_val * S * I – gamma * I – mu * I.

dR = gamma * I.

dD = mu * I.

return np.array([dS, dI, dR, dD]).

def g(X, sigma1, sigma2, sigma3, sigma4):

return np.array([[sigma1 * X[0], 0, 0, 0],

           [0, sigma2 * X[1], 0, 0],

           [0, 0, sigma3 * X[2], 0],

           [0, 0, 0, sigma4 * X[3]]])

for i in range (1, n):

# Euler-Maruyama.

X = np.array([S[i – 1], I[i – 1], R[i – 1], D[i – 1]]).

dX = f(*X, lambda_val, gamma, mu) * dt +.

g(X, sigma1, sigma2, sigma3, sigma4).dot(dW[i,:])

  X_new = X + dX.

  S[i], I[i], R[i], D[i] = X_new.

# Plotting results.

plt.figure (figsize = (12, 6)).

# Stochastic SIRD using Euler-Maruyama.

plt.plot(np.arange (0, T + dt, dt), [S[0]] + list (S),

label = ‘Susceptible (Stochastic)’, linestyle = ‘dashed’, color = ‘blue’).

plt.plot(np.arange(0, T + dt, dt), [I[0]] + list(I),

label = ‘Infected (Stochastic)’, linestyle = ‘dashed’, color = ‘orange’).

plt.plot(np.arange(0, T + dt, dt), [R[0]] + list(R),

label = ‘Recovered (Stochastic)’, linestyle = ‘dashed’, color = ‘green’).

plt.plot(np.arange(0, T + dt, dt), [D[0]] + list(D),

label = ‘Deceased (Stochastic)’, linestyle = ‘dashed’, color = ‘red’).

plt.legend(fontsize = ‘small’).

plt.grid(True).

plt.xlabel(‘Time’).

plt.ylabel(‘Population’).

plt.title(‘Deterministic vs. Euler-Maruyama’).

plt.show().

We can observe in Figure 3 that by using the values of ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.05$, the curve of recovered individuals overestimates that of the deterministic case. This phenomenon is evident when conducting an individual variable analysis in Figure 4 and when performing 50 simulations in Figure 5. We attribute this to the convergence order of the method, p = 0.5.

Figure 3.

Simulation results of the SIRD model using Euler–Maruyama method.

Figure 4.

Plots of the join three methods.

Figure 5.

50 simulations of each variable of the stochastic SIRD model using Euler–Maruyama method.

4.2.3 Milstein code

import numpy as np.

import matplotlib.pyplot as plt.

# Model parameters.

T0 = 0, T = 365, n = 36,500, dt = (T - T0) / n.

lambda_val = 0.00000000424404.

gamma = 0.0457600061.

mu = 0.0019902810.

sigma1 = 0.05, sigma2 = 0.05, sigma3 = 0.05, sigma4 = 0.05.

# Milstein method for the SIRD model.

S_milstein = np.zeros(n).

I_milstein = np.zeros(n).

R_milstein = np.zeros(n).

D_milstein = np.zeros(n).

dW = np.zeros(n).

S_milstein[0] = 32,625,831.

I_milstein[0] = 116.

R_milstein[0] = 0.

D_milstein[0] = 1.

dW[0] = np.sqrt(dt) * np.random.randn().

for i in range(1, n):

dW[i] = np.sqrt(dt) * np.random.randn().

S_milstein[i] = S_milstein[i - 1] - lambda_val * S_milstein[i - 1].

* I_milstein[i - 1] * dt + sigma1 * S_milstein[i - 1]* dW[i - 1].

+ 0.5 * (sigma1 ** 2) * S_milstein[i - 1] * (dW[i - 1] ** 2 - dt).

I_milstein[i] = I_milstein[i - 1] + (lambda_val * S_milstein[i - 1].

* I_milstein[i - 1] - gamma * I_milstein[i - 1] - mu * I_milstein[i - 1]).

+ sigma2 * I_milstein[i - 1] * dW[i - 1] + 0.5 * (sigma2 ** 2).

* I_milstein[i - 1] * (dW[i - 1] ** 2 - dt).

R_milstein[i] = R_milstein[i - 1] + gamma * I_milstein[i - 1] * dt.

+ sigma3 * R_milstein[i - 1] * dW[i - 1] + 0.5 * (sigma3 ** 2).

* R_milstein[i - 1] * (dW[i - 1] ** 2 - dt).

D_milstein[i] = D_milstein[i - 1] + mu * I_milstein[i - 1] * dt +.

sigma4 * D_milstein[i - 1] * dW[i - 1] + 0.5 * (sigma4 ** 2).

* D_milstein[i - 1] * (dW[i - 1] ** 2 - dt).

# Plotting results.

plt.figure(figsize = (12, 6)).

plt.plot(np.arange(T0, T, dt), [S_milstein[0]] + list(S_milstein[1:]),

linestyle = ‘dashed’, label = ‘Susceptible (Milstein)’).

plt.plot(np.arange(T0, T, dt), [I_milstein[0]] + list(I_milstein[1:]),

linestyle = ‘dashed’, label = ‘Infected (Milstein)’).

plt.plot(np.arange(T0, T, dt), [R_milstein[0]] + list(R_milstein[1:]),

linestyle = ‘dashed’, label = ‘Recovered (Milstein)’).

plt.plot(np.arange(T0, T, dt), [D_milstein[0]] + list(D_milstein[1:]),

linestyle = ‘dashed’, label = ‘Deceased (Milstein)’).

plt.legend(fontsize = ‘small’).

plt.grid(True).

plt.xlabel(‘Time’).

plt.ylabel(‘Population’).

plt.title(‘SIRD Model Simulation with Deterministic and Milstein Methods’).

plt.show().

In the case of Milstein with the same values for ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.05$, we observe in Figure 6 an underestimation in the recovered individuals. In the individual analysis depicted in Figure 4, it is noted that during a segment of its trajectory, an overestimation occurs, while in another segment, an underestimation is observed. Finally, in Figure 7, conducting 50 simulations reveals that, collectively, the trajectories are underestimated compared to the deterministic case. We attribute this phenomenon to the inherent convergence order of the method, p = 1.

Figure 6.

Simulation results of the SIRD model using Milstein method.

Figure 7.

50 simulations of each variable of the stochastic SIRD model using Milstein method.

4.2.4 Runge: Kutta 4 code

import numpy as np.

import matplotlib.pyplot as plt.

# Set random seed for reproducibility.

np.random.seed(100).

# Model parameters for the SIRD system.

T0 = 0, T = 365, n = 36,500, dt = (T - T0) / n.

lambda_val = 0.00000000424404.

gamma = 0.0457600061.

mu = 0.0019902810.

sigma1 = 0.05, sigma2 = 0.05, sigma3 = 0.05, sigma4 = 0.05.

# Initialize arrays to store S, I, R, D.

S = np.zeros(n).

I = np.zeros(n).

R = np.zeros(n).

D = np.zeros(n).

# Initialize stochastic increments.

dW = np.sqrt(dt) * np.random.randn(4, n) # Now dW has 4 components.

# Initial conditions.

S[0] = 32,625,831, I[0] = 116, R[0] = 0, D[0] = 1.

# Define deterministic part of the system.

def deterministic_part(S, I, R, D):

  return np.array([−lambda_val * S * I,

          lambda_val * S * I - gamma * I - mu * I,

          gamma * I,

          mu * I]).

# Define stochastic part of the system.

def stochastic_part(S, I, R, D):

return np.array([sigma1 * S, sigma2 * I, sigma3 * R, sigma4 * D]).

# Stochastic Runge–Kutta method.

def stochastic_runge_kutta(X, f, g, dt, dW):

     stages = 4.

     A = np.array([[0, 0, 0, 0],

           [0.5, 0, 0, 0],

          [0, 0.5, 0, 0],

          [0, 0, 1, 0]])

     Gamma = np.array([[0, 0, 0, 0],

              [0.5, 0, 0, 0],

              [0, 0.5, 0, 0],

              [0, 0, 1, 0]])

     b = np.array([1 / 6, 1 / 3, 1 / 3, 1 / 6])

     q = np.array([1/6, 1/3, 1/3, 1/6]).

     c = np.array([0, 0.5, 0.5, 1]).

     X_hat = np.zeros((stages, len(X))).

     for i in range(stages):

    sum_f = sum(A[i, j] * f(*X_hat[j]) for j in range(stages)).

    sum_g = sum(Gamma[i, j] * g(*X_hat[j]) for j in range(stages)).

    X_hat[i] = X + dt * sum_f + dW[i] * np.sqrt(dt) * sum_g.

     sum_b = sum(b[j] * f(*X_hat[j]) for j in range(stages)).

     sum_q = sum(q[j] * g(*X_hat[j]) for j in range(stages)).

     X_new = X + dt * sum_b + dW[−1] * np.sqrt(dt) * sum_q.

     return X_new.

# Time-stepping loop.

for i in range(1, n):

     X = np.array([S[i - 1], I[i - 1], R[i - 1], D[i - 1]]).

     X_new = stochastic_runge_kutta(X, deterministic_part, stochastic_part,                dt, dW[:, i - 1]).

     S[i], I[i], R[i], D[i] = X_new.

# Plot results.

plt.plot(np.arange(0, T + dt, dt), [S[0]] + list(S), label = ‘Susceptible’).

plt.plot(np.arange(0, T + dt, dt), [I[0]] + list(I), label = ‘Infected’).

plt.plot(np.arange(0, T + dt, dt), [R[0]] + list(R), label = ‘Recovered’).

plt.plot(np.arange(0, T + dt, dt), [D[0]] + list(D), label = ‘Deceased’).

plt.legend().

plt.grid(True).

plt.xlabel(‘Time’).

plt.ylabel(‘Population’).

plt.title(‘Stochastic Runge–Kutta Method for SIRD Model’).

plt.show().

With the same values of ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.05$, in this case, we observe in Figure 8 that all four variables follow the deterministic path. This observation holds true both in the individual analysis depicted in Figure 4 and in the 50 simulations presented in Figure 9. Analyzing the three methods, we can conclude that the fourth-order Runge–Kutta method best fits the deterministic case. We attribute this to the stability and convergence order of the method, p = 1.5, [7].

Figure 8.

Simulation results of the SIRD model using Runge Kutta 4 stochastic method.

Figure 9.

50 simulations of each variable of the stochastic SIRD model using RK4 method.

In Figure 4, we present a comparison for each variable of the stochastic SIRD model using the three methods explained in this chapter: Euler–Maruyama (yellow dashed lines), Milstein (green dashed lines), Runge–Kutta (red dashed lines), and the deterministic simulation (solid blue line).

In Figures 5–9, we conducted 50 simulations for each method to examine the average behavior of the simulations. It is worth noting that in some cases, stochastic models may provide a better description of the variable over extended periods, as seen in the case of recovered, for example.

5. Conclusion

5.1 Theoretical insights and stochastic methods

This section delves into theoretical results for stochastic methods, including discussions on the strong and weak convergence of the Euler–Maruyama and Milstein methods. Additionally, it explores specific conditions for the convergence of SRK methods, contingent on the deterministic Runge–Kutta methods’ convergence.

5.2 Application to COVID-19 modeling in Peru

The study applies the developed methodologies to a specific case: modeling the spread of COVID-19 in Peru. Official health reports are utilized to calculate parameters for the SIR model. The section presents deterministic and stochastic simulations of the SIRD model, offering insights into population dynamics over time.

5.3 Implementation and comparative analysis

Providing practical implementation, this part offers Python code snippets for various stochastic methods, including Euler–Maruyama, Milstein, and stochastic Runge–Kutta. The concluding part includes comparative visualizations between deterministic and stochastic simulations, emphasizing the impact of stochastic components on the model’s outcomes.

This comprehensive summary covers theoretical foundations, practical applications in COVID-19 modeling, and detailed code snippets, offering a thorough overview of the research.

Conflict of interest

The authors declare no conflict of interest.