About me

Introduction

Brownian motion is a mathematical description of the random motion of a “large” particle immersed in a fluid and which is not subject to any other interaction than shocks with the “small” molecules of the surrounding fluid. However, this theory is only used to describe a random motion of a large particle. The latter is very useful in finance, for example the famous Black-Scholes model which evaluates the price of an asset over continuous time. Our report consists of five parts, the first one is reserved to the historical side of the Brownian motion, in the second part we will define mathematically the Brownian motion, as well as its construction in the third part with the help of a random walk or by a Gaussian process. Fourthly, we will see some properties of Brownian motion. And finally the arithmetic and geometric Brownian motion.

Importing libraries

The following R packages are required for this study.

## Load Packages
# Importing the required packageS
library(Metrics)
library(zoo,warn.conflicts=FALSE)
library(lubridate,warn.conflicts=FALSE)
library(mgcv,warn.conflicts=FALSE)
library(rugarch,warn.conflicts=FALSE)
# visualization
suppressPackageStartupMessages(library(ggplot2)) 
# calculating returns
suppressPackageStartupMessages(library(PerformanceAnalytics)) 
library(tidyverse)  # ggplot(), %>%, mutate(), and friends
library(tidyr) #For ggplot
library(broom)  # Convert models to data frames
library(scales)  # Format numbers with functions like comma(), percent(), and dollar()
library(modelsummary)  # Create side-by-side regression tables
# To read csv and xlsx files 
library(readr)
library(here)
library(readxl)
#library(nlme)
#Transform data to time series
library(zoo)
library(dplyr)
library(ggrepel)
library(ggpubr)
library(forecast)
library(pastecs)
#For Maximum Likelihood method
library(bbmle)
#For Garch model
library(fBasics)
library(fGarch)

Construction of the Brownian Motion

Construction by a Gaussian Process

La fonction BMN permet de simuler un mouvement brownien standard \(\{B_t, t \geq 0 \}\) dans l’intervalle de temps \([t_0, T]\) avec un pas \(\Delta t=\frac{T-t_0}{N}\).

library(Sim.DiffProc)

Brownian motion is very easy to simulate. To start off, let’s simulate a single instance of Brownian motion for 100 generations of discrete time in which the variance of the diffusion process is \(\sigma^2 = 0.01\) per generation.

t <- 0:100  # time
sig2 <- 0.01
## first, simulate a set of random deviates
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
## now compute their cumulative sum
x <- c(0, cumsum(x))
plot(t, x, type = "l", ylim = c(-2, 2))

The trick here is that we draw random normal deviates for the change over each t time intervals; and then to get the state of our chain at each time interval we just compute the cumulative sum of all the individual changes using `cumsum. We can do this because the distribution of the changes under Brownian motion is invariant and does not depend on the state of the chain.

We can also easily do a whole bunch of simulations like this at once, using the same conditions:

nsim <- 100
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2)), nsim, length(t) - 1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
plot(t, X[1, ], xlab = "time", ylab = "phenotype", ylim = c(-2, 2), type = "l")
apply(X[2:nsim, ], 1, function(x, t) lines(t, x), t = t)

## NULL

To see how the outcome depends on \(\sigma^2\), let’s compate the result when we divide sig2 by 10:

X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2/10)), nsim, length(t) - 
    1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
plot(t, X[1, ], xlab = "time", ylab = "phenotype", ylim = c(-2, 2), type = "l")
apply(X[2:nsim, ], 1, function(x, t) lines(t, x), t = t)

## NULL

Cool.

There are a number of different ways we could’ve done this. For example, instead of using apply, which economizes on code, we could have used a for loop as follows:

X <- matrix(0, nsim, length(t))
for (i in 1:nsim) X[i, ] <- c(0, cumsum(rnorm(n = length(t) - 1, sd = sqrt(sig2))))
plot(t, X[1, ], xlab = "time", ylab = "phenotype", ylim = c(-2, 2), type = "l")
for (i in 1:nsim) lines(t, X[i, ])

As mentioned above, the expected variance under Brownian motion is just \(\sigma^2 \times t\). To see this easiest, we can just do the following. Here I will use 10,000 simulations for 100 generations under the same conditions to “smooth out” our result:

nsim <- 10000
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2)), nsim, length(t) - 
    1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
v <- apply(X, 2, var)
plot(t, v, type = "l", xlab = "time", ylab = "variance among simulations")

var(X[, length(t)])  # this should be about 1.00
## [1] 1.008685
par(mfrow=c(1,1))
plot(BM(N =1000,M=1,x0=0,t0=0,T=1),main="Brownian trajectory")

X <-BM(N =10000,M=100,x0=0,t0=0,T=10)
plot(X,plot.type="single",main="Flow of Brownian trajectories")
lines(as.vector(time(X)),rowMeans(X),col="red")

Simple Brownian Motion Code in One Dimension

N = 1000
dis = rnorm(N, 0, 1);
dis = cumsum(dis);
plot(dis, type= "l",main= "Brownian Motion in One Dimension", xlab="time", ylab="displacement")

Simple Brownian Motion Code in Two

N = 1000; xdis = rnorm(N, 0 ,1);
ydis = rnorm(N, 0 ,1);
xdis = cumsum(xdis);
ydis = cumsum(ydis);
plot(xdis, ydis, type="l", main ="Brownian Motion in Two Dimension", xlab="x displacement", 
     ylab = "y displacement");

Brownian Motion with Displacement Code in One Dimension

N = 1000;
dis = rnorm(N, 0, 1);
at = rpois(N,1)
for(i in 1:N){
if(at[i] != 0){
dis[i]= dis[i]*at[i];
}}
dis = cumsum(dis)
plot(dis, type= "l",main= "Brownian Motion in One Dimension with Poisson Arrival Process",
     xlab="time", ylab="displacement")

Simulation using Brownian motion

Now let’s try to simulate using Brownian motion up the branches of a tree. Sticking with discrete time, we first need a discrete time phylogeny, which we can simulate using pbtree in the “phytools” package:

library(phytools)
t <- 100  # total time
n <- 30  # total taxa
b <- (log(n) - log(2))/t
tree <- pbtree(b = b, n = n, t = t, type = "discrete")
## simulating with both taxa-stop (n) and time-stop (t) is
## performed via rejection sampling & may be slow
## 
##   101 trees rejected before finding a tree
plotTree(tree)

Now, to simulate on all the branches of the tree we just need to simulate on each branch separately, and then “shift” each daughter branch by the end state of it’s parent. We can do this, remember, because the outcome of Brownian evolution at each time step is independent of all other time steps.

## simulate evolution along each edge
X <- lapply(tree$edge.length, function(x) c(0, cumsum(rnorm(n = x, sd = sqrt(sig2)))))
## reorder the edges of the tree for pre-order traversal
cw <- reorder(tree)
## now simulate on the tree
ll <- tree$edge.length + 1
for (i in 1:nrow(cw$edge)) {
    pp <- which(cw$edge[, 2] == cw$edge[i, 1])
    if (length(pp) > 0) 
        X[[i]] <- X[[i]] + X[[pp]][ll[pp]] else X[[i]] <- X[[i]] + X[[1]][1]
}
## get the starting and ending points of each edge for plotting
H <- nodeHeights(tree)
## plot the simulation
plot(H[1, 1], X[[1]][1], ylim = range(X), xlim = range(H), xlab = "time", ylab = "phenotype")
for (i in 1:length(X)) lines(H[i, 1]:H[i, 2], X[[i]])
## add tip labels if desired
yy <- sapply(1:length(tree$tip.label), function(x, y) which(x == y), y = tree$edge[, 
    2])
yy <- sapply(yy, function(x, y) y[[x]][length(y[[x]])], y = X)
text(x = max(H), y = yy, tree$tip.label)

There are probably many more economical ways to do that, but this is one.

In reality, most simulations of Brownian motion are conducted using continuous rather than discrete time. This is because the additivity of Brownian motion means that the expected variances among & covariances between species are the same in whether we simulate \(t\) steps each with variance \(\sigma^2\), or one big step with variance \(\sigma^2_t\).

Here is an example of a continuous time simulation & visualization using canned functions.

## simulate Brownian evolution on a tree with fastBM
x <- fastBM(tree, sig2 = sig2, internal = TRUE)
## visualize Brownian evolution on a tree
phenogram(tree, x, spread.labels = TRUE, spread.cost = c(1, 0))

Construction by a limit of a random walk

n <-10 ;T <-1;
t <-seq (0,T, length =100)
S <-cumsum (2*( runif (n ) >0.5) -1)
W <-sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0))
W <-as.numeric (W)/sqrt (n)
plot (t,W, type ="l",xlab=expression(W[t]),ylim=c(-2,2),las=1)
lines (t,W, lty =3)
n <-100
S <-cumsum (2*( runif (n ) >0.5) -1)
W <-sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0))
W <-as.numeric (W)/sqrt (n)
lines (t,W, lty =2)
n <-1000
S <-cumsum (2*( runif (n ) >0.5) -1)
W <-sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0))
W <-as.numeric (W)/sqrt (n)
lines (t,W, lty =3)
legend("topleft",c("n=10","n=100","n=1000"),lty=c(1,2,3),lwd=1,cex=0.9)

Brownian trajectory as a limit of a random walk.

Arithmetic Brownian motion

par(mfrow=c(1,1))
X<-ABM(M=1,N=1000,theta=1,sigma=1)
plot(X,main="Trajectory of a Brownian motion Arithmetic")

FL<-ABM(N=1000,M = 10,theta=1,sigma=1)
plot(FL,plot.type="single",main="Flow of Brownian Motion Arithmetic ")
lines(as.vector(time(FL)),rowMeans(FL),col="red")

Geometric Brownian motion

par(mfrow=c(1,1))
X<-GBM(M=1,N=1000,theta=1,sigma=1)
plot(X,main="Trajectory of a Geometric Brownian Motion")

FLG<-GBM(N=1000,M = 10,theta=1,sigma=1)
plot(FLG,plot.type="single",main="Flow of Geometric Brownian Motion ")
lines(as.vector(time(FLG)),rowMeans(FLG),col="red")

N = 1000;
dis = rnorm(N, 0 ,1);
dis = cumsum(dis);
plot(dis, type = "l", main="Pseudo-fractal simulation with Brownian motion",xlim=c(0,1000),
     ylim=c(- 70,70))
fnorm = matrix(nrow= 5, ncol = N)
for(i in 1:5){
for(j in 1:N){
fnorm[i,j] = rnorm(1,0,1)
}}
node = c(100,300,500,700,900);
fdis = matrix(nrow= 5, ncol = N)
for(i in 1:5){
fdis[i,] = cumsum(fnorm[i,]) + dis[node[i]]
}
for(i in 1:5){
lines(node[i]:N,fdis[i,1:(N-node[i]+1)])
}

BMF = function(start, end, total,var){ N = end - start + 1;
dis = rep(NA, total);
replace = BM(N,var)
for(i in start:end){
dis[i] = replace[i + 1 - start];
}
return(dis)
}

total = 1000;
xdis = BMF(1,1000,total,1)
plot(xdis, type = "l", main="Pseudo-Fractal Simulation with Brownian motion with 2 iteration",
     ylim = c(-100, 100), xlab = "steps",ylab="scale", lwd= 1.4)

References