# Numerically solving a system of SDE's with Levy noise?

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

\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}]]}]

• Could you add code with parameters to reproduce Figure 1? Nov 28, 2022 at 15:52
• @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.
– Math
Nov 28, 2022 at 16:50
• You don't understand. Just add line of code with parameters from Table 1, for example A=0.9; ... Nov 28, 2022 at 16:57
• 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
– bmf
Jan 15, 2023 at 13:20
• 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? Jan 17, 2023 at 10:10

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}]


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"}]}]


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}]]}]



• This is very good Sir! But one question; what values were taken for the plots attached in your answer?
– Math
Jan 19, 2023 at 12:52
• 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?
– Math
Jan 19, 2023 at 13:22
• @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. Jan 19, 2023 at 14:56
• 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?
– Math
Jan 19, 2023 at 15:08
• @Math It is just parameters from Table 1. For plot above we use first column in a form rul[1]. Jan 19, 2023 at 15:14