Resources on Modelling SDE with Parameter Switches

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

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

1 Like

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)



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))


Yes. So you don’t need dt.

That’s type-stable.

1 Like

Thank you!