We want to solve the SDE with the following features:
- Start from an arbitrary initial condition of
N
particles
μ
is the drift of the Brownian motion
σ
is the variance
β
is the arrival rate of jumps
- Upon an arrival of a jump, the go to
max(z, z')
where z'
is drawn from the distribution.
To simulate, we will "solve" a jump-diffusion SDE. See http://docs.juliadiffeq.org/latest/tutorials/jump_diffusion.html
In the simplest version, we can solve this with constants for all of the parameters. Later we can have these a function of the state to simulate a "poor man's lucas-moll".
To get this working, see http://docs.juliadiffeq.org/latest/tutorials/jump_diffusion.html#Jump-Diffusion-1 for more details and http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1 for solving the SDE on its own.
But it would be something like the following to setup the SDE (without the jumps):
using DifferentialEquations, Distributions, StochasticDiffEq, DiffEqJump
function μ_SDE(du,u,p,t)
du .= p.μ
end
function σ_SDE(du,u,p,t)
du .= p.σ
end
p = (μ = 0.01, σ = 0.1, β = 0.1, N = 5) # if all constant
T = 10.0 # maximum time length
dist_iv = Normal() # initial condition distribution, though could hardcode
x_iv = rand(dist_iv, p.N) # just draws from the inital condition
prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, T), p)
Now, to turn this into a jump problem it will take a little more work figuring out the interface. The issue is that we want N
independent jump procesess. So maybe something like:
function affect!(integrator, n)
integrator.u[n] = max(integrator.u[n], integrator.u[rand(1:integrator.p.N)]) # choose random integrator form 1 to N
end
jumps = [ConstantRateJump(0.2, (integrator) -> affect!(integrator, n) for n in 1:p.N)] # add one for each n
jump_prob = JumpProblem(prob, Direct(), jumps) # might be wrong...
sol = solve(jump_prob)
plot(sol) # Can't plot if too many points, instead use ensemble.
The thing I am unclear of here is the jumps
setup, since we need N
of them. Should that be an array? Should it be a JumpSet
? Don't know, but ask for help on slack if having issues.
After you get the basics working, you should be able to crank up N
, use http://docs.juliadiffeq.org/latest/features/ensemble.html if you want to get fancy, etc. to get fancy, etc.