Solving Stochastic Differential Equations in Different Time Windows

zekeriya.sari · September 11, 2018, 11:58am

I want to solve a stochastic differential equation in different time windows under Wiener process, e.g., instead of solving the system for a time interval of (0, 2) seconds, I want to solve the system for a time duration of (0, 1) first, and then for a time duration of (1, 2). However, although I feed the problem in the second time window with the properly modified Wiener process (I update the initial time and initial value of the wiener process in the second time window with the last time and value of the wiener process in the first time window), the solutions of the system solved in two time windows and the that of the system solved at once do not coinside.

Below are the script I wrote for this problem and the plot of the solutions I obtained. Any ideas what I am doing wrong?

using DifferentialEquations
using Random
using Plots; plotly()

# Define the problem
f(dx, x, u, t) = (dx .= x)
g(dx, x, u, t) = (dx .= x)
dt = 2^(-6)

# Solve for 0-2 seconds
Random.seed!(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0.)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
Random.seed!(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0.)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end])
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

plt1 = plot(sol_full.t, vcat(sol_full.u...),
    linewidth=2, label="Full solution", linecolor=:red)
plot!(sol_part1.t, vcat(sol_part1.u...),
    linewidth=2, label="Solution part 1", linecolor=:green)
plot!(sol_part2.t, vcat(sol_part2.u...),
    linewidth=2, label="Solution part 2", linecolor=:blue)

plt2 = plot(sol_full.W.t, vcat(sol_full.W.u...),
    linewidth=2, label="Full solution Wiener", linecolor=:magenta)
plot!(sol_part1.W.t, vcat(sol_part1.W.u...),
    linewidth=2, label="Solution part 1 Wiener", linecolor=:cyan)
plot!(sol_part2.W.t, vcat(sol_part2.W.u...),
    linewidth=2, label="Solution part 2 Wiener", linecolor=:orange)

plot(plt1, plt2, layout=(2, 1))
gui()

solving_sde_in_time_windows

ChrisRackauckas · September 11, 2018, 3:18pm

You’re using global RNG seeding while the noise processes themselves have local RNGs.

github.com

SciML/DiffEqDocs.jl/blob/master/docs/src/features/noise_process.md#direct-construction-of-a-noiseprocess

# [Noise Processes](@id noise_process)

Noise processes are essential in continuous stochastic modeling. The `NoiseProcess`
types are distributionally-exact, meaning they are not solutions of
stochastic differential equations and instead are directly generated according
to their analytical distributions. These processes are used as the noise term
in the SDE and RODE solvers. Additionally, the noise processes themselves can
be simulated and solved using the DiffEq common interface (including the Monte
Carlo interface).

For more details, see [DiffEqNoiseProcess.jl](https://docs.sciml.ai/DiffEqNoiseProcess/)

The docs have been updated to reflect this. So what you’ll want to do is something like:

using RandomNumbers
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng = rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end],rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

to have it continue with the same RNG.

zekeriya.sari · September 11, 2018, 3:41pm

Thank you so much for your reply. But, I think that SDEProblem type does not accept rng as an argument and I think I should write

rng = Xorshifts.Xoroshiro128Plus(0)
noise = WienerProcess(0., 0., rng=rng)

since NoiseProcess has rng field. Athough I applied your suggestion, it did not change much, still the solutions do not coinside.

using DifferentialEquations
using RandomNumbers
using Plots; plotly()

# Define the problem
f(dx, x, u, t) = (dx .= x)
g(dx, x, u, t) = (dx .= x)
dt = 2^(-6)

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

plt1 = plot(sol_full.t, vcat(sol_full.u...),
    linewidth=2, label="Full solution", linecolor=:red)
plot!(sol_part1.t, vcat(sol_part1.u...),
    linewidth=2, label="Solution part 1", linecolor=:green)
plot!(sol_part2.t, vcat(sol_part2.u...),
    linewidth=2, label="Solution part 2", linecolor=:blue)

plt2 = plot(sol_full.W.t, vcat(sol_full.W.u...),
    linewidth=2, label="Full solution Wiener", linecolor=:magenta)
plot!(sol_part1.W.t, vcat(sol_part1.W.u...),
    linewidth=2, label="Solution part 1 Wiener", linecolor=:cyan)
plot!(sol_part2.W.t, vcat(sol_part2.W.u...),
    linewidth=2, label="Solution part 2 Wiener", linecolor=:orange)

plot(plt1, plt2, layout=(2, 1))
gui()

solving_sde_in_time_windows_updated

ChrisRackauckas · September 11, 2018, 4:01pm

Yeah, that was a typo in my post. Edited to be fixed.

ChrisRackauckas · September 11, 2018, 4:02pm

You need to reset the RNG for the second run, like:

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

zekeriya.sari · September 11, 2018, 5:03pm

Resetting the RNG in the second run did not solve the problem, either. Here is the case.

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)  # Reset the RNG!
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

solving_sde_in_time_windows_updated2

ChrisRackauckas · September 11, 2018, 6:12pm

Oh you might need to set reseed=false.

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)  # Reset the RNG!
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

That should be it. It automatically reseeds so that way giving the same noise process can give different solves which is normally wanted with a Monte Carlo, but here that’s not what you want.

zekeriya.sari · September 11, 2018, 6:29pm

reseed=false made the solutions coincide, but this time I lost the WienerProcess. Here is the case,

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(0)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

solving_sde_in_time_windows_updated3

ChrisRackauckas · September 11, 2018, 7:06pm

Yeah you’re running into a bug here, but I’m not exactly sure what it is. I’ll need to sit down and really hammer through this. I’ll see if I find time later today.

zekeriya.sari · September 11, 2018, 7:19pm

I am really grateful for your efforts. I am going through the documentation and source codes of the ecosystem to find a clue to cause of the problem. Thank you for all.

mschauer · September 11, 2018, 7:42pm

If you have to handle the Wiener processes explicitly than you can look into using GitHub - mschauer/Bridge.jl: A statistical toolbox for diffusion processes and stochastic differential equations. Named after the Brownian Bridge. :

using Bridge
using Random
using Plots
# Define the problem
f(t, x) = x
g(t, x) = x
dt = 2^(-6)

# Solve for 0-2 seconds
Random.seed!(0)
x0 = 1/2
t = 0.:dt:2.
W = sample(t, Wiener())
X = solve(Euler(), x0, W, (f, g))


Random.seed!(0)
x0 = 1/2
t1 = 0.:dt:1.
t2 = 1.:dt:2.
W1 = sample(t1, Wiener())
W2 = sample(t2, Wiener())
X1 = solve(EulerMaruyama(), x0, W1, (f, g))
X2 = solve(EulerMaruyama(),  X1.yy[end]+0.1, W2, (f, g)) # shift by 0.1 to make path visible

plot(X)
plot!(X1)
plot!(X2)

zekeriya.sari · September 11, 2018, 8:37pm

Bridge.jl did the job. This is the first time I tried Bridge.jl and it seems a very good tool to deal with stochastic differential equations.

The thing is that I need a syntax for drift and diffusion terms of the SDE as f(dx, x, u, t), g(dx, x, u, t) or f(x, u, t) , g(x, u, t). Here x refers to the state, dx refers to the derivative of the state and u which is a function of t refers to the input and t is the current time of the state of a dynamical system modelled by stochastic differential equations. Currently, I do not know how to apply this syntax to the one required by Bridge.jl which, I think, requires drift and diffusion functions of the form

Bridge.b(t, x, P::OrnsteinUhlenbeck) = -P.β * x
Bridge.σ(t, x, P::OrnsteinUhlenbeck) = P.σ
Bridge.a(t, x, P::OrnsteinUhlenbeck) = P.σ^2

as its documentation says. Maybe I need more reading about Bridge.jl. Thank you for the advice.

mschauer · September 11, 2018, 9:41pm

I am not sure what you need. Maybe the following is meaningful for you:

# Controlled example

struct ControlledOU{S} <: ContinuousTimeProcess{Float64}
    t0::Float64
    dt::Float64
    u::Matrix{S}
end

ttoi(t, P) = 1 + round(Int, (t-P.t0)*P.dt)
function Bridge.b!(t, y, out, P::ControlledOU) 
    for i in 1:2
        out[i] = -0.2 * y[i] + P.u[i, ttoi(t, P)]
    end
    out
end
Bridge.σ!(t, y, dw, out, P::ControlledOU) = out .= y.*dw 


x0 = [1/2, 0.5]
dt = 2^(-6)


t = 0.:dt:10.
n = length(t)
u = rand(2,n)

P = ControlledOU(t[1], dt, u)

W = VSamplePath(t, zeros(2, n))
X = VSamplePath(t, zeros(2, n))


sample!(W, Wiener())
solve!(EulerMaruyama!(), X, x0, W, P)
plot(X)

zekeriya.sari · September 12, 2018, 8:28am

@mschauer I read the docs and studied your code. I think the example you gave solves a differential equation of the form.

dx1 = (-0.2 *x1 + u1)dt + x1dW
dx2 = (-0.2 *x2 + u2)dt + x2dW

But, if I am not wrong the input function u is pre-generated as a vector in the with the line

u = rand(2, n)

and then the value of u(t) is then indexed using the function ttoi. This is not the case I wanted. But,I do not want the input function to be pre-generated but be a function that can be called at any time without being indexed by considering the initial time t0 and time step dt of the Wiener process P. But, this is a very good example for a new starter of Bridge.jl. Thank you very much for the example.
By the way, I cannot find the docs for VSamplePath in the docs but it is in the source code. I wanted to let you know in case you missed it.

mschauer · September 13, 2018, 4:01pm

Thank you for the feedback! I am not sure if I understand but you can of course play a bit with the design:

struct ControlledOU{S,T} <: ContinuousTimeProcess{Float64}
    t0::Float64
    dt::Float64
    u::T
end

function Bridge.b!(t, y, out, P::ControlledOU) 
    u = P.u(t)
    for i in 1:2
        out[i] = -0.2 * y[i] + u[i]
    end
    out
end

P = ControlledOU(t[1], dt, sin)

…

ChrisRackauckas · September 30, 2018, 6:56am

I traced the issue back to being an upstream issue with RandomNumbers.jl on v1.0. If you use a seed of zero you always get zeros out of randn:

using RandomNumbers
rng = Xorshifts.Xoroshiro128Plus(0)
randn(rng,Float64) # 0.0

That issue was upstreamed: zero seeds give zeros · Issue #48 · JuliaRandom/RandomNumbers.jl · GitHub . So if you don’t use a seed of zero and you set reseed=false, then the random process is not reseeded and you repeat the values:

using StochasticDiffEq, RandomNumbers, DiffEqNoiseProcess

# Define the problem
f(dx, x, u, t) = (dx .= x)
g(dx, x, u, t) = (dx .= x)
dt = 2^(-6)

# Solve for 0-2 seconds
rng = Xorshifts.Xoroshiro128Plus(1)
x0 = [1 / 2]
tspan = (0., 2.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_full = solve(prob, EM(), dt=dt)

# Solve for 0-1 and 1-2 seconds
rng = Xorshifts.Xoroshiro128Plus(1)
x0 = [1 / 2]
tspan = (0., 1.)
noise = WienerProcess(0., 0., rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part1 = solve(prob, EM(), dt=dt)

x0 = sol_part1.u[end]
tspan = (sol_part1.t[end], 2.)
noise = WienerProcess(sol_part1.W.t[end], sol_part1.W.u[end], rng=rng, reseed=false)
prob = SDEProblem(f, g, x0, tspan, noise=noise)
sol_part2 = solve(prob, EM(), dt=dt)

using Plots
plt1 = plot(sol_full.t, vcat(sol_full.u...),
    linewidth=2, label="Full solution", linecolor=:red)
plot!(sol_part1.t, vcat(sol_part1.u...),
    linewidth=2, label="Solution part 1", linecolor=:green)
plot!(sol_part2.t, vcat(sol_part2.u...),
    linewidth=2, label="Solution part 2", linecolor=:blue)

plt2 = plot(sol_full.W.t, vcat(sol_full.W.u...),
    linewidth=2, label="Full solution Wiener", linecolor=:magenta)
plot!(sol_part1.W.t, vcat(sol_part1.W.u...),
    linewidth=2, label="Solution part 1 Wiener", linecolor=:cyan)
plot!(sol_part2.W.t, vcat(sol_part2.W.u...),
    linewidth=2, label="Solution part 2 Wiener", linecolor=:orange)

plot(plt1, plt2, layout=(2, 1))

res

Hope that helps. Sorry that took so long.

zekeriya.sari · September 30, 2018, 9:54am

Thank you for your time and all your efforts @ChrisRackauckas . This solved the problem.

ymardoukhi · November 24, 2023, 9:16pm

@ChrisRackauckas I was trying to reproduce this plot by just rerunning your script. But it seems there is again a mismatch between the Wiener process in the second window compared to the full solution. This is what I get.

Notice that the Wiener processes in the first and second windows follow the same trend, which is odd. I tried different seeds, and the behaviour is still the same. Here is the plot when the seed is 2.

I am using Julia 1.9.3; here is the version of the packages I use.

BenchmarkTools v1.3.2
DiffEqCallbacks v2.33.1
DiffEqNoiseProcess v5.19.0
DifferentialEquations v7.11.0
Distributions v0.25.103
ForwardDiff v0.10.36
Interpolations v0.14.7
RandomNumbers v1.5.3
StaticArrays v1.6.5
StochasticDiffEq v6.63.2

ChrisRackauckas · November 24, 2023, 10:14pm

What is the script you’re running? The script has changed in the last 5 years, this was pre Julia v1.0. Check the docs on DiffEqNoiseProcess.

ymardoukhi · November 26, 2023, 5:59pm

The script is the same one you posted as the solution.