Modeling with Stochasticity

All previous differential equations tutorials deal with deterministic Systems. In this tutorial, we add randomness. In particular, we show how to represent a stochastic differential equation as a SDESystem.

Note

The high level @mtkmodel macro used in the getting started tutorial is not yet compatible with SDESystem. We thus have to use a lower level interface to define stochastic differential equations. For an introduction to this interface, read the programmatically generating Systems tutorial.

Let's take the Lorenz equation and add noise to each of the states. To show the flexibility of ModelingToolkit, we do not use homogeneous noise, with constant variance, but instead use heterogeneous noise, where the magnitude of the noise scales with (0.3 times) the magnitude of each of the states:

\[\begin{aligned} \frac{dx}{dt} &= (\sigma (y-x)) &+ 0.3x\frac{dB}{dt} \\ \frac{dy}{dt} &= (x(\rho-z) - y) &+ 0.3y\frac{dB}{dt} \\ \frac{dz}{dt} &= (xy - \beta z) &+ 0.3z\frac{dB}{dt} \\ \end{aligned}\]

Where $B$, is standard Brownian motion, also called the Wiener process. We use notation similar to the Langevin equation, often used in physics. By "multiplying" the equations by $dt$, the notation used in probability theory can be recovered.

We use this Langevin-like notation because it allows us to extend MTK modeling capacity from ODEs to SDEs, using only a single new concept, @brownians variables, which represent $\frac{dB}{dt}$ in the above equation.

using ModelingToolkit, StochasticDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using Plots

@parameters σ=10.0 ρ=2.33 β=26.0
@variables x(t)=5.0 y(t)=5.0 z(t)=1.0
@brownians B
eqs = [D(x) ~ σ * (y - x) + 0.3x * B,
    D(y) ~ x * (ρ - z) - y + 0.3y * B,
    D(z) ~ x * y - β * z + 0.3z * B]

@mtkcompile de = System(eqs, t)

\[ \begin{align} \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} &= x\left( t \right) y\left( t \right) - z\left( t \right) \beta \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= - y\left( t \right) + x\left( t \right) \left( - z\left( t \right) + \rho \right) \\ \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= \left( - x\left( t \right) + y\left( t \right) \right) \sigma \end{align} \]

Even though we did not explicitly use SDESystem, ModelingToolkit can still infer this from the equations.

typeof(de)

System

We continue by solving and plotting the SDE.

prob = SDEProblem(de, [], (0.0, 100.0))
sol = solve(prob, SRIW1())
plot(sol, idxs = [(1, 2, 3)])

The noise present in all 3 equations is correlated, as can be seen on the below figure. The figure also shows the multiplicative nature of the noise. Because states x and y generally take on larger values, the noise also takes on a more pronounced effect on these states compared to the state z.

plot(sol)

If you want uncorrelated noise for each equation, multiple @brownians variables have to be declared.

@brownians Bx By Bz
eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx,
    D(y) ~ x * (ρ - z) - y + 0.3y * By,
    D(z) ~ x * y - β * z + 0.3z * Bz]
@mtkcompile de = System(eqs, t)
prob = SDEProblem(de, [], (0.0, 100.0))
sol = solve(prob, SRIW1())
plot(sol)