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