# Jump Diffusion Equations

This tutorial assumes you have read the Ordinary Differential Equations tutorial.

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.

## Defining a ConstantRateJump Problem

To start, let's solve an ODE with constant rate jumps. A jump is defined as being "constant rate" if the rate is only dependent on values from other constant rate jumps, meaning that its rate must not be coupled with time or the solution to the differential equation. However, these types of jumps are cheaper to compute.

(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]
end
prob = ODEProblem(f,[0.2],(0.0,10.0))
```

Notice that, even though our equation is on 1 number, 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`

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:

`affect!(integrator) = (integrator.u[1] = integrator.u[1]/2)`

Then we build our jump:

`jump = ConstantRateJump(rate,affect!)`

Next, we extend our `ODEProblem`

to a `JumpProblem`

by attaching the jump:

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

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 which is coupled to the differential equation. Let's let the rate be the current value of the solution, that is:

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

Using the same `affect!`

`affect!(integrator) = (integrator.u[1] = integrator.u[1]/2)`

we build a `VariableRateJump`

:

`jump = VariableRateJump(rate,affect!)`

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

`jump2 = deepcopy(jump)`

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

:

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

which we once again solve using an ODE solver:

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

## 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]
end
prob = SDEProblem(f,g,[0.2],(0.0,10.0))
```

and couple it to the jumps:

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

We then solve it using an SDE algorithm:

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

## Coupling Jump Problems

In many applications one is interested in coupling two stochastic processes. This has applications in Monte Carlo simulations and sensitivity analysis, for example. Currently, the coupling that is implemented for jump processes is known as the split coupling. The split coupling couples two jump processes by coupling the underlying Poisson processes driving the jump components.

Suppose `prob`

and `prob_control`

are two problems we wish to couple. Then the coupled problem is obtained by

`prob_coupled = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map)`

Here, `coupling_map`

specifies which jumps to couple. If `(j,i)`

is in `coupling_map`

, then the `i`

th jump in `prob`

will be coupled to the `j`

th jump in `prob_control`

. Note that currently `SplitCoupledJumpProblem`

is only implemented for constant rate jump problems.

As an example, consider a doubly stochastic Poisson process, that is, a Poisson process whose rate is itself a stochastic process. In particular, we will take the rate to randomly switch between zero and `10`

at unit rates:

```
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(u0,tspan)
jump_prob = JumpProblem(prob,Direct(),jump1,jump2,jump3)
```

The doubly stochastic poisson process has two sources of randomness: one due to the Poisson process, and another due to random evolution of the rate. This is typical of many multiscale stochastic processes appearing in applications, and it is often useful to compare such a process to one obtained by removing one source of randomness. In present context, this means looking at an ODE with constant jump rates, where the deterministic evolution between jumps is given by the expected value of the Poisson process:

```
function f(du,u,p,t)
du[1] = u[2]*10
du[2] = 0.
du[3] = 0.
end
prob_control = ODEProblem(f,u0,tspan)
jump_prob_control = JumpProblem(prob_control,Direct(),jump2,jump3)
```

Let's couple the two problems by coupling the jumps corresponding the switching of the rate:

```
coupling_map = [(2,1),(3,2)]
prob_coupled = SplitCoupledJumpProblem(jump_prob,jump_prob_control,Direct(),coupling_map)
```

Now `prob_coupled`

will be dealt with like any other `JumpProblem`

:

`sol = solve(prob_coupled,Tsit5())`