####################################################
# Function definitions for Chapter 1


####################################################
# estimation of pi
estPi <- function(N){
  x <- 2*runif(N)-1
  y <- 2*runif(N)-1
  h <- sum(x^2+y^2<=1)
  return(4*h/N)
}


####################################################
# simulation of an AR process
AR <- function(X0,a,mu,sigma,T){
  X <- rep(0,T+1)
  X[1] <- X0
  for (t in 1:T){
    X[t+1] <- mu+a*X[t]+rnorm(1,0,sigma)
  }
  plot(0:T,X,pch=15)
  lines(0:T,X)
}

####################################################
# simulation of geometric Brownian motion
simGBM <- function (S0,mu,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  muT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  pathMatrix <- matrix(nrow=numRepl,ncol=numSteps+1)
  pathMatrix[,1] <- S0
  for (i in 1:numRepl){
    for (j in 2:(numSteps+1)){
      pathMatrix[i,j] <- 
        pathMatrix[i,j-1]*exp(rnorm(1,muT,sigmaT))
    }
  }
  return(pathMatrix)
}

####################################################
# simulation of an M/M/1 queue
MM1_Queue <- 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(mean(W))
}

