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?)?