Piecewise Deterministic Markov Processes and Jump Diffusion Equations

Jump Diffusion equations are stochastic differential equations with discontinuous jumps. These can be written as:

\[du = f(u,p,t)dt + \sum_{j}g_j(u,p,t)dW_j(t) + \sum_{i}h_i(u,p,t)dN_i(t)\]

where $N_i$ is a Poisson-counter which denotes jumps of size $h_i$. In this tutorial we will show how to solve problems with even more general jumps. In the special case that $g_j = 0$ for all $j$, we'll call the resulting jump-ODE a piecewise deterministic Markov process.

Before running this tutorial please install the following packages if they are not already installed

using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")

DifferentialEquations.jl will install DiffEqJump, along with the needed ODE and SDE solvers.

We then load these packages, and set some plotting defaults, as

using DifferentialEquations, Plots
default(; lw = 2)

Defining a ConstantRateJump Problem

To start, let's solve an ODE that is coupled to a ConstantRateJump. A jump is defined as being "constant rate" if the rate is only dependent on values from other ConstantRateJumps or MassActionJumps (a special type of ConstantRateJump). This means that its rate must not be coupled with time, the solution to the differential equation, or a solution component that is changed by a VariableRateJump. ConstantRateJumps are cheaper to compute than VariableRateJumps, and so should be preferred when mathematically appropriate.

(Note: if your rate is only "slightly" dependent on the solution of the differential equation, then it may be okay to use a ConstantRateJump. Accuracy loss will be related to the percentage that the rate changes over the jump intervals.)

Let's solve the following problem. We will have a linear ODE with a Poisson counter of rate 2 (which is the mean and variance), where at each jump the current solution will be halved. To solve this problem, we first define the ODEProblem:

function f(du,u,p,t)
    du[1] = u[1]
    nothing
end

prob = ODEProblem(f, [0.2], (0.0, 10.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 0.2

Notice that, even though our equation is scalar, we define it using the in-place array form. Variable rate jump equations will require this form. Note that for this tutorial we solve a one-dimensional problem, but the same syntax applies for solving a system of differential equations with multiple jumps.

Now we define our rate equation for our jump. Since it's just the constant value 2, we do:

rate(u, p, t) = 2

rate (generic function with 1 method)

Now we define the affect! of the jump. This is the same as an affect! from a DiscreteCallback, and thus acts directly on the integrator. Therefore, to make it halve the current value of u, we do:

function affect!(integrator)
    integrator.u[1] = integrator.u[1] / 2
    nothing
end

affect! (generic function with 1 method)

Then we build our jump:

jump = ConstantRateJump(rate, affect!)

ConstantRateJump{typeof(Main.rate), typeof(Main.affect!)}(Main.rate, Main.affect!)

Next, we extend our ODEProblem to a JumpProblem by attaching the jump:

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

JumpProblem with problem ODEProblem with aggregator Direct
Number of constant rate jumps: 1
Number of variable rate jumps: 0
Number of mass action jumps: 0

We can now solve this extended problem using any ODE solver:

sol = solve(jump_prob, Tsit5())
plot(sol)

Variable Rate Jumps

Now let's define a jump with a rate that is dependent on the differential equation via the solution vector. Let's set the rate to be the current value of the solution, that is:

rate(u,p,t) = u[1]

rate (generic function with 1 method)

Using the same affect! we build a VariableRateJump:

jump = VariableRateJump(rate, affect!)

VariableRateJump{typeof(Main.rate), typeof(Main.affect!), Nothing, Float64, Int64}(Main.rate, Main.affect!, nothing, true, 10, (true, true), 1.0e-12, 0)

To make things interesting, let's copy this jump:

jump2 = deepcopy(jump)

VariableRateJump{typeof(Main.rate), typeof(Main.affect!), Nothing, Float64, Int64}(Main.rate, Main.affect!, nothing, true, 10, (true, true), 1.0e-12, 0)

so that way we have two independent jump processes. We now couple these jumps to the ODEProblem:

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

JumpProblem with problem ODEProblem with aggregator Direct
Number of constant rate jumps: 0
Number of variable rate jumps: 2
Number of mass action jumps: 0

which we once again solve using an ODE solver:

sol = solve(jump_prob, Tsit5())
plot(sol)

In this way we have solve a mixed jump-ODE, i.e. a piecewise deterministic Markov process.

Jump Diffusion

Now we will finally solve the jump diffusion problem. The steps are the same as before, except we now start with a SDEProblem instead of an ODEProblem. Using the same drift function f as before, we add multiplicative noise via:

function g(du, u, p, t)
  du[1] = u[1]
  nothing
end

prob = SDEProblem(f, g, [0.2], (0.0, 10.0))

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 0.2

and couple it to the jumps:

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

JumpProblem with problem SDEProblem with aggregator Direct
Number of constant rate jumps: 0
Number of variable rate jumps: 2
Number of mass action jumps: 0

We then solve it using an SDE algorithm:

sol = solve(jump_prob, SRIW1())
plot(sol)