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

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

####################################################
#   Chunk 1: Estimation of pi
set.seed(55555)
estPi(1000)
estPi(1000)
estPi(1000)
estPi(1000000)
pi


####################################################
#   Chunk 2: simulation of an AR process
set.seed(55555)
AR(40,.8,8,1,50)


####################################################
#   Chunk 3: simulation of geometric Brownian motion
set.seed(55555)
S0 <- 40
mu <- 0.25
sigma1 <- 0.1
sigma2 <- 0.4
T <- 1
numSteps <- 50
numRepl <- 10
Paths1 <- simGBM(S0,mu,sigma1,T,numSteps,numRepl)
Paths2 <- simGBM(S0,mu,sigma2,T,numSteps,numRepl)
maxVal <- max(Paths1,Paths2)
minVal <- min(Paths1,Paths2)
labx <- bquote(italic(t))    # how to generate plot labels 
laby <- bquote(italic(S(t))) # with italic fonts 
par(mfrow = c(1,2))
plot(Paths1[1,],type="l",ylab=laby,xlab=labx,
     main="Vol. 10%",ylim=c(minVal,maxVal))
for (k in 2:10) lines(Paths1[k,])
plot(Paths2[1,],type="l",ylab=laby,xlab=labx,
     main="Vol. 40%",ylim=c(minVal,maxVal))
for (k in 2:10) lines(Paths2[k,])
par(mfrow = c(1,1))


####################################################
#   Chunk 4: simulation of an M/M/1 queue
set.seed(55555)    # ensure repeatability
lambda = 1; mu = 1.1
rho = lambda/mu    # utilization is 90.90%
rho/mu/(1-rho)     # exact value
MM1_Queue(lambda,mu,100000)
MM1_Queue(lambda,mu,100000)
MM1_Queue(lambda,mu,100000)
# change parameters
lambda = 1
mu = 2
rho = lambda/mu    # utilization is 50% 
rho/mu/(1-rho)     #exact value
MM1_Queue(lambda,mu,100000)
MM1_Queue(lambda,mu,100000)
MM1_Queue(lambda,mu,100000)


####################################################
#   Chunk 5: define and plot the Rastrigin function
rastrigin <- function(x,y) 20+x^2-10*cos(2*pi*x) +
                            y^2-10*cos(2*pi*y) 
x <- seq(-2,2,length=100)
y <- x
z <- outer(x,y,rastrigin)
persp(x,y,z,theta=30, phi=35, expand=0.6,
      col='lightgreen', shade=0.1, ltheta=120,
      ticktype='detailed')
contour(x,y,z)

