solve SDEProblem with same driving noise · helpdesk (published)

Using DifferentialEquations.jl, I would like to solve several SDEs with the same driving Wiener noise. I am surprised to see that sol.u and sol2.u are different. Any suggestions how come?

function f2(du,u,p,t)
  du[1] = p[1] * (u[2] - u[1])
  du[2] = p[2] * (u[1] - u[2])
end

function g2(du,u,p,t)
  du[1] = 0.1
  du[2] = 0.1
end

tspan = (0.0, 2.0)
u0 = [0.0, 0.0]
p = [3.0, 1.0]
Z = WienerProcess(0.0, [0.0, 0.0])

prob = SDEProblem(f2, g2, u0, tspan, p,  noise=Z)
sol = solve(prob, EM(false), dt=0.01)

prob2 = SDEProblem(f2, g2, u0, tspan, p, noise=Z)
sol2 = solve(prob2, EM(false), dt=0.01)

sol.W == sol2.W # true
sol.u == sol2.u # false
sol.t == sol2.t # true

Frank Schäfer (Apr 06 2021 at 09:10):

To reproduce the same noise values with a NoiseProcess, you'll need to fix the seed:

using StochasticDiffEq, DiffEqNoiseProcess, Random
function f2(du,u,p,t)
  du[1] = p[1] * u[1]
end

function g2(du,u,p,t)
  du[1] = 0.1
end

tspan = (0.0, 2.0)
u0 = [0.0]
p = [0.98]
Z = WienerProcess(0.0, [0.0])

Random.seed!(10)
prob = SDEProblem(f2, g2, u0, tspan, p, noise=Z)
sol = solve(prob, EM(false), dt=0.01, save_noise=true)
Random.seed!(10)
prob2 = SDEProblem(f2, g2, u0, tspan, p, noise=Z)
sol2 = solve(prob2, EM(false), dt=0.01, save_noise=true)

sol.W.W == sol2.W.W # true
sol.u == sol2.u # true
sol.t == sol2.t # true

using Plots

plot(sol)
plot!(sol2)

However, I am also a bit surprised that sol.W.W in your example shows true.. I think that this is actually not the case and a new independent driving noise is drawn.

The correct way to re-use the noise values is to use NoiseWrapper from DiffEqNoiseProcess:

Zrep = NoiseWrapper(sol.W)
prob2= SDEProblem(f2, g2, u0, tspan, p, noise=Zrep)
sol2 = solve(prob2, EM(false), dt=0.01, save_noise=false)

plot(sol)
plot!(sol2)

Frank van der Meulen (Apr 06 2021 at 10:01):

That is confusing: you pass the noise to the function call, so why would you ever need to fix the seed? No random numbers need to be generated. I can work around this now of course, but the logic is not clear to me.

Looking at the "correct way" to do this, without save_noise=false something would change?

Frank Schäfer (Apr 06 2021 at 12:14):

I opened an issue regarding sol.W == sol2.W # true https://github.com/SciML/StochasticDiffEq.jl/issues/412. In your example, I'd have expected that in sol2 a new, independent noise process is simulated. So it's weird that the same noise values are displayed.

The random seed is just to fix the generation of the Wiener noise values (assuming that the noise values are fully regenerated).

I wouldn't expect a difference with save_noise=false; except that the noise values might grow a little because of floating point errors, i.e., if you step not exactly to the values of the time grid, a noise value will be sampled using the bridge distribution close to the time points of the grid, and this new noise value -- which will be very similar to the noise value on the grid -- would probably be stored as well.

Stream: helpdesk (published)

Topic: solve SDEProblem with same driving noise

Frank van der Meulen (Apr 04 2021 at 10:26):

Frank Schäfer (Apr 06 2021 at 09:10):

Frank van der Meulen (Apr 06 2021 at 10:01):

Frank Schäfer (Apr 06 2021 at 12:14):