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}];
(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}]
(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
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]
The next difficult step is to tell ItoProcess
to take GaussianProcess
as a legitimate argument.