 # PQN #008: How to simulate stock prices with Python ## How to simulate stock prices with Python

In today’s issue, I’m going to show you how to simulate stock prices using Geometric Brownian Motion (GBM).

Simulating prices is fundamental for pricing derivatives. In the case of GBM, it is the key part of pricing equity options using Black-Scholes. Understanding how to simulate stock prices is foundational for any quantitative finance work. In fact, the math behind GBM was a good portion of my master’s degree.

The good news?

You don’t need a master’s degree to build your own stock price simulator in Python.

I’m going to show you how to do it.

Step by step.

## But first: A quiz

Which stock price chart is real?

(Answer at the bottom.)

## Step 1: What is Geometric Brownian Motion?

This is an abstract concept so I want to explain what GBM is all about.

Brownian motion comes from physics. It describes the random movement of particles in a substance. A Wiener process is a one-dimentional Brownian motion. It’s named after Norbert Wiener who won a Nobel Prize studying one-dimentional Brownian motions.

The Wiener process features prominently in quantitative finance because of some useful mathemetical properties.

The GBM is a continuous-time stochastic process where where the log of the random variable follows the Wiener process with drift.

What?

It’s a data series that trends up or down through time with a defined level of volatility.

And it’s perfect for simulating stock prices.

## Step 2: Import the libraries

I start by importing the libraries.

```import numpy as np
import matplotlib.pyplot as plt```

Since I’m simulating stock prices, I don’t need to download data.

## Step 3: Set the input parameters

To simulate stock prices, we need some input parameters.

```# setup params for brownian motion
s0 = 131.00
sigma = 0.25
mu = 0.35

# setup the simulation
paths = 1000
delta = 1.0/252.0
time = 252 * 5```

Start by defining the initial stock price, s0. Then sigma which is the percentage volatility. Finally, mu (drift), which determine the overall trend.

Setup the simulation with 1,000 simulation paths. delta refers to the time step. I want to model a new stock price every day so I use 1/252 to represent 1 day in 252 trading days. The length of the simulation is 5 years.

## Step 4: Build the functions

The first step is to build a function that returns a Wiener process. This function returns a 2-dimentional array with 1,260 rows and 1,000 columns. Each row is a day and each column is a simulation path.

```def wiener_process(delta, sigma, time, paths):
"""Returns a Wiener process

Parameters
----------
delta : float
The increment to downsample sigma
sigma : float
Percentage volatility
time : int
Number of samples to create
paths : int
Number of price simulations to create

Returns
-------
wiener_process : np.ndarray

Notes
-----
This method returns a Wiener process.
The Wiener process is also called Brownian
Wiener process check out the Wikipedia
page: http://en.wikipedia.org/wiki/Wiener_process
"""

# return an array of samples from a normal distribution
return sigma * np.random.normal(loc=0, scale=np.sqrt(delta), size=(time, paths))```

Next, I define a function that creates the GBM returns.

```def gbm_returns(delta, sigma, time, mu, paths):
"""Returns from a Geometric brownian motion

Parameters
----------
delta : float
The increment to downsample sigma
sigma : float
Percentage volatility
time : int
Number of samples to create
mu : float
Percentage drift
paths : int
Number of price simulations to create

Returns
-------
gbm_returns : np.ndarray

Notes
-----
This method constructs random Geometric Brownian
Motion (GBM).
"""
process = wiener_process(delta, sigma, time, paths)
return np.exp(
process + (mu - sigma**2 / 2) * delta
)```

Finally, I prepend a row of 1s to the returns array and multiply the starting stock price by the cumulative product of the GBM returns to produce the price paths.

```def gbm_levels(s0, delta, sigma, time, mu, paths):
"""Returns price paths starting at s0

Parameters
----------
s0 : float
The starting stock price
delta : float
The increment to downsample sigma
sigma : float
Percentage volatility
time : int
Number of samples to create
mu : float
Percentage drift
paths : int
Number of price simulations to create

Returns
-------
gbm_levels : np.ndarray
"""
returns = gbm_returns(delta, sigma, time, mu, paths)

stacked = np.vstack([np.ones(paths), returns])
return s0 * stacked.cumprod(axis=0)```

## Step 4: Visualize the results

I want to demonstrate two examples.

The first example simulates 1,000 price paths over 5 years. It starts at a price of 131 with 25% annualized volatility and 35% average return per year. This is the volatility and return of Apple in 2021.

```price_paths = gbm_levels(s0, delta, sigma, time, mu, paths)
plt.plot(price_paths, linewidth=0.25)
plt.show()```

As you might expect, a 35% drift causes most price paths to increase from from the initial price. In fact we can test exactly how many have increased from the initial price.

`len(price_paths[-1, price_paths[-1, :] > s0])`

999 of the 1,000 samples have increased.

In the second example I set the drift to 0.0.

```price_paths = gbm_levels(s0, delta, sigma, time, 0.0, paths)
plt.plot(price_paths, linewidth=0.25)
plt.show()```

And I get a much different picture.

Only 402 prices end up higher than the original price.

Spend some time and play around with the variables. What happens if you double volatility? What happens if you set mu to a negative number?

Quiz answer: The first stock chart is AAPL. The second is a simulation.

Well, that’s it for today. I hope you enjoyed it.

See you again next week.