# Simulating a stochastic volatility model with jumps using DifferentialEquations.jl and JumpProcesses.jl
**URL:** https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534
**Category:** Modelling & Simulations
**Tags:** differentialequation
**Created:** [April 26, 2024, 1:43pm UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534 "2024-04-26T13:43:52Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@colintbowers](https://discourse.julialang.org/u/colintbowers)
#### Post date: [April 26, 2024, 1:43pm UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534/1 "2024-04-26T13:43:52Z")
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:
```julia
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:
```julia
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:
```julia
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
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 26, 2024, 3:41pm UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534/2 "2024-04-26T15:41:49Z")
Is there a reason why your initial condition changes to being length 3, `[log(30.0), 0.0, 0.0]`?
---
### Author:  [@colintbowers](https://discourse.julialang.org/u/colintbowers)
#### Post date: [April 27, 2024, 12:22am UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534/3 "2024-04-27T00:22:12Z")
Ahhhh crap. Yep, that’ll about do it. Sorry for wasting time with a trivial error on my part. Many thanks, somehow I’d gone around in circles for an hour without spotting that mistake.
On the plus side, does this mean my understanding above about how to add jump processes is correct? That is, does the position of `cjump_1` and `cjump_2` in the call to `JumpProblem` indicate which equation they apply to?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 27, 2024, 1:55pm UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534/4 "2024-04-27T13:55:06Z")
> [@colintbowers](#):
>
> That is, does the position of `cjump_1` and `cjump_2` in the call to `JumpProblem` indicate which equation they apply to?
Nope, just the affect! definition does that, the othering doesn’t matter. You can have 20 jumps if you want, but whether it changes value 1 or 2, or both, is just based on the definition of the affect!. So your code is right, but not because of the ordering here.
---
### Author:  [@colintbowers](https://discourse.julialang.org/u/colintbowers)
#### Post date: [April 28, 2024, 12:46pm UTC](https://discourse.julialang.org/t/simulating-a-stochastic-volatility-model-with-jumps-using-differentialequations-jl-and-jumpprocesses-jl/113534/5 "2024-04-28T12:46:29Z")
That makes perfect sense to me now that you’ve said it. Thanks again for taking the time to respond.