Brownian Motion Simulation

At the moment of pricing options, the indisputable benchmark is the Black Scholes Merton (BSM) model presented in 1973 at the Journal of Political Economy. In the paper, they derive a mathematical formula to price options based on a stock that follows a Geometric Brownian Motion.

​

In this article, we will replicate that valuation, but instead of deriving the value of the option mathematically, we are going to generate 1000 paths through Monte Carlo simulation.

​

The libraries that we are going to use are:
- Math to transform the data with square root (x) and exponential (ex)
- Numpy to manipulate the data.
- Matplotlib to create the charts
 

import numpy as np                                                                                                                   

import math                                                                                                                                

import matplotlib.pyplot as plt                                                                                                 

​

The BSM model describes the market with one risky asset called the stock (S) and one riskless asset known as the bond (represented as a short term interest rate r). For this analysis, the environment that we are defining consist in an initial stock value S0=100  and a constant interest rate r=0.10. The volatility of the stock is represented by σ=0.2 and the time period for the simulation is one year t=1

​

S0 = 100.                                                                                                                                     

r = 0.1                                                                                                                                          

sigma = 0.20                                                                                                                               

T = 1                                                                                                                                            

​

For the simulation we will use 10,000 different paths (I=10000), and the year will be divided into 250 sub-periods (M=250).

​

M = 250                                                                                                                                      

I = 10000                                                                                                                                    

dt = T / M                                                                                                                                   

​

To create the different paths, we begin by utilizing the function np.random.standard_normal that draw  (M+1)×I samples from a standard Normal distribution. To ensure that the mean is 0 and the standard deviation is 1 we adjust the generated values with a technique called moment matching.

​

z = np.random.standard_normal((M+1,I))                                                                                

z -= z.mean()                                                                                                                               

z /= z.std()                                                                                                                                   

​

Once we have the normally distributed sample matrix, we can proceed to calculate the differents paths utilizing the difference equation of the Geometric Brownian Motion stochastic process. We substitute the values of the equation with the previously defined variables and replace the z term with the previously generated normal matrix.

​

S = np.zeros((M+1,I))                                                                                                                 

S [ 0 ] = S0                                                                                                                                  

for i in range (1, M+1):                                                                                                               

          S[ i ] = S[ i-1 ] * np.exp((r - sigma ** 2/2) * dt + sigma * math.sqrt(dt) * z[ i ])             

​

​

​

​