Simulating a stopping time in SDE without using a fixed time span?

Hi everyone,

I am trying to simulate a stochastic process on [0,1] that will hit the boundary with probability 1 in a finite amount of time. The package I am using is the SDE solver in DifferentialEquation.jl. However, all the recipes I found need a fixed time span for the solver. i.e. I can only specify a common time frame, say, 0\le t\le T for all the sample paths in the ensemble. But I want the solver to stop after hitting the boundary, or continue to solve for longer time if the path has not hit the boundary before T, is it possible in DifferentialEquations.jl?

My cumbersome solution now is to do a while-loop for each solving attempt with a fixed time frame, and break out if at the end of each iteration the process has hit the boundary, otherwise continue to solve using the end state of the process at the last iteration. But I think it is not quite efficient.

Thanks!

https://diffeq.sciml.ai/stable/features/callback_functions/#Example-2:-Terminating-an-Integration

4 Likes

I thought I mention here that if someone is interested in the inverse problem - that is observing the hitting time of the boundary and then trying to say something about the underlying process - please ping me.

1 Like

Thanks Chris, it is exactly what I was looking for!

Hi Chris,

I tried this Callback with SDE and it works great, but when I used in in Discrete Problems, there seems to be some instability in the performance. For instance, I try the following toy example:

f(u,p,t) = u+1

prob = DiscreteProblem(f,3,(0.0,20.0))

rate(u,p,t) = 1
affect!(integrator) = (integrator.u = integrator.u + 1)
jump = ConstantRateJump(rate,affect!)

condition(u,t,integrator) = (u >= 10)
function affect_stop!(integrator)
    println("hit!") 
    terminate!(integrator)
end
cb = DiscreteCallback(condition,affect_stop!)

jump_prob = JumpProblem(prob,Direct(),jump)

sol = solve(jump_prob,callback=cb)

The solution is supposed to stop when it reached 10 and print out “hit!”

But sometimes I get solutions like this:

hit!
hit!
hit!
hit!
retcode: Terminated
Interpolation: left-endpoint piecewise constant
t: 22-element Array{Float64,1}:
 0.0
 1.0
 1.2177921614384977
 1.2177921614384977
 1.3483018817001802
 1.3483018817001802
 1.9142430896065223
 1.9142430896065223
 2.2600182450752926
 2.2600182450752926
 2.2600182450752926
 2.2600182450752926
 2.3244255770166746
 2.3244255770166746
 2.3244255770166746
 2.3244255770166746
 2.9412940817863067
 2.9412940817863067
 2.9412940817863067
 2.9412940817863067
 3.9412940817863067
 3.9412940817863067
u: 22-element Array{Int64,1}:
  3
  4
  5
  6
  7
  8
  9
 10
 11
 11
 11
 12
 13
 13
 13
 14
 15
 15
 15
 16
 17
 17

The jump process went on even after it was terminated. Then I thought maybe it is because I added a continuous-time jump process to the iteration, and I should use the ContinuousCallback() function, but switching to that also doesn’t work.

Thanks for your feedback!

I just discovered that if I add

SSAStepper()

In the solve() function, then the result can be quite accurately attained. No clue why it works.

I turned this into an issue and explained why it happens in https://github.com/SciML/DifferentialEquations.jl/issues/691. Track that issue for a solution down the road.