
####################################################
#    Function definitions for Chapter 11
####################################################


####################################################
# AUXILIARY
####################################################

### scalar Geometric Brownian motion (vectorized version)
simvGBM <- function (S0,mu,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  nuT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  normMatrix <- matrix(rnorm(numSteps*numRepl,nuT,sigmaT),nrow=numRepl)
  logIncrements <- cbind(log(S0),normMatrix)
  logPath <- t(apply(logIncrements,1,cumsum))
  return(exp(logPath))
}

### MUltivariate GBM
simMVGBM <- function (S0Vet,muVet,sigmaVet,corrMatrix,T,numSteps,numRepl){
  # precompute invariant quantities
  dt <- T/numSteps
  nuT <- (muVet-sigmaVet^2/2)*dt
  sigmaT <- sqrt(dt)*sigmaVet
  lowerCholFactor <- t(chol(corrMatrix))
  numFactors <- length(muVet)
  # prepare three-dimensional array
  pathArray <- array(dim=c(numRepl,numSteps+1,numFactors))
  for (k in 1:numFactors) pathArray[,1,k] <- S0Vet[k]
  for (i in 1:numRepl){
    for (j in 2:(numSteps+1)){
      # generate independent standard normals and then correlate them
      zVet <- rnorm(numFactors)
      epsVet <- lowerCholFactor %*% zVet
      # new sample vector
      pathArray[i,j,]=pathArray[i,j-1,]*exp(nuT+sigmaT*epsVet)
    }
  }
  return(pathArray)
}

### GBM paths using low-discrepancy sequences
lowDiscrGBM <- function (S0,mu,sigma,T,numSteps,numRepl,
                         type="sobol",scrambling=0){
  dt <- T/numSteps
  nuT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  if (toupper(type) == "HALTON")
    uSeq <- halton(numRepl,dim=numSteps)
  else
    uSeq <- sobol(numRepl,dim=numSteps,scrambling=scrambling)
  normMatrix <- nuT + sigmaT*qnorm(uSeq)
  logIncrements <- cbind(log(S0),normMatrix)
  logPath <- t(apply(logIncrements,1,cumsum))
  return(exp(logPath))
}

### European-style 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)  
}

### European-style 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)  
}

### The Vasicek model
simVasicek<-function(r0,rbar,gamma,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  # precompute time invariant quantities
  expg <- exp(-gamma*dt)
  auxc <- rbar*(1-expg)
  vol <- sigma*sqrt((1-expg^2)/(2*gamma))
  paths <- matrix(0,nrow=numRepl,ncol=numSteps+1)
  paths[,1] <- r0
  epsMat <- matrix(rnorm(numRepl*numSteps),nrow=numRepl)
  for (i in 1:numRepl){
    for (t in 1:numSteps)
      paths[i,t+1] <- paths[i,t]*expg+auxc+vol*epsMat[i,t]
  }
  return(paths)
}

### THE CIR model: Exact approach
simCIR <- function(r0,rbar,gamma,alpha,T,
                   numSteps,numRepl){
  dt <- T/numSteps
  # precompute time invariant quantities
  expg <- alpha*(1-exp(-gamma*dt))/(4*gamma)
  nu <- 4*rbar*gamma/alpha
  lambdaC <- exp(-gamma*dt)/expg
  paths <- matrix(0,nrow=numRepl,ncol=numSteps+1)
  paths[,1] <- r0
  for (i in 1:numRepl){
    for (t in 1:numSteps){
      chi <- rchisq(1,df=nu,ncp=lambdaC*paths[i,t])
      paths[i,t+1] <- expg*chi}}
  return(paths)
}


####################################################
# European-style exchange/spread option
####################################################

### Exact pricing formula
EuExchange <- function(V0,U0,sigmaV,sigmaU,rho,T,r){
  sigmahat <- sqrt(sigmaU^2+sigmaV^2-2*rho*sigmaU*sigmaV)
  d1 <- (log(V0/U0) + 0.5*T*sigmahat^2)/(sigmahat*sqrt(T))
  d2 <- d1 - sigmahat*sqrt(T)
  return(V0*pnorm(d1) - U0*pnorm(d2))
}

### Naive Monte Carlo 
EuSpreadNaiveMC <- function(V0,U0,K,sigmaV,sigmaU,rho,T,r,
                            numRepl,c.level=0.95){
  # generate sample paths
  S0Vet <- c(V0,U0); muVet <- rep(r,2) 
  sigmaVet <- c(sigmaV,sigmaU)
  corrMatrix <- matrix(c(1,rho,rho,1),nrow=2) 
  numSteps <- 1
  pathArray <- simMVGBM(S0Vet,muVet,sigmaVet,corrMatrix,
                        T,numSteps,numRepl)
  # collect discounted  payoffs
  dpayoff <- exp(-r*T)* 
    pmax(0,pathArray[,2,1]-pathArray[,2,2]-K)
  # in pathArray, the first index refers to sample paths,
  # the second one to time, and the last one to the asset
  aux <- t.test(dpayoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

### Stratified Monte Carlo (single stratification)
### NOTE: There is an error in the book
EuSpreadStrat1 <- function(V0,U0,K,sigmaV,sigmaU,rho,T,r,
                           numRepl,numStrats,c.level=0.95){
  nuTV <- (r - 0.5*sigmaV^2)*T
  sigmaTV <- sigmaV*sqrt(T)
  nuTU <- (r - 0.5*sigmaU^2)*T
  sigmaTU <- sigmaU*sqrt(T)
  rho2 <- sqrt(1-rho^2)
  numPerStrat <- ceiling(numRepl/numStrats)
  meanStrat <- numeric(numStrats)
  varStrat <- numeric(numStrats)
  for (k in 1:numStrats) {
    U <- (k-1)/numStrats+runif(numPerStrat)/numStrats
    epsV <- qnorm(U)
    epsU <- rnorm(numPerStrat,rho*epsV,rho2)
    VT <- V0*exp(nuTV+sigmaTV*epsV)
    UT <- U0*exp(nuTU+sigmaTU*epsU)
    dPayoff <- exp(-r*T)*pmax(0,VT-UT-K)
    meanStrat[k] <- mean(dPayoff)
    varStrat[k] <- var(dPayoff)
  }
  meanTot <- mean(meanStrat)
  varTot <- mean(varStrat)/numRepl
  alpha <- 1-c.level
  half <- qnorm(1-alpha/2)*sqrt(varTot)
  ci <- c(meanTot-half,meanTot+half)
  return(list(estimate=meanTot,conf.int=ci))
}


### Stratified Monte Carlo (double stratification)
### NOTE: There is an error in the book
EuSpreadStrat2 <- function(V0,U0,K,sigmaV,sigmaU,rho,T,r,
                           numRepl,numStrats,c.level=0.95){
  nuTV <- (r - 0.5*sigmaV^2)*T
  sigmaTV <- sigmaV*sqrt(T)
  nuTU <- (r - 0.5*sigmaU^2)*T
  sigmaTU <- sigmaU*sqrt(T)
  rho2 <- sqrt(1-rho^2)
  numPerStrat <- ceiling(numRepl/numStrats)
  meanStrat <- matrix(0,numStrats,numStrats)
  varStrat <- matrix(0,numStrats,numStrats)
  for (k in 1:numStrats) {
    for (j in 1:numStrats) {
      U1 <- (k-1)/numStrats+runif(numPerStrat)/numStrats
      U2 <- (j-1)/numStrats+runif(numPerStrat)/numStrats
      epsV <- qnorm(U1)
      epsU <- rho*epsV + rho2*qnorm(U2)
      VT <- V0*exp(nuTV+sigmaTV*epsV)
      UT <- U0*exp(nuTU+sigmaTU*epsU)
      dPayoff <- exp(-r*T)*pmax(0,VT-UT-K)
      meanStrat[k,j] <- mean(dPayoff)
      varStrat[k,j] <- var(dPayoff)
    }
  }
  meanTot <- mean(meanStrat)
  varTot <- sum(varStrat)/(numRepl*numStrats^2)
  alpha <- 1-c.level
  half <- qnorm(1-alpha/2)*sqrt(varTot)
  ci <- c(meanTot-half,meanTot+half)
  return(list(estimate=meanTot,conf.int=ci))
}


####################################################
# European-style barrier options
####################################################

### Down-and-out put: Exact price
EuDOPut <- function(S0,K,T,r,sigma,Sb){
  a <- (Sb/S0)^(-1 + (2*r / sigma^2))
  b <- (Sb/S0)^(1 + (2*r / sigma^2)) 
  d1 <- (log(S0/K) + (r+sigma^2 / 2)* T) / (sigma*sqrt(T))
  d2 <- (log(S0/K) + (r-sigma^2 / 2)* T) / (sigma*sqrt(T))
  d3 <- (log(S0/Sb) + (r+sigma^2 / 2)* T) / (sigma*sqrt(T))
  d4 <- (log(S0/Sb) + (r-sigma^2 / 2)* T) / (sigma*sqrt(T))
  d5 <- (log(S0/Sb) - (r-sigma^2 / 2)* T) / (sigma*sqrt(T))
  d6 <- (log(S0/Sb) - (r+sigma^2 / 2)* T) / (sigma*sqrt(T))
  d7 <- (log(S0*K/Sb^2) - (r-sigma^2 / 2)* T) / (sigma*sqrt(T))
  d8 <- (log(S0*K/Sb^2) - (r+sigma^2 / 2)* T) / (sigma*sqrt(T))
  P <- K*exp(-r*T)*(pnorm(d4)-pnorm(d2)-a*(pnorm(d7)-pnorm(d5))) -
    S0*(pnorm(d3)-pnorm(d1)-b*(pnorm(d8)-pnorm(d6)))
  return(P)
}

### Down-and-out put: Naive MC
DownOutPutNaive <- function (S0,K,T,r,sigma,Sb,numSteps,
                             numRepl,c.level=0.95){
  payoff <- numeric(numRepl)
  numCrossed <- 0
  for (k in 1:numRepl){
    path <- simvGBM(S0,r,sigma,T,numSteps,1)
    if (any(path <= Sb)){
      payoff[k] <- 0
      numCrossed <- numCrossed+1
    } else
      payoff[k] <- max(0, K-path[numSteps+1])
  }
  aux <- t.test(exp(-r*T)*payoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci,
              numCrossed=numCrossed))
}

### Down-and-in put: Importance sampling
DownInPutIS <- function(S0,K,T,r,sigma,Sb,numSteps,
                        numRepl,b,c.level=0.95){
  # invariant quantities and vector allocation
  dt <- T/numSteps
  nudt <- (r-0.5*sigma^2)*dt
  bdt <- b*dt
  sidt <- sigma*sqrt(dt)
  numCrossed <- 0
  payoff <- numeric(numRepl)
  ISRatio <- numeric(numRepl)
  # Generate asset paths and payoffs  
  for (k in 1:numRepl){
    vetZ <- nudt-bdt+sidt*rnorm(numSteps)
    logPath <- cumsum(c(log(S0), vetZ))
    path <- exp(logPath)
    if (any(path <= Sb)){
      numCrossed <- numCrossed+1
      payoff[k] <- max(0, K-path[numSteps+1])
      ISRatio[k] <-    # evaluate likelihood ratio
        exp(numSteps*bdt^2/2/sigma^2/dt +
              bdt/sigma^2/dt*sum(vetZ) - 
              numSteps*bdt/sigma^2*(r - sigma^2/2))
    } 
  }
  ISvet <- exp(-r*T)*payoff*ISRatio
  aux <- t.test(ISvet,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci,
              numCrossed=numCrossed))
}


####################################################
# European-style Asian options
####################################################

### Exact geometric average price
EuAsianGeometric <- function(S0,K,T,r,sigma,numSamples){
  dT <- T/numSamples
  nu <- r-sigma^2/2
  a <- log(S0)+nu*dT+0.5*nu*(T-dT)
  b <- sigma^2*dT+sigma^2*(T-dT)*
    (2*numSamples-1)/6/numSamples
  x <- (a-log(K)+b)/sqrt(b)
  value <- exp(-r*T)*
    (exp(a+b/2)*pnorm(x)-K*pnorm(x-sqrt(b)))
  return(value)  
}

### Naive MC for arithmetic average
EuAsianNaive <- function(S0,K,T,r,sigma,numSamples,
                         numRepl,c.level=0.95){
  payoff <- numeric(numRepl)
  for (k in 1:numRepl){
    path <- simvGBM(S0,r,sigma,T,numSamples,1)
    payoff[k] <- max(0, mean(path[-1])-K)
  }
  aux <- t.test(exp(-r*T)*payoff,conf.level=c.level)
  value=as.numeric(aux$estimate)
  ci=as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

### Generic function to use control variates
EuAsianCV <- function(S0,K,T,r,sigma,numSamples,numRepl,
                      SampleCV,expCV,numPilot,c.level=0.95){
  # precompute discount factor 
  df <- exp(-r*T)
  # pilot replications to set control parameter
  cVarsP <- numeric(numPilot)
  arithPrices <- numeric(numPilot)
  for (k in 1:numPilot){
    path <- simvGBM(S0,r,sigma,T,numSamples,1)
    cVarsP[k] <- SampleCV(path)
    arithPrices[k] <- df*max(0,mean(path[-1]) - K)
  }
  covar <- cov(cVarsP, arithPrices)
  cstar <- -covar/var(cVarsP)
  # MC run
  CVestimate <- numeric(numRepl)
  for (k in 1:numRepl){
    path <- simvGBM(S0,r,sigma,T,numSamples,1)
    cVar <- SampleCV(path)
    arithPrice <- df*max(0,mean(path[-1]) - K)
    CVestimate[k] <- arithPrice + cstar*(cVar-expCV)
  }
  aux <- t.test(CVestimate,conf.level=c.level)
  value=as.numeric(aux$estimate)
  ci=as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

### Auxiliary function to compare correlations in CV
EuAsianCompareCorr <- function(S0,K,T,r,sigma,numSamples,
                               numRepl){
  df <- exp(-r*T)
  averagePrices <- numeric(numRepl)
  vanillaPrices <- numeric(numRepl)
  geomPrices <- numeric(numRepl)
  arithPrices <- numeric(numRepl)
  for (k in 1:numRepl){
    path <- simvGBM(S0,r,sigma,T,numSamples,1)
    averagePrices[k] <- mean(path[-1])
    vanillaPrices[k] <- df*max(0,path[numSamples+1] - K)
    geomPrices[k] <- df*max(0,
                            (prod(path[-1]))^(1/numSamples) - K)
    arithPrices[k] <- df*max(0,mean(path[-1]) - K)
  }
  cat('Correlations between arithmetic Asian and\n')
  aux <- cor(averagePrices, arithPrices)
  cat('   average stock price ', aux, '\n')
  aux <- cor(vanillaPrices, arithPrices)
  cat('   vanilla call ', aux, '\n')
  aux <- cor(geomPrices, arithPrices)
  cat('   geometric Asian ', aux, '\n')
}

### Asian option with low-discepancy sequences
EuAsianLOWD <- function(S0,K,T,r,sigma,numSamples,
                        numRepl,type="sobol",scrambling=0){
  paths <- lowDiscrGBM(S0,r,sigma,T,numSamples,
                       numRepl,type,scrambling)
  ave <- rowMeans(paths[,-1])
  payoff <- pmax(0, ave-K)
  return(exp(-r*T)*mean(payoff))
}

####################################################
# American-style options: Broadie-Glasserman trees
####################################################

sampleBGMTree <- function(S0,drift,sigma,T,b,numSteps){
  # preallocate memory
  numNodes <- sum(b^(1:(numSteps-1)))
  tree <- numeric(numNodes+1) 
  # precompute invariant quantities once
  dt <- T/numSteps
  muT = (drift - 0.5*sigma^2)*dt
  sigmaT = sigma*sqrt(dt)
  # fill tree
  tree[1] <- S0
  for(i in 0:(numNodes - b^(numSteps-1))){
    tree[(i*b+1):((i+1)*b)+1] <- tree[i+1]*e
    xp(rnorm(b,muT,sigmaT))
  }
  return(tree)
}

getboundsBG <- function(tree,b,numSteps,T,r,payoffFun){
  # This function computes the upper and lower bound on 
  # a Bermudan-style option based on a single BG tree
  # It should be called repeatedly in order to build a 
  # confidence interval; the two bounds are collected 
  # into a vector c(lbound,ubound)
  
  # compute payoff on the last time layer 
  # (nodes from beginLast to numNodes)
  numNodes <- length(tree)
  beginLast <- numNodes - b^(numSteps-1)
  payoffValues <- payoffFun(tree[beginLast:numNodes])
  # evaluate discount factor
  disc = exp(-r*T/numSteps)
  # initialize tree of upper bounds
  nodesToGo <- beginLast-1 
  ubounds <- numeric(numNodes)
  ubounds[beginLast:numNodes] = payoffValues
  # trace back the tree and compute upper bounds
  for(i in nodesToGo:0){
    continValue <- disc*mean(ubounds[i*b+1:b+1])
    intrinsicValue <- f(tree[i+1])
    ubounds[i+1] <- max(continValue,intrinsicValue)
  }
  # now do the same with lower bounds: 
  # initialize and trace back
  lbounds <- numeric(numNodes)
  lbounds[beginLast:numNodes] = payoffValues
  aux = numeric(b)
  for(i in nodesToGo:0){
    intrinsicValue <- payoffFun(tree[i+1])
    for(j in 1:b){
      aux[j] <- disc*mean(lbounds[i*b+(1:b)[-j]+1])
      aux[j] <- ifelse(intrinsicValue>aux[j], 
                       intrinsicValue, lbounds[i*b+j+1])
    }
    lbounds[i+1] <- mean(aux)
  }
  return(c(lbounds[1], ubounds[1]))
}

####################################################
# American-style options: Linear regression and ADP
####################################################

genericLS <- function(paths,T,r,payoffFun,basisFuns){
  # get number of paths, exercise opportunities, etc.
  size <-dim(paths)
  numPaths <-size[1]
  numSteps <-size[2]-1 # paths include initial state
  S0 <- paths[1,1] # first column holds initial price
  paths <- paths[,-1] # get rid of first column
  dt <-T/numSteps # time step
  df <-exp(-r*dt*(1:numSteps)) # discount factors
  numBasis <-length(basisFuns) # number of basis functions
  # 
  cashFlows <- sapply(paths[,numSteps],payoffFun)
  exerciseTime <- numSteps*rep(1,numPaths)
  for (step in seq(from=numSteps-1, to=1, by=-1)){
    intrinsicValue <- sapply(paths[,step],payoffFun)
    inMoney <- which(intrinsicValue > 0)
    xData <- paths[inMoney,step]
    regrMat <- matrix(nrow=length(xData), ncol=numBasis)
    for (k in 1:numBasis) regrMat[,k] <- 
      sapply(xData,basisFuns[[k]])
    yData <- cashFlows[inMoney]*df[exerciseTime[inMoney]-step]
    alpha <- as.numeric(lm.fit(regrMat,yData)$coefficients)
    continuationValue <- regrMat%*%alpha
    index <- which(intrinsicValue[inMoney]>continuationValue)
    exercisePaths <- inMoney[index]
    cashFlows[exercisePaths] <- intrinsicValue[inMoney[index]]
    exerciseTime[exercisePaths] <- step
  }
  price <- max(payoffFun(S0),mean(cashFlows*df[exerciseTime]))
  return(price)
}

####################################################
# Call option with stochastic volatility (Heston model)
####################################################

CallStochVol <- function(S0,K,T,r,V0,alpha,Vbar,xi,
                         rho,numRepl,numSteps,c.level=0.95){
  # invariant quantities
  dt <- T/numSteps
  alphadt <- alpha*dt
  xidt <- xi*sqrt(dt)
  oneplusrdt <- 1+r*dt
  rho2 <- sqrt(1-rho^2)
  VPath <- numeric(numSteps+1)
  VPath[1] <- V0
  SPath <- numeric(numSteps+1)
  SPath[1] <- S0
  payoff <- numeric(numRepl)
  for (k in 1:numRepl){
    eps1 <- rnorm(numSteps)
    eps2 <- rho*eps1+rho2*rnorm(numSteps)
    for (t in 1:numSteps){
      VPath[t+1] <- max(0,VPath[t]+alphadt*(Vbar-VPath[t])+
                          xidt*sqrt(VPath[t])*eps1[t])
      SPath[t+1] <- max(0,SPath[t]*(oneplusrdt+
                                      sqrt(VPath[t]*dt)*eps2[t]))
    }
    payoff[k] <- max(0, SPath[numSteps+1]-K)
  }
  aux <- t.test(exp(-r*T)*payoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))  
}

####################################################
# The Vasicek model for the short rate
####################################################

### Exact ZC bond price under Vasicek
ZCBVasicek <- function(r0,rbar,gamma,sigma,T,F=100){
  B <- (1-exp(-gamma*T))/gamma
  A <- (B-T)*(rbar-sigma^2/(2*gamma^2))-sigma^2*B^2/(4*gamma)
  return(F*exp(A-B*r0))
}

### ZCB price under Vasicek with MC
ZCBVasicekMC <- function(r0,rbar,gamma,sigma,T,mumRepl,
                         numSteps,F=100,c.level=0.95){
  paths <- simVasicek(r0,rbar,gamma,sigma,T,numSteps,numRepl)
  dt <- T/numSteps
  # approximate integral with sum (paths are on rows!)
  sums <- rowSums(paths[,1:numSteps]*dt)
  prices <- F*exp(-sums)
  aux <- t.test(prices,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))  
}

### Exact ZC bond call option price under Vasicek
ZCBCallVasicek <- function(r0,rbar,gamma,sigma,TB,
                           K,TO,F=100){
  Tmat <- TB-TO
  Kn <- K/F
  B <- (1-exp(-gamma*Tmat))/gamma
  ZTB <- ZCBVasicek(r0,rbar,gamma,sigma,TB,1)
  ZTO <- ZCBVasicek(r0,rbar,gamma,sigma,TO,1)
  SZ <- B*sqrt(sigma^2*(1-exp(-2*gamma*TO))/2/gamma)
  d1 <- log(ZTB/Kn/ZTO)/SZ + SZ/2
  d2 <- d1 - SZ
  unitPrice <- ZTB*pnorm(d1)-Kn*ZTO*pnorm(d2)
  return(F*unitPrice)
}

### Monte Carlo for ZC bond call option price under Vasicek
ZCBCallVasicekMC <- function(r0,rbar,gamma,sigma,TB,K,
                             TO,mumRepl,numSteps,F=100,c.level=0.95){
  # find discount factors and rates at option maturities
  dt <- TO/numSteps
  paths <- simVasicek(r0,rbar,gamma,sigma,TO,
                      numSteps,numRepl)
  sums <- rowSums(paths[,1:numSteps]*dt)
  df <- exp(-sums)
  # find bond prices for each path
  Tmat <- TB-TO
  bondPrices <- ZCBVasicek(paths[,numSteps+1],rbar,
                           gamma,sigma,Tmat,F)
  # find discounted payoff
  dPayoff <- df*pmax(0, bondPrices-K)
  aux <- t.test(dPayoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))  
}

####################################################
# The CIR model for the short rate
####################################################

### Exact ZC bond price under CIR
ZCBCIR <- function(r0,rbar,gamma,alpha,T,F=100){
  psi <- sqrt(gamma^2+2*alpha)
  aux <- (gamma+psi)*(exp(psi*T)-1)+2*psi
  B <- 2*(exp(psi*T)-1)/aux
  A <- 2*rbar*gamma/alpha*log(2*psi*exp((psi+gamma)*T/2)/aux)
  return(F*exp(A-B*r0))
}

### ZCB price under CIR with MC
ZCBCIRMC <- function(r0,rbar,gamma,alpha,T,mumRepl,
                     numSteps,F=100,c.level=0.95){
  paths <- simCIR(r0,rbar,gamma,alpha,T,numSteps,numRepl)
  dt <- T/numSteps
  # approximate integral with sum (paths are on rows!)
  sums <- rowSums(paths[,1:numSteps]*dt)
  prices <- F*exp(-sums)
  aux <- t.test(prices,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))  
}
