Piecewise Deterministic Markov Processes and Jump Diffusion Equations

Jump Diffusion equations are stochastic differential equations (SDEs) 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)\]

The change in state vector $u$ depends on three concurrent dynamics, (1) deterministic drift $dt$ with size $f(u, p, t)$, (2) stochastic diffusions $dW_j(t)$ with sizes $g_j(u, p, t)$, and, (3) stochastic jumps $dN_i(t)$ with sizes $h_i(u, p, t)$. By diffusion we mean a continuous stochastic process which is usually represented as Gaussian white noise (i.e. $W_j(t)$ is a Brownian Motion). By jump we mean a discrete stochastic process which is usually represented by a Poisson process.

In this tutorial, we will show how to solve problems with jumps more general than the homogeneous Poisson process. 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 JumpProcesses, 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
Non-trivial mass matrix: false
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:

const_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(), const_jump)

JumpProblem with problem ODEProblem with aggregator Direct
Number of jumps with discrete aggregation: 1
Number of jumps with continuous aggregation: 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)

Note that every time a jump occurs the function peaks and then decreases by half following our specification.

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:

var_jump = VariableRateJump(rate, affect!)

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

We now couple the jump to the ODEProblem:

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

JumpProblem with problem ODEProblem with aggregator Direct
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 1
Number of mass action jumps: 0

which we once again solve using an ODE solver:

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

Again, when a jump occurs the function peaks and then decreases by half. In contrast to our previous example, as the function value decreases the rate of jumps also decreases.

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

Multiple jumps

In this section we allow for multiple jumps in the SDE. Let's assume that the dynamics of our SDE is affected by the two jumps we defined in the previous sections. It is easy to re-use all the elements we defined above and initialize a JumpProblem with multiple jumps.

jump_prob = JumpProblem(prob, Direct(), const_jump, var_jump)

JumpProblem with problem ODEProblem with aggregator Direct
Number of jumps with discrete aggregation: 1
Number of jumps with continuous aggregation: 1
Number of mass action jumps: 0

which we solve in the same manner using an ODE solver:

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

The constant jump ensures that the function jumps at more or less regular intervals, while the variable jump increases the frequency of jumps when the function reaches higher levels.

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
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 0.2

and couple it to the jumps:

jump_prob = JumpProblem(prob, Direct(), const_jump, var_jump)

JumpProblem with problem SDEProblem with aggregator Direct
Number of jumps with discrete aggregation: 1
Number of jumps with continuous aggregation: 1
Number of mass action jumps: 0

We then solve it using an SDE algorithm:

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

The function now behaves different than before due to the diffusion component which drifts randomly around zero. Again, note that jumps are more frequent when the function obtains a higher value and the function sharply decreases by half anytime a jump takes place.