Mathematical validity of using SDE solver schemes for custom NoiseProcess(es) in DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

As we know, DifferentialEquations.jl has the amazing feature of allowing us to define our own NoiseProcess object and use all the solving interface with it. My question here is more related to the mathematical validity of using the solving schemes available for SDEProblems with custom noise.

I’m currently reading the famous book “Platen, Bruti-Liberati (2010) - Numerical Solution of Stochastic Differential Equations with Jumps in Finance” to see if I can find any hints on this (haven’t so far), but I’m creating this thread simultaneously because I think this might be an interesting topic.

So, we know many of the higher order (strong) convergence schemes (i.e., better than 0.5 or Euler-Maruyama) rely on properties of Brownian motion that simplify the computations in the Taylor expansions for the numerical schemes. The Milstein scheme is a basic example of this.

Now, consider we want to define our own NoiseProcess that defines a pure jump Lévy process (for example). For some processes, we have access to methods for sampling the increment distribution directly; a very good example for this is the alpha stable process. The common approach for simulating SDE’s when we have a means of sampling the increment distribution is to use the Euler-Maruyama scheme with the noise term sampling Noise(t + delta) - Noise(t); which is equal in distribution to Noise(delta) in the case of Lévy processes, just like Brownian motion.

Now, the formula for strong order 1.0 including jumps that I found in Platen’s book at section 6.2 “Strong Order 1.0 Taylor Scheme” (formula 6.2.1 up to 6.2.3) assumes we are using a compound Poisson process description. In practice, this might be undesirable, because some processes have infinite jump activity around jumps of size zero. Once again, the alpha stable process is a good example of this: we cannot simulate it exactly as a compound Poisson process, but we can simulate its increments exactly.

Up until now, I have been using a common technique which involves decomposing the compound Poisson representation into a small jump component and a big jump component, to avoid the infinite jumps around zero. However, I have now been wanting to simulate SDE’s driven by other noise processes whose small/big jump decomposition is very problematic, and hence using the direct simulation method for the increments is desirable. Using a custom NoiseProcess is a perfect way to achieve this. However, can we safely use the Runge-Kutta Milstein schemes defined for SDEProblems with custom noise processes and not worry that we are applying something where it should not be applied? Or are we forced to use only the lower order convergence schemes based on the Euler one? And what about the methods that deal with stiffness (are they the implicit ones?)?

2

Well, I just went to take a look at https://github.com/SciML/StochasticDiffEq.jl/blob/64a581c1dd268f4031caebdd07066360460dfeac/src/perform_step/low_order.jl

And found an example for one of the versions of the RKMil(stein) algorithm at

@muladd function perform_step!(integrator,cache::RKMilConstantCache)
  @unpack t,dt,uprev,u,W,p,f = integrator
  du1 = integrator.f(uprev,p,t)
  K = @.. uprev + dt * du1
  L = integrator.g(uprev,p,t)
  mil_correction = zero(u)
  if alg_interpretation(integrator.alg) == :Ito
    utilde =  K + L*integrator.sqdt
    ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
    mil_correction = ggprime.*(W.dW.^2 .- abs(dt))./2
  elseif alg_interpretation(integrator.alg) == :Stratonovich
    utilde = uprev + L*integrator.sqdt
    ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
    mil_correction = ggprime.*(W.dW.^2)./2
  else
    error("Alg interpretation invalid. Use either :Ito or :Stratonovich")
  end
  u = K+L.*W.dW+mil_correction

  if integrator.opts.adaptive
    du2 = integrator.f(K,p,t+dt)
    Ed = dt*(du2 - du1)/2
    En = W.dW.^3 .* ((du2-L)/(integrator.sqdt)).^2 / 6

    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

The assumption that we are using Brownian motion seems to be present at

ggprime.*(W.dW.^2 .- abs(dt))./2

Which I believe is the correction term
g(t_n,X)g'(t_n,X) \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} dW(t) dW(s) = g(t_n,X)g'(t_n,X) \frac {[W(t_{n+1}) - W(t_{n})]^2 - (t_{n+1} - t_n)}{2}

Which won’t apply for a pure jump process (the case I’m discussing here), where that term would be
g(t_n,X)g'(t_n,X) \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} dJ(t) dJ(s)
And I still haven’t found how to calculate this iterated integral, anyway.

TLDR: It seems we can’t use the higher-order schemes with a custom NoiseProcces because the calculations of iterated integral terms assume we are using Brownian motion.

3

They are safe, but indeed you get order loss. It’s not hard to show that the other schemes will converge at the rate of the Holder continuity which is what Euler-Maryama will give on arbitrary noise processes. So it is safe, but not higher order so you should just use the Euler methods when going to some of the other cases. That’s a bit hard to test for automatically since things like NoiseGrid are commonly used with Brownian processes though.

Indeed. That could in theory get generalized but I don’t think anyone has even written a paper on that. The most common thing to do is just rewrite the noise process as an SDE. Things like pink noise and such all have SDE forms which can be used and then you get a formulation for which higher order is easy again.

4

It’s not hard to show that the other schemes will converge at the rate of the Holder continuity which is what Euler-Maryama will give on arbitrary noise processes.

That’s pretty interesting, and unexpected to me. Since some of these processes have no continuous (stochastic) component, I assumed the schemes would actually have higher error rate (and lower convergence) than the EM one, due to wrong formulae.

Indeed. That could in theory get generalized but I don’t think anyone has even written a paper on that. The most common thing to do is just rewrite the noise process as an SDE. Things like pink noise and such all have SDE forms which can be used and then you get a formulation for which higher order is easy again.

I have also done a small literature review on the subject and couldn’t find anything as well. There seem to have been plenty of people deriving product formulae for iterated integrals of Lévy processes, so it’s probably something that could be computed for specific cases where you already know the Lévy measure of the process you want to integrate on, albeit I imagine the computations would be a nightmare to do.
I honestly felt tempted to give it a try to derive the second order correction terms for alpha-stable processes and one of their tempered versions (in 1 dim), since it would be just a computing a single iterated integral and a quadratic covariation term, but I was intimidated by the technicality of these papers that derive the product formulae for iterated integrals.

5

It would come from the results on Random Ordinary Differential Equations (RODEs), which is probably a better formulation for arbitrary noise processes.

Indeed it would be nice to generalize the interface for this, but then we’d also have to work out the math for every case too so it’s not too much point right now. But it’s really not to hard to setup if the literature begins to exist :wink: .