# Piecewise Deterministic Markov Processes and Jump Diffusion Equations

This tutorial assumes you have read the Ordinary Differential Equations tutorial in `DifferentialEquations.jl`

.

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 `ConstantRateJump`

s or `MassActionJump`

s (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 `VariableRateJump`

s, 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)
```