####################################################
#    Scripts for Chapter 6 of
#    P. Brandimarte
#    Handbook in Monte Carlo Simulation
#    Applications in Financial Engineering,
#    Risk Management, and Economics
#    Wiley, 2014

#  Source the function file
source("Chapter06_Functs.R")

####################################################
#   Chunk 1: Sampling SCALAR geometric Brownian motion
####################################################

### plotting sample paths for different volatilities
set.seed(55555)
paths1 <- simGBM(50,0.1,0.1,1,300,5)
paths2 <- simGBM(50,0.1,0.4,1,300,5)
yBounds <- c(min(min(paths1),min(paths2)),
             max(max(paths1),max(paths2)))
par(mfrow=c(1,2))
plot(1:301,paths1[1,],ylim=yBounds,type='l',
     xlab="sigma=10%",ylab='')
for (j in 2:5)
  lines(1:301,paths1[j,])
plot(1:301,paths2[1,],ylim=yBounds,type='l',
     xlab="sigma=40%",ylab='')
for (j in 2:5)
  lines(1:301,paths2[j,])
par(mfrow=c(1,1))

### comparing vectorized and non-vectorized code
system.time(simGBM(50,0.1,0.1,1,3000,100))
system.time(simvGBM(50,0.1,0.1,1,3000,100))


####################################################
#   Chunk 2: Pricing a vanilla call by crude MC
####################################################

S0 <- 50
K <- 60
T <- 1
r <- 0.04
sigma <- 0.4
numRepl <- 50000
EuVanillaCall(S0,K,T,r,sigma)
set.seed(55555)
out <- BlsMCNaive(S0,K,T,r,sigma,numRepl)
out$estimate
out$conf.int
100*(out$conf.int[2]-out$conf.int[1])/2/out$estimate


####################################################
#   Chunk 3: Multidimensional GBM
####################################################

### Plotting bidimensional GBM
set.seed(55555)
S0 <- c(30,50)
mus <- c(0.1,0.1)
sigmas <- c(0.3,0.3)
paths1 <- sim2GBM(S0,mus,sigmas,0.3,1,300)
paths2 <- sim2GBM(S0,mus,sigmas,0.9,1,300)
yBounds <- c(min(min(paths1),min(paths2)),
             +              max(max(paths1),max(paths2)))
par(mfrow=c(1,2))
plot(1:301,paths1[1,],ylim=yBounds,type='l',
     xlab="rho=0.4",ylab='')
lines(1:301,paths1[2,])
plot(1:301,paths2[1,],ylim=yBounds,type='l',
     xlab="rho=0.9",ylab='')
lines(1:301,paths2[2,])
par(mfrow=c(1,1))

####################################################
#   Chunk 4: Path generation by the Brownian brdge
####################################################

### Checking path generation for the standard Wiener process 
### by a Brownian bridge.}\label{fig:bbridge1}
set.seed(55555)
numRepl <- 100000
T <- 1
numSteps <- 4
Wpaths <- matrix(0, nrow=numRepl, ncol=1+numSteps)
for (i in 1:numRepl)
  Wpaths[i,] <- WienerBridge(T, numSteps)
m <- colMeans(Wpaths[,-1]);m
sdev <- apply(Wpaths[,-1],2,sd);sdev
sqrt((1:numSteps)*T/numSteps)


####################################################
#   Chunk 5: Path generation for the Vasicek 
#            short rate model
####################################################

### sample and plot for the Vasicek model
r0 <- 0.02; rbar <- 0.05; gamma <- 0.1 
sigma <- 0.015; T <- 1
numSteps <- 250; numRepl <- 40
set.seed(55555)
par(mfrow=c(2,1))
paths <- simVasicek(r0,rbar,gamma,sigma,T,numSteps,numRepl)
plot(paths[1,], ylim=range(paths,rbar), type='l', 
     col='mediumorchid2',ylab='',xlab=paste('gamma=',gamma))
for(i in 2:numRepl)
  lines(paths[i,], col=colors()[floor(runif(1,1,657))] )
abline( h=rbar, col='black',lwd=2)
text(0, rbar+0.005, paste("long-term rate =",rbar),adj=0)
gamma <- 5 
set.seed(55555)
paths <- simVasicek(r0,rbar,gamma,sigma,T,numSteps,numRepl)
plot(paths[1,], ylim=range(paths,rbar), type='l', 
     col='mediumorchid2',ylab='',xlab=paste('gamma=',gamma))
for(i in 2:numRepl)
  lines(paths[i,], col=colors()[floor(runif(1,1,657))] )
abline( h=rbar, col='black',lwd=2)
text(0, rbar+0.005, paste("long-term rate =",rbar),adj=0)
par(mfrow=c(1,1))


####################################################
#   Chunk 6: Path generation for the CIR
#            short rate model
####################################################

### get exact density and plot against estimate
### based on Euler scheme (using kernel estimates)
rbar <- 0.05
r0 <- 0.04
gamma <- 0.2
alpha <- 0.1^2
par(mfrow=c(2,1))
# case T=0.25
set.seed(55555)
exact1 <- simCIR(r0,rbar,gamma,alpha,0.25,1,100000)
euler1 <- simEulerCIR(r0,rbar,gamma,alpha,0.25,1,100000)
plot(density(exact1[,2]),ylab='',xlab='rates',
     main='T=0.25',xlim=range(exact1[,2],euler1[,2]))
lines(density(euler1[,2]),lty='dashed')
# case T=1
set.seed(55555)
exact2 <- simCIR(r0,rbar,gamma,alpha,1,1,100000)
euler2 <- simEulerCIR(r0,rbar,gamma,alpha,1,1,100000)
plot(density(exact2[,2]),ylab='',xlab='rates',
     main='T=1',xlim=range(exact2[,2],euler2[,2]))
lines(density(euler2[,2]),lty='dashed')
par(mfrow=c(1,1))


####################################################
#   Chunk 7: Simulating a jump--diffusion process
####################################################

### generalized Wiener process plus jumps
set.seed(55555)
T <- 1; numSteps <- 1000; h <- T/numSteps
# parameters of the jump component
lambda <- 7; nu <- 0; xi <- 2
# parameters of the diffusion component
beta <- 10;  sigma <- 2
# sample the diffusion increments 
dX <- beta*h + sigma*sqrt(h)*rnorm(numSteps)
# sample and add the jump component
N <- rpois(1, lambda*T)
jumpPos <- ceiling(runif(N)*T*numSteps) 
jumpSize <- rnorm(N,nu,xi)
dX[jumpPos] <- dX[jumpPos] + jumpSize
# add components and cumulate sample path
X <- c(0, cumsum(dX))
if (N==0){
  plot(0:numSteps,X,type='l')
}else{
  plot(0:numSteps,X,type='l')
  points(jumpPos,X[jumpPos])
  points(jumpPos,X[jumpPos+1],pch=19)
}
