SDEProblem driven by two different NoiseProcess(es)

Specific Domains Modelling & Simulations
,
1

How can we pass two different driving noise processes for an SDEProblem? The goal is to simulate a simple system of SDE’s:

dX(t) = f(X, Y, t)dt + g(X, Y, t)dZ_1(t) + h(X, Y, t)dZ_2(t)
dY(t) = \alpha(Y,t)dt + \beta(Y,t)dZ_2(t)

To exemplify with a concrete example of a process with stochastic volatility:

dX(t) = bdt +Y(t)dZ_1(t) - \rho dZ_2(t)
dY(t) = \eta (\mu - Y(t))dt + dZ_2(t)
Where Z_2(t) is some subordinator (always positive and increasing).

Let’s assume we are using custom noise processes declared using the interface from DiffEqNoiseProcess, so no possibility of exploring the multivariate properties of Wiener processes. I am using my own noise processes here, let’s call them MyProcess1 for Z_1 and MyProcess2 for Z_2.

Both have completely different properties so I don’t see how to do this without creating a single NoiseProcess equivalent to Z(t) = [Z_1(t), Z_2(t)], which seems to imply a tremendous amount of juggling to define dist and bridge functions that act in different ways on each coordinate.

function f(du,u,p,t)
du[1] = b
du[2] = eta * (mu - u[2])
end

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

And then we would define our NoiseProcess in the way that I said above (quite troublesome in practice).
Is there any handier way of doing this? Something similar in principle to:

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

function g2(du, u, p ,t)
du[1] = -rho
du[2] = 1
end

SDEProblem(f, g1, g2, u0, tspan; noise=(MyProcess1, MyProcess2))
2

There isn’t. That’s a cool idea but right now it cannot be done.

3

Got you. For anyone curious, here’s a simple MWE on how to implement this the laborious way, based on the system in the first post and using a difference of gamma’s for the first process and a gamma subordinator for the second. I think this a generalization of a variance gamma process with independent gamma-driven OU volatility.

using DifferentialEquations
using DiffEqNoiseProcess
using Distributions
using Plots

tspan = (0.0, 1.0)

#just to have less boring parameter values
addnoise(x) = x * (1 + 0.1*(rand() - 0.5))
αplus, αminus, αsub = [addnoise(10) for i=1:3]
θplus, θminus, θsub = [addnoise(1) for i=1:3]


function TestDualNoiseProcess!()
    function multidist!(dX, X, dt, u, p, t, rng)
        dX[1] = rand(rng, Gamma(αplus*dt, θplus)) - rand(rng, Gamma(αminus*dt, θminus))
        dX[2] = rand(rng, Gamma(αsub*dt, θsub))
        return dX
    end
    noiseprocess = NoiseProcess{true}(0.0, zeros(2), nothing,
                                      multidist!, nothing;
                                      continuous=false
                                      )
    return noiseprocess
end

testnoise = TestDualNoiseProcess!()
testnoise_prob = NoiseProblem(testnoise, tspan)
testnoise_solve = solve(testnoise_prob, alg=EM(), dt=0.001)
# plot(testnoise_solve.t, transpose(reduce(hcat, testnoise_solve.u)))


function f!(du, u, p, t)
du[1] = 0.0 # b
du[2] = 1 - u[2]
end

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


testsde_prob = SDEProblem(f!, g!, [0.0, 1.0], tspan;
                          noise_rate_prototype = zeros(2,2),
                          noise=testnoise
                          )

testsde_solve = solve(testsde_prob, alg=EM(), dt=0.001)

plot(testsde_solve)

Here’s a plot of a realization for this particular system:

Screenshot from 2023-04-19 10-29-27
Screenshot from 2023-04-19 10-29-271533×998 47.7 KB

To somewhat automatize this, I’m trying to create a function that creates a vector noise from independent marginals of a specific kind, fed as a collection of these marginals. Or making it even more general, by building different NoiseProcesses with independent components by making their dists/bridges act separately on each corresponding coordinate. I’ll try to post a draft here. It would still require the user to define the noise integrand g(du,u,p,t) in a consistent manner.

2 Likes