Linear system with all implemented noises

This time we consider linear equations in two different ways. Either a series of scalar equations each with a different noise or a system of equations with a vector-valued noise composed of all the implemented noises.

The equation

More precisely, we consider the RODE

\[ \begin{cases} \displaystyle \frac{\mathrm{d}{X}_t}{\mathrm{d} t} = - \|{Y}_t\|^2 {X}_t + {Y}_t, \qquad 0 \leq t \leq T, \\ \left. {X}_t \right|_{t = 0} = {X}_0, \end{cases}\]

where $\{{Y}_t\}_{t\geq 0}$ is either a scalar noise with one of the many implemented noises or a vector-valued noise with all the noises, each as an independent component of the vector process.

Again, the target solution is construced by solving the system via Euler method at a much higher resolution.

Numerical approximation

Setting up the problem

First we load the necessary packages

using Plots
using Measures
using LinearAlgebra
using Random
using Distributions
using RODEConvergence

Then we set up some variables, as in the first example. First, the RNG:

rng = Xoshiro(123)

The time interval is given by the end points

t0, tf = 0.0, 1.0

and the mesh parameters are set to

ntgt = 2^18
ns = 2 .^ (6:9)

4-element Vector{Int64}:
  64
 128
 256
 512

and

nsample = ns[[1, 2, 3]]

3-element Vector{Int64}:
  64
 128
 256

The number of simulations for the Monte Carlo estimate is set to

m = 400

Now we define all the noise parameters

y0 = 0.2

ν = 0.3
σ = 0.5

μ = 0.3

μ1 = 0.5
μ2 = 0.3
ω = 3π
primitive_a = t -> μ1 * t - μ2 * cos(ω * t) / ω
primitive_b2 = t -> σ^2 * ( t/2 - sin(ω * t) * cos(ω * t) / 2ω )

λ = 5.0
dylaw = Exponential(0.5)

steplaw = Uniform(0.0, 1.0)

λ₀ = 3.0
a = 2.0
δ = 3.0

nr = 6
transport(t, r) = mapreduce(ri -> cbrt(sin(ri * t)), +, r) / length(r)
ylaw = Gamma(7.5, 2.0)

hurst = 0.6

The noise is, then, defined as a (vectorial) ProductProcess, where each coordinate is one of the implemented noise types:

noise = ProductProcess(
    WienerProcess(t0, tf, 0.0),
    OrnsteinUhlenbeckProcess(t0, tf, y0, ν, σ),
    GeometricBrownianMotionProcess(t0, tf, y0, μ, σ),
    HomogeneousLinearItoProcess(t0, tf, y0, primitive_a, primitive_b2),
    CompoundPoissonProcess(t0, tf, λ, dylaw),
    PoissonStepProcess(t0, tf, λ, steplaw),
    ExponentialHawkesProcess(t0, tf, λ₀, a, δ, dylaw),
    TransportProcess(t0, tf, ylaw, transport, nr),
    FractionalBrownianMotionProcess(t0, tf, y0, hurst, ntgt)
)

Both the Wiener and the Orsntein-Uhlenbeck processes are additive noises so the strong order 1 convergence for them is not a surprise based on previous work since the Euler method coincides with the Euler-Maruyama method which is known to be of order 1 for additive noises.

The geometric Brownian motion noise and the time-dependent linear Itô diffusion noise are multiplicative noises and were also known to yield strong order 1 convergence since the Euler method for the RODE coincides also with the Milstein method for SDEs, which is known to be of order 1.

Another way of addressing the geometric Brownian motion process as a noise for a RODE is to observe that it, say $\{G_t\}_t,$ is given explicitly in the form $G_t = g(t, W_t)$ for a Wiener process $\{W_t\}_{t}$, so the Euler method for the associated RODE coincides with the Euler method for an associated RODE with additive noise $\{W_t\}_t.$ However, the corresponding nonlinear term does not have global Lipschitz bound, so one needs to be careful with this interpretation. Our results, however, apply without further assumptions.

All the other noises, however, would be thought to have a lower order of convergence but our results prove they are also of order 1. Hence, their combination is also expected to be of order 1, as illustrated here.

Notice we chose the hurst parameter of the fractional Brownian motion process to be between 1/2 and 1, so that the strong convergence is also of order 1, just like for the other types of noises in noise. Previous knowledge would expect a strong convergence of order $H$, with $1/2 < H < 1,$ instead.

The system with all noises combined

Now we define the right hand side of the system of equations, in the in-place version, for the sake of performance. Here, the vector variable dx is updated with the right hand side of the system of equations. The norm squared of the noise y at a given time t is computed via sum(abs2, y).

f!(dx, t, x, y, p) = (dx .= .- sum(abs2, y) .* x .+ y)

params = nothing

The initial condition for the system takes into account the number of equations in the system, which is determined by the dimension of the vector-valued noise.

x0law = MvNormal(zeros(length(noise)), I)

IsoNormal(
dim: 9
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Σ: [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]
)

We now add some information about the simulation, for the caption of the convergence figure.

info = (
    equation = "\$\\mathrm{d}{X}_t/\\mathrm{d}t = - \\| {Y}_t\\|^2 {X}_t + {Y}_t\$",
    noise = "vector-valued noise \$\\{{Y}_t\\}_t\$ with all the implemented noises",
    ic = "\${X}_0 \\sim \\mathcal{N}({0}, \\mathrm{I})\$"
)

The method for which want to estimate the rate of convergence is, naturally, the Euler method, which is also used to compute the target solution, at the finest mesh:

target = RandomEuler(length(noise))
method = RandomEuler(length(noise))

RandomEuler{Float64, Distributions.Multivariate}([3.5e-323, 2.243e-321, 2.15e-321, 2.085e-321, 1.996e-321, 1.937e-321, 1.853e-321, 1.8e-321, 5.377940751268117e-299])

Order of convergence for the system with all the noises

With all the parameters set up, we build the ConvergenceSuite:

suite = ConvergenceSuite(t0, tf, x0law, f!, noise, params, target, method, ntgt, ns, m)

ConvergenceSuite{Float64, Distributions.IsoNormal, ProductProcess{Float64, Tuple{WienerProcess{Float64}, OrnsteinUhlenbeckProcess{Float64}, GeometricBrownianMotionProcess{Float64}, HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}, CompoundPoissonProcess{Float64, Distributions.Exponential{Float64}}, PoissonStepProcess{Float64, Distributions.Uniform{Float64}}, ExponentialHawkesProcess{Float64, Distributions.Exponential{Float64}}, TransportProcess{Float64, typeof(Main.var"Main".transport), Distributions.Gamma{Float64}, 1}, FractionalBrownianMotionProcess{AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, false, 1, UnitRange{Int64}}, Float64}, FFTW.cFFTWPlan{ComplexF64, -1, false, 1, Tuple{Int64}}}}}, Nothing, typeof(Main.var"Main".f!), 2, 2, RandomEuler{Float64, Distributions.Multivariate}, RandomEuler{Float64, Distributions.Multivariate}}(0.0, 1.0, IsoNormal(
dim: 9
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Σ: [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]
)
, Main.var"Main".f!, ProductProcess{Float64, Tuple{WienerProcess{Float64}, OrnsteinUhlenbeckProcess{Float64}, GeometricBrownianMotionProcess{Float64}, HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}, CompoundPoissonProcess{Float64, Distributions.Exponential{Float64}}, PoissonStepProcess{Float64, Distributions.Uniform{Float64}}, ExponentialHawkesProcess{Float64, Distributions.Exponential{Float64}}, TransportProcess{Float64, typeof(Main.var"Main".transport), Distributions.Gamma{Float64}, 1}, FractionalBrownianMotionProcess{AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, false, 1, UnitRange{Int64}}, Float64}, FFTW.cFFTWPlan{ComplexF64, -1, false, 1, Tuple{Int64}}}}}((WienerProcess{Float64}(0.0, 1.0, 0.0), OrnsteinUhlenbeckProcess{Float64}(0.0, 1.0, 0.2, 0.3, 0.5), GeometricBrownianMotionProcess{Float64}(0.0, 1.0, 0.2, 0.3, 0.5), HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}(0.0, 1.0, 0.2, Main.var"Main".var"#2#3"(), Main.var"Main".var"#5#6"()), CompoundPoissonProcess{Float64, Distributions.Exponential{Float64}}(0.0, 1.0, 5.0, Distributions.Exponential{Float64}(θ=0.5)), PoissonStepProcess{Float64, Distributions.Uniform{Float64}}(0.0, 1.0, 5.0, Distributions.Uniform{Float64}(a=0.0, b=1.0)), ExponentialHawkesProcess{Float64, Distributions.Exponential{Float64}}(0.0, 1.0, 3.0, 2.0, 3.0, Distributions.Exponential{Float64}(θ=0.5)), TransportProcess{Float64, typeof(Main.var"Main".transport), Distributions.Gamma{Float64}, 1}(0.0, 1.0, Distributions.Gamma{Float64}(α=7.5, θ=2.0), Main.var"Main".transport, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), FractionalBrownianMotionProcess{AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, false, 1, UnitRange{Int64}}, Float64}, FFTW.cFFTWPlan{ComplexF64, -1, false, 1, Tuple{Int64}}}(0.0, 1.0, 0.2, 0.6, 262144, [0.0, 1.0877824e-316, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im], ComplexF64[0.0 + 4.6609763e-317im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im], 1.9073486328125e-6 * FFTW backward plan for 524288-element array of ComplexF64
(dft-ct-dit/128
  (dftw-genericbuf/64-128-4096
    (dft-direct-128-x64 "n1fv_128_avx2"))
  (dft-indirect-before
    (dft-vrank>=1-x128/1
      (dft-ct-dit/8
        (dftw-direct-8/56 "t2bv_8_avx2")
        (dft-ct-dif/8
          (dftw-directsq-8/28-x8 "q1bv_8_avx2")
          (dft-vrank>=1-x8/1
            (dft-direct-64-x8 "n1bv_64_avx2")))))
    (dft-r2hc-1
      (rdft-rank0-tiled/2-x4096-x128)))), FFTW forward plan for 524288-element array of ComplexF64
(dft-ct-dit/128
  (dftw-genericbuf/64-128-4096
    (dft-direct-128-x64 "n1fv_128_avx2"))
  (dft-indirect-before
    (dft-vrank>=1-x128/1
      (dft-ct-dit/8
        (dftw-direct-8/56 "t2fv_8_avx2")
        (dft-ct-dif/8
          (dftw-directsq-8/28-x8 "q1fv_8_avx2")
          (dft-vrank>=1-x8/1
            (dft-direct-64-x8 "n1fv_64_avx2")))))
    (dft-r2hc-1
      (rdft-rank0-tiled/2-x4096-x128))))))), nothing, RandomEuler{Float64, Distributions.Multivariate}([3.5e-323, NaN, 6.9450912992678e-310, 6.9452102012903e-310, 6.94508999420486e-310, 6.94521268529986e-310, 6.94509117603044e-310, 3.177e-321, NaN]), RandomEuler{Float64, Distributions.Multivariate}([3.5e-323, 2.243e-321, 2.15e-321, 2.085e-321, 1.996e-321, 1.937e-321, 1.853e-321, 1.8e-321, 5.377940751268117e-299]), 262144, [64, 128, 256, 512], 400, [1, 1, 1, 1], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 2.4e-322; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 6.50432e-319 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [2.413417796e-315 2.177419237e-315 … NaN NaN; 2.343794223e-315 2.33639881e-315 … -7.962215409160601 2.22724724846113e-310; … ; 0.0 -4.909382477812359 … 2.46850319e-315 2.558617644914287e-309; 2.336398574e-315 -5.222434485176767 … 2.1186010639637e-310 2.56404995413916e-309])

Then we are ready to compute the errors via solve:

@time result = solve(rng, suite)

152.600349 seconds (5.46 G allocations: 81.490 GiB, 3.20% gc time, 1.33% compilation time)

The computed strong error for each resolution in ns is stored in result.errors, and a raw LaTeX table can be displayed for inclusion in the article:

table = generate_error_table(result, suite, info)

    \begin{center}
        \begin{tabular}[htb]{|r|l|l|l|}
            \hline N & dt & error & std err \\
            \hline \hline
            64 & 0.0156 & 0.205 & 0.0093 \\
            128 & 0.00781 & 0.0993 & 0.00451 \\
            256 & 0.00391 & 0.049 & 0.00223 \\
            512 & 0.00195 & 0.0243 & 0.00111 \\
            \hline
        \end{tabular}
    \end{center}

    \bigskip

    \caption{Mesh points (N), time steps (dt), strong error (error), and standard error (std err) of the Euler method for $\mathrm{d}{X}_t/\mathrm{d}t = - \| {Y}_t\|^2 {X}_t + {Y}_t$ for each mesh resolution $N$, with initial condition ${X}_0 \sim \mathcal{N}({0}, \mathrm{I})$ and vector-valued noise $\{{Y}_t\}_t$ with all the implemented noises, on the time interval $I = [0.0, 1.0]$, based on $M = 400$ sample paths for each fixed time step, with the target solution calculated with $262144$ points. The order of strong convergence is estimated to be $p = 1.024$, with the 95\% confidence interval $[0.8926, 1.1556]$.}

The calculated order of convergence is given by result.p:

println("Order of convergence `C Δtᵖ` with p = $(round(result.p, sigdigits=2)) and 95% confidence interval ($(round(result.pmin, sigdigits=3)), $(round(result.pmax, sigdigits=3)))")

Order of convergence `C Δtᵖ` with p = 1.0 and 95% confidence interval (0.893, 1.16)

Plots

We illustrate the rate of convergence with the help of a plot recipe for ConvergenceResult:

plt_result = plot(result)

For the sake of illustration, we plot a sample of the norms of a sequence of approximations of a target solution, along with the norm of the target:

plts = [plot(suite, ns=nsample, xshow=i, resolution=2^4, title="Coordinate $i", titlefont=8) for i in axes(suite.xt, 2)]

plot(plts..., legend=false)

We can also visualize the noises associated with this sample solution, both individually, as they enter the non-homogenous term,

plt_noises = plot(suite, xshow=false, yshow=true, linecolor=:auto, label=["W" "OU" "gBm" "hlp" "cP" "sP" "H" "T" "fBm"], legend=:topright)

and combined, with their sum squared, as it enters the homogenous term,

plot(suite, xshow=false, yshow= y -> sum(abs2, y), label="\$\\left\\|\\left\\|{Y}_t\\right\\|\\right\\|^2\$")

We finally combine all the convergence plot and the noises into a single plot, for the article:

plt_noises = plot(suite, xshow=false, yshow=true, linecolor=:auto, legend=nothing)

plt_combined = plot(plt_result, plt_noises, legendfont=6, size=(800, 240), title=["(a) non-homogeneous linear system" "(b) sample paths of all noises"], titlefont=10, bottom_margin=5mm, left_margin=5mm)

Scalar equations with the individual noises

Now we simulate a series of Random ODEs, each with one of the noises above, instead of the system with all combined noises.

In the univariate case, the right hand side of the equation becomes

f(t, x, y, p) = - y^2 * x + y

f (generic function with 1 method)

The initial condition is also univariate

eachx0law = Normal()

Distributions.Normal{Float64}(μ=0.0, σ=1.0)

and so is the Euler method

eachtarget = RandomEuler()
eachmethod = RandomEuler()

RandomEuler{Float64, Distributions.Univariate}(Float64[])

Now we compute the error for each noise and gather the order of convergence in a vector.

ps = Float64[result.p]
pmins = Float64[result.pmin]
pmaxs = Float64[result.pmax]
noises = String["all noises combined"]

for eachnoise in noise.processes
    eachsuite = ConvergenceSuite(t0, tf, eachx0law, f, eachnoise, params, eachtarget, eachmethod, ntgt, ns, m)

    @time eachresult = solve(rng, eachsuite)

    @info "noise = $(typeof(eachnoise)) => p = $(eachresult.p) ($(eachresult.pmin), $(eachresult.pmax))"

    push!(noises, string(nameof(typeof(eachnoise)))[begin:end-7])
    push!(ps, eachresult.p)
    push!(pmins, eachresult.pmin)
    push!(pmaxs, eachresult.pmax)
end

  1.526121 seconds (314.50 k allocations: 14.395 MiB, 26.47% compilation time)
[ Info: noise = WienerProcess{Float64} => p = 1.000500529941167 (0.8601306728605727, 1.1360999792781403)
  1.597391 seconds (299.86 k allocations: 14.615 MiB, 24.33% compilation time)
[ Info: noise = OrnsteinUhlenbeckProcess{Float64} => p = 0.9868297993393622 (0.8581873477311118, 1.115474089463054)
  2.016133 seconds (281.92 k allocations: 12.705 MiB, 18.53% compilation time)
[ Info: noise = GeometricBrownianMotionProcess{Float64} => p = 0.9907240757370498 (0.8063472331122628, 1.1773787722720688)
107.270371 seconds (5.45 G allocations: 81.263 GiB, 3.83% gc time, 0.37% compilation time)
[ Info: noise = HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"} => p = 0.989188275284611 (0.8257025635784274, 1.1529600780970857)
  1.478755 seconds (301.26 k allocations: 13.677 MiB, 26.69% compilation time)
[ Info: noise = CompoundPoissonProcess{Float64, Distributions.Exponential{Float64}} => p = 1.0192335984025374 (0.7424546579566343, 1.2722975913381867)
  2.008487 seconds (296.81 k allocations: 13.480 MiB, 19.21% compilation time)
[ Info: noise = PoissonStepProcess{Float64, Distributions.Uniform{Float64}} => p = 0.9861495954720401 (0.8352250747130328, 1.1372961112150073)
  2.019327 seconds (316.24 k allocations: 14.402 MiB, 20.17% compilation time)
[ Info: noise = ExponentialHawkesProcess{Float64, Distributions.Exponential{Float64}} => p = 1.021105203493439 (0.9107408153334632, 1.1314144614894124)
 13.869921 seconds (289.55 k allocations: 13.087 MiB, 2.77% compilation time)
[ Info: noise = TransportProcess{Float64, typeof(Main.var"Main".transport), Distributions.Gamma{Float64}, 1} => p = 1.067703318808653 (0.9786497787079639, 1.1569229440927251)
 21.432310 seconds (809.26 k allocations: 40.730 MiB, 2.72% compilation time)
[ Info: noise = FractionalBrownianMotionProcess{AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, false, 1, UnitRange{Int64}}, Float64}, FFTW.cFFTWPlan{ComplexF64, -1, false, 1, Tuple{Int64}}} => p = 0.9834867097957234 (0.8402845335481829, 1.1259490480462113)

We print them out for inclusing in the paper:

noises_short = ["all"; "W"; "OU"; "gBm"; "hlp"; "cP"; "sP"; "H"; "T"; "fBm"]

for (noisej, noiseshortj, pj, pminj, pmaxj) in zip(noises, noises_short, ps, pmins, pmaxs)
    println("$noisej ($noiseshortj) & $(round(pj, sigdigits=6)) & $(round(pminj, sigdigits=6)) & $(round(pmaxj, sigdigits=6)) \\\\")
end

all noises combined (all) & 1.02409 & 0.892626 & 1.15557 \\
Wiener (W) & 1.0005 & 0.860131 & 1.1361 \\
OrnsteinUhlenbeck (OU) & 0.98683 & 0.858187 & 1.11547 \\
GeometricBrownianMotion (gBm) & 0.990724 & 0.806347 & 1.17738 \\
HomogeneousLinearIto (hlp) & 0.989188 & 0.825703 & 1.15296 \\
CompoundPoisson (cP) & 1.01923 & 0.742455 & 1.2723 \\
PoissonStep (sP) & 0.98615 & 0.835225 & 1.1373 \\
ExponentialHawkes (H) & 1.02111 & 0.910741 & 1.13141 \\
Transport (T) & 1.0677 & 0.97865 & 1.15692 \\
FractionalBrownianMotion (fBm) & 0.983487 & 0.840285 & 1.12595 \\

The following plot helps in visualizing the result.

plt_eachnoise = plot(ylims=(-0.1, 1.5), ylabel="\$p\$", guidefont=10, legend=:bottomright)
scatter!(plt_eachnoise, noises_short, ps, yerror=(ps .- pmins, pmaxs .- ps), xrotation = 30, label="computed")
hline!(plt_eachnoise, [1.0], linestyle=:dash, label="theory",bottom_margin=5mm, left_margin=5mm)

Strong order $p$ of convergence of the Euler method for $\mathrm{d}X_t/\mathrm{d}t = - Y_t^2 X_t + Y_t$ for a series of different noise $\{Y_t\}_t$ (scattered dots: computed values; dashed line: expected $p = 1;$ with 95% confidence interval).

Combined with the noise sample paths:

plt_combined = plot(plt_eachnoise, plt_noises, legendfont=6, size=(800, 240), title=["(a) non-homogeneous linear system" "(b) sample paths of all noises"], titlefont=10, bottom_margin=5mm, left_margin=5mm)

This page was generated using Literate.jl.