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) 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?
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)
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?
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).
Excellent, I’m glad to hear that. The interface to all of these methods is remarkably natural and easy to learn. Thanks!
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.
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
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.
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:
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
condition(t,u,integrator)
rand(Poisson(rate(integrator)) > 0
end
for readability, or
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 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!
2 Likes
I’m solving more or less the same problem in 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. 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 in integrate_mcwf() and jump(). Maybe this code can be of some use…
No. If I am interpreting it correctly, he’s asking for a 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 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.
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 which was written to state the inaccuracies of this paper. 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.
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.
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.
Thanks for the clarification - I got confused by the name similarity to the 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.
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.
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.
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.
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?)
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!
No in the interval, when the norm is equal to the random number.
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!
Yes exactly. I guess this works only for a very specialized set of problems.
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:
I’ll have some easy to use implementations of those jump callbacks together in a few days and start to get everything tagged.