# SDE simulations
# STAT 753
# April 23, 2020 (updated 4-29-20)

###########################################################
# Simulate various stochastic differential equations (SDEs)
###########################################################

# Install the 'sde' package, then call the sde library
library(sde)


#########################################################################################
# Built in functions for Brownian motion, Brownian bridge, and Geometric Brownian motion
#########################################################################################

######### Parameters used in functions below ############
# x - initial value of the process (default is 0)
# y - can specify terminal value of process, for instance 0 at time 1 for Brownian bridge 
# t0 - initial time (default is 0)
# T - final time
# N - number of intervals in which to split time vector [t0,T]
# r - interest rate of GBM
# sigma - volatility of GBM


# Brownian motion
B1 <- BM(x=0, t0=0, T=1, N=1000)
plot(B1, main = "Brownian Motion", ylab="Value")


# Brownian bridge
N=100
#N=1000
B2 <- BBridge(x=0, y=0, t0=0, T=1, N=N)
plot(B2, main = "Brownian Bridge", ylab="Value")
abline(0, 0, col="red")


# Geometric Brownian motion
B3 <- GBM(x=1, r=1, sigma=1, T=1, N=1000)
plot(B3, main = "Geometric Brownian Motion", ylab="Value")



#############################
# Brownian motion with drift 
#############################

# Default parameter values
# t0 = 0      - initial time
# T = 1       - final time
# X0 = 1      - initial state
# N = 100     - number of simulation steps
# M = 1       - number of trajectories
# delta = *   - time step of simulation *(if unspecified, it's (T-t0)/N; if specified, then T=N*delta)

# 1 sample path (trajectory) of BM with drift
set.seed(123)
d <- expression(5)  #drift term: mu = 5
s <- expression(3.5)   #scaling term: sigma = 3.5
sx <- expression(0)
BM <- sde.sim(X0=0, drift=d, sigma=s, sigma.x=sx)   #method="euler" by default
plot(BM, main="Brownian Motion with Drift Sample Path", col="forestgreen", lwd=2)

# Multiple trajectories of the OU process
set.seed(123)
BM2 <- sde.sim(X0=0, drift=d, sigma=s, M=3)
plot(BM2, main="Multiple trajectories of BM with Drift", col="orange2", lwd=2)


# Set final time T=10 - drifts upward
BM <- sde.sim(T=10, X0=10, drift=d, sigma=s, sigma.x=sx)
plot(BM, main="Brownian Motion with Drift Sample Path", col="forestgreen", lwd=2)


#############################
# Ornstein-Uhlenbeck process
#############################

# Same set of parameters as above

# 1 sample path (trajectory) of the OU process
set.seed(123)
d <- expression(-5*x)  #drift term: theta = 5
s <- expression(3.5)   #diffusion term: sigma = 3.5
sx <- expression(0)
X <- sde.sim(T=10, X0=10, drift=d, sigma=s, sigma.x=sx)
plot(X, main="Ornstein-Uhlenbeck Sample Path", col="tomato", lwd=2)

# Multiple trajectories of the OU process
set.seed(123)
X2 <- sde.sim(X0=10, drift=d, sigma=s, M=3)
plot(X2, main="Multiple trajectories of OUP", col="darkturquoise", lwd=2)


# Set final time T=10 - reverts to mean 0
X <- sde.sim(T=10, X0=10, drift=d, sigma=s, sigma.x=sx)
plot(X, main="Ornstein-Uhlenbeck Sample Path", col="tomato", lwd=2)



# Exact simulation method - periodic drift
set.seed(123)
d <- expression(sin(x))
d.x <- expression(cos(x)) 
A <- function(x) 1-cos(x)
X3 <- sde.sim(method="EA", delta=1/20, X0=0, N=500, 
        drift=d, drift.x = d.x, A=A)
plot(X3, main="Periodic Drift", col="darkorchid2")



#########################################################
# Also, install the 'Sim.DiffProc' package, call library
# Name comes from: simulate diffusion process
#########################################################

library(Sim.DiffProc)

## 1-dim SDE - Ito 
## dX(t) = 2 * (3 - X(t)) * dt + 2 * X(t) * dW(t)

# Set drift and diffusion terms
f <- expression(2*(3-x))
g <- expression(2*x)

# The function `snssde1d' simulates the numerical solution to 1-dim SDE
run1 <- snssde1d(drift=f,diffusion=g,M=10,x0=1,N=1000)
run1

# Plot sample paths
plot(run1, col="magenta3")



## 1-dim SDE - Stratonovich 
run2 <- snssde1d(drift=f,diffusion=g,M=10,x0=1,N=1000,type="str")
run2

# Plot sample paths
plot(run2, col="limegreen")



## 2-dim SDE - Ito
## dX(t) = 4*(-1-X(t))*Y(t) dt + 0.2 dW1(t)
## dY(t) = 4*(1-Y(t))*X(t) dt + 0.2 dW2(t)
set.seed(1234)

# Set drift and diffusion terms
fx <- expression(4*(-1-x)*y , 4*(1-y)*x )
gx <- expression(0.25*y,0.2*x)

# The function `snssde2d' simulates the numerical solution to 2-dim SDE
mod2d1 <- snssde2d(drift=fx,diffusion=gx,x0=c(x0=1,y0=-1),M=1000)
mod2d1
summary(mod2d1)

# Plot
#dev.new()
plot(mod2d1,type="n")
mx <- apply(mod2d1$X,1,mean)
my <- apply(mod2d1$Y,1,mean)
lines(time(mod2d1),mx,col=1)
lines(time(mod2d1),my,col=2)
legend("topright",c(expression(E(X[t])),expression(E(Y[t]))),lty=1,inset = 0.01,col=c(1,2),cex=0.95)

# Plot 
#dev.new()
plot2d(mod2d1) ## in plane (O,X,Y)
lines(my~mx,col=2)



##########################
##########################
### Exercises

# Work through the practice problems found at:
# http://www.iam.fmph.uniba.sk/institute/stehlikova/fd16/ex/ex02.html