# Tips for Estimating Non-linear Latent SDE

**URL:** <https://discourse.julialang.org/t/tips-for-estimating-non-linear-latent-sde/104326>\
**Category:** Statistics\
**Tags:** inference, turing, nonlinear, differentialequation\
**Created:** [September 27, 2023, 7:45pm UTC](https://discourse.julialang.org/t/tips-for-estimating-non-linear-latent-sde/104326 "2023-09-27T19:45:52Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![dleather](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dleather/32/43762_2.png) [@dleather](https://discourse.julialang.org/u/dleather)\
**Post date:** [September 27, 2023, 7:45pm UTC](https://discourse.julialang.org/t/tips-for-estimating-non-linear-latent-sde/104326/1 "2023-09-27T19:45:52Z")

</div>

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](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.

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 7, 2024, 1:42pm UTC](https://discourse.julialang.org/t/tips-for-estimating-non-linear-latent-sde/104326/2 "2024-04-07T13:42:20Z")

</div>

> [@dleather](#):
>
> EDIT: After coming upon this ([https://www.sas.upenn.edu/~jesusfv/hmc\_dssm.pdf](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.

Yes, I would suggest HMC as a good way to do this. Doing one solve of the SDE as the likelihood is an unbiased estimate of the likelihood, so doing HMC over that is a quasi-likelihood formulation that is fast and pretty robust.