# 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: 
#   1999 - 2004, Diethelm Wuertz, GPL
#   Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
#   info@rmetrics.org
#   www.rmetrics.org
# for the code accessed (or partly included) from other R-ports:
#   see R's copyright and license files
# for the code accessed (or partly included) from contributed R-ports
# and other sources
#   see Rmetrics's copyright file


################################################################################
# FUNCTION:                  DESCRIPTION:
#  fOPTION                    S4 Class Representation
# 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
# S3 METHODS:                DESCRIPTION:
#  print.option               Print Method
#  summary.otion              Summary Method
################################################################################


setClass("fOPTION", 
    representation(
        call = "call",
        parameters = "list",
        price = "numeric",     
        title = "character",
        description = "character"
    ) 
)
   

################################################################################


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, title = NULL,
description = NULL)
{   # 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) 
        
    # Parameters:
    param = list()
    param$TypeFlag = TypeFlag
    param$S = S
    param$X = X
    param$Time = Time
    param$r = r
    param$b = b
    param$sigma = sigma
    
    # Add title and description:
    if (is.null(title)) title = "Black Scholes Option Valuation"
    if (is.null(description)) description = as.character(date())
    
    # Return Value:
    new("fOPTION", 
        call = match.call(),
        parameters = param,
        price = result, 
        title = title,
        description = description
        )        
}


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


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, title = NULL,
description = NULL)
{   # 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]
    
    # Result:
    result = GBSOption(TypeFlag = TypeFlag, S = FT, X = X, Time = Time, 
        r = r, b = 0, sigma = sigma)@price
    
    # Parameters:
    param = list()
    param$TypeFlag = TypeFlag
    param$FT = FT
    param$X = X
    param$Time = Time
    param$r = r
    param$sigma = sigma
    
    # Add title and description:
    if (is.null(title)) title = "Black 76 Option Valuation"
    if (is.null(description)) description = as.character(date())
    
    # Return Value:
    new("fOPTION", 
        call = match.call(),
        parameters = param,
        price = result, 
        title = title,
        description = description
        )            
}


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


MiltersenSchwartzOption = 
function (TypeFlag = c("c", "p"), Pt, FT, X, time, Time, sigmaS, sigmaE, 
sigmaF, rhoSE, rhoSF, rhoEF, KappaE, KappaF, title = NULL,
description = NULL) 
{   # 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)) }
    
    # Parameters:
    param = list()
    param$TypeFlag = TypeFlag
    param$Pt = Pt
    param$FT = FT
    param$X = X
    param$time = time
    param$Time = Time
    param$sigmaS = sigmaS
    param$sigmaE = sigmaE
    param$sigmaF = sigmaF
    param$rhoSE = rhoSE
    param$rhoSF = rhoSF
    param$rhoEF = rhoEF
    param$KappaE = KappaE
    param$KappaF = KappaF
    
    # Add title and description:
    if (is.null(title)) title = "Miltersen Schwartz Option Valuation"
    if (is.null(description)) description = as.character(date())
    
    # Return Value:
    new("fOPTION", 
        call = match.call(),
        parameters = param,
        price = result, 
        title = title,
        description = description
        )        
}  


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


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, ...)
}


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


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

    # Description:
    #   Print method for objects of class "fOPTION".
    
    # FUNCTION:
    
    # Print Method:
    object = x
    Parameter = unlist(x@parameters)
    Names = names(Parameter)
    Parameter = cbind(as.character(Parameter))
    rownames(Parameter) = paste("", Names)
    colnames(Parameter) = "Value:"
    
    # Title:
    cat("\nTitle:\n ")
    cat(x@title, "\n")
    
    # Call:
    cat("\nCall:", paste("", deparse(x@call)), "", sep = "\n")
    
    # Parameters:
    cat("Parameters:\n")
    print(Parameter, quote = FALSE)
    
    # Price:
    cat("\nOption Price:\n ")
    cat(x@price, "\n")
    
    # Description:
    cat("\nDescription:\n ")
    cat(x@description, "\n\n")  
    
    # Return Value:
    invisible()
}


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


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

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


################################################################################