14
$\begingroup$

Question

Is it possible to extend the function ItoProcess so that it takes correlated noise?

I.e would like to be able to write

 eqn = Thread[{x'[t], y'[t]} == 1/10 {-y[t], x[t]} + {n1[t], n2[t]}];
 proc = ItoProcess[eqn, {x[t], y[t]}, {{x, y}, {1, 1}}, t, {n1 \[Distributed] GaussianProcess[],  n2 \[Distributed] GaussianProcess[]}];

as a means of solving a Stochastic Differential Equation sourced by correlated random noise. So I would like to have a function called GaussianProcess and that ItoProcess understands it. (If I understand correctly the Documentation WienerProcess is uncorrelated).

Is it quite possible this functionality exists with a different name?

Attempt

I know how to do this using the GaussianRandomField (defined in this thread) and NDSolve

 noise = Interpolation[#][t] & /@ Table[GaussianRandomField[256, 1, Function[k, k^-2]] // Chop, {2}];
 Plot[noise, {t, 1, 256}];

enter image description here

(Note the level of correlation in the noise)

After interpolation, I can integrate using NDSolveValue

eqn = Thread[{x'[t], y'[t]} == 1/10 {-y[t], x[t]} + noise];
eqn = Join[eqn, {x[1] == 1, y[1] == 1}];
sol = NDSolveValue[eqn, {x, y}, {t, 1, 256}];
ParametricPlot[#[t] & /@ sol, {t, 1, 256}]

enter image description here

(amusing random plot!)

Yet, it would be great IMHO if a function like GaussianProcess existed and could be fed to the ItoProcess framework.

Motivation

One advantage of such solution would be to inherit all the wrapping that ItoProcess has, while being able to specify the exact PowerSpectrum of the Gaussian noise.

PS: For clarity, let me replicate here the GaussianRandomField function that @HenrikSchumacher wrote

GaussianRandomField[
  size : (_Integer?Positive) : 256, dim : (_Integer?Positive) : 1, 
  Pk_: Function[k, k^-3]] := Module[{Pkn, fftIndgen, noise, amplitude, s2},
  s2 = Quotient[size, 2];
  fftIndgen = N@ArrayPad[Range[0, s2], {0, s2 - 1}, "ReflectedNegation"];
  amplitude = Sqrt[Outer[Plus, Sequence@@ ConstantArray[fftIndgen^2, dim],dim]];
  amplitude[[Sequence @@ ConstantArray[1, dim]]] = 1.;
  amplitude = Pk[amplitude];
  amplitude[[Sequence @@ ConstantArray[1, dim]]] = 0.;
  noise = Fourier[RandomVariate[NormalDistribution[], ConstantArray[size, dim]]];
  Re[InverseFourier[noise amplitude]]
  ]

Note that this code provides the opportunity to generate correlated Gaussian Random Processes,

 tt = GaussianRandomField[128, 2];tt[[;; , ;; 8]] // Transpose // ListLinePlot

enter image description here

which could be handy.

Complement

As a possible wrapper to GaussianRandomField one could define

Clear[gaussianProcess, GaussianProcess];
gaussianProcess[R : (_?Positive) : 1, L : (_?Positive) : 10, 
  dx : (_?Positive) : 1/100, nb : (_Integer?Positive) : 1] := 
 Module[{tt, k, nn}, 
  If[dx > R/2, Print["Insufficient Sampling"]; Abort[]];
  tt = Table[tt = GaussianRandomField[nn = Round[L/dx];
      nn = If[OddQ[nn], nn + 1, nn], 1, 
      Function[k, Exp[-1/2 R^2 (2 Pi k/L)^2]]];
    tt /= StandardDeviation[tt], {nb}];
  If[nb == 1, tt = tt[[1]]];
  TemporalData[tt, {0., L}, 
   ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1, 
     Method -> "Spline"}]];    
   GaussianProcess /: 
   RandomFunction[ GaussianProcess[R_], {0, t1_, dt_}, n_: 1] := 
    gaussianProcess[R, t1, dt, n]

Then RandomFunction produces the timelines as it does with WienerProcess.

  dat = RandomFunction[GaussianProcess[0.1], {0, 20, 0.01}, 5]

enter image description here

The next difficult step is to tell ItoProcess to take GaussianProcess as a legitimate argument.

$\endgroup$
1
  • 2
    $\begingroup$ It seems an interesting problem. But I guess it would be better if some notes briefing the theoretical fundamentals are provided. $\endgroup$ Apr 23, 2020 at 8:10

1 Answer 1

4
+75
$\begingroup$

It is not a solution, but some comparison two different approach to the same problem (oscillator with random force). What we expecting from ItoProcess in this case? We looking for model of force, and there are several possibilities for ItoProcess as option for process. There are also several method of solution. Nevertheless there is one combination similar to the possible solution:

ito = ItoProcess[{\[DifferentialD]x[
      t] == -1/10 y[t] \[DifferentialD]t + \[DifferentialD]w1[t]/
      15, \[DifferentialD]y[t] == 
    1/10 x[t] \[DifferentialD]t + \[DifferentialD]w2[t]/15}, {x[t], 
   y[t]}, {{x, y}, {1, 1}}, 
  t, {w1 \[Distributed] OrnsteinUhlenbeckProcess[0, 1, .03, 1], 
   w2 \[Distributed] OrnsteinUhlenbeckProcess[0, 1, .03, 1]}];
 path = RandomFunction[ito, {1., 256, .05}, 1, Method -> "Milstein"];
 dat = Flatten[Transpose@path["ValueList"], 1];
 {ListLinePlot[path], 
 Graphics[Table[{Hue[i/Length[dat]], Point[dat[[i]]]}, {i, 
    Length[dat]}], AspectRatio -> Automatic, Frame -> True, 
  FrameTicks -> Automatic]}    

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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