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

#  Required packages

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

####################################################
#   Chunk 1: Estimation of average waiting time
#            in an M/M/1 queue
####################################################

### check the variability of the sample mean
set.seed(55555)
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())
mean(MM1_Queue_V())

### find naive (WRONG!) confidence intervals 
set.seed(55555)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)
as.numeric(t.test(MM1_Queue_V())$conf.int)

### see the effect of utilization on autocorrelation
set.seed(55555)
W1 <- MM1_Queue_V()
W2 <- MM1_Queue_V(lambda=1,mu=2)
par(mfrow=c(1,2))
acf(W1, main="utilization = 90.91%")
acf(W2, main="utilization = 50%")
par(mfrow=c(1,1))

####################################################
#   Chunk 2: Pricing a chooser option
####################################################

### Find the exact price and compare with vanillas
S0 <- 50; K <- 50; r <- 0.05
T1 <- 2/12; T2 <- 7/12; sigma <- 0.4
EuVanillaCall(S0,K,T2,r,sigma)
EuVanillaPut(S0,K,T2,r,sigma)
AYLI(S0,K,r,T1,T2,sigma)

### Find the price by naive MC (WRONG confidence intervals!!!)
NRepl1 <- 1000
NRepl2 <- 1000
set.seed(55555)
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int
AYLIMC(S0,K,r,T1,T2,sigma,NRepl1,NRepl2)$conf.int


