2
$\begingroup$

A jump-diffusion process is a stochastic process where both continuous noise (in my case complex Wiener noise $dZ,dZ^*$ such that $dZ^2=dZ^{*2}=0,|dZ|^2=dt$) and discrete Jumps (in my case Poissonian $dN,dM$ such that $dN^2=dN$, $dM^2=dM$, $\langle dN\rangle=A dt$, $\langle dM \rangle=B dt$ ). My equations are in the Itô sense.

Using an update rule with a deterministic evolution term $dt$ as well as the noises, I can evolve a complex vector (or ket-if you like) $v$ through time. Coeficients $A$ and $B$ depend on $v$, but not so much that a $dN>1$ would cause problems.

I am currently solving this system with the plain forward Euler method in MATLAB, generating Wiener nose with randn and Poisson-noise with poissrnd. In addition, after each timestep I do an additional renormalization $v\rightarrow v/||v||$ (I know the norm is a constant of motion).

However, in order to make everything converge properly, I have to set the timestep $dt$ so small that the solving of the equations takes a long time. Are there good practical solvers (e.g. Runge-Kutta) with more performance that Euler?

When I try to look around what is commonly used for jump-diffusion processes, e.g in DifferentialEquations.jl or QuantumOptics.jl, the expected time between two jumps is typically much larger than the timestep( $dN,dM\approx 0$) so that one can distinguish deterministic evolution from jumps taking place at random times. This assumption is not fulfilled in my case, the jumprate is so high (it dominates the dynamics) that jumps may coincide the same timestep. On the other hand, it is not so macroscopically large that we can just use the central-limit-theorem and approximate the poisson-terms by gaussian noise.

Improve this question
$\endgroup$
2
  • $\begingroup$ Is it essential that your renormalisation happens at every time step or is just to avoid numerical overflows or similar? $\endgroup$
    – Wrzlprmft
    Commented Mar 8, 2018 at 16:53
  • $\begingroup$ @Wrzlprmft the second case- there is also a simpeler version of the equation for which the renormalisation ís essential, possible advantage of the latter is that depends only on dZ and not on dZ*, but I expect it 'd be harder for the latter to have an higher-order method. $\endgroup$
    – Wouter
    Commented Mar 8, 2018 at 17:43

1 Answer 1

Reset to default
2
$\begingroup$

The setup that is used in DifferentialEquations.jl and QuantumOptics.jl is what's known as time-adaptive jumping. It's nice because it allows for jump events to do things like change the number of DEs, and the jumps are computed exactly. However, it does have the limitation that jumps are computed exactly, so if you have a high jump rate then this slows down.

So then the other option is to use regular jumps instead of time-adapted jumps. The easiest of course is the Euler method, also known as tau-leaping in systems biology circles. The best resource on this topic for jump diffusions is probably Numerical Solution of Stochastic Differential Equations with Jumps in Finance. This book includes many high order schemes for regular jumping, but they do require lots of derivatives except for a few Runge-Kutta type schemes. There are many schemes for pure-jump problems (i.e. where the drift and diffusion terms are zero) that are not mentioned in that book but only in systems biology literature (such as binomial leaping), but of course if you always have drift and diffusion terms that doesn't matter. This is still a very active area of research.

As for implementations that are out there, I do not know of any that are not for pure-jump problems (and for pure-jump problems the ones I know of are for biological problems like StochSS). FWIW DifferentialEquations.jl is restructuring to allow for implementing the regular jump integrators (i.e. going through and implementing most of Platen's, along with using SciCompDSL.jl to automatically build the Jacobians for the user), but I wouldn't expect that for a few months.

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks! Now, I start to realize that the regime where the numerical problems occur that forced me to make the timestep so small may be the same regime where the jumprate is the lowest, so that the time-adaptive jumping in that regime (and euler as before with larger timestep elsewhere) could be a way out after all. I see on the docs of DifferentialEquations.jl that complex numbers are compatible with the SDE-solvers. Does this include problems that depend both on dZ and dZ* or only the ones that only depend on dZ and not dZ*? $\endgroup$
    – Wouter
    Commented Mar 8, 2018 at 18:14
  • $\begingroup$ The default use of complex numbers matches what you're describing. The real and imaginary parts of the number in the du vector are for the respective real and complex changes. The complex random numbers are drawn with the variance you describe. $\endgroup$
    – Chris Rackauckas
    Commented Mar 8, 2018 at 18:25
  • $\begingroup$ Thanks, will accept this answer tomorrow if it all seems to work out. $\endgroup$
    – Wouter
    Commented Mar 8, 2018 at 18:46
  • $\begingroup$ I've been reading a bit in the link from finance you posted and it seems very interesting, especially the part on strong jump-adapted schemes. But I'm still a bit struggling with the transition from a real to a complex process. If I write e.g. dZ=:(dW_x+idW_p)/sqrt(2) with real dW_x,dW_p, there remains a difference between adding an complex prefactor before the dW_x or putting the imaginary part before the dW_p. I see that in principle one could redefine a general complex noise dW_gen at each time such that dW_gen is absent, but then dW_gen^2 will be neither dt or zero but time-dependent? $\endgroup$
    – Wouter
    Commented Mar 9, 2018 at 10:44
  • $\begingroup$ Also in e.g. the SRIW1 algorithm you have this Milstein-like 'dW²-dt' part, this looks like it should only work for real dW? $\endgroup$
    – Wouter
    Commented Mar 9, 2018 at 10:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.