The Euler–Maruyama Method

Why Stochastic?¶

Many systems we want to model are not purely deterministic:

• Finance: Stock prices fluctuate due to unpredictable market forces. The Black–Scholes model for asset prices is an SDE, specifically geometric Brownian motion.

• Biology: Populations are subject to environmental noise. Deterministic logistic growth $dX = rX(1 - X/K)\,dt$ becomes stochastic when we add $\sigma X\,dW$ to model random environmental fluctuations.

• Physics: A pollen grain suspended in water gets buffeted by water molecules, the original Brownian motion observed by Robert Brown in 1827. An ODE for the grain’s position would predict it stays still (no net force). The actual path is a random, jittery curve.

Brownian Motion¶

The Drunkard’s Walk¶

Imagine a person standing at the origin on the real line. At each tick of a clock (spaced $\delta t$ apart), they step right by $+\delta$ or left by $-\delta$, each with probability $\tfrac{1}{2}$. Let $\xi_j \in \{-\delta,\,+\delta\}$ denote the $j$-th step. The steps are independent, and by symmetry each has

Here $\mathbb{E}[X]$ denotes the expected value (or mean) of a random variable $X$: the average over all possible outcomes, weighted by their probabilities. Since $\mathbb{E}[\xi_j] = 0$, the variance is $\operatorname{Var}(\xi_j) = \mathbb{E}[\xi_j^2] = \delta^2$.

After $n$ steps the position is $S_n = \sum_{j=1}^n \xi_j$. Since $S_{n+1} = S_n + \xi_{n+1}$, we can ask: if the walker is currently at position $s$, what is the expected position after one more step? In general, the conditional expectation of a random variable $Y$ given that another random variable $X$ takes the value $x$ is

the average of $Y$ restricted to only those outcomes where $X = x$. In our case $Y = S_{n+1}$, $X = S_n$, and $x = s$. Given $S_n = s$, the only two possible values of $S_{n+1}$ are $s + \delta$ and $s - \delta$, each with conditional probability $\tfrac{1}{2}$. So:

This says: no matter where the walker currently stands, the expected position after one more step is unchanged. To get the unconditional expectation, we average over all possible values of $S_n$, weighting each by its probability. Write $p_k = \mathbb{P}(S_n = k\delta)$ for the probability that the walker is at position $k\delta$ after $n$ steps. Then

In the second equality we used the conditional result above: $\mathbb{E}[S_{n+1} \mid S_n = k\delta] = k\delta$. The final sum is the definition of $\mathbb{E}[S_n]$: recall that $S_n$ takes the values $k\delta$ with probabilities $p_k$, so $\mathbb{E}[S_n] = \sum_k (k\delta)\,p_k$. (This argument is an instance of the law of total expectation.) Since $S_0 = 0$, unrolling gives $\mathbb{E}[S_n] = \mathbb{E}[S_{n-1}] = \cdots = \mathbb{E}[S_0] = 0$.

For the mean-square displacement $\mathbb{E}[S_n^2]$ (the average squared distance from the origin, taken over many realizations), square $S_{n+1} = S_n + \xi_{n+1}$ and average over the two cases:

The cross term $2\,S_n\,\delta$ cancels between the two cases. Taking expectations gives the recursion $\mathbb{E}[S_{n+1}^2] = \mathbb{E}[S_n^2] + \delta^2$, so starting from $S_0 = 0$:

After time $t = n\,\delta t$ we have $n = t/\delta t$, so the mean-square displacement is

For this to have a well-defined limit as the discretization is refined, we need $\delta^2 / \delta t$ to remain constant, i.e. $\delta \propto \sqrt{\delta t}$. The natural choice is $\delta = \sqrt{\delta t}$, giving $\mathbb{E}[S_n^2] = t$.

This is precisely what is observed experimentally: if you track a pollen grain suspended in water under a microscope, the mean-square displacement grows linearly with time. Einstein (1905) explained this by the argument above, connecting the microscopic randomness of molecular collisions to the macroscopic diffusion rate.

The Drunkard’s Walk (General Version)¶

Imagine a person standing at the origin who, at regular time intervals $\delta t$, takes a step of random size. Let $\xi_j$ denote the $j$-th step, drawn from a probability distribution with density $p(\xi)$ (so that $\mathbb{P}(a \leq \xi_j \leq b) = \int_a^b p(\xi)\,d\xi$). Suppose the steps are independent with mean zero and variance $\sigma^2$:

Here $\mathbb{E}[X]$ denotes the expected value (or mean) of a random variable $X$: the average over all possible outcomes, weighted by their probabilities.

After $n$ steps the position is $S_n = \sum_{j=1}^n \xi_j$. What are the mean and variance of $S_n$? Writing out the expected value against the joint density $p(\xi_1, \ldots, \xi_n)$:

The integral is linear in the sum. Since the steps are independent, the joint density factors as $p(\xi_1, \ldots, \xi_n) = p(\xi_1)\cdots p(\xi_n)$, and the cross-integrals collapse: $\int p(\xi_k)\,d\xi_k = 1$ for $k \neq j$. This gives

For the variance, independence is essential. Recall that $\operatorname{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$. Since $\mathbb{E}[S_n] = 0$, the variance reduces to $\operatorname{Var}(S_n) = \mathbb{E}[S_n^2]$, the mean-square displacement: the average of the squared distance from the origin, taken over many realizations of the walk. Expanding $S_n^2$:

Independence means $\mathbb{E}[\xi_i\,\xi_j] = \mathbb{E}[\xi_i]\,\mathbb{E}[\xi_j] = 0$ for $i \neq j$, so all cross terms vanish and

After time $t = n\,\delta t$ the mean-square displacement is $\mathbb{E}[S_n^2] = (t / \delta t)\,\sigma^2$. For this to depend only on $t$ and not on the discretization $\delta t$, we need $\sigma^2 \propto \delta t$, i.e. each step has standard deviation proportional to $\sqrt{\delta t}$.

This is precisely what is observed experimentally: if you track a pollen grain suspended in water under a microscope, the mean-square displacement grows linearly with time. Einstein (1905) explained this by the argument above, connecting the microscopic randomness of molecular collisions to the macroscopic diffusion rate.

Why Gaussian? The Central Limit Theorem¶

Nothing in the drunkard’s walk argument used the shape of the step distribution — only that the steps are independent with mean zero and finite variance. The $\pm\delta$ coin flip, a uniform distribution, or any other zero-mean, finite-variance law all give the same linear growth of mean-square displacement. So why does Brownian motion have Gaussian increments?

The answer is the Central Limit Theorem (CLT): if $\xi_1, \xi_2, \ldots$ are i.i.d. with mean 0 and variance $\sigma^2$, then the normalized sum converges in distribution to a Gaussian,

Equivalently, $S_n \xrightarrow{d} \mathcal{N}(0, n\sigma^2)$. In our setting $n = t/\delta t$ steps of variance $\sigma^2 = \delta t$ give $S_n \xrightarrow{d} \mathcal{N}(0, t)$, regardless of the step distribution. The Gaussian emerges as the universal limit of many small independent kicks — a coin-flip walk, a uniform walk, or molecular collisions all produce the same Gaussian in the continuum limit.

This universality is made precise by Donsker’s theorem (the functional CLT): the rescaled random walk converges to Brownian motion as a continuous process, not just at a single time. Brownian motion is the unique scaling limit of every finite-variance random walk, in the same way the Gaussian is the universal limit for sums of independent random variables. Using Gaussian increments from the start simply gives exact Brownian increments at every scale, rather than only in the limit.

From the Random Walk to Continuous Noise¶

Setting $\sigma^2 = \delta t$ and taking the steps to be Gaussian $\xi_j = \sqrt{\delta t}\,Z_j$ with $Z_j \sim \mathcal{N}(0,1)$, the random walk updates as

The $\sqrt{\delta t}$ scaling is precisely what Einstein’s argument requires. When $\delta t \to 0$, this random walk converges to a continuous random function $W(t)$ called Brownian motion (or a Wiener process).

Computing with Stochastic Processes¶

A stochastic process like $W(t)$ gives a different function on every realization. The best we can do is describe statistics of the process: expected values, variances, and correlations across realizations. When we solve an SDE numerically, the solution $X_n$ at each time step is itself a random variable that depends on all the random increments $\Delta W_0, \Delta W_1, \ldots, \Delta W_{n-1}$.

In practice, we rarely know the density of $X_n$ explicitly, so we cannot evaluate $\mathbb{E}[X_n]$ as an integral. Instead, we approximate it by Monte Carlo sampling: run $M$ independent simulations and average the results,

This is a quadrature rule for the integral $\int x\,p(x)\,dx$, where the sample points $X^{(i)}$ are drawn from the distribution $p$ rather than placed on a deterministic grid. See the companion notebook Monte Carlo Methods for Monte Carlo integration examples, importance sampling, and the Metropolis--Hastings algorithm.

From Euler’s Method to Euler–Maruyama¶

The Integral Form of an ODE (Duhamel’s Principle)¶

Start from $\frac{dx}{dt} = f(t,x)$. Integrate both sides from $t_n$ to $t_{n+1}$:

This is Duhamel’s principle. Euler’s method approximates the integral by freezing the integrand at the left endpoint (a rectangle rule):

The SDE and Its Integral Form¶

This notation is shorthand for an integral equation, Duhamel’s principle with two integrals. Over one time step $[t_n, t_{n+1}]$:

We approximate both integrals by the rectangle-rule logic:

• Deterministic integral: $\int_{t_n}^{t_{n+1}} a(s, X(s))\,ds \approx a(t_n, X_n) \cdot h$. This is Euler’s approximation.

• Stochastic integral: $\int_{t_n}^{t_{n+1}} b(s, X(s))\,dW(s) \approx b(t_n, X_n) \cdot \Delta W_n$. We freeze $b$ at the left endpoint; what remains is $\int_{t_n}^{t_{n+1}} dW(s) = W(t_{n+1}) - W(t_n) = \Delta W_n$.

The Euler–Maruyama Scheme¶

Convergence of Euler–Maruyama¶

Both $X_n$ (numerical) and $X(t_n)$ (exact) are random variables. On any particular run, the error $|X_n - X(t_n)|$ depends on the Brownian path. We need convergence concepts that account for randomness. Now that we have the expected value at our disposal, there are two natural ways to measure the error.

Strong Convergence: Mean of the Error¶

Strong convergence asks: how close is each numerical path to the true path, on average? We compute the pathwise error $|X_n - X(t_n)|$ on each realization, then take the expected value across all realizations. This matters when individual trajectories are important, for example when simulating a specific stock price path or a particular particle trajectory.

Weak Convergence: Error of the Means¶

Strong convergence is demanding: it requires every individual path to be accurate. Often we only care about statistics of the solution, for instance the expected payoff of a financial derivative or the mean concentration of a chemical species. Weak convergence asks a different question: how well does the method reproduce the expected value? Rather than averaging the error, we look at the error of the average: $|\mathbb{E}[X_n] - \mathbb{E}[X(t_n)]|$.

Stochastic Calculus and Itô’s Formula¶

To understand where the exact GBM solution comes from, and in particular the $-\frac{1}{2}\mu^2$ correction, we need a stochastic version of the chain rule. This requires first making sense of integration with respect to Brownian motion.

The Itô Stochastic Integral¶

We have already seen that the integral form of an SDE involves an expression $\int_0^T H(s)\,dW(s)$, where $H$ is some process. The Itô integral defines this as the limit of left-endpoint Riemann sums:

The left-endpoint choice is what makes this the Itô integral (as opposed to a midpoint or right-endpoint convention). It is also the choice that Euler–Maruyama naturally implements.

Two properties of the Itô integral are essential for everything that follows:

Itô’s Formula¶

Now we can derive the stochastic chain rule. Suppose $X(t)$ satisfies the SDE $dX = f(X)\,dt + g(X)\,dW$, and let $\varphi$ be a twice continuously differentiable function. We want the SDE satisfied by $Y(t) = \varphi(X(t))$.

Summary¶

• Euler–Maruyama extends Euler’s method to SDEs.

• Numerical solutions are random variables; each simulation gives one sample path.

• The $\sqrt{h}$ scaling of Brownian increments is why $(dW)^2 = dt$, why the chain rule gains a correction (Itô’s formula), and why strong convergence is order $1/2$ instead of 1.

• Strong order $1/2$, weak order 1.