Stochastic Differential Equations

Example 1: Scalar SDEs

In this example, we will solve the equation

\[du = f(u,p,t)dt + g(u,p,t)dW,\]

where $f(u,p,t)=αu$ and $g(u,p,t)=βu$. We know via Stochastic calculus that the solution to this equation is

\[u(t,Wₜ)=u₀\exp\left[\left(α-\frac{β^2}{2}\right)t+βWₜ\right].\]

To solve this numerically, we define a stochastic problem type using SDEProblem by specifying f(u, p, t), g(u, p, t), and the initial condition:

using DifferentialEquations
α = 1
β = 1
u₀ = 1 / 2
f(u, p, t) = α * u
g(u, p, t) = β * u
dt = 1 // 2^(4)
tspan = (0.0, 1.0)
prob = SDEProblem(f, g, u₀, tspan)

SDEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 0.5

The solve interface is then the same as ODEs. Here, we will use the classic Euler-Maruyama algorithm EM and plot the solution:

sol = solve(prob, EM(), dt = dt)
using Plots
plot(sol)

Using Higher Order Methods

One unique feature of DifferentialEquations.jl is that higher-order methods for stochastic differential equations are included. To illustrate it, let us compare the accuracy of the EM() method and a higher-order method SRIW1() with the analytical solution. This is a good way to judge the accuracy of a given algorithm, and is also useful to test convergence of new methods being developed. To setup our problem, we define u_analytic(u₀, p, t, W) and pass it to the SDEFunction as:

u_analytic(u₀, p, t, W) = u₀ * exp((α - (β^2) / 2) * t + β * W)
ff = SDEFunction(f, g, analytic = u_analytic)
prob = SDEProblem(ff, u₀, (0.0, 1.0))

SDEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 0.5

We can now compare the EM() solution with the analytic one:

sol = solve(prob, EM(), dt = dt)
plot(sol, plot_analytic = true)

Now, we choose a higher-order solver SRIW1() for better accuracy. By default, the higher order methods are adaptive. Let's first switch off adaptivity and compare the numerical and analytic solutions :

sol = solve(prob, SRIW1(), dt = dt, adaptive = false)
plot(sol, plot_analytic = true)

Now, let's allow the solver to automatically determine a starting dt. This estimate at the beginning is conservative (small) to ensure accuracy.

sol = solve(prob, SRIW1())
plot(sol, plot_analytic = true)

We can instead start the method with a larger dt by passing it to solve:

sol = solve(prob, SRIW1(), dt = dt)
plot(sol, plot_analytic = true)

Ensemble Simulations

Instead of solving single trajectories, we can turn our problem into a EnsembleProblem to solve many trajectories all at once. This is done by the EnsembleProblem constructor:

ensembleprob = EnsembleProblem(prob)

EnsembleProblem with problem SDEProblem

The solver commands are defined at the Parallel Ensemble Simulations page. For example, we can choose to have 1000 trajectories via trajectories=1000. In addition, this will automatically parallelize using Julia native parallelism if extra processes are added via addprocs(), but we can change this to use multithreading via EnsembleThreads(). Together, this looks like:

sol = solve(ensembleprob, EnsembleThreads(), trajectories = 1000)

EnsembleSolution Solution of length 1000 with uType:
RODESolution{Float64, 1, Vector{Float64}, Vector{Float64}, Dict{Symbol, Float64}, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 1, Float64, Float64, Float64, Vector{Float64}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Float64}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), typeof(Main.g), LinearAlgebra.UniformScaling{Bool}, typeof(Main.u_analytic), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, typeof(Main.g), Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, Nothing}, StochasticDiffEq.SOSRI, StochasticDiffEq.LinearInterpolationData{Vector{Float64}, Vector{Float64}}, SciMLBase.DEStats, Nothing, Nothing}

Many more controls are defined at the Ensemble simulations page, including analysis tools. A very simple analysis can be done with the EnsembleSummary, which builds mean/var statistics and has an associated plot recipe. For example, we can get the statistics at every 0.01 timesteps and plot the average + error using:

using DifferentialEquations.EnsembleAnalysis
summ = EnsembleSummary(sol, 0:0.01:1)
plot(summ, labels = "Middle 95%")
summ = EnsembleSummary(sol, 0:0.01:1; quantiles = [0.25, 0.75])
plot!(summ, labels = "Middle 50%", legend = true)

Additionally, we can easily calculate the correlation between the values at t=0.2 and t=0.7 via

timepoint_meancor(sol, 0.2, 0.7) # Gives both means and then the correlation coefficient

(0.604757885368387, 0.9681210751065629, 0.4647732320187854)

Example 2: Systems of SDEs with Diagonal Noise

In general, a system of SDEs

\[du = f(u,p,t)dt + g(u,p,t)dW,\]

where g is now a matrix of values, is numerically integrated in the same way as ODEs. A common scenario is when we have diagonal noise, which is the default for DifferentialEquations.jl. Physically this means that every variable in the system gets a different random kick. Consequently, g is a diagonal matrix and we can handle this in a simple manner by defining the deterministic part f(du,u,p,t) and the stochastic part g(du2,u,p,t) as in-place functions.

Consider for example a stochastic variant of the Lorenz equations, where we introduce a simple additive noise to each of x,y,z, which is simply 3*N(0,dt). Here, N is the normal distribution and dt is the time step. This is done as follows:

using DifferentialEquations
using Plots
function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

function σ_lorenz!(du, u, p, t)
    du[1] = 3.0
    du[2] = 3.0
    du[3] = 3.0
end

prob_sde_lorenz = SDEProblem(lorenz!, σ_lorenz!, [1.0, 0.0, 0.0], (0.0, 10.0))
sol = solve(prob_sde_lorenz)
plot(sol, idxs = (1, 2, 3))

Note that it's okay for the noise function to mix terms. For example

function σ_lorenz!(du, u, p, t)
    du[1] = sin(u[3]) * 3.0
    du[2] = u[2] * u[1] * 3.0
    du[3] = 3.0
end

σ_lorenz! (generic function with 1 method)

is a valid noise function, which will once again give diagonal noise by du2.*W.

Example 3: Systems of SDEs with Scalar Noise

In this example, we'll solve a system of SDEs with scalar noise. This means that the same noise process is applied to all SDEs. We need to define a scalar noise process using the Noise Process interface. Since we want a WienerProcess that starts at 0.0 at time 0.0, we use the command W = WienerProcess(0.0,0.0,0.0) to define the Brownian motion we want, and then give this to the noise option in the SDEProblem. For a full example, let's solve a linear SDE with scalar noise using a high order algorithm:

using DifferentialEquations
using Plots
f(du, u, p, t) = (du .= u)
g(du, u, p, t) = (du .= u)
u0 = rand(4, 2)

W = WienerProcess(0.0, 0.0, 0.0)
prob = SDEProblem(f, g, u0, (0.0, 1.0), noise = W)
sol = solve(prob, SRIW1())
plot(sol)

Example 4: Systems of SDEs with Non-Diagonal Noise

In the previous examples we had diagonal noise, that is a vector of random numbers dW whose size matches the output of g where the noise is applied element-wise, and scalar noise where a single random variable is applied to all dependent variables. However, a more general type of noise allows for the terms to linearly mixed via g being a matrix.

Note that nonlinear mixings are not SDEs but fall under the more general class of random ordinary differential equations (RODEs) which have a separate set of solvers.

Let's define a problem with four Wiener processes and two dependent random variables. In this case, we will want the output of g to be a 2x4 matrix, such that the solution is g(u,p,t)*dW, the matrix multiplication. For example, we can do the following:

using DifferentialEquations
f(du, u, p, t) = du .= 1.01u
function g(du, u, p, t)
    du[1, 1] = 0.3u[1]
    du[1, 2] = 0.6u[1]
    du[1, 3] = 0.9u[1]
    du[1, 4] = 0.12u[1]
    du[2, 1] = 1.2u[2]
    du[2, 2] = 0.2u[2]
    du[2, 3] = 0.3u[2]
    du[2, 4] = 1.8u[2]
end
prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

In our g we define the functions for computing the values of the matrix. We can now think of the SDE that this solves as the system of equations

\[du_1 = f_1(u,p,t)dt + g_{11}(u,p,t)dW_1 + g_{12}(u,p,t)dW_2 + g_{13}(u,p,t)dW_3 + g_{14}(u,p,t)dW_4 \\ du_2 = f_2(u,p,t)dt + g_{21}(u,p,t)dW_1 + g_{22}(u,p,t)dW_2 + g_{23}(u,p,t)dW_3 + g_{24}(u,p,t)dW_4\]

meaning that for example du[1,1] and du[2,1] correspond to stochastic changes with the same random number in the first and second SDEs.

The matrix itself is determined by the keyword argument noise_rate_prototype in the SDEProblem constructor. This is a prototype for the type that du will be in g. This can be any AbstractMatrix type. Thus, we can define the problem as

# Define a sparse matrix by making a dense matrix and setting some values as not zero
using SparseArrays
A = zeros(2, 4)
A[1, 1] = 1
A[1, 4] = 1
A[2, 4] = 1
A = sparse(A)

# Make `g` write the sparse matrix values
function g(du, u, p, t)
    du[1, 1] = 0.3u[1]
    du[1, 4] = 0.12u[2]
    du[2, 4] = 1.8u[2]
end

# Make `g` use the sparse matrix
prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = A)

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

and now g(u,p,t) writes into a sparse matrix, and g(u,p,t)*dW is sparse matrix multiplication.

Example 4: Colored Noise

Colored noise can be defined using the Noise Process interface. In that portion of the docs, it is shown how to define your own noise process my_noise, which can be passed to the SDEProblem

SDEProblem(f, g, u0, tspan, noise = my_noise)

Note that general colored noise problems are only compatible with the EM and EulerHeun methods. This is discussed in the SDE solvers page.

Example: Spatially-Colored Noise in the Heston Model

Let's define the Heston equation from financial mathematics:

\[dS = μSdt + \sqrt{v}SdW_1 \\ dv = κ(Θ-v)dt + σ\sqrt{v}dW_2 \\ dW_1 dW_2 = ρ dt\]

In this problem, we have a diagonal noise problem given by:

function f(du, u, p, t)
    du[1] = μ * u[1]
    du[2] = κ * (Θ - u[2])
end
function g(du, u, p, t)
    du[1] = √u[2] * u[1]
    du[2] = σ * √u[2]
end

g (generic function with 1 method)

However, our noise has a correlation matrix for some constant ρ. Choosing ρ=0.2:

ρ = 0.2
Γ = [1 ρ; ρ 1]

2×2 Matrix{Float64}:
 1.0  0.2
 0.2  1.0

To solve this, we can define a CorrelatedWienerProcess which starts at zero (W(0)=0) via:

tspan = (0.0, 1.0)
heston_noise = CorrelatedWienerProcess!(Γ, tspan[1], zeros(2), zeros(2))

t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.0, 0.0]

This is then used to build the SDE:

SDEProblem(f, g, ones(2), tspan, noise = heston_noise)

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

Of course, to fully define this problem, we need to define our constants. Constructors for making common models like this easier to define can be found in the modeling toolkits. For example, the HestonProblem is pre-defined as part of the financial modeling tools.