Heston SDE exact methods

Specific Domains Numerics
, ,
1

I am new to DifferentialEquations.jl and I tried out some SDE simulations with EM.
I implemented an exact simulation of Heston SV model using the paper of Broadie and Kaya (https://www.researchgate.net/publication/220243902_Exact_Simulation_of_Stochastic_Volatility_and_Other_Affine_Jump_Diffusion_Processes)
I used the sampling method they propose in the dist function, following the guide: Abstract Noise Processes · DiffEqNoiseProcess.jl.

Is this the proper way to go?

How are almost-exact methods supported? For example commonly the volatility process is exactly simulated or approximated with QE and then the price process is approximated with Euler scheme.
I am not sure if a blueprint for such a scheme is already there or not.

By the way, I am relatively new to Julia, and I find the resources on SCIML in general very nice and well documented👍

2

Yeah this part isn’t so well documented.

You can just define the SDE to be solved. Some of them are found here:

But this space needs a bit of work.

3

Hi, thanks for that.
I went through it and through the NoiseProcess documentation too.
Two considerations. Heston process can be simulated using the transition distribution directly. I created a sampler and a NoiseProcess taking the existing GeometricBrownianMotionProcess as a blueprint.
It looks like this:

# the first component of W is the log of price in heston, the second its vol
function HestonNoise2D(μ, κ, θ, σ, ρ, t0, W0, Z0 = nothing; kwargs...)
    @inline function Heston_dist(DW, W, dt, u, p, t, rng)
        S, V = exp(W[end][1]), W[end][2]
        heston_dist_at_t = HestonDistribution2D(S, V, κ, θ, σ, ρ, μ, t + dt)
        return @fastmath rand(rng, heston_dist_at_t; kwargs...)  # Calls exact Heston sampler
    end

    return NoiseProcess{false}(t0, W0, Z0, Heston_dist, nothing)
end

The one that is on DiffEqFinancial works only with EM or similar discretization methods, not distributionally exact, if I understand correctly.

Also, that one has the issue that variance can become negative for some sets of parameters, making the sqrt() fail. I solved it like this, maybe there’s a better way using built-in functions?

function HestonProblem(μ, κ, Θ, σ, ρ, u0, tspan; seed = UInt64(0), kwargs...)
    f = function (u, p, t)
        adj_var = max(u[2], 0)
        return @. [μ * u[1], κ * (Θ - adj_var)]
    end
    g = function (u, p, t)
        adj_var = sqrt(max(u[2], 0))
        return @. [adj_var * u[1], σ * adj_var]
    end
    Γ = [1 ρ; ρ 1] 

    noise = CorrelatedWienerProcess(Γ, tspan[1], zeros(2))

    sde_f = SDEFunction(f, g)
    return SDEProblem(sde_f, u0, (tspan[1], tspan[2]), noise=noise, seed=seed, kwargs...)
end

Here I just “absorbed” vol using 0 if it becomes negative. I did not actually set it to zero. There are multiple strategies to handle this. Would one have to create a different SDEProblem for each strategy (Abs, Max, etc)?

My curiosity was also about the following: what if I have an exact sampler for one component of my 2D process and I want to discretize the other component?