R/GpdFit.R

In fExtremes: Rmetrics - Modelling Extreme Events in Finance

# 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


################################################################################
# FUNCTION:               GPD SIMULATION:
#  gpdSim                  Simulates a GPD distributed process
# FUNCTION:               GPD PARAMETER ESTIMATION:
#  'fGPDFIT'               S4 class representation
#  gpdFit                  Fits Parameters of GPD distribution
#  .gpdpwmFit              Fits GPD with probability weighted moments
#  .gpdmleFit              Fits GPD with max log-likelihood approach
#   .gpdLLH                 Computes GPD log-likelihood function
################################################################################


gpdSim <-
    function(model = list(xi = 0.25, mu = 0, beta = 1), n = 1000, seed = NULL)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Generates random variates from a GPD distribution

    # FUNCTION:

    # Seed:
    if (is.null(seed)) seed = NA else set.seed(seed)

    # Simulate:
    ans = rgpd(n = n, xi = model$xi, mu = model$mu, beta = model$beta)
    ans = as.ts(ans)

    # Control:
    attr(ans, "control") =
        data.frame(t(unlist(model)), seed = seed, row.names = "")

    # Return Value:
    ans
}


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


setClass("fGPDFIT",
    representation(
        call = "call",
        method = "character",
        parameter = "list",
        data = "list",
        fit = "list",
        residuals = "numeric",
        title = "character",
        description = "character"
    )
)


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


gpdFit <-
    function(x, u = quantile(x, 0.95), type = c("mle", "pwm"),
    information = c("observed", "expected"),
    title = NULL, description = NULL, ...)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits a generalized Pareto model to excesses

    # Arguments:

    # Details:
    #   Returns an object of class "fGPDFIT" representing the fit of a
    #   generalized Pareto model to excesses over a high threshold

    # Notes:
    #   This is a wrapper to EVIR's 'gpd' function.

    # FUNCTION:

    # Settings:
    call = match.call()
    type = match.arg(type)
    information = match.arg(information)

    # Check Type and Convert:
    X = x
    xClass = class(x)
    if (xClass == "timeSeries") stopifnot(isUnivariate(x))
    x = as.vector(x)
    N = length(x)

    # Compute Exceedances:
    exceedances = x[x > u]
    Names = as.character((1:N)[x > u])
    exceedances = as.vector(exceedances)
    names(exceedances) = Names

    # Estimate Parameters:
    if (type == "mle") {
        fit = .gpdmleFit(x, u, information)
        fit$llh = fit$fit$value
        fit$convergence = fit$fit$convergence
    } else if (type == "pwm") {
        fit = .gpdpwmFit(x, u)
        fit$llh = NA
        fit$convergence = NA
    }
    fit$p.less.thresh = fit$prob = 1 - length(x[x > u]) / length(x)
    fit$threshold = u
    fit$data = x

    # Compute Residuals:
    xi = fit$par.ests["xi"]
    beta = fit$par.ests["beta"]
    residuals = log(1 + (xi * (as.vector(exceedances)-u))/beta)/xi

    # Add title and description:
    if (is.null(title)) title = "GPD Parameter Estimation"
    if (is.null(description)) description = description()

    # Compose Parameter List:
    parameter = list(u = u, type = type)
    if (information == "mle") parameter$information = information

    # Return Value:
    new("fGPDFIT",
        call = call,
        method = c("gpd", type),
        parameter = parameter,
        data = list(x = X, exceedances = exceedances),
        fit = fit,
        residuals = residuals,
        title = title,
        description = description)
}


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


.gpdmleFit <-
    function (x, u, information = c("observed", "expected"), ...)
{
    # A Copy from Evir

    data = x
    threshold = u

    exceedances <- data[data > threshold]
    excess <- exceedances - threshold
    Nu <- length(excess)
    xbar <- mean(excess)

    s2 <- var(excess)
    xi0 <- -0.5 * (((xbar * xbar)/s2) - 1)
    beta0 <- 0.5 * xbar * (((xbar * xbar)/s2) + 1)
    theta <- c(xi0, beta0)

    negloglik <- function(theta, tmp)
    {
        xi <- theta[1]
        beta <- theta[2]
        cond1 <- beta <= 0
        cond2 <- (xi <= 0) && (max(tmp) > (-beta/xi))
        if (cond1 || cond2) {
            f <- 1e+06
        } else {
            y <- logb(1 + (xi * tmp)/beta)
            y <- y/xi
            f <- length(tmp) * logb(beta) + (1 + xi) * sum(y)
        }
        f
    }

    fit <- optim(theta, negloglik, hessian = TRUE, ..., tmp = excess)
    names(fit$par) = c("xi", "beta")

    if (fit$convergence) warning("Optimization may not have been succeeded.")
    par.ests <- fit$par
    converged <- fit$convergence
    nllh.final <- fit$value

    information <- match.arg(information)
    if (information == "observed")
        varcov <- solve(fit$hessian)
    if (information == "expected") {
        one <- (1 + par.ests[1])^2/Nu
        two <- (2 * (1 + par.ests[1]) * par.ests[2]^2)/Nu
        cov <- -((1 + par.ests[1]) * par.ests[2])/Nu
        varcov <- matrix(c(one, cov, cov, two), 2)
    }

    par.ses <- sqrt(diag(varcov))
    names(par.ses) = c("xi", "beta")

    list(par.ests = par.ests, par.ses = par.ses, fit = fit, varcov = varcov)
}


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


.gpdmleFitCheck =
    function(x, u = quantile(x, 0.95),
    information = c("observed", "expected"), ...)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits GPD with max log-likelihood approach

    # FUNCTION:

    x = as.vector(x)
    excess = x[x > u] - u
    theta = .gpdpwmFit(x = x, u = u)$par.ests

    # Parameter Estimation:
    fit = optim(theta, .gpdLLH, hessian = TRUE, excess = excess)
    names(fit$par.ests) = c("xi", "beta")

    # Error Estimates:
    if (information[1] == "observed") {
        varcov = solve(fit$hessian)
    }
    if (information[1] == "expected") {
        Nu = length(excess)
        one = (1 + fit$par[1])^2/Nu
        two = (2 * (1 + fit$par[1]) * fit$par[2]^2)/Nu
        cov = -((1 + fit$par[1]) * fit$par[2])/Nu
        varcov = matrix(c(one, cov, cov, two), 2)
    }
    par.ses = sqrt(diag(varcov))
    names(par.ses) = c("xi", "beta")

    # Return Value:
    list(par.ests = fit$par, par.ses = par.ses, fit = fit, varcov = varcov)
}


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


.gpdpwmFit <-
    function (x, u)
{
    # A Copy from Evir

    data = x
    threshold = u

    data <- as.numeric(data)
    n <- length(data)
    exceedances <- data[data > threshold]
    excess <- exceedances - threshold
    Nu <- length(excess)
    xbar <- mean(excess)

    a0 <- xbar
    gamma <- -0.35
    delta <- 0
    pvec <- ((1:Nu) + gamma)/(Nu + delta)
    a1 <- mean(sort(excess) * (1 - pvec))
    xi <- 2 - a0/(a0 - 2 * a1)
    beta <- (2 * a0 * a1)/(a0 - 2 * a1)
    par.ests <- c(xi, beta)
    names(par.ests) = c("xi", "beta")

    denom <- Nu * (1 - 2 * xi) * (3 - 2 * xi)
    if (xi > 0.5) {
        denom <- NA
        warning("Asymptotic Standard Errors not available for PWM when xi>0.5.")
    }
    one <- (1 - xi) * (1 - xi + 2 * xi^2) * (2 - xi)^2
    two <- (7 - 18 * xi + 11 * xi^2 - 2 * xi^3) * beta^2
    cov <- beta * (2 - xi) * (2 - 6 * xi + 7 * xi^2 - 2 * xi^3)
    varcov <- matrix(c(one, cov, cov, two), 2)/denom
    information <- "expected"
    converged <- NA
    nllh.final <- NA
    par.ses <- sqrt(diag(varcov))
    names(par.ses) = c("xi", "beta")

    list(par.ests = par.ests, par.ses = par.ses, fit = NA, varcov = NA)
}



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


.gpdpwmFitCheck <-
    function(x, u = quantile(x, 0.95))
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Fits GPD with probability weighted moments

    # FUNCTION:

    # PWM Fit:
    x = as.vector(x)
    excess = x[x > u] - u
    Nu = length(excess)
    a0 = mean(excess)
    pvec = ((1:Nu) - 0.35)/Nu
    a1 = mean(sort(excess) * (1 - pvec))
    xi = 2 - a0/(a0 - 2 * a1)
    beta = (2 * a0 * a1)/(a0 - 2 * a1)
    par.ests = c(xi = xi, beta = beta)
    names(par.ests) = c("xi", "beta")
    par.ses = c(xi = NA, beta = NA)
    names(par.ses) = c("xi", "beta")

    # Return Value:
    list(par.ests = par.ests, par.ses = par.ses, fit = NA, varcov = NA)
}


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


.gpdLLH =
function(theta, excess)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes GPD log-likelihood function

    # FUNCTION:

    # LLH:
    xi = theta[1]
    beta = theta[2]
    cond = (beta <= 0) || ((xi <= 0) && (max(excess) > (-beta/xi)))
    if (cond) {
        func = NA
    } else {
        y = log(1+(xi*excess)/beta) / xi
        func = length(excess) * log(beta) + (1+xi)*sum(y)
    }

    # Return Value:
    func
}


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