# 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:
#  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]
                
    # Integrate:
    call1 = integrate(.fstarHN, 0, Inf, const = 1, model = model, 
        S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
    # For SPlus Compatibility:
    if (is.null(call1$value)) call1$value = call1$integral
    call2 = integrate(.fstarHN, 0, Inf, const = 0, model = model, 
        S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
    # For SPlus Compatibility:
    if (is.null(call2$value)) call2$value = call2$integral
        
    # 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
}


.fstarHN <- 
function(phi, const, model, S, X, Time.inDays, r.daily) 
{          
    # Internal Function:

    # 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 
}


# ------------------------------------------------------------------------------
    
    
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]
            
    # Delta:
    if (Selection == "Delta") {
        # Integrate:
        delta1 = integrate(.fdeltaHN, 0, Inf, const = 1, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        # For SPlus Compatibility:
        if (is.null(delta1$value)) delta1$value = delta1$integral   
        delta2 = integrate(.fdeltaHN, 0, Inf, const = 0, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily) 
        # For SPlus Compatibility:
        if (is.null(delta2$value)) delta2$value = delta2$integral   
        # 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") {
        # Integrate:    
        gamma1 = integrate(.fgammaHN, 0, Inf, const = 1, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        # For SPlus Compatibility:
        if (is.null(gamma1$value)) gamma1$value = gamma1$integral    
        gamma2 = integrate(.fgammaHN, 0, Inf, const = 0, model = model, 
            S = S, X = X, Time.inDays = Time.inDays, r.daily = r.daily)
        # For SPlus Compatibility:
        if (is.null(gamma2$value)) gamma2$value = gamma2$integral   
        # 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
}


.fdeltaHN <- 
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 * 
        .fHN(phi, const, model, S, X, Time.inDays, r.daily) / S
    
    # Return Value:
    Re(fdelta) 
}


.fgammaHN <- 
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 ) *
        .fHN(phi, const, model, S, X, Time.inDays, r.daily) / S^2
    
    # Return Value:
    Re(fgamma) 
}



.fHN <- 
function(phi, const, model, S, X, Time.inDays, r.daily) 
{           
    # Internal Function:
    
    # 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) }
    fun = exp(-cphi0*log(X)+cphi*log(S)+a+b*sigma2)/cphi0/pi
    
    # Return Value:
    fun 
}


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


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


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