Specify jumps in a Heston-like model?

Specific Domains Finance and Economics
1

I am working with a model similar to the Heston model, and would like to add jumps. I would like to have the jump rate which depend on the ordinary volatility as
rate1(u,p,t) = λ0 .+ λ1.*exp(u[2]/2.0)
I understand that this should be specified as a variable rate jump. However, I can’t seem to get the code to run if I do that. It does run if I specify it as a constant rate jump, but I’m doubting if the solution is correct when I do that. The code for the jump part follows. Any pointers would be welcome!

(this post is substantially similar to and earlier post, but with the diffeq tag added, sorry for the noise)

## jump in price
rate1(u,p,t) = λ0 .+ λ1.*exp(u[2]/2.0) # volatility dependent jump rate
# jump is normal with st. dev. equal to λ1 times current st. dev.
affect1!(integrator) = (integrator.u[1] = integrator.u[1].+randn(size(integrator.u[1])).*λ2.*exp(integrator.u[2]./2.0))

# this works:
jump1 = ConstantRateJump(rate1,affect1!)
# this does not
#jump1 = VariableRateJump(rate1,affect1!)
jump_prob = JumpProblem(prob,Direct(), jump1)
sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)
2

You may also want to look at https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl for simulating this, using a trivial flow. At least, you will have a way to compare to another numerical solution.

2 Likes
3

That doesn’t support SDEs?

@mcreel I was planning to look at your post today. Sorry about that!

1 Like
4

He did not give the flow in between jumps, I assumed it is constant

5

It’s from his other post, and the solver is SRIW1

6

Here’s the whole file, if it would help. As is, it runs and gives output that looks as expected, I’m just not sure if the jumps are occurring at the proper rate. I also don’t know how to determine when a jump has occurred.

using DifferentialEquations, Plots

function MyProblem(μ0,μ1,κ,α,σ,ρ,u0,tspan)
    f = function (du,u,p,t)
        du[1] = μ0 + μ1*(u[2]-α)/σ # drift in log prices
        du[2] = κ*(α-u[2]) # mean reversion in shocks
    end
    g = function (du,u,p,t)
        du[1] = exp(u[2]/2.0)
        du[2] = σ
    end
    Γ = [1.0 ρ;ρ 1.0] # Covariance Matrix
    noise = CorrelatedWienerProcess!(Γ,tspan[1],zeros(2),zeros(2))
    sde_f = SDEFunction{true}(f,g)
    SDEProblem(sde_f,g,u0,tspan,noise=noise)
end


function main()
# trading days, closing, etc

# assume trading period is 1/3 of day (8 hours) 
# but that latent price evolves continuously
# observed return is daily difference of log price at
# closing time
TradingDays = 1000 # total days in sample
Days = Int(TradingDays*7/5)  # calendar days
MinPerDay = 1440 # minutes per day
MinPerTic = 5 # minutes between tics, lower for better accuracy
tics = Int(MinPerDay/MinPerTic) # number of tics in day
dt = 1/tics # divisions per day
closing = Int(floor(tics/3)) # closing tic: closing happens after 1/3 of day

# parameters
μ0 = 0.0
μ1 = 0.0
κ = 0.1
α = 0.15
σ = 0.15
ρ = -0.7
λ0 = 1.0 # constant in jump rate
λ1 = 1.0 # slope in jump rate
λ2 = 3.0 # size of jumps
σme = 0.05 # standard dev of measurement error in returns
u0 = [0;α]
prob = MyProblem(μ0, μ1, κ, α, σ, ρ, u0, (0.0,Days))

## jump in price
rate(u,p,t) = λ0 .+ λ1.*exp(u[2]/2.0) # volatility dependent jump rate
# jump is normal with st. dev. equal to λ1 times current st. dev.
affect1!(integrator) = (integrator.u[1] = integrator.u[1].+randn(size(integrator.u[1])).*λ2.*exp(integrator.u[2]./2.0))

# this works:
jump = ConstantRateJump(rate,affect1!)
# this does not
#jump = VariableRateJump(rate,affect1!)
jump_prob = JumpProblem(prob,Direct(), jump)
sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)

# get log price at end of trading days
z = zeros(TradingDays+1)
global j = 0 # counter for day of week
global k = 0 # counter for trading days
for i = 0:(Days)
    # set day of week, and record if it's a trading day
    j +=1
    if j<6
        k +=1 # advance trading day
        z[k]=(sol.u)[i*tics+closing][1]
    end
    if j==7 # restart the week if Sunday
        j = 0
    end   
end
z[2:end]-z[1:end-1] + σme*randn(TradingDays) # returns are diff of log price
end
z = main()

plot(z)
7

This one is going to be quite tough. The first thing to work out is that VariableRateJumps grow the size of the system, so the problem was that you needed to define a 3 noise process for this to work:

using StochasticDiffEq, DiffEqJump, DiffEqNoiseProcess, Plots

function MyProblem(μ0,μ1,κ,α,σ,ρ,u0,tspan)
    f = function (du,u,p,t)
        du[1] = μ0 + μ1*(u[2]-α)/σ # drift in prices
        du[2] = κ*(α-u[2]) # mean reversion in shocks
    end
    g = function (du,u,p,t)
        du[1] = exp(u[2]/2.0)
        du[2] = σ
    end
    Γ = [1 ρ 0;ρ 1 0;0 0 0] # Covariance Matrix
    noise = CorrelatedWienerProcess!(Γ,tspan[1],zeros(3),zeros(3))
    sde_f = SDEFunction{true}(f,g)
    SDEProblem(sde_f,g,u0,tspan,noise=noise)
end

μ0 = 0.0
μ1 = 0.0
κ = 0.05
α = 0.2
σ = 0.2
ρ = -0.7
λ0 = 2.0
λ1 = 3.0
u0 = [0;α]
dt = 0.01
prob = MyProblem(μ0, μ1, κ, α, σ, ρ, u0, (0.0,1000.0))

## jump in price
#rate1(u,p,t) = λ0 # constant jump rate
rate1(u,p,t) = λ0.*exp(u[2]) # volatility dependent jump rate
# jump is normal with st. dev. equal to λ1 times current st. dev.
affect1!(integrator) = (integrator.u[1] = integrator.u[1].+randn(size(integrator.u[1])).*λ1.*exp(integrator.u[2]./2.0))

# this works:
#jump1 = ConstantRateJump(rate1,affect1!)
# this does not
jump1 = VariableRateJump(rate1,affect1!)
jump_prob = JumpProblem(prob,Direct(), jump1)
sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)

However, this now errors at the jump stage because continuous callbacks require pulling back the noise process, so it requires the definition of a bridging distribution. If you’re willing to work out the bridging distribution for the correlated noise process then we’re done. Basically, we need the distribution for if you know W0 and Wh, what’s the distribution of Wt in the middle.

1 Like
8

Thanks for the explanation. I will fall back to using a constant rate jump, which will be good enough for my purpose, which is simply to illustrate an econometric estimator. I do have one remaining question, though, how can one see the realization of jumps is a simulation, their timing and size? Thanks again.

9

If you look for non-unique time points you can see the jump at those spots.

10

Thanks, that does the trick. Here’s a plot of returns with the jumps highlighted in red:
rets_and_jumps

In this example, a jump is a normally distributed r.v. with mean zero, so not all jumps are outliers.

If anyone would like to see how to find jumps, here’s the code:
check jumps