# Euler–Maruyama method

In Itô calculus, the Euler–Maruyama method (also called the Euler method) is a method for the approximate numerical solution of a stochastic differential equation (SDE). It is a simple generalization of the Euler method for ordinary differential equations to stochastic differential equations. It is named after Leonhard Euler and Gisiro Maruyama. Unfortunately, the same generalization cannot be done for any arbitrary deterministic method [1].

Consider the stochastic differential equation (see Itô calculus)

${\displaystyle \mathrm {d} X_{t}=a(X_{t})\,\mathrm {d} t+b(X_{t})\,\mathrm {d} W_{t},}$

with initial condition X0 = x0, where Wt stands for 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:

• partition the interval [0, T] into N equal subintervals of width ${\displaystyle \Delta t>0}$:
${\displaystyle 0=\tau _{0}<\tau _{1}<\cdots <\tau _{N}=T{\text{ and }}\Delta t=T/N;}$
• set Y0 = x0;
• recursively define Yn for 1 ≤ n ≤ N by
${\displaystyle \,Y_{n+1}=Y_{n}+a(Y_{n})\,\Delta t+b(Y_{n})\,\Delta W_{n},}$
where
${\displaystyle \Delta W_{n}=W_{\tau _{n+1}}-W_{\tau _{n}}.}$

The random variables ΔWn are independent and identically distributed normal random variables with expected value zero and variance ${\displaystyle \Delta t}$.

## Example

### Numerical simulation

Gene expression modelled as stochastic process

An area that has benefited significantly from SDE is biology or more precisely mathematical biology. Here the number of publications on the use of stochastic model grew, as most of the models are nonlinear, demanding numerical schemes.

The graphic depicts a stochastic differential equation being solved using the Euler Scheme. The deterministic counterpart is shown as well.

### Computer implementation

The following Python code implements the Euler–Maruyama method and uses it to solve the Ornstein–Uhlenbeck process defined by

${\displaystyle dY_{t}=\theta \cdot (\mu -Y_{t})\,{\mathrm {d} }t+\sigma \,{\mathrm {d} }W_{t}}$
${\displaystyle Y_{0}=Y_{\mathrm {init} }.}$

The random numbers for ${\displaystyle {\mathrm {d} }W_{t}}$ are generated using the NumPy mathematics package.

 1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import matplotlib.pyplot as plt
4
5 num_sims = 5  # Display five runs
6
7 t_init = 3
8 t_end  = 7
9 N      = 1000  # Compute 1000 grid points
10 dt     = float(t_end - t_init) / N
11 y_init = 0
12
13 c_theta = 0.7
14 c_mu    = 1.5
15 c_sigma = 0.06
16
17 def mu(y, t):
18     """Implement the Ornstein–Uhlenbeck mu."""  # = \theta (\mu-Y_t)
19     return c_theta * (c_mu - y)
20
21 def sigma(y, t):
22     """Implement the Ornstein–Uhlenbeck sigma."""  # = \sigma
23     return c_sigma
24
25 def dW(delta_t):
26     """Sample a random number at each call."""
27     return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))
28
29 ts    = np.arange(t_init, t_end, dt)
30 ys    = np.zeros(N)
31
32 ys[0] = y_init
33
34 for _ in range(num_sims):
35     for i in range(1, ts.size):
36         t = (i-1) * dt
37         y = ys[i-1]
38         ys[i] = y + mu(y, t) * dt + sigma(y, t) * dW(dt)
39     plt.plot(ts, ys)
40
41 plt.show()


The following is simply the translation of the above code into the MATLAB (R2019b) programming language:

 1 %% Initialization and Utility
2 close all;
3 clear all;
4
5 numSims = 5;            % display five runs
6 tBounds = [3 7];        % The bounds of t
7 N      = 1000;          % Compute 1000 grid points
8 dt     = (tBounds(2) - tBounds(1)) / N ;
9 y_init = 1;             % Initial y condition
10
11
12 pd = makedist('Normal',0,sqrt(dt)); % Initialize the probability distribution for our
13                          % random variable with mean 0 and
14                          % stdev of sqrt(dt)
15
16 c = [0.7, 1.5, 0.06];   % the initial Theta, Mu, and Sigma, respectively
17
18 ts    = linspace(tBounds(1), tBounds(2), N); % From t0-->t1 with N points
19 ys    = zeros(1,N);     % 1xN Matrix of zeros
20
21 ys(1) = y_init;
22 %% Computing the Process
23 for j = 1:numSims
24     for i = 2:numel(ts)
25         t = (i-1) .* dt;
26         y = ys(i-1);
27         mu      = c(1) .* (c(2) - y);
28         sigma   = c(3);
29         dW      = random(pd);
30
31         ys(i) = y + mu .* dt + sigma .* dW;
32     end
33     figure(69)
34     hold on;
35     plot(ts, ys, 'o')
36 end