
####################################################
#   Function definitions for Chapter 7


################################################
# Simulate an M/M/1 queue and return the vector
# of the waiting times
MM1_Queue_V <- function(lambda=1, mu=1.1, howmany=10000){
  W = numeric(howmany)
  for (j in 2:howmany) {
    intTime = rexp(1,rate=lambda)
    servTime = rexp(1,rate=mu)
    W[j] = max(0, W[j-1] + servTime - intTime)
  }
  return(W)
}

################################################
# Pricing of a chooser option (also uses vanilla
# calls and puts for exact valuation)

EuVanillaCall <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  d2 <- d1-sqrt(T)*sigma
  value <- S0*pnorm(d1)-K*exp(-r*T)*pnorm(d2)
  return(value)  
}

EuVanillaPut <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  d2 <- d1-sqrt(T)*sigma
  value <- K*exp(-r*T)*pnorm(-d2) - S0*pnorm(-d1)
  return(value)  
}

# Exact chooser price
AYLI <- function(S0,K,r,T1,T2,sigma){
  Kc <- K*exp(-r*(T2-T1))
  value <- EuVanillaCall(S0,K,T2,r,sigma) + 
    EuVanillaPut(S0,Kc,T1,r,sigma)
  return(value)
}

# Naive MC to price the chooser
AYLIMC <- function(S0,K,r,T1,T2,sigma,NRepl1,NRepl2){
  # compute auxiliary quantities outside the loop
  DeltaT <- T2-T1
  muT1 <- (r-sigma^2/2)*T1
  muT2 <- (r-sigma^2/2)*(T2-T1)
  siT1 <- sigma*sqrt(T1)
  siT2 <- sigma*sqrt(T2-T1)
  # vector to contain payoffs
  DiscountedPayoffs <- numeric(NRepl1*NRepl2)
  # sample at time T1
  Samples1 <- rnorm(NRepl1)
  PriceT1 <- S0*exp(muT1 + siT1*Samples1)
  for (k in 1:NRepl1){
    Samples2 <- rnorm(NRepl2)
    PriceT2 <- PriceT1[k]*exp(muT2 + siT2*Samples2)
    ValueCall <- exp(-r*DeltaT)*mean(pmax(PriceT2-K, 0))
    ValuePut <- exp(-r*DeltaT)*mean(pmax(K-PriceT2, 0))
    if (ValueCall > ValuePut)
      DiscountedPayoffs[(1+(k-1)*NRepl2):(k*NRepl2)] <- 
      exp(-r*T2)*pmax(PriceT2-K, 0)
    else
      DiscountedPayoffs[(1+(k-1)*NRepl2):(k*NRepl2)] <-
      exp(-r*T2)*pmax(K-PriceT2, 0)
  }
  aux <- t.test(DiscountedPayoffs)
  price <- as.numeric(aux$estimate)
  conf.int <- as.numeric(aux$conf.int)
  return(list(price=price, conf.int=conf.int))
}

