$$ \newcommand{\mcl}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\avg}[1]{\langle#1 \rangle} \newcommand{\pth}[1]{\left( #1 \right)} \newcommand{\bck}[1]{\left\{ #1 \right\}} \newcommand{\sbck}[1]{\left[ #1 \right]} \newcommand{\bsbck}[1]{\Big[ #1 \Big]} \newcommand{\deriv}[1]{\frac{\mrm{d} #1}{\mrm{d}t}} $$
I recently learned that stochastic differential equations (SDEs) had a whole lot of theory behind them and I am genuinely surprised at the complexity that a "simple" noise term seems to introduce in the numerical strategies employed to solve such equations.
I found this paper by D.Higham very instructive, and learned about Ito and Stratonovich frameworks for SDEs, but I don't understand how to get started on my own set of equations. My system is of the form:
$$ \deriv{X_k} = -X_k + f\left( \sum_i \alpha_{i,k} X_i + \eta \right) $$ where $\alpha_{i,j}$ are real constants, $f$ is a smooth non-linear function and $\eta$ is normally distributed.
From Higham's paper, I assume that the "Stochastic Chain Rule" would be relevant here, but it is unclear to me how to reduce this set of equations to an Ito or any other suitable form. Most numerical solvers I encountered so far require an Ito form, but clearly here the noise term is not directly additive.
I don't know how to go about integrating this system with noise. Without noise, I was using a simple Runge-Kutta 4 method but from what I read I assume that "adapting" this with noise is not going to be straightforward at all?
EDIT:
Based on the recommendations of @LutzL, a sensible approach seems to be as follows:
- Define the non-linear feedback term $$Z_k = \sum_i \alpha_{i,k} X_i + \eta$$ the differential of which is $$d Z_k = \sum_i \alpha_{i,k} d X_i + \sigma_k d W_k$$
- We know from the original equation that $$d X_k = ( -X_k + f(Z_k) )dt$$ so we can rewrite $$d Z_k = \sum_i \alpha_{i,k} ( f(Z_k)-X_k )dt + \sigma_k d W_k$$
- These equations are now in Ito form and can be solved using e.g. an Euler-Maruyama method. Note that if we started off with a system of $N$ equations, the resulting Ito form now comprises $2N$ equations.
Questions that persist for me are:
- @LutzL mentioned that the noise timestep $\Delta t$ has to be larger than the solver's time-step; in practice this causes the noise time-course to have unexpected spectral properties (ie, a sub-sampled white-noise is no longer white). Why is this necessary? Is there a reference in the literature about this specific problem?
- Is it ok to use a mixed method, Euler-Maruyama for $Z_k$, and RK4 for $X_k$?
- If so, should intermediary steps of $X_k$ each sample a value of $Z_k$, or should there be only one sample of $Z_k$ for all intermediary steps?
An intermediate step in making sense of your description is to choose some time step $Δt$ and holding the random input constant over intervals of size $Δt$. Obviously, this assumes that $Δt$ is somewhat larger than the step size of the numerical integrator. In this setup, i.e., where $η(t)$ is constant on the intervals $[m·Δt, (m+1)·Δt)$, the given differential equation makes perfect sense.
The obvious theoretical problem then would be how to scale this approach for increasingly smaller $Δt$ such that the resulting solutions have some kind of sensible limit in at least the probabilistic distributions, even better if some kind of path convergence could be achieved.
Note that even in the simple case of replacing $η(t)$ constant over intervals of length $Δt$ with $\tilde η(t)$ constant over intervals of half length $Δt/2$ so that locally mean and variation stay constant, if $η(m·Δt)$ is $N(0,σ^2·Δt)$ then $\tilde η(n·Δt/2)$ needs to be $N(0,σ^2·Δt/2)$ which neatly leads back to the limit of the summation processes being a Brownian motion, i.e., $η(t)=σ·ΔB(t)=σ·(B(t+Δt)-B(t))$.
In sum this leads to replacing the argument of $f$ with a random variable $Z$ with $dZ=\sum α_i·dX_i+σ·dW$ (so that it reduces to the non-random model for $σ=0$). Inserting $dX_k=(-X_k+f(Z))·dt$ leads to $$ dZ=\sum α_i·(-X_i+f(Z))·dt+σ·dW $$
Since all of $X_k,Z$ are truly random resp. proper stochastic processes, none of the SDE's are smooth enough to justify the partial use of a Runge-Kutta method. They have to be solved jointly by a stochastic integrator.