Expectation of Process Noise

SciMLExpectations.jl is able to calculate the average trajectory of a stochastic differential equation. This done by representing the Wiener process using the Kosambi–Karhunen–Loève theorem.

using SciMLExpectations
using Cuba
using Integrals
using StochasticDiffEq
using DiffEqNoiseProcess
using Distributions

f(du, u, p, t) = (du .= u)
g(du, u, p, t) = (du .= u)
u0 = collect(1:4)

W = WienerProcess(0.0, 0.0, 0.0)
prob = SDEProblem(f, g, u0, (0.0, 1.0), noise = W)
sm = ProcessNoiseSystemMap(prob, 8, LambaEM(), abstol = 1e-3, reltol = 1e-3)
cov(x, u, p) = x, p
observed(sol, p) = sol[:, end]
exprob = ExpectationProblem(sm, observed, cov)
sol1 = solve(exprob, Koopman(), ireltol = 1e-3, iabstol = 1e-3, batch = 64,
    quadalg = CubaDivonne())
sol1.u

4-element Vector{Float64}:
  4.366804481354138
  8.733608962708276
 13.10041344406242
 17.467217925416552

sol2 = solve(exprob, MonteCarlo(1_000_000))
sol2.u

4-element Vector{Float64}:
  4.383326100977312
  8.766652201954624
 13.149978302931933
 17.53330440390925

In theory, any numerical integration method from Integrals.jl is supported, but in practice many of the techniques struggle with correctly calculating the expected value. We got the best results with the Divonne algorithm from Cuba.jl.