Stochastic filtration problem

Hi all,

I’m working with a physical system whose state u(t) follows an SDE
du=f(u,t)dt+g(u,t)dW_t

The system produces a measured signal Y_t given by
dY_t = h(u,t)dt+dW_t
where dW_t is the same Wiener increment in both equations.

In my experiment I know the initial state u_0 and I measure Y_t at discrete times. Using the observation equation, I can formally substitute dW_t = dY_t - h(u,t)dt into the equation, giving
du = f(u,t)dt + g(u,t)\bigl(dY_t - h(u,t)dt\bigr) = (f(u,t) - g(u,t)h(u,t))dt + g(u,t)dY_t.

My goal: to numerically propagate u(t) using the measured signal Y(t), ideally with something better than fixed-step Euler–Maruyama.
Trying to improve the integration method, I ran into two problems:

  1. dY_t is sampled at discrete timesteps. To derive it midstep, we would need to interpolate Y_t. This may require corrections due to Itô vs. Stratonovich interpretations, which I’m not sure how to derive.
  2. I’m not sure how to implement the numerical integration. ODE methods seem inappropriate because dY_t scales like \sqrt{dt}, not dt. For SDE methods, I need to somehow substitute the innovation term dW_t = dY_t - h(u,t)dt as the noise process. I considered defining a custom AbstractNoiseProcess that stores the timeseries of Y and derives dW_t based on that. However, these methods typically require a pseudoprocess dZ which I’m not sure if and how I can similarly emulate.

Any help would be greatly appreciated.

p.s. Extra structure in my equations: f and h are linear in u, while g is quadratic.

Have you considered using something like an extended or unscented Kalman filter that supports correlation between process and measurement noise to do this?

If you want to be particular, you could linearize f and g at each time step and sample the covariance matrix \nabla g (\nabla g)^T assuming the linearized dynamics (implemented in ControlSystemsBase.c2d)

Here’s a tutorial that shows how to use one particular filtering package with correlated process and measurement noise