Imagine that a pollen particle is suspended in a glass of water. If we were to observe and record the vertical position of the particle over time, we would find that its movements were random. And if we were to plot this position, we’d get a jagged path through time (Figure $1$, left). This path would be just one of many possible paths, and if we were to repeat this observational experiment many times, we would not expect to see the same path again.

Given this randomness, how can we reason about this phenomenon? Can we say anything interesting or useful about the particle? For most of human history, this was a seemingly impossible task. A key insight, a conceptual pillar in probability theory, is to separate what actually happened (Figure $1$, left) from other possible outcomes (Figure $1$, right). This approach allows us to reason about the world through counterfactuals: what are all the possible paths the pollen could have taken? How likely is each path? What can we say about the distribution of outcomes?

This understanding of the pollen particle as a random process is a deep idea, and it took many decades and scientists to understand. The phenomenon was first observed in the 1830s by the Scottish botanist Robert Brown. Brown used a microscope to observe pollen particles suspended in water, and to his surprise, he saw the particles moving! At first, he thought this meant that the pollen particles were alive, but he tested and then rejected this hypothesis by observing the same effect with particles that he was convinced were inanimate, such as glass powder, minerals, and even pulverized fragments of the Egyptian Sphinx (Góra, 2006)! For roughly half-a-century, the phenomenon remained a mystery, although it became known as Brownian motion.

Then starting in 1905, Albert Einstein published a series of papers in which he hythesized that the pollen particles were moving because they were being bombarded by invisible molecules in the liquid (Einstein, 1905). In the following year, the Polish physicist Marian Smoluchowski independently published essentially the same theory (Von Smoluchowski, 1906). At the time, this theory was controversial, because the idea of molecules was not yet widely accepted. However, using statistical mechanics, Einstein and Smoluchowski were able to make testable predictions about the behavior of the particles, and another scientist, Jean Baptiste Perrin, verified the model a few years later (Perrin, 1909). And since Einstein’s breakthrough work, Brownian motion has been widely studied and more deeply understood. In the mathematical community, Brownian motion was formalized by Norbert Wiener (Wiener, 1923), and thus Brownian motion is often referred to as a Wiener process, particularly by mathematicians.

The goal of this post is to better understand Brownian motion. Brownian motion is an important concept because it can be used to model many phenomenon, from particles suspended in liquids to the prices of stocks. Ultimately, we’ll reconstruct the marginal distribution of our pollen particle at any given point in time. As we will see this, this is the normal distribution. This deep connection means that we can make mathematically precise probabilistic statements about a completely random process.

Random walks

Let’s begin with a simplified model of our pollen particle in discrete time. This is a stochastic process called a random walk. In the next section, we’ll extend this to continuous time, which is Brownian motion.

Imagine we can discretize time and then observe a single discrete “tick” on the clock. What happens to the pollen particle during this one tick? In our simple model of the world, we’re going to imagine that we flip a coin, not necessarily fair, and that the pollen particle moves up or down the same amount based on the outcome of that coin toss. The coin toss models the fact that the pollen particle is being randomly bombarded by water molecules and thus its position at the next time point is random. So the pollen particle cannot stay in place; after one tick of the clock, it moves up or down.

Formally, let $S_0$

Now consider the position $S_n$

$S_n = u Z_1 + u Z_2 + \dots + u Z_n. \tag{1}$

Since each $Z_i$

Clearly, as we repeatedly flip our coin, the set of possible locations of the pollen particle expands linearly with $n$. We can visualize all these possible locations as a directed graph or tree, sometimes called a binomial tree (Figure $3$, left)—we’ll explain the name in a moment. The tree layers (vertical slices) are zero-indexed, and so the root node occurs at time $n=0$

To help us identify nodes, let’s introduce the counting number $k$, which indexes the leaf nodes, taking values in $k \in \{0, 1, \dots, n\}$

Now that we understand this simple, discrete-time model for our pollen particle, let’s tackle our motivating question: which outcomes (leaf nodes) are most likely? Any given path is random, but can we say something about the distribution of outcomes?

To start, let’s compute the probability of arriving at the highlighted leaf node in Figure $4$. This is really the probability of arriving at a given node $(n, k)$

$\mathbb{P}\left(\{ \text{two heads and one tails} \}\right) = p^2 q. \tag{2}$

However, there are three ways flip two heads in three coin tosses,

$\{ HHT, HTH, THH \}, \tag{3}$

which is another way of saying that there are three paths to the highlighted node. Since each path is a mutually exclusive outcome, we compute our desired probability by summing the probability of all outcomes in Equation $2$ by the number of paths:

$\mathbb{P}\left(\{ \text{arriving at node $(3, 2)$} \}\right) = \mathbb{P}(K_3 = 2) = 3 p^2 q. \tag{4}$

For example, if $p=1/2$

To compute this probability in general, we just need a way to compute the number of ways to get $k$ successes or heads in $n$ trials. Since order matters, the number of ways to pick $k$ heads from $n$ coin tosses is

$n (n-1) (n-2) \dots (n-k+1) = \frac{n!}{(n-k)!}. \tag{5}$

First, we can choose any of $n$ coin tosses to be a heads. Then we can pick any of $n-1$

However, this overcounts the possible paths. For example, this does not distinguish between $H_1 H_2$

$\frac{\text{ordered ways to pick $k$ heads from $n$ tosses}}{\text{permutations of $k$ heads}} \;\; = \;\; \frac{n(n-1)(n-2) \dots (n-k+1)}{k(k-1)(k-2) \dots 1}. \tag{6}$

This number in Equation $6$ is often called the binomial coefficient, pronounced “$n$ choose $k$”, and is denoted as

${n \choose k} \triangleq \frac{n!}{k! (n-k)!} \tag{7}$

Putting it all together, the probability of arriving at the $k$-th node in the $n$-th layer of a binomial tree is

$\mathbb{P}(\{K_n = k\}) = {n \choose k} p^k q^{n-k}. \tag{8}$

The fact that these probabilities sum to one is just a trivial application of the binomial theorem. See A1.

Computing the mean and variance of $K_n$

$K_n = \frac{Z_1 + 1}{2} + \frac{Z_2 + 1}{2} + \dots + \frac{Z_n + 1}{2}. \tag{9}$

The mean of each $Z_i$

$\begin{aligned} \mathbb{E}[K_n] &= \sum_{i=1}^n \mathbb{E}\!\left[\frac{Z_i + 1}{2}\right] = np, \\ \mathbb{V}[K_n] &= \sum_{i=1}^{n} \mathbb{V}\!\left[\frac{Z_i + 1}{2}\right] = npq. \end{aligned} \tag{10}$

For the variance calculation, we use Bienaymé’s identity and the fact that $(Z_i+1)/2$

Of course, $K_n$

$\begin{aligned} S_n &= u \left( Z_1 + Z_2 + \dots + Z_n \right) \\ &= u \left(K_n - (n - K_n)\right) \\ &= u \left(2 K_n - n\right). \end{aligned} \tag{11}$

Since $n$ and $u$ are non-random, the events $\{S_n = u(2k-n)\}$

$\mathbb{P}(\{S_n = u(2k-n)\}) = \mathbb{P}(\{K_n = k\}) = {n \choose k} p^k q^{n-k}. \tag{12}$

So the location of our pollen particle by a given layer $n$ is determined by the distribution of a random variable $K_n$

$\begin{aligned} \mathbb{E}[S_n] &= \mathbb{E}[u(2K_n - n)] = 2nu(p-1/2), \\ \mathbb{V}[S_n] &= \mathbb{V}[u(2K_n - n)] = \mathbb{V}[2 u K_n] = 4 u^2 npq. \end{aligned} \tag{13}$

In the special case in which $p = 1/2$

$\begin{aligned} \mathbb{E}[S_n] &= 0, \\ \mathbb{V}[S_n] &= n u^2. \end{aligned} \tag{14}$

We can explore this distribution by plotting the function for various parameterizations (Figure $5$). Note that while $K_n$

And another way to visualize this is to imagine larger and larger binomial trees (Figure $6$). The distribution for the locations $S_n$

To summarize so far, we have done something remarkable. We have modeled the motion of a completely random particle, and yet we can say something concrete and precise about its distribution of locations over time.

To do this, however, we had to assume that time was discrete. So the natural next question is: what’s the distribution of our process in the continuous-time limit? At this point in our story, it is not far-fetched to guess that it’s the normal distribution. De Moivre proved the De Moivre–Laplace theorem, the earliest version of a central limit theorem (CLT), in 1738, so roughly a hundred years before Robert Brown observed Brownian motion. So scientists and mathematicians already knew that a sum of independent and identically distributed random variables converge to a Gaussian. The key insight in the development of Brownian motion was to realize that the bombardment of a pollen particle could be modeled as such as a sum.

Convergence of a rescaled random walk

So now let’s imagine what happens when the molecular bombardments on our pollen particle increase in number but decrease proportionally in impact. So we have more bombardments but they move the pollen particle less per bombardment. This rescaling is critical, or else the variance of our process would explode. Put in physical terms, if we increased the number of bombardments of our pollen particle but did not scale down the size of the move, the pollen particle’s moves would grow implausibly large.

To formalize this, let’s first fix $p = 1/2$

$B_t^{(n)} = u Z_1 + u Z_2 + \dots + u Z_{\lfloor tn \rfloor}, \quad u := \frac{1}{\sqrt{n}}. \tag{15}$

The notation ${\lfloor tn \rfloor}$ just indicates flooring to an integer since $t$ is a positive real number. And we need $u$ to scale with $n$, and so we set $u = 1 / \sqrt{n}$

$B_t = \lim_{n \rightarrow \infty} B_t^{(n)}. \tag{16}$

So again, we hold physical time $t$ fixed, and we make our binomial tree finer and finer (larger $n$ for fixed $t$). If we remove the grid of the binomial tree which clutters the visualization, and just visualize paths for finer and finer $n$, we can create visualizations similar to Figure $6$ but for much larger $n$ (Figure $7$).

Now we can ask the same question we asked in the discrete-time case: after physical time $t$, what is the distribution of our pollen particle’s position? As we observed above, it must be a normal distribution! Here, the insight is not that the binomial distribution converges to the normal distribution—again, this was known a hundred years before Robert Brown’s observations. The insight is that by modeling the continuous-time limit of a random walk as in Equation $15$, this rescaled random walk $B_t^{(n)}$

Let’s see this a bit more formally. The De Moivre–Laplace theorem states that a properly standardized binomial random variable converges to the normal distribution. In our notation, $K_n \sim \text{binom}(n, p)$

$\frac{K_n - \mathbb{E}[K_n]}{\sqrt{\mathbb{V}[K_n]}} = \frac{K_n - np}{\sqrt{npq}} = \frac{K_n - n/2}{\sqrt{n/4}} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(0, 1). \tag{17}$

Now observe that $B_t^{(n)}$

$\begin{aligned} B_t^{(n)} &= u Z_1 + u Z_2 + \dots + u Z_{\lfloor tn \rfloor} \\ &= \frac{1}{\sqrt{n}} \left( Z_1 + Z_2 + \dots + Z_{\lfloor tn \rfloor} \right) \\ &= \frac{1}{\sqrt{n}} \left(2 K_{\lfloor tn \rfloor} - {\lfloor tn \rfloor} \right) \\ &= \frac{1}{\sqrt{n}} \frac{\sqrt{\lfloor tn \rfloor}}{\sqrt{\lfloor tn \rfloor}} \; 2 \left(K_{\lfloor tn \rfloor} - {\lfloor tn \rfloor}/2 \right) \\ &= \sqrt{\frac{\lfloor tn \rfloor}{n}} \frac{K_{\lfloor tn \rfloor} - {\lfloor tn \rfloor}/2}{\sqrt{\lfloor tn \rfloor / 4}}. \end{aligned} \tag{18}$

By De Moivre–Laplace, we can say:

$\frac{K_{\lfloor tn \rfloor} - {\lfloor tn \rfloor}/2}{\sqrt{\lfloor tn \rfloor / 4}} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(0, 1). \tag{19}$

And as $n \rightarrow \infty$

$\sqrt{\frac{\lfloor tn \rfloor}{n}} \;\rightarrow\; \sqrt{t}. \tag{20}$

Since the standardized binomial converges in distribution to $\mathcal{N}(0, 1)$

$B_t^{(n)} \;\stackrel{d}{\rightarrow}\; \sqrt{t} \; \mathcal{N}(0, 1) = \mathcal{N}(0, t). \tag{21}$

That’s it! As an aside, I think that in a modern treatment, we would invoke Slutsky’s theorem to arrive at Equation $21$. Slutsky’s theorem states that if a sequence of random variables converges in distribution and is multiplied by a sequence converging to a constant, then the product converges in distribution to the constant times the limit.

The geometric interpretation of this is that the marginal distribution after time $t$ is simply the normal distribution $\mathcal{N}(0, t)$

Now that we see the simplest version of the derivation in its entirety, we can make two important adjustments. First, notice that our bombardment factor $u = 1/\sqrt{n}$

$u = \frac{\sigma}{\sqrt{n}}. \tag{22}$

It’s easy to see that this will flow through the derivation in Equation $18$ and give us

$\sigma B_t^{(n)} \;\stackrel{d}{\rightarrow}\; \sigma \sqrt{t} \; \mathcal{N}(0, 1) = \mathcal{N}(0, \sigma^2 t). \tag{23}$

But I think the more interesting adjustment is adding a drift parameter $\mu$. Of course, we could just shift our Brownian motion directly:

$\mu + \sigma B_t^{(n)} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(\mu, \sigma^2 t). \tag{24}$

But this has no physical meaning for our process. It’s just an arbitrary shift, not a drift. A richer way to approach this is to encode it directly into the bias of our coin flip. Intuitively, if we flip a biased coin (so $p \neq 1/2$

However, there’s a problem with this approach: since $p$ is constrained to $[0, 1]$

$\begin{aligned} \mathbb{E}[Z_i] &= 2p - 1, \\ \mathbb{E}\!\left[B_t^{(n)}\right] &= \frac{1}{\sqrt{n}} \lfloor tn \rfloor \mathbb{E}[Z_1]. \end{aligned} \tag{25}$

A more elegant approach is to make $p$ a function $\mu$. However, we cannot naively do this, since our drift could explode as $n \rightarrow \infty$

$p_n = \frac{1}{2} + \frac{\mu}{2 \sigma \sqrt{n}}. \tag{26}$

Intuitively, the factor $\mu /(2 \sigma \sqrt{n})$

$\mathbb{E}[Z_i] = 2 p_n - 1 = \frac{\mu}{\sigma \sqrt{n}}, \tag{27}$

and so the mean of our process—let’s denote it as $X_n^{(n)}$

$\mathbb{E}\!\left[X_t^{(n)}\right] = \frac{\sigma}{\sqrt{n}} \lfloor tn \rfloor \mathbb{E}[Z_1] = \frac{\sigma}{\sqrt{n}} \lfloor tn \rfloor \frac{\mu}{\sigma \sqrt{n}} \;\rightarrow\; \mu t. \tag{28}$

Putting these two adjustments together—one for the drift and one for the size of the bombardment—we can see that the general result is non-standard Brownian motion:

$X_t^{(n)} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(\mu t, \sigma^2 t). \tag{29}$

Alternatively, we could simply rewrite the main derivation (Equation $18$) using $u$ and $p_n$

Note that this isn’t a proof that the rescaled random walk converges to Brownian motion as a process. That requires more advanced mathematics such as Donsker’s theorem. Rather, it’s a claim about its marginal distribution at any fixed time $t$. But I think this provides amazing intuition for what Brownian motion really is without requiring much beyond elementary probability.

Conclusion

I still remember sitting in class for a course on probability and random process and watching the professor churn through the algebra to produce the insight in Equation $19$. It felt surprising and then obvious. The normal distribution is everywhere precisely because it is the limiting distribution for sums of independent and identically distributed random variables. We can shift or scale our random walk. We can make it asymmetric. It doesn’t really matter. We’ll still converge to a normal. And in my mind, this derivation builds good intuition for other properties of Brownian motion. For example, we can say that Brownian motion is a martingale or that it has stationary Gausian increments. The mathematics needed to make these claims precise might require some work, but the basic intuition is encoded in the derivations and visualizations above.

Appendix

A1. Binomial theorem

The binomial theorem is the following identity, which holds for any non-negative integer power $n$:

$(x + y)^n = \sum_{k=0}^n {n \choose k} x^k y^{n-k}. \tag{A1.1}$

This is easy to prove by induction. One can trivially check that the base case holds. And the inductive step is as follows:

$\begin{aligned} (x + y)^n (x + y) &= \sum_{k=0}^n {n \choose k} x^k y^{n-k} (x + y) \\ &= \sum_{k=0}^n {n \choose k} x^{k+1} y^{n-k} + \sum_{k=0}^n {n \choose k} x^k y^{n-k+1} \\ &:= A + B. \end{aligned} \tag{A1.2}$

If we write each sum $A$ and $B$ explicitly, it’s clear that we have $n-1$

$\begin{aligned} A &= {n \choose 0} x^1 y^n + {n \choose 1} x^2 y^{n-1} + \dots + {n \choose n-1} x^n y^1 + {n \choose n} x^{n+1} y^0, \\ \\ B &= {n \choose 0} x^0 y^{n+1} + {n \choose 1} x^1 y^n + \dots + {n \choose n-1} x^{n-1} y^2 + {n \choose n} x^n y^1. \end{aligned} \tag{A1.3}$

Collecting the $n-1$

$\begin{aligned} A+B &= \left[{n \choose 0} + {n \choose 1}\right] x^1 y^n + \dots + \left[{n \choose n-1} + {n \choose n}\right] x^n y^1 \\ &\quad + {n \choose 0} x^0 y^{n+1} + {n \choose n} x^{n+1} y^0. \end{aligned} \tag{A1.4}$

Finally, we can use the following identity to collapse bracketed binomial coefficients:

${n \choose k} = {n-1 \choose k-1} + {n-1 \choose k}. \tag{A1.5}$

And we can rewrite the non-overlapping terms in terms of $n+1$

${n \choose 0} = {n+1 \choose 0} = {n \choose n} = {n+1 \choose n+1} = 1. \tag{A1.6}$

This completes the inductive step:

$\begin{aligned} &(x+y)^{n+1} \\ &= {n+1 \choose 1} x^1 y^n + \dots + {n+1 \choose n} x^n y^1 + {n+1 \choose 0} x^0 y^{n+1} + {n+1 \choose n+1} x^{n+1} y^0 \\ &= \sum_{k=1}^{n+1} {n+1 \choose k} x^k y^{n+1-k}. \end{aligned} \tag{A1.7}$

Finally, the fact that the binomial distribution normalizes—discussed around Equation $8$—is simply a direct application of the binomial theorem for $x = p$

A2. Convergence with mean-centering and $\sigma$ scaling

Let $p_n$

$p_n = \frac{1}{2} + \frac{\mu}{2 \sigma \sqrt{n}}. \tag{A3.1}$

Then clearly

$\begin{aligned} \mathbb{E}[Z_i] &= 2 p_n - 1 = \frac{\mu}{\sigma \sqrt{n}}, \\\\ \mu_{K} := \mathbb{E}[K_{\lfloor tn \rfloor}] &= {\lfloor tn \rfloor} p_n \\ &= {\lfloor tn \rfloor} \left( \frac{1}{2} + \frac{\mu}{2 \sigma \sqrt{n}}\right), \\\\ \sigma_{K}^2 := \mathbb{V}[K_{\lfloor tn \rfloor}] &= {\lfloor tn \rfloor} p_n (1 - p_n) \\ &= {\lfloor tn \rfloor} \left( \frac{1}{2} + \frac{\mu}{2 \sigma \sqrt{n}}\right) \left( \frac{1}{2} - \frac{\mu}{2 \sigma \sqrt{n}}\right) \\ &= {\lfloor tn \rfloor} \left( \frac{1}{4} - \frac{\mu^2}{4 \sigma^2 n}\right). \end{aligned} \tag{A3.2}$

Let’s redefine $X_t^{(n)}$

$X_t^{(n)} = u Z_1 + u Z_2 + \dots + u Z_{\lfloor tn \rfloor}, \quad u := \frac{\sigma}{\sqrt{n}}. \tag{A3.3}$

We can write this as:

$\begin{aligned} X_t^{(n)} &= \frac{\sigma}{\sqrt{n}} \left[ Z_1 + Z_2 + \dots + Z_{\lfloor tn \rfloor} \right] \\ &= \frac{\sigma}{\sqrt{n}} \left[ 2 K_{\lfloor tn \rfloor} - {\lfloor tn \rfloor} \right] \\ &= \frac{\sigma}{\sqrt{n}} \; 2 \left[ K_{\lfloor tn \rfloor} - \lfloor tn \rfloor \frac{1}{2} \right] \\ &= \frac{\sigma}{\sqrt{n}} \; 2 \left[ K_{\lfloor tn \rfloor} - \lfloor tn \rfloor \frac{1}{2} - {\lfloor tn \rfloor} \left( \frac{\mu}{2 \sigma \sqrt{n}}\right) + {\lfloor tn \rfloor} \left( \frac{\mu}{2 \sigma \sqrt{n}}\right) \right] \\ &= \frac{\sigma}{\sqrt{n}} \; 2 \left[ K_{\lfloor tn \rfloor} - \lfloor tn \rfloor \left( \frac{1}{2} - \frac{\mu}{2 \sigma \sqrt{n}} \right) \right] + {\lfloor tn \rfloor} \frac{\mu}{n} \\ &= \frac{\sigma}{\sqrt{n}} \; 2 \left[ K_{\lfloor tn \rfloor} - \mu_K \right] + {\lfloor tn \rfloor} \frac{\mu}{n} \\ &= \frac{\sigma}{\sqrt{n}} \sqrt{4 \frac{\lfloor tn \rfloor}{\lfloor tn \rfloor} \frac{1 - \frac{\mu}{\sigma^2 n}}{1 - \frac{\mu}{\sigma^2 n}}} \left[ K_{\lfloor tn \rfloor} - \mu_K \right] + {\lfloor tn \rfloor} \frac{\mu}{n} \\ &= \frac{\sigma}{\sqrt{n}} \sqrt{\lfloor tn \rfloor \left( 1 - \frac{\mu}{\sigma^2 n} \right)} \left[ \frac{K_{\lfloor tn \rfloor} - \mu_K}{\sqrt{\frac{\lfloor tn \rfloor \left(1 - \frac{\mu}{\sigma^2 n}\right)}{4}}} \right] + {\lfloor tn \rfloor} \frac{\mu}{n} \\ &= \frac{\sigma}{\sqrt{n}} \sqrt{\lfloor tn \rfloor \left( 1 - \frac{\mu}{\sigma^2 n} \right)} \left[ \frac{K_{\lfloor tn \rfloor} - \mu_K}{\sigma_K} \right] + {\lfloor tn \rfloor} \frac{\mu}{n} \end{aligned} \tag{A3.4}$

Finally, it’s clear that the prefactor converges to $\sigma \sqrt{t}$

$\frac{\sigma}{\sqrt{n}} \sqrt{\lfloor tn \rfloor \left( 1 - \frac{\mu}{\sigma^2 n} \right)} \;\rightarrow\; \sigma \sqrt{t} \tag{A3.5}$

while the last term converges to $\mu t$ as $n \rightarrow \infty$

${\lfloor tn \rfloor} \frac{\mu}{n} \;\rightarrow\; \mu t. \tag{A3.6}$

And since

$\frac{K_{\lfloor tn \rfloor} - \mu_K}{\sigma_K} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(0, 1), \tag{A3.7}$

then again by Slutsky’s theorem, we know

$X_t^{(n)} \;\stackrel{d}{\rightarrow}\; \mathcal{N}(\mu t, \sigma^2 t). \tag{A3.8}$

So $X_t^{(n)}$