# Some inconsistencies in an ensemble of stochastic differential equations
**URL:** https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639
**Category:** General Usage
**Tags:** question, differentialequation
**Created:** [June 12, 2022, 1:59pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639 "2022-06-12T13:59:26Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [June 12, 2022, 1:59pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/1 "2022-06-12T13:59:26Z")
I am solving a stochastic differential equation using `EnsembleProblem`.
When I then simply plot a summary of the solution, I get a reasonable result.
However, when I calculate the mean myself using `timeseries_point_meanvar` I find that the _same_ solution is unstable.
Here is the ensemble summary
And here is the calculated mean with `timeseries_point_meanvar`
The code I am running for the 1st plot is:
```julia
using DifferentialEquations
using Plots
function rayleigh(du,u,p,t)
T, M, m, L, num_bath_particles = p
Utag(x) = -T * (L - 2* x) / (x*(L-x))
Gamma(x) = num_bath_particles * (1/x + 1/(L-x)) * sqrt(8*π * m * T)
du[1] = u[2]
du[2] = (-Utag(u[1]) - Gamma(u[1]) * u[2]) / M
end
function σ_rayleigh(du,u,p,t)
T, M, m, L, num_bath_particles = p
Gamma(x) = num_bath_particles * (1/x + 1/(L-x)) * sqrt(8*π * m * T)
du[1] = 0
du[2] = sqrt(2* abs(Gamma(u[1])) * T) / M
end
T = 1.
M = 10^6
m = 1.
L = 10.
num_bath_particles = 1
p = [T, M, m, L, num_bath_particles]
x0 = 2.5
v0 = 0.0
prob_sde_rayleigh = SDEProblem(rayleigh, σ_rayleigh, [x0, v0], (0.0,10^6), p)
ensembleprob = EnsembleProblem(prob_sde_rayleigh)
dt = 1000.
sim = solve(ensembleprob, EnsembleThreads(), trajectories=1000, saveat=dt)
using DifferentialEquations.EnsembleAnalysis
summ = EnsembleSummary(sim, 0:dt:10^6)
plot(summ,labels="Middle 95%", idxs=1)
```
and for the 2nd plot is
```julia
dt = 1000
t_to_run = 10^6
t_arr = 1:dt:t_to_run
m_series, v_series = timeseries_point_meanvar(sim, 1:dt:t_to_run)
x_mean= zeros(length(dt:dt:t_to_run))
x_var= zeros(length(dt:dt:t_to_run))
v_mean= zeros(length(dt:dt:t_to_run))
v_var = zeros(length(dt:dt:t_to_run))
for i=1:length(dt:dt:t_to_run)
x_var[i] = v_series.u[i][1]
x_mean[i] = m_series.u[i][1]
v_var[i] = v_series.u[i][2]
v_mean[i] = m_series.u[i][2]
end
plot(x_mean)
```
Why is there such inconsistency?
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 12, 2022, 8:21pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/2 "2022-06-12T20:21:26Z")
I think the top plot is displaying the median and the bottom plot is displaying the arithmetic mean. For multiplicative noise processes these tend to diverge.
cmp
> **[Ergodicity Breaking in Geometric Brownian Motion](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.110.100603)**
>
> Geometric Brownian motion (GBM) is a model for systems as varied as financial instruments and populations. The statistical properties of GBM are complicated by nonergodicity, which can lead to ensemble averages exhibiting exponential growth while any...
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 12, 2022, 9:31pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/3 "2022-06-12T21:31:34Z")
Yes, IIRC the default for the summary plot is to use quantiles, so it’ll be showing the median and the 95% quantiles, which in some cases of course can be very different from the mean.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [June 13, 2022, 7:16am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/4 "2022-06-13T07:16:04Z")
@josuagrw thanks very much. I did not realize the summary gives the median.
@ChrisRackauckas what is IIRC?
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 13, 2022, 8:54am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/5 "2022-06-13T08:54:25Z")
IIRC — if I recall correctly
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 13, 2022, 8:58am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/6 "2022-06-13T08:58:50Z")
Just as a sidenote: If your model is correct, this could be a quite interesting scientific result. It would mean that there are very rare cases which have such intensity that they overshadow the more frequent “normal” cases. Nassim Taleb writes extensively on the risk implications for such systems.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [June 13, 2022, 9:16am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/7 "2022-06-13T09:16:17Z")
Yes, that’s a good point. Since you told me that I was actually looking at the median, I am trying to figure out if I have a problem with my model.
The coordinate I am looking at should be bounded in some range (in this case between 0 and 10), and this is modeled by an infinite force at x=0, 10 which rises in a continuous though fast manner. I suspect that the solver manages to jump over this barrier and then go wild.
I am looking at very large times, so the solver might take too coarse timesteps, even though I was thinking that it should find appropriate timesteps.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 13, 2022, 9:18am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/8 "2022-06-13T09:18:49Z")
You might need Stratonovich integration to describe this process correctly…
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [June 13, 2022, 9:21am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/9 "2022-06-13T09:21:43Z")
thanks for the tip! I will try to figure it out.
---
### Author:  [@patrick](https://discourse.julialang.org/u/patrick)
#### Post date: [June 13, 2022, 4:07pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/10 "2022-06-13T16:07:02Z")
I don’t think this is the case. For any Itô SDE, one can write down an equivalent Stratonovich SDE, and similarly for the reverse (perhaps assuming mild technical hypotheses). For a brief reference, see [here](https://www.robots.ox.ac.uk/~lsgs/posts/2018-09-30-ito-strat.html).
Maybe the modeling would be more intuitive/practical in the Stratonovich interpretation, however.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 13, 2022, 4:43pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/11 "2022-06-13T16:43:05Z")
What I mean is, there might be a mismatch between the form of the model (might be stratonovich, common in physics) and the numerical solver (which is ito by default).
---
### Author:  [@patrick](https://discourse.julialang.org/u/patrick)
#### Post date: [June 13, 2022, 5:12pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/12 "2022-06-13T17:12:19Z")
My apologies, I misunderstood your comment!
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [June 14, 2022, 8:35am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/13 "2022-06-14T08:35:55Z")
I tried some of the Stratonovich solvers, like `SRA` and `SOSRA` etc., (without putting too much thought into it yet) and some did yield a better result. This means that for some ensembles the mean explodes and some do not. I will further investigate this.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [June 14, 2022, 8:45am UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/14 "2022-06-14T08:45:33Z")
That’s interesting. I want to stress @patrick s point: Switching between Ito and Stratonovich integration changes the _meaning_ of the parameters in your model. So you need to make sure beforehand in which form your SDE is formulated and then choose the correct type of solver for that formulation.
The article @patrick linked to is indeed extremely enlightening on that topic.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 14, 2022, 6:40pm UTC](https://discourse.julialang.org/t/some-inconsistencies-in-an-ensemble-of-stochastic-differential-equations/82639/15 "2022-06-14T18:40:21Z")
> [@roi.holtzman](#):
>
> I tried some of the Stratonovich solvers, like `SRA` and `SOSRA` etc., (without putting too much thought into it yet) and some did yield a better result. This means that for some ensembles the mean explodes and some do not. I will further investigate this.
SRA and SOSRA aren’t not for Ito or Stratonovich necessarily, just additive noise, and additive noise is the same in both interpretations.