# Driving a system of differential equations with an AR1-Process

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?

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]


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]


• Works very well, thank you Chris! Oct 24, 2018 at 14:31