# Basic Gaussian Process simulation/ white noise - massive problem

I have a certain problem with Mathematica or more precisely, myself handling it. It is correlated with a mathematical problem and its realization/simulation in Mathematica.

I want to start with the mathematical part: I am given a gaussian distributed random variable $\epsilon(t)$ with $E\left[\epsilon(t)\right]=0$ and correlation function $E\left[\epsilon(t)\epsilon(s)\right]=\sigma_{\epsilon}^2\delta(t-s)$, so gaussian white noise. The thing i want to calculate is $E\left[\exp\left(-\frac{1}{2}\int_{0}^{t}\epsilon(t')\mathrm{d}t'\right)\right]$. We see that the integral of a gaussian-distributed random variable is again gaussian distributed, so from the moment generating function of the gaussian distribution we get: \begin{equation*} E\left[\exp\left(k X\right)\right]=G(k)=\exp\left(\mu_X k+\frac{1}{2}k^2\sigma_X^2 \right) \end{equation*}

With $X=\int_{0}^{t}\epsilon(t')\mathrm{d}t'$ it follows that $\mu_X=E\left[X\right]=0$ and $E\left[X^2\right]=\sigma_{\epsilon}^2t$.

So it follows: \begin{equation*} G\left(-\frac{1}{2}\right)=E\left[\exp\left(-\frac{1}{2}\int_{0}^{t}\epsilon(t')\mathrm{d}t'\right)\right]=\exp\left(\frac{1}{8}\sigma^2_\epsilon t\right) \end{equation*}

So far so good, this is the thing I want to 'get as a result' theoretically. Now I wanted to verify it with a little simulation in Mathematica.

I really hope you guys can help me, I can't work any further without having solved that certain problem...and I'm getting desperate. I also seems that the result depends on the time-steps I choose... But the plots really don't match.

(simulated with )

(theoretically)

Standarddeviation=5;
(*define standard deviation of the gaussian random variable epsilon *)

DeltaT=0.04;
(* time steps with which i want to evaluate the integral *)

Repetitions=5000;
(* in order to get a good average, i take 5000 reps*)

Time = 20;
(*maximal time *)

AverageExactly[t_] := Exp[1/8*Standarddeviation^2*t];
(*exact average theoretically*)

(*It follows the loop for the 5000 repetitions *)
For[z = 0, z < Repetitions, z++,
GaussianDistribution=RandomVariate[NormalDistribution[0,Standarddeviation],
Round[Time/DeltaT + 1]];
(*generate Gaussian numbers, there shall be enough that all
time is covered*)
SimulatedValue[t_] := Exp[-1/2*Deltat*Sum[GaussianDistribution[[i]],{i,1,t/Deltat}]];
(* Evaluate the integral in form of a sum*)
ListSimulatedValues =
Table[SimulatedValue[t], {t, 0, Time, DeltaT}];
(*make a table out of the simulated values*)
ListAxis = Table[t, {t, 0, Time, Deltat}];
(*define the axis-values*)

ListSimulatedValuesPlot[z] =
Transpose@{# & @@@ ListAxis, # & @@@
ListSimulatedValues}; (*merge both tables*)

]

SimulatedAveraged = List[1/Repetitions*Sum[ListSimulatedValuesPlot[z],{z, 0,  Repetitions - 1}]]
(*average over all 5000 repetitions *)
ListPlot[SimulatedAveraged =     List[1/Repetitions*Sum[ListSimulatedValuesPlot[z], {z, 0, Repetitions - 1}]],
PlotRange -> {{0, 20}, {0.9, 2}}] (*Plot simulated values *)
Plot[AverageExactly[t], {t, 0, 20}](*Plot theoretical values values *)

• "I tried to add comments in my code..." - bless you, dear sir. – J. M.'s ennui Aug 11 '17 at 14:52
• Without having tried the code myself, are you sure the plots are different? Look at the y-axes of both plots; they're completely different so you really can't tell if the lines are the same or not. (Also: if you're expecting exponential behaviour, use ListLogPlot to make the plot more readable) – Sjoerd Smit Aug 11 '17 at 15:58
• Thanks for the answers! But yes, sadly I am sure....and can't find the solution anywhere... – Martin Aug 11 '17 at 16:12

See this for example. The simulation is band-limited, requiring a correction to the variance of the distribution. If you do that and stick to low values for time the agreement is okay.

sd = 5;
dt = 0.001;
reps = 10000;
time = 1;

exact[t_] := Exp[1/8*sd^2*t];

gd = RandomVariate[NormalDistribution[0, sd/Sqrt[dt]], {Round[time/dt + 1], reps}];

sim = Mean[Exp[-(dt/2) Transpose@Accumulate[gd]]];

Show[
Plot[exact[t], {t, 0, time}, PlotLegends -> {"Exact"}],
ListLinePlot[sim, DataRange -> {0, time}, PlotStyle -> Orange,
PlotLegends -> {"Simulation"}]]


• Thanks a lot! I don't really understand why you have to scale the standard deviation with $\delta t$ yet, but i hope i will understand it later. The problem is: how would i do this for other processes (not gaussian distributed, but for example the Ornstein-Uhlenbeck-process). I wouldn't even have any 'theoretical' value to test my simulations – Martin Aug 11 '17 at 21:43
• Another question: for my certain problem I need the simulation to be accurate for 'longer times'. in fact I have a function that is multiplied with this integral/sum, and I want it to be accurate up to time about 20. – Martin Aug 12 '17 at 0:08
• @Martin I don't fully understand the maths either. The simulation stays accurate up to longer times if you increase the number of reps - so how far you can go will largely depend on how long you are prepared to wait. – Simon Woods Aug 12 '17 at 8:04
• Okay, thank you very much for your help! :) – Martin Aug 12 '17 at 10:22