SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
141 stars 35 forks source link

Initial time points with jumps are double saved #58

Closed rcalxrc08 closed 3 years ago

rcalxrc08 commented 6 years ago

I am solving a Jump Diffusion SDE with the following code:

using DifferentialEquations
Nsim=10000;
Nsteps=100;
#Starting Point
r=0.02
sigma=0.2
T=1.0;
d=0.01;
u0=0.0;
#Drift part
f(u,p,t) = (r-d-sigma^2/2.0)
#Diffusion part
g(u,p,t) = sigma
#Time step
Dt =T/Nsteps;
#Time Window
tspan = (0.0,T)
#Definition of the SDE
prob = SDEProblem(f,g,u0,tspan)
#Definition of the Jump
rate(u,p,t) = 2
affect!(integrator) = (integrator.u = integrator.u/2)
jump = ConstantRateJump(rate,affect!)
#Definition of the Jump Problem
jump_prob = JumpProblem(prob,Direct(),jump)
#Definition of the MonteCarlo Problem
monte_prob = MonteCarloProblem(jump_prob)
sol = solve(monte_prob,SOSRI(),num_monte=Nsim,parallel_type=:threads,dt=Dt,adaptive=false)

When I look into the elements of the paths of the solution I get:

lengthT=length(sol[1].t) # =109 != Nsteps+1 and random
lengthUniqueT=length(unique(sol[1].t)) # =103 != Nsteps+1 and random

The problem is not only the repetitions of the elements (for example the initial time is repeated twice) but also the time discretization step that is not uniform.

The problem shows up also if I don't put the flag on the parallel type parallel_type=:threads, but I don't know which is the default argument for such flag.

ChrisRackauckas commented 5 years ago

Wow, I think I missed this one. Sorry!

There are duplicate values for the left and right continuous values at the point of a jump discontinuity. So each jump time will be doubled. Since the jumps are at random times, that leads to the non-uniform times since dt will decrease to hit the jump times exactly.

I do see a slight issue here though. The first time shouldn't be doubled. This is because the callback initialization for the jump is hitting the code for the right-left discontinuity handling in saving, and then causing a second save. That should get fixed.

isaacsas commented 3 years ago

This code no longer runs, but the doubling of the initial value should be closed by https://github.com/SciML/DiffEqJump.jl/pull/165

As Chris said, the nonuniform times are because the jump times are saved too. We can reopen if there are still issues when using EnsembleProblems.