How do I solve stochastic differential equations in Julia?

Question

I try to understand how to solve stochastic differential equations (SDEs) numerically (I have no experience in any language, but for some reasons I chose Julia). As a starting model, I decided to use Lotka-Volterra equations. I read manual and tutorial for DifferentialEquations.jl. While my model is a simple system of ODE:

everything works fine:

using DifferentialEquations
using Plots
function lotka(du,u,p,t);
    ????, ????, ????, ???? = p; 
    du[1] = ????*u[1]-????*u[1]*u[2];
    du[2] = ????*u[1]*u[2]-????*u[2];
end
u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1);
tspan = (0.0,50.0);
prob = ODEProblem(lotka,u0,tspan,p);
sol = solve(prob);
plot(sol,vars = (1,2))

Now I want to add Ornstein-Uhlenbeck noise:

The stupid straightforward solution is:

using DifferentialEquations
using Plots
function lotka(du,u,p,t);
    ????, ????, ????, ????, ????, ???? = p; 
    du[1] = ????*u[1]-????*u[1]*u[2]+u[3];
    du[2] = ????*u[1]*u[2]-????*u[2]+u[4];
    du[3] = -u[3]/????+sqrt((2.0*????^2.0/????))*randn();
    du[4] =  -u[4]/????+sqrt((2.0*????^2.0/????))*randn();
end
u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1,0.1,0.1);
tspan = (0.0,50.0);
prob = ODEProblem(lotka,u0,tspan,p);
sol = solve(prob);

But, as expected, it did not work, as solver is not for such SDE problem.

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/jrun/.julia/packages/DiffEqBase/ujEgN/src/integrator_interface.jl:116

I tried to read SDE Julia documentation, but without good example I could not understand how to deal with it. Moreover, my math background is poor, and it seems I did not understand notation correctly. How can I make this work for SDEs?

UPDATE:

Finally, I have the following code:

using DifferentialEquations,Plots;

function lotka(du,u,p,t);
    ????, ????, ????, ????, ????, ???? = p; 
    du[1] = ????*u[1]-????*u[1]*u[2];
    du[2] = ????*u[1]*u[2]-????*u[2];
end
function stoch(du,u,p,t)
  ????, ????, ????, ????, ????, ???? = p; 
  du[1] = 1 
  du[2] = 1
end
u0 = [15.0,1.0];
p = (0.3,0.05,0.7,0.1,0.9,0.1);
????, ????, ????, ????, ????, ???? = p; 
OU = OrnsteinUhlenbeckProcess(1/????, 0.0, ????, 0.0, 0.0);
tspan = (0.0,100.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p,noise = OU);
sol = solve(prob,EM(),dt = 1e-3, adaptive=false);

To close this topic I have two last questions:

1. Is it correct code? (I am afraid, that my noise is not diagonal in this case)

2. May I have different initial noise (W0) for each of two equations?

Answer 1

It seems you have an SDE where the first two terms are driven by the second two which are stochastic? In that case, you'd make lotka the deterministic terms:

function lotka(du,u,p,t);
    𝛼, 𝛽, 𝛾, 𝛿, 𝜎, 𝜏 = p; 
    du[1] = 𝛼*u[1]-𝛽*u[1]*u[2]+u[3];
    du[2] = 𝛿*u[1]*u[2]-𝛾*u[2]+u[4];
    du[3] = -u[3]/𝜏
    du[4] = -u[4]/𝜏
end

and then stoch the stochastic part:

function stoch(du,u,p,t)
  𝛼, 𝛽, 𝛾, 𝛿, 𝜎, 𝜏 = p; 
  du[1] = 0 
  du[2] = 0
  du[3] = sqrt((2.0*𝜎^2.0/𝜏))
  du[4] = sqrt((2.0*𝜎^2.0/𝜏))
end

This is now written in the form du = f(u,p,t)dt + g(u,p,t)dW. Notice that, just like how you don't write dt, you don't write the randn() since that's handled (quite more intricately) in the solver. Using these you define an SDEProblem and solve:

u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1,0.1,0.1);
tspan = (0.0,50.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p);
sol = solve(prob);

(You do need to be careful with this model though since it's not guaranteed to stay positive, so it may go unstable if there is too much noise. Not just the integrator, but the solution itself)

Quick Explanation of Why Using an ODE Solver Will Not Work

Just for clarity, the reason why what you were doing above will not work is for two reasons. For one, adaptivity makes assumptions based on solution properties. For example, a standard 5th order integrator assumes the solution is 5 times differentiable and uses that in its error estimates for choosing stepsizes. An SDE is 0 times differentiable: it's only epsilon differentiable for every epsilon<1/2. So already, the error estimates and time step choices are going to be wildly off when directly using an ODE method for SDEs.

But secondly, ODE solvers adapt using rejection sampling. They will choose a dt, give it a try, and if that fails, then reduce the dt. Here you are taking a random number within your derivative, and a 5th order integrator will call your derivative function 7 times, getting 7 different values for what it thinks the derivative is, calculate the error estimate, realize it's wildly off (since it's not even differentiable so the whole algorithm made bad assumptions), reduce the time step, and then take completely different random numbers so that the smaller dt ends up being a completely different trajectory. As you can tell, this whole thing is wildly unstable and doesn't actually solve the real SDE that we had in the first place.

You can get around this by being a lot smarter with how you take said random numbers, how you define the error estimates, and using high order methods that assume "Ito differentiability" (i.e. "differentiability" in terms of the SDE's components). This is described in the paper which is the foundation of the current crop of adaptive SDE solvers.

Answering Updates

For the code you have

using DifferentialEquations,Plots;

function lotka(du,u,p,t);
    𝛼, 𝛽, 𝛾, 𝛿, 𝜎, 𝜏 = p; 
    du[1] = 𝛼*u[1]-𝛽*u[1]*u[2];
    du[2] = 𝛿*u[1]*u[2]-𝛾*u[2];
end
function stoch(du,u,p,t)
  𝛼, 𝛽, 𝛾, 𝛿, 𝜎, 𝜏 = p; 
  du[1] = 1 
  du[2] = 1
end
u0 = [15.0,1.0];
p = (0.3,0.05,0.7,0.1,0.9,0.1);
𝛼, 𝛽, 𝛾, 𝛿, 𝜎, 𝜏 = p; 
OU = OrnsteinUhlenbeckProcess(1/𝜏, 0.0, 𝜎, 0.0, 0.0);
tspan = (0.0,100.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p,noise = OU);
sol = solve(prob,EM(),dt = 1e-3, adaptive=false);

This has 1 OU process added to both variables. If you want them to be diagonal, i.e. two independent OU processes, you use:

OU = OrnsteinUhlenbeckProcess(1/𝜏, 0.0, 𝜎, 0.0, [0.0,0.0]);

More generally, [0.0,0.0] are the starting points, and you change these to solve (2). For a small performance optimization you can choose:

OU = OrnsteinUhlenbeckProcess!(1/𝜏, 0.0, 𝜎, 0.0, [0.0,0.0]);

Both this formulation and the formulation using Brownian SDEs work. The SDE formulation work with adaptive methods but is a larger system, while the OUProcess formulation is a smaller system but only will work well with EM() and fixed time steps. Which is better will be problem-dependent, but I'd probably prefer the SDE in most cases. The OUProcess form will be much better on RODEs though.