# Code for simulating jump diffusion?
**URL:** https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154
**Category:** General Usage
**Tags:** diffeq
**Created:** [December 26, 2016, 5:14pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154 "2016-12-26T17:14:58Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [December 26, 2016, 5:14pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/1 "2016-12-26T17:14:58Z")
I would like to simulate from a jump diffusion model. I have written C++ code for use with Octave ([NeuralNetsForIndirectInference.jl/CTSVmodel.cc at master · mcreel/NeuralNetsForIndirectInference.jl · GitHub](https://github.com/mcreel/NeuralNetsForIndirectInference.jl/blob/master/CTSV/Code/Common/CTSVmodel.cc)) and I could use that as the basis for Julia code, but I’m hoping that someone has already done this. Any resources out there?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 26, 2016, 7:14pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/2 "2016-12-26T19:14:08Z")
When the events and callbacks are added to StochasticDiffEq.jl (DifferentialEquations.jl) I’ll write an example for jump diffusions. That shouldn’t be more than 2 or so weeks off (I need the callbacks done by then for other reasons, but once they’re done it’s a small step to get jump diffusions)
---
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [January 25, 2017, 8:03am UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/3 "2017-01-25T08:03:05Z")
I have been looking at DifferentialEquations.jl, impressive work! I understand that a callback could be used to insert jumps into the solution. Would it be possible that the probability that a jump occurs in a given interval after t depend on the values of the forcing or noise functions at t? Or would jumps have to be exogenous?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 8:11am UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/4 "2017-01-25T08:11:06Z")
> [@mcreel](#):
>
> Would it be possible that the probability that a jump occurs in a given interval after t depend on the values of the forcing or noise functions at t? Or would jumps have to be exogenous?
It will be possible once the integrator+callbacks are available for the SDE solvers. I’m actually working on that right now and will probably have it done by tomorrow. I’ll post with an example when it’s ready on master (and tag it soon after).
---
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [January 25, 2017, 8:31am UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/5 "2017-01-25T08:31:02Z")
Excellent, I’m glad to hear that. The interface to all of these methods is remarkably natural and easy to learn. Thanks!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 2:53pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/6 "2017-01-25T14:53:10Z")
Alright, so the adding of the integrator + callback interfaces to StochasticDiffEq.jl has finished its first draft. I won’t be releasing it quite yet since I’ll be adding a bunch of unit tests and making sure the edges are all robust + profiling and optimizing. However, I suspect that from where I am at I’ll release either tonight or tomorrow (I’m going full-time on this until it’s done).
A little bit of discussion on the math here. Note that I will be adding callbacks for making this more user friendly (should be able to create these callbacks just by giving a rate function), but this discussion will let you know what can currently be done without error and what has what amount of discretization error. Given the method from your link, it looks like you’re fine with having a good amount of discretization error due to interval size, but I thought I’d let you know where I am trying to fix that.
> [@mcreel](#):
>
> Would it be possible that the probability that a jump occurs in a given interval after t depend on the values of the forcing or noise functions at t? Or would jumps have to be exogenous?
Both are possible. Let’s split the discussion since there’s interesting things for both.
### Constant Rate Jumps
First of all, if the jumps have a constant rate (this includes the case where the rate is constant between jumps as well), we can use the `DiscreteCallback` framework to model it exactly. This can be done via what’s known in systems biology as Gillespie’s method. Let’s say the jump rate is `lambda`. Then, given the jumps are Poisson distributed, we know that the time to the next jump is exponentially distributed where the rate parameter there is also `lambda` (easy to prove). Thus given the previous jump was at a time `t`, the next jump’s distribution is given by a random number of the form
```julia
rand(Exponential(1/lambda))
```
(using `Distributions.jl` for the exponentially distributed random number, using `1/lambda` since it internally uses the scale instead of the rate parameter). So here’s how we’d solve this. We can overload the call on a type to make this our `DiscreteCallback` condition.
```julia
type ConstantJump{T}
next_jump::T
end
(p::ConstantJump{T})(t,u,integrator) = (p.next_jump == t)
```
Notice that this function is true when `t == next_jump`. Next we need the affect function. Our affect is to take a new jump, and add it to `tstops` as well. This will force the solver to hit the timestep of `next_jump` exactly, and cause the next condition. We do this by the overload:
```julia-auto
function (p::ConstantJump{T}(integrator)
p.next_jump = integrator.t + rand(Exponential(1/lambda))
add_tstop!(integrator,p.next_jump)
u_modified!(integrator,false) # Stops re-evaluation checks
end
```
We’d start by taking a time from `tspan[1]+rand(Exponential(1/lambda))` which would be the for the first jump. That would be part of the constructor for the callback, and we’d also have to start with this value in `tstops`.
This format makes us hit every time for the constant rate jumps exactly, so there’s no discretization error for the jumps. If the jump rate is small, we’d want to replace this with a tau-leaping method for more efficiency (handles multiple jumps at a time, but we’d have to modify it to account for the continuous updates, so that’s a story for another day).
### Non-Constant Rate Jumps : Naive Approach
For this case we cannot determine values for the jump in advance, so instead we use the Poisson distribution directly on the interval. Let’s start with a naive approach. From the looks of your code, I think this is similar to what you had. The probability that a Poisson random variable has a count in an interval of length `dt` with rate `lambda` is `(1-exp(-dt*lambda)`. Here, we can evaluate lambda at either endpoint, or we can do some sophisticated interpolation (use the average of the endpoints), So let `lambda = rate(integrator)` where `rate(integrator)` can use `integrator(t)` the interpolation of the value throughout the whole current interval to determine some value for lambda.
Given that we can calculate this, we can take a random number to determine if an event occurred in the interval. Once again, we can use a `DiscreteCallback` to implement this, where now we have the condition that
```julia
condition(t,u,integrator)
rand(Poisson(rate(integrator)) > 0
end
```
for readability, or
```julia
condition(t,u,integrator)
rand() > exp(-integrator.dt*rate(integrator))
end
```
should be slightly more efficient but have the same distribution. Now, since we aren’t constraining the timestep, we can define any `affect!(integrator)` function and be on our merry way.
This will work, but has discretization error since it does not get the timepoint of the jump exactly, and instead will apply it at `t+dt`. I believe this would converge with 1/2-order accuracy.
### Second Non-Constant Jump Rate Approach: Wong-Zakai
Instead of using a discrete callback, we can use a `ContinuousCallback`. Here we take a uniform random number to determine the timepoint in the interval the event occured, and can use a linear interpolation and backtrack to that timepoint. This should also be 1st-order accurate. There’s a caveat here: to get strong-order accuracy, we have to smartly handle how we change the Brownian motion. [The adaptive timestepping scheme can account for this](http://chrisrackauckas.com/assets/Papers/ChrisRackauckas-AdaptiveSRK.pdf) so we’re good, though this will require some setup and testing (which is what I am doing today).
### Higher Order Accuracy
To get higher order accuracy, we need to be even smarter. What we need to do solve a first passage time problem and get an approximate distribution for both the probability that an event occurred and a distribution for the time of the event. I haven’t found any literature on this, so this may be one of my next research projects.
### tl;dr
Using `DiscreteCallbacks` will do this easy. It can do constant time jumps exactly. It can do non-constant rates decently, but using a linear interpolation and the internal stacks from the adaptive method we can do even better. However, the holy grail requires solving a first-passage problem. I’ll have functions which build these callbacks as part of the Callback Library to make this all easier, but this is all for transparency about the methods.
Hope this helps!
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 4:48pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/7 "2017-01-25T16:48:41Z")
I’m solving more or less the same problem in [QuantumOptics.jl](https://github.com/bastikr/QuantumOptics.jl). As ode solver I use an adaptive Runge-Kutta Dormand-Prince 4(5) solver with event handling and dense output that can be found at [https://github.com/bastikr/QuantumOptics.jl/blob/master/src/ode\_dopri.jl](https://github.com/bastikr/QuantumOptics.jl/blob/master/src/ode_dopri.jl). The function of interest is in this case the `ode_event()`. An example of how it can be used to solve a stochastic equation can be found in [https://github.com/bastikr/QuantumOptics.jl/blob/master/src/mcwf.jl](https://github.com/bastikr/QuantumOptics.jl/blob/master/src/mcwf.jl) in `integrate_mcwf()` and `jump()`. Maybe this code can be of some use…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 5:24pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/8 "2017-01-25T17:24:31Z")
> [@bastikr](#):
>
> As ode solver I use an adaptive Runge-Kutta Dormand-Prince 4(5) solver with event handling and dense output that can be found at [https://github.com/bastikr/QuantumOptics.jl/blob/master/src/ode\_dopri.jl](https://github.com/bastikr/QuantumOptics.jl/blob/master/src/ode_dopri.jl). The function of interest is in this case the ode\_event().
No. If I am interpreting it correctly, he’s asking for a [jump diffusion](https://en.wikipedia.org/wiki/Jump_diffusion), meaning a stochastic differential equation with jumps. It’s a very bad idea to use a deterministic method in this case. From Kloden Platen Schwartz / Rossler we know that this method will converge with at most order 1/2, but even then I am not sure this method will converge at all. This is a bad idea and will be very inaccurate.
In any case, if it’s just a deterministic equation with jumps, then [the DifferentialEquations.jl callback interface](http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-1:-Bouncing-Ball-1) will do what he wants (it’ll actually handle your case too). Even then, you likely want to use a more efficient method than dopri, like [Tsit5](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/NonStiffODE/Linear%20Work-Precision%20Diagrams.ipynb).
> [@bastikr](#):
>
> An example of how it can be used to solve a stochastic equation can be found in [https://github.com/bastikr/QuantumOptics.jl/blob/master/src/mcwf.jl](https://github.com/bastikr/QuantumOptics.jl/blob/master/src/mcwf.jl) in integrate\_mcwf() and jump(). Maybe this code can be of some use…
Your approach is the above “Second Non-Constant Jump Rate Approach: Wong-Zakai”, but you use the full interpolation from the dopri method which again, is unjustified if the equation is an SDE and will result in significant order loss. Wong-Zakai’s proof shows you can actually do this at low accuracy with a linear interpolation, but again a first-passage problem has to be approximated for higher accuracy.
Just to clarify, I don’t mean to be mean here, but I just want to clarify that deterministic methods do not achieve appropriate accuracy on stochastic differential equations like jump diffusions. Many times they will not strong converge, and their weak convergence will be extremely low order. One paper which states this very nicely is [this rebuttal](http://personal.strath.ac.uk/d.j.higham/papers/pre_nsde.pdf) which was written to state the inaccuracies of [this paper](https://arxiv.org/pdf/quant-ph/0407039.pdf). If you don’t properly address the subtleties, then the best method is just to use Euler-Maruyama. But in Julia we already have freely available and well-tested high strong order strongly adaptive methods, so we might as well build tools off of them as I demonstrated above.
---
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [January 25, 2017, 5:53pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/9 "2017-01-25T17:53:06Z")
Thanks for the interesting discussion. To be clear, I am interested in SDEs with jumps, as in your detailed message above. For my purposes, simulation based parameter estimation, speed is more important than great accuracy, but having options to go either way would be great. My C++ code for Octave is fairly quick, I believe, but it may have poor accuracy, I am not very knowledgeable about these methods.
I will be happy to try the new code out in the near future, when I get a few hours free.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 6:09pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/10 "2017-01-25T18:09:53Z")
> [@mcreel](#):
>
> For my purposes, simulation based parameter estimation, speed is more important than great accuracy, but having options to go either way would be great. My C++ code for Octave is fairly quick, I believe, but it may have poor accuracy, I am not very knowledgeable about these methods.
Then the Wong-Zukai approach should do well since the overhead over the naive approach should be minimal but it should allow the adaptive timestepping method to take as large steps as possible without worrying too much about errors due to large intervals. I’ll make sure to have an example illustrate the difference.
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 6:16pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/11 "2017-01-25T18:16:04Z")
Thanks for the clarification - I got confused by the name similarity to the [Quantum jump method](https://en.wikipedia.org/wiki/Quantum_jump_method). My knowledge of SDEs is very limited. As far as I’m aware this method can not be written as an SDE but I might be wrong about that. I will try to understand if these objections are also valid for this method. Thank you for taking the time to provide many interesting links! It would be definitely be advantageous to use specialized packages instead of using my own implementation. However, when I started writing it I could not find any acceptable alternatives.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 6:33pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/12 "2017-01-25T18:33:03Z")
> [@bastikr](#):
>
> As far as I’m aware this method can not be written as an SDE but I might be wrong about that. I will try to understand if these objections are also valid for this method.
Oh I see. You can write that as a jump diffusion, but the diffusion term would be trivial. The objections are not valid for the quantum jump method because there’s no continuous stochasticity, so as long as you use event handling to appropriately determine the discontinuities you’re fine. That seems to be what you’re doing.
Though may I ask, how are you determining the location in the interval for the jump? It seems the jumps are Poisson with rates are dependent on the wave function, so I am curious if there’s an algorithm determining the “right” random number (with some order of accuracy) for determining the location, since the approach I have is 1st order (uniform throughout the interval: same as the assumption that the rate is constant on the interval). This makes sense for SDEs since it’s a high enough accuracy, but for ODEs with jumps it would kill the (strong) accuracy.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 6:40pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/13 "2017-01-25T18:40:13Z")
> [@bastikr](#):
>
> It would be definitely be advantageous to use specialized packages instead of using my own implementation.
Your implementation is pretty solid, though I did notice you’re not using a PI-stepsize controller which is immensely helpful with Dormand-Prince methods (see Hairer Solving Ordinary Differential Equations II for details on that). DifferentialEquations.jl’s (well, technically OrdinaryDiffEq.jl’s) event handling and callback framework should be mature enough to handle implementing the quantum jump method (though a `ContinuousCallback`), and it has a dopri implementation (`DP5`, though `Tsit5` is more modern and usually more efficient) with PI-stepsize controls if you want to give it a try. Just open an issue if you want to discuss more.
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 6:50pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/14 "2017-01-25T18:50:23Z")
> [@ChrisRackauckas](#):
>
> Though may I ask, how are you determining the location in the interval for the jump? It seems the jumps are Poisson with rates are dependent on the wave function, so I am curious if there’s an algorithm determining the “right” random number (with some order of accuracy) for determining the location, since the approach I have is 1st order (uniform throughout the interval: same as the assumption that the rate is constant on the interval).
In this case the probability that NO jump has happened is just the norm of the state vector. So I choose a random number in the interval [0,1] (with a uniform distribution) and perform the jump when the norm equals this number. The statistical average of these trajectories is then the same as simulating the corresponding quantum master equation.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 6:58pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/15 "2017-01-25T18:58:02Z")
> [@bastikr](#):
>
> perform the jump when the norm equals this number
Yeah, but where in time do you perform the jump? At the end of the interval? (I assume you mean when the norm is greater than or equal to this number?)
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 6:58pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/16 "2017-01-25T18:58:26Z")
> [@ChrisRackauckas](#):
>
> Your implementation is pretty solid, though I did notice you’re not using a PI-stepsize controller which is immensely helpful with Dormand-Prince methods (see Hairer Solving Ordinary Differential Equations II for details on that).
For my implementation i used Hairer’s Solving Ordinary Differential Equations I, so I never encountered PI-stepsize controllers 😀.[quote=“ChrisRackauckas, post:13, topic:1154”]
DifferentialEquations.jl’s (well, technically OrdinaryDiffEq.jl’s) event handling and callback framework should be mature enough to handle implementing the quantum jump method (though a ContinuousCallback), and it has a dopri implementation (DP5, though Tsit5 is more modern and usually more efficient) with PI-stepsize controls if you want to give it a try. Just open an issue if you want to discuss more.
[/quote]
Thanks for the kind offer! I’m really tempted to switch to a better solver. If I find some time this weekend I will take a look!
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 6:59pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/17 "2017-01-25T18:59:46Z")
> [@ChrisRackauckas](#):
>
> Yeah, but where in time do you perform the jump? At the end of the interval? (I assume you mean when the norm is greater than or equal to this number?)
No in the interval, when the norm is equal to the random number.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2017, 7:10pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/18 "2017-01-25T19:10:02Z")
Oh, I see. Yeah, I wouldn’t think that the scheme is strongly convergent, though as you say it’s likely weakly (statistically) convergent (I wouldn’t know how to prove it but sounds like it’s something known in the field?). I was hoping to grab some intuition from it but sadly it seems very dependent on the fact that the norm is not preserved (it uses the monotonicity to make sure the jump occurs sometime in the future, right?).
Thanks for sharing, it’s a cool idea!
---
### Author:  [@bastikr](https://discourse.julialang.org/u/bastikr)
#### Post date: [January 25, 2017, 7:19pm UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/19 "2017-01-25T19:19:01Z")
> [@ChrisRackauckas](#):
>
> I was hoping to grab some intuition from it but sadly it seems very dependent on the fact that the norm is not preserved (it uses the monotonicity to make sure the jump occurs sometime in the future, right?).
Yes exactly. I guess this works only for a very specialized set of problems.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 30, 2017, 9:27am UTC](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154/20 "2017-01-30T09:27:28Z")
An update:
I have the full event handling framework working for SDEs on the `Roots` branch, so callbacks to implement jumps can be done using the methods I described. However, I am contemplating whether this release can really be considered v0.5-compatible, since it introduces a performance degredation due to a v0.5 bug:
> [@Handling #18465 on v0.5](https://discourse.julialang.org/t/handling-18465-on-v0-5/1773):
>
> [https://github.com/JuliaLang/julia/issues/18465](https://github.com/JuliaLang/julia/issues/18465) This issue causes some really weird performance regressions. I have a potential one related to using Roots.jl in a new release of the StochasticDiffEq.jl module for DifferentialEquations.jl that I would like to do. However, there doesn’t seem to be any Julia v0.5.x coming and no posted workaround, so I was wondering what I should do in terms of package release. I could just make a last v0.5 release pre-using Roots.jl, and then have a version wi…
I’ll have some easy to use implementations of those jump callbacks together in a few days and start to get everything tagged.
[Next page](https://discourse.julialang.org/t/code-for-simulating-jump-diffusion/1154.md?page=2)