#  Krishna Kumar 
#  kriskumar@earthlink.net
#  Last modified 10/07/2004
#  http://financemonk.com
#  Have fun!. No rights reserved.
#
#


BinTreeParams<-function(imod, S, X, r, q, tyr, sigma, nstep)
{
  #Returns Vector of Binomial Tree Parameters (imod = 0 for JR, 1 for CRR, 2 for LR)
  #Returns as row vector with u, d then p
  #BSDTwo fn
  #PPNormInv fn
  # variables: delt, nu, ermqdt, rnmudt, u, d, p, d2, pdash
   if (imod == 2) 
   {
   nstep <- nstep+ (nstep %% 2)
   }
   delt <- tyr / nstep
   nu <- (r - q - 0.5 * sigma ^ 2)
   ermqdt <- exp((r - q) * delt)
   if (imod == 0) {
   u <- exp(nu * delt + sigma * sqrt(delt))
   d <- exp(nu * delt - sigma * sqrt(delt))
   p <- 0.5
   }
   else if (imod == 1) {
   u <- exp(sigma * sqrt(delt))
   d <- 1 / u
   p <- 0.5 + (nu * sqrt(delt)) / (2 * sigma)
   }
   else {
   d2 <- Dtwo(S, X, r, q, tyr, sigma)
   p <- PPNormInv(d2, nstep)
   pdash <- PPNormInv(d2 + sigma * sqrt(tyr), nstep)
   u <- ermqdt * pdash / p
   d <- ermqdt * (1 - pdash) / (1 - p)
   }
   BinTreeParams <- c(u, d, p)
 }

optionvalue.euro<-function(imod=2, iopt=1, S=80, X=110, r=0.05, q=0, tyr=1, sigma=0.3, nstep=100)
{
   #Returns Binomial European Option Value (imod<-0 for JR, 1 for CRR, 2 for LR)
   #Uses BinTreeParamVec fn
   # Variables : u, d, p, pdash, a, bind1, bind2, v,udpvec 
#Debug cat("imod",imod,"\n")
   if (imod == 2)   nstep <- nstep + nstep%%2
   if (S > 0 && X > 0 && tyr > 0 && sigma > 0 )
   {
   udpvec <- BinTreeParams(imod, S, X, r, q, tyr, sigma, nstep)
   u <- udpvec[1]
   d <- udpvec[2]
   p <- udpvec[3]
   pdash <- p * u / exp((r - q) * tyr / nstep)
   a <- max(0, 1 + floor(log(X / (S * d ^ nstep)) / log(u / d)))
   bind1 <- 1 - pbinom(a - 1, nstep, pdash)
   bind2 <- 1 - pbinom(a - 1, nstep, p)
 if (a <= nstep) { 
  	 v <- S * exp(-q * tyr) * bind1 - X * exp(-r * tyr) * bind2
   	if (iopt == -1) v <- v - S * exp(-q * tyr) + X * exp(-r * tyr)
}
   else
  {
   v <- 0
   }
	}
   else {
   v <- -1 
	}

optionvalue.euro<- v
} 



greek.binomial<- function(imod=1, TypeFlag="ce", S=100, X=110, r=0.05, q=0, tyr=1, sigma=0.3, nstep=100)
{
# greek.binomial(imod=1,TypeFlag= "pa", S = 50, X = 50, tyr=5/12, r = 0.1, q = 0, sigma = 0.4, nstep = 50)
# the extended CRR tree with node extensions!.
   if (imod == 2)   nstep <- nstep + nstep%%2
    mstep<-nstep+2
    TypeFlag = TypeFlag[1]
    
    z = NA
   
    if (TypeFlag == "ce" || TypeFlag == "ca") 
        z <- +1
    if (TypeFlag == "pe" || TypeFlag == "pa") 
        z <- -1
    if (is.na(z)) stop("TypeFlag misspecified: ce|ca|pe|pa")
   dtval = tyr/nstep
   udpvec <- BinTreeParams(imod, S, X, r, q, tyr, sigma, nstep)
   u <- udpvec[1]
   d <- udpvec[2]
   p <- udpvec[3]
   Df = exp(-r * dtval)
   OptionValue = z * (S * u^(0:mstep) * d^(mstep:0) - X)
   OptionValue = (abs(OptionValue) + OptionValue)/2
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for (j in seq(from = mstep - 1, to = 2, by = -1)) 
		for (i in 0:j) 
		{
	{ OptionValue[i + 1] = (p * OptionValue[i + 2] + (1 - p) * OptionValue[i + 
            1]) * Df }
			}
    }
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for (j in seq(from = mstep - 1, to = 2, by = -1)) for (i in 0:j) OptionValue[i + 
            1] = max((z * (S * u^i * d^(abs(i - j)) - X)), (p * 
            OptionValue[i + 2] + (1 - p) * OptionValue[i + 1]) * 
            Df)
    }
#cat(OptionValue[1]," ", OptionValue[3]," ",OptionValue[2],"\n")
    greek.binomial<-(OptionValue[1]-OptionValue[3])/(S*u^2-S*d^2)
  # cat("greek binomial :", greek.binomial,"\n")

}



recursive.binomial<- function(imod=2, TypeFlag="ce", S=100, X=110, r=0.05, q=0, tyr=1, sigma=0.3, nstep=100)
{
# recursive.binomial(imod=2,TypeFlag= "pa", S = 50, X = 50, tyr=5/12, r = 0.1, q = 0, sigma = 0.4, nstep = 5)
#

    TypeFlag = TypeFlag[1]
    z = NA
    if (TypeFlag == "ce" || TypeFlag == "ca") 
        z = +1
    if (TypeFlag == "pe" || TypeFlag == "pa") 
        z = -1
    if (is.na(z)) stop("TypeFlag misspecified: ce|ca|pe|pa")
   dtval = tyr/nstep
   udpvec <- BinTreeParams(imod, S, X, r, q, tyr, sigma, nstep)
   u <- udpvec[1]
   d <- udpvec[2]
   p <- udpvec[3]
   Df = exp(-r * dtval)
   OptionValue = z * (S * u^(0:nstep) * d^(nstep:0) - X)
   OptionValue = (abs(OptionValue) + OptionValue)/2
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for (j in seq(from = nstep - 1, to = 0, by = -1)) for (i in 0:j) OptionValue[i + 
            1] = (p * OptionValue[i + 2] + (1 - p) * OptionValue[i + 
            1]) * Df
    }
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for (j in seq(from = nstep - 1, to = 0, by = -1)) for (i in 0:j) OptionValue[i + 
            1] = max((z * (S * u^i * d^(abs(i - j)) - X)), (p * 
            OptionValue[i + 2] + (1 - p) * OptionValue[i + 1]) * 
            Df)
    }
cat(OptionValue[1])
    recursive.binomial<-OptionValue[1]
}


optionvalue.amer<-function(imod=2, iopt=1, S=100, X=110, r=0.05, q=0, tyr=1, sigma=0.3, nstep=100)
{
   #Returns Binomial European Option Value (imod<-0 for JR, 1 for CRR, 2 for LR)
   #Uses BinTreeParamVec fn
   # Variables : u, d, p, pdash, a, bind1, bind2, v,udpvec 
   
   #Debug cat("imod",imod,"\n")
   dtval<-tyr/nstep
   if (imod == 2)   nstep <- nstep + nstep%%2
   if (S > 0 && X > 0 && tyr > 0 && sigma > 0 )
   {
   udpvec <- BinTreeParams(imod, S, X, r, q, tyr, sigma, nstep)
   u <- udpvec[1]
   d <- udpvec[2]
   p <- udpvec[3]
   a <- max(0, 1 + floor(log(X / (S * d ^ nstep)) / log(u / d)))
 if (a <= nstep) {
dpowers <- d^(t(c(nstep:0)))
upowers <- u^(t(c(0:nstep)))
scale1 <- p*exp(-r*dtval)
scale2 <- (1-p)*exp(-r*dtval)

# Option values at time T
v <- pmax(iopt*(S*dpowers*upowers-X),0)
# Re-trace to get option value at time zero
for (i in nstep:1)
{
    Si <- S*dpowers[(nstep-i+2):(nstep+1)]*upowers[1:i]
    v <- pmax(pmax(iopt*(Si-X),0),scale1*v[2:(i+1)] + scale2*v[1:i])
}
#cat("\n Time Step=",nstep," Option Value =",v,"\n")
}
   else
  {
   v <- 0
   }
	}
   else {
   v <- -1 
	}

optionvalue.amer<- v
}


optionvalue.bbs<-function(imod=2, TypeFlag="c", S=100, X=110, r=0.05, q=0, tyr=1, sigma=0.3, nstep=100)
{
   # optionvalue.bbs(imod=2,TypeFlag= "c", S = 50, X = 50, tyr=5/12, r = 0.1, q = 0, sigma = 0.4, nstep = 5)

    TypeFlag = TypeFlag[1]
    z = NA
    if (TypeFlag == "c") 
        z = +1
    if (TypeFlag == "p") 
        z = -1
    if (is.na(z)) stop("TypeFlag misspecified: ce|ca|pe|pa")
   dtval = tyr/nstep
   udpvec <- BinTreeParams(imod, S, X, r, q, tyr, sigma, nstep)
   u <- udpvec[1]
   d <- udpvec[2]
   p <- udpvec[3]
   Df = exp(-r * dtval)
   OptionValue = z * (S * u^(0:nstep) * d^(nstep:0) - X)
   OptionValue = (abs(OptionValue) + OptionValue)/2

  for (i in 0:nstep-1)
{
OptionValue[i]<-GBSOption(TypeFlag,S=S,X=X,r=r,Time=dtval,b=r-q,sigma=sigma)$price

}
        for (j in seq(from = nstep - 1, to = 0, by = -1)) 
   {  
 for (i in 0:j) 

 { 
 OptionValue[i + 1] = max((z * (S * u^i * d^(abs(i - j)) - X)), (p * OptionValue[i + 2] + (1 - p) * OptionValue[i + 1]) * Df) 

  }
 }   

optionvalue.bbs<- OptionValue[1]
}
#######################
















PPNormInv<-function(z, n){
   #Returns the Peizer-Pratt Inversion
   #Only defined for n odd
   #Used in LR Binomial Option Valuation
   n <- n + n%%2
   c1 <- exp(-((z / (n + 1 / 3 + 0.1 / (n + 1))) ^ 2) * (n + 1 / 6))
   PPNormInv <- 0.5 + sign(z) * sqrt((0.25 * (1 - c1)))
}

Dtwo<-function(S, X, r, q, tyr, sigma)
{
#Returns the Black-Scholes d2 value
   BSDTwo <- (log(S / X) + (r - q - 0.5 * sigma ^ 2) * tyr) / (sigma * sqrt(tyr))
}
