Simulating a stochastic volatility model with jumps using DifferentialEquations.jl and JumpProcesses.jl

Specific Domains Modelling & Simulations
1

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

1 Like
2

Is there a reason why your initial condition changes to being length 3, [log(30.0), 0.0, 0.0]?

3

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?

4

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.

5

That makes perfect sense to me now that you’ve said it. Thanks again for taking the time to respond.