Riding the waves of randomness: using Brownian motion for stock price simulations in Python
Introduction
As an ubiquitous phenomenon in nature, Brownian motion refers to the random motion of particles within a medium, usually a gas or a liquid. Having been discovered by the botanist Robert Brown in 1827 through a microscopy analysis on pollen seeds immersed in water, it took about further eighty years to establish a profound quantitative understanding of this phenomenon. Among other scientists, noteworthy are hereby the contributions by Albert Einstein that provided a theoretical foundation of Brownian motion associated with diffusion [1, 2].
In simple terms, the idea behind Brownian motion is that a trajectory follows a random and unpredictable path. This idea has been widely utilized across several scientific and application-oriented disciplines, among which also includes quantitative finance for describing the random behavior of stock price movements. In this story, we are going to use Python in order to carry out stock price simulations based on Brownian motion.
Disclaimer: This article is only for educational purposes and shall not be taken as an investment, financial, or other professional advice.
Standard Brownian motion (Wiener process)
The stochastic nature of stock prices relies on the empirical observation that daily returns nearly follow a normal distribution. The driving force behind this part is the standard Brownian motion, also known as Wiener process. Specifically, for an increment {W(t) : t ≥ 0} accounting for a change in stock price at time t, the Wiener process exhibits the following characteristics:
- It starts at zero, i.e., W(0) = 0
- For any timestep dt, the difference between the subsequent and current increment follows a normal distribution: W(t+dt) - W(t) ~ N(0, dt)
- The evolution of W(t) is continuous
An example for implementing W(t) in Python is given in the script below, along with one plot of a possible evolution
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
#Function for calculating the temporal evolution of W
#for n=1000 timesteps
def wiener_process(dt=0.1, x0=0, n=1000):
# Initialize W(t)
W = np.zeros(n+1)
# Create n+1 timesteps: t=0,1,2,3...n
t = np.linspace(x0, n, n+1)
# Use the cumulative sum for calculating W at every timestep
W[1:n+1] = np.cumsum(npr.normal(0, np.sqrt(dt), n))
return t, W
# Function for plotting W
def plot_process(t, W):
plt.plot(t, W)
plt.xlabel('Time(t)')
plt.ylabel('W(t)')
plt.title('Wiener process')
plt.show()
[t, W] = wiener_process()
plot_process(t, W)
Stock price behavior described by arithmetic Brownian motion
Given the characteristics of the Wiener process stated above, it can theoretically represent the daily stock price fluctuations without considering any bullish or bearish behavior (i.e., growing or declining trend in price). Because, however, such case does not occur for actual stocks, we also need to take into account a deterministic component, i.e., a drift, that is positive (negative) in case of a bullish (bearish) trend.
With this in mind, the stock price, denoted S(t) at time t, experiences a change dSₜ during a timestep dt that is given by
which represents an arithmetic Brownian motion. From eq. (1), we see that two fundamental quantities enter in the expressions of the deterministic and stochastic components:
- The average return μ that is responsible for the drift
- The risk/volatility σ that accounts for the stochastic behavior associated with the Wiener process
A qualitative representation of S(t) given by eq. (1) is shown in Figure 2
A drawback with eq. (1), however, is that stock prices within this model can generally attain negative values. Because in reality this never happens, we need to encounter a different version of eq. (1), in which the stock is not subjected to a normal, but a log-normal distribution instead. This is achieved using geometric Brownian motion.
Stock price simulation with Geometric Brownian motion
The geometric Brownian motion is characterized by
that has the following solution:
In this form, since S(t) follows a log-normal distribution, the non-negativity of the stock price is always met. We now can make use of eq. (3) to simulate the evolution of a stock price with the following parameters
- S(0) = S₀ = 10 $
- μ = 0,1
- σ = 0,05
within the timespan [0, T] (using T = 1 by default), based on the below .py script:
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
def simulate_geometric_brownian_motion(S0, T=1, N=1000, mu=0.1, sigma=0.05):
#timestep dt
dt = T/N
# Initialize W(t)
W = np.zeros(N+1)
# Create N timesteps
t = np.linspace(0, T, N+1)
# Use the cumulative sum for calculating W at every timestep
W[1:N+1] = np.cumsum(npr.normal(0, np.sqrt(dt), N))
#Calculate S(t)
S = S0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * W)
return t, S
def plot_simulation(t, S):
plt.plot(t, S)
plt.xlabel('Time (t)')
plt.ylabel('Stock Price S(t)')
plt.title('Geometric Brownian Motion')
plt.show()
#Calculate S with an initial price S0 = 10 $
[t, S] = simulate_geometric_brownian_motion(10)
plot_simulation(t, S)
In order to infer whether S(t) truly follows a log-normal distribution, one would need to conduct many trajectories of geometric Brownian motions, and then plot the distributions of S(t). This can be verified using Aleatory, one powerful Python library for simulating and visualising stochastic processes [3].
Simulating Geometric Brownian motion with Aleatory
The simulation is pretty straightforward, since all the process conducted in the above .py script, is implemented in Aleatory within a class called GBM. Using any development environment, e.g., Jupyter Notebook, an instance of GBM can be created as
# Install aleatory, if necessary
pip install aleatory
# Import the geometric Brownian motion class GBM
from aleatory.processes import GBM
# Create an instance of GBM
geometric_brownian_motion = GBM(drift=1, volatility=0.5, initial=1.0, T=1.0, rng=None)with the following parameters (using again T = 1 by default)
- S₀ = 1 $
- μ = 1
- σ = 0.5
The last entry refers to the type of random number generator, which is given as None by default. Given this, the plot can be drawn by using the draw method,
geometric_brownian_motion.draw(n=100, N=200)where n corresponds to the number of steps in each trajectory, and N to the number of trajectories to simulate.
As shown in Figure 4, the simulation yields several stock price trajectories, with a log-normal distribution, as indicated by its distinct skewness.
Conclusion
In summary, using Brownian motion to simulate stock prices can provide a deeper insight into market dynamics, as well as into potential future trends of stocks. While Python is equipped with the extensive Aleatory library to carry out a variety of such simulations, it is also important to note that Brownian motion offers a simplified representation of stock movements, and that it does not capture other factors that can have a critical impact under real world market conditions (economic events, interest rates, etc.). Nevertheless, when complemented with further analyses and empirical data, Brownian motion is the tool of choice for informed decision-making in the financial market.
Sources
[1] A. Einstein, Ann. Phys. 17, 549 (1905).
[2] A. Einstein, Zeitschrift fuer Elektrochemie und angewandte physikalische Chemie 13, 41 (1907).
