Solving Stochastic Differential Equations in Different Time Windows

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

1 Like

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

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.

1 Like

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

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

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)

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

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.

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

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.

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.

If you have to handle the Wiener processes explicitly than you can look into using https://github.com/mschauer/Bridge.jl :

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)

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.

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)

@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.

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)

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: https://github.com/sunoru/RandomNumbers.jl/issues/48 . 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.

3 Likes

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