7
$\begingroup$

Consider this system from the following paper titled: The long-time behaviour of a stochastic SIR epidemic model with distributed delay and multidimensional L´evy jumps

https://arxiv.org/pdf/2003.08219.pdf

Their system reads:

\begin{align} dS&=\big[A-\mu_1 S(t)-\beta S(t) D(t)\big]dt+\sigma_1 S(t) dW_1(t)+\int_U \lambda_1(u)S(t-)\tilde{N}(dt,du)\\[1ex] dI&=\big[\beta S(t) D(t)-(\mu_2+\gamma)I(t)\big]dt+\sigma_2 I(t) dW_2(t)+\int_U \lambda_2(u)I(t-)\tilde{N}(dt,du)\\[1ex] dD&=\big[\eta(I(t)-D(t)]dt+\sigma_4 D(t) dW_4(t)+\int_U \lambda_4(u)D(t-)\tilde{N}(dt,du) \end{align}

where $S(t-), I(t-)$ and $D(t-)$ are left limits of $S(t), I(t)$ and $D(t)$. $W_i$ are standard Brownian motions with $\sigma_i>0$. $N$ is a Poisson counting measure with compensating martingale $\tilde{N}$ and characteristic measure $\nu$ on a measurable subset $U$ of $(0,\infty)$, satisfying $\nu(U)<\infty$. We assume $\nu$ is a Levy measure such that $\tilde{N}=N(dt, du)-\nu(du)dt$ and assume $\lambda_i:Z\times\Omega\to \mathbb{R}$ are bounded and continuous.

How would I plot Figures 1-4 using the nominal values provided in the paper? I will add a bounty when applicable as this question is interesting to me.

Edit:

From Table 1, the value of parameters:

Figure 1 :

A = 0.9; \[Mu]1 = 0.3; \[Beta] = 0.07; \[Gamma] = 0.05; \[Mu]2 = 0.5; \
\[Eta] = 0.09; \[Sigma]1 = 0.15; \[Sigma]2 = 0.25; \[Sigma]4 = 0.27; \
\[Lambda]1 = 0.2; \[Lambda]2 = 0.23; \[Lambda]4 = 0.1;

Figure 2 :

A = 0.3; \[Mu]1 = 0.3; \[Beta] = 1.3; \[Gamma] = 0.05; \[Mu]2 = 0.5; \
\[Eta] = 0.09; \[Sigma]1 = 0.15; \[Sigma]2 = 0.25; \[Sigma]4 = 0.27; \
\[Lambda]1 = 0.2; \[Lambda]2 = 0.23; \[Lambda]4 = 0.1;

Figure 3 :

A = 0.6; \[Mu]1 = 0.4; \[Beta] = 0.35; \[Gamma] = 0.2; \[Mu]2 = 0.3; \
\[Eta] = 0.7; \[Sigma]1 = 0.2; \[Sigma]2 = 0.15; \[Sigma]4 = 0.13; \
\[Lambda]1 = 0.5; \[Lambda]2 = 0.7; \[Lambda]4 = 0.3;

Figure 4 :

A = 0.6; \[Mu]1 = 0.4; \[Beta] = 0.8; \[Gamma] = 0.3; \[Mu]2 = 0.3; \
\[Eta] = 0.2; \[Sigma]1 = 0.169; \[Sigma]2 = 0.15; \[Sigma]4 = 0.13; \
\[Lambda]1 = 0.5; \[Lambda]2 = 0.7; \[Lambda]4 = 0.3;

Deterministic model as in Figure 3:

tmax = 600;
\[Beta] = 0.35;
A = 0.6;
\[Mu]1 = 0.4;
\[Mu]2 = 0.3;
\[Eta] = 0.7;
\[Gamma] = 0.2;
SExDd = NDSolveValue[{
    S'[t] == A - \[Beta] *S[t]*Dd[t] - \[Mu]1*S[t],
    Ex'[t] == \[Beta] *S[t]*Dd[t] - (\[Mu]2 + \[Gamma])*Ex[t],
    Dd'[t] == \[Eta] (Ex[t] - Dd[t]),
    S[0] == 0.2,
    Ex[0] == 0.3,
    Dd[0] == 0.4},
   {S, Ex, Dd},
   {t, 0, tmax}];
{f1, f2, f3} = SExDd;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Red, Orange}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Red, Orange}, {"S(t)", "I(t)", "D(t)"}, 
    LegendFunction -> Framed], {0.85, 0.35}], ImageSize -> 550]

Deterministic model as in Figure 4:

tmax = 300;
\[Beta] = 0.8;
A = 0.6;
\[Mu]1 = 0.4;
\[Mu]2 = 0.3;
\[Eta] = 0.2;
\[Gamma] = 0.3;
SExDd = NDSolveValue[{
    S'[t] == A - \[Beta] *S[t]*Dd[t] - \[Mu]1*S[t],
    Ex'[t] == \[Beta] *S[t]*Dd[t] - (\[Mu]2 + \[Gamma])*Ex[t],
    Dd'[t] == \[Eta] (Ex[t] - Dd[t]),
    S[0] == 0.2,
    Ex[0] == 0.3,
    Dd[0] == 0.4},
   {S, Ex, Dd},
   {t, 0, tmax}];
{f1, f2, f3} = SExDd;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Red, Orange}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Red, Orange}, {"S(t)", "I(t)", "D(t)"}, 
    LegendFunction -> Framed], {0.85, 0.25}], ImageSize -> 550]

My question is; how do we extend the deterministic code to a stochastic one that includes the Levy term?

EDIT 2:

When we simulate this paper: https://www.researchgate.net/profile/Dianli-Zhao/publication/332241627_Threshold_dynamics_of_the_stochastic_epidemic_model_with_jump-diffusion_infection_force/links/5ca85d24a6fdcca26d013e72/Threshold-dynamics-of-the-stochastic-epidemic-model-with-jump-diffusion-infection-force.pdf, our plots, although similar, do not show the effects of Levy noise as what paper shows, any idea why?

The code for new paper:

µ(*Natural mortality rate of S,I,*)= {0.05, 0.05, 0.05, 0.4};
\[Beta] (*Transmission rate*)= {0.3, 0.3, 0.3, 0.8};
\[Delta](*Transmission rate*)= {0.05, 0.05, 0.05, 0.8};
\[Gamma] (*Recovered rate*)= {0.1, 0.1, 0.1, 0.3};
\[Sigma]1 (*Intensity of W1(t)*)= {0.1, 0.1, 0, 0.169};
\[Lambda]1 (*Intensity of W2(t)*)= {0.2, 0, 0, 0.15};


tmax = 801; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1], {0, tmax}];
ListStepPlot[{pL1}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;

eq1 = -s'[t] + (mu  - mu s[t] - beta s[t] i[t]) - 
   sigma1 s[t] i[t] dW1[t] + lambda1  L1[t];
eq2 = -i'[t] + (beta s[t] i[t] - (mu + delta + gamma) i[t]) + 
   sigma1 s[t] i[t] dW1[t] + lambda1  L1[t];

ic = {s[0] == 0.5, i[0] == 0.1};
rul[j_] := {beta -> \[Beta][[j]], gamma -> \[Gamma][[j]], 
   mu -> µ[[j]], delta -> \[Delta][[j]], sigma1 -> \[Sigma]1[[j]], 
   lambda1 -> \[Lambda]1[[j]]};

eqn[j_] := {eq1, eq2} /. rul[j];

sol[j_] := NDSolve[{eqn[j] == {0, 0}, ic}, {s, i}, {t, 0, tmax - 1}];

With[{sol = sol[1], sol2 = sol[2], 
  sol3 = sol[3]}, {Plot[
   Evaluate[s[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.85, .85}]], 
  Plot[Evaluate[i[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.8, .8}]]}]
Improve this question
$\endgroup$
21
  • $\begingroup$ Could you add code with parameters to reproduce Figure 1? $\endgroup$
    – Alex Trounev
    Nov 28, 2022 at 15:52
  • $\begingroup$ @AlexTrounev I tried coding in R, I can code the deterministic system and the system with Brownian motion but not with Levy noise. I cannot reproduce figure 1-4 unfortunately, hence my question. I did my other plots in my document using Mathematica(they were ODe's) but I don't know how to solve a SDE system. $\endgroup$
    – Math
    Nov 28, 2022 at 16:50
  • 3
    $\begingroup$ You don't understand. Just add line of code with parameters from Table 1, for example A=0.9; ... $\endgroup$
    – Alex Trounev
    Nov 28, 2022 at 16:57
  • 3
    $\begingroup$ I think that what @AlexTrounev is trying to point to you is that you should include in latex form the numerical values that we have to use in order to solve the differential equations. Don't code them, just write them down so we don't have to read the full paper $\endgroup$
    – bmf
    Jan 15 at 13:20
  • 1
    $\begingroup$ First of all you have to do the typing. Try to add some Mathematica code that does at least something. It will be easier for others to build on the existing code. You mention that you can code deterministic system. Why not in Mathematica? Is your post about Mathematical at all? $\endgroup$
    – yarchik
    Jan 17 at 10:10

1 Answer 1

Reset to default
8
+100
$\begingroup$

To simulate Levy jumps we can use numerical model described in the paper Analysis of a stochastic SEIS epidemic model with the standard Brownian motion and Lévy jump. In this model we use WhiteNoiseProcess[] and PoissonProcess[] to simulate diffusion and jumps as follows

SeedRandom[1234];
A (*Recruitment rate*)= { 0.9, 0.3, 0.6, 0.6};
µ1 (*Natural mortality rate of S*)= { 0.3, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*) = {0.07, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= { 0.05, 0.05, 0.2, 0.3};
µ2 (*General mortality of I*)= { 0.5, 0.5, 0.3, 0.3};
\[Eta] (*Exponentially fading memory rate*)= {0.09, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= { 0.15, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= { 0.25, 0.25, 0.15, 0.15};
\[Sigma]4 (*Intensity of W4(t)*)= { 0.27, 0.27, 0.13, 0.13};
\[Lambda]1 (*Jump intensity of S*)= { 0.2, 0.2, 0.5, 0.5};
\[Lambda]2 (*Jump intensity of I*)= { 0.23, 0.23, 0.3, 0.3};
\[Lambda]4 (*Jump intensity of D*)= { 0.1, 0.1, 0.7, 0.7};

tmax = 301; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}]; pL2 = 
 RandomFunction[PoissonProcess[1.1], {0, tmax}]; pL4 = 
 RandomFunction[PoissonProcess[.9], {0, tmax}]; ListStepPlot[{pL1, 
  pL2, pL4}]

Figure 1

Using data shown above we define 3 functions

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences; dpL2 = 
 pL2["SliceData", Range[tmax]] // First // Differences; dpL4 = 
 pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax; 
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax; 
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

With these functions the SDE model can be written in a form

eq1 = -s'[t] + (a - mu1 s[t] - beta s[t] d[t]) + sigma1 s[t] dW1[t] + 
   lambda1 s[t] L1[t];
eq2 = -i'[t] + (beta s[t] d[t] - (mu2 + gamma) i[t]) + 
   sigma2 i[t] dW2[t] + lambda2 i[t] L2[t];
eq3 = -r'[t] + (gamma i[t] - mu3 r[t]);
eq4 = -d'[t] + eta (i[t] - d[t]) + sigma4 d[t] dW4[t] + 
   lambda4 d[t] L4[t];

ic = {s[0] == 0.6, i[0] == 0.3, d[0] == 0.05}; 
rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], 
  gamma -> \[Gamma][[j]], eta -> \[Eta][[j]], mu1 -> µ1[[j]], 
  mu2 -> µ2[[j]], sigma1 -> \[Sigma]1[[j]], sigma2 -> \[Sigma]2[[j]], 
  sigma4 -> \[Sigma]4[[j]], lambda1 -> \[Lambda]1[[j]], 
  lambda2 -> \[Lambda]2[[j]], lambda4 -> \[Lambda]4[[j]]};

Numerical solution

eqn[j_]: = {eq1, eq2, eq4} /. rul[j];

sol[j_]: = NDSolve[{eqn[j] == {0, 0, 0}, ic}, {s, i, d}, {t, 0, tmax - 1}];

Visualization

With[{sol = sol[1]}, {Plot[Evaluate[s[t] /. sol], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> Green, 
   FrameLabel -> {"Time t (Days)", "S"}], 
  Plot[Evaluate[i[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Blue, 
   FrameLabel -> {"Time t (Days)", "I"}], 
  Plot[Evaluate[d[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Red, 
   FrameLabel -> {"Time t (Days)", "D"}]}]

Figure 2

Update 1. To compute model described in the paper we use Euler method and data

SeedRandom[123];

µ(*Natural mortality rate of S,I,*)= {0.05, 0.05, 0.05, 0.4};
\[Beta] (*Transmission rate*)= {0.3, 0.3, 0.3, 0.8};
\[Delta](*Transmission rate*)= {0.05, 0.05, 0.05, 0.8};
\[Gamma] (*Recovered rate*)= {0.1, 0.1, 0.1, 0.3};
\[Sigma]1 (*Intensity of W1(t)*)= {0.1, 0.1, 0, 0.169};
\[Lambda]1 (*Intensity of W2(t)*)= {0.1, 0, 0, 0.15};


tmax = 301; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1], {0, tmax}];



dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]];

eq1 = -s'[t] + (mu - mu s[t] - beta s[t] i[t]) - 
   s[t] i[t] (sigma1 dW1[t] - lambda1 L1[t]);
eq2 = -i'[t] + (beta s[t] i[t] - (mu + delta + gamma) i[t]) + 
   s[t] i[t] (sigma1 dW1[t] - lambda1 L1[t]);

ic = {s[0] == 0.5, i[0] == 0.1};
rul[j_] := {beta -> \[Beta][[j]], gamma -> \[Gamma][[j]], 
   mu -> µ[[j]], delta -> \[Delta][[j]], sigma1 -> \[Sigma]1[[j]], 
   lambda1 -> \[Lambda]1[[j]]};

eqn[j_] := {eq1, eq2} /. rul[j];

sol[j_] := 
  NDSolve[{eqn[j] == {0, 0}, ic}, {s, i}, {t, 0, tmax - 1}, 
   MaxSteps -> 10^6, Method -> "ExplicitEuler", 
   StartingStepSize -> 0.001];

With[{sol = sol[1], sol2 = sol[2], 
  sol3 = sol[3]}, {Plot[
   Evaluate[s[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.85, .85}]], 
  Plot[Evaluate[i[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.8, .8}]]}]

Figure 3

Improve this answer
$\endgroup$
11
  • $\begingroup$ This is very good Sir! But one question; what values were taken for the plots attached in your answer? $\endgroup$
    – Math
    Jan 19 at 12:52
  • $\begingroup$ And why is the output of this simulation much more "smooth" than the one in the paper: arxiv.org/pdf/2003.08219.pdf#page=19&zoom=100,56,369? $\endgroup$
    – Math
    Jan 19 at 13:22
  • $\begingroup$ @Math We don't know what pseudorandom generator state has been used in the paper while in my answer it is fixed as SeedRandom[1234]. All other parameters are taken from Table 1, Figure 1. $\endgroup$
    – Alex Trounev
    Jan 19 at 14:56
  • $\begingroup$ Is there a reason why there is a set for each parameter in the code? I presume the first value will be taken while the others are there for other figures, correct? $\endgroup$
    – Math
    Jan 19 at 15:08
  • $\begingroup$ @Math It is just parameters from Table 1. For plot above we use first column in a form rul[1]. $\endgroup$
    – Alex Trounev
    Jan 19 at 15:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.