Didn’t create an issue on the repo because this might be just my interpretation being wrong.
I am trying to run code based on the tutorial at Jump Diffusion page.
It mentions that this is supposed to be a Cox process (doubly stochastic) whose rate switches between 10 and 0, but I’m having trouble understanding how this code is supposed to achieve this.
The system that seems to be implied is something like
x(t) = N_a(y, t)
y(t) = N_c(t) - N_b(t)
z(t) = N_b(t) - N_c(t)
Where
N_a (y, t) ~ Poisson(10yt)
N_b(y, t) ~ Poisson(yt)
N_c(z, t) ~ Poisson(zt)
First of all, the code fails if I try to solve this system with almost any u0 except (1,1,1) or (1,0,0) with the error message
“Tried to add a tstop that is behind the current time. This is strictly forbidden”
And here is what I get if I initialize to (1,1,1):
The code used, below:
rate(u,p,t) = u[2]*10
affect!(integrator) = integrator.u[1] += 1.
jump1 = ConstantRateJump(rate,affect!)
rate(u,p,t) = u[2]
affect!(integrator) = (integrator.u[2] -= 1.;integrator.u[3] += 1.)
jump2 = ConstantRateJump(rate,affect!)
rate(u,p,t) = u[3]
affect!(integrator) = (integrator.u[2] += 1.;integrator.u[3] -= 1.)
jump3 = ConstantRateJump(rate,affect!)
prob = DiscreteProblem(ones(1),(0.0,10.0))
jump_prob = JumpProblem(prob,Direct(),jump1,jump2,jump3)
plot(solve(jump_prob))