Use Sobol sequence in DifferentialEquations.jl / StochasticDiffEq.jl

Specific Domains Modelling & Simulations
,
1

Hi everyone,

I’m trying to use the Sobol sequence to simulate Stochastic Differential Equations with StochasticDiffEq.jl / DifferentialEquations.jl. I’m trying to use the Sobol sequence to generate the random numbers because of the improved convergence rates that this sequence introduces in the simulation of this type of equations.

To that end I defined the following:

using DifferentialEquations 
using StatsFuns
using Sobol
using Plots

abstract type AbstractSobolSequence <: RandomNumbers.AbstractRNG{UInt32} end 

struct SobolSequence{T} <: AbstractSobolSequence 
    s::SobolSeq
end

function Random.randn(rng::SobolSequence, ::Union{Type{Float16}, Type{Float32}, Type{Float64}})    
    return norminvcdf(rng)
end


S0 = 111.0
r  = 0.0
σ  = 0.1
p  = (μ = r, σ = σ)

tspan = (0.0, 1.0)

sbl_rng  = SobolSequence(1, Float64)

W_sbl  = WienerProcess(0.0, 1.0, 1.0; rng=sbl_rng, reseed=false)

sde_sbl  = SDEProblem(f, g, S0, tspan, p; noise=W_sbl)

s1 = solve(EnsembleProblem(sde_sbl), trajectories=1000)
plot(s1)

which produces the following:

paths

which is clearly not correct.

I think that this issue could be solved by generating a SobolSequence of dimension n, where n is the number of time steps, but as DifferentialEquations.jl uses adaptive algorithms the number of time steps is not defined beforehand and I can’t do that.

To summarize, my questions :

  • Is it possible to use a Sobol sequence with DifferentialEquations.jl?

  • Is it possible to use a Sobol sequence with adaptive discretization algorithms?

  • If it’s possible which could be the strategy to generate the random numbers?

Thanks in advance!

2

Maybe my math is a bit rusty, but the Sobol sequence is primarily used in Quasi Monte-Carlo-type applications to generate advantageous pseudorandom numbers in a hypercube. Usually to estimate expectations/integrals such as \int_{[0,1]^d} f(x) \, \mathrm{d} x\approx \frac{1}{N} \sum_{j=1}^N f(x_j). For such types of applications Sobol numbers can give good convergence. In particular since the points have low discrepancy (cover the whole space quite uniformly everywhere).

But for time-stepping these numbers are not useful as they are correlated (exactly because of the low discrepancy property). That’s why you see the periodic patterns.

3

Thanks for your answer, but I’m not trying to time step the sequence. I’m trying to:

1- Generate numbers using the Sobol sequence
2- Use Sobol numbers to generate random normal numbers
3- Make DifferentialEquations.jl to simulate a SDE with these numbers instead of the default random numbers generator.

4

The point is that those random numbers are not independent, identically distributed. Hence the generated noise does not satisfy the axioms of a Wiener process.

1 Like
5

My tip is to use GitHub - SciML/DiffEqNoiseProcess.jl: A library of noise processes for stochastic systems like stochastic differential equations (SDEs) and other systems that are present in scientific machine learning (SciML) to define your own noise process. That way it’s easy to define the rng function without type piracy or creating a new method on randn.
If the numerical results are the same, the problem might lie elsewhere. Or there might not be a problem at all because of what @SteffenPL mentioned. Simulating an sde trajectory is time stepping. If you just want an integral over the distribution of terminal values, that’s a different thing and it is true that you can calculate it by Monte Carlo simulation of multiple sde trajectories and taking the samples of the desired time point; but there are other ways as well.

1 Like
6

Thanks for answer.

I’m actually doing that when I define the WiennerProcess and pass the SobolSequence as rng. In addition, I defined the randn function which is called in DiffEqNoiseProcess.jl when sampling the random numbers to generate the WienerProcess.

7

My suggestion was for you to define a different process with its own random sampling, not to replace the Wiener’s random sampling function. Check out Abstract Noise Processes · DiffEqNoiseProcess.jl

It should be equivalent to what you did but it could be a sanity check.