Brownian Motions and Quantifying Randomness in Physical Systems

Stochastic calculus has come a long way since Robert Brown described the motion of pollen through a microscope in 1827. It’s now a key player in data science, quant finance, and mathematical biology. This article is drawn from notes I wrote for an undergraduate statistical physics course a few months ago. There won’t be any mathematical rigor.

Brownian processes (and a Wetherspoons customer)

In 1d, the term Brownian motion is reserved for continuous functions ##W(t)## that satisfy three key conditions:

1. The motion starts at zero: ##W(0) = 0##
2. Its increments ##W(t_{k+1}) – W(t_k)## for ##0 \leq t_1 \leq \dots \leq t_n## are independent of each other
3. Its increments are normally distributed ##W(t_{k+1}) – W(t_k) \sim N(0, t_{k+1} – t_k)## with variance equal to the time increment.

You might like to verify that the following well-known properties are true: (i) ##E(W(t)) = 0##, (ii) ##E(W(t)^2) = t## and (iii) ##W(t) \sim N(0,t)##, (iv) the continuous time Markov property:

$$P(X(t) = j | X(t_1) = k_1, \dots, X(t_n) = k_n) = P(X(t) = j | X(t_n) = k_n)$$Also, (v) Brownian motions are martingales, which means:

$$E(W(t+s) | W(t)) = W(t)$$

Martingales have cool properties. One is that a martingale stopped at a stopping time is a martingale (what does that mean?).

Consider a random walk consisting of several independent and identically distributed (IID) random variables ##X_i## which take values ##\pm 1## with probabilities ##p_{\pm} = 1/2##. Define ##S_n = X_1 + \dots + X_n##. This ##S_n## is a symmetric random walk: notice ##E(S_n) = 0##.

The quantity ##S_n^2 – n## is also a martingale. If ##S_n## is known, then ##S_{n+1}## is either ##S_{n} + 1## or ##S_{n}-1## with equal probabilities of ##1/2##, so:

$$E(S_{n+1}^2 – (n+1)) = \frac{1}{2}[(S_{n}+1)^2 + (S_{n}-1)^2] – (n+1) = S_{n}^2 – n$$

This fact will be useful in the next example:

Example: A drunken martingale

To give an example of why martingale-ness is a useful property, consider the following problem. A drunk man starts at the ##n^{\mathrm{th}}## pavement tile on a 100-tile street. He moves either ##+1## or ##-1## tile each step, with a probability of ##1/2## to go forward or backward. What is the probability that he ends up at the ##100^{\mathrm{th}}## tile (as opposed to tile ##0##), and what’s his expected number of steps?

You can solve this with a more direct probability argument, but the martingale approach is faster. Transform the coordinates so that his current position is defined as tile zero so that he ends at either tile ##-n## or tile ##100-n##. His path is now a symmetric random walk starting at zero and ending at one of these points.

Let ##N## be the number of steps until he reaches one of the endpoints. Because ##S_N## and ##S_N^2 – N## are martingales stopped at stopping time (with the stopping time being ##N##), these random variables are also martingales, hence ##E(S_N) = S_0 = 0## and ##E(S_N^2 – N) = S_0^2 – 0 = 0##.

If ##p_0## and ##p_{100}## are the probabilities of arriving at tile 0 and 100 respectively, then you can also figure out directly from the definition of the expectation that

$$E(S_N) = -n p_{0} + (100-n) p_{100}= 0$$

and

$$E(S_N^2 – N) = (-n)^2 p_{0} + (100-n)^2 p_{100} – E(N) = 0$$

You can check that, after putting ##p_{0} + p_{100} = 1## and solving, you obtain the probability of reaching the ##100^{\mathrm{th}}## tile as ##p_{100} = n/100## and the expected number of steps ##E(N) = n(100-n)##.

Stochastic differential equations (SDEs)

A typical stochastic differential equation for a random process ##X(t)## is

$$dX(t) = \mu(X, t)dt + \sigma(X,t) dW$$In this equation ##\mu(X,t) dt## is a drift term, and ##\sigma(X,t) dW## is a diffusive term which captures the stochastic character of the process. Mathematicians will tell you that the SDE above is an informal statement of the integral equation

$$X(t+\Delta t) – X(t) = \int_t^{t+\Delta t} \mu(X(t’),t’) dt’ + \int_t^{t+\Delta t} \sigma(X(t’), t’) dW(t’)$$ Stochastic processes obey their own slightly unique calculus (stochastic calculus). A famous result is the Itô Lemma, which describes how to take the differential of a function ##f(X)## of a stochastic process.

Kiyosi Itô

To arrive at the result, use a Taylor expansion:

$$df = f(X+dX, t+dt) – f(x,t) = \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial X} dX + \frac{1}{2} \frac{\partial^2 f}{\partial X^2} dX^2 + \dots$$

The ##dX^2## term is retained because, from the SDE, it is first order in ##dt## (because ##dW^2 = dt##, so ##dX^2 = \sigma^2 dt + \dots##). The Itô Lemma results from inserting the SDE into the Taylor expansion:

$$df = \left( \frac{\partial f}{\partial t} + \mu \frac{\partial f}{\partial X} + \frac{1}{2} \sigma^2 \frac{\partial^2 f}{\partial X^2} \right)dt + \sigma \frac{\partial f}{\partial X} dW$$

Example: Geometric Brownian Motion (GBM)

A famous SDE is

$$dX = \mu X dt + \sigma X dW$$

for constants ##\mu## and ##\sigma##. The SDE can be solved by considering the transformation ##f(X) = \log X##, so that ##\partial f/\partial X = 1/X## and ##\partial^2 f/\partial X^2 = -1/X^2##. Obviously ##\partial f/\partial t=0##. The Itô Lemma gives the following result:

$$df = \left( \mu – \frac{1}{2} \sigma^2\right)dt + \sigma dW$$

Integrating and exponentiating, the behavior of ##X(t)## is exponential growth multiplied by a stochastic factor ##e^{\sigma W(t)}##.

$$X = X_0 \exp \left[ \left( \mu – \frac{1}{2} \sigma^2\right)t + \sigma W(t) \right]$$

Simulation of GBM for 20 random paths.

Kramers-Moyal Expansion

Now for a change of perspective. Consider instead the full probability ##P(X,t)## that a stochastic process takes a value ##X## at time ##t##. How can we calculate the probability density at time ##t + \Delta t##? The famous Chapman-Kolmogorov theorem is helpful:

$$P(X, t+\Delta t) = \int G(X, t+\Delta t | X’, t) P(X’, t) dX’$$

The function ##G## is the propagator (the conditional probability for each path). It’s a similar concept to the propagator in a path integral. Make a sneaky substitution ##\Delta X := X-X’##, and write

$$P(X, t+\Delta t) = \int G(X + \Delta X – \Delta X, t+\Delta t | X – \Delta X, t) P(X-\Delta X, t) d(\Delta X)$$(Eagle-eyed readers will have noticed that the negative sign in ##d(-\Delta X)## has been canceled by an un-shown reversal of the limits of integration.) This form is convenient because you can do a Taylor expansion in ##-\Delta X##,

$$P(X, t+\Delta t) = \int \sum_{n=0}^{\infty} \frac{(-\Delta X)^n}{n!} \frac{\partial^n}{\partial X^n} \left[ G(X+\Delta X, t+\Delta t | X, t)P(X,t) \right] d(\Delta X)$$You can make an even sneakier improvement by moving the ##\Delta X## integration (and the ##(\Delta X)^n## term) inside the partial derivative,

$$P(X, t+ \Delta t) = \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} \frac{\partial^n}{\partial X^n} \left[ \int (\Delta X)^n G(X+\Delta X, t+\Delta t | X,t) d(\Delta X) \cdot P(X,t)\right]$$But what is ##\int (\Delta X)^n G(X+\Delta X, t+\Delta t | X,t) d(\Delta X)##? It’s the expectation ##E((\Delta X)^n)##. What you end up with is called the Kramers-Moyal expansion for the probability,

$$P(X, t+\Delta t) – P(X,t) = \sum_{n=1}^{\infty} \frac{(-1)^n}{n!} \frac{\partial^n}{\partial X^n} \left[ E((\Delta X)^n) P(X,t) \right]$$

You can turn this into a PDE by dividing through by ##\Delta t## and taking the limit ##\Delta t \rightarrow 0##.

$$\frac{\partial P}{\partial t} = \sum_{n=1}^{\infty} (-1)^n \frac{\partial^n}{\partial X^n} \left[ D^{(n)}(X) P(X,t) \right]$$

This is the famous Kramers-Moyal expansion with Kramers-Moyal coefficients:

$$D^{(n)}(X) := \lim_{\Delta t\rightarrow 0} \frac{1}{n!} \frac{E((\Delta X)^n) }{\Delta t}$$

The Fokker-Planck equation

If you cut off the Kramers-Moyal expansion at ##n=2##, then the result is the Fokker-Planck equation:

$$\frac{\partial P}{\partial t} = -\frac{\partial}{\partial X} [D^{(1)}(X) P] + \frac{\partial^2}{\partial X^2}[D^{(2)}(X) P]$$

Retaining only the first two terms works great a lot of the time. It can fail if there are pathological jump conditions.

Physicists are more familiar with ##\mu := D^{(1)}(x)## and ##\tfrac{1}{2} \sigma^2 = D = D^{(2)}(X)##, where ##D## is Einstein’s diffusion constant:

$$\frac{\partial P}{\partial t} = – \frac{\partial}{\partial X}(\mu P) + \frac{1}{2} \frac{\partial^2}{\partial X^2} (\sigma^2 P)$$

Example: Force-free Brownian motion

What do you notice in the case of zero drift, ##\mu = 0##? The Fokker-Planck equation becomes the diffusion equation

$$\frac{\partial P}{\partial t} = D \frac{\partial^2 P}{\partial X^2}$$

with probability current ##J = -D \partial P/\partial X##.

Example: Constant-force Brownian motion

In 1d, a typical (Langevin) equation of motion of a particle in a potential and subject to thermal noise is

$$\ddot{X} = -\gamma \dot{X} + F(X) + \sqrt{2kT\gamma} \frac{dW}{dt}$$

In this equation, ##\gamma## is the friction constant, ##F(X)## is the force field. The stochastic forcing term ##\sqrt{2kT\gamma} \tfrac{dW}{dt}## is sometimes called the Gaussian white noise. Restrict to the case where the particle is over-damped, ##|\ddot{X}| \ll \gamma |\dot{X}|##, so that you can all but ignore the ##\ddot{X}## term.

A mathematician would rather write the SDE:

$$dX = \frac{1}{\gamma} F(X) dt + \sigma dW$$

Simulation of diffusion under a constant force for 20 random paths.

For thermal stochastic noise, the volatility parameter ##\sigma = \sqrt{2kT/\gamma}## has a simple temperature dependence. The Fokker-Planck equation for this system is

$$\frac{\partial P}{\partial t} = -\frac{\partial}{\partial X} \left[ \frac{1}{\gamma} F(X) P \right] + \frac{1}{2} \sigma^2 \frac{\partial^2 P}{\partial X^2}$$

For a constant force ##F(X) = F_0##, the solution is analytic: it’s a Gaussian with mean position ##\bar{x}(t) := E(X(t)) = (F_0/\gamma)t## and a variance ##\sigma^2 t## increasing linearly with time:

$$P(X,t) = \frac{1}{\sqrt{2\pi \sigma^2 t}} \exp\left[ -\frac{(X-\bar{x}(t))^2}{2\sigma^2 t}\right]$$

Imagine what this probability density looks like as time proceeds. The peak moves to the right with speed ##F_0/\gamma##, and its amplitude gets squashed down asymptotically as it spreads out around the mean.

Ornstein-Uhlenbeck processes

So much for constant force fields. What about a harmonic force ##F(X) = -k(X-X_0)##? The corresponding SDE represents an Ornstein-Uhlenbeck process

$$dX = -\frac{k}{\gamma} (X-X_0) dt + \sigma dW$$

Simulation of an Ornstein-Uhlenbeck process with 20 random paths.

You can imagine a particle attached to a spring, being pulled back towards its equilibrium position ##X_0## while being disturbed by some stochastic noise. The process’s mean will approach the equilibrium position asymptotically, but the variance approaches a non-zero constant value ##\sigma^2/(2k)##. Physically: the stochastic noise prevents the system from ever settling!

You can find the full probability density ##P(X,t)## with a neat trick. Define new stochastic variable ##Y := e^{kt} X##. If you write down ##dY## using the Itô Lemma and insert it into the Ornstein-Uhlenbeck SDE, you get:

$$dY = kX_0 e^{kt} dt + e^{kt} \sigma dW$$

Now with the initial condition ##X(0) = 0##. After converting back to ##X##, $$X = X_0(1-e^{-kt}) + \sigma \int_0^t e^{k(t’-t)} dW(t’)$$The expectation is ##E(X(t)) := \bar{x}(t) = X_0(1-e^{-kt})##, because ##E(dW) = 0##. The variance is marginally more fiddly to expand out, but after dropping small terms you get

$$\mathrm{Var}(X) = E(X^2) – E(X)^2 = E \left[\left(  \sigma \int_0^t e^{k(t’-t)} dW(t’) \right)^2 \right]$$

It helps to re-write the right-hand side as a double integral$$\mathrm{Var}(X) = \sigma^2 \int_0^t \int_0^t e^{k(t’-t)} e^{k(t”-t)} E(dW(t’) dW(t”))$$The useful relationship for differential Brownian motions ##E(dW(t’) dW(t”)) = \delta(t’-t”)## means that this simplifies to

$$\mathrm{Var}(X) = \frac{\sigma^2}{2k}(1-e^{-2kt})$$In analogy to the result obtained for constant forcing, the full probability density in a harmonic potential is also a Gaussian, but with a very different time dependence. The function ##P(x,t)## is given by

$$P(x,t) = \frac{1}{\sqrt{\pi \frac{\sigma^2}{k}(1-e^{-2kt}) t}} \exp \left(- \frac{(X-\bar{x}(t))^2}{\frac{\sigma^2}{k}(1-e^{-2kt}) t}\right)$$with a central position ##E(X(t)) := \bar{x}(t) = X_0(1-e^{-kt})##.