I’ve been using DifferentialEquations.jl to simulate the following stochastic volatility model:

du_t = 0 dt + v_t ( -0.576 dW_{1,t} + \sqrt{1 - 0.576^2} dW_{2,t})

dv_t^2 = 0.035 (0.636 - v_t^2) dt + 0 . 144 v_t^2 dW_{1,t}

The code I’m using to do this follows:

```
function f_mu!(du, u, p, t)
du[1] = 0.0
du[2] = 0.035 * (0.636 - u[2])
end
function f_sigma!(du, u, p, t)
vscaled = sqrt(u[2])
du[1,1] = -0.576 * vscaled
du[1,2] = sqrt(1 - 0.576^2) * vscaled
du[2,1] = 0.144 * u[2]
du[2,2] = 0.0
end
numperiod = 500
dt = 1 // 1000
prob = SDEProblem(f_mu!, f_sigma!, [log(30.0), 0.1], (0.0, Float64(numperiod)), noise_rate_prototype=fill(0.1, 2, 2))
sol = solve(prob, EulerHeun(), dt=dt)
```

This all works nicely. But now I want to add a jump component to both equations, ie:

du_t = 0 dt + v_t ( -0.576 dW_{1,t} + \sqrt{1 - 0.576^2} dW_{2,t}) + dJ_{1,t}

dv_t^2 = 0.035 (0.636 - v_t^2) dt + 0 . 144 v_t^2 dW_{1,t} + dJ_{2,t}

In both cases, the jump components can be expressed as `ConstantRateJump`

. I’ve read through the docs, and it seems to me that they are saying I can do this using the following code:

```
prob = SDEProblem(f_mu!, f_sigma!, [log(30.0), 0.0, 0.0], (0.0, Float64(numperiod)), noise_rate_prototype=fill(0.001, 2, 2));
f_rate_1 = ((u, p, t) -> 1.0)
f_rate_2 = ((u, p, t) -> 1.5)
f_affect_1! = ((integrator) -> (integrator.u[1] += (rand([-1, 1]) * 4.0)))
f_affect_2! = ((integrator) -> (integrator.u[2] += 1.6))
cjump_1 = ConstantRateJump(f_rate_1, f_affect_1!)
cjump_2 = ConstantRateJump(f_rate_2, f_affect_2!)
jump_prob = JumpProblem(prob, Direct(), cjump_1, cjump_2)
sol = solve(jump_prob, EulerHeun(), dt=dt)
```

where `cjump_1`

is supposed to correspond to dJ_{1,t} and `cjump_2`

is supposed to correspond to dJ_{2,t}. The resulting error:

```
ERROR: DimensionMismatch: first dimension of A, 2, does not match length of y, 3
```

suggests I actually don’t know what I’m doing. When I look back at the code, the most likely culprit (to me) is that it isn’t clear how I’m telling the solver that `cjump_1`

is the jump component of the first equation, and `cjump_2`

is the jump component in the second equation. Initially, I thought that the solver inferred which equation each jump process belonged to by their argument position in the `JumpProblem`

call, but the more I think about it, the less likely this seems. Nonetheless, that’s what I took the docs to imply.

I’d be very grateful if anyone could point out where my understanding has gone wrong here?

Cheers and thanks,

Colin