Speed up jump diffusion solution for estimation?

General Usage
1

I am interested in estimating the parameters of a jump diffusion using simulation-based methods, and filtering statistics through a neural net. The basic idea is working quite well. Now, I need to make the simulations as fast as possible, while keeping decent accuracy, so as not to loose identification of the parameters. I need to generate many solutions at different parameter vectors. Currently, I’m using the SRIW1 algorithm (in DifferentialEquations.jl):

sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)

I am wondering if a different algorithm would be more better, to increase speed. I tried a few others, and did not see very large differences in solution time. I have experimented with julia threads and OPENBLAS threads, but I don’t see that the SRIW1 method is using either of them. I don’t have access to a useful GPU. Is there any recommendation for changing algorithms, tolerances, or means of using parallelism?

For reference, here’s the non-jump portion of the model:

function Diffusion(μ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
1 Like
2

How often are the jumps? I think it’s going to be difficult to beat what you have there (unless the jumps are rapid, in which case some of the new regular jump diffusion solvers could be applicable), and focusing on the parallelism might be the best choice.

3

What do your jumps look like? Are they mass action jumps? Where is the profile spending most of its time?

4

The jumps are infrequent (λ0 is between 0.0 and 0.05).

## jump in log price
rate(u,p,t) = ## jump in log price
rate(u,p,t) = λ0
# jump is random sign times  λ1 times current st. dev.
affect1!(integrator) = (integrator.u[1] = integrator.u[1].+rand([-1.0,1.0]).*λ1.*exp(integrator.u[2]./2.0))
jump = ConstantRateJump(rate,affect1!)
jump_prob = JumpProblem(prob,Direct(), jump)
sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)

The time span is long: (0,10000). I have parallel processing working ok, once the solution is done, but the solution itself seems to be done using only a single thread. That’s where I was hoping to find a way to parallelize. I’m using this in a MCMC algorithm, at the moment, but perhaps I will have to go to SMC, instead.

5

I’m not sure classical methods can do any better here (long time spans), but we do keep trying to optimize this more and more. saveat at longer intervals probably helps, since I assume you don’t need every step in there?

BTW, I think for estimation, directly using the PDE formulation might be more efficient if you’re fitting via distributional information. That might be something to look into.

6

I am computing realized volatility measures, so I need the solution at every regular short interval. I’m fitting via statistics, as in ABC or GMM. The solution is already pretty quick, I just wanted to make sure I was not missing some obvious opportunity to speed it up. I think my best option for speeding this up is to use posterior sampler that can sample in parallel, rather than the simple MCMC I’m currently using. Thanks!

7

Yeah, that’s probably a good way to go.

8

In the end, I have used multiple MCMC chains, which seems to work pretty well. Estimating the model using data generated at known parameter values recovers the true parameters with pretty good precision, as can be seen here https://github.com/mcreel/SNM/tree/master/examples/JD/SimulationEstimation. Thanks for the help.

1 Like
9

Perhaps too late but conditional on the jumps the parameters of the drift part are conjugate — you can just read their posterior from the trajectory with proposition 5.6 in https://arxiv.org/pdf/1712.03807.pdf if you have a continuous time observation

1 Like
10

@ChrisRackauckas Do you have conjugate analysis of drift parameters somewhere in the ecosystem?

11

The example I’m working with is observed in discrete time (end of trading day), so I believe that this proposition would not apply. Thanks for the reference, it has some pretty interesting animations.

12

R has a package called yuima that can estimate jump and other stochastic processes. I am looking for a Julia counterpart.

13

Yeah, you could still think of using the result in a data augmentation setting where you sample the joint posterior of parameters and latent path, but that’ll become a story of its own (my paper linked is concerned with that)

1 Like
14

@ChrisRackauckas Do you have conjugate analysis of drift parameters somewhere in the ecosystem?

Work in progress in Regression by frankschae · Pull Request #15 · mschauer/MitosisStochasticDiffEq.jl · GitHub