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.377446336690395
  8.75489267338079
 13.13233901007118
 17.50978534676158

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

4-element Vector{Float64}:
  4.384131037074277
  8.768262074148554
 13.15239311122283
 17.536524148297108

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.