I was hoping to use Turing.jl alongside DifferentialEquations.jl to estimate a partially observable SDE via quasi-likelihood methods. Being new to continuous-time models, I was hoping to get some suggestions on the most efficient way to proceed.
In the model below, z(t), π(t), and r(t) are assumed to be observed at discrete time-intervals without measurement error. While zg(t) and zπ(t) are not observed by the statistician, and thus are referred to as “latent factors”. These latent factors impact the drift of g(t) and π(t) respectively, and non-linearly impact the effect of g(t) and π(t) on D(r(t)). Lastly, note that the diffusion matrix is assumed to be diagonal, and all variables have constant drift with the exception of r(t), whose diffusion coefficient is proportional to sqrt(r(t)), resembling a Cox-Ingeroll-Ross interest rate model.
My guess is this kind of estimation requires an extended Kalman filter or a particle filter in order to compute the likliehood function, but I am not sure which is more suitable to my problem and easiest to implement in the SciML and Turing frameworks.
EDIT: After coming upon this (https://www.sas.upenn.edu/~jesusfv/hmc_dssm.pdf), maybe HMC is the way to go while sampling the latent states as part of the estimation.
using ModelingToolkit
using StochasticDiffEq
using Plots
using NaNMath
@variables t g(t) π(t) r(t) zg(t) zπ(t)
@parameters α_g α_π α_r κ_g κ_π κ_r κ_zg κ_zπ ϕ_g ϕ_π ζ_0g ζ_1g ζ_2g ζ_0π ζ_1π ζ_2π ζ_gπ σ_g σ_π σ_r σ_zg σ_zπ
D = Differential(t)
eqs = [D(g) ~ α_g + zg - κ_g * g - ϕ_g * r + ϕ_g * π,
D(π) ~ α_π + zπ - κ_π * π + ϕ_π * g,
D(r) ~ α_r - κ_r * r + (ζ_0g + ζ_1g * zg + ζ_2g * zg^2) * g +
(ζ_0π + ζ_1π * zπ + ζ_2π * zπ^2) * π + ζ_gπ * zg * zπ,
D(zg) ~ - κ_zg * zg,
D(zπ) ~ - κ_zπ * zπ
]
noiseeqs = [σ_g,
σ_π,
σ_r * NaNMath.sqrt(r),
σ_zg,
σ_zπ
]
@named macro_mdl = SDESystem(eqs, noiseeqs, t, [g, π, r, zg, zπ], [α_g, α_π, α_r, κ_g, κ_π, κ_r, κ_zg, κ_zπ, ϕ_g, ϕ_π, ζ_0g, ζ_1g, ζ_2g, ζ_0π, ζ_1π, ζ_2π, ζ_gπ, σ_g, σ_π, σ_r, σ_zg, σ_zπ])
u0map = [g => 0.0, π => 0.0075, r => 0.02, zg => 0.0, zπ => 0.0]
parammap = [α_g => 0.005, α_π => 0.005, α_r => 0.01, κ_g => 0.5, κ_π => 0.5, κ_r => 0.5, κ_zg => 0.5, κ_zπ => 0.5, ϕ_g => 0.5, ϕ_π => 0.5, ζ_0g => 0.01, ζ_1g => 0.1, ζ_2g => 0.2, ζ_0π => 0.0, ζ_1π => 0.3, ζ_2π => 0.2, ζ_gπ => 0.1, σ_g => 0.01, σ_π => 0.01, σ_r => 0.01, σ_zg => 0.01, σ_zπ => 0.01]
prob = SDEProblem(macro_mdl, u0map, (0.0,10000.0), parammap)
sol = solve(prob, SOSRI())
plot(sol)