Diagonal noise with specified noise process in solve (Ornstein Uhlenbeck)

Specific Domains Modelling & Simulations
,
1

Hello, I am trying to simulate a SDE problem with diagonal noise specifically Orstein Uhlenbeck noise, however I am only able to get scalar noise when specifying the noise process in SDEProblem (i.e. one noise process instead of independent noise processes for each DE). I have tried to resolve this issue by going through the examples in the docs (Stochastic Differential Equations · DifferentialEquations.jl). Specifically I have tried to extend Example 2 in the tutorial by adding a noise process in the SDEProblem formulation both a WienerProcess and the desired OrnsteinUhlenbeckProcess.

using DifferentialEquations
using Plots
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

function σ_lorenz(du, u, p, t)
    du[1] = 3.0
    du[2] = 3.0
    du[3] = 3.0
end

prob_sde_lorenz = SDEProblem(lorenz, σ_lorenz, [1.0, 0.0, 0.0], (0.0, 10.0))
sol = solve(prob_sde_lorenz,  save_noise=true)
println(size(sol.W))

W = WienerProcess(0.0,0.0, 0.0)
prob_sde_lorenz2 = SDEProblem(lorenz, σ_lorenz, [1.0, 0.0, 0.0], (0.0, 10.0), noise=W)
sol2= solve(prob_sde_lorenz2, save_noise=true)
println(size(sol2.W))


OU = OrnsteinUhlenbeckProcess(2.,0.,1.,0.,0.)
prob_sde_lorenz3 = SDEProblem(lorenz, σ_lorenz, [1.0, 0.0, 0.0], (0.0, 10.0), noise=OU)
sol3 = solve(prob_sde_lorenz3,EM(), adaptive=false, dt=0.001, save_noise=true)
println(size(sol3.W))

where the output is:

(3, 19852)
(857,)
(10001,)

How can I specify the noise process in such a way that I get an independent noise process for each equation with an Ornstein Uhlenbeck noise process?

2

The best way to do this is to represent OU via SDE terms and use those SDE terms in the other equations.

3

Thanks for the answer Chris. I assume that you mean using a matrix of g such as the non-diagonal noise example in the docs (Example 4 in Stochastic Differential Equations · DifferentialEquations.jl), but have g such that it is a diagonal matrix.

However, I am having issues when trying to specify the noise process in that example and I’m not really sure why

using DifferentialEquations
f(du, u, p, t) = du .= 1.01u

# Define a sparse matrix by making a dense matrix and setting some values as not zero
using SparseArrays
A = zeros(2, 4)
A[1, 1] = 1
A[1, 4] = 1
A[2, 4] = 1
A = sparse(A)

# Make `g` write the sparse matrix values
function g(du, u, p, t)
    du[1, 1] = 0.3u[1]
    du[1, 4] = 0.12u[2]
    du[2, 4] = 1.8u[2]
end

# Make `g` use the sparse matrix
prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = A)
sol = solve(prob,  save_noise=true)
println(size(sol.W))


# specify noise process as Wiener Process
W = WienerProcess(0.0,0.0, 0.0)
prob2 = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = A, noise=W)
sol2 = solve(prob2,  save_noise=true)
println(size(sol2.W))

With the follow output and stack trace:

(4, 408)
ERROR: DimensionMismatch: first array has length 2 which does not match the length of the second, 8.
Stacktrace:
  [1] generic_mul!(C::Vector{Float64}, X::SparseMatrixCSC{Float64, Int64}, s::Float64, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra C:\Users\nilsk\AppData\Local\Programs\julia-1.8.2\share\julia\stdlib\v1.8\LinearAlgebra\src\generic.jl:113
  [2] mul!
    @ C:\Users\nilsk\AppData\Local\Programs\julia-1.8.2\share\julia\stdlib\v1.8\LinearAlgebra\src\generic.jl:144 [inlined]
  [3] mul!
    @ C:\Users\nilsk\AppData\Local\Programs\julia-1.8.2\share\julia\stdlib\v1.8\LinearAlgebra\src\matmul.jl:276 [inlined]
  [4] perform_step!(integrator::StochasticDiffEq.SDEIntegrator{LambaEM{true}, true, Vector{Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SparseMatrixCSC{Float64, Int64}}, LambaEM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats, Nothing}, StochasticDiffEq.LambaEMCache{Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, Float64}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Float64, Nothing, Nothing}, cache::StochasticDiffEq.LambaEMCache{Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, Float64})
    @ StochasticDiffEq C:\Users\nilsk\.julia\packages\StochasticDiffEq\0PuS7\src\perform_step\lamba.jl:57
  [5] solve!(integrator::StochasticDiffEq.SDEIntegrator{LambaEM{true}, true, Vector{Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SparseMatrixCSC{Float64, Int64}}, LambaEM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats, Nothing}, StochasticDiffEq.LambaEMCache{Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, Float64}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Float64, Nothing, Nothing})
    @ StochasticDiffEq C:\Users\nilsk\.julia\packages\StochasticDiffEq\0PuS7\src\solve.jl:613
  [6] __solve(prob::SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SparseMatrixCSC{Float64, Int64}}, alg::LambaEM{true}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Nothing, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:default_set, :second_time, :save_noise), Tuple{Bool, Bool, Bool}}})
    @ StochasticDiffEq C:\Users\nilsk\.julia\packages\StochasticDiffEq\0PuS7\src\solve.jl:7
  [7] __solve(::SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SparseMatrixCSC{Float64, Int64}}, ::Nothing; default_set::Bool, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:second_time, :save_noise), Tuple{Bool, Bool}}})
    @ DifferentialEquations C:\Users\nilsk\.julia\packages\DifferentialEquations\yS3VA\src\default_solve.jl:9
  [8] #__solve#50
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\WXn2i\src\solve.jl:1097 [inlined]
  [9] #solve_call#26
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\WXn2i\src\solve.jl:473 [inlined]
 [10] #solve_up#32
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\WXn2i\src\solve.jl:839 [inlined]
 [11] solve(::SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEFunction{true, SciMLBase.FullSpecialize, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SparseMatrixCSC{Float64, Int64}}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:save_noise,), Tuple{Bool}}})
    @ DiffEqBase C:\Users\nilsk\.julia\packages\DiffEqBase\WXn2i\src\solve.jl:802
 [12] top-level scope
    @ REPL[23]:1
4

No, you can do this with diagonal noise.

function f(du,u,p,t)
  du[1] = -a * u[1]
  du[2] = 1.01u[1] + u[2] # u[2] is the OU term
end

function g(du,u,p,t)
  du[1] = b
  du[2] = 0
end