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:
- 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.
- 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
AbstractNoiseProcessthat 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.