# Stochastic model building: assistance with the DiffEq ecosystem

**URL:** <https://discourse.julialang.org/t/stochastic-model-building-assistance-with-the-diffeq-ecosystem/46611>\
**Category:** Numerics\
**Tags:** diffeq, sde\
**Created:** [September 14, 2020, 9:37pm UTC](https://discourse.julialang.org/t/stochastic-model-building-assistance-with-the-diffeq-ecosystem/46611 "2020-09-14T21:37:28Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [September 14, 2020, 9:37pm UTC](https://discourse.julialang.org/t/stochastic-model-building-assistance-with-the-diffeq-ecosystem/46611/1 "2020-09-14T21:37:29Z")

</div>

So, it’s time. I need to put some noise in my nice, clean, simple ODEs. Luckily DifferentialEquations.jl has me covered—the documentation is extensive… But perhaps too extensive: I don’t know where to start.

I _think_, ultimately, I’ll need to construct a [Jump Diffusion](https://diffeq.sciml.ai/v6.16/tutorials/jump_diffusion/#Jump-Diffusion) problem, but need some baby steps and assistance before I get there.

My system is a coupled ODE with three equations. One of those contains the term r \frac{N}{N + H} B\_t, where B\_t is the dependent variable `u[1]`, `r` and `H` are constants. `N` is also a constant, but is altered currently via callback at discrete time points.

Ignoring that last point for the moment, I’d like to add some noise to `N`. I can do this naively, changing `N::Float64` to `N::Normal(n, sigma)`, where `n` is the original value of `N`, then call `rand(N)` within my ODE function.

I’ve also tried something like the [scalar SDE example](https://diffeq.sciml.ai/v6.16/tutorials/sde_example/#Example-1:-Scalar-SDEs), where I set up the problem like:

```julia
function f(du, u, p, t)
   r, N, H = p #N is a Float64 here.
   du[1] = r*(N/(N+H))*u[1] + ...
  du[2] = ...
  du[3] = ...
end
function g(du, u, p, t)
   r, N, H = p
   du[1] = r*(N/(N+H))*u[1]    
   du[2] = 0.0    
   du[3] = 0.0                                                               
end
prob = SDEProblem(f,g,u0,tspan,p)

```

Here are the three results for `u[1]`:  
 ![stochastic](https://global.discourse-cdn.com/julialang/original/3X/8/9/89f670c444ed2906d08723e5df94143d3b9c3f6d.png)

At the moment, I don’t need to be too rigorous about the form of this noise, but I like the look of the SDEProblem’s solution better. However, the interpretation of what I’ve done in `g` is troubling me. Have I introduced a Weiner process with an amplitude of `r*(N/(N+H))*u[1]`? If I’m just looking to add some stochasticity to `N`, is this going about it the wrong way? If I set `du[1]` in `g` to just `N` or `N*u[1]` the results are non-nonsensical.

So, if I’ve got this part correct, is it a better solution than the `rand(N)` inside the ODE? As alternatives I tried using [Random ODEs](https://diffeq.sciml.ai/v6.16/tutorials/rode_example/#Example-1:-Scalar-RODEs) and generating my own [noise process](https://diffeq.sciml.ai/v6.16/features/noise_process/), but I don’t think either of them make sense for my problem. Basically I’m getting bogged down in the options.

I think I’m on the right path with the SDEProblem though, just need some assistance with a plain language interpretation of what I did there.

Next: assuming I have that right, I’d like to add in discrete jumps via a callback. Say `N=0.7` at one time point, 10 steps later a jump may change it to `N=1.2`, then perhaps 20 steps after that, down to `N=0.9`. I won’t know the size or direction of these jumps before time.

I’d like the Weiner process to take over in the steps between jumps however. This is why I suspect my system will ultimately become a [Jump Diffusion](https://diffeq.sciml.ai/v6.16/tutorials/jump_diffusion/#Jump-Diffusion) problem. Based on what I’ve said here, do you agree that this is the correct path for my process, or am I way off the mark?