# Brownian Motion simulations 
# STAT 753
# April 21, 2020


#########################################################################
# Simulate a standard Brownian motion (i.e. sigma = 1) in 1 dimension
#########################################################################

# First generate N=1000 random normal variates
N = 1000
xs <- rnorm(N,0,1)
plot(xs, pch=19, col="blue", ylab="Value", main="Random Normal Variates")

# Plot density as a histogram
hist(xs, breaks=25, col="tomato", xlab= "Value", main="Histogram of 1000 random normal values")

# Generate a sample path of a standard Brownian motion
samplepath <- cumsum(xs)

# Need the process to start at 0 at time 0
samplepath <- c(0,cumsum(xs))

# Plot
plot(samplepath, type="l", xlab="Time", ylab="Value", main="Sample Path of Standard Brownian Motion in 1D")



###########################################
# Simulate Brownian motion with drift (1D)
###########################################

# Standard BM
tmax = 100
t = seq(from=0, to=tmax, by=0.01)  
x = rnorm(length(t))  
xt = cumsum(x) 
plot(xt, type="l", xlab="Time", ylab="Value", main="Sample Path of Standard Brownian Motion (1D)")


# BM with drift
sig = 3 
mu = 3
drift.xt = sig*xt + mu*t

# Plot sample path for BM with drift
plot(t, xt, type='l', ylim = c(min(drift.xt, xt), max(drift.xt, xt)), xlab="Time", ylab="Value",
      main="Sample Path of Brownian Motion with Drift (1D)")

# Optional - add linear trend line using linear regression function (lm)
lines(t, drift.xt, col = "red" )
abline(lm(drift.xt~t), col = "red" )
abline(lm(xt~t), col = "black" )
grid()
legend('topleft', c("BM with Drift", "Standard BM"), col=c("red","black"), lty=c(1,1))



####################################################
# Simulate standard Brownian motion in 2 dimensions
####################################################

# Generate N random normal variates for each dimension (now in a plane so for x and y)
N = 1000
xs <- rnorm(N,0,1)
ys <- rnorm(N,0,1)

# Cumulative sum, ensure the process starts at 0
xdis <- c(0,cumsum(xs))
ydis <- c(0,cumsum(ys))

# Plot
plot(xdis, ydis, type="l", col="blue", xlab="X-Value", ylab="Y-Value", 
     main="Sample Path of Standard Brownian Motion in 2D")

# Optional - color the starting point (origin) and ending point with different colors
points(xdis[1], ydis[1], pch=19, col="green")
points(xdis[length(xdis)], ydis[length(ydis)], pch=19, col="red")



########################################################
# Simulate 1D Brownian motion with reflecting boundary -- next time!
########################################################

# Set upper bound (ub) and lower bound (lb) for the boundary
ub = 10
lb = -10

# Generate N random normal variates
N = 1000
xs <- rnorm(N,0,1)

# Cumulative sum, start process at 0
samplepath <- c(0,cumsum(xs))

# Didn't finish this....