In Itô calculus, the Euler–Maruyama method (also simply called the Euler method) is a method for the approximate numerical solution of a stochastic differential equation (SDE). It is an extension of the Euler method for ordinary differential equations to stochastic differential equations named after Leonhard Euler and Gisiro Maruyama. The same generalization cannot be done for any arbitrary deterministic method.[1]
Consider the stochastic differential equation (see Itô calculus)
dXt=a(Xt,t)dt+b(Xt,t)dWt,
with initial condition X0 = x0, where Wt denotes the Wiener process, and suppose that we wish to solve this SDE on some interval of time [0, ''T'']. Then the Euler–Maruyama approximation to the true solution X is the Markov chain Y defined as follows:
\Deltat>0
0=\tau0<\tau1< … <\tauN=Tand\Deltat=T/N;
Yn=Yn+a(Yn,\taun)\Deltat+b(Yn,\taun)\DeltaWn,
where
\DeltaWn=
| W | |
| \taun |
-
| W | |
| \taun |
.
The random variables ΔWn are independent and identically distributed normal random variables with expected value zero and variance Δt.
An area that has benefited significantly from SDEs is mathematical biology. As many biological processes are both stochastic and continuous in nature, numerical methods of solving SDEs are highly valuable in the field.
The graphic depicts a stochastic differential equation solved using the Euler-Maruyama method. The deterministic counterpart is shown in blue.
The following Python code implements the Euler–Maruyama method and uses it to solve the Ornstein–Uhlenbeck process defined by
dYt=\theta ⋅ (\mu-Yt){d}t+\sigma{d}Wt
Y0=Yinit.
The random numbers for
{d}Wt
import numpy as npimport matplotlib.pyplot as plt
class Model: """Stochastic model constants.""" THETA = 0.7 MU = 1.5 SIGMA = 0.06
def mu(y: float, _t: float) -> float: """Implement the Ornstein–Uhlenbeck mu.""" return Model.THETA * (Model.MU - y)
def sigma(_y: float, _t: float) -> float: """Implement the Ornstein–Uhlenbeck sigma.""" return Model.SIGMA
def dW(delta_t: float) -> float: """Sample a random number at each call.""" return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))
def run_simulation: """ Return the result of one full simulation.""" T_INIT = 3 T_END = 7 N = 1000 # Compute at 1000 grid points DT = float(T_END - T_INIT) / N TS = np.arange(T_INIT, T_END + DT, DT) assert TS.size
Y_INIT = 0
ys = np.zeros(TS.size) ys[0] = Y_INIT for i in range(1, TS.size): t = T_INIT + (i - 1) * DT y = ys[i - 1] ys[i] = y + mu(y, t) * DT + sigma(y, t) * dW(DT)
return TS, ys
def plot_simulations(num_sims: int): """ Plot several simulations in one image.""" for _ in range(num_sims): plt.plot(*run_simulation)
plt.xlabel("time") plt.ylabel("y") plt.show
if __name__
The following is simply the translation of the above code into the MATLAB (R2019b) programming language:
numSims = 5; % display five runstBounds = [3 7]; % The bounds of tN = 1000; % Compute at 1000 grid pointsdt = (tBounds(2) - tBounds(1)) / N ;y_init = 0; % Initial y condition
% Initialize the probability distribution for our% random variable with mean 0 and stdev of sqrt(dt)pd = makedist('Normal',0,sqrt(dt));
c = [0.7, 1.5, 0.06]; % Theta, Mu, and Sigma, respectively
ts = linspace(tBounds(1), tBounds(2), N); % From t0-->t1 with N pointsys = zeros(1,N); % 1xN Matrix of zeros
ys(1) = y_init;%% Computing the Processfor j = 1:numSims for i = 2:numel(ts) t = tBounds(1) + (i-1) .* dt; y = ys(i-1); mu = c(1) .* (c(2) - y); sigma = c(3); dW = random(pd); ys(i) = y + mu .* dt + sigma .* dW; end figure; hold on; plot(ts, ys, 'o')end