Resources on Modelling SDE with Parameter Switches

Specific Domains Modelling & Simulations
,
1

I’m exploring the SciMl ecosystem for the possibility of modelling something like an OU process where the drift, autoregressive, and volatility parameters are subject to discrete switches governed by a continuous-time Markov chain.

For example, let s_t \in \{1,2,\ldots,N\} be governed by the intensity matrix U, such that the transition matrix of s_t between period t and t + \tau is \exp(\tau U).

The process I’d like to model is then dX_t = \theta(s_t) (\mu(s_t) - X_t) dt + \sigma(s_t) dW_t .

As someone who traditionally works in discrete-time, I’m having a hard time mapping this problem into ModellingToolkit.jl.

2 Likes
2

You’re looking at a jump diffusion equation. See the jump diffusion tutorial: Jump Diffusion Equations · DifferentialEquations.jl

1 Like
3

Give the rate function, rate(u, p ,t) = ..., does not specify dt as a possible argument, is it correct to assume giving the proper rate for dt=1 will automatically induce the scaling invoked by \exp(\tau Q), where \tau in the previously post is equivalent to dt?

For example, with p = (Q) representing a Continuous-Time Markov chain which jumps between two states s_t = \{1,2\} would it be correct to code up the jump as:

function rate(u, p, t)
    P = exp(Q)
    if u[1]==1
        return P[1,1] 
    else
        return P[2,2]
    end
end

function affect!(integrator)
    if integrator.u[1] == 1
        integrator.u[1] += 1
    else
        integrator.u[1] -= 1
    end
end

jump = ConstantRateJump(rate, affect!)

Q = [-0.1 0.1; 0.3 -0.3]
p = Q
prob = DiscreteProblem([1], (0.0, 100.0), (Q))

jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, SSAStepper())

plot(sol)
4

My current solution for the OU process where parameters \theta(s_t), \mu(s_t), \sigma(s_t), are subject to discrete switches governed by discrete 2-element Continuous-time Markov chain, s_t \in \{1,2\}, which is parameterized by intensity matrix Q such that p(s_{t+dt} | s_t) \sim \exp(d t Q) is below.

I do not believe the solution is ideal, as I treat the discrete process, s_t, as the second element of an SDEProblem, rather than a mixture between a SDEProblem and a DiscreteProblem. Of particular concern is the call to Int(u[2]) in the rate(du, u, p, t) function, which I can’t imagine is good for performance (and may cause type-instability?).

using DifferentialEquations
using JumpProcesses
using Plots

Q = [-0.1 0.1; 0.3 -0.3]
θ = [0.4, 0.2]
μ = [0.9, -0.9]
σ = [0.1, 0.5]
p = (θ, μ, σ, Q)

function rate(u, p, t)
    P = exp(p[4])
    if u[2]==1
        return P[1,1] 
    else
        return P[2,2]
    end
end

function affect!(integrator)
    if integrator.u[2] == 1
        integrator.u[2] += 1
    else
        integrator.u[2] -= 1
    end
end

jump = ConstantRateJump(rate, affect!)

function f(du, u, p ,t) #dt
    θ = p[1]
    μ = p[2]
    du[1] = θ[Int(u[2])] * (μ[Int(u[2])] - u[1])
    du[2] = 0.0
    nothing
end

function g(du, u, p ,t) #dW
    σ = p[3]
    du[1] = σ[Int(u[2])]
    du[2] = 0.0
    nothing
end

prob = SDEProblem(f, g, [0.0, 1.0], (0.0, 1000.0), p)

jump_prob = JumpProblem(prob, Direct(), jump)

sol = solve(jump_prob, SRIW1())

plot(sol[1:200], layout = (2,1))
5

Yes. So you don’t need dt.

That’s type-stable.

1 Like
6

Thank you!