Simulating Brownian Motion in Python with Numpy

Sat 21 January 2017
In [1]:
import random
import math
import numpy as np
from functools import partial
from bokeh.io import show, output_notebook
from bokeh.layouts import row
from bokeh.plotting import figure
from bokeh.embed import notebook_div
import plotly.plotly as py
from plotly.graph_objs import *
random.seed(10)
output_notebook(hide_banner=True)
In [2]:
def brownian_path(N):
    Δt_sqrt = math.sqrt(1 / N)
    Z = np.random.randn(N)
    Z[0] = 0
    B = np.cumsum(Δt_sqrt * Z)
    return B
In [3]:
N = 365
T = [x for x in range(N)]
B = brownian_path(365)
p = figure(title='Sample path of a brownian motion')
r = p.line(T, B)
t = show(p)

Quadratic Variation and Total Variation of Brownian Motion

In [4]:
def ssq(B):
    """sum of squares"""
    ΔB = np.diff(B)
    ΔBsq = np.square(ΔB)
    ΔBqv = np.sum(ΔBsq)
    return ΔBqv
In [5]:
def ssd(B):
    """sum of abs"""
    ΔB = np.diff(B)
    ΔBabs = np.abs(ΔB)
    ΔBtv = np.sum(ΔBabs)
    return ΔBtv
In [6]:
N_max = 100000
N_seq = np.arange(100, N_max, step=100)
ΔBqv_seq = np.empty(N_seq.shape)
ΔBtv_seq = np.empty(N_seq.shape)
for i, n in enumerate(N_seq):
    B = brownian_path(n)
    ΔBqv_seq[i] = ssq(B)
    ΔBtv_seq[i] = ssd(B)
In [7]:
opts = dict(plot_width=450, plot_height=450, min_border=0)

p_qv = figure(**opts, title='Sum of squares sequence')
r_qv = p_qv.line(N_seq, ΔBqv_seq)

p_tv = figure(**opts, title='Sum of absolute changes sequence')
r_tv = p_tv.line(N_seq, ΔBtv_seq)

t_v = show(row(p_qv, p_tv))

Geometric Brownian Motion

In [8]:
gbm = lambda μ, σ, x0, t, bt: math.exp(np.log(x0) + (μ - 1 / 2 * σ ** 2) * t + σ * bt)
In [9]:
μ = 0.0003
σ = 0.025
x0 = 1
B = brownian_path(365)
GB = []
for t, bt in enumerate(B):
    gbt = gbm(μ, σ, x0, t, bt)
    GB.append(gbt)
pg = figure(title='Sample path of a geometric brownian motion')
rg = pg.line(range(len(GB)), GB)
tg = show(pg)

What do a brownian motion and geometric brownian motion with the same brownian sample path look like side by side?

In [10]:
μ = 1 / 2
σ = 1
x0 = 1
B = brownian_path(365)
GB = []
for t, bt in enumerate(B):
    gbt = gbm(μ, σ, x0, t, bt)
    GB.append(gbt)

opts = dict(plot_width=450, plot_height=450, min_border=0)
pg2 = figure(**opts, title='Geometric brownian motion')
rg2 = pg2.line(range(len(GB)), GB)

pb2 = figure(**opts, title='Brownian motion')
rb2 = pb2.line(range(len(B)), B)

tgb = show(row(pb2, pg2))