.packageName <- "fOptions"

# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                  DESCRIPTION:
#  NDF                        Normal distribution function
#  CND                        Cumulative normal distribution function
#  CBND                       Cumulative bivariate normal distribution  
# FUNCTION:                  DESCRIPTION:
#  GBSOption                  Computes Option Price from the GBS Formula
#  GBSCharacteristics         Computes Option Price and all Greeks of GBS Model
#   BlackScholesOption         Synonyme Function Call to GBSOption
#  GBSGreeks                  Computes one of the Greeks of the GBS formula
# FUNCTION:                  DESCRIPTION:
#  Black76Option              Computes Prices of Options on Futures
# FUNCTION:                  DESCRIPTION:
#  MiltersenSchwartzOption    Pricing a Miltersen Schwartz Option
# METHODS:                   DESCRIPTION:
#  print.option               Print Method
#  summary.otion              Summary Method
################################################################################


NDF = 
function(x) 
{   # A function implemented by Diethelm Wuertz
     
    # Description:
    #   Calculate the normal distribution function.
    
    # FUNCTION:
    
    # Compute:
    result = exp(-x*x/2)/sqrt(8*atan(1))
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


CND = 
function(x)
{   # A function implemented by Diethelm Wuertz
     
    # Description:
    #   Calculate the cumulated normal distribution function.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
        
    # FUNCTION:
    
    # Compute:
    k  = 1 / ( 1 + 0.2316419 * abs(x) )      
    a1 =  0.319381530; a2 = -0.356563782; a3 =  1.781477937 
    a4 = -1.821255978; a5 =  1.330274429 
    result = NDF(x) * (a1*k + a2*k^2 + a3*k^3 + a4*k^4 + a5*k^5) - 0.5
    result = 0.5 - result*sign(x)
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


CBND = 
function(x1, x2, rho) 
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Calculate the cumulative bivariate normal distribution function.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Compute:
    # Take care for the limit rho = +/- 1
    a = x1
    b = x2
    if (abs(rho) == 1) rho = rho - (1e-12)*sign(rho)
    # cat("\n a - b - rho :"); print(c(a,b,rho))
    X = c(0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334)
    y = c(0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604)
    a1 = a / sqrt(2 * (1 - rho^2))
    b1 = b / sqrt(2 * (1 - rho^2))   
    if (a <= 0 && b <= 0 && rho <= 0) { 
       Sum1 = 0
       for (I in 1:5) {
            for (j in 1:5) {
            Sum1 = Sum1 + X[I] * X[j] * 
              exp(a1*(2*y[I]-a1) + b1*(2*y[j]-b1) + 
              2*rho*(y[I]-a1)*(y[j]-b1)) } }
       result = sqrt(1 - rho^2) / pi * Sum1 
       return(result) }
    if (a <= 0 && b >= 0 && rho >= 0) {
        result = CND(a) - CBND(a, -b, -rho)
        return(result) }
    if (a >= 0 && b <= 0 && rho >= 0) {
        result = CND(b) - CBND(-a, b, -rho)
        return(result) }
    if (a >= 0 && b >= 0 && rho <= 0) {
        result = CND(a) + CND(b) - 1 + CBND(-a, -b, rho)
        return(result) }
    if (a * b * rho >= 0 ) { 
        rho1 = (rho*a - b) * sign(a) / sqrt(a^2 - 2*rho*a*b + b^2)
        rho2 = (rho*b - a) * sign(b) / sqrt(a^2 - 2*rho*a*b + b^2)
        delta = (1 - sign(a) * sign(b)) / 4
        result = CBND(a, 0, rho1) + CBND(b, 0, rho2) - delta 
        return(result) }
    
    # Return Value:
    invisible()
}


# ******************************************************************************


GBSOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Calculate the Generalized Black-Scholes option
    #   price either for a call or a put option.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
    d2 = d1 - sigma*sqrt(Time)
    if (TypeFlag == "c")
        result = S*exp((b-r)*Time)*CND(d1) - X*exp(-r*Time)*CND(d2) 
    if (TypeFlag == "p")
        result = X*exp(-r*Time)*CND(-d2) - S*exp((b-r)*Time)*CND(-d1) 
    
    # Return Value:
    option = list(
        price = result, 
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


GBSCharacteristics = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz
 
    # Description:
    #   Calculate the Options Characterisitics (Premium 
    #   and Greeks for a Generalized Black-Scholes option 
    #   either for a call or a put option.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Premium and Function Call to all Greeks
    TypeFlag = TypeFlag[1]
    premium = GBSOption(TypeFlag, S, X, Time, r, b, sigma)$price  
    delta = GBSGreeks("Delta", TypeFlag, S, X, Time, r, b, sigma)  
    theta = GBSGreeks("Theta", TypeFlag, S, X, Time, r, b, sigma)
    vega = GBSGreeks("Vega", TypeFlag, S, X, Time, r, b, sigma)
    rho = GBSGreeks("Rho", TypeFlag, S, X, Time, r, b, sigma)
    lambda = GBSGreeks("Lambda", TypeFlag, S, X, Time, r, b, sigma)  
    gamma = GBSGreeks("Gamma", TypeFlag, S, X, Time, r, b, sigma)  
    # Return Value:
    list(premium = premium, delta = delta, theta = theta, 
        vega = vega, rho = rho, lambda = lambda, gamma = gamma) 
} 


# ------------------------------------------------------------------------------


BlackScholesOption = 
function(...) 
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   A synonyme for GBSOption
    
    # FUNCTION:
    
    # Return Value:
    GBSOption(...)
}



# ******************************************************************************


GBSGreeks = 
function(Selection = c("Delta", "Theta", "Vega", "Rho", "Lambda", "Gamma",
"CofC"), TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz
 
    # Description:
    #   Calculate the Options Greeks for a Generalized  
    #   Black-Scholes option either for a call or a put 
    #   option.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    Selection = Selection[1]
    
    # Internal Functions:
    GBSDelta = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        if (TypeFlag == "c") result = exp((b-r)*Time)*CND(d1)
        if (TypeFlag == "p") result = exp((b-r)*Time)*(CND(d1)-1)
        result }
    GBSTheta = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        d2 = d1 - sigma*sqrt(Time)
        Theta1 = -(S*exp((b-r)*Time)*NDF(d1)*sigma)/(2*sqrt(Time))  
        if (TypeFlag == "c") result = Theta1 - 
            (b-r)*S*exp((b-r)*Time)*CND(+d1) - r*X*exp(-r*Time)*CND(+d2) 
        if (TypeFlag == "p") result = Theta1 + 
            (b-r)*S*exp((b-r)*Time)*CND(-d1) + r*X*exp(-r*Time)*CND(-d2) 
        result }
    GBSVega = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        result = S*exp((b-r)*Time)*NDF(d1)*sqrt(Time) # Call,Put
        result }
    GBSRho = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        d2 = d1 - sigma*sqrt(Time)
        CallPut = GBSOption(TypeFlag, S, X, Time, r, b , sigma)$price
        if (TypeFlag == "c") {
            if (b != 0) {result =  Time * X * exp(-r*Time)*CND( d2)} 
            else {result = -Time * CallPut } }
        if (TypeFlag == "p") {
            if (b != 0) {result = -Time * X * exp(-r*Time)*CND(-d2)}
            else { result = -Time * CallPut } }
        result }
    GBSLambda = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        CallPut = GBSOption(TypeFlag,S,X,Time,r,b,sigma)$price
        if (TypeFlag == "c") result = exp((b-r)*Time)* CND(d1)*S / CallPut
        if (TypeFlag == "p") result = exp((b-r)*Time)*(CND(d1)-1)*S / CallPut
        result }        
    GBSGamma = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        result = exp((b-r)*Time)*NDF(d1)/(S*sigma*sqrt(Time)) # Call,Put
        result }
    GBSCofC = function(TypeFlag, S, X, Time, r, b, sigma) {
        d1 = ( log(S/X) + (b+sigma*sigma/2)*Time ) / (sigma*sqrt(Time))
        if (TypeFlag == "c") result = Time*S*exp((b-r)*Time)*CND(d1)
        if (TypeFlag == "p") result = -Time*S*exp((b-r)*Time)*CND(-d1)
        result }
    
    # Function Call to all Greeks via selection parameter 
    result = NA
    if (Selection == "Delta" || Selection == "delta")
            result = GBSDelta (TypeFlag, S, X, Time, r, b, sigma)  
    if (Selection == "Theta" || Selection == "theta")
            result = GBSTheta (TypeFlag, S, X, Time, r, b, sigma)
    if (Selection == "Vega" || Selection == "vega")
            result = GBSVega  (TypeFlag, S, X, Time, r, b, sigma)
    if (Selection == "Rho" || Selection == "rho")
            result = GBSRho   (TypeFlag, S, X, Time, r, b, sigma)
    if (Selection == "Lambda" || Selection == "lambda")
            result = GBSLambda(TypeFlag, S, X, Time, r, b, sigma)  
    if (Selection == "Gamma" || Selection == "gamma")
            result = GBSGamma (TypeFlag, S, X, Time, r, b, sigma)  
    if (Selection == "CofC" || Selection == "cofc")
            result = GBSCofC  (TypeFlag, S, X, Time, r, b, sigma)
    
    # Return Value:
    result    
} 


# ******************************************************************************


Black76Option = 
function(TypeFlag = c("c", "p"), FT, X, Time, r, sigma)
{   # A function implemented by Diethelm Wuertz
 
    # Description:
    #   Calculate Options Price for Black (1977) Options 
    #   on futures/forwards

    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Return Value:
    GBSOption(TypeFlag, S = FT, X = X, Time = Time, r = r, b = 0, sigma = sigma)
}


# ******************************************************************************


MiltersenSchwartzOption = 
function (TypeFlag = c("c", "p"), Pt, FT, X, time, Time, sigmaS, sigmaE, 
sigmaF, rhoSE, rhoSF, rhoEF, KappaE, KappaF) 
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Miltersen Schwartz (1997) commodity option model.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    TyoeFlag = TypeFlag[1]
    
    # Compute:
    vz = sigmaS^2*time+2*sigmaS*(sigmaF*rhoSF*1/KappaF*(time-1/KappaF*
      exp(-KappaF*Time)*(exp(KappaF*time)-1))-sigmaE*rhoSE*1/KappaE*
      (time-1/KappaE*exp(-KappaE*Time)*(exp(KappaE*time)-1)))+sigmaE^2*
      1/KappaE^2*(time+1/(2*KappaE)*exp(-2*KappaE*Time)*(exp(2*KappaE*time)-
      1)-2*1/KappaE*exp(-KappaE*Time)*(exp(KappaE*time)-1))+sigmaF^2*
      1/KappaF^2*(time+1/(2*KappaF)*exp(-2*KappaF*Time)*(exp(2*KappaF*time)-
      1)-2*1/KappaF*exp(-KappaF*Time)*(exp(KappaF*time)-1))-2*sigmaE*
      sigmaF*rhoEF*1/KappaE*1/KappaF*(time-1/KappaE*exp(-KappaE*Time)*
      (exp(KappaE*time)-1)-1/KappaF*exp(-KappaF*Time)*(exp(KappaF*time)-
      1)+1/(KappaE+KappaF)*exp(-(KappaE+KappaF)*Time)*(exp((KappaE+KappaF)*
      time)-1))
    vxz = sigmaF*1/KappaF*(sigmaS*rhoSF*(time-1/KappaF*(1-exp(-KappaF*
      time)))+sigmaF*1/KappaF*(time-1/KappaF*exp(-KappaF*Time)*(exp(KappaF*
      time)-1)-1/KappaF*(1-exp(-KappaF*time))+1/(2*KappaF)*exp(-KappaF*
      Time)*(exp(KappaF*time)-exp(-KappaF*time)))-sigmaE*rhoEF*1/KappaE*
      (time-1/KappaE*exp(-KappaE*Time)*(exp(KappaE*time)-1)-1/KappaF*(1-
      exp(-KappaF*time))+1/(KappaE+KappaF)*exp(-KappaE*Time)*
      (exp(KappaE*time)-exp(-KappaF*time))))
    vz = sqrt(vz)
    d1 = (log(FT/X)-vxz+vz^2/2)/vz
    d2 = (log(FT/X)-vxz-vz^2/2)/vz
    
    # Call/Put:
    if (TypeFlag == "c") {
        result = Pt*(FT*exp(-vxz)*CND(d1)-X*CND(d2)) }
    if (TypeFlag == "p") {
        result = Pt*(X*CND(-d2)-FT*exp(-vxz)*CND(-d1)) }
    
    # Return Value:
    result
}  


# ******************************************************************************


GBSVolatility = function(price, TypeFlag = c("c", "p"), S, X, Time, r, b, 
tol = .Machine$double.eps, maxiter = 10000)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Compute implied volatility
    
    # Example:
    #   sigma = GBSVolatility(price=10.2, "c", S=100, X=90, Time=1/12, r=0, b=0)
    #   sigma
    #   GBSOption("c", S=100, X=90, Time=1/12, r=0, b=0, sigma=sigma)$price
    
    # FUNCTION:
    
    # Option Type:
    TypeFlag = TypeFlag[1]
    
    # Internal Function:
    f = function(x, price, TypeFlag, S, X, Time, r, b, ...) {
        GBS = GBSOption(TypeFlag = TypeFlag, S = S, X = X, Time = Time, 
            r = r, b = b, sigma = x)$price 
        price - GBS}
    
    # Search for Root:
    volatility = uniroot(f, interval = c(-10,10), price = price, 
        TypeFlag = TypeFlag, S = S, X = X, Time = Time, r = r, b = b, 
        tol = tol, maxiter = maxiter)$root
        
    # Return Value:
    volatility
}


# ******************************************************************************


print.option = 
function(x, ...)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Print method for objects of class "option".
    
    # FUNCTION:
    
    # Print Method:
    object = x
    cat("\nCall:", deparse(object$call), "", sep = "\n")
    cat("Option Price:\n")
    cat(object$price, "\n")
}


# ------------------------------------------------------------------------------


summary.option = 
function(object, ...)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Summary method for objects of class "option".
    
    # FUNCTION:
    
    # Summary Method:
    print(object, ...)
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyright (C) 1998-2003 by Diethelm Wuertz for this R-Port


################################################################################
# FUNCTION:                  DESCRIPTION:
#  RollGeskeWhaleyOption      Roll-Geske-Whaley Calls on Dividend Paying Stocks
#  BAWAmericanApproxOption    Barone-Adesi and Whaley Approximation
#  BSAmericanApproxOption     Bjerksund and Stensland Approximation
################################################################################


RollGeskeWhaleyOption = 
function(S, X, time1, Time2, r, D, sigma) 
{   # A function implemented by Diethelm Wuertz
 
    # Description:
    #   Calculates the option price of an American call on a stock
    #   paying a single dividend with specified time to divident
    #   payout. The option valuation formula derived by Roll, Geske 
    #   and Whaley is used.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    big = 100000000
    eps = 1.0e-5
    t1 = time1
    T2 = Time2
    
    # Compute:
    Sx = S - D * exp(-r * t1)
    if (D <= X * (1 - exp(-r*(T2-t1)))) {         
        result = GBSOption("c", Sx, X, T2, r, b=r, sigma)$price
        cat("\nWarning: Not optimal to exercise\n")
        return(result) }
    ci = GBSOption("c", S, X, T2-t1, r, b=r, sigma)$price
    HighS = S
    while ( ci-HighS-D+X > 0 && HighS < big ) {
        HighS = HighS * 2
        ci = GBSOption("c", HighS, X, T2-t1, r, b=r, sigma)$price }
    if (HighS > big) {
        result = GBSOption("c", Sx, X, T2, r, b=r, sigma)$price
        stop()}
    LowS = 0
    I = HighS * 0.5
    ci = GBSOption("c", I, X, T2-t1, r, b=r, sigma)$price 
    # Search algorithm to find the critical stock price I
    while ( abs(ci-I-D+X) > eps && HighS - LowS > eps ) {
         if (ci-I-D+X < 0 ) { HighS = I }
        else { LowS = I }
        I = (HighS + LowS) / 2
        ci = GBSOption("c", I, X, T2-t1, r, b=r, sigma)$price }
    a1 = (log(Sx/X) + (r+sigma^2/2)*T2) / (sigma*sqrt(T2))
    a2 = a1 - sigma*sqrt(T2)
    b1 = (log(Sx/I) + (r+sigma^2/2)*t1) / (sigma*sqrt(t1))
    b2 = b1 - sigma*sqrt(t1)
    result = Sx*CND(b1) + Sx*CBND(a1,-b1,-sqrt(t1/T2)) -
        X*exp(-r*T2)*CBND(a2,-b2,-sqrt(t1/T2)) - 
            (X-D)*exp(-r*t1)*CND(b2)
    
    # Return Value:
    result
}


# ******************************************************************************


BAWAmericanApproxOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz
 
    # Description:
    #   Calculates the option price of an American call or put
    #   option on an underlying asset for a given cost-of-carry rate.
    #   The quadratic approximation method by Barone-Adesi and
    #   Whaley is used.

    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Internal Function - The Call:
    BAWAmCallApproxOption <<- function(S, X, Time, r, b, sigma) {
        # Newton Raphson Algorithm:
        Kc <<- function(X, Time, r, b, sigma) {   
            # Newton Raphson algorithm to solve for the critical commodity 
            # price for a Call.
            # Calculation of seed value, Si
            n = 2*b/sigma^2
            m = 2*r/sigma^2
            q2u = (-(n-1)+sqrt((n-1)^2+4*m))/2
            Su = X/(1-1/q2u)
            h2 = -(b*Time+2*sigma*sqrt(Time))*X/(Su-X)
            Si = X+(Su-X)*(1-exp(h2))
            K = 2*r/(sigma^2*(1-exp(-r*Time)))
            d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
            Q2 = (-(n-1)+sqrt((n-1)^2+4*K))/2
            LHS = Si-X
            RHS = GBSOption("c", Si, X, Time, r, b, sigma)$price + 
                (1-exp((b-r)*Time)*CND(d1))*Si/Q2
            bi = exp((b-r)*Time)*CND(d1)*(1-1/Q2) +
                (1-exp((b-r)*Time)*CND(d1)/(sigma*sqrt(Time)))/Q2
            E = 0.000001
            # Newton Raphson algorithm for finding critical price Si
            while (abs(LHS-RHS)/X > E) {
                Si = (X+RHS-bi*Si)/(1-bi)
                d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
                LHS = Si-X
                RHS = GBSOption("c", Si, X, Time, r, b, sigma)$price + 
                    (1-exp((b-r)*Time)*CND(d1))*Si/Q2
                bi = exp((b-r)*Time)*CND(d1)*(1-1/Q2) + 
                (   1-exp((b-r)*Time)*CND(d1)/(sigma*sqrt(Time)))/Q2 }
            # Return Value:
            Si}
        # Compute:
        if (b >= r) {
            result = GBSOption("c", S, X, Time, r, b, sigma)$price }
        else {
            Sk = Kc(X, Time, r, b, sigma)
            n = 2*b/sigma^2
            K = 2*r/(sigma^2*(1-exp(-r*Time)))
            d1 = (log(Sk/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
            Q2 = (-(n-1)+sqrt((n-1)^2+4*K))/2
            a2 = (Sk/Q2)*(1-exp((b-r)*Time)*CND(d1))
            if (S < Sk) {
                result = GBSOption("c", S, X, Time, r, b, sigma)$price +
                    a2*(S/Sk)^Q2 }
            else {
                result = S-X } }
        # Return Value:
        result }

    # Internal Function - The Put:
    BAWAmPutApproxOption <<- function(S, X, Time, r, b, sigma) {
        # Internal Function:
        Kp <<- function(X, Time, r, b, sigma) {   
            # Newton Raphson algorithm to solve for the critical commodity 
            # price for a Put.
            # Calculation of seed value, Si
            n = 2*b/sigma^2
            m = 2*r/sigma^2
            q1u = (-(n-1)-sqrt((n-1)^2+4*m))/2
            Su = X/(1-1/q1u)
            h1 = (b*Time-2*sigma*sqrt(Time))*X/(X-Su)
            Si = Su+(X-Su)*exp(h1) 
            K = 2*r/(sigma^2*(1-exp(-r*Time)))
            d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
            Q1 = (-(n-1)-sqrt((n-1)^2+4*K))/2
            LHS = X-Si
            RHS = GBSOption("p", Si, X, Time, r, b, sigma)$price -
                (1-exp((b-r)*Time)*CND(-d1))*Si/Q1
            bi = -exp((b-r)*Time)*CND(-d1)*(1-1/Q1) -
                (1+exp((b-r)*Time)*CND(-d1)/(sigma*sqrt(Time)))/Q1
            E = 0.000001
            # Newton Raphson algorithm for finding critical price Si
            while (abs(LHS-RHS)/X > E ) {
                Si = (X-RHS+bi*Si)/(1+bi)
                d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
                LHS = X-Si
                RHS = GBSOption("p", Si, X, Time, r, b, sigma)$price -
                    (1-exp((b-r)*Time)*CND(-d1))*Si/Q1
                bi = -exp((b-r)*Time)*CND(-d1)*(1-1/Q1) -
                    (1+exp((b-r)*Time)*CND(-d1)/(sigma*sqrt(Time)))/Q1 }
            # Return Value:
            Si}
        # Compute:
        Sk = Kp(X, Time, r, b, sigma)
        n = 2*b/sigma^2
        K = 2*r/(sigma^2*(1-exp(-r*Time)))
        d1 = (log(Sk/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
        Q1 = (-(n-1)-sqrt((n-1)^2+4*K))/2
        a1 = -(Sk/Q1)*(1-exp((b-r)*Time)*CND(-d1))
        if (S > Sk) {
            result = GBSOption("p", S, X, Time, r, b, sigma)$price + 
                a1*(S/Sk)^Q1 }
        else {
            result = X-S }  
        # Return Value:
        result}
    
    # Compute:
    if (TypeFlag == "c") {
        result = BAWAmCallApproxOption(S, X, Time, r, b, sigma) }
    if (TypeFlag == "p") {      
        result = BAWAmPutApproxOption(S, X, Time, r, b, sigma) }
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


BSAmericanApproxOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Calculates the option price of an American call or 
    #   put option stocks, futures, and currencies. The 
    #   approximation method by Bjerksund and Stensland is used.
    
    # References:
    #   Haug E.G., The Complete Guide to Option Pricing Formulas
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Utility function phi:
    phi <<- function(S, Time, gamma, H, I, r, b, sigma) {
        lambda = (-r + gamma*b + 0.5*gamma * (gamma-1)*sigma^2) * Time
        d = -(log(S/H) + (b + (gamma-0.5)*sigma^2)*Time) / 
            (sigma*sqrt(Time))
        kappa = 2 * b / (sigma^2) + (2*gamma - 1)
        result = exp(lambda)*S^gamma * 
        (CND(d)-(I/S)^kappa*CND(d-2*log(I/S)/(sigma*sqrt(Time))))
        result }
    # Call Approximation:
    BSAmericanCallApprox <<- function(S, X, Time, r, b, sigma) { 
        if (b >= r) { 
        # Never optimal to exersice before maturity
        result = list(
        Premium=GBSOption("c", S, X, Time, r, b, sigma)$price,
          TriggerPrice=NA)}
      else {
        Beta = (1/2 - b/sigma^2) + sqrt((b/sigma^2 - 1/2)^2 + 2*r/sigma^2)
        BInfinity = Beta/(Beta-1) * X
        B0 = max(X, r/(r-b) * X)
        ht = -(b*Time + 2*sigma*sqrt(Time)) * B0/(BInfinity-B0)
        # Trigger Price I:
        I = B0 + (BInfinity-B0) * (1 - exp(ht))
        alpha = (I-X) * I^(-Beta)
        if (S >= I) { 
          result = list(
           Premium=S-X, TriggerPrice=I) }
        else {
          result = list(
            Premium=alpha*S^Beta - alpha*phi(S,Time,Beta,I,I,r,b,sigma) + 
              phi(S,Time,1,I,I,r,b,sigma) - phi(S,Time,1,X,I,r,b,sigma) - 
              X*phi(S,Time,0,I,I,r,b,sigma) + 
            X*phi(S,Time,0,X,I,r,b,sigma), TriggerPrice=I) } }
      result}
    
    # The Bjerksund and Stensland (1993) American approximation:
    if (TypeFlag == "c") {
      result = BSAmericanCallApprox(S, X, Time, r, b, sigma) }
    if (TypeFlag == "p") {
      # Use the Bjerksund and Stensland put-call transformation
      result = BSAmericanCallApprox(X, S, Time, r - b, -b, sigma) }
    
    # Return Value:
    result
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                 DESCRIPTION:
#  CRRBinomialTreeOption     Cox-Ross-Rubinstein Binomial Tree Option Model
#  JRBinomialTreeOption      JR Modfication to the Binomial Tree Option
#  TIANBinomialTreeOption    Tian's Modification to the Binomial Tree Option
# FUNCTION:
#  BinomialTreeOption        CRR Binomial Tree Option with Cost of Carry Term
#  BinomialTreePlot          Plots results from the CRR Option Pricing Model
################################################################################


CRRBinomialTreeOption = 
function(TypeFlag = c("ce", "pe", "ca", "pa"), S, X, Time, r, b, sigma, n)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Cox-Ross-Rubinstein Binomial Tree Option Model
    
    # FUNCTION:
    
    # Check Flags:
    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")
  
    # Parameters:
    dt = Time/n
    u  = exp(sigma*sqrt(dt))
    d  = 1/u
    p  = (exp(b*dt)-d)/(u-d)
    Df = exp(-r*dt)
    
    # Iteration:
    OptionValue = z*(S*u^(0:n)*d^(n:0) - X)
    OptionValue = (abs(OptionValue)+OptionValue)/2
    
    # European Option:
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for ( j in seq(from=n-1, to=0, by=-1) ) 
            for ( i in 0:j )         
                OptionValue[i+1] = 
                (p*OptionValue[i+2] + (1-p)*OptionValue[i+1]) * Df }
    
    # American Option:
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for ( j in seq(from=n-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) }
    
    # Return Value:
    OptionValue[1]
}


# ------------------------------------------------------------------------------


JRBinomialTreeOption = 
function(TypeFlag = c("ce", "pe", "ca", "pa"), S, X, Time, r, b, sigma, n)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   JR Modfication to the Binomial Tree Option
    
    # FUNCTION:
    
    # Check Flags:
    TypeFlag = TypeFlag[1]
    if (TypeFlag == "ce" || TypeFlag == "ca") z = +1
    if (TypeFlag == "pe" || TypeFlag == "pa") z = -1
    
    # Parameters:
    dt = Time/n
    u = exp( (r-sigma^2/2)*dt+sigma*sqrt(dt) )
    d = exp( (r-sigma^2/2)*dt-sigma*sqrt(dt) )
    p = 1/2
    Df = exp(-r*dt)
    
    # Iteration:
    OptionValue = z*(S*u^(0:n)*d^(n:0) - X)
    OptionValue = (abs(OptionValue)+OptionValue)/2
    
    # European Option:
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for ( j in seq(from=n-1, to=0, by=-1) ) 
            for ( i in 0:j )         
                OptionValue[i+1] = 
                (p*OptionValue[i+2] + (1-p)*OptionValue[i+1]) * Df }
    
                # American Option:
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for ( j in seq(from=n-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) }
    
    # Return Value:
    OptionValue[1]
}


# ------------------------------------------------------------------------------


TIANBinomialTreeOption = 
function(TypeFlag = c("ce", "pe", "ca", "pa"), S, X, Time, r, b, sigma, n)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Tian's Modification to the Binomial Tree Option
    
    # FUNCTION:
    
    # Check Flags:
    TypeFlag = TypeFlag[1]
    if (TypeFlag == "ce" || TypeFlag == "ca") z = +1
    if (TypeFlag == "pe" || TypeFlag == "pa") z = -1  
    
    # Parameters:
    dt = Time/n 
    M = exp ( b*dt )
    V = exp ( sigma^2 * dt )
    u = (M*V/2) * ( V + 1 + sqrt(V*V + 2*V - 3) )
    d = (M*V/2) * ( V + 1 - sqrt(V*V + 2*V - 3) )
    p = (M-d)/(u-d)
    Df = exp(-r*dt)
    
    # Iteration:
    OptionValue = z*(S*u^(0:n)*d^(n:0) - X)
    OptionValue = (abs(OptionValue)+OptionValue)/2
    
    # European Option:
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for ( j in seq(from=n-1, to=0, by=-1) ) 
            for ( i in 0:j )         
                OptionValue[i+1] = 
                (p*OptionValue[i+2] + (1-p)*OptionValue[i+1]) * Df }
    
    # American Option:
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for ( j in seq(from=n-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) }
                    
    # Return Value:
    OptionValue[1]
}


# ******************************************************************************


BinomialTreeOption = 
function(TypeFlag = c("ce", "pe", "ca", "pa"), S, X, Time, r, b, sigma, n)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Calculates option prices from the Cox-Ross-Rubinstein
    #   Binomial tree model.
    
    # Note:
    #   The model described here is a version of the CRR Binomial
    #   Tree model. Including a cost of carry term b, the model can
    #   used to price European and American Options on
    #     b=r       stocks
    #     b=r-q     stocks and stock indexes paying a continuous  
    #               dividend yield q
    #     b=0       futures
    #     b=r-rf    currency options with foreign interst rate rf
    
    # Example:
    #   par(mfrow=c(1,1))
    #   Tree = BinomialTree("pa", 100, 95, 0.5, 0.08, 0.08, 0.3, 5)
    #   print(round(Tree, digits=3))
    #   BinomialTreePlot(Tree, main="American Put Option")
    #
    # Reference:
    #   E.G. Haug, The Complete Guide to Option Pricing Formulas, 
    #   1997, Chapter 3.1.1
    
    # FUNCTION:
    
    # Check Flags:
    TypeFlag = TypeFlag[1]
    if (TypeFlag == "ce" || TypeFlag == "ca") z = +1
    if (TypeFlag == "pe" || TypeFlag == "pa") z = -1    
    
    # Parameters:
    dt = Time / n
    u  = exp(sigma*sqrt(dt))
    d  = 1 / u
    p  = (exp(b*dt) - d) / (u - d)
    Df = exp(-r*dt)
    
    # Algorithm:
    OptionValue = z*(S*u^(0:n)*d^(n:0) - X)
    offset = 1
    Tree = OptionValue = (abs(OptionValue)+OptionValue)/2   
    
    # European Type:
    if (TypeFlag == "ce" || TypeFlag == "pe") {
        for (j in (n-1):0) {
            Tree <-c(Tree, rep(0, times=n-j))
            for (i in 0:j) {         
                OptionValue[i+offset] = 
                    (p*OptionValue[i+1+offset] + 
                (1-p)*OptionValue[i+offset]) * Df 
                Tree = c(Tree, OptionValue[i+offset]) } } }
                
    # American Type:
    if (TypeFlag == "ca" || TypeFlag == "pa") {
        for (j in (n-1):0) { 
            Tree <-c(Tree, rep(0, times=n-j))
            for (i in 0:j) { 
                OptionValue[i+offset] = 
                max((z * (S*u^i*d^(abs(i-j)) - X)), 
                        (p*OptionValue[i+1+offset] + 
                (1-p)*OptionValue[i+offset]) * Df ) 
                Tree = c(Tree, OptionValue[i+offset]) } } } 
                
    # Tree-Matrix of form (here n=4):
    # x x x x
    # . x x x
    # . . x x
    # . . . x
    Tree = matrix(rev(Tree), byrow = FALSE, ncol = n+1)
    
    # Return Value:
    Tree
}


# ------------------------------------------------------------------------------


BinomialTreePlot = 
function(BinomialTreeValues, dx = -0.025, dy = 0.4, cex = 1, digits = 2, ...) 
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Plots the binomial tree of the Cox-Ross-Rubinstein
    #   binomial tree model.
    
    # Example:
    #   par(mfrow=c(1,1))
    #   Tree = BinomialTree("a", "p", 100, 95, 0.5, 0.08, 0.08, 0.3, 5)
    #   print(round(Tree, digits=3))
    #   BinomialTreePlot(Tree, main="American Put Option")

    # FUNCTION:
    
    # Tree:
    Tree = round(BinomialTreeValues, digits=digits)
    depth = ncol(Tree)
    plot(x=c(1,depth), y=c(-depth+1, depth-1), type="n", col=0, ...)
    points(x=1, y=0)
    text(1+dx, 0+dy, deparse(Tree[1,1]), cex=cex)
    for (i in 1:(depth-1) ) {
        y = seq(from=-i, by=2, length=i+1)
        x = rep(i, times=length(y))+1
        points(x, y, col=1) 
        for (j in 1:length(x))
            text(x[j]+dx, y[j]+dy, deparse(Tree[length(x)+1-j,i+1]), cex=cex)   
        y = (-i):i
        x = rep(c(i+1,i), times=2*i)[1:length(y)]
        lines(x, y, col=2)
    }
    
    # Return Value:
    invisible()
}


# --- 3.1.2 --------------------------------------------------------------------


# Options on a Stock Paying a Known Dividend Yield
# not yet implemented


# --- 3.1.3 --------------------------------------------------------------------


# BarrierBinomialTree
# not yet implemented


# --- 3.1.4 --------------------------------------------------------------------


# ConvertibleBond
# not yet implemented


# --- 3.2 ----------------------------------------------------------------------


# TrinomialTree
# not yet implemented


# --- 3.3 ----------------------------------------------------------------------


# ThreeDimensionalBinomialTree
# PayoffFunction
# not yet implemented


# --- 3.4.1 --------------------------------------------------------------------


# ImpliedBinomialTree
# not yet implemented


# --- 3.4.2 --------------------------------------------------------------------


# ImpliedTrinomialTree
# not yet implemented


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                       DESCRIPTION:
#  Multiple Exercies Options:           
#   ExecutiveStockOption            Executive Stock Option
#   ForwardStartOption              Forward Start Option
#   RatchetOption                   Ratchet [Compound] Option
#   TimeSwitchOption                Time Switch Option
#   SimpleChooserOption             Simple Chooser Option
#   ComplexChooserOption            Complex Chooser Option
#   OptionOnOption                  Options On Options
#   HolderExtendibleOption          Holder Extendible Option
#   WriterExtendibleOption          Writer Extendible Option
################################################################################


ExecutiveStockOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma, lambda)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Executive stock options

    # References:
    #   Jennergren and Naslund (1993) 
    #   Haug, Chapter 2.1
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Calculate Price:
    ExecutiveStock = exp (-lambda * Time) * 
        GBSOption(TypeFlag = TypeFlag, S = S, X = X, Time = Time, 
            r = r, b = b, sigma = sigma)$price

    # Return Value:
    option = list(
        price = ExecutiveStock, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


ForwardStartOption =  
function(TypeFlag = c("c", "p"), S, alpha, time1, Time2, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Forward Start Options

    # References:
    #   Rubinstein (1990)
    #   Haug, Chapter 2.2
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Compute Settings:
    Time = time1
    time = Time2
    
    # Compute Price:
    ForwardStart = S * exp ((b - r) * time ) *
        GBSOption(TypeFlag, S = 1, X = alpha, Time = Time-time, 
            r = r, b = b, sigma = sigma)$price
    
    # Return Value:
    option = list(
        price = ForwardStart, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


RatchetOption = 
function(TypeFlag = c("c", "p"), S, alpha, time1, Time2, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Ratchet Option, 
    #   other names are MovingStrikeOption or CliquetOption

    # References:
    #   Haug, Chapter 2.3
    
    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Calculate Price
    Ratchet = 0 
    for ( i in 1:length(Time2) ) {
        Ratchet = Ratchet +
            ForwardStartOption(TypeFlag = TypeFlag, 
                S = S, alpha = alpha, time1 = time1[i], Time2 = Time2[i], 
                r = r, b = b, sigma = sigma)$price }
            
    # Return Value:
    option = list(
        price = Ratchet, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


TimeSwitchOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma, A, m, dt)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Discrete time switch options

    # References:
    #   Pechtl (1995)
    #   Haug, Chapter 2.4

    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Compute Settings:
    n = Time / dt
    Sum = 0
    if (TypeFlag == "c") Z = +1
    if (TypeFlag == "p") Z = -1
    
    # Calculate Price:
    Sum = 0
    for (I in (1:n)) {
        d = (log(S/X) + (b - sigma^2/2) * I * dt) / (sigma * sqrt(I * dt))
        Sum = Sum + CND (Z * d) * dt }
    TimeSwitch = A * exp (-r * Time) * Sum + dt * A * exp(-r * Time) * m
            
    # Return Value:
    option = list(
        price = TimeSwitch, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


SimpleChooserOption = 
function(S, X, time1, Time2, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Simple Chooser Options

    # References:
    #   Rubinstein (1991)
    #   Haug, Chapter 2.5.1

    # FUNCTION:
    
    # Compute Settings:
    d = (log(S/X) + (b + sigma ^ 2 / 2) * Time2) / (sigma * sqrt(Time2))
    y = (log(S/X) + b * Time2 + sigma ^ 2 * time1 / 2) / 
        (sigma * sqrt(time1))
    
    # Calculate Price:
    SimpleChooser = S * exp ((b - r) * Time2) * CND(d) -
        X * exp(-r * Time2) * CND(d - sigma * sqrt(Time2)) -
        S * exp ((b - r) * Time2) * CND(-y) +
        X * exp(-r * Time2) * CND(-y + sigma * sqrt(time1))   
            
    # Return Value:
    option = list(
        price = SimpleChooser, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


ComplexChooserOption = 
function(S, Xc, Xp, Time, Timec, Timep, r, b, sigma, doprint = FALSE)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Complex Chooser Options

    # References:
    #   Haug, Chapter 2.5.2
    
    # FUNCTION:
    
    # Compute Settings:
    Tc = Timec
    Tp = Timep
    
    # Calculate Price:
    CriticalValueChooser = 
    function(S, Xc, Xp, Time, Tc, Tp, r, b, sigma){
        Sv = S
        ci = GBSOption("c", Sv, Xc, Tc - Time, r, b, sigma)$price
        Pi = GBSOption("p", Sv, Xp, Tp - Time, r, b, sigma)$price
        dc = GBSGreeks("Delta", "c", Sv, Xc, Tc - Time, r, b, sigma)
        dp = GBSGreeks("Delta", "p", Sv, Xp, Tp - Time, r, b, sigma)
        yi = ci - Pi
        di = dc - dp
        epsilon = 0.001
        # Newton-Raphson:
        while (abs(yi) > epsilon) {
            Sv = Sv - (yi) / di
            ci = GBSOption("c", Sv, Xc, Tc - Time, r, b, sigma)$price
            Pi = GBSOption("p", Sv, Xp, Tp - Time, r, b, sigma)$price
            dc = GBSGreeks("Delta", "c", Sv, Xc, Tc - Time, r, b, sigma)
            dp = GBSGreeks("Delta", "p", Sv, Xp, Tp - Time, r, b, sigma)
            yi = ci - Pi
            di = dc - dp }
        result = Sv
        result}
        
    # Complex chooser options:
    I = CriticalValueChooser (S, Xc, Xp, Time, Tc, Tp, r, b, sigma)
    if (doprint) {
        cat("\nCritical Value:\n")
        print(I) 
        cat ("\n")}
    d1 = (log(S / I) + (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
    d2 = d1 - sigma * sqrt (Time)
    y1 = (log(S / Xc) + (b + sigma ^ 2 / 2) * Tc) / (sigma * sqrt(Tc))
    y2 = (log(S / Xp) + (b + sigma ^ 2 / 2) * Tp) / (sigma * sqrt(Tp))
    rho1 = sqrt (Time / Tc)
    rho2 = sqrt (Time / Tp)
    
    ComplexChooser = S * exp ((b - r) * Tc) * CBND(d1, y1, rho1) -
        Xc * exp(-r * Tc) * CBND(d2, y1 - sigma * sqrt(Tc), rho1) -
        S * exp((b - r) * Tp) * CBND(-d1, -y2, rho2) +
        Xp * exp(-r * Tp) * CBND(-d2, -y2 + sigma * sqrt(Tp), rho2)
            
    # Return Value:
    option = list(
        price = ComplexChooser, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


OptionOnOption = 
function(TypeFlag = c("cc", "cp", "pc", "pp"), S, X1, X2, time1, Time2, r, 
b, sigma, doprint = FALSE)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Option on Option

    # References:
    #   Geske (1977), Geske (1979b), Hodges and Selby (1987), 
    #   Rubinstein (1991a) et al.
    #   Haug, Chpater 2.6

    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Compute Settings:
    Time = time1
    time = Time2
    
    # Internal Function:
    CriticalValueOptionOnOption = 
    function(TypeFlag, X1, X2, Time, r, b, sigma) { 
        # Calculation of critical price options on options
        Si = X1
        ci = GBSOption(TypeFlag, Si, X1, Time, r, b, sigma)$price
        di = GBSGreeks("Delta", TypeFlag, Si, X1, Time, r, b, sigma)
        epsilon = 0.000001
        # Newton-Raphson algorithm:
        while (abs(ci - X2) > epsilon) {
                Si = Si - (ci - X2) / di
                ci = GBSOption(TypeFlag, Si, X1, Time, r, b, sigma)$price
                di = GBSGreeks("Delta", TypeFlag, Si, X1, Time, r, b, sigma) }
        result = Si
        result }
            
    # Option On Option:
    T2 = Time
    t1 = time
    TypeFlag2 = "p"
    if (TypeFlag == "cc" || TypeFlag == "pc") TypeFlag2 = "c"
    I = CriticalValueOptionOnOption(TypeFlag2, X1, X2, T2-t1, r, b, sigma)
    if (doprint) { cat("\nCriticalValue: ", I, "\n") }
    rho = sqrt (t1 / T2)
    y1 = (log(S / I) + (b + sigma ^ 2 / 2) * t1) / (sigma * sqrt(t1))
    y2 = y1 - sigma * sqrt (t1)
    z1 = (log(S / X1) + (b + sigma ^ 2 / 2) * T2) / (sigma * sqrt(T2))
    z2 = z1 - sigma * sqrt (T2)
    if (TypeFlag == "cc") 
    OptionOnOption = S * exp ((b - r) * T2) * CBND(z1, y1, rho) -
        X1 * exp(-r * T2) * CBND(z2, y2, rho) - X2 * exp(-r * t1) * 
        CND(y2)
    if (TypeFlag == "pc") 
    OptionOnOption = X1 * exp (-r * T2) * CBND(z2, -y2, -rho) -
        S * exp((b - r) * T2) * CBND(z1, -y1, -rho) + X2 * 
        exp(-r * t1) * CND(-y2)
    if (TypeFlag == "cp") 
    OptionOnOption = X1 * exp (-r * T2) * CBND(-z2, -y2, rho) -
        S * exp((b - r) * T2) * CBND(-z1, -y1, rho) - X2 * 
        exp(-r * t1) * CND(-y2) 
    if (TypeFlag == "pp") 
    OptionOnOption = S * exp ((b - r) * T2) * CBND(-z1, y1, -rho) -
        X1 * exp(-r * T2) * CBND(-z2, y2, -rho) + exp(-r * t1) *
        X2 * CND(y2) 
            
    # Return Value:
    option = list(
        price = OptionOnOption, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


HolderExtendibleOption = 
function(TypeFlag = c("c", "p"), S, X1, X2, time1, Time2, r, b, sigma, A)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Options that can be extended by the Holder

    # References:
    #   Haug, Chapter 2.7.1

    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    
    # Calculate Price:
    HolderExtendible = NA
    if (TypeFlag == "c") {
        HolderExtendible = max(c(S-X1, GBSOption(TypeFlag = "c", S = S, 
            X = X2, Time = Time2-time1, r = r, b = b, sigma = sigma)$price - 
            A, 0)) }
    if (TypeFlag == "p") {
        HolderExtendible = max(c(X1-S, GBSOption(TypeFlag = "p", S = S, 
            X = X2, Time = Time2-time1, r = r, b = b, sigma = sigma)$price - 
            A, 0)) }
            
    # Return Value:
    option = list(
        price = HolderExtendible, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


WriterExtendibleOption = 
function(TypeFlag = c("c", "p"), S, X1, X2, time1, Time2, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Writer Extendible Options

    # References:
    #   Haug, Chapter 2.7.2

    # FUNCTION:
    
    # Settings:
    TypeFlag = TypeFlag[1]
    rho = sqrt (time1 / Time2)
    z1 = (log(S/X2) + (b + sigma^2 / 2) * Time2) / (sigma * sqrt(Time2))
    z2 = (log(S/X1) + (b + sigma^2 / 2) * time1) / (sigma * sqrt(time1))
        
    # Calculate Price:
    if (TypeFlag == "c") 
        WriterExtendible = 
            GBSOption(TypeFlag, S, X1, time1, r, b, sigma)$price +
                S * exp((b - r) * Time2) * CBND(z1, -z2, -rho) - 
                X2 * exp(-r * Time2) * CBND(z1 - sqrt(sigma^2 * Time2), 
                -z2 + sqrt(sigma^2 * time1), -rho)
    if (TypeFlag == "p") 
        WriterExtendible = 
            GBSOption(TypeFlag, S, X1, time1, r, b, sigma)$price +
                X2 * exp(-r * Time2) * CBND(-z1 + sqrt(sigma^2 * Time2), 
                z2 - sqrt(sigma^2 * time11), -rho) - 
                S * exp((b - r) * Time2) * CBND(-z1, z2, -rho)
            
    # Return Value:
    option = list(
        price = WriterExtendible, 
        call = match.call() )
    class(option) = "option"
    option
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                       DESCRIPTION:
# Multiple Asset Options:
#   TwoAssetCorrelationOption       Two Asset Correlation Option
#   ExchangeOneForAnotherOption     Exchane One For Another Option  
#   ExchangeOnExchangeOption        Exchange Exchange Option
#   EuropeanExchangeOption          European Exchange Optionn
#   AmericanExchangeOption          American Exchange Option
#   TwoRiskyAssetsOption            Option On The MinMax
#   SpreadApproxOption              Spread Approximated Option              
################################################################################


TwoAssetCorrelationOption = 
function(TypeFlag = c("c", "p"), S1, S2, X1, X2, Time, r, b1, b2, 
sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Two asset correlation options

    # References:
    #   Haug, Chapter 2.8.1

    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    y1 = (log(S1/X1) + (b1 - sigma1^2 / 2) * Time) / (sigma1*sqrt(Time))
    y2 = (log(S2/X2) + (b2 - sigma2^2 / 2) * Time) / (sigma2*sqrt(Time))
    
    # Calculate Call and Put:
    if (TypeFlag == "c") 
    TwoAssetCorrelation = S2 * exp ((b2 - r) * Time) *
        CBND(y2 + sigma2 * sqrt(Time), y1 + rho * sigma2 * sqrt(Time), rho) -
        X2 * exp (-r * Time) * CBND(y2, y1, rho) 
    if (TypeFlag == "p") 
    TwoAssetCorrelation = X2 * exp (-r * Time) * CBND(-y2, -y1, rho) -
        S2 * exp ((b2 - r) * Time) *
        CBND(-y2 - sigma2 * sqrt(Time), -y1 - rho * sigma2 * sqrt(Time), rho)
    
    # Return Value:
    option = list(
        price = TwoAssetCorrelation, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


EuropeanExchangeOption = 
function(S1, S2, Q1, Q2, Time, r, b1, b2, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           
     
    # Description:
    #   Exchange-One-Asset-for-Another-Asset options -
    #   European option to exchange one asset for another
     
    # References:
    #   Haug, Chapter 2.8.2 (European)

    # FUNCTION:
    
    # Compute Settings:
    sigma = sqrt (sigma1 ^ 2 + sigma2 ^ 2 - 2 * rho * sigma1 * sigma2)
    d1 = ((log(Q1*S1/(Q2 * S2)) + (b1-b2+sigma^2/2)*Time)/(sigma*sqrt(Time)))
    d2 = d1 - sigma * sqrt (Time)
    
    # calculate Price:
    EuropeanExchange = Q1 * S1 * exp ((b1 - r) * Time) * CND(d1) -
        Q2 * S2 * exp((b2 - r) * Time) * CND(d2)
    
    # Return Value:
    option = list(
        price = EuropeanExchange, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


AmericanExchangeOption = 
function(S1, S2, Q1, Q2, Time, r, b1, b2, sigma1, sigma2, rho, doprint = FALSE)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Exchange-One-Asset-for-Another-Asset options -
    #   American option to exchange one asset for another

    # References:
    #   Haug, Chapter 2.8.2 (American)

    # FUNCTION:
    
    # Compute Settings:
    sigma = sqrt(sigma1^2 + sigma2^2 - 2 * rho * sigma1 * sigma2)
    
    # Calculate Price:
    AmericanExchange = BSAmericanApproxOption("c", Q1*S1, Q2*S2, 
        Time, r-b2, b1-b2, sigma)
        
    # Print Trigger Price:
    if (doprint) {cat("\nTriggerPrice: ", AmericanExchange$TriggerPrice, "\n")}
    
    # Return Value:
    option = list(
        price = AmericanExchange$Premium, 
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


ExchangeOnExchangeOption = 
function(TypeFlag = c("1", "2", "3", "4"), S1, S2, Q, time1, Time2, r, 
b1, b2, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           
  
    # Description:
    #   Exchange-One-Asset-for-Another-Asset options -

    # References:
    #   Haug, Chapter 2.8.3 
    
    # FUNCTION:
    
    # Define Functions:
    TypeFlag = TypeFlag[1]
    q = Q
    
    # Third:
    # To run under SPlus we require "<<-"
    CriticalPart3 <<- function(id, I, time1, Time2, v) {    
        if (id == 1) {
            z1 = (log(I)+v^2/2*(Time2 - time1))/(v*sqrt(Time2-time1))
            z2 = (log(I)-v^2/2*(Time2 - time1))/(v*sqrt(Time2-time1))
            CriticalPart3 = I * CND(z1) - CND(z2) }     
        if (id == 2) {
            z1 = (-log(I)+v^2/2*(Time2-time1))/(v*sqrt(Time2-time1))
            z2 = (-log(I)-v^2/2*(Time2-time1))/(v*sqrt(Time2-time1))
            CriticalPart3 = CND(z1) - I * CND(z2) }
        CriticalPart3 }   
        
    # Second:
    CriticalPart2 <<- function(id, I, time1, Time2, v) {
        if (id == 1) {
            z1 = (log(I)+v^2/2*(Time2-time1))/(v*sqrt(Time2-time1))
            CriticalPart2 = CND(z1) }      
        if (id == 2) {
            z2 = (-log(I)-v^2/2*(Time2-time1))/(v*sqrt(Time2-time1))
            CriticalPart2 = -CND(z2) }
        CriticalPart2 }
        
    # Numerical search algorithm to find critical price I
    CriticalPrice = function(id, I1, time1, Time2, v, q) {
        Ii = I1
        yi = CriticalPart3(id, Ii, time1, Time2, v)
        # cat("\nCriticalPart3: ", yi)
        di = CriticalPart2(id, Ii, time1, Time2, v)
        # cat("\nCriticalPart2: ", di)
        epsilon = 0.00001
        while (abs(yi - q) > epsilon) {
            Ii = Ii - (yi - q) / di
            yi = CriticalPart3(id, Ii, time1, Time2, v)
            # cat("\nCriticalPart3: ", yi)
            di = CriticalPart2(id, Ii, time1, Time2, v)
            # cat("\nCriticalPart2: ", di) 
            }
        CriticalPrice = Ii
        CriticalPrice }
    
    # Compute:
    v = sqrt(sigma1 ^ 2 + sigma2 ^ 2 - 2 * rho * sigma1 * sigma2)
    I1 = S1 * exp((b1 - r) * (Time2 - time1)) / 
        (S2 * exp((b2 - r) * (Time2 - time1)))   
    if (TypeFlag == "1" || TypeFlag == "2") {
        id = 1 }
    else {
        id = 2 }  
    I = CriticalPrice(id, I1, time1, Time2, v, q)
    
    d1 = (log(S1 / (I * S2)) + (b1 - b2 + v ^ 2 / 2) * time1) / 
        (v * sqrt(time1))
    d2 = d1 - v * sqrt(time1)
    d3 = (log((I * S2) / S1) + (b2 - b1 + v ^ 2 / 2) * time1) / 
        (v * sqrt(time1))
    d4 = d3 - v * sqrt(time1)
    y1 = (log(S1 / S2) + (b1 - b2 + v ^ 2 / 2) * Time2) / (v * sqrt(Time2))
    y2 = y1 - v * sqrt(Time2)
    y3 = (log(S2 / S1) + (b2 - b1 + v ^ 2 / 2) * Time2) / (v * sqrt(Time2))
    y4 = y3 - v * sqrt(Time2)
    
    # Calculate Price:
    if (TypeFlag == "1")
        ExchangeOnExchange = -S2 * exp((b2 - r) * Time2) * 
            CBND(d2, y2, sqrt(time1/Time2)) + S1 * exp((b1-r) * Time2) * 
            CBND(d1, y1, sqrt(time1/Time2)) - q * S2 * exp((b2-r) * time1) * 
            CND(d2)
    if (TypeFlag == "2")
        ExchangeOnExchange = S2 * exp((b2 - r) * Time2) * 
            CBND(d3, y2, -sqrt(time1/Time2)) - S1 * exp((b1-r) * Time2) * 
            CBND(d4, y1, -sqrt(time1/Time2)) + q * S2 * exp((b2 - r) * time1) * 
            CND(d3)
    if (TypeFlag == "3")
        ExchangeOnExchange = S2 * exp((b2 - r) * Time2) * 
            CBND(d3, y3, sqrt(time1/Time2)) - S1 * exp((b1-r) * Time2) * 
            CBND(d4, y4, sqrt(time1/Time2)) - q * S2 * exp((b2-r) * time1) * 
            CND(d3)
    if (TypeFlag == "4")
        ExchangeOnExchange = -S2 * exp((b2 - r) * Time2) * 
            CBND(d2, y3, -sqrt(time1/Time2)) + S1 * exp((b1-r) * Time2) * 
            CBND(d1, y4, -sqrt(time1/Time2)) + q * S2 * exp((b2-r) * time1) * 
            CND(d2)   
    
    # Return Value:
    option = list(
        price = ExchangeOnExchange, 
        call = match.call() )
    class(option) = "option"
    option  
}


# ------------------------------------------------------------------------------


TwoRiskyAssetsOption = 
function(TypeFlag = c("cmin", "cmax", "pmin", "pmax"), S1, S2, X, Time, 
r, b1, b2, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Option on two risky assets
    
    # References:
    #   Haug, Chapter 2.8.4 
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    v = sqrt(sigma1 ^ 2 + sigma2 ^ 2 - 2 * rho * sigma1 * sigma2)
    rho1 = (sigma1 - rho * sigma2) / v
    rho2 = (sigma2 - rho * sigma1) / v
    d = (log(S1 / S2) + (b1 - b2 + v ^ 2 / 2) * Time) / (v * sqrt(Time))
    y1 = (log(S1 / X) + (b1 + sigma1 ^ 2 / 2) * Time) / (sigma1 * sqrt(Time))
    y2 = (log(S2 / X) + (b2 + sigma2 ^ 2 / 2) * Time) / (sigma2 * sqrt(Time))
    
    # Calculate Price:
    OnTheMaxMin = NA
    if (TypeFlag == "cmin")
        OnTheMaxMin = S1 * exp((b1 - r) * Time) * 
            CBND(y1, -d, -rho1) + S2 * exp((b2 - r) * Time) * 
            CBND(y2, d - v * sqrt(Time), -rho2) - X * exp(-r * Time) * 
            CBND(y1 - sigma1 * sqrt(Time), y2 - sigma2 * sqrt(Time), rho)
    if (TypeFlag == "cmax")
        OnTheMaxMin = S1 * exp((b1 - r) * Time) * 
            CBND(y1, d, rho1) + S2 * exp((b2 - r) * Time) * 
            CBND(y2, -d + v * sqrt(Time), rho2) - X * exp(-r * Time) * 
            (1 - CBND(-y1 + sigma1*sqrt(Time), -y2 + sigma2 * sqrt(Time), rho))
    if (TypeFlag == "pmin")
        OnTheMaxMin = X * exp(-r * Time) - S1 * exp((b1 - r) * Time) + 
            EuropeanExchangeOption(S1, S2, 1, 1, Time, r, b1, b2, 
                sigma1, sigma2, rho)$price + 
            TwoRiskyAssetsOption("cmin", S1, S2, X, Time, r, b1, b2, 
                sigma1, sigma2, rho)$price
    if (TypeFlag == "pmax")
        OnTheMaxMin = X * exp(-r * Time) - S2 * exp((b2 - r) * Time) - 
            EuropeanExchangeOption(S1, S2, 1, 1, Time, r, b1, b2, 
                sigma1, sigma2, rho)$price + 
            TwoRiskyAssetsOption("cmax", S1, S2, X, Time, r, b1, b2, 
                sigma1, sigma2, rho)$price   
    
    # Return Value:
    option = list(
        price = OnTheMaxMin, 
        call = match.call() )
    class(option) = "option"
    option    
}


# ------------------------------------------------------------------------------


SpreadApproxOption = 
function(TypeFlag = c("c", "p"), S1, S2, X, Time, r, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Spread Option Approximation

    # References:
    #   Haug, Chapter 2.8.5
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    F1 = S1
    F2 = S2
    sigma = sqrt(sigma1 ^ 2 + (sigma2 * F2 / (F2 + X)) ^ 2 - 2 * rho * 
        sigma1 * sigma2 * F2 / (F2 + X))
    FF = F1 / (F2 + X) 
    
    # Calculate Price
    SpreadApproximation = 
        GBSOption(TypeFlag, FF, 1, Time, r, 0, sigma)$price * (F2 + X)   
    
    # Return Value:
    option = list(
        price = SpreadApproximation, 
        call = match.call() )
    class(option) = "option"
    option
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                           DESCRIPTION:
# Lookback Options:
#   FloatingStrikeLookbackOption        Floating Strike Lookback Option
#   FixedStrikeLookbackOption           Fixed Strike Lookback Option
#   PTFloatingStrikeLookbackOption      Partial Floating Strike LB Option
#   PTFixedStrikeLookbackOption         Partial Fixed Strike LB Option  
#   ExtremeSpreadOption                 Extreme Spread Option
################################################################################


FloatingStrikeLookbackOption = 
function(TypeFlag = c("c", "p"), S, SMinOrMax, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Floating strike lookback options

    # References:
    #   Haug, Chapter 2.9.1

    # FUNCTION:
    
    # Comute Settungs:
    TypeFlag = TypeFlag[1]
    if (TypeFlag == "c") m = SMinOrMax # Min
    if (TypeFlag == "p") m = SMinOrMax # Max
    a1 = (log(S / m) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    a2 = a1 - sigma * sqrt (Time)
    
    # Calculate Call and Put:
    if (TypeFlag == "c") 
        FloatingStrikeLookback = S * exp ((b - r) * Time) * CND(a1) -
            m * exp(-r * Time) * CND(a2) + exp (-r * Time) *
            sigma^2 / (2 * b) * S * ((S / m)^(-2 * b / sigma^2) *
            CND(-a1 + 2 * b / sigma * sqrt(Time)) - exp(b * Time) * CND(-a1))
    if (TypeFlag == "p") 
        FloatingStrikeLookback = m * exp (-r * Time) * CND(-a2) -
            S * exp((b - r) * Time) * CND(-a1) + exp (-r * Time) * 
            sigma^2 / (2 * b) * S * (-(S / m)^(-2 * b / sigma^2) *
            CND(a1 - 2 * b / sigma * sqrt(Time)) + exp(b * Time) * CND(a1)) 
    
    # Return Value:
    option = list(
        price = FloatingStrikeLookback, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


FixedStrikeLookbackOption = 
function(TypeFlag = c("c", "p"), S, SMinOrMax, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Fixed strike lookback options

    # References:
    #   Haug, Chapter 2.9.2

    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    if (TypeFlag == "c") m = SMinOrMax
    if (TypeFlag == "p") m = SMinOrMax
    d1 = (log(S / X) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    d2 = d1 - sigma * sqrt (Time)
    e1 = (log(S / m) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    e2 = e1 - sigma * sqrt (Time)
    
    # Calculate Call and Put:
    if (TypeFlag == "c" && X > m) 
            FixedStrikeLookback = S * exp ((b - r) * Time) * CND(d1) -
                X * exp(-r * Time) * CND(d2) +
                S * exp (-r * Time) * sigma^2 / (2 * b) * 
            (-(S / X)^(-2 * b / sigma^2) *
                CND(d1 - 2 * b / sigma * sqrt(Time)) + exp(b * Time) * CND(d1))
    if (TypeFlag == "c" && X <= m)  
            FixedStrikeLookback = exp (-r * Time) * (m - X) +
                S * exp((b-r) * Time) * CND(e1) - exp(-r * Time) * m * CND(e2) +
                S * exp (-r * Time) * sigma^2 / (2 * b) * 
            (-(S / m)^(-2 * b / sigma^2) *
                CND(e1 - 2 * b / sigma * sqrt(Time)) + exp(b * Time) * CND(e1))
    if (TypeFlag == "p" && X < m) 
            FixedStrikeLookback = -S * exp ((b - r) * Time) * CND(-d1) +
                X * exp(-r * Time) * CND(-d1 + sigma * sqrt(Time)) +
                S * exp (-r * Time) * sigma^2 / (2 * b) * 
            ((S / X)^(-2 * b / sigma^2) *
                CND(-d1 + 2 * b / sigma * sqrt(Time)) - exp(b*Time) * CND(-d1))
    if (TypeFlag == "p" && X >= m) 
            FixedStrikeLookback = exp (-r * Time) * (X - m) -
                S * exp((b - r) * Time) * CND(-e1) + 
            exp(-r * Time) * m * CND(-e1 + sigma * sqrt(Time)) +
                exp (-r * Time) * sigma^2 / (2 * b) * S * 
            ((S / m)^(-2 * b / sigma^2) *
                CND(-e1 + 2 * b / sigma * sqrt(Time)) - exp(b*Time) * CND(-e1))
    
    # Return Value:
    option = list(
        price = FixedStrikeLookback, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


PTFloatingStrikeLookbackOption = 
function(TypeFlag = c("c", "p"), S, SMinOrMax, time1, Time2, r, b, 
sigma, lambda)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Partial-time floating strike lookback options

    # References:
    #   Haug, Chapter 2.9.3
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    T2 = Time2
    t1 = time1
    if (TypeFlag == "c") m = SMinOrMax
    if (TypeFlag == "p") m = SMinOrMax
    d1 = (log(S / m) + (b + sigma^2 / 2) * T2) / (sigma * sqrt(T2))
    d2 = d1 - sigma * sqrt (T2)
    e1 = (b + sigma^2 / 2) * (T2 - t1) / (sigma * sqrt(T2 - t1))
    e2 = e1 - sigma * sqrt (T2 - t1)
    f1 = (log(S / m) + (b + sigma^2 / 2) * t1) / (sigma * sqrt(t1))
    f2 = f1 - sigma * sqrt (t1)
    g1 = log (lambda) / (sigma * sqrt(T2))
    g2 = log (lambda) / (sigma * sqrt(T2 - t1))
        
    # Calculate Call and Puts:
    if (TypeFlag == "c") {
        part1 = S * exp ((b - r) * T2) * CND(d1 - g1) -
            lambda * m * exp(-r * T2) * CND(d2 - g1)
        part2 = exp (-r * T2) * sigma^2 /
            (2 * b) * lambda * S * ((S / m)^(-2 * b / sigma^2) *
            CBND(-f1 + 2*b*sqrt(t1) / sigma, -d1 + 2 * b * sqrt(T2) / sigma -
            g1, sqrt(t1 / T2)) - exp (b * T2) * lambda^(2 * b / sigma^2) *
            CBND(-d1 - g1, e1 + g2, -sqrt(1 - t1 / T2))) +
            S * exp ((b - r)*T2) * CBND(-d1 + g1, e1 - g2, -sqrt(1 - t1 / T2))
        part3 = exp (-r*T2) * lambda * m * CBND(-f2, d2 - g1, -sqrt(t1 / T2)) -
            exp (-b * (T2-t1)) * exp((b - r)*T2) * (1 + sigma^2 / (2 * b)) *
            lambda * S * CND(e2 - g2) * CND(-f1) }
    if (TypeFlag == "p") {
        part1 = lambda * m * exp (-r * T2) * CND(-d2 + g1) -
            S * exp((b - r) * T2) * CND(-d1 + g1)
        part2 = -exp (-r * T2) * sigma^2 /
            (2 * b) * lambda * S * ((S / m)^(-2 * b / sigma^2) *
            CBND(f1 - 2 * b * sqrt(t1) / sigma, d1 - 2 * b * sqrt(T2) / sigma +
            g1, sqrt(t1 / T2)) - exp (b * T2) * lambda^(2 * b / sigma^2) *
            CBND(d1 + g1, -e1 - g2, -sqrt(1 - t1 / T2))) -
            S * exp ((b - r)*T2) * CBND(d1 - g1, -e1 + g2, -sqrt(1 - t1 / T2))
        part3 = -exp (-r*T2) * lambda*m * CBND(f2, -d2 + g1, -sqrt(t1 / T2)) +
            exp (-b * (T2-t1)) * exp((b - r)*T2) * (1 + sigma^2 / (2 * b)) *
            lambda * S * CND(-e2 + g2) * CND(f1) }
    PartialFloatLookback = part1 + part2 + part3
    
    # Return Value:
    option = list(
        price = PartialFloatLookback, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


PTFixedStrikeLookbackOption = 
function(TypeFlag = c("c", "p"), S, X, time1, Time2, r, b, sigma) 
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Partial Time Fixed Strike Lookback Option 

    # References:
    #   Haug, Chapter 2.9.4
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    d1 = (log(S / X) + (b + sigma^2 / 2) * Time2) / 
        (sigma * sqrt(Time2))
    d2 = d1 - sigma * sqrt(Time2)
    e1 = ((b + sigma^2 / 2) * (Time2 - time1)) / 
        (sigma * sqrt(Time2 - time1))
    e2 = e1 - sigma * sqrt(Time2 - time1)
    f1 = (log(S / X) + (b + sigma^2 / 2) * time1) / (sigma * sqrt(time1))
    f2 = f1 - sigma * sqrt(time1)
    
    # Calculate Call and Put:
    if (TypeFlag == "c") {
        PartialFixedLB = S * exp((b - r) * Time2) * 
            CND(d1) - exp(-r * Time2) * X * 
            CND(d2) + S * exp(-r * Time2) * sigma^2 / (2 * b) * 
                (-(S / X)^(-2 * b / sigma^2) * 
            CBND(d1 - 2 * b * sqrt(Time2) / sigma, -f1 + 2 * b * sqrt(time1) / 
                sigma, -sqrt(time1 / Time2)) + exp(b * Time2) * CBND(e1, d1, 
                sqrt(1 - time1 / Time2))) - S * exp((b - r) * Time2) * 
            CBND(-e1, d1, -sqrt(1 - time1 / Time2)) - 
                X * exp(-r * Time2) * 
            CBND(f2, -d2, -sqrt(time1 / Time2)) + exp(-b * (Time2 - time1)) * 
                (1 - sigma^2 / (2 * b)) * S * exp((b - r) * Time2) * 
            CND(f1) * 
            CND(-e2) }
   if (TypeFlag == "p") {
        PartialFixedLB = X * exp(-r * Time2) * 
            CND(-d2) - S * exp((b - r) * Time2) * 
            CND(-d1) + S * exp(-r * Time2) * sigma^2 / (2 * b) * 
                ((S / X)^(-2 * b / sigma^2) * CBND(-d1 + 2 * b * 
                sqrt(Time2) / sigma, f1 - 2 * b * sqrt(time1) / sigma, 
                -sqrt(time1 / Time2)) - exp(b * Time2) * 
            CBND(-e1, -d1, sqrt(1 - time1 / Time2))) + 
                S * exp((b - r) * Time2) * 
            CBND(e1, -d1, -sqrt(1 - time1 / Time2)) + X * exp(-r * Time2) * 
            CBND(-f2, d2, -sqrt(time1 / Time2)) - exp(-b * (Time2 - time1)) * 
                (1 - sigma^2 / (2 * b)) * S * exp((b - r) * Time2) * 
            CND(-f1) * 
            CND(e2) }
    
    # Return Value:
    option = list(
        price = PartialFixedLB, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


ExtremeSpreadOption = 
function(TypeFlag = c("c", "p", "cr", "pr"), S, SMin, SMax, time1, Time2, 
r, b, sigma)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Extreme Spread Option 

    # References:
    #   Haug, Chapter 2.9.5
   
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    v = sigma
    Time = Time2
    if (TypeFlag == "c"  || TypeFlag == "cr") { eta = +1 }
    if (TypeFlag == "p"  || TypeFlag == "pr") { eta = -1 }
    if (TypeFlag == "c"  || TypeFlag == "p")  { kappa = +1 }
    if (TypeFlag == "cr" || TypeFlag == "pr") { kappa = -1 }          
    if (kappa * eta == +1) { Mo = SMax }
    if (kappa * eta == -1) { Mo = SMin }    
    mu1 = b - v^2 / 2
    mu = mu1 + v^2
    m = log(Mo/S)   
    ExtremeSpread = NA
    
    # Extreme Spread Option:
    if (kappa == 1) {
    ExtremeSpread = eta * (S * exp((b - r) * Time) * (1 + v^2 / (2 * b)) * 
        CND(eta * (-m + mu * Time) / (v*sqrt(Time))) - 
            exp(-r * (Time - time1)) * S * exp((b - r) * Time) * 
            (1 + v^2 / (2 * b)) * 
        CND(eta * (-m + mu * time1) / (v*sqrt(time1))) + exp(-r * Time) * 
            Mo * CND(eta * (m - mu1 * Time) / (v*sqrt(Time))) - 
            exp(-r * Time) * Mo * v^2 / (2 * b) * 
            exp(2 * mu1 * m / v^2) * 
        CND(eta * (-m - mu1 * Time) / (v*sqrt(Time))) - exp(-r * Time) * Mo * 
        CND(eta * (m - mu1 * time1) / (v*sqrt(time1))) + exp(-r * Time) * 
            Mo * v^2 / (2 * b) * exp(2 * mu1 * m / v^2) * 
        CND(eta * (-m - mu1 * time1) / (v*sqrt(time1)))) }
    
    # Reverse Extreme Spread Option:  
    if (kappa == -1) {  
    ExtremeSpread = -eta * (S * exp((b - r) * Time) * (1 + v^2 / (2 * b)) * 
        CND(eta * (m - mu * Time) / (v*sqrt(Time))) + exp(-r * Time) * Mo * 
        CND(eta * (-m + mu1 * Time) / (v*sqrt(Time))) - exp(-r * Time) * 
            Mo * v^2 / (2 * b) * exp(2 * mu1 * m / v^2) * 
        CND(eta * (m + mu1 * Time) / (v*sqrt(Time))) - S * 
            exp((b - r) * Time) * (1 + v^2 / (2 * b)) * 
        CND(eta * (-mu * (Time - time1)) / (v*sqrt(Time - time1))) - 
            exp(-r * (Time - time1)) * S * exp((b - r) * Time) * 
            (1 - v^2 / (2 * b)) * 
        CND(eta * (mu1 * (Time - time1)) / (v*sqrt(Time - time1)))) }      
    
    # Return Value:
    option = list(
        price = ExtremeSpread, 
        call = match.call() )
    class(option) = "option"
    option              
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                     DESCRIPTION:
# Barrier Options:
#   StandardBarrierOption         Standard Barrier Option
#   DoubleBarrierOption           Double Barrier Option
#   PTSingleAssetBarrierOption    Partial Time Barrier Option
#   TwoAssetBarrierOption         Two Asset Barrier
#   PTTwoAssetBarrierOption       Partial Time TwoAsset Barrier Option
#   LookBarrierOption             Look Barrier Option
#   DiscreteBarrierOption         Discrete Adjusted Barrier Option
#   SoftBarrierOption             Soft Barrier
################################################################################


StandardBarrierOption = 
function(TypeFlag = c("cdi", "cui", "pdi", "pui", "cdo", "cuo", "pdo", "puo"), 
S, X, H, K, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Standard Barrier Options

    # References:
    #   Haug, Chapter 2.10.1

    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    StandardBarrier = NA
    mu = (b - sigma ^ 2 / 2) / sigma ^ 2
    lambda = sqrt (mu ^ 2 + 2 * r / sigma ^ 2)
    X1 = log (S / X) / (sigma * sqrt(Time)) + (1 + mu) * sigma * sqrt(Time)
    X2 = log (S / H) / (sigma * sqrt(Time)) + (1 + mu) * sigma * sqrt(Time)
    y1 = log (H ^ 2 / (S * X)) / (sigma * sqrt(Time)) + (1 + mu) * sigma * 
        sqrt(Time)
    y2 = log (H / S) / (sigma * sqrt(Time)) + (1 + mu) * sigma * sqrt(Time)
    Z  = log (H / S) / (sigma * sqrt(Time)) + lambda * sigma * sqrt(Time)
    if (TypeFlag == "cdi" || TypeFlag == "cdo") { eta = +1; phi = +1 }
    if (TypeFlag == "cui" || TypeFlag == "cuo") { eta = -1; phi = +1 }
    if (TypeFlag == "pdi" || TypeFlag == "pdo") { eta = +1; phi = -1 }
    if (TypeFlag == "pui" || TypeFlag == "puo") { eta = -1; phi = -1 }
    f1 = phi * S * exp ((b - r) * Time) * CND(phi * X1) -
        phi * X * exp(-r * Time) * CND(phi * X1 - phi * sigma * sqrt(Time))
    f2 = phi * S * exp ((b - r) * Time) * CND(phi * X2) -
        phi * X * exp(-r * Time) * CND(phi * X2 - phi * sigma * sqrt(Time))
    f3 = phi * S * exp ((b - r) * Time) * (H / S) ^ (2 * (mu + 1))  *
        CND(eta * y1) - phi * X * exp(-r * Time) * (H / S) ^ (2 * mu) *
        CND(eta * y1 - eta * sigma * sqrt(Time))
    f4 = phi * S * exp ((b - r) * Time) * (H / S) ^ (2 * (mu + 1)) *
        CND(eta * y2) - phi * X * exp(-r * Time) * (H / S) ^ (2 * mu) *
        CND(eta * y2 - eta * sigma * sqrt(Time))
    f5 = K * exp (-r * Time) * (CND(eta * X2 - eta * sigma * sqrt(Time)) -
        (H / S) ^ (2 * mu) * CND(eta * y2 - eta * sigma * sqrt(Time)))
    f6 = K * ((H / S) ^ (mu + lambda) * CND(eta * Z) + (H / S)^(mu - lambda) * 
        CND(eta * Z - 2 * eta * lambda * sigma * sqrt(Time)))
    if (X >= H) {
        if (TypeFlag == "cdi") StandardBarrier = f3 + f5
        if (TypeFlag == "cui") StandardBarrier = f1 + f5
        if (TypeFlag == "pdi") StandardBarrier = f2 - f3 + f4 + f5
        if (TypeFlag == "pui") StandardBarrier = f1 - f2 + f4 + f5
        if (TypeFlag == "cdo") StandardBarrier = f1 - f3 + f6
        if (TypeFlag == "cuo") StandardBarrier = f6
        if (TypeFlag == "pdo") StandardBarrier = f1 - f2 + f3 - f4 + f6
        if (TypeFlag == "puo") StandardBarrier = f2 - f4 + f6 }
    if (X < H) {
        if (TypeFlag == "cdi") StandardBarrier = f1 - f2 + f4 + f5
        if (TypeFlag == "cui") StandardBarrier = f2 - f3 + f4 + f5
        if (TypeFlag == "pdi") StandardBarrier = f1 + f5
        if (TypeFlag == "pui") StandardBarrier = f3 + f5
        if (TypeFlag == "cdo") StandardBarrier = f2 + f6 - f4
        if (TypeFlag == "cuo") StandardBarrier = f1 - f2 + f3 - f4 + f6
        if (TypeFlag == "pdo") StandardBarrier = f6
        if (TypeFlag == "puo") StandardBarrier = f1 - f3 + f6 }
    
    # Return Value:
    option = list(
        price = StandardBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


DoubleBarrierOption = 
function(TypeFlag = c("co", "ci", "po", "pi"), S, X, L, U, Time, r, b, 
sigma, delta1, delta2)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Double barrier options
 
    # References:
    #   Haug, Chapter 2.10.2

    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    DoubleBarrier = NA
    FU = U * exp (delta1 * Time)
    E = L * exp (delta1 * Time)
    Sum1 = Sum2 = 0
    
    # Call:
    if (TypeFlag == "co" || TypeFlag == "ci") {
        for (n in -5:5) {
            d1 = (log(S * U ^ (2 * n) / (X * L ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d2 = (log(S * U ^ (2 * n) / (FU * L ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d3 = (log(L ^ (2 * n + 2) / (X * S * U ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d4 = (log(L ^ (2 * n + 2) / (FU * S * U ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            mu1 = 2 * (b - delta2 - n * (delta1 - delta2)) / sigma^2 + 1
            mu2 = 2 * n * (delta1 - delta2) / sigma^2
            mu3 = 2 * (b - delta2 + n * (delta1 - delta2)) / sigma^2 + 1
            Sum1 = Sum1 + (U^n / L ^ n) ^ mu1 * (L / S) ^ mu2 *
                (CND(d1) - CND(d2)) - (L^(n + 1) / (U ^ n * S)) ^ mu3 *
                (CND(d3) - CND(d4))
            Sum2 = Sum2 + (U^n / L ^ n) ^ (mu1 - 2) * (L/S)^mu2 *
                (CND(d1 - sigma * sqrt(Time)) - CND(d2 - sigma * sqrt(Time))) -
                (L^(n + 1) / (U ^ n * S))^(mu3 - 2) *
                (CND(d3 - sigma * sqrt(Time)) - CND(d4 - sigma * sqrt(Time))) }
        OutValue = S * exp ((b-r)*Time) * Sum1 - X * exp(-r*Time) * Sum2 }
    
    # Put:
    if (TypeFlag == "po" || TypeFlag == "pi") {
        for ( n in (-5:5) ) {
            d1 = (log(S * U ^ (2 * n) / (E * L ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d2 = (log(S * U ^ (2 * n) / (X * L ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d3 = (log(L ^ (2 * n + 2) / (E * S * U ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            d4 = (log(L ^ (2 * n + 2) / (X * S * U ^ (2 * n))) +
                (b + sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
            mu1 = 2 * (b - delta2 - n * (delta1 - delta2)) / sigma ^ 2 + 1
            mu2 = 2 * n * (delta1 - delta2) / sigma ^ 2
            mu3 = 2 * (b - delta2 + n * (delta1 - delta2)) / sigma ^ 2 + 1
            Sum1 = Sum1 + (U^n / L^n)^mu1 * (L / S) ^ mu2 *
                (CND(d1) - CND(d2)) -
                (L ^ (n + 1) / (U ^ n * S)) ^ mu3 *
                (CND(d3) - CND(d4))
            Sum2 = Sum2 + (U ^n / L^n)^(mu1 - 2) * (L/S)^mu2 *
                (CND(d1 - sigma * sqrt(Time)) - CND(d2 - sigma * sqrt(Time))) -
                (L^(n + 1) / (U ^ n * S))^(mu3 - 2) *
                (CND(d3 - sigma * sqrt(Time)) - CND(d4 - sigma * sqrt(Time))) }
         OutValue = X * exp (-r*Time) * Sum2 - S * exp((b - r)*Time) * Sum1 }
    
    # Final Values:
    if (TypeFlag == "co" || TypeFlag == "po") 
        DoubleBarrier = OutValue
    if (TypeFlag == "ci") 
        DoubleBarrier = 
            GBlackScholes("c", S, X, Time, r, b, sigma)$price - OutValue
    if (TypeFlag == "pi") 
        DoubleBarrier = 
            GBlackScholes("p", S, X, Time, r, b, sigma)$price - OutValue
    
    # Return Value:
    option = list(
        price = DoubleBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}



# ------------------------------------------------------------------------------


PTSingleAssetBarrierOption = 
function(TypeFlag = c("cdoA", "cuoA", "pdoA", "puoA", "coB1", "poB1", 
"cdoB2", "cuoB2"), S, X, H, time1, Time2, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Partial-time single asset barrier options
     
    # References:
    #   Haug, Chapter 2.10.3

    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    PartialTimeBarrier = NA
    t1 = time1
    T2 = Time2
    if (TypeFlag == "cdoA") eta = 1
    if (TypeFlag == "cuoA") eta = -1
    
    # Continue:
    d1 = (log(S/X) + (b + sigma^2/2) * T2) / (sigma * sqrt(T2))
    d2 = d1 - sigma * sqrt (T2)
    f1 = (log(S/X) + 2 * log(H/S) + (b + sigma^2/2) * T2) / (sigma * sqrt(T2))
    f2 = f1 - sigma * sqrt (T2)
    e1 = (log(S / H) + (b + sigma ^ 2 / 2) * t1) / (sigma * sqrt(t1))
    e2 = e1 - sigma * sqrt (t1)
    e3 = e1 + 2 * log (H / S) / (sigma * sqrt(t1))
    e4 = e3 - sigma * sqrt (t1)
    mu = (b - sigma ^ 2 / 2) / sigma ^ 2
    rho = sqrt (t1 / T2)
    g1 = (log(S / H) + (b + sigma ^ 2 / 2) * T2) / (sigma * sqrt(T2))
    g2 = g1 - sigma * sqrt (T2)
    g3 = g1 + 2 * log (H / S) / (sigma * sqrt(T2))
    g4 = g3 - sigma * sqrt (T2)
    z1 = CND (e2) - (H / S) ^ (2 * mu) * CND(e4)
    z2 = CND (-e2) - (H / S) ^ (2 * mu) * CND(-e4)
    z3 = CBND (g2, e2, rho) - (H / S) ^ (2 * mu) * CBND(g4, -e4, -rho)
    z4 = CBND (-g2, -e2, rho) - (H / S) ^ (2 * mu) * CBND(-g4, e4, -rho)
    z5 = CND (e1) - (H / S) ^ (2 * (mu + 1)) * CND(e3)
    z6 = CND (-e1) - (H / S) ^ (2 * (mu + 1)) * CND(-e3)
    z7 = CBND (g1, e1, rho) - (H / S) ^ (2 * (mu + 1)) * CBND(g3, -e3, -rho)
    z8 = CBND (-g1, -e1, rho) - (H / S) ^ (2 * (mu + 1)) * CBND(-g3, e3, -rho)
    
    if (TypeFlag == "cdoA" || TypeFlag == "cuoA") { 
        # call down-and out and up-and-out type A
        PartialTimeBarrier = 
            S * exp ((b - r) * T2) * (CBND(d1, eta * e1, eta * rho) -
            (H / S) ^ (2 * (mu + 1)) * CBND(f1, eta * e3, eta * rho)) -
            X * exp (-r * T2) * (CBND(d2, eta * e2, eta * rho) -
            (H / S) ^ (2 * mu) * CBND(f2, eta * e4, eta * rho)) }
    
    if (TypeFlag == "cdoB2" && X < H) {  
        # call down-and-out type B2
        PartialTimeBarrier = 
            S * exp ((b - r) * T2) * (CBND(g1, e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(g3, -e3, -rho)) 
            - X * exp (-r * T2) * (CBND(g2, e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(g4, -e4, -rho)) }
    
    if (TypeFlag == "cdoB2" && X > H) {
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("coB1", 
                S, X, H, t1, T2, r, b, sigma)$price }
    
    if (TypeFlag == "cuoB2" && X < H) {  
        # call up-and-out type B2
        PartialTimeBarrier = 
            S * exp ((b - r) * T2) * (CBND(-g1, -e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(-g3, e3, -rho)) 
            - X * exp (-r * T2) * (CBND(-g2, -e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(-g4, e4, -rho)) 
            - S * exp ((b - r) * T2) * (CBND(-d1, -e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(e3, -f1, -rho)) 
            + X * exp (-r * T2) * (CBND(-d2, -e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(e4, -f2, -rho))}
    
    if (TypeFlag == "coB1" && X > H) {  
        # call out type B1
        PartialTimeBarrier = 
            S * exp ((b - r) * T2) * (CBND(d1, e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(f1, -e3, -rho)) 
            - X * exp (-r * T2) * (CBND(d2, e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(f2, -e4, -rho)) }
    
    if (TypeFlag == "coB1" && X < H) {
        PartialTimeBarrier = 
            S * exp ((b - r) * T2) * (CBND(-g1, -e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(-g3, e3, -rho)) 
            - X * exp (-r * T2) * (CBND(-g2, -e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(-g4, e4, -rho)) 
            - S * exp ((b - r) * T2) * (CBND(-d1, -e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(-f1, e3, -rho)) 
            + X * exp (-r * T2) * (CBND(-d2, -e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(-f2, e4, -rho)) 
            + S * exp ((b - r) * T2) * (CBND(g1, e1, rho) 
            - (H / S) ^ (2 * (mu + 1)) * CBND(g3, -e3, -rho)) 
            - X * exp (-r * T2) * (CBND(g2, e2, rho) 
            - (H / S) ^ (2 * mu) * CBND(g4, -e4, -rho)) }
    
    if (TypeFlag == "pdoA") {  
        # put down-and out and up-and-out type A
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("cdoA", 
                S, X, H, t1, T2, r, b, sigma)$price - 
                S * exp((b - r) * T2) * z5 + X * exp(-r * T2) * z1 }
                
    if (TypeFlag == "puoA") {
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("cuoA", 
                S, X, H, t1, T2, r, b, sigma)$price -
                S * exp((b - r) * T2) * z6 + X * exp(-r * T2) * z2 }
                
    if (TypeFlag == "poB1") {  
        # put out type B1
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("coB1", 
                S, X, H, t1, T2, r, b, sigma)$price -
                S * exp((b - r) * T2) * z8 + X * exp(-r * T2) * z4 -
                S * exp((b - r) * T2) * z7 + X * exp(-r * T2) * z3 }
                
    if (TypeFlag == "pdoB2") {  
        # put down-and-out type B2
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("cdoB2", 
                S, X, H, t1, T2, r, b, sigma)$price - 
                S * exp((b - r) * T2) * z7 + X * exp(-r * T2) * z3 }
                
    if (TypeFlag == "puoB2") {  
        # put up-and-out type B2
        PartialTimeBarrier = 
            PTSingleAssetBarrierOption("cuoB2", 
                S, X, H, t1, T2, r, b, sigma)$price - 
                S * exp((b - r) * T2) * z8 + X * exp(-r * T2) * z4 }
    
    # Return Value:
    option = list(
        price = PartialTimeBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


TwoAssetBarrierOption = 
function(TypeFlag = c("cuo", "cui", "cdo", "cdi", "puo", "pui", "pdo", "pdi"), 
S1, S2, X, H, Time, r, b1, b2, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Partial-time singel asset barrier options
 
    # References:
    #   Haug, Chapter 2.10.4

    # FUNCTION:

    # Compute:
    TypeFlag = TypeFlag[1]
    v1 = sigma1
    v2 = sigma2
    mu1 = b1 - v1 ^ 2 / 2
    mu2 = b2 - v2 ^ 2 / 2
    d1 = (log(S1 / X) + (mu1 + v1 ^ 2 / 2) * Time) / (v1 * sqrt(Time))
    d2 = d1 - v1 * sqrt (Time)
    d3 = d1 + 2 * rho * log (H / S2) / (v2 * sqrt(Time))
    d4 = d2 + 2 * rho * log (H / S2) / (v2 * sqrt(Time))
    e1 = (log(H / S2) - (mu2 + rho * v1 * v2) * Time) / (v2 * sqrt(Time))
    e2 = e1 + rho * v1 * sqrt (Time)
    e3 = e1 - 2 * log (H / S2) / (v2 * sqrt(Time))
    e4 = e2 - 2 * log (H / S2) / (v2 * sqrt(Time))
    
    # Make Decisions:
    if (TypeFlag == "cuo" || TypeFlag == "cui") {
        eta = 1
        phi = 1 }       
    if (TypeFlag == "cdo" || TypeFlag == "cdi") {
        eta = 1
        phi = -1 }      
    if (TypeFlag == "puo" || TypeFlag == "pui") {
        eta = -1
        phi = 1 }
    if (TypeFlag == "pdo" || TypeFlag == "pdi") {
        eta = -1
        phi = -1 }

    # Calculate Knock Out Value:        
    KnockOutValue = 
        eta * S1 * exp ((b1 - r) * Time) * 
            (CBND ( eta * d1, phi * e1, -eta * phi * rho) - 
            exp (2 * (mu2 + rho * v1 * v2) * 
            log(H / S2) / v2 ^ 2) *
            CBND(eta * d3, phi * e3, -eta * phi * rho)) - 
            eta * exp(-r * Time) * X *
            (CBND(eta * d2, phi * e2, -eta * phi * rho) - 
            exp (2 * mu2 * log(H / S2) / v2 ^ 2) *
            CBND(eta * d4, phi * e4, -eta * phi * rho))
    
    # Calculate Two Asset Barrier:
    if (TypeFlag == "cuo" || TypeFlag == "cdo" || 
        TypeFlag == "puo" || TypeFlag == "pdo") 
        TwoAssetBarrier = 
            KnockOutValue           
    if (TypeFlag == "cui" || TypeFlag == "cdi")
        TwoAssetBarrier = 
            GBlackScholes("c", S1, X, Time, r, b1, v1)$price - KnockOutValue        
    if (TypeFlag == "pui" || TypeFlag == "pdi")
        TwoAssetBarrier = 
            GBlackScholes("p", S1, X, Time, r, b1, v1)$price - KnockOutValue
    
    # Return Value:
    option = list(
        price = TwoAssetBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


"PTTwoAssetBarrierOption" = 
function(TypeFlag = c("cdo", "pdo", "cdi", "pdi", "cuo", "puo", "cui", "pui"), 
S1, S2, X, H, time1, Time2, r, b1, b2, sigma1, sigma2, rho)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Partial-time two asset barrier options
     
    # References:
    #   Haug, Chapter 2.10.5

    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    t1 = time1
    T2 = Time2
    v1 = sigma1
    v2 = sigma2
    if (TypeFlag == "cdo" || TypeFlag == "pdo" || 
        TypeFlag == "cdi" || TypeFlag == "pdi") {
        phi = -1 }
    else {
        phi = 1 }
    if (TypeFlag == "cdo" || TypeFlag == "cuo" || 
        TypeFlag == "cdi" || TypeFlag == "cui") {
        eta = 1 }
    else {
        eta = -1 }
    
    mu1 = b1 - v1 ^ 2 / 2
    mu2 = b2 - v2 ^ 2 / 2
    d1 = (log(S1 / X) + (mu1 + v1 ^ 2) * T2) / (v1 * sqrt(T2))
    d2 = d1 - v1 * sqrt (T2)
    d3 = d1 + 2 * rho * log (H / S2) / (v2 * sqrt(T2))
    d4 = d2 + 2 * rho * log (H / S2) / (v2 * sqrt(T2))
    e1 = (log(H / S2) - (mu2 + rho * v1 * v2) * t1) / (v2 * sqrt(t1))
    e2 = e1 + rho * v1 * sqrt (t1)
    e3 = e1 - 2 * log (H / S2) / (v2 * sqrt(t1))
    e4 = e2 - 2 * log (H / S2) / (v2 * sqrt(t1))
    
    OutBarrierValue = 
        eta * S1 * exp ((b1 - r) * T2) *
        (CBND(eta * d1, phi * e1, -eta * phi * rho * sqrt(t1 / T2)) -
        exp(2 * log(H / S2) * (mu2 + rho * v1 * v2) / (v2 ^ 2)) *
        CBND (eta * d3, phi * e3, -eta * phi * rho * sqrt(t1 / T2))) -
        eta * exp (-r * T2) * X *
        (CBND(eta * d2, phi * e2, -eta * phi * rho * sqrt(t1 / T2)) -
        exp(2 * log(H / S2) * mu2 / (v2 ^ 2)) *
        CBND (eta * d4, phi * e4, -eta * phi * rho * sqrt(t1 / T2)))
    
    if (TypeFlag == "cdo" || TypeFlag == "cuo" || 
        TypeFlag == "pdo" || TypeFlag == "puo") 
        PartialTimeTwoAssetBarrier = 
            OutBarrierValue
    if (TypeFlag == "cui" || TypeFlag == "cdi") 
        PartialTimeTwoAssetBarrier = 
            GBlackScholes("c", S1, X, T2, r, b1, v1)$price - OutBarrierValue
    if (TypeFlag == "pui" || TypeFlag == "pdi") 
        PartialTimeTwoAssetBarrier = 
            GBlackScholes("p", S1, X, T2, r, b1, v1)$price - OutBarrierValue
    
    # Return Value:
    option = list(
        price = PartialTimeTwoAssetBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


LookBarrierOption = 
function(TypeFlag = c("cuo", "cui", "pdo", "pdi"), S, X, H, time1, Time2, 
r, b, sigma)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Look-Barrier Options

    # References:
    #   Haug, Chapter 2.10.6

    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    t1 = time1
    T2 = Time2
    # Take care of the limit t1 -> T2
    if (T2 == t1) t1 = t1*(1-1.0e-12)
    v = sigma
    hh = log(H/S)
    K = log(X/S)
    mu1 = b - sigma^2/2
    mu2 = b + sigma^2/2
    rho = sqrt(t1/T2)
    
    # Make Decisions - Settings:
    if (TypeFlag == "cuo" || TypeFlag == "cui") {
        eta = +1
        m = min (hh, K) }
    if (TypeFlag == "pdo" || TypeFlag == "pdi") {
        eta = -1
        m = max (hh, K) }
    
    # Compute the g Values:
    g1 = (CND(eta * (hh - mu2 * t1) / (sigma * sqrt(t1))) -
        exp(2 * mu2 * hh / sigma ^ 2) *
        CND(eta * (-hh - mu2 * t1) / (sigma * sqrt(t1)))) -
        (CND(eta * (m - mu2 * t1) / (sigma * sqrt(t1))) -
        exp(2 * mu2 * hh / sigma ^ 2) *
        CND(eta * (m - 2 * hh - mu2 * t1) / (sigma * sqrt(t1))))
    g2 = (CND(eta * (hh - mu1 * t1) / (sigma * sqrt(t1))) -
        exp(2 * mu1 * hh / sigma ^ 2) *
        CND(eta * (-hh - mu1 * t1) / (sigma * sqrt(t1)))) -
        (CND(eta * (m - mu1 * t1) / (sigma * sqrt(t1))) -
        exp(2 * mu1 * hh / sigma ^ 2) *
        CND(eta * (m - 2 * hh - mu1 * t1) / (sigma * sqrt(t1))))
    
    # Needed by Out Value:
    part1 = S * exp((b-r)*T2) * (1+v^2/(2*b)) * (
            CBND(
                eta*(+m-mu2*t1)/(v*sqrt(t1)), 
                eta*(-K+mu2*T2)/(v*sqrt(T2)), 
                -rho) - exp(2*mu2*hh/v^2) * 
            CBND(
                eta*(m-2*hh-mu2*t1)/(v*sqrt(t1)), 
                eta*(2*hh-K+mu2*T2)/(v*sqrt(T2)), 
                -rho) )
    part2 = - X * exp(-r*T2) * (
            CBND(
                eta*(+m-mu1*t1)/(v*sqrt(t1)), 
                eta*(-K+mu1*T2)/(v*sqrt(T2)), 
                -rho) - exp(2*mu1*hh/v^2) * 
            CBND(
                eta*(m-2*hh-mu1*t1)/(v*sqrt(t1)), 
                eta*(2*hh-K+mu1*T2)/(v*sqrt(T2)), 
                -rho) )
    part3 = -exp(-r*T2) * v^2/(2*b) * (
            S*(S/X)^(-2*b/v^2) * 
                CBND(
                    eta * (m + mu1 * t1) / (v * sqrt(t1)), 
                    eta * (-K - mu1 * T2) / (v * sqrt(T2)), 
                    -rho) - 
            H*(H/X)^(-2*b/v^2) * 
                CBND(
                    eta*(m - 2 * hh + mu1 * t1) / (v * sqrt(t1)), 
                    eta * (2 * hh - K - mu1 * T2) / (v * sqrt(T2)), 
                    -rho) )
    part4 = S * exp((b-r)*T2) * ((1+v^2/(2 * b)) * 
        CND(eta*mu2*(T2-t1)/(v*sqrt(T2-t1))) + 
        exp(-b*(T2-t1))*(1-v^2/(2*b)) * 
            CND(eta*(-mu1*(T2-t1))/(v*sqrt(T2-t1))))*g1 - 
        exp(-r*T2)*X*g2
    
    # Calculate Out Value:
    OutValue = eta * (part1 + part2 + part3 + part4)
    
    # Option Price:
    if (TypeFlag == "cuo" || TypeFlag == "pdo") 
        LookBarrier = 
            OutValue
    if (TypeFlag == "cui") 
        LookBarrier = PTFixedStrikeLookbackOption("c", S, X, t1, T2, 
            r, b, sigma) - OutValue
    if (TypeFlag == "pdi") 
        LookBarrier = PTFixedStrikeLookbackOption("p", S, X, t1, T2, 
            r, b, sigma) - OutValue
    
    # Return Value:
    option = list(
        price = LookBarrier, 
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


DiscreteBarrierOption = 
function(S, H, sigma, dt) 
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Discrete Barrier Options
     
    # References:
    #   Haug, Chapter 2.10.

    # FUNCTION:
    
    # Compute:
    DiscreteBarrier = NA
    if (H > S) { 
        DiscreteBarrier = H * exp(0.5826 * sigma * sqrt(dt)) }
    if (H < S) {
        DiscreteBarrier = H * exp(-0.5826 * sigma * sqrt(dt)) }
    
    # Return Value:
    option = list(
        price = DiscreteBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


SoftBarrierOption = 
function(TypeFlag = c("cdi", "cdo", "pdi", "pdo"), S, X, L, U, Time , 
r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Soft Barrier Option

    # References:
    #   Haug, Haug Chapter 2.10.8


    # FUNCTION:
    
    # Compute:
    TypeFlag = TypeFlag[1]
    v = sigma
    
    # Make Decisions - Settings:
    if (TypeFlag == "cdi" || TypeFlag == "cdo") {
        eta = 1}
    else {
        eta = -1 }
        
    # Continue:  
    mu = (b + v ^ 2 / 2) / v ^ 2
    lambda1 = exp(-1 / 2 * v ^ 2 * Time * (mu + 0.5) * (mu - 0.5))
    lambda2 = exp(-1 / 2 * v ^ 2 * Time * (mu - 0.5) * (mu - 1.5))
    d1 = log(U ^ 2 / (S * X)) / (v * sqrt(Time)) + mu * v * sqrt(Time)
    d2 = d1 - (mu + 0.5) * v * sqrt(Time)
    d3 = log(U ^ 2 / (S * X)) / (v * sqrt(Time)) + (mu - 1) * v * sqrt(Time)
    d4 = d3 - (mu - 0.5) * v * sqrt(Time)
    e1 = log(L ^ 2 / (S * X)) / (v * sqrt(Time)) + mu * v * sqrt(Time)
    e2 = e1 - (mu + 0.5) * v * sqrt(Time)
    e3 = log(L ^ 2 / (S * X)) / (v * sqrt(Time)) + (mu - 1) * v * sqrt(Time)
    e4 = e3 - (mu - 0.5) * v * sqrt(Time)  
    
    # Compute Value:
    Value = eta * 1 / (U - L) * (S * exp((b - r) * Time) * S ^ (-2 * mu) * 
        (S * X) ^ (mu + 0.5) / (2 * (mu + 0.5)) * 
        ((U ^ 2 / (S * X)) ^ (mu + 0.5) * CND(eta * d1) - 
        lambda1 * CND(eta * d2) - (L ^ 2 / (S * X)) ^ (mu + 0.5) * 
        CND(eta * e1) + lambda1 * CND(eta * e2)) - X * exp(-r * Time) * 
            S ^ (-2 * (mu - 1)) * (S * X) ^ (mu - 0.5) / (2 * (mu - 0.5)) * 
            ((U ^ 2 / (S * X)) ^ (mu - 0.5) * CND(eta * d3) - lambda2 * 
        CND(eta * d4) - (L ^ 2 / (S * X)) ^ (mu - 0.5) * 
        CND(eta * e3) + lambda2 * 
        CND(eta * e4)))
    ### print(Value)
    
    # Continue: 
    if (TypeFlag == "cdi" || TypeFlag == "pui") {
        SoftBarrier = 
            Value }
    if (TypeFlag == "cdo") {
        SoftBarrier = 
            GBSOption("c", S, X, Time, r, b, v)$price - Value }
    if (TypeFlag == "puo") {
        SoftBarrier = 
            GBSOption("p", S, X, Tome, r, b, v)$price - Value }
    
    # Return Value:
    option = list(
        price = SoftBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                     DESCRIPTION:
# Binary Options:
#  GapOption                     Gap Option
#  CashOrNothingOption           Cash Or Nothing Option
#  TwoAssetCashOrNothingOption   Two Asset Cash-Or Nothing Option
#  AssetOrNothingOption          Asset Or Nothing Option
#  SuperShareOption              Super Share Option
#  BinaryBarrierOption           Binary Barrier Option
################################################################################


GapOption = 
function(TypeFlag = c("c", "p"), S, X1, X2, Time, r, b, sigma)
{   # A function imlemented by Diethelm Wuertz           
    
    # Description:
    #   Gap Options
 
    # References:
    #   Haug, Haug Chapter 2.11.1

    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    d1 = (log(S/X1) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    d2 = d1 - sigma*sqrt (Time)
    if (TypeFlag == "c") 
        GapOption = 
            S*exp((b-r)*Time)*CND(d1) - X2*exp(-r*Time)*CND(d2)
    if (TypeFlag == "p") 
        GapOption =  
            X2*exp(-r*Time)*CND(-d2) - S*exp((b-r)*Time)*CND(-d1) 
    
    # Return Value:
    option = list(
        price = GapOption, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


CashOrNothingOption = 
function(TypeFlag = c("c", "p"), S, X, K, Time, r, b, sigma)
{   # A function imlemented by Diethelm Wuertz           

    # Description:
    #   Cash-Or-Nothing Options

    # References:
    #   Haug, Cahpter 2.11.2

    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    d = (log(S / X) + (b - sigma ^ 2 / 2) * Time) / (sigma * sqrt(Time))
    if (TypeFlag == "c") 
        CashOrNothing = K * exp (-r * Time) * CND(d)
    if (TypeFlag == "p") 
        CashOrNothing = K * exp (-r * Time) * CND(-d)
    
    # Return Value:
    option = list(
        price = CashOrNothing, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


TwoAssetCashOrNothingOption = 
function(TypeFlag = c("c", "p", "ud", "du"), S1, S2, X1, X2, K, Time, r, 
b1, b2, sigma1, sigma2, rho)
{   # A function imlemented by Diethelm Wuertz           

    # Description:
    #   Two Asset Cash-Or-Nothing Options

    # References:
    #   Haug, Chapter 2.11.3

    # Arguments:
    #   1: Asset One, 2: Asset Two
    #   TypeFlag
    #       1   Call
    #       2   Put
    #       3   Up-Down
    #       4   Down-Up
    #   S=c(S1,S2)      Asset Prices
    #   K           Payout
    #   X=c(X1,X2)      Strikes
    #   b=c(b1,b2)      Cost-of-Carry
    #   sigma=c(sigma1,sigma2)  Volatilities
    #   rho         Correlation
    #
    
    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    d11 = (log(S1/X1) + (b1 - sigma1^2/2) * Time) / 
        (sigma1*sqrt(Time))
    d22 = (log(S2/X2) + (b2 - sigma2^2/2) * Time) / 
        (sigma2*sqrt(Time))
    # Select:
    if (TypeFlag == "c") 
        TwoAssetCashOrNothing = K * exp (-r * Time) * 
            CBND( d11,  d22,  rho)
    if (TypeFlag == "p") 
        TwoAssetCashOrNothing = K * exp (-r * Time) * 
            CBND(-d11, -d22,  rho)
    if (TypeFlag == "ud")     
        TwoAssetCashOrNothing = K * exp (-r * Time) * 
            CBND( d11, -d22, -rho)
    if (TypeFlag == "du") 
        TwoAssetCashOrNothing = K * exp (-r * Time) * 
            CBND(-d11,  d22, -rho)
    
    # Return Value:
    option = list(
        price = TwoAssetCashOrNothing, 
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


AssetOrNothingOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function imlemented by Diethelm Wuertz           

    # Description:
    #   Asset-or-Nothing Options

    # Reference:
    #   Cox Rubinstein (1985)
    #   Haug, Chapter 2.11.4
    
    
    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    d = (log(S/X) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    if (TypeFlag == "c") 
            AssetOrNothing = S * exp ((b - r) * Time) * CND( d)
    if (TypeFlag == "p") 
            AssetOrNothing = S * exp ((b - r) * Time) * CND(-d)
    
    # Return Value:
    option = list(
        price = AssetOrNothing, 
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


SuperShareOption = 
function(S, XL, XH, Time, r, b, sigma)
{   # A function imlemented by Diethelm Wuertz           

    # Description:
    #   Supershare Options

    # Reference:
    #   Hakansson (1976)
    #   Haug, Chapter 2.11.5

    # FUNCTION:
    
    # Compute Price:
    d1 = (log(S/XL) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    d2 = (log(S/XH) + (b + sigma^2 / 2) * Time) / (sigma * sqrt(Time))
    SuperShare = (S * exp((b-r)*Time) / XL) * (CND(d1) - CND(d2))
    
    # Return Value:
    option = list(
        price = SuperShare, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


BinaryBarrierOption = 
function(TypeFlag = as.character(1:28), S, X, H, K, Time, r, b, sigma, 
eta, phi)
{   # A function imlemented by Diethelm Wuertz           
    
    # Description:
    #   Binary Barrier Options
    
    # Reference:
    #   Reiner and Rubinstein (1991)
    #   Haug, Chapter 2.11.6

    # FUNCTION:  
    
    # Compute Price:
    TypeFlag = as.integer(TypeFlag[1])
    eta = rep(c(+1,-1), 14)[TypeFlag]
    #         1  2  3  4  5  6  7  8  9 10 11 12 13 14 
    #        15 16 17 18 19 20 21 22 23 24 25 26 27 28
    phi = c(+0,+0,+0,+0,-1,+1,-1,+1,+1,-1,+1,-1,+1,+1,
             +1,+1,-1,-1,-1,-1,+1,+1,+1,+1,-1,-1,-1,-1)[TypeFlag]
    v = sigma
    mu = (b - v ^ 2 / 2) / v ^ 2
    lambda = sqrt(mu ^ 2 + 2 * r / v ^ 2)
    X1 = log(S / X) / (v * sqrt(Time)) + (mu + 1) * v * sqrt(Time)
    X2 = log(S / H) / (v * sqrt(Time)) + (mu + 1) * v * sqrt(Time)
    y1 = log(H ^ 2 / (S * X)) / (v * sqrt(Time)) + (mu + 1) * v * sqrt(Time)
    y2 = log(H / S) / (v * sqrt(Time)) + (mu + 1) * v * sqrt(Time)
    Z = log(H / S) / (v * sqrt(Time)) + lambda * v * sqrt(Time)
    
    # Values:  
    a1 = S * exp((b - r) * Time) * CND(phi * X1)
    b1 = K * exp(-r * Time) * CND(phi * X1 - phi * v * sqrt(Time))
    a2 = S * exp((b - r) * Time) * CND(phi * X2)
    b2 = K * exp(-r * Time) * CND(phi * X2 - phi * v * sqrt(Time))
    a3 = S * exp((b - r) * Time) * (H / S) ^ (2 * (mu + 1)) * 
        CND(eta * y1)
    b3 = K * exp(-r * Time) * (H / S) ^ (2 * mu) * 
        CND(eta * y1 - eta * v * sqrt(Time))
    a4 = S * exp((b - r) * Time) * (H / S) ^ (2 * (mu + 1)) * 
        CND(eta * y2)
    b4 = K * exp(-r * Time) * (H / S) ^ (2 * mu) * 
        CND(eta * y2 - eta * v * sqrt(Time))
    a5 = K * ((H / S) ^ (mu + lambda) * 
        CND(eta * Z) + (H / S) ^ (mu - lambda) * 
        CND(eta * Z - 2 * eta * lambda * v * sqrt(Time)))
    # Select:
    BinaryBarrier = NA
    if (X > H) {
        if (TypeFlag ==  1) BinaryBarrier = a5
        if (TypeFlag ==  2) BinaryBarrier = a5
        if (TypeFlag ==  3) BinaryBarrier = a5
        if (TypeFlag ==  4) BinaryBarrier = a5      
        if (TypeFlag ==  5) BinaryBarrier = b2 + b4
        if (TypeFlag ==  6) BinaryBarrier = b2 + b4     
        if (TypeFlag ==  7) BinaryBarrier = a2 + a4
        if (TypeFlag ==  8) BinaryBarrier = a2 + a4     
        if (TypeFlag ==  9) BinaryBarrier = b2 - b4
        if (TypeFlag == 10) BinaryBarrier = b2 - b4 
        if (TypeFlag == 11) BinaryBarrier = a2 - a4
        if (TypeFlag == 12) BinaryBarrier = a2 - a4     
        if (TypeFlag == 13) BinaryBarrier = b3
        if (TypeFlag == 14) BinaryBarrier = b3
        if (TypeFlag == 15) BinaryBarrier = a3
        if (TypeFlag == 16) BinaryBarrier = a1
        if (TypeFlag == 17) BinaryBarrier = b2 - b3 + b4
        if (TypeFlag == 18) BinaryBarrier = b1 - b2 + b4
        if (TypeFlag == 19) BinaryBarrier = a2 - a3 + a4
        if (TypeFlag == 20) BinaryBarrier = a1 - a2 + a3
        if (TypeFlag == 21) BinaryBarrier = b1 - b3
        if (TypeFlag == 22) BinaryBarrier = 0
        if (TypeFlag == 23) BinaryBarrier = a1 - a3
        if (TypeFlag == 24) BinaryBarrier = 0
        if (TypeFlag == 25) BinaryBarrier = b1 - b2 + b3 - b4
        if (TypeFlag == 26) BinaryBarrier = b2 - b4
        if (TypeFlag == 27) BinaryBarrier = a1 - a2 + a3 - a4
        if (TypeFlag == 28) BinaryBarrier = a2 - a4 }
    # Continue:
    if (X < H) {
        if (TypeFlag ==  1) BinaryBarrier = a5
        if (TypeFlag ==  2) BinaryBarrier = a5
        if (TypeFlag ==  3) BinaryBarrier = a5
        if (TypeFlag ==  4) BinaryBarrier = a5      
        if (TypeFlag ==  5) BinaryBarrier = b2 + b4
        if (TypeFlag ==  6) BinaryBarrier = b2 + b4     
        if (TypeFlag ==  7) BinaryBarrier = a2 + a4
        if (TypeFlag ==  8) BinaryBarrier = a2 + a4     
        if (TypeFlag ==  9) BinaryBarrier = b2 - b4
        if (TypeFlag == 10) BinaryBarrier = b2 - b4     
        if (TypeFlag == 11) BinaryBarrier = a2 - a4
        if (TypeFlag == 12) BinaryBarrier = a2 - a4     
        if (TypeFlag == 13) BinaryBarrier = b1 - b2 + b4
        if (TypeFlag == 14) BinaryBarrier = b2 - b3 + b4
        if (TypeFlag == 15) BinaryBarrier = a1 - a2 + a4
        if (TypeFlag == 16) BinaryBarrier = a2 - a3 + a4
        if (TypeFlag == 17) BinaryBarrier = b1
        if (TypeFlag == 18) BinaryBarrier = b3
        if (TypeFlag == 19) BinaryBarrier = a1
        if (TypeFlag == 20) BinaryBarrier = a3
        if (TypeFlag == 21) BinaryBarrier = b2 - b4
        if (TypeFlag == 22) BinaryBarrier = b1 - b2 + b3 - b4
        if (TypeFlag == 23) BinaryBarrier = a2 - a4
        if (TypeFlag == 24) BinaryBarrier = a1 - a2 + a3 - a4
        if (TypeFlag == 25) BinaryBarrier = 0
        if (TypeFlag == 26) BinaryBarrier = b1 - b3
        if (TypeFlag == 27) BinaryBarrier = 0
        if (TypeFlag == 28) BinaryBarrier = a1 - a3 }
    
    # Return Value:
    option = list(
        price = BinaryBarrier, 
        call = match.call() )
    class(option) = "option"
    option
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                         DESCRIPTION:
# Asian Options:
#   GeometricAverageAsianOption       Geometric Average Rate Option
# Arithmetic Average-Rater Options:
#   TurnbullWakemanAsianApproxOption  Turnbull-Wakeman Approximated Asian Option
#   LevyAsianApproxOption             Levy Approximated Asian Option
################################################################################


GeometricAverageRateOption = 
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Geometric Average Rate Options
     
    # References:
    #   Kemma and Vorst (1990)
    #   Haug, Chapter 2.12.1

    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    b.A = 0.5 * (b - sigma^2 / 6)
    sigma.A = sigma / sqrt (3)
    GeometricAverageRate = 
        GBSOption (TypeFlag = TypeFlag, S = S, X = X, Time = Time, 
            r = r, b = b.A, sigma = sigma.A)$price
    
    # Return Value:
    option = list(
        price = GeometricAverageRate,
        call = match.call() )
    class(option) = "option"
    option 
}


# ------------------------------------------------------------------------------


TurnWakeAsianApproxOption = 
function(TypeFlag = c("c", "p"), S, SA, X, Time, time, tau, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Arithmetic average rate options
    #   Turnbull-Wakeman's Approximation
 
    # References:
    #   Haug, Chapter 2.12.2

    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    m1 = (exp(b * Time) - exp(b * tau)) / (b * (Time - tau))
    m2 = 2 * exp((2 * b + sigma^2) * Time) / ((b + sigma^2) * 
        (2*b + sigma^2) * (Time - tau)^2) + 2 * exp((2 * b + sigma^2) *
        tau) / (b * (Time - tau)^2) * (1/(2 * b + sigma^2) - 
        exp(b * (Time - tau)) / (b + sigma^2))
    b.A = log(m1) / Time
    sigma.A = sqrt(log(m2) / Time - 2*b.A)
    t1 = Time - time
    if (t1 > 0) { 
        X = Time/time * X - t1/time * SA
        TurnbullWakemanAsianApprox = 
            GBSOption(TypeFlag, S, X, time, r, b.A, sigma.A)$price *
            time/Time }
    else {
        TurnbullWakemanAsianApprox = 
            GBSOption(TypeFlag, S, X, time, r, b.A, sigma.A)$price }
    
    # Return Value:
    option = list(
        price = TurnbullWakemanAsianApprox, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


LevyAsianApproxOption = 
function(TypeFlag = c("c", "p"), S, SA, X, Time, time, r, b, sigma)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Arithmetic average rate options
    #   Levy's Approximation
    
    # References:
    #   Haug, Chapter 2.12.2

    # FUNCTION:
    
    # Compute Price:
    TypeFlag = TypeFlag[1]
    SE = S / (Time*b) * (exp((b-r)*time) - exp(-r*time))
    m = 2 * S ^ 2 / (b + sigma ^ 2) * ((exp((2 * 
        b + sigma^2) * time) - 1) / (2 * b + sigma^2) - 
        (exp(b * time) - 1) / b)
    d = m / (Time^2)
    Sv = log (d) - 2 * (r * time + log(SE))
    XStar = X - (Time - time) / Time * SA
    d1 = 1 / sqrt (Sv) * (log(d) / 2 - log(XStar))
    d2 = d1 - sqrt (Sv)
    if (TypeFlag == "c") {
        LevyAsianApprox = SE * CND (d1) - XStar * exp(-r*time) * 
            CND(d2)}
    if (TypeFlag == "p") {
        LevyAsianApprox = (SE * CND(d1) - XStar * exp(-r*time) * 
            CND(d2)) - SE + XStar * exp (-r*time) }
    
    # Return Value:
    option = list(
        price = LevyAsianApprox, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


#CurranAsianApproxOption = 
#function()
#{  # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Arithmetic average rate option
    #   Curran's Approximation
    
    # References:
    #   Haug, Chapter 2.12.2

    # FUNCTION:
    
    # Compute Price:
    #   CurranAsianApprox = NA
    
    # Return Value:
#   CurranAsianApprox
#}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                     DESCRIPTION:
#  Currency Translated Options:
#   FEInDomesticFXOption          Foreign Exchang In Domestic Currency
#   QuantoOption                  Quanto Option
#   EquityLinkedFXOption          EquityLinked FX Option
#   TakeoverFXOption              Takeover FX Option
################################################################################


FEInDomesticFXOption = 
function(TypeFlag = c("c", "p"), S, E, X, Time, r, q, sigmaS, sigmaE, rho)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Foreign equity option struck in domestic currency
    
    # References:
    #   Haug, Chapter 2.13.1
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    sigma = sqrt(sigmaE^2 + sigmaS^2 + 2*rho*sigmaE*sigmaS)
    d1 = (log(E*S/X) + (r-q+sigma^2/2) * Time) / (sigma*sqrt(Time))
    d2 = d1 - sigma * sqrt(Time)
    
    # Calculate Call and Put:
    if (TypeFlag == "c") {
        ForeignEquityInDomesticFX = 
            E * S * exp(-q*Time)*CND(d1) - X * exp(-r*Time)*CND(d2) }
    if (TypeFlag == "p") {
        ForeignEquityInDomesticFX = 
            X * exp(-r*Time)*CND(-d2) - E * S * exp(-q*Time)*CND(-d1) }
    
    # Return Value:
    option = list(
        price = ForeignEquityInDomesticFX, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


QuantoOption = 
function(TypeFlag = c("c", "p"), S, Ep, X, Time, r, rf, q, sigmaS, 
sigmaE, rho)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Fixed exchange rate foreign equity options

    # References:
    #   Haug, Chapter 2.13.2

    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    d1 = (log(S/X) + (rf-q-rho*sigmaS*sigmaE + sigmaS^2/2) * Time) / 
        (sigmaS*sqrt(Time))
    d2 = d1 - sigmaS*sqrt (Time)
    
    # Calculate Call and Put:
    if (TypeFlag == "c") {
        Quanto = Ep*(S*exp((rf-r-q-rho*sigmaS*sigmaE)*Time) *
            CND(d1) - X*exp(-r*Time)*CND(d2)) }
    if (TypeFlag == "p") {
        Quanto = Ep*(X*exp(-r*Time)*CND(-d2) - 
            S*exp((rf-r-q-rho*sigmaS*sigmaE)* Time)*CND(-d1)) }
    
    # Return Value:
    option = list(
        price = Quanto, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


EquityLinkedFXOption = 
function(TypeFlag = c("c", "p"), E, S, X, Time, r, rf, q, sigmaS, 
sigmaE, rho)
{   # A function implemented by Diethelm Wuertz           
    
    # Description:
    #   Equity Linked FX Option -
    
    # References:
    #   Haug, Chapter 2.13.3
    
    # FUNCTION:
    
    # Compute Settings:
    TypeFlag = TypeFlag[1]
    vS = sigmaS
    vE = sigmaE
    d1 = (log(E / X) + (r - rf + rho * vS * vE + vE ^ 2 / 2) * Time) / 
        (vE * sqrt(Time))
    d2 = d1 - vE * sqrt(Time)
    
    # Calculate Call and Put:
    if (TypeFlag == "c") {
        EquityLinkedFXO = E * S * exp(-q * Time) * CND(d1) - 
            X * S * exp((rf - r - q - rho * vS * vE) * Time) * CND(d2) }
    if (TypeFlag == "p") {
        EquityLinkedFXO = X * S * exp((rf - r - q - rho * vS * vE) * Time) * 
            CND(-d2) - E * S * exp(-q * Time) * CND(-d1) }
    
    # Return Value:
    option = list(
        price = EquityLinkedFXO, 
        call = match.call() )
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------


TakeoverFXOption = 
function(V, B, E, X, Time, r, rf, sigmaV, sigmaE, rho)
{   # A function implemented by Diethelm Wuertz           

    # Description:
    #   Takeover FX  Option -
 
    # References:
    #   Haug, Chapter 2.13.4
    
    # FUNCTION:
    
    # Compute Settings:
    v = V
    b = B
    vV = sigmaV
    vE = sigmaE
    a1 = (log(v / b) + (rf - rho * vE * vV - vV ^ 2 / 2) * Time) / 
        (vV * sqrt(Time))
    a2 = (log(E / X) + (r - rf - vE ^ 2 / 2) * Time) / 
        (vE * sqrt(Time))
    
    # Calculate:
    TakeoverFX = b * (E * exp(-rf * Time) * 
        CBND(a2 + vE * sqrt(Time), -a1 - rho * vE * sqrt(Time), -rho) - 
            X * exp(-r * Time) * CBND(-a1, a2, -rho))           
    
    # Return Value:
    option = list(
        price = TakeoverFX, 
        call = match.call() )
    class(option) = "option"
    option
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:             DESCRIPTION:
#  hngarchSim            Simulates an HN-GARCH(1,1) Time Series Process
#  hngarchFit            Fits a HN-GARCH model by Gaussian Maximum Likelihood
#  print.hngarch         Print method, reports results
#  summary.hngarch       Summary method, diagnostic analysis
#  hngarchStats          Computes Unconditional Moments of a HN-GARCH Process
################################################################################


hngarchSim = 
function(model = list(lambda = 4, omega = 4*0.0002, alpha = 0.3*0.0002, 
beta = 0.3, gamma = 0, rf = 0), n = 1000, innov = NULL, n.start = 100, 
start.innov = NULL, rand.gen = rnorm, ...)
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Simulates a HN-GARCH time series with user supplied innovations.
    
    # Details:
    #   The function simulates a Heston Nandi Garch(1,1) process with 
    #   structure parameters specified through the list 
    #   `model(lambda, omega, alpha, beta, gamma, rf)'
    #   The function returns the simulated time series points
    #   neglecting those from the first "start.innov" innovations. 
    
    # Example:
    #   x = hngarch()
    #   plot(100*x, type="l", xlab="Day numbers", 
    #     ylab="Daily Returns %", main="Heston Nandi GARCH")
    #   S0 = 1
    #   plot(S0*exp(cumsum(x)), type="l", xlab="Day Numbers", 
    #     ylab="Daily Prices", main="Heston Nandi GARCH") }

    # FUNCTION:
    
    # Innovations:
    if (is.null(innov)) innov = rand.gen(n, ...)
    if (is.null(start.innov)) start.innov = rand.gen(n.start, ...)
    
    # Parameters:
    lambda = model$lambda
    omega = model$omega
    alpha = model$alpha
    beta = model$beta
    gamma = model$gamma
    rfr = model$rf
    
    # Start values:
    x = h = Z = c(start.innov, innov)
    nt = n.start + n
   
    # Recursion:
    h[1] = ( omega + alpha )/( 1 - alpha*gamma*gamma - beta )
    x[1] = rfr + lambda*h[1] + sqrt(h[1]) * Z[1]
    for (i in 2:nt) {
        h[i] = omega + alpha*(Z[i-1] - gamma*sqrt(h[i-1]))^2 + beta*h[i-1]
        x[i] = rfr + lambda*h[i] + sqrt(h[i]) * Z[i] } 
             
    # Time Series:
    x = x[-(1:n.start)]
        
    # Return Value:
    x
}


# ------------------------------------------------------------------------------


hngarchFit = 
function(x, model=list(lambda = -0.5, omega = var(x), alpha = 0.1*var(x), 
beta = 0.1, gamma = 0, rf = 0), symmetric = TRUE, trace = FALSE, ...)
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Fits Heston-Nandi Garch(1,1) time series model
        
    # FUNCTION:
    
    # Parameters:
    rfr = model$rf; rfr <<- rfr
    lambda = model$lambda
    omega = model$omega
    alpha = model$alpha
    beta = model$beta
    gamma = model$gamma; gamma <<- gamma
    
    # Continue:
    x <<- x
    trace <<- trace
    symmetric <<- symmetric
    opt = list()
    
    # Log-Likelihood Function:
    llh = function(par) {
        # h = sigma^2
        h = Z = x
        lambda = par[1]
        # Transform - to keep them between 0 and 1:
        omega = 1/(1+exp(-par[2]))
        alpha = 1/(1+exp(-par[3]))
        beta = 1/(1+exp(-par[4]))
        # Add gamma if selected:
        if (!symmetric) gamma = par[5]      
        # HN Garch Filter:
        h[1] = ( omega + alpha )/( 1 - alpha*gamma*gamma - beta)
        Z[1] = ( x[1] - rfr - lambda*h[1] ) / sqrt(h[1])
        for ( i in 2:length(Z) ) {
            h[i] = omega + alpha * ( Z[i-1] - gamma * sqrt(h[i-1]) )^2 +
                beta * h[i-1]
            Z[i] = ( x[i] - rfr - lambda*h[i] ) / sqrt(h[i])  }     
        # Calculate Log - Likelihood for Normal Distribution:       
        llh = -sum(log( dnorm(Z)/sqrt(h) ))
        if (trace > 0) {
            cat("Parameter Estimate\n")
            print(c(lambda, omega, alpha, beta, gamma)) }
        Z <<- Z
        h <<- h
        llh}
    
    # Transform Parameters and Calculate Start Parameters:
    par.omega = -log((1-omega)/omega)  # for 2
    par.alpha = -log((1-alpha)/alpha)  # for 3
    par.beta = -log((1-beta)/beta)     # for 4
    par.start = c(lambda, par.omega, par.alpha, par.beta)
    if(!symmetric) par.start = c(par.start, gamma)
    
    # Initial Log Likelihood:
    opt$value = llh(par = par.start)
    opt$estimate = par.start
    if (trace) {
        print(c(lambda, omega, alpha, beta, gamma))
        print(opt$value)}
     
    # Estimate parameters:
    opt = nlm(llh, par.start, ...)
    
   # Log-Likelihood:
    opt$minimum = -opt$minimum + length(x)*sqrt(2*pi)
        
    # Backtransform estimated parameters:
    lambda = opt$estimate[1]
    omega = opt$estimate[2] = (1/(1+exp(-opt$estimate[2])))
    alpha = opt$estimate[3] = (1/(1+exp(-opt$estimate[3])))
    beta = opt$estimate[4] = (1/(1+exp(-opt$estimate[4])))
    if (symmetric) opt$estimate[5] = 0
    gamma = opt$estimate[5] 
    
    # Add to Output:
    opt$model = list(lambda = lambda, omega = omega, alpha = alpha,
        beta = beta, gamma = gamma, rf = rfr)
    opt$h = h
    opt$residuals = Z
    opt$call = match.call()
    
    # Statistics - Printing:
    opt$persistence = beta + alpha*gamma*gamma
    opt$sigma2 = ( omega + alpha ) / ( 1 - opt$persistence )
    
    # Print Estimated Parameters:
    if (trace > 0) print(opt$estimate)
                
    # Return Value:
    class(opt) = "hngarch"
    opt
}


# ------------------------------------------------------------------------------


print.hngarch = 
function(x, ...)
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Print method for the  HN-GARCH time series model. 
    
    # FUNCTION:

    # Print:
    object = x
    if (!inherits(object, "hngarch")) 
        stop("method is only for garch objects")
    
    cat("\nCall:\n", deparse(object$call), "\n", sep = "")
    
    cat("\nCoefficients: lambda, omega, alpha, beta, gamma\n")
    print.default(format(object$estimate, ...), print.gap = 2, 
        quote = FALSE)
    
    cat("\nLog-Likelihood:\n")
    print.default(object$minimum)
    
    cat("\nPersistence and Variance:\n")
    print.default(c(object$persistence, object$sigma2))
    
    # Return value:
    invisible(object)
}


# ------------------------------------------------------------------------------


summary.hngarch = 
function(object, ...)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Summary method,
    #   Computes diagnostics for a HN-GARCH time series model. 

    # FUNCTION:
    
    # Print Object:
    print(object)
        
    # Create Graphs:
    plot(x, type = "l", xlab = "Days", ylab = "log-Returns", 
        main="Log-Returns", ...)
    plot(sqrt(object$h), type = "l", xlab = "Days", ylab = "sqrt(h)", 
        main = "Conditional Standard Deviations", ...)
    plot(object$residuals, type = "l", xlab = "Days", ylab = "Z", 
        main = "Innovations", ...)
    
    # Return Value:
    invisible()
}


# ******************************************************************************


hngarchStats =  
function(model)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Calculates the first 4 moments of the unconditional log
    #   return distribution for a stationary HN GARCH(1,1) process  
    #   with standard normally distributed innovations. The function   
    #   returns a list with the theoretical values for the mean, the 
    #   variance, the skewness and the kurtosis} of the (unconditional)
    #   log return distribution. We have also access to the persistence  
    #   of the corresponding HN GARCH(1,1) process and to the values 
    #   for E[sigma^2], E[sigma^4], E[sigma^6], and E[sigma^8], which are
    #   needed for the computation of the moments of the unconditional 
    #   log return distribution. The only arguments are the risk free
    #   interest rate r and the structure parameters of the HN GARCH(1,1) 
    #   process, which are specified in the model list model=list(alpha, 
    #   beta, omega, gamma, lambda)}.
    
    # Reference:
    #   A function originally written by Reto Angliker
    #   License: GPL

    # FUNCTION:
    
    # Check:
    if (model$alpha < 0) {
        warning("Negative value for the parameter alpha")}  
    if (model$beta < 0) 
        {warning("Negative value for the parameter beta") }  
    if (model$omega < 0) 
        {warning("Negative value for the parameter omega")}  
    
    # Short:
    lambda = model$lambda
    omega = model$omega
    alpha = model$alpha
    beta = model$beta
    gamma = model$gamma
    
    # Moments of the Normal Distribution        
    expect2 = 1
    expect4 = 3
    expect6 = 15
    expect8 = 105
    
    # Symmetric Case:
    if(model$gamma == 0) {
        persistence = beta  
        meansigma2 = (omega+alpha) /(1-beta)
        meansigma4 = (omega^2 + 2*omega*alpha + 2*omega*beta*meansigma2 +
            3*alpha^2 + 2*alpha*beta*meansigma2) / (1 - beta^2)
        meansigma6 = (omega^3 + 3*omega^2*alpha + 3*omega^2*beta*meansigma2 +
            9*omega*alpha^2 + 6*omega*alpha*beta*meansigma2 +
            3*omega*beta^2*meansigma4 + 15*alpha^3 +
            9*alpha^2*beta*meansigma2 + 3*alpha*beta^2*meansigma4) / (1-beta^3)
        meansigma8 = 
            (omega^4 + expect8*alpha^4 + 12*omega^2*alpha*beta*meansigma2 +
            60*alpha^3*beta*meansigma2 + 18*alpha^2*beta^2*meansigma4 +
            4*alpha*beta^3*meansigma6 + 36*omega*alpha^2*beta*meansigma2 +
            12*omega*alpha*beta^2*meansigma4 + 4*omega^3*alpha +
            4*omega^3*beta*meansigma2 + 18*omega^2*alpha^2 +
            6*omega^2*beta^2*meansigma4 + 60*omega*alpha^3 +
            4*omega*beta^3*meansigma6)/ (1 - beta^4) }  
    
    # Asymmetric Case:
    if(gamma != 0) {
        persistence = beta + alpha*gamma^2
        meansigma2 = (omega+alpha) / (1-beta-alpha*gamma^2)
        meansigma4 = (omega^2 + 2*omega*beta*meansigma2 + 
            alpha^2*expect4 + 2*beta*meansigma2*alpha*expect2 + 
            6*alpha^2*expect2*gamma^2*meansigma2 +
            2*omega*alpha*gamma^2*meansigma2 + 2*omega*alpha*expect2) / 
            (1 - beta^2 - 2*beta*alpha*gamma^2 - alpha^2*gamma^4)
        meansigma6 = 
            (3*omega*alpha^2*expect4 + 3*omega^2*alpha*gamma^2*meansigma2 +
            3*beta*meansigma2*alpha^2*expect4 + 
            3*beta^2*meansigma4*alpha*expect2 +
            15*alpha^3*expect4*gamma^2*meansigma2 +
            15*alpha^3*expect2*gamma^4*meansigma4 +
            3*omega*alpha^2*gamma^4*meansigma4 +
            3*omega^2*beta*meansigma2 + 3*omega^2*alpha*expect2 +
            3*omega*beta^2*meansigma4 + omega^3 + alpha^3*expect6 +
            18*beta*meansigma4*alpha^2*expect2*gamma^2 +
            6*omega*beta*meansigma2*alpha*expect2 +
            6*omega*beta*meansigma4*alpha*gamma^2 +
            18*omega*alpha^2*expect2*gamma^2*meansigma2) /
            (1 - 3*beta^2*alpha*gamma^2 - 3*beta*alpha^2*gamma^4 -
            alpha^3*gamma^6 - beta^3)
        meansigma8 = (omega^4 + alpha^4*expect8 +
            6*omega^2*alpha^2*expect4 + 4*omega^3*beta*meansigma2 +
            4*omega^3*alpha*expect2 + 6*omega^2*beta^2*meansigma4 +
            4*omega*beta^3*meansigma6 + 4*omega*alpha^3*expect6 +
            12*omega^2*beta*meansigma2*alpha*expect2 +
            12*omega^2*beta*meansigma4*alpha*gamma^2 +
            36*omega^2*alpha^2*expect2*gamma^2*meansigma2 +
            4*omega^3*alpha*gamma^2*meansigma2 +
            6*omega^2*alpha^2*gamma^4*meansigma4 +
            6*beta^2*meansigma4*alpha^2*expect4 +
            4*beta^3*meansigma6*alpha*expect2 +
            4*beta*meansigma2*alpha^3*expect6 +
            28*alpha^4*expect6*gamma^2*meansigma2 +
            70*alpha^4*expect4*gamma^4*meansigma4 +
            28*alpha^4*expect2*gamma^6*meansigma6 +
            4*omega*alpha^3*gamma^6*meansigma6 +
            60*beta*meansigma4*alpha^3*expect4*gamma^2 +
            60* beta*meansigma6*alpha^3*expect2*gamma^4 +
            36*beta^2*meansigma6*alpha^2*expect2*gamma^2 +
            12*omega*beta*meansigma2*alpha^2*expect4 +
            12*omega*beta^2*meansigma4*alpha*expect2 +
            12*omega*beta^2*meansigma6*alpha*gamma^2 + 
            12*omega*beta*meansigma6*alpha^2*gamma^4 +
            60*omega*alpha^3*expect4*gamma^2*meansigma2 +
            60*omega*alpha^3*expect2*gamma^4*meansigma4 +
            72*omega*beta*meansigma4*alpha^2*expect2*gamma^2) /
            (1 - beta^4 - alpha^4*gamma^8 - 4*beta^3*alpha*gamma^2 -
            6*beta^2*alpha^2*gamma^4 - 4*beta*alpha^3*gamma^6 ) }
    if (persistence > 1) { warning(paste(
        "The selected HN GARCH model is not stationary and",
        "the expressions for the moments are no more valid")) }
    
    # Leverage:
    leverage = -2*alpha*gamma*meansigma2
    
    # Unconditional Values:
    uc.mean = lambda*meansigma2
    uc.variance = lambda^2*(meansigma4 - meansigma2^2) + meansigma2
    uc.skewness = (3*lambda*meansigma4 - 3*lambda*meansigma2^2 +
        lambda^3*meansigma6 - 3*lambda^3*meansigma2*meansigma4 +
        2*lambda^3*meansigma2^3 ) / sqrt(uc.variance)^3
    uc.kurtosis = (meansigma4*3 + 6*lambda^2*meansigma6 -
        12*lambda^2*meansigma2*meansigma4 + 6*lambda^2*meansigma2^3 +
        lambda^4*meansigma8 - 4*lambda^4*meansigma2*meansigma6 +
        6*lambda^4*meansigma2^2*meansigma4 - 
        3*lambda^4*meansigma2^4 ) / uc.variance^2
    
    # Return Value:             
    list(mean = uc.mean, variance = uc.variance, skewness = uc.skewness, 
    kurtosis = uc.kurtosis, persistence = persistence, leverage = leverage,
        meansigma2 = meansigma2, meansigma4 = meansigma4, meansigma6 = 
        meansigma6, meansigma8 = meansigma8)
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:           DESCRIPTION:
#  HNGOption           Computes Option Price from the HN-GARCH Formula
#  HNGGreeks           Calculates one of the Greeks of the HN-GARCH Formula
#  HNGCharacteristics  Computes Option Price and all Greeks of HN-GARCH Model
################################################################################

    
HNGOption = 
function(TypeFlag = c("c", "p"), model, S, X, Time.inDays, r.daily) 
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Calculates the price of a HN-GARCH option.
    
    # Details:
    #   The function calculates the price of a Heston-Nandi GARCH(1,1)
    #   call or put option. 

    # FUNCTION:
    
    # Option Type:
    TypeFlag = TypeFlag[1]
    
    # Internal Function:
    f.star = function(phi, const, model, S, X, Time.inDays, r.daily) {          
        # Model Parameters:
        lambda = -1/2 
        omega = model$omega
        alpha = model$alpha 
        gamma = model$gamma + model$lambda + 1/2
        beta = model$beta   
        sigma2 = (omega + alpha)/(1 - beta - alpha * gamma^2)   
        # Function to be integrated:
        cphi0 = phi*complex(real = 0, imag = 1) 
        cphi = cphi0 + const   
        a = cphi * r.daily
        b = lambda*cphi + cphi*cphi/2           
        for (i in 2:Time.inDays) {
            a = a + cphi*r.daily + b*omega - log(1-2*alpha*b)/2
            b = cphi*(lambda+gamma) - gamma^2/2 + beta*b + 
                0.5*(cphi-gamma)^2/(1-2*alpha*b) }
        f = Re(exp(-cphi0*log(X)+cphi*log(S)+a+b*sigma2 )/cphi0)/pi 
        # Return Value:
        f }
                
    # Integrate:
    call1 = integrate(f.star, 0, Inf, const = 1, model = model, 
        S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
    call2 = integrate(f.star, 0, Inf, const = 0, model = model, 
        S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        
    # Compute Call Price:
    call.price = S/2 + exp(-r.daily*Time.inDays) * call1$value - 
        X * exp(-r.daily*Time.inDays) * ( 1/2 + call2$value )
    
    # Select Option Price:
    price = NA
    if (TypeFlag == "c" ) price = call.price
    if (TypeFlag == "p" ) price = call.price + X*exp(-r.daily*Time.inDays) - S
        
    # Return Value:
    option = list(
        price = price, 
        call = match.call())
    class(option) = "option"
    option
}


# ------------------------------------------------------------------------------
    
    
HNGGreeks = 
function(Selection = c("Delta", "Gamma"), TypeFlag = c("c", "p"), model, 
S, X, Time.inDays, r.daily) 
{   # A function implemented by Diethelm Wuertz
    
    # Description:
    #   Calculates the Greeks of a HN-GARCH option.
    
    # Details:
    #   The function calculates the delta and gamma Greeks of 
    #   a Heston Nandi GARCH(1,1) call or put option. 
    
    # FUNCTION:
    
    # Type Flags:
    Selection = Selection[1]
    TypeFlag = TypeFlag[1]
        
    # Internal Function:
    f = function(phi, const, model, S, X, Time.inDays, r.daily) {           
        # Model Parameters:
        lambda = -1/2 
        omega = model$omega
        alpha = model$alpha 
        gamma = model$gamma + model$lambda + 1/2
        beta = model$beta   
        sigma2 = (omega + alpha)/(1 - beta - alpha * gamma^2)   
        # Function to be integrated:
        cphi0 = phi*complex(real = 0, imag = 1) 
        cphi = cphi0 + const   
        a = cphi * r.daily
        b = lambda*cphi + cphi*cphi/2           
        for (i in 2:Time.inDays) {
            a = a + cphi*r.daily + b*omega - log(1-2*alpha*b)/2
            b = cphi*(lambda+gamma) - gamma^2/2 + beta*b + 
                0.5*(cphi-gamma)^2/(1-2*alpha*b) }
        f = exp(-cphi0*log(X)+cphi*log(S)+a+b*sigma2)/cphi0/pi
        # Return Value:
        f }
            
    # Delta:
    if (Selection == "Delta") {
        fdelta = function(phi, const, model, S, X, Time.inDays, r.daily) {          
            # Function to be integrated:
            cphi0 = phi * complex(real = 0, imag = 1) 
            cphi = cphi0 + const
            fdelta = cphi * 
                f(phi, const, model, S, X, Time.inDays, r.daily) / S
            # Return Value:
            Re(fdelta) }
        # Integrate:
        delta1 = integrate(fdelta, 0, Inf, const = 1, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        delta2 = integrate(fdelta, 0, Inf, const = 0, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily) 
        # Compute Call and Put Delta :
        greek = 1/2 + exp(-r.daily*Time.inDays) * delta1$value - 
            X * exp(-r.daily*Time.inDays) * delta2$value 
        if (TypeFlag == "p") greek = greek - 1 }
    
    # Gamma:
    if (Selection == "Gamma") {
        fgamma = function(phi, const, model, S, X, Time.inDays, r.daily) {          
            # Function to be integrated:
            cphi0 = phi * complex(real = 0, imag = 1) 
            cphi = cphi0 + const
            fgamma = cphi * ( cphi - 1 ) *
                f(phi, const, model, S, X, Time.inDays, r.daily) / S^2
            # Return Value:
            Re(fgamma) }
        # Integrate:    
        gamma1 = integrate(fgamma, 0, Inf, const = 1, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        gamma2 = integrate(fgamma, 0, Inf, const = 0, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        # Compute Call and Put Gamma :
        greek = put.gamma = exp(-r.daily*Time.inDays) * gamma1$value - 
            X * exp(-r.daily*Time.inDays) * gamma2$value }
        
    # Return Value: 
    greek
}


# ------------------------------------------------------------------------------


HNGCharacteristics = 
function(TypeFlag = c("c", "p"), model, S, X, Time.inDays, r.daily)
{   # A function implemented by Diethelm Wuertz
     
    # Description:
    #   The function calculates the option price for the Heston 
    #   Nandi Garch(1,1) option model together with the delta 
    #   and gamma option sensitivies.

    # FUNCTION: 
    
    # Premium and Function Call to all Greeks
    TypeFlag = TypeFlag[1]
    premium = HNGOption(TypeFlag, model, S, X, Time.inDays, r.daily)  
    delta = HNGGreeks("Delta", TypeFlag, model, S, X, Time.inDays, r.daily)  
    gamma = HNGGreeks("Gamma", TypeFlag, model, S, X, Time.inDays, r.daily)  
    
    # Return Value:
    list(premium = premium, delta = delta, gamma = gamma)   
} 


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# for the Rmetrics port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:             DESCRIPTION:                      
#  runif.pseudo           Uniform Pseudo Random number sequence
#  rnorm.pseudo           Normal Pseudo Random number sequence
#  runif.halton           Uniform Halton low discrepancy sequence
#  rnorm.halton           Normal Halton low discrepancy sequence
#  runif.sobol            Uniform Sobol low discrepancy sequence
#  rnorm.sobol            Normal Sobol low discrepancy sequence
###############################################################################


runif.pseudo = 
function(n, dimension, init = NULL) 
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Uniform Pseudo Random number sequence   
    
    # FUNCTION:
    
    # Deviates:
    result = matrix(runif(n*dimension), ncol = dimension)
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


rnorm.pseudo = 
function(n, dimension, init = TRUE) 
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Normal Pseudo Random number sequence    
    
    # FUNCTION:
    
    # Deviates:
    result = matrix(rnorm(n*dimension), ncol = dimension)
    
    # Return Value:
    result
}


# -----------------------------------------------------------------------------
    

runif.halton = 
function (n, dimension, init = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Uniform Halton Low Discrepancy Sequence

    # Details: 
    #   DIMENSION : dimension <= 200
    #       N : LD numbers to create

    # FUNCTION:
    
    # Restart Settings:
    if (init) {
        runif.halton.seed <<- list()
        runif.halton.seed$base <<- rep(0, dimension)
        runif.halton.seed$offset <<- 0 }
    
    # Generate:
    qn = rep(0, n*dimension)
    
    # SUBROUTINE HALTON(QN, N, DIMEN, BASE, OFFSET, INIT, TRANSFORM)
    result = .Fortran("halton",  
        as.double(qn),
        as.integer(n),
        as.integer(dimension),
        as.integer(runif.halton.seed$base),
        as.integer(runif.halton.seed$offset),
        as.integer(init),
        as.integer(0),
        PACKAGE = "fOptions")
    
    # For the next numbers save:    
    runif.halton.seed$base <<- result[[4]]
    runif.halton.seed$offset <<- result[[5]]
    
    # Deviates:
    result = matrix(result[[1]], ncol = dimension)
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


rnorm.halton = 
function (n, dimension, init = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Normal Halton Low Discrepancy Sequence

    # Details: 
    #   DIMENSION : dimension <= 200
    #       N : LD numbers to create

    # FUNCTION:
    
    # Restart Settings:
    if (init) {
        rnorm.halton.seed <<- list()
        rnorm.halton.seed$base <<- rep(0, dimension)
        rnorm.halton.seed$offset <<- 0 }
    
    # Generate:
    qn = rep(0, n*dimension)
    
    # SUBROUTINE HALTON(QN, N, DIMEN, BASE, OFFSET, INIT, TRANSFORM)
    result = .Fortran("halton",  
        as.double(qn),
        as.integer(n),
        as.integer(dimension),
        as.integer(rnorm.halton.seed$base),
        as.integer(rnorm.halton.seed$offset),
        as.integer(init),
        as.integer(1),
        PACKAGE = "fOptions")
        
    # For the next numbers save:    
    rnorm.halton.seed$base <<- result[[4]]
    rnorm.halton.seed$offset <<- result[[5]]
    
    # Deviates:
    result = matrix(result[[1]], ncol = dimension)
    
    # Return Value:
    result
}


# -----------------------------------------------------------------------------


runif.sobol = 
function (n, dimension, init = TRUE, scrambling = 0, seed = 4711)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Uniform Sobol Low Discrepancy Sequence

    # Details: 
    #   DIMENSION : dimension <= 200
    #           N : LD numbers to create
    #  SCRAMBLING : One of the numbers 0,1,2,3
    #

    # FUNCTION:
    
    # Restart Settings:
    if (init) {
        runif.sobol.seed <<- list()
        runif.sobol.seed$quasi <<- rep(0, dimension)
        runif.sobol.seed$ll <<- 0
        runif.sobol.seed$count <<- 0
        runif.sobol.seed$sv <<- rep(0, dimension*30)
        runif.sobol.seed$seed <<- seed }
    
    # Generate:
    qn = rep(0.0, n*dimension)
    
    # SSOBOL(QN,N,DIMEN,QUASI,LL,COUNT,SV,IFLAG,SEED,INIT,TRANSFORM)
    result = .Fortran("sobol",  
        as.double(qn),
        as.integer(n),
        as.integer(dimension),
        as.double (runif.sobol.seed$quasi),
        as.integer(runif.sobol.seed$ll),
        as.integer(runif.sobol.seed$count),
        as.integer(runif.sobol.seed$sv),
        as.integer(scrambling),
        as.integer(runif.sobol.seed$seed),
        as.integer(init),
        as.integer(0),
        PACKAGE = "fOptions")
        
    # For the next numbers save:    
    runif.sobol.seed$quasi <<- result[[4]]
    runif.sobol.seed$ll <<- result[[5]]
    runif.sobol.seed$count <<- result[[6]]
    runif.sobol.seed$sv <<- result[[7]]
    runif.sobol.seed$seed <<- result[[9]]
    
    # Deviates:
    result = matrix(result[[1]], ncol = dimension)
    
    # Return Value:
    result
}


# ------------------------------------------------------------------------------


rnorm.sobol = 
function (n, dimension, init = TRUE, scrambling = 0, seed = 4711)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Normal Sobol Low Discrepancy Sequence

    # Details: 
    #   DIMENSION : dimension <= 200
    #           N : LD numbers to create
    #  SCRAMBLING : One of the numbers 0,1,2,3

    # FUNCTION:
    
    # Restart Settings:
    if (init) {
        rnorm.sobol.seed <<- list()
        rnorm.sobol.seed$quasi <<- rep(0, dimension)
        rnorm.sobol.seed$ll <<- 0
        rnorm.sobol.seed$count <<- 0
        rnorm.sobol.seed$sv <<- rep(0, dimension*30)
        rnorm.sobol.seed$seed <<- seed  }

    # Generate:
    qn = rep(0.0, n*dimension)      
    
    # SSOBOL(QN,N,DIMEN,QUASI,LL,COUNT,SV,IFLAG,SEED,INIT,TRANSFORM)
    result = .Fortran("sobol",  
        as.double(qn),
        as.integer(n),
        as.integer(dimension),
        as.double (rnorm.sobol.seed$quasi),
        as.integer(rnorm.sobol.seed$ll),
        as.integer(rnorm.sobol.seed$count),
        as.integer(rnorm.sobol.seed$sv),
        as.integer(scrambling),
        as.integer(rnorm.sobol.seed$seed),
        as.integer(init),
        as.integer(1),
        PACKAGE = "fOptions")   
                
    # For the next numbers save:    
    rnorm.sobol.seed$quasi <<- result[[4]]
    rnorm.sobol.seed$ll <<- result[[5]]
    rnorm.sobol.seed$count <<- result[[6]]
    rnorm.sobol.seed$sv <<- result[[7]]
    rnorm.sobol.seed$seed <<- result[[9]]
    
    # Deviates:
    result = matrix(result[[1]], ncol = dimension) 
    
    # Return Value:
    result
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


################################################################################
# FUNCTION:                  DESCRIPTION:
#  MonteCarloOption           Valuate Options by Monte Carlo Simulation
################################################################################


MonteCarloOption = function(delta.t, pathLength, mcSteps, mcLoops, 
init = TRUE, innovations.gen, path.gen, payoff.calc, antithetic = TRUE, 
standardization = FALSE, trace = TRUE, ...)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Valuates Options by Monte Carlo Simulation

    # Arguments:
    #   delta.t    - The length of the time interval, by default one day.
    #   pathLength - Number of Time Intervals which add up to the path.
    #   mcSteps    - The number of Monte Carlo Steps performed in one loop.
    #   mcLoops    - The number of Monte Carlo Loops
    #   init       - Should the random number generator be initialized ?
    #                This variable must appear in the argument list of the
    #                random number generator, even it will not ne used
    #   innovations.gen 
    #              - the generator function for the innovations
    #   path.gen   - the generator for the MC paths
    #   payoff.calc
    #              - the payoff claculator function
    #   antithetic - if TRUE, antithetic paths are used in the MC simulation
    #   standardization
    #              - if TRUE, the random numbers will be standardized so that
    #                their mean is zero and their variance is zero
    #   trace      - a logical, should the iteration path be traced ?
    #   ...        - additional parameters passed to innovations.gen.
    
    # Notes:
    #   Global Variables:
    #     The options parameter must be globally available. 
    #     For a Black-Scholes or a simple Asian Option these are:
    #     TypeFlag, S, X, Time, r, b, sigma
    #   Required Functions:
    #     The user must specify the following functions:
    #     innovations.gen()
    #     path.gen()
    #     payoff.calc()
    
    # FUNCTION
        
    # Monte Carlo Simulation:
    delta.t <<- delta.t
    if (trace) cat("\nMonte Carlo Simulation Path:\n\n")
        iteration = rep(0, length = mcLoops)
        # MC Iteration Loop:
        cat("\Loop:\t", "No\t")
        for ( i in 1:mcLoops ) {
            if ( i > 1) init = FALSE
            # Generate Innovations:
              eps =  innovations.gen(mcSteps, pathLength, init = init, ...)
            # Use Antithetic Variates if requested:
              if (antithetic) 
                  eps = rbind(eps, -eps)
            # Standardize Variates if requested:
              if (standardization) eps = 
                 (epsilon-mean(epsilon))/sqrt(var(as.vector(epsilon)))
            # Calculate for each path the option price:
              path = t(path.gen(eps))
              payoff = NULL
              for (j in 1:dim(path)[1]) 
                  payoff = c(payoff, payoff.calc(path[, j])) 
              iteration[i] = mean(payoff)
            # Trace the Simualtion if desired:
              if (trace) 
                 cat("\nLoop:\t", i, "\t:", iteration[i], sum(iteration)/i ) 
        }
        if (trace) cat("\n")

    # Return Value:
    iteration
}


# ******************************************************************************


# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA 02111-1307 USA

# Copyrights (C)
# for this R-port: 
#   Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
#   date: Terry Therneau <therneau@mayo.edu>
#     R port by Th. Lumley <thomas@biostat.washington.edu>  K. Halvorsen 
#       <khal@alumni.uv.es>, and Kurt Hornik <Kurt.Hornik@R-project.org>
#   ts: Collected by Brian Ripley. See SOURCES
#   tseries: Compiled by Adrian Trapletti <a.trapletti@bluewin.ch>
# for ical:
#   libical: Libical is an Open Source implementation of the IETF's 
#     iCalendar Calendaring and Scheduling protocols. (RFC 2445, 2446, 
#     and 2447). It parses iCal components and provides a C API for 
#     manipulating the component properties, parameters, and subcomponents.
#   Olsen's VTIMEZONE: These data files are released under the GNU 
#     General Public License, in keeping with the license options of 
#     libical. 
# for the holiday database:
#   holiday information collected from the internet and governmental 
#   sources obtained from a few dozens of websites


################################################################################
# FUNCTION:             DESCRIPTION:
#  xmpfOptions           Popups the example menu
################################################################################


xmpfOptions =
function() 
{   # A function implemented by Diethelm WUertz

    # Description:
    #   Popups the example menu
    
    # FUNCTION:
    
    # Popup:    
    path = paste(.Library,"/fOptions", sep = "") 
    entries = read.00Index (file.path(path, "demo", "00Index")) 
    example = select.list(entries[,1])
    selected = 0
    for (i in 1:length(entries[,1])) {
        if (example == entries[i,1]) selected = i}
    if (example == "") {
        cat("\nNo demo selected\n")}
    else    {
        cat("\nLibrary: ", "fOptions", "\nExample: ", 
            entries[selected,1], 
        "\nTitle:   ", entries[selected,2], "\n")
        source(paste(path, "/demo/", example, ".R", sep=""), 
            echo = TRUE, verbose = FALSE)}
    if(TRUE) cat("\n")
    
    # Return Value:
    invisible()
}
    
    
#*******************************************************************************
# fOptions - A SOFTWARE COLLECTION FOR FINANCIAL ENGINEERS
# PART IV: Pricing and Hedging of Options
#
# collected by Diethelm Wuertz
#    
#*******************************************************************************
                                                        

# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C) 
# this R-port: 
#   by Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# for the code accessed (or partly included) from other R-ports:
#   R: see R's copyright and license file
# for Haug's Option Pricing Formulas:
#   Formulas are implemented along the book and the Excel spreadsheets of 
#     E.G. Haug, "The Complete Guide to Option Pricing"; documentation
#     is partly taken from www.derivicom.com which implements
#     a C Library based on Haug. For non-academic and commercial use 
#     we recommend the professional software from "www.derivicom.com".  


#*******************************************************************************

# Default Settings:

    xmpOptions = function(prompt = "") {invisible(prompt)}
 

.First.lib = 
function(lib, pkg)
{   # A function implemented by D. Wuertz
    
    # Package:
    cat("\nfOptions:   Valuation of Options ")
    
    # Requires:
    # DEBUG <- FALSE
    # sink("@sink@")           
    # ...
    # sink()
    # unlink("@sink@")

    # Load dll: 
    library.dynam("fOptions", pkg, lib)
}

