
Disclaimer: The inspiration for this article stemmed from Georgia Tech’s Online Masters in Analytics (OMSA) program study material. I am proud to pursue this excellent Online MS program. You can also check the details here.
What is Brownian motion?
Physical origin and properties
Brownian motion, or pedesis, is the randomized motion of molecular-sized particles suspended in a fluid. It results from the stochastic collisions of the particles with the fast-moving molecules in the fluid (energized due to the internal thermal energy). The aforementioned fluid is supposed to be at the so-called thermal equilibrium, where no preferential direction of flow exists (as opposed to various transport phenomena).
Discovery and early explanations
Brownian motion is named after the Scottish botanist Robert Brown, who first described the phenomenon in 1827 while observing pollens (from the Clarkia pulchella plant) immersed in water, through a microscope.
For a long time, the scientific community did not think much of it. Then, in 1905, a 26-year old Swiss patent clerk changed the world of physics by analyzing the phenomena with the help of the laws of thermodynamics.

Albert Einstein published a seminal paper where he modeled the motion of the pollen, influenced by individual water molecules, and depending on the thermal energy of the fluid. Although this paper is relatively less celebrated than his other 1905 papers, it is one of his most cited publications. In fact, Einstein’s explanation of Brownian motion served as the first mathematically sound evidence that molecules exist.
You can read this enjoyable article commemorating the 100-year of Einstein’s paper.
Brownian motion as a stochastic process
Many-body interaction
The many-body interactions, that yield the intricate yet beautiful pattern of Brownian motion, cannot be solved by a first-principle model that accounts for the detailed motion of the molecules. Consequently, only probabilistic macro-models applied to molecular populations can be employed to describe it.
This is the reasoning behind the description of Brownian motion mostly as a purely stochastic process in its modern form. Almost all practical application also adopts this approach.
Wiener process
Mathematical properties of the one-dimensional Brownian motion was first analyzed American mathematician Norbert Wiener. The resulting formalism is a real-valued continuous-time stochastic process, called the Wiener process
It is one of the best known stochastic processes with attractive properties like stationarity and independent increments. Consequently, it finds frequent applications in a wide range of fields covering pure and applied mathematics, quantitative finance, economic modeling, quantum physics, and even evolutionary biology.
Common applications
The concept and mathematical models of Brownian motion play a vital role in modern mathematics such as stochastic calculus and diffusion processes. The Wiener process is also used to represent the integral of a white noise Gaussian process, which, often acts as a ubiquitous model of noise in electrical and electronics engineering. The same idea is also applied to model noise and errors in filter design and random/unknown forces in control theory.

In quantum physics, diffusion phenomena related to the Fokker-Planck and Langevin equations are studied with the help of Brownian motion. It also underlies the formation of the rigorous path integral formulation of quantum mechanics. For example, using the Feynman-Kac formula, a solution to the famous Schrodinger equation can be represented in terms of the Wiener process. The model of eternal inflation in physical cosmology takes inspiration from the Brownian motion dynamics.
In the world of finance and econometric modeling, Brownian motion holds a mythical status. It features prominently in almost all major mathematical theories of finance. In particular, the famous Black-Scholes option pricing model, for which Myron Scholes received 1997 Nobel prize in economics, depends on this formalism.

Python implementation
A rather simple equation
The core equation at the heart of generating data points following a Brownian motion dynamics is rather simple,

where Yi could be a basic stochastic process like Random Walk or sample from a Normal distribution.
A Brownian
class
The Jupyter notebook for the implementation can be found here.
In the Python code below, we define a class Brownian
with a few useful methods,
gen_random_walk()
: Generates motion from the Random Walk processgen_normal()
: Generates motion by drawing from the Normal Distributionstock_price()
: Models a stock price using the so-called ‘Geometric Brownian Motion’ formula
class Brownian():
"""
A Brownian motion class constructor
"""
def __init__(self,x0=0):
"""
Init class
"""
assert (type(x0)==float or type(x0)==int or x0 is None), "Expect a float or None for the initial value"
self.x0 = float(x0)
def gen_random_walk(self,n_step=100):
"""
Generate motion by random walk
Arguments:
n_step: Number of steps
Returns:
A NumPy array with `n_steps` points
"""
# Warning about the small number of steps
if n_step < 30:
print("WARNING! The number of steps is small. It may not generate a good stochastic process sequence!")
w = np.ones(n_step)*self.x0
for i in range(1,n_step):
# Sampling from the Normal distribution with probability 1/2
yi = np.random.choice([1,-1])
# Weiner process
w[i] = w[i-1]+(yi/np.sqrt(n_step))
return w
def gen_normal(self,n_step=100):
"""
Generate motion by drawing from the Normal distribution
Arguments:
n_step: Number of steps
Returns:
A NumPy array with `n_steps` points
"""
if n_step < 30:
print("WARNING! The number of steps is small. It may not generate a good stochastic process sequence!")
w = np.ones(n_step)*self.x0
for i in range(1,n_step):
# Sampling from the Normal distribution
yi = np.random.normal()
# Weiner process
w[i] = w[i-1]+(yi/np.sqrt(n_step))
return w
def stock_price(
self,
s0=100,
mu=0.2,
sigma=0.68,
deltaT=52,
dt=0.1
):
"""
Models a stock price S(t) using the Weiner process W(t) as
`S(t) = S(0).exp{(mu-(sigma^2/2).t)+sigma.W(t)}`
Arguments:
s0: Iniital stock price, default 100
mu: 'Drift' of the stock (upwards or downwards), default 1
sigma: 'Volatility' of the stock, default 1
deltaT: The time period for which the future prices are computed, default 52 (as in 52 weeks)
dt (optional): The granularity of the time-period, default 0.1
Returns:
s: A NumPy array with the simulated stock prices over the time-period deltaT
"""
n_step = int(deltaT/dt)
time_vector = np.linspace(0,deltaT,num=n_step)
# Stock variation
stock_var = (mu-(sigma**2/2))*time_vector
# Forcefully set the initial value to zero for the stock price simulation
self.x0=0
# Weiner process (calls the `gen_normal` method)
weiner_process = sigma*self.gen_normal(n_step)
# Add two time series, take exponent, and multiply by the initial stock price
s = s0*(np.exp(stock_var+weiner_process))
return s
Process with an initial value of zero and using random walk
We can use a basic stochastic process such as Random Walk, to generate the data points for Brownian motion.

b = Brownian()
for i in range(4):
plt.plot(b.gen_random_walk(1000))
plt.show()

Process with an initial value of 20 and using Normal distribution
We can generate Brownian motion data by drawing from Normal distribution.
b = Brownian(20)
for i in range(4):
plt.plot(b.gen_normal(1000))
plt.show()

Stock price simulation
We implemented the Geometric Brownian Motion model in the class as a method.

In the demo, we simulate multiple scenarios with for 52 time periods (imagining 52 weeks a year). Note, all the stock prices start at the same point but evolve randomly along different trajectories.
Note that the dynamics are controlled by the mean and variance parameters of the underlying Normal distribution. This somehow emulates the growth trend and the ‘volatility’ of the stock.
For example, a stock with a positive growth trend will have a positive mean. For this particular simulation, the choice of the mean (mu) is 0.2 and the choice of standard deviation (square root of the variance) is 0.68.
We define a utility function for plotting first.
def plot_stock_price(mu,sigma):
"""
Plots stock price for multiple scenarios
"""
plt.figure(figsize=(9,4))
for i in range(5):
plt.plot(b.stock_price(mu=mu,
sigma=sigma,
dt=0.1))
plt.legend(['Scenario-'+str(i) for i in range(1,6)],
loc='upper left')
plt.hlines(y=100,xmin=0,xmax=520,
linestyle='--',color='k')
plt.show()
Note that, although the scenarios look sufficiently stochastic, they have a downward trend. This is because even with a positive mean, we have a slightly high spread or volatility.
plot_stock_price(mu=0.2,sigma=0.7)

A slight change in the volatility
We simulate the stock price again with slightly less volatility (but with the same mean as before) and get a completely different outcome this time. The trend looks neutral i.e. some scenario shows an increase in the stock price after the 52-week time period, whereas some show a decrease.
plot_stock_price(mu=0.2,sigma=0.65)

Lowering the volatility further down
If we lower the volatility, even more, we get a clear positive trend.
plot_stock_price(mu=0.2,sigma=0.6)

Two-dimensional plot
In the following example, we show a two-dimensional Brownian motion much like the actually suspended particle in the fluid medium goes through.
b1 = Brownian()
b2 = Brownian()
x = b1.gen_normal(1000)
y = b2.gen_normal(1000)
plt.plot(x,y,c='b')
xmax,xmin,ymax,ymin = x.max(),x.min(),y.max(),y.min()
scale_factor = 1.25
xmax,xmin,ymax,ymin = xmax*scale_factor,xmin*scale_factor,ymax*scale_factor,ymin*scale_factor
plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
plt.show()

Summary
We showed how to generate random datasets corresponding to the Brownian motion in one and two dimensions. We also showed an application of the idea in stock price simulation using the Geometric Brownian motion model.
Having a ready-made Python implementation for this important stochastic process is extremely important because of its ubiquitousness in various real-life applications.
For example, Data Science practitioners can readily take this implementation and integrate it with their model of a stochastic process when they are analyzing a quantitative finance or physics model.
Again, the Jupyter notebook for the implementation can be found here. Readers are welcome to fork it and extend as per their requirements.
Also, you can check the author’s GitHub repositories for code, ideas, and resources in Machine Learning and data science. If you are, like me, passionate about AI/machine learning/data science, please feel free to add me on LinkedIn or follow me on Twitter.
Tirthajyoti Sarkar – Sr. Principal Engineer – Semiconductor, AI, Machine Learning – ON…