2

Looking for reference on doing numerical solving of a system of SDE's in python. I am trying to model a stochastic SIR model:

$$\begin{cases} dS = -\beta SIdt - \sigma SIdW \\ dI = (\beta SI - \gamma I)dt + \sigma SIdW \\ dR = \gamma Idt \end{cases}$$

I have been using a python package called SDEint, however, I am running into errors of getting negative sample paths which should not occur.

I have posted my code here: Stochastic SIR using SDEint python package

I have thought of trying to build my own solver but I am not terrible proficient in python. Any references to examples or literature that will help me in constructing my own solver is appreciated.

Additionally, perhaps there is a way to parametrize the system to avoid getting negative values in the sample paths. Such a parametrization is mentioned here How is the simulation done of the black scholes model?

Rodrigo de Azevedo
  • 1
oliverjones
  • 4,199
  • 2
    I have had the same problem using something like Euler-Maruyama in the setting of an SDE that has a decay in both the drift and the diffusion which prevents the true sample path from becoming negative but fails to do so with a numerical sample path. One fix is to directly track the evolution of the logarithm of the variable that you know remains positive. – Ian Dec 31 '21 at 21:02
  • @Ian could you explain how I could implement that fix? I believe I understand but don't see how I would immediately do that. – oliverjones Dec 31 '21 at 21:18
  • Look at $\ln(S)$ and $\ln(I)$ (you could look at $\ln(R)$ if you want but it is not really important) and write down SDEs for them using the Ito formula. – Ian Dec 31 '21 at 21:19
  • By the way, it's important here that the BMs in the first two SDEs are the same BM, not two independent BMs. This ensures that you still have conservation of population. – Ian Dec 31 '21 at 21:20
  • @Ian oh okay, yes I understand, I was thinking in terms of fixing my code not in of changing the SDE. Yes I will try that. Yeah, took a minute to figure out how to get the same BM in the package I am using. – oliverjones Dec 31 '21 at 21:21
  • 2
    Yeah, changing the algorithm to deal with this is actually rather subtle; doing things like reducing the time step when you get near the edge of the "admissible domain" in a naive way causes cross-correlations that you don't want. It is much easier to make a change of variable to something with "no wall". – Ian Dec 31 '21 at 21:26
  • @Ian I am still getting negative values, even when trying many different $\beta$ and $\gamma$ values. If I figure out a further fix I will post it. – oliverjones Jan 01 '22 at 00:51

1 Answers1

4

You can really ignore $R$ in this setting; if it is desired it can be reconstructed as $R=S_0+I_0+R_0-S-I$.

With that in mind, you can look at $f(S,I)=\ln(S)$ and $g(S,I)=\ln(I)$. You can use Ito's formula on each of these. For $f$ (assuming I didn't make any mistakes) you have

$$df=\left ( S^{-1} (-\beta S I) - \frac{1}{2} S^{-2} \sigma^2 S^2 I^2 \right ) dt - S^{-1} (\sigma S I) dW \\ =\left ( -\beta I - \frac{1}{2} \sigma^2 I^2 \right ) dt - \sigma I dW$$

and for $g$ you have

$$dg=\left ( I^{-1} (\beta S I) - I^{-1} (\gamma I) - \frac{1}{2} I^{-2} \sigma^2 S^2 I^2 \right ) dt + I^{-1} (\sigma S I) dW \\ =\left ( \beta S - \gamma - \frac{1}{2} \sigma^2 S^2\right ) dt + \sigma S dW$$

and now you can replace $S$ with $e^f$ and $I$ with $e^g$.

It is important here that $S$ and $I$ and thus $f$ and $g$ are driven by the same Brownian motion, so that the noise term conserves the total population.

With this formulation there is no more barrier to be concerned about crossing. The reconstructed solutions for $S$ and $I$ will not cross $0$ by definition since they will be written as $e^f$ and $e^g$ respectively.

Ian
  • 101,645