# Stability of the numerical methods for SDE

I've been figuring out with the methods for integrating of stochastic differential equations in Mathematica. I've considered the one-dimensional system: $$dx=-x dt+\sigma x dw$$ with some initial condition. I've hoped to see how the trajectory goes to zero while diminishing the initial condition till zero due to the theorem which states that if the deterministic system is stable then the stochastic system of the form given above is either stable for any value of the noise intensity $$\sigma$$.

However, I've found that the trajectory bursts close to the end of the integration interval. Moreover, integrating of the same system for different time intervals always finishes with bursting close to the end of the interval. Can it be the problem of the numerical methods or its realization in Mathematica? Why does the solution lose its stability?

I've tried each numerical method for the system both in Ito and Sratonovich forms.

Code example:

\[Sigma] = 20;
Sc = 10000;
proc = StratonovichProcess[\[DifferentialD]x[t] == -x[t] \[DifferentialD]t + (\[Sigma] x[t] ) \[DifferentialD]w[t],x[t], {x, 1/Sc}, t, w \[Distributed] WienerProcess[]];
rf = RandomFunction[proc, {0., 20., 0.01},
Method -> "KloedenPlatenSchurz"];
ListLinePlot[rf, Filling -> Axis, PlotRange -> Full]


The strength of noise σ = 20 seems large so that it can overwhelm the mean trend. Luckily Mathematica can give some analytical insights into this stochastic process.

Clear[σ]
proc = StratonovichProcess[\[DifferentialD]x[t] ==
-x[t] \[DifferentialD]t + (σ x[t]) \[DifferentialD]w[t],
x[t], {x, x0}, t, w \[Distributed] WienerProcess[]];
Mean[proc[t]]
Variance[proc[t]]

(* E^(1/2 t (-2 + σ^2)) x0 *)
(* E^(-2 t + t σ^2) (-1 + E^(t σ^2)) x0^2 *)


Looks like the mean goes to zero only if σ<Sqrt[2] and the variance blows up if σ>1.

ItoProcess gives slightly different results:

proc = ItoProcess[\[DifferentialD]x[t] ==
-x[t] \[DifferentialD]t + (σ x[t]) \[DifferentialD]w[t],
x[t], {x, x0}, t, w \[Distributed] WienerProcess[]];
Mean[proc[t]]
Variance[proc[t]

(* E^-t x0 *)
(* E^(-2 t) (-1 + E^(t σ^2)) x0^2 *)


Here the mean goes to zero in all cases, but the variance blows up if σ>Sqrt[2].

• Thanks for showing how the moments behave. Indeed, these moments are solution of the system of ordinary differential equation which is of closed form for any linear SDE. I see, the system is stable in mean-square sense only when $\sigma<Sqrt[2]$ for the SDE in both forms. However, originally I meant stochastic stability (stability in probability) which is weaker than that in p-moment sense for any p>0. The second part of the question should be neglected, I suppose, because it's hardly possible to infer something about the probability measure by the values of distinct random trajectories. – Artem Zefirov Nov 4 '18 at 17:08
• This is already outside my area of expertise, but PDF[proc[t]][x] might also be useful. – Chris K Nov 4 '18 at 17:17
• The density function $p(t,x)$ goes to $0$ uniformly by $x$ when $t \rightarrow \infty$ and has its peak at $x=0$. It follows that $X(t)=0$ is stable in distribution for any $\sigma$ as it should be. – Artem Zefirov Nov 4 '18 at 17:42
• I've found that decreasing of the step size seems to make the computed trajectories more stable. The method is explicit so, as my intuition said to me at the beginning, one have to take into account stability of the numerical method which is conditional for explicit methods and depends on the step size. – Artem Zefirov Nov 5 '18 at 19:31