6
$\begingroup$

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.

Here is what I did in Mathematica (I tried to add comments in my code, please ask if there are questions/ambiguities)

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.

simulation (simulated with )

theory (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 *)
$\endgroup$
  • 3
    $\begingroup$ "I tried to add comments in my code..." - bless you, dear sir. $\endgroup$ – J. M. is away Aug 11 '17 at 14:52
  • 1
    $\begingroup$ 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) $\endgroup$ – Sjoerd Smit Aug 11 '17 at 15:58
  • $\begingroup$ Thanks for the answers! But yes, sadly I am sure....and can't find the solution anywhere... $\endgroup$ – Martin Aug 11 '17 at 16:12
3
$\begingroup$

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"}]]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – Martin Aug 11 '17 at 21:43
  • $\begingroup$ 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. $\endgroup$ – Martin Aug 12 '17 at 0:08
  • $\begingroup$ @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. $\endgroup$ – Simon Woods Aug 12 '17 at 8:04
  • $\begingroup$ Okay, thank you very much for your help! :) $\endgroup$ – Martin Aug 12 '17 at 10:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.