
####################################################
#   Function definitions for Chapter 3


################################################
### Function to check the memoryless property of 
#   the exponential dirtsibution empirically
CheckMemory <- function(lambda,t,s,numRepl){
  times <- rexp(numRepl, rate=lambda)
  num_t <- sum(times>t)
  cat('P(T>t) = ', num_t/numRepl,'\n')
  cat('P(T>s) = ', sum(times>s)/numRepl,'\n')
  times_t = times[times>t]
  cat('P(T>t+s|T>t) = ', sum(times_t>(t+s))/num_t,'\n')
}

### Trying different correlation coefficients
tryCorrelation <- function(){
  rho <- 0.8
  numObs <- 10000
  X1 <- rnorm(numObs)
  X2 <- rho*X1 + sqrt(1-rho^2)*rnorm(numObs)
  F1 <- exp(X1)
  cat("Correlations",'\n')
  RXX <- cor(X1,X2,method='pearson')
  RFX <- cor(F1,X2,method='pearson')
  cat("  Pearson:  (X1,X2)=", RXX, " (F1,X2)=" , RFX, '\n')
  RXX <- cor(X1,X2,method='spearman')
  RFX <- cor(F1,X2,method='spearman')
  cat("  Spearman: (X1,X2)=", RXX, " (F1,X2)=" , RFX, '\n')
  RXX <- cor(X1,X2,method='kendall')
  RFX <- cor(F1,X2,method='kendall')
  cat("  Kendall:  (X1,X2)=", RXX, " (F1,X2)=" , RFX, '\n')
}

################################################
### Sample paths of a standard Wiener process
simWiener <- function (T,numSteps,numRepl){
  dt <- T/numSteps
  # generate vector of standard normals 
  epsVet <- rnorm(numSteps*numRepl)
  # transform vector to matrix of increments and cumulate 
  incrMatrix <- sqrt(dt)*cbind(numeric(numRepl),
                               matrix(epsVet,numRepl,numSteps))
  pathMatrix <- t(apply(incrMatrix,1,cumsum))
  return(pathMatrix)
}

### Sample paths of a geometric Brownian motion
simGBM <- function (S0,mu,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  nuT <- (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,nuT,sigmaT))
    }
  }
  return(pathMatrix)
}

################################################
### Black-Scholes-Merton formula for vanilla call
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)  
}

### Black-Scholes-Merton formula for vanilla put
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)  
}

### Black-Scholes-Merton formulas for vanilla call/put greeks
BSMdelta <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(pnorm(d1))  
}

BSMgamma <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(dnorm(d1)/(S0*sigma*sqrt(T)))  
}

BSMvega <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(dnorm(d1)*S0*sqrt(T))
}

BSMthetacall <- 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*dnorm(d1)*sigma/(2*sqrt(T))-
    r*K*exp(-r*T)*pnorm(d2)
  return(value)  
}

BSMthetaput <- 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*dnorm(d1)*sigma/(2*sqrt(T))+
    r*K*exp(-r*T)*pnorm(-d2)
  return(value)  
}

