Stochastic differential equations

How to add noise to a differential equation without angering the mathematicians: an ensemble of sample paths, the Euler–Maruyama scheme, and the odd rule that noise scales with the square root of the time step.

Almost every model I work with has noise in it. Synaptic input arrives at irregular times, spontaneous activity in the cortex fluctuates, and even a tidy rate model gets a noise term the moment it is compared with data. The usual way to write this down is to take an ordinary differential equation and add a random forcing term. That looks harmless, and in practice it mostly is, but the step from “add noise” to “a stochastic differential equation” hides a few genuine surprises. In 2022 I gave a short talk about these at our lab’s coffee meeting under the title “How to avoid angering the mathematicians”. The first figure grows a bundle of sample paths while you watch, and the second lets you break the simulation on purpose by picking the wrong time step.

From a spring to a random walk

My talk started with a spring, because that is where the trouble starts. A mass on a damped spring, pushed by an external force \(w(t)\), obeys

\[\frac{d^2x}{dt^2} + \gamma\,\frac{dx}{dt} + \nu^2 x = w(t),\]

with \(x\) the position, \(\gamma\) the damping and \(\nu\) the natural frequency. This is a second-order linear ordinary differential equation, and everything about it is well understood: it has exactly one solution for a given start, and you can find it by hand or with a standard numerical scheme such as Euler or Runge–Kutta. Now let the force be random. The natural choice is white noise, a process \(w(t)\) that has zero mean, is Gaussian, and whose values at any two different times are independent. The spring is now a stochastic differential equation, and the textbook guarantees no longer apply.

The reason is that white noise is a strange object. Independence at every pair of times means it is discontinuous everywhere, and it takes arbitrarily large values in any interval, however short. The theory of ordinary differential equations asks for a continuous right-hand side, so the usual existence and uniqueness theorem cannot even be stated, and the integral \(\int w(t)\,dt\) that would give the solution does not exist as an ordinary integral because \(w\) is unbounded. The way out is to stop talking about \(w\) and talk instead about its running sum, Brownian motion \(\beta(t)\). Brownian motion is continuous, so integrals against its increments \(d\beta\) can be made to work. Formally \(w = d\beta/dt\), but \(\beta\) has no derivative anywhere, which is exactly why \(w\) is so badly behaved. This is the origin of Itô calculus, and it is also why its chain rule has an extra term that the ordinary one does not: a tiny increment \(d\beta\) is not negligible when squared, because \((d\beta)^2\) behaves like \(dt\). All of this is why a mathematician winces when \(w(t)\) is written into an equation as if it were an ordinary function. The mathematics allows the equation, as long as you say what you mean by the integral.

The name “white noise”, by the way, comes from the spectrum. By the Wiener–Khinchin theorem (Wiener 1930; Khintchine 1934), the power spectrum of a stationary process is the Fourier transform of its autocorrelation, and the delta-shaped autocorrelation of \(w\) transforms into a flat spectrum, all frequencies present in equal measure, as in white light. Brownian motion, the running sum of white noise, has a spectrum that falls as \(1/f^2\) and is called brown noise, after the motion rather than the colour, and the \(1/f\) spectrum in between is pink. That deserves a page of its own, and it is on the list.

The two equations on this page are the ones you are most likely to meet. Written in the Itô form, they are

\[\text{Brownian motion:}\qquad dx = D\, d\beta\]

\[\text{Ornstein–Uhlenbeck:}\qquad dx = -\theta\, x\, dt + D\, d\beta .\]

The first is a random walk: \(x\) wanders with no preference for any value and spreads out forever, which is also the equation for the position of a diffusing particle. The second adds a spring that pulls \(x\) back towards zero with rate \(\theta\), so the wandering is kept in check. The Ornstein–Uhlenbeck process (Uhlenbeck and Ornstein 1930) was invented for the velocity of a Brownian particle, and it is the standard model for anything that fluctuates around a set point, from membrane potentials to interest rates. The third option in the figure is the noisy spring itself, with \(\nu = 2\) fixed and the slider setting the damping \(\gamma\). Written as a first-order system in position and velocity it is

\[\text{Noisy spring:}\qquad dx = v\, dt, \qquad dv = (-\gamma\, v - \nu^2 x)\, dt + D\, d\beta .\]

Because these equations are linear, their solutions can be described completely. The whole ensemble of possible paths is Gaussian at every time, with a mean that follows the noise-free equation and a variance that grows in a known way. Starting from \(x_0\), Brownian motion has mean \(x_0\) and variance \(D^2 t\), so its spread grows without limit as \(\sqrt{t}\). The Ornstein–Uhlenbeck process has mean \(x_0 e^{-\theta t}\) and variance \(\frac{D^2}{2\theta}\left(1 - e^{-2\theta t}\right)\), which settles at \(D^2/(2\theta)\). The first figure draws these formulas over the simulated ensemble, so you can check them.

Solving them numerically

Most of the time nobody solves these equations by hand. The simplest scheme is Euler–Maruyama, which is the Euler method with one change. For an equation \(dx = f(x)\,dt + L(x)\,d\beta\) it reads

\[x_{k+1} = x_k + f(x_k)\,\Delta t + L(x_k)\,\Delta\beta_k, \qquad \Delta\beta_k \sim \mathcal N(0, \Delta t),\]

so each step adds a Gaussian random number with variance \(\Delta t\), or equivalently a standard normal number multiplied by \(\sqrt{\Delta t}\). That square root is the whole point. If you treat the noise like any other term and multiply it by \(\Delta t\), the random part shrinks faster than the deterministic part as the step gets smaller, and in the limit it vanishes. Put differently, multiplying the noise by \(\Delta t\) is the same as running Euler–Maruyama with a noise strength \(D\sqrt{\Delta t}\): the scheme is internally consistent, but it solves a different equation, and which one depends on the time step. The second figure shows this directly. Higher-order schemes exist, Milstein being the most common, but for the models on this page Euler–Maruyama with a small step is all you need.

Reading the figure

In the first figure, the panel on the left shows twelve sample paths in orange, all started at \(x_0 = 2\) and simulated with Euler–Maruyama at \(\Delta t = 0.01\). Behind them, the blue line is the mean of a larger ensemble of 300 paths and the blue band is that ensemble’s mean plus and minus one standard deviation, as the legend under the panel says. The dashed black lines are the formulas for the Brownian and Ornstein–Uhlenbeck cases, mean plus and minus one standard deviation, so when the simulation is right the band and the dashed lines should agree up to sampling noise. The noisy spring has no dashed lines; its ensemble statistics would need a longer formula, and the band speaks for itself.

The panel on the right is a histogram of where the 300 paths are at the current moment, turned on its side so that it shares the vertical axis with the paths. The dashed curve is the Gaussian predicted by the formulas, or, for the spring, the Gaussian with the ensemble’s own mean and standard deviation. The three numbers above it are the ensemble mean, the ensemble standard deviation and the predicted standard deviation. The histogram is rescaled to fit its panel, so compare widths along the vertical axis, not bar lengths. Changing a slider or the model restarts the simulation from \(t = 0\).

The second figure runs the same random walk, \(dx = d\beta\) with \(D = 1\), twice with the time step you choose. On the left each step adds a standard normal number times \(\sqrt{\Delta t}\), on the right the same number times \(\Delta t\). The dashed lines in both panels are the same thing, the exact spread \(\sqrt{t}\) of the equation being solved, and the number under each panel is the measured standard deviation of 500 paths at \(t = 10\). The third panel collects that number for every available time step, blue for Euler–Maruyama and orange for naive Euler, against the exact value \(\sqrt{10} \approx 3.16\) as a dashed line, with the current time step drawn larger.

Things to try

Start with the Ornstein–Uhlenbeck process at the default values. The mean decays from 2 towards 0 with time constant \(1/\theta\), and the band widens until it reaches a standard deviation of \(D/\sqrt{2\theta}\), which is 0.71 for the defaults. Once the mean has settled, the histogram on the right stops changing shape even though every individual path keeps moving. That is the stationary distribution, and it is what people mean when they say that an Ornstein–Uhlenbeck process “has” a variance. Raise \(\theta\) and the band narrows and settles faster. Lower it and the process starts to look like Brownian motion for a while before the pull makes itself felt.

Switch to Brownian motion. Now nothing pulls \(x\) back, the mean stays at 2 and the band keeps growing as \(\sqrt t\), reaching a standard deviation of \(D\sqrt{10} \approx 3.2\) for \(D = 1\) at the end of the run. The histogram just keeps flattening. Set \(D\) to 2 and watch how many paths leave the plot: the standard deviation at \(t = 10\) is now 6.3, well beyond the axis.

Switch to the noisy spring. With \(\gamma = 1\) the mean oscillates a few times before it dies away, and the band is narrower than for the Ornstein–Uhlenbeck process at the same \(D\), because the noise drives the velocity rather than the position and the spring filters it. Set the damping to 0.1 and the paths ring at the natural frequency for the whole run, still roughly in step with each other, while the mean decays only slowly and the band widens a little with every swing. Ten time units is too short to see where this ends up, but the trend is the point: the noise slowly knocks the paths out of phase with one another, so over a long enough run the mean fades away while every single path keeps oscillating. The ensemble mean is then a poor description of any member of the ensemble, which is worth remembering whenever a model is compared with trial-averaged data.

In the second figure, begin with \(\Delta t = 1\). The two path panels are indistinguishable, because \(\sqrt{1} = 1\), and the two dots in the third panel coincide on the dashed line. Go the other way first, to \(\Delta t = 2\) and then 5. With only five steps in the whole run the paths are crude in both panels, but the blue spread stays near 3.2 while the orange one climbs to 4.5 and then 7.1: a coarse step overstates the noise in the naive scheme just as a fine step understates it. Now reduce the step. The left panel stays the same apart from looking smoother, and its spread stays near 3.2 at every step: the blue dots sit on the line. The right panel shrinks. Its paths still spread as \(\sqrt t\), since it is a random walk like any other, only one with a weaker noise, and the spread at \(t = 10\) drops to 1.0 at \(\Delta t = 0.1\) and 0.32 at \(\Delta t = 0.01\). The orange dots fall away from the line as \(\sqrt{\Delta t}\). A numerical scheme is supposed to give an answer that stops depending on \(\Delta t\) once the step is small enough, and the naive one never does. In the limit of a vanishing step its noise would disappear altogether. This is the practical face of everything in the section on white noise. A random force that leaves a finite mark in a finite time must have an infinite amplitude in the limit, and the \(\sqrt{\Delta t}\) in the Euler–Maruyama scheme is the bookkeeping that makes that work on a computer. If a simulation of a noisy model gives results that change when you refine the time step, this is the first thing to check.

This page follows my 2022 coffee-meeting talk, which drew mostly on the lecture notes by Särkkä (2012), since published as a book (Särkkä and Solin 2019). Higham (2001) is the friendliest introduction to the numerical side and the source of the Euler–Maruyama figure in the talk. A disclaimer is in order: I once tried studying mathematics alongside biophysics and gave it up in favour of the physics part of my degree, so the mathematicians in the title may be angered by this page after all. If you are one of them, Särkkä and Solin (2019) has the proper version.

Numerical details

Both figures use the Euler–Maruyama scheme exactly as written above, with standard normal numbers from the Box–Muller transform. The first figure integrates 300 paths with \(\Delta t = 0.01\) from \(x_0 = 2\) (and \(v_0 = 0\) for the spring) to \(T = 10\), draws 12 of them, and computes the mean and standard deviation across all 300 at every step; the spring uses \(\nu = 2\). The second figure integrates 500 paths of \(dx = d\beta\) with \(D = 1\) to \(T = 10\) at the chosen \(\Delta t\) and draws 15 of them; the convergence panel repeats the 500-path run for each of the nine available time steps with both schemes. The naive Euler scheme, with the noise increment multiplied by \(\Delta t\) instead of \(\sqrt{\Delta t}\), is used nowhere except in the right panel of that figure, where it is the point.

References

Higham DJ (2001). An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review 43: 525–546. DOI · PDF

Itô K (1944). Stochastic integral. Proceedings of the Imperial Academy 20: 519–524.

Khintchine A (1934). Korrelationstheorie der stationären stochastischen Prozesse. Mathematische Annalen 109: 604–615. DOI

Särkkä S (2012). Applied Stochastic Differential Equations, lecture notes, Aalto University. PDF

Särkkä S, Solin A (2019). Applied Stochastic Differential Equations. Cambridge University Press. DOI

Uhlenbeck GE, Ornstein LS (1930). On the theory of the Brownian motion. Physical Review 36: 823–841. DOI

Vatiwutipong P, Phewchean N (2019). Alternative way to derive the distribution of the multivariate Ornstein–Uhlenbeck process. Advances in Difference Equations 2019: 276. DOI · PDF

Wiener N (1930). Generalized harmonic analysis. Acta Mathematica 55: 117–258.

Lecture notes on stochastic systems, Institute for Dynamic Systems and Control, ETH Zürich. PDF