5
$\begingroup$

I have the following system of differential equations:

v[t_] := RandomVariate[NormalDistribution[]];
sol = NDSolve[{x1'[t] == -0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 
      1.22 x4[t] + v[t], 
    x2'[t] == 
     0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] + 0.28 x4[t] + v[t], 
    x3'[t] == 
     0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t] + v[t], 
    x4'[t] == 
     0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t] + v[t], 
    x1[0] == 0, x2[0] == 0, x3[0] == 0, x4[0] == 0}, {x1, x2, x3, 
    x4}, {t, 0, 60}];
graphAll = 
 Plot[Evaluate[{x1[t], x2[t], x3[t], x4[t]} /. sol], {t, 0, 60}, 
  PlotRange -> All]

I want to drive the system with an AR1-Process. So far I was only able to specify some random variate from a normal distribution, but it does not seem to work properly.

Any ideas how to do this?

$\endgroup$

1 Answer 1

6
$\begingroup$

You need to use RandomFunction[ItoProcess[]] instead of NDSolve. The syntax is a frustratingly a little different than NDSolve. Driven by a WienerProcess as in your example:

sol = RandomFunction[ItoProcess[{
  \[DifferentialD]x1[t] == \[DifferentialD]v[t] + 
    (-0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 1.22 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x2[t] == \[DifferentialD]v[t] + 
    (0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] +0.28 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x3[t] == \[DifferentialD]v[t] +
    (0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x4[t] == \[DifferentialD]v[t] +
    (0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t]) \[DifferentialD]t},
  {x1[t], x2[t], x3[t], x4[t]}, {{x1, x2, x3, x4}, {0, 0, 0, 0}}, t,
  v \[Distributed] WienerProcess[0, 1]], {0, 100, 0.01}];
ListLinePlot[sol, PlotRange -> All]

Mathematica graphics

For the AR-1 part, I think you need to add a first-order decay equation for v[t] driven by a WienerProcess like:

ϕ = 1;
σ = 1;
sol = RandomFunction[ItoProcess[{
  \[DifferentialD]x1[t] == (v[t] - 0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 1.22 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x2[t] == (v[t] + 0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] + 0.28 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x3[t] == (v[t] + 0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x4[t] == (v[t] + 0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t]) \[DifferentialD]t,
  \[DifferentialD]v[t] == -ϕ v[t] \[DifferentialD]t + σ \[DifferentialD]W[t]},
  {x1[t], x2[t], x3[t], x4[t], v[t]}, {{x1, x2, x3, x4, v}, {0, 0, 0, 0, 0}}, t,
  W \[Distributed] WienerProcess[0, 1]], {0, 100, 0.01}];
ListLinePlot[sol]

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ Works very well, thank you Chris! $\endgroup$
    – holistic
    Oct 24, 2018 at 14:31

Your Answer

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

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