.packageName <- "fExtremes"

# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# BUILIN - PACKAGE DESCRIPTION:
#  Package: evir
#  Version: 1.1
#  Date: 2004-05-05
#  Title: Extreme Values in R
#  Author: S original (EVIS) by Alexander McNeil 
#    <mcneil@math.ethz.ch>, R port by Alec 
#    Stephenson <alec_stephenson@hotmail.com>.
#  Maintainer: Alec Stephenson <alec_stephenson@hotmail.com>
#  Depends: R (>= 1.5.0)
#  Description: Functions for extreme value theory, which may be 
#    divided into the following groups; exploratory data analysis, 
#    block maxima, peaks over thresholds (univariate and bivariate), 
#    point processes, gev/gpd distributions.
#  License: GPL (Version 2 or above)
#  URL: http://www.maths.lancs.ac.uk/~stephena/
#  Packaged: Wed May  5 15:29:24 2004; stephena
################################################################################
# BUILIN - PACKAGE DESCRIPTION:
#  Package: ismev
#  Version: 1.1
#  Date: 2003/11/25
#  Title: An Introduction to Statistical Modeling of Extreme Values  
#  Author: Original S functions by Stuart Coles
#    <Stuart.Coles@bristol.ac.uk>, R port and R documentation files 
#    by Alec Stephenson <a.stephenson@lancaster.ac.uk>.
#  Maintainer: Alec Stephenson <a.stephenson@lancaster.ac.uk>
#  Depends: R (>= 1.5.0)
#  Description: Functions to support the computations carried out in
#    `An Introduction to Statistical Modeling of Extreme Values' by
#    Stuart Coles. The functions may be divided into the following 
#    groups; maxima/minima, order statistics, peaks over thresholds
#    and point processes.  
#  License: GPL (Version 2 or above)
#  URL: http://www.maths.lancs.ac.uk/~stephena/
################################################################################


# This file contains the following functions:
# gev.fit  gev.diag  gev.pp  gev.qq  gev.rl  gev.his
# gevf  gevq  gev.dens  gev.profxi  gev.prof

"gev.fit"<-
function(xdat, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, 
mulink = identity, siglink = identity, shlink = identity, show = TRUE, 
method = "Nelder-Mead", maxit = 10000, ...)
{
#
# obtains mles etc for gev distn
#
	z <- list()
        npmu <- length(mul) + 1
        npsc <- length(sigl) + 1
        npsh <- length(shl) + 1
	z$trans <- FALSE	# if maximization fails, could try
# changing in1 and in2 which are 
# initial values for minimization routine
	in2 <- sqrt(6 * var(xdat))/pi
	in1 <- mean(xdat) - 0.57722 * in2
	if(is.null(mul)) {
		mumat <- as.matrix(rep(1, length(xdat)))
		muinit <- in1
	}
	else {
		z$trans <- TRUE
		mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
		muinit <- c(in1, rep(0, length(mul)))
	}
	if(is.null(sigl)) {
		sigmat <- as.matrix(rep(1, length(xdat)))
		siginit <- in2
	}
	else {
		z$trans <- TRUE
		sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
		siginit <- c(in2, rep(0, length(sigl)))
	}
	if(is.null(shl)) {
		shmat <- as.matrix(rep(1, length(xdat)))
		shinit <- 0.1
	}
	else {
		z$trans <- TRUE
		shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
		shinit <- c(0.1, rep(0, length(shl)))
	}
	z$model <- list(mul, sigl, shl)
	z$link <- deparse(substitute(c(mulink, siglink, shlink)))
	init <- c(muinit, siginit, shinit)
        gev.lik <- function(a) {
        # computes neg log lik of gev model
        mu <- mulink(mumat %*% (a[1:npmu]))
        sc <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
	xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
	y <- (xdat - mu)/sc
	y <- 1 + xi * y
	if(any(y <= 0) || any(sc <= 0)) return(10^6)
	sum(log(sc)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
        }
	x <- optim(init, gev.lik, hessian = TRUE, method = method,
                   control = list(maxit = maxit, ...))
	z$conv <- x$convergence
        mu <- mulink(mumat %*% (x$par[1:npmu]))
	sc <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
	xi <- shlink(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
	z$nllh <- x$value
	z$data <- xdat
	if(z$trans) {
		z$data <-  - log(as.vector((1 + (xi * (xdat - mu))/sc)^(
			-1/xi)))
	}
	z$mle <- x$par
        z$cov <- solve(x$hessian)
	z$se <- sqrt(diag(z$cov))
	z$vals <- cbind(mu, sc, xi)
        if(show) {
	    if(z$trans)
		print(z[c(2, 3, 4)])
	    else print(z[4])
	    if(!z$conv)
                print(z[c(5, 7, 9)])
	}
	invisible(z)
}

"gev.diag"<-
function(z)
{
#
# produces diagnostic plots for output of
# gev.fit stored in z
#
	n <- length(z$data)
	x <- (1:n)/(n + 1)
	if(z$trans) {
      		oldpar <- par(mfrow = c(1, 2))
       		plot(x, exp( - exp( - sort(z$data))), xlab = 
       			"Empirical", ylab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Probability Plot")
       		plot( - log( - log(x)), sort(z$data), ylab = 
       			"Empirical", xlab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Quantile Plot (Gumbel Scale)")
       	}
       	else {
       		oldpar <- par(mfrow = c(2, 2))
       		gev.pp(z$mle, z$data)
       		gev.qq(z$mle, z$data)
       		gev.rl(z$mle, z$cov, z$data)
       		gev.his(z$mle, z$data)
       	}
       	par(oldpar)
       	invisible()
}

"gev.pp"<-
function(a, dat)
{
#
# sub-function for gev.diag
# produces probability plot
#
	plot((1:length(dat))/length(dat), gevf(a, sort(dat)), xlab = 
		"Empirical", ylab = "Model", main = "Probability Plot")
	abline(0, 1, col = 4)
}

"gev.qq"<-
function(a, dat)
{
#
# function called by gev.diag
# produces quantile plot
#
	plot(gevq(a, 1 - (1:length(dat)/(length(dat) + 1))), sort(dat), ylab = 
		"Empirical", xlab = "Model", main = "Quantile Plot")
	abline(0, 1, col = 4)
}

"gev.rl"<-
function(a, mat, dat)
{
#
# function called by gev.diag
# produces return level curve and 95 % confidence intervals
# on usual scale
#
	eps <- 1e-006
	a1 <- a
	a2 <- a
	a3 <- a
	a1[1] <- a[1] + eps
	a2[2] <- a[2] + eps
	a3[3] <- a[3] + eps
	f <- c(seq(0.01, 0.09, by = 0.01), 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 
		0.8, 0.9, 0.95, 0.99, 0.995, 0.999)
	q <- gevq(a, 1 - f)
	d1 <- (gevq(a1, 1 - f) - q)/eps
	d2 <- (gevq(a2, 1 - f) - q)/eps
	d3 <- (gevq(a3, 1 - f) - q)/eps
	d <- cbind(d1, d2, d3)
	v <- apply(d, 1, q.form, m = mat)
	plot(-1/log(f), q, log = "x", type = "n", xlim = c(0.1, 1000), ylim = c(
		min(dat, q), max(dat, q)), xlab = "Return Period", ylab = 
		"Return Level")
	title("Return Level Plot")
	lines(-1/log(f), q)
	lines(-1/log(f), q + 1.96 * sqrt(v), col = 4)
	lines(-1/log(f), q - 1.96 * sqrt(v), col = 4)
	points(-1/log((1:length(dat))/(length(dat) + 1)), sort(dat))
}

"gev.his"<-
function(a, dat)
{
#
# Plots histogram of data and fitted density
# for output of gev.fit stored in z
#
	h <- hist(dat, prob = TRUE, plot = FALSE)
	if(a[3] < 0) {
		x <- seq(min(h$breaks), min(max(h$breaks), (a[1] - a[2]/a[3] - 
			0.001)), length = 100)
	}
	else {
		x <- seq(max(min(h$breaks), (a[1] - a[2]/a[3] + 0.001)), max(h$
			breaks), length = 100)
	}
	y <- gev.dens(a, x)
	hist(dat, prob = TRUE, ylim = c(0, max(y)), xlab = "z", ylab = "f(z)", 
		main = "Density Plot")
	points(dat, rep(0, length(dat)))
	lines(x, y)
}

"gevf"<-
function(a, z)
{
#
# ancillary function calculates gev dist fnc
#
	if(a[3] != 0) exp( - (1 + (a[3] * (z - a[1]))/a[2])^(-1/a[3])) else 
			gum.df(z, a[1], a[2])
}

"gevq"<-
function(a, p)
{
	if(a[3] != 0)
		a[1] + (a[2] * (( - log(1 - p))^( - a[3]) - 1))/a[3]
	else gum.q(p, a[1], a[2])
}

"gev.dens"<-
function(a, z)
{
#
# evaluates gev density with parameters a at z
#
	if(a[3] != 0) (exp( - (1 + (a[3] * (z - a[1]))/a[2])^(-1/a[3])) * (1 + (
			a[3] * (z - a[1]))/a[2])^(-1/a[3] - 1))/a[2] else {
		gum.dens(c(a[1], a[2]), z)
	}
}

"gev.profxi"<-
function(z, xlow, xup, conf = 0.95, nint = 100)
{
#
# plots profile log-likelihood for shape parameter
# in gev model
#
	cat("If routine fails, try changing plotting interval", fill = TRUE)
	v <- numeric(nint)
	x <- seq(xup, xlow, length = nint)
	sol <- c(z$mle[1], z$mle[2])
        gev.plikxi <- function(a) {
        # computes profile neg log lik
        if (abs(xi) < 10^(-6)) {
                y <- (z$data - a[1])/a[2]
                if(a[2] <= 0) l <- 10^6
                else l <- length(y) * log(a[2]) + sum(exp(-y)) + sum(y)
        }
        else {
		y <- (z$data - a[1])/a[2]
		y <- 1 + xi * y
		if(a[2] <= 0 || any(y <= 0))
			l <- 10^6
		else l <- length(y) * log(a[2]) + sum(y^(-1/xi)) + sum(log(y
			)) * (1/xi + 1)
	}
	l
        }
	for(i in 1:nint) {
		xi <- x[i]
		opt <- optim(sol, gev.plikxi)
		sol <- opt$par ; v[i] <- opt$value
	}
	plot(x,  - v, type = "l", xlab = "Shape Parameter", ylab = 
		"Profile Log-likelihood")
	ma <-  - z$nllh
	abline(h = ma, col = 4)
	abline(h = ma - 0.5 * qchisq(conf, 1), col = 4)
	invisible()
}

"gev.prof"<-
function(z, m, xlow, xup, conf = 0.95, nint = 100)
{
#
# plots profile log likelihood for m 'year' return level
# in gev model
#
        if(m <= 1) stop("`m' must be greater than one")
	cat("If routine fails, try changing plotting interval", fill = TRUE)
	p <- 1/m
	v <- numeric(nint)
	x <- seq(xlow, xup, length = nint)
	sol <- c(z$mle[2], z$mle[3])
        gev.plik <- function(a) {
        # computes profile neg log lik
        if (abs(a[2]) < 10^(-6)) {
                mu <- xp + a[1] * log(-log(1 - p))
                y <- (z$data - mu)/a[1]
                if(is.infinite(mu) || a[1] <= 0) l <- 10^6
                else l <- length(y) * log(a[1]) + sum(exp(-y)) + sum(y)
        }
	else {
                mu <- xp - a[1]/a[2] * (( - log(1 - p))^( - a[2]) - 1)
		y <- (z$data - mu)/a[1]
		y <- 1 + a[2] * y
                if(is.infinite(mu) || a[1] <= 0 || any(y <= 0))
			l <- 10^6
		else l <- length(y) * log(a[1]) + sum(y^(-1/a[2])) + sum(log(
				y)) * (1/a[2] + 1)
	}
	l
        }
        for(i in 1:nint) {
                xp <- x[i]
		opt <- optim(sol, gev.plik)
		sol <- opt$par ; v[i] <- opt$value 
	}
	plot(x,  - v, type = "l", xlab = "Return Level", ylab = 
		" Profile Log-likelihood")
	ma <-  - z$nllh
	abline(h = ma, col = 4)
	abline(h = ma - 0.5 * qchisq(conf, 1), col = 4)
	invisible()
}







# This file contains the following functions:
# gpd.fitrange  gpd.fit  gpd.diag  gpd.pp  gpd.qq  gpd.rl
# gpd.his  gpdf  gpdq  gpdq2  gpd.dens  gpd.profxi  gpd.prof

"gpd.fitrange"<-
function(data, umin, umax, nint = 10, show = FALSE)
{
#
# computes mle's in gpd model, adjusted for threshold, 
# over range of threshold choices.
#
	m <- s <- up <- ul <- matrix(0, nrow = nint, ncol = 2)
	u <- seq(umin, umax, length = nint)
	for(i in 1:nint) {
		z <- gpd.fit(data, u[i], show = show)
		m[i,  ] <- z$mle
		m[i, 1] <- m[i, 1] - m[i, 2] * u[i]
		d <- matrix(c(1,  - u[i]), ncol = 1)
		v <- t(d) %*% z$cov %*% d
		s[i,  ] <- z$se
		s[i, 1] <- sqrt(v)
		up[i,  ] <- m[i,  ] + 1.96 * s[i,  ]
		ul[i,  ] <- m[i,  ] - 1.96 * s[i,  ]
	}
	names <- c("Modified Scale", "Shape")
	oldpar <- par(mfrow = c(2, 1))
	for(i in 1:2) {
		um <- max(up[, i])
		ud <- min(ul[, i])
		plot(u, m[, i], ylim = c(ud, um), xlab = "Threshold", ylab = 
			names[i], type = "b")
		for(j in 1:nint)
			lines(c(u[j], u[j]), c(ul[j, i], up[j, i]))
	}
        par(oldpar)
        invisible()
}

"gpd.fit"<-
function(xdat, threshold, npy = 365, ydat = NULL, sigl = NULL, shl = NULL, 
siglink = identity, shlink = identity, show = TRUE, method = "Nelder-Mead", 
maxit = 10000, ...)
{
# 
# obtains mles etc for gpd model
#
	z <- list()
        npsc <- length(sigl) + 1
	npsh <- length(shl) + 1
        n <- length(xdat)
	z$trans <- FALSE
	if(is.function(threshold))
            stop("`threshold' cannot be a function")
	u <- rep(threshold, length.out = n)
        if(length(unique(u)) > 1) z$trans <- TRUE
	xdatu <- xdat[xdat > u]
	xind <- (1:n)[xdat > u]
	u <- u[xind]
	in2 <- sqrt(6 * var(xdat))/pi
	in1 <- mean(xdat, na.rm = TRUE) - 0.57722 * in2
	if(is.null(sigl)) {
		sigmat <- as.matrix(rep(1, length(xdatu)))
		siginit <- in2
	}
	else {
		z$trans <- TRUE
		sigmat <- cbind(rep(1, length(xdatu)), ydat[xind, sigl])
		siginit <- c(in2, rep(0, length(sigl)))
	}
	if(is.null(shl)) {
		shmat <- as.matrix(rep(1, length(xdatu)))
		shinit <- 0.1
	}
	else {
		z$trans <- TRUE
		shmat <- cbind(rep(1, length(xdatu)), ydat[xind, shl])
		shinit <- c(0.1, rep(0, length(shl)))
	}
	init <- c(siginit, shinit)
	z$model <- list(sigl, shl)
	z$link <- deparse(substitute(c(siglink, shlink)))
        z$threshold <- threshold
	z$nexc <- length(xdatu)
	z$data <- xdatu	
        gpd.lik <- function(a) {
        # calculates gpd neg log lik
	sc <- siglink(sigmat %*% (a[seq(1, length = npsc)]))
	xi <- shlink(shmat %*% (a[seq(npsc + 1, length = npsh)]))
	y <- (xdatu - u)/sc
	y <- 1 + xi * y
	if(min(sc) <= 0)
		l <- 10^6
	else {
		if(min(y) <= 0)
			l <- 10^6
		else {
			l <- sum(log(sc)) + sum(log(y) * (1/xi + 1))
		}
	}
	l
        }
        x <- optim(init, gpd.lik, hessian = TRUE, method = method,
                   control = list(maxit = maxit, ...))
	sc <- siglink(sigmat %*% (x$par[seq(1, length = npsc)]))
	xi <- shlink(shmat %*% (x$par[seq(npsc + 1, length = npsh)]))
	z$conv <- x$convergence
	z$nllh <- x$value
	z$vals <- cbind(sc, xi, u)
	if(z$trans) {
		z$data <-  - log(as.vector((1 + (xi * (xdatu - u))/sc)^(-1/xi))
			)
	}
	z$mle <- x$par
	z$rate <- length(xdatu)/n
        z$cov <- solve(x$hessian)
	z$se <- sqrt(diag(z$cov))
	z$n <- n
	z$npy <- npy
	z$xdata <- xdat
        if(show) {
	    if(z$trans)
		print(z[c(2, 3)])
	    if(length(z[[4]]) == 1)
		print(z[4])
	    print(z[c(5, 7)])
	    if(!z$conv)
		print(z[c(8, 10, 11, 13)])
        }
	invisible(z)
}

"gpd.diag"<-
function(z)
{
#
# produces diagnostic plots for gpd model
# estimated using gpd.fit with output stored in z
#
	n <- length(z$data)
	x <- (1:n)/(n + 1)
       	if(z$trans) {
       		oldpar <- par(mfrow = c(1, 2))
       		plot(x, 1 - exp( - sort(z$data)), xlab = "Empirical", 
       			ylab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Probability Plot")
       		plot( - log(1 - x), sort(z$data), ylab = "Empirical", 
       			xlab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Quantile Plot (Exptl. Scale)")
       	}
       	else {
       		oldpar <- par(mfrow = c(2, 2))
       		gpd.pp(z$mle, z$threshold, z$data)
       		gpd.qq(z$mle, z$threshold, z$data)
       		gpd.rl(z$mle, z$threshold, z$rate, z$n, z$npy, z$cov, z$
       			data, z$xdata)
       		gpd.his(z$mle, z$threshold, z$data)
       	}
        par(oldpar)
       	invisible()
}

"gpd.pp"<-
function(a, u, dat)
{
# 
# function called by gpd.diag
# produces probability plot for gpd model
#
	plot((1:length(dat))/length(dat), gpdf(a, u, sort(dat)), xlab = 
		"Empirical", ylab = "Model", main = "Probability Plot")
	abline(0, 1, col = 4)
}

"gpd.qq"<-
function(a, u, dat)
{
#
# function called by gpd.diag
# produces quantile plot for gpd model
#
	plot(gpdq(a, u, 1 - (1:length(dat)/(length(dat) + 1))), sort(dat), ylab
		 = "Empirical", xlab = "Model", main = "Quantile Plot")
	abline(0, 1, col = 4)
}

"gpd.rl"<-
function(a, u, la, n, npy, mat, dat, xdat)
{
#
# function called by gpd.diag
# produces return level curve and 95% confidence intervals
# for fitted gpd model
	a <- c(la, a)
	eps <- 1e-006
	a1 <- a
	a2 <- a
	a3 <- a
	a1[1] <- a[1] + eps
	a2[2] <- a[2] + eps
	a3[3] <- a[3] + eps
	jj <- seq(-1, 3.75 + log10(npy), by = 0.1)
	m <- c(1/la, 10^jj)
	q <- gpdq2(a[2:3], u, la, m)
	d1 <- (gpdq2(a1[2:3], u, la, m) - q)/eps
	d2 <- (gpdq2(a2[2:3], u, la, m) - q)/eps
	d3 <- (gpdq2(a3[2:3], u, la, m) - q)/eps
	d <- cbind(d1, d2, d3)
	mat <- matrix(c((la * (1 - la))/n, 0, 0, 0, mat[1, 1], mat[1, 2], 0, 
		mat[2, 1], mat[2, 2]), nc = 3)
	v <- apply(d, 1, q.form, m = mat)
	plot(m/npy, q, log = "x", type = "n", xlim = c(0.1, max(m)/npy), ylim
		 = c(u, max(xdat, q[q > u - 1] + 1.96 * sqrt(v)[q > u - 1])), 
		xlab = "Return period (years)", ylab = "Return level", main = 
		"Return Level Plot")
	lines(m[q > u - 1]/npy, q[q > u - 1])
	lines(m[q > u - 1]/npy, q[q > u - 1] + 1.96 * sqrt(v)[q > u - 1], col
		 = 4)
	lines(m[q > u - 1]/npy, q[q > u - 1] - 1.96 * sqrt(v)[q > u - 1], col
		 = 4)
	nl <- n - length(dat) + 1
	sdat <- sort(xdat)
	points((1/(1 - (1:n)/(n + 1))/npy)[sdat > u], sdat[sdat > u])	
	#	points(1/(1 - (1:n)/(n + 1))/npy, 
#		sort(xdat))
#	abline(h = u, col = 3)
}

"gpd.his"<-
function(a, u, dat)
{
#
# function called by gpd.diag
# produces histogram and density plot
#
	h <- hist(dat, prob = TRUE, plot = FALSE)
	x <- seq(u, max(h$breaks), length = 100)
	y <- gpd.dens(a, u, x)
	hist(dat, prob = TRUE, ylim = c(0, max(y)), xlab = "x", ylab = "f(x)", 
		main = "Density Plot")
	lines(x, y, col = 4)
}

"gpdf"<-
function(a, u, z)
{
#
# ancillary function
# calculates gpd distribution function
#
	1 - (1 + (a[2] * (z - u))/a[1])^(-1/a[2])
}

"gpdq"<-
function(a, u, p)
u + (a[1] * (p^( - a[2])	#
# ancillary function
# computes gpd quantiles
#
 - 1))/a[2]

"gpdq2"<-
function(a, u, la, m)
{
#
# ancillary function
# calculates quantiles of gpd model
#
	u + (a[1] * ((m * la)^(a[2]) - 1))/a[2]
}

"gpd.dens"<-
function(a, u, z)
{
#
# ancillary function computes gpd density
#
	(1 + (a[2] * (z - u))/a[1])^(-1/a[2] - 1)/a[1]
}

"gpd.profxi"<-
function(z, xlow, xup, conf = 0.95, nint = 100)
{
#
# plots profile log likelihood for shape parameter
# in gpd model
#
	cat("If routine fails, try changing plotting interval", fill = TRUE)
	xdat <- z$data ; u <- z$threshold
	v <- numeric(nint)
	x <- seq(xup, xlow, length = nint)
	sol <- z$mle[1]
        gpd.plikxi <- function(a) {
        # calculates profile log lik
	if(abs(xi) < 10^(-4)) l <- length(xdat) * log(a) + sum(xdat - u)/a
		 else {
		y <- (xdat - u)/a
		y <- 1 + xi * y
                if(any(y <= 0) || a <= 0)
			l <- 10^6
		else l <- length(xdat) * log(a) + sum(log(y)) * (1/xi + 1)
	}
	l
        }
	for(i in 1:nint) {
		xi <- x[i]
		opt <- optim(sol, gpd.plikxi, method = "BFGS")
		sol <- opt$par ; v[i] <- opt$value 
	}
	plot(x,  - v, type = "l", xlab = "Shape Parameter", ylab = 
		"Profile Log-likelihood")
	ma <-  - z$nllh
	abline(h = ma, lty = 1)
	abline(h = ma - 0.5 * qchisq(conf, 1), lty = 1)
	invisible()
}



"gpd.prof"<-
function(z, m, xlow, xup, npy = 365, conf = 0.95, nint = 100)
{
#
# plots profile log-likelihood for m-year return level
# in gpd model
#
	cat("If routine fails, try changing plotting interval", fill = TRUE)
        xdat <- z$data ; u <- z$threshold ; la <- z$rate
	v <- numeric(nint)
	x <- seq(xlow, xup, length = nint)
        m <- m * npy
	sol <- z$mle[2]
        gpd.plik <- function(a) {
        # calculates profile neg log lik
        if(m != Inf) sc <- (a * (xp - u))/((m * la)^a - 1) else sc <- (u - xp)/
			a
	if(abs(a) < 10^(-4))
		l <- length(xdat) * log(sc) + sum(xdat - u)/sc
	else {
		y <- (xdat - u)/sc
		y <- 1 + a * y
                if(any(y <= 0) || sc <= 0)
			l <- 10^6
		else l <- length(xdat) * log(sc) + sum(log(y)) * (1/a + 1)
	}
	l
        }
	for(i in 1:nint) {
		xp <- x[i]
		opt <- optim(sol, gpd.plik, method = "BFGS")
		sol <- opt$par ; v[i] <- opt$value
	}
	plot(x,  - v, type = "l", xlab = "Return Level", ylab = 
		"Profile Log-likelihood")
	ma <-  - z$nllh
	abline(h = ma)
	abline(h = ma - 0.5 * qchisq(conf, 1))
	invisible()
}





# This file contains the following functions:
# gum.fit  gum.diag  gum.rl  gum.df  gum.q  gum.dens

"gum.fit"<-
function(xdat, ydat = NULL, mul = NULL, sigl = NULL, mulink = identity, 
siglink = identity, show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
#
# finds mles etc for gumbel model
#
	z <- list()
        npmu <- length(mul) + 1
        npsc <- length(sigl) + 1
	z$trans <- FALSE
	in2 <- sqrt(6 * var(xdat))/pi
	in1 <- mean(xdat) - 0.57722 * in2
	if(is.null(mul)) {
		mumat <- as.matrix(rep(1, length(xdat)))
		muinit <- in1
	}
	else {
		z$trans <- TRUE
		mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
		muinit <- c(in1, rep(0, length(mul)))
	}
	if(is.null(sigl)) {
		sigmat <- as.matrix(rep(1, length(xdat)))
		siginit <- in2
	}
	else {
		z$trans <- TRUE
		sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
		siginit <- c(in2, rep(0, length(sigl)))
	}
	z$model <- list(mul, sigl)
	z$link <- c(deparse(substitute(mulink)), deparse(substitute(siglink)))
	init <- c(muinit, siginit)
        gum.lik <- function(a) {
        # calculates neg log lik of gumbel model
	mu <- mulink(mumat %*% (a[1:npmu]))
	sc <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
        if(any(sc <= 0)) return(10^6)
	y <- (xdat - mu)/sc
        sum(log(sc)) + sum(y) + sum(exp( - y))
        }
	x <- optim(init, gum.lik, hessian = TRUE, method = method,
                   control = list(maxit = maxit, ...))
	z$conv <- x$convergence
        if(!z$conv) {
                mu <- mulink(mumat %*% (x$par[1:npmu]))
	        sc <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
	        z$nllh <- x$value
	        z$data <- xdat
	        if(z$trans) {
		        z$data <- as.vector((xdat - mu)/sc)
	        }
	        z$mle <- x$par
                z$cov <- solve(x$hessian)
	        z$se <- sqrt(diag(z$cov))
	        z$vals <- cbind(mu, sc)
        }
        if(show) {
	    if(z$trans)
		print(z[c(2, 3, 4)])
	    else print(z[4])
	    if(!z$conv)
                print(z[c(5, 7, 9)])
        }
	invisible(z)
}

"gum.diag"<-
function(z)
{
#
# produces diagnostic plots for output of
# gum.fit stored in z
#
	z$mle <- c(z$mle, 0)
	n <- length(z$data)
	x <- (1:n)/(n + 1)
	if(z$trans) {
	        oldpar <- par(mfrow = c(1, 2))
	        plot(x, exp( - exp( - sort(z$data))), xlab = "empirical",
                     ylab = "model")
	       	abline(0, 1, col = 4)
	       	title("Residual Probability Plot")
	       	plot( - log( - log(x)), sort(z$data), xlab = 
	       		"empirical", ylab = "model")
	       	abline(0, 1, col = 4)
	       	title("Residual Quantile Plot (Gumbel Scale)")
	}
       	else {
       		oldpar <- par(mfrow = c(2, 2))
       		gev.pp(z$mle, z$data)
       		gev.qq(z$mle, z$data)
       		gum.rl(z$mle, z$cov, z$data)
       		gev.his(z$mle, z$data)
       	}
       	par(oldpar)
       	invisible()
}

"gum.rl"<-
function(a, mat, dat)
{
#
# function called by gum.diag
# produces return level curve and 95 % confidence intervals
# on usual scale for gumbel model
#
	eps <- 1e-006
	a1 <- a
	a2 <- a
	a1[1] <- a[1] + eps
	a2[2] <- a[2] + eps
	f <- c(seq(0.01, 0.09, by = 0.01), 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 
		0.8, 0.9, 0.95, 0.99, 0.995, 0.999)
	q <- gevq(a, 1 - f)
	d1 <- (gevq(a1, 1 - f) - q)/eps
	d2 <- (gevq(a2, 1 - f) - q)/eps
	d <- cbind(d1, d2)
	v <- apply(d, 1, q.form, m = mat)
	plot(-1/log(f), q, log = "x", type = "n", xlim = c(0.1, 1000), ylim = c(
		min(dat, q), max(dat, q)), xlab = "Return Period", ylab = 
		"Return Level")
	title("Return Level Plot")
	lines(-1/log(f), q)
	lines(-1/log(f), q + 1.96 * sqrt(v), col = 4)
	lines(-1/log(f), q - 1.96 * sqrt(v), col = 4)
	points(-1/log((1:length(dat))/(length(dat) + 1)), sort(dat))
}

"gum.df"<-
function(x, a, b)
{
#
# ancillary function calculates dist fnc of gumbel model
#
	exp( - exp( - (x - a)/b))
}

"gum.q"<-
function(x, a, b)
{
#
# ancillary routine
# calculates quantiles of gumbel distn
#
	a - b * log( - log(1 - x))
}

"gum.dens"<-
function(a, x)
{
#
# ancillary function calculates density for gumbel model
#
	y <- (x - a[1])/a[2]
	(exp( - y) * exp( - exp( - y)))/a[2]
}






# This file contains the following functions:
# identity  q.form  mrl.plot

"identity"<-
function(x)
x

"q.form"<-
function(d, m)
{
#
# ancillary routine
# evaluates quadratic forms
#
	t(as.matrix(d)) %*% m %*% as.matrix(d)
}

"mrl.plot"<-
function(data, umin = min(data), umax = max(data) - 0.1, conf = 0.95, nint = 
	100)
{
#
# function to produce empirical mean residual life plot
# as function of threshold.
# confidence intervals included as well.
#
	x <- xu <- xl <- numeric(nint)
	u <- seq(umin, umax, length = nint)
	for(i in 1:nint) {
		data <- data[data > u[i]]
		x[i] <- mean(data - u[i])
		sdev <- sqrt(var(data))
		n <- length(data)
		xu[i] <- x[i] + (qnorm((1 + conf)/2) * sdev)/sqrt(n)
		xl[i] <- x[i] - (qnorm((1 + conf)/2) * sdev)/sqrt(n)
	}
	plot(u, x, type = "l", xlab = "u", ylab = "Mean Excess", ylim = c(min(
		xl[!is.na(xl)]), max(xu[!is.na(xu)])))
	lines(u[!is.na(xl)], xl[!is.na(xl)], lty = 2)
	lines(u[!is.na(xu)], xu[!is.na(xu)], lty = 2)
}

# This file contains the following functions:
# pp.fitrange  pp.fit  pp.diag  pp.pp  pp.qq
# ppf  ppq  ppp

"pp.fitrange"<-
function(data, umin, umax, npy = 365, nint = 10, show = FALSE)
{
#
# produces estimates and 95% confidence intervals
# for point process model across range of thresholds
#
        m <- s <- up <- ul <- matrix(0, nrow = nint, ncol = 3)
	u <- seq(umin, umax, length = nint)
	for(i in 1:nint) {
		z <- pp.fit(data, u[i], npy, show = show)
		m[i,  ] <- z$mle
		s[i,  ] <- z$se
		up[i,  ] <- z$mle + 1.96 * z$se
		ul[i,  ] <- z$mle - 1.96 * z$se
	}
	names <- c("Location", "Scale", "Shape")
	oldpar <- par(mfrow = c(1, 3))
	for(i in 1:3) {
		um <- max(up[, i])
		ud <- min(ul[, i])
		plot(u, m[, i], ylim = c(ud, um), xlab = "Threshold", ylab = 
			names[i], type = "b")
		for(j in 1:nint)
			lines(c(u[j], u[j]), c(ul[j, i], up[j, i]))
	}
        par(oldpar)
        invisible()
}

"pp.fit"<-
function(xdat, threshold, npy = 365, ydat = NULL, mul = NULL, sigl = NULL, 
shl = NULL, mulink = identity, siglink = identity, shlink = identity, 
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
	z <- list()
        npmu <- length(mul) + 1
	npsc <- length(sigl) + 1
	npsh <- length(shl) + 1
        n <- length(xdat)
	z$trans <- FALSE
	if(is.function(threshold)) 
            stop("`threshold' cannot be a function")
	u <- rep(threshold, length.out = n)
	if(length(unique(u)) > 1) z$trans <- TRUE
	xdatu <- xdat[xdat > u]
	xind <- (1:n)[xdat > u]
	u <- u[xind]
	in2 <- sqrt(6 * var(xdat))/pi
	in1 <- mean(xdat) - 0.57722 * in2
	if(is.null(mul)) {
		mumat <- as.matrix(rep(1, length(xdatu)))
		muinit <- in1
	}
	else {
		z$trans <- TRUE
		mumat <- cbind(rep(1, length(xdatu)), ydat[xind, mul])
		muinit <- c(in1, rep(0, length(mul)))
	}
	if(is.null(sigl)) {
		sigmat <- as.matrix(rep(1, length(xdatu)))
		siginit <- in2
	}
	else {
		z$trans <- TRUE
		sigmat <- cbind(rep(1, length(xdatu)), ydat[xind, sigl])
		siginit <- c(in2, rep(0, length(sigl)))
	}
	if(is.null(shl)) {
		shmat <- as.matrix(rep(1, length(xdatu)))
		shinit <- 0.1
	}
	else {
		z$trans <- TRUE
		shmat <- cbind(rep(1, length(xdatu)), ydat[xind, shl])
		shinit <- c(0.1, rep(0, length(shl)))
	}
	init <- c(muinit, siginit, shinit)
	z$model <- list(mul, sigl, shl)
	z$link <- deparse(substitute(c(mulink, siglink, shlink)))
        z$threshold <- threshold
	z$npy <- npy
	z$nexc <- length(xdatu)
	z$data <- xdatu
        pp.lik <- function(a) {
	mu <- mulink(mumat %*% (a[1:npmu]))
	sc <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
	xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
        if(any(sc <= 0)) return(10^6)
	if(min(1 + ((xi * (u - mu))/sc)) < 0) {
		l <- 10^6
	}
	else {
		y <- (xdatu - mu)/sc
		y <- 1 + xi * y
		if(min(y) <= 0)
			l <- 10^6
		else l <- sum(log(sc)) + sum(log(y) * (1/xi + 1)) + n/npy * 
				mean((1 + (xi * (u - mu))/sc)^(-1/xi))
	}
	l
        }
	x <- optim(init, pp.lik, hessian = TRUE, method = method,
                   control = list(maxit = maxit, ...))
        mu <- mulink(mumat %*% (x$par[1:npmu]))
	sc <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
	xi <- shlink(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
	z$conv <- x$convergence
	z$nllh <- x$value
	z$vals <- cbind(mu, sc, xi, u)
	z$gpd <- apply(z$vals, 1, ppp, npy)
	if(z$trans) {
		z$data <- as.vector((1 + (xi * (xdatu - u))/z$gpd[2,  ])^(-1/xi
			))
	}
	z$mle <- x$par
        z$cov <- solve(x$hessian)
	z$se <- sqrt(diag(z$cov))
        if(show) {
	    if(z$trans)
		print(z[c(2, 3)])
	    if(length(z[[4]]) == 1)
		print(z[4])
	    print(z[c(5, 6, 8)])
	    if(!z$conv)
		print(z[c(9, 12, 14)])
        }
        invisible(z)
}

"pp.diag"<-
function(z)
{
	n <- length(z$data)
	x <- (1:n)/(n + 1)
	if(z$trans) {
		oldpar <- par(mfrow = c(1, 2))
		plot(x, sort(z$data), xlab = "empirical", ylab = "model")
		abline(0, 1, col = 3)
		title("Residual Probability Plot")
		plot( - log(1 - x),  - log(1 - sort(z$data)), ylab = 
			"empirical", xlab = "model")
		abline(0, 1, col = 3)
		title("Residual quantile Plot (Exptl. Scale)")
	}
	else {
		oldpar <- par(mfrow = c(1, 2), pty = "s")
		pp.pp(z$mle, z$threshold, z$npy, z$data)
		pp.qq(z$mle, z$threshold, z$npy, z$data)
	}
	par(oldpar)
	invisible()
}

"pp.pp"<-
function(a, u, npy, dat)
{
#
# function called by pp.diag
# produces probability plot
#
	y <- apply(as.matrix(sort(dat)), 1, ppf, a = a, u = u, npy = npy)
	plot((1:length(dat))/length(dat), y, xlab = "empirical", ylab = "model",
		main = "Probability plot")
	abline(0, 1, col = 4)
}

"pp.qq"<-
function(a, u, npy, dat)
{
#
# function called by pp.diag
# computes quantile plot
#
	y <- apply(as.matrix((length(dat):1/(length(dat) + 1))), 1, ppq, a = a, 
		u = u, npy = npy)
	plot(y, sort(dat), ylab = "empirical", xlab = "model", main = 
		"Quantile Plot")
	abline(0, 1, col = 4)
}

"ppf"<-
function(a, z, u, npy)
{
#
# ancillary function
# calculates distribution function in point process model
#
	b <- ppp(c(a, u), npy)
	1 - (1 + (b[3] * (z - u))/b[2])^(-1/b[3])
}

"ppq"<-
function(a, u, npy, p)
{
#
# ancillary routine
# finds quantiles in point process model
#
	b <- ppp(c(a, u), npy)
	u + (b[2] * (((p))^( - b[3]) - 1))/b[3]
}

"ppp"<-
function(a, npy)
{
	u <- a[4]
	la <- 1 - exp( - (1 + (a[3] * (u - a[1]))/a[2])^(-1/a[3])/npy)
	sc <- a[2] + a[3] * (u - a[1])
	xi <- a[3]
	c(la, sc, xi)
}









# This file contains the following functions:
# rlarg.fit  rlarg.diag  rlarg.pp  rlarg.qq
# rlargf  rlargq  rlargq2

"rlarg.fit"<-
function(xdat, r = dim(xdat)[2], ydat = NULL, mul = NULL, sigl = NULL, 
shl = NULL, mulink = identity, siglink = identity, shlink = identity, 
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
#
# calculates mles etc for rlargest order statistic model
#
	z <- list()
        npmu <- length(mul) + 1
        npsc <- length(sigl) + 1
        npsh <- length(shl) + 1
        z$trans <- FALSE
	in2 <- sqrt(6 * var(xdat[, 1]))/pi
	in1 <- mean(xdat[, 1]) - 0.57722 * in2
	if(is.null(mul)) {
		mumat <- as.matrix(rep(1, dim(xdat)[1]))
		muinit <- in1
	}
	else {
		z$trans <- TRUE
		mumat <- cbind(rep(1, dim(xdat)[1]), ydat[, mul])
		muinit <- c(in1, rep(0, length(mul)))
	}
	if(is.null(sigl)) {
		sigmat <- as.matrix(rep(1, dim(xdat)[1]))
		siginit <- in2
	}
	else {
		z$trans <- TRUE
		sigmat <- cbind(rep(1, dim(xdat)[1]), ydat[, sigl])
		siginit <- c(in2, rep(0, length(sigl)))
	}
	if(is.null(shl)) {
		shmat <- as.matrix(rep(1, dim(xdat)[1]))
		shinit <- 0.1
	}
	else {
		z$trans <- TRUE
		shmat <- cbind(rep(1, dim(xdat)[1]), ydat[, shl])
		shinit <- c(0.1, rep(0, length(shl)))
	}
        xdatu <- xdat[, 1:r, drop = FALSE]
        init <- c(muinit, siginit, shinit)
	z$model <- list(mul, sigl, shl)
	z$link <- deparse(substitute(c(mulink, siglink, shlink)))
        u <- apply(xdatu, 1, min, na.rm = TRUE)
        rlarg.lik <- function(a) {
        # calculates neg log lik
	mu <- mulink(drop(mumat %*% (a[1:npmu])))
	sc <- siglink(drop(sigmat %*% (a[seq(npmu + 1, length = npsc)])))
	xi <- shlink(drop(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)])))
        if(any(sc <= 0)) return(10^6)
        y <- 1 + xi * (xdatu - mu)/sc
	if(min(y, na.rm = TRUE) <= 0)
		l <- 10^6
	else {
                y <- (1/xi+1) * log(y) + log(sc)
                y <- rowSums(y, na.rm = TRUE)
                l <- sum((1 + xi * (u - mu)/sc)^(-1/xi) + y)
        }
	l
        }
	x <- optim(init, rlarg.lik, hessian = TRUE, method = method,
                   control = list(maxit = maxit, ...))
        mu <- mulink(drop(mumat %*% (x$par[1:npmu])))
	sc <- siglink(drop(sigmat %*% (x$par[seq(npmu + 1, length = npsc)])))
	xi <- shlink(drop(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)])))
	z$conv <- x$convergence
	z$nllh <- x$value
	z$data <- xdat
	if(z$trans) {
		for(i in 1:r)
			z$data[, i] <-  - log((1 + (as.vector(xi) * (xdat[, i] - 
				as.vector(mu)))/as.vector(sc))^(-1/as.vector(xi
				)))
	}
	z$mle <- x$par
        z$cov <- solve(x$hessian)
	z$se <- sqrt(diag(z$cov))
	z$vals <- cbind(mu, sc, xi)
	z$r <- r
        if(show) {
	    if(z$trans)
		print(z[c(2, 3)])
	    print(z[4])
	    if(!z$conv)
		print(z[c(5, 7, 9)])
        }
	invisible(z)
}

"rlarg.diag"<-
function(z, n = z$r)
{
#
# takes output from rlarg.fit
# produces probability and quantile plots for
# each order statistic
#
	z2 <- z
	z2$data <- z$data[, 1]
        oldpar <- par(ask = TRUE, mfcol = c(2, 2))
	if(z$trans) {
		for(i in 1:n) {
			rlarg.pp(c(0, 1, 0), z$data[, 1:z$r], i)
			rlarg.qq(c(0, 1, 0), z$data[, 1:z$r], i)
		}
	}
	else {
		gev.diag(z2)
		for(i in 1:n) {
			rlarg.pp(z$mle, z$data, i)
			rlarg.qq(z$mle, z$data, i)
		}
	}
	par(oldpar)
	invisible()
}

"rlarg.pp"<-
function(a, dat, k)
{
#
# ancillary function
# calculates probability plot in r largest model
#
	da <- dat[!is.na(dat[, k]), k]
	plot((1:length(da))/length(da), rlargf(a, sort(da), k), xlab = "", ylab
		 = "")
	title(paste("k=", k, sep = ""), cex = 0.7)
	abline(0, 1, col = 4)
}

"rlarg.qq"<-
function(a, dat, k)
{
#
# ancillary function
# calculates quantile plot in r largest model
#
	da <- dat[!is.na(dat[, k]), k]
	plot(rlargq(a, 1 - (1:length(da)/(length(da) + 1)), k, da), sort(da), 
		xlab = "", ylab = "")
	title(paste("k=", k, sep = ""), cex = 0.7)
	abline(0, 1, col = 4)
}

"rlargf"<-
function(a, z, k)
{
#
# ancillary function
# calculates dist fnc in r largest model
#
	eps <- 10^(-6)
	res <- NULL
	if(abs(a[3]) < eps)
		tau <- exp( - (z - a[1])/a[2])
	else tau <- (1 + (a[3] * (z - a[1]))/a[2])^(-1/a[3])
	for(i in 1:length(tau)) {
		if(is.na(tau[i]))
			res[i] <- 1
		else res[i] <- exp( - tau[i]) * sum(tau[i]^(0:(k - 1))/gamma(1:(
				k)))
	}
	res
}

"rlargq"<-
function(a, p, k, dat)
{
#
# ancillary routine 
# for finding quantiles in r largest model
	res <- NULL
	for(i in 1:length(p)) {
		inter <- c(min(dat) - 1, max(dat) + 1)
		res[i] <- uniroot(rlargq2, inter, a = a, kk = k, p = p[i])$root
	}
	res
}

"rlargq2"<-
function(x, a, kk, p)
{
#
# ancillary routine
# for finding quantiles in r largest model
#
	res <- rlargf(a, x, kk) - (1 - p)
	res
}



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


"gev" <- 
function(data, block = NA, ...)
{
    n.all <- NA
    if(!is.na(block)) {
        n.all <- length(data)
        if(is.character(block)) {
            times <- as.POSIXlt(attributes(data)$times)
            if(block %in% c("semester", "quarter")) {
                sem <- quart <- times$mon
                sem[sem %in% 0:5] <- quart[quart %in% 0:2] <- 0
                sem[sem %in% 6:11] <- quart[quart %in% 3:5] <- 1
                quart[quart %in% 6:8] <- 2
                quart[quart %in% 9:11] <- 3
            }
            grouping <- switch(block,
                semester = paste(times$year, sem),
                quarter = paste(times$year, quart),
                month = paste(times$year, times$mon),
                year = times$year,
                stop("unknown time period"))
            data <- tapply(data, grouping, max)
        }
        else {
            data <- as.numeric(data)
            nblocks <- (length(data) %/% block) + 1
            grouping <- rep(1:nblocks, rep(block, nblocks))[1:length(data)]
            data <- tapply(data, grouping, max)
        }
    }
    data <- as.numeric(data)
    n <- length(data)
    sigma0 <- sqrt(6 * var(data))/pi
    mu0 <- mean(data) - 0.57722 * sigma0
    xi0 <- 0.1
    theta <- c(xi0, sigma0, mu0)
    negloglik <- function(theta, tmp)
    {
      	y <- 1 + (theta[1] * (tmp - theta[3]))/theta[2]
       	if((theta[2] < 0) || (min(y) < 0))
       	    out <- 1e+06
       	else {
       	    term1 <- length(tmp) * logb(theta[2])
       	    term2 <- sum((1 + 1/theta[1]) * logb(y))
       	    term3 <- sum(y^(-1/theta[1]))
       	    out <- term1 + term2 + term3
       	}
       	out
    }
    fit <- optim(theta, negloglik, hessian = TRUE, ..., tmp = data)
    if(fit$convergence)
        warning("optimization may not have succeeded")
    par.ests <- fit$par
    varcov <- solve(fit$hessian)
    par.ses <- sqrt(diag(varcov))
    out <- list(n.all = n.all, n = n, data = data, block = block, par.ests
       	 = par.ests, par.ses = par.ses, varcov = varcov, converged = 
       	fit$convergence, nllh.final = fit$value)
    names(out$par.ests) <- c("xi", "sigma", "mu")
    names(out$par.ses) <- c("xi", "sigma", "mu")
    class(out) <- "gev"
    out
}

"gumbel" <- 
function(data, block = NA, ...)
{
	n.all <- NA
	data <- as.numeric(data)
        if(!is.na(block)) {
	  n.all <- length(data)
	  if(fg <- n.all %% block) {
              data <- c(data, rep(NA, block - fg))
              warning(paste("final group contains only", fg, "observations"))
          }
          data <- apply(matrix(data, nrow = block), 2, max, na.rm = TRUE)
	}
	n <- length(data)
	sigma0 <- sqrt(6 * var(data))/pi
	mu0 <- mean(data) - 0.57722 * sigma0
	theta <- c(sigma0, mu0)
	negloglik <- function(theta, tmp)
	{
		y <- (tmp - theta[2])/theta[1]
		if(theta[1] < 0)
			out <- 1e+06
		else {
			term1 <- length(tmp) * logb(theta[1])
			term2 <- sum(y)
			term3 <- sum(exp( - y))
			out <- term1 + term2 + term3
		}
		out
	}
        fit <- optim(theta, negloglik, hessian = TRUE, ..., tmp = data)
        if(fit$convergence)
            warning("optimization may not have succeeded")
	par.ests <- fit$par
	varcov <- solve(fit$hessian)
	par.ses <- sqrt(diag(varcov))
	out <- list(n.all = n.all, n = n, data = data, block = block, par.ests
		 = par.ests, par.ses = par.ses, varcov = varcov, converged = 
		fit$convergence, nllh.final = fit$value)
	names(out$par.ests) <- c("sigma", "mu")
	names(out$par.ses) <- c("sigma", "mu")
	class(out) <- "gev"
	out
}

"plot.gev" <- 
function(x, ...)
{
	par.ests <- x$par.ests
	mu <- par.ests["mu"]
	sigma <- par.ests["sigma"]
	if(!("xi" %in% names(par.ests)))
	    xi <- 0
	else xi <- par.ests["xi"]
	if(xi != 0)
	    residuals <- (1 + (xi * (x$data - mu))/sigma)^(-1/xi)
	else residuals <- exp( - exp( - (x$data - mu)/sigma))
	choices <- c("Scatterplot of Residuals", "QQplot of Residuals")
	tmenu <- paste("plot:", choices)
	pick <- 1
	while(pick > 0) {
	    pick <- menu(tmenu, title =
                         "\nMake a plot selection (or 0 to exit):")
	    switch(pick,
		   {
		       plot(residuals, ylab = "Residuals",
                            xlab = "Ordering", ...)
		       lines(lowess(1:length(residuals), residuals))
		   },
		   qplot(residuals, ...))
	}
}

"rlevel.gev" <- 
function(out, k.blocks = 20, add = FALSE, ...)
{
	par.ests <- out$par.ests
	mu <- par.ests["mu"]
	sigma <- par.ests["sigma"]
	if(!("xi" %in% names(par.ests)))
	    stop("Use this function after a GEV rather than a Gumbel fit")
	else xi <- par.ests["xi"]
	pp <- 1/k.blocks
	v <- qgev((1 - pp), xi, mu, sigma)
	if(add) abline(h = v)
	data <- out$data
        overallmax <- out$nllh.final
	sigma0 <- sqrt(6 * var(data))/pi
	xi0 <- 0.01
	theta <- c(xi0, sigma0)
	parloglik <- function(theta, tmp, pp, rli)
	{
		mu <- rli + (theta[2] * (1 - ( - logb(1 - pp))^( - theta[
			1])))/theta[1]
		y <- 1 + (theta[1] * (tmp - mu))/theta[2]
		if((theta[2] < 0) | (min(y) < 0))
			out <- 1e+06
		else {
			term1 <- length(tmp) * logb(theta[2])
			term2 <- sum((1 + 1/theta[1]) * logb(y))
			term3 <- sum(y^(-1/theta[1]))
			out <- term1 + term2 + term3
		}
		out
	}
	parmax <- NULL
	rl <- v * c(0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1, 1.1, 1.2,
                    1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 4.5)
	for(i in 1:length(rl)) {
		fit <- optim(theta, parloglik, hessian = FALSE, tmp = data,
                             pp = pp, rli = rl[i])
		parmax <- rbind(parmax, fit$value)
	}
	parmax <-  - parmax
	overallmax <-  - overallmax
	crit <- overallmax - qchisq(0.9999, 1)/2
	cond <- parmax > crit
	rl <- rl[cond]
	parmax <- parmax[cond]
	smth <- spline(rl, parmax, n = 200)
	aalpha <- qchisq(0.95, 1)
	if(!add) {
	    plot(rl, parmax, type = "p", ...)
	    abline(h = overallmax - aalpha/2)
	    abline(v = v)
	    lines(smth)
	}
        ind <- smth$y > overallmax - aalpha/2
	ci <- range(smth$x[ind])
	if(add) {
	    abline(h = ci[1], lty = 2, col = 2)
	    abline(h = ci[2], lty = 2, col = 2)
	}
	as.numeric(c(ci[1], v, ci[2]))
}


"gpdbiv" <- 
function(data1 = NA, data2 = NA, u1 = NA, u2 = NA, ne1 = NA,
    ne2 = NA, global = FALSE, method = "BFGS", ...)
{
    data1 <- as.numeric(data1)
    data2 <- as.numeric(data2)
    
    Zfunc <- function(y, u, lambda, xi, sigma)
        (lambda^-1) * (1 + (xi * pmax((y - u), 0))/sigma)^(1/xi)
    Kfunc <- function(y, u, lambda, xi, sigma)
        -lambda^(-xi) * (sigma^-1) * (Zfunc(y, u, lambda, xi, sigma))^(1 - xi)
    Vfunc <- function(x, y, alpha)
        (x^(-1/alpha) + y^(-1/alpha))^alpha
    Vfunc1 <- function(x, y, alpha)
        -x^(-(1/alpha) - 1) * (x^(-1/alpha) + y^(-1/alpha))^(alpha - 1)
    Vfunc2 <- function(x, y, alpha)
        -(alpha - 1) * (alpha^-1) * (x * y)^(-(1/alpha) - 1) *
          (x^(-1/alpha) + y^(-1/alpha))^(alpha - 2)
    fun <- list(Z = Zfunc, K = Kfunc, V = Vfunc, V1 = Vfunc1, V2 = Vfunc2) 
   
    if(is.na(ne1) && is.na(u1))
        stop(paste("Enter either a threshold or",
                   "the number of upper extremes for margin 1"))
    if(!is.na(ne1) && !is.na(u1))
        stop(paste("Enter EITHER a threshold or",
                   "the number of upper extremes for margin 1"))
    if(is.na(ne2) && is.na(u2))
        stop(paste("Enter either a threshold or",
                   "the number of upper extremes for margin 2"))
    if(!is.na(ne2) && !is.na(u2))
        stop(paste("Enter EITHER a threshold or",
                   "the number of upper extremes for margin 2"))

    out1 <- gpd(data1, threshold = u1, ne = ne1)
    par.ests1 <- out1$par.ests
    par.ses1 <- out1$par.ses
    
    out2 <- gpd(data2, threshold = u2, ne = ne2)
    par.ests2 <- out2$par.ests
    par.ses2 <- out2$par.ses

    uu <- c(out1$threshold, out2$threshold)
    ne <- c(out1$n.exceed, out2$n.exceed)
    mpar <- c(par.ests1, par.ests2) 
    
    delta1 <- as.numeric(data1 > uu[1])
    delta2 <- as.numeric(data2 > uu[2])
    lambda1 <- sum(delta1)/length(data1)
    lambda2 <- sum(delta2)/length(data2)

    theta <- 0.8
    if(global) {
        theta <- c(theta, mpar)
        mpar <- NULL
    }
	
    negloglik <- function(theta, data1, data2, uu, delta1, delta2,
        lambda1, lambda2, mpar, fun)
    {
      	alpha <- theta[1]
	if(is.null(mpar)) {
            xi1 <- theta[2] ; sigma1 <- theta[3]
	    xi2 <- theta[4] ; sigma2 <- theta[5]
	}
        else {
            xi1 <- mpar[1] ; sigma1 <- mpar[2]
	    xi2 <- mpar[3] ; sigma2 <- mpar[4]
        }
	cond1 <- (alpha <= 0) | (alpha >= 1)
	cond2 <- sigma1 <= 0
	cond3 <- sigma2 <= 0
	if(cond1 || cond2 || cond3)
	   out <- 1e+06
	else {
	    term4 <- (1 - delta1) * (1 - delta2) * logb(1 -
                fun$V(lambda1^-1, lambda2^-1, alpha))
	    term3 <- delta1 * (1 - delta2) * logb(fun$K(data1, uu[1], lambda1,
                xi1, sigma1) * fun$V1(fun$Z(data1, uu[1], lambda1, xi1,
                sigma1), lambda2^-1, alpha))
	    term2 <- delta2 * (1 - delta1) * logb(fun$K(data2, uu[2], lambda2,
                xi2, sigma2) * fun$V1(fun$Z(data2, uu[2], lambda2, xi2,
                sigma2), lambda1^-1, alpha))
	    term1 <- delta1 * delta2 * logb(fun$K(data1, uu[1], lambda1, xi1,
                sigma1) * fun$K(data2, uu[2], lambda2, xi2, sigma2) *
                fun$V2(fun$Z(data1, uu[1], lambda1, xi1, sigma1), fun$Z(data2,
                uu[2], lambda2, xi2, sigma2), alpha))
	    allterm <- term1 + term2 + term3 + term4
	    out <-  - sum(allterm)
	}
	out
    }
    fit <- optim(theta, negloglik, hessian = TRUE, method = method, ...,
                 data1 = data1, data2 = data2, uu = uu,
                 delta1 = delta1, delta2 = delta2, lambda1 = lambda1,
                 lambda2 = lambda2, mpar = mpar, fun = fun)
    if(fit$convergence)
        warning("optimization may not have succeeded")
    par.ests <- fit$par
    varcov <- solve(fit$hessian)
    par.ses <- sqrt(diag(varcov))
    alpha <- par.ests[1]
    alpha.se <- par.ses[1]
    if(global) {
        par.ests1 <- c(par.ests[2], par.ests[3])
        names(par.ests1) <- c("xi", "beta")
        par.ses1 <- c(par.ses[2], par.ses[3])
        par.ests2 <- c(par.ests[4], par.ests[5])
        names(par.ests2) <- c("xi", "beta")
        par.ses2 <- c(par.ses[4], par.ses[5])
    }
    out <- list(data1 = data1[delta1 == 1], delta1 = (delta1 ==
      	1 & delta2 == 1)[delta1 == 1], data2 = data2[
       	delta2 == 1], delta2 = (delta1 == 1 & delta2 == 1)[delta2 ==
       	1], u1 = uu[1], ne1 = ne[1], lambda1 = lambda1, u2 = uu[2],
        ne2 = ne[2], lambda2 = lambda2, alpha = alpha, alpha.se = alpha.se, 
       	par.ests1 = par.ests1, par.ses1 = par.ses1, par.ests2 = 
       	par.ests2, par.ses2 = par.ses2, converged = fit$convergence,
       	nllh.final = fit$value, dependence = "logistic", 
       	dep.func = Vfunc)
    class(out) <- "gpdbiv"
    out
}

"interpret.gpdbiv" <- 
function(out, x, y)
{
    Vfuncf <- out$dep.func
    newfunc <- function(x, y, alpha, u1, lambda1, xi1, sigma1, u2, lambda2,
	           xi2, sigma2, vfunc)
    {
        Zfunc <- function(y, u, lambda, xi, sigma)
	    (lambda^-1) * (1 + (xi * pmax((y - u), 0))/sigma)^(1/xi)
	1 - vfunc(Zfunc(x, u1, lambda1, xi1, sigma1), Zfunc(y, u2,
		  lambda2, xi2, sigma2), alpha)
    }
    marg <- function(x, u1, lambda1, xi1, sigma1)
    {
        1 - lambda1 * (1 + (xi1 * (x - u1))/sigma1)^(-1/xi1)
    }
    newfunc2 <- function(x, y, alpha, u1, lambda1, xi1, sigma1, u2, lambda2,
		    xi2, sigma2, marg, newfunc, vfunc)
    {
        1 - marg(x, u1, lambda1, xi1, sigma1) - marg(y, u2, lambda2, xi2,
        sigma2) + newfunc(x, y, alpha, u1, lambda1, xi1, sigma1, u2,
        lambda2, xi2, sigma2, vfunc)
    }
    
    if(out$u1 > x) stop("Point below x threshold")
    if(out$u2 > y) stop("Point below y threshold")
    p1 <- 1 - marg(x, out$u1, out$lambda1, out$par.ests1[1], out$
		par.ests1[2])
    p2 <- 1 - marg(y, out$u2, out$lambda2, out$par.ests2[1], out$
		par.ests2[2])
    p12 <- newfunc2(x, y, out$alpha, out$u1, out$lambda1, out$par.ests1[1],
                    out$par.ests1[2], out$u2, out$lambda2, out$par.ests2[1],
                    out$par.ests2[2], marg, newfunc, Vfuncf)
    
    cat("Thresholds:", out$u1, out$u2, "\n")
    cat("Extreme levels of interest (x,y):", x, y, "\n")
    cat("P(X exceeds x)", p1, "\n")
    cat("P(Y exceeds y)", p2, "\n")
    cat("P(X exceeds x AND Y exceeds y)", p12, "\n")
    cat("P(X exceeds x) * P(Y exceeds y)", p1 * p2, "\n")
    cat("P(Y exceeds y GIVEN X exceeds x)", p12/p1, "\n")
    cat("P(X exceeds x GIVEN Y exceeds y)", p12/p2, "\n")
    invisible(as.numeric(c(p1, p2, p12, p1 * p2, p12/p1, p12/p2)))
}

"plot.gpdbiv" <- 
function(x, extend = 1.1, n.contours = 15, ...)
{
    Zfunc <- function(y, u, lambda, xi, sigma)
        (lambda^-1) * (1 + (xi * pmax((y - u), 0))/sigma)^(1/xi)

    joint <- function(xx, y, alpha, u1, lambda1, xi1, sigma1, u2, lambda2,
		xi2, sigma2, Vfunc)
    {
        1 - Vfunc(Zfunc(xx, u1, lambda1, xi1, sigma1),
                  Zfunc(y, u2, lambda2, xi2, sigma2), alpha)
    }
    marg <- function(xx, u1, lambda1, xi1, sigma1)
    {
        1 - lambda1 * (1 + (xi1 * (xx - u1))/sigma1)^(-1/xi1)
    }
    survivor <- function(xx, y, alpha, u1, lambda1, xi1, sigma1, u2,
                    lambda2, xi2, sigma2, marg, joint, Vfunc)
    {
	1 - marg(xx, u1, lambda1, xi1, sigma1) - marg(y, u2, lambda2,
	    xi2, sigma2) + joint(xx, y, alpha, u1, lambda1, xi1,
	    sigma1, u2, lambda2, xi2, sigma2, Vfunc)
    }
    
    xx <- seq(from = x$u1, to = extend * max(x$data1), length = 200)
    y <- seq(from = x$u2, to = extend * max(x$data2), length = 200)
    choices <- c("Exceedance data", 
		 "Contours of Bivariate Distribution Function", 
		 "Contours of Bivariate Survival Function",
                 "Tail of Marginal 1", "Tail of Marginal 2")
    tmenu <- paste("plot:", choices)
    pick <- 1
    while(pick > 0) {
        par(mfrow = c(1, 1))
	pick <- menu(tmenu, title =
                     "\nMake a plot selection (or 0 to exit):")
	if(pick == 1) {
	    par(mfrow = c(2, 1))
	    plot(x$data1, main = "Marginal1", type = "n", ...)
	    points((1:length(x$data1))[x$delta1 == 0],
                   x$data1[x$delta1 == 0])
	    points((1:length(x$data1))[x$delta1 == 1],
                   x$data1[x$delta1 == 1], col = 2)
	    plot(x$data2, main = "Marginal2", type = "n", ...)
	    points((1:length(x$data2))[x$delta2 == 0],
                   x$data2[x$delta2 == 0])
	    points((1:length(x$data2))[x$delta2 == 1],
                   x$data2[x$delta2 == 1], col = 2)
	}
	if(pick == 4) {
	    x$name <- "Marginal1"
	    x$par.ests <- x$par.ests1
	    x$data <- x$data1
	    x$threshold <- x$u1
	    x$p.less.thresh <- 1 - x$lambda1
	    tailplot(x, ...)
	}
	if(pick == 5) {
	    x$name <- "Marginal2"
	    x$par.ests <- x$par.ests2
	    x$data <- x$data2
	    x$threshold <- x$u2
	    x$p.less.thresh <- 1 - x$lambda2
	    tailplot(x, ...)
	}
	if(pick == 2) {
	    z <- outer(xx, y, joint, alpha = x$alpha, u1 = x$u1,
                       lambda1 = x$lambda1, xi1 = x$par.ests1[1],
                       sigma1 = x$par.ests1[2], u2 = x$u2, lambda2 =
                       x$lambda2, xi2 = x$par.ests2[1], sigma2 =
                       x$par.ests2[2], Vfunc = x$dep.func)
	    par(xaxs = "i", yaxs = "i")
	    contour(xx, y, z, nlevels = n.contours, main = "Joint", ...)
	}
	if(pick == 3) {
	    z2 <- outer(xx, y, survivor, alpha = x$alpha, u1 = x$u1,
                        lambda1 = x$lambda1, xi1 = x$par.ests1[1],
                        sigma1 = x$par.ests1[2], u2 = x$u2, lambda2 =
                        x$lambda2, xi2 = x$par.ests2[1], sigma2 =
                        x$par.ests2[2], marg = marg, joint = joint,
                        Vfunc = x$dep.func)
	    level.thresh <- x$lambda1 + x$lambda2 - (x$lambda1^(1/x$alpha) +
                x$lambda2^(1/x$alpha))^x$alpha
	    contour(xx, y, z2, nlevels = n.contours, main = "Survival", ...)
	}
    }
}

"emplot" <- 
function(data, alog = "x", labels = TRUE, ...)
{
    data <- sort(as.numeric(data))
    ypoints <- 1 - ppoints(data)
    plot(data, ypoints, log = alog, xlab = "", ylab = "", ...)
    if(labels) {
        xxlab <- "x"
	yylab <- "1 - F(x)"
	if(alog != "")
	    xxlab <- paste(xxlab, "(on log scale)")
	if(alog == "xy" || alog == "yx")
	    yylab <- paste(yylab, "(on log scale)")
	 title(xlab = xxlab, ylab = yylab)
    }
    invisible(list(x = data, y = ypoints))
}

"exindex" <- 
function(data, block, start = 5, end = NA, reverse = FALSE,
    auto.scale = TRUE, labels = TRUE, ...)
{
    sorted <- rev(sort(as.numeric(data)))
    n <- length(sorted)
    if(is.character(block)) {
        times <- as.POSIXlt(attributes(data)$times)
        if(block %in% c("semester", "quarter")) {
            sem <- quart <- times$mon
            sem[sem %in% 0:5] <- quart[quart %in% 0:2] <- 0
            sem[sem %in% 6:11] <- quart[quart %in% 3:5] <- 1
            quart[quart %in% 6:8] <- 2
            quart[quart %in% 9:11] <- 3
        }
        grouping <- switch(block,
            semester = paste(times$year, sem),
            quarter = paste(times$year, quart),
            month = paste(times$year, times$mon),
            year = times$year,
            stop("unknown time period"))
        b.lengths <- as.numeric(tapply(data, grouping, length))
        b.maxima <- as.numeric(tapply(data, grouping, max))
    }
    else {
	data <- as.numeric(data)
	nblocks <- (length(data) %/% block) + 1
	grouping <- rep(1:nblocks, rep(block, nblocks))[1:length(data)]
	b.lengths <- tapply(data, grouping, length)
	b.maxima <- tapply(data, grouping, max)
    }
    b.lengths <- b.lengths[!is.na(b.lengths)]
    b.maxima <- rev(sort(b.maxima[!is.na(b.maxima)]))
    if(is.numeric(block)) r <- block
    else r <- round(mean(b.lengths[2:(length(b.lengths) - 1)]))
    k <- round(n/r)
    un <- unique(b.maxima)[-1]
    K <- match(un, b.maxima) - 1
    N <- match(un, sorted) - 1
    if(is.na(end)) end <- k
    cond <- (K < end) & (K >= start)
    un <- un[cond]
    K <- K[cond]
    N <- N[cond]
    theta2 <- K/N
    theta <- logb(1 - K/k)/(r * logb(1 - N/n))
    out <- cbind(N, K, un, theta2, theta)
    yrange <- range(theta)
    index <- K
    if(reverse)	index <-  - K
    if(auto.scale)
        plot(index, theta, ylim = yrange, type = "l", xlab = "", ylab = "",
             axes = FALSE, ...)
    else plot(index, theta, type = "l", xlab = "", ylab = "", axes =
              FALSE, ...)
    axis(1, at = index, lab = paste(K), tick = FALSE)
    axis(2)
    axis(3, at = index, lab = paste(format(signif(un, 3))), tick = FALSE)
    box()
    if(labels) {
      	ylabel <- paste("theta (", k, " blocks of size ", r, ")", sep = "")
	title(xlab = "K", ylab = ylabel)
	mtext("Threshold", side = 3, line = 3)
    }
    invisible(out)
}

"hill" <- 
function(data, option = c("alpha","xi","quantile"), start = 15, end = NA,
    reverse = FALSE, p = NA, ci = 0.95, auto.scale = TRUE, labels = TRUE, ...)
{
    data <- as.numeric(data)
    ordered <- rev(sort(data))
    ordered <- ordered[ordered > 0]
    n <- length(ordered)
    option <- match.arg(option)
    if((option == "quantile") && (is.na(p)))
        stop("Input a value for the probability p")
    if((option == "quantile") && (p < 1 - start/n)) {
	cat("Graph may look strange !! \n\n")
	cat(paste("Suggestion 1: Increase `p' above",
                  format(signif(1 - start/n, 5)), "\n"))
	cat(paste("Suggestion 2: Increase `start' above ",
                  ceiling(length(data) * (1 - p)), "\n"))
    }
    k <- 1:n
    loggs <- logb(ordered)
    avesumlog <- cumsum(loggs)/(1:n)
    xihat <- c(NA, (avesumlog - loggs)[2:n])
    alphahat <- 1/xihat
    y <- switch(option,
	    alpha = alphahat,
	    xi = xihat,
	    quantile = ordered * ((n * (1 - p))/k)^(-1/alphahat))
    ses <- y/sqrt(k)
    if(is.na(end)) end <- n
    x <- trunc(seq(from = min(end, length(data)), to = start))
    y <- y[x]
    ylabel <- option
    yrange <- range(y)
    if(ci && (option != "quantile")) {
       	qq <- qnorm(1 - (1 - ci)/2)
       	u <- y + ses[x] * qq
       	l <- y - ses[x] * qq
       	ylabel <- paste(ylabel, " (CI, p =", ci, ")", sep = "")
       	yrange <- range(u, l)
    }
    if(option == "quantile") ylabel <- paste("Quantile, p =", p)
    index <- x
    if(reverse) index <-  - x
    if(auto.scale)
        plot(index, y, ylim = yrange, type = "l", xlab = "", ylab = "",
	     axes = FALSE, ...)
    else plot(index, y, type = "l", xlab = "", ylab = "", axes = FALSE, ...)
    axis(1, at = index, lab = paste(x), tick = FALSE)
    axis(2)
    threshold <- findthresh(data, x)
    axis(3, at = index, lab = paste(format(signif(threshold, 3))),
         tick = FALSE)
    box()
    if(ci && (option != "quantile")) {
       	lines(index, u, lty = 2, col = 2)
       	lines(index, l, lty = 2, col = 2)
    }
    if(labels) {
       	title(xlab = "Order Statistics", ylab = ylabel)
       	mtext("Threshold", side = 3, line = 3)
    }
    invisible(list(x = index, y = y))
}

"meplot" <- 
function(data, omit = 3, labels = TRUE, ...)
{
    data <- as.numeric(data)
    n <- length(data)
    myrank <- function(x, na.last = TRUE)
    {
        ranks <- sort.list(sort.list(x, na.last = na.last))
	if(is.na(na.last))
	     x <- x[!is.na(x)]
	for(i in unique(x[duplicated(x)])) {
	    which <- x == i & !is.na(x)
	    ranks[which] <- max(ranks[which])
	}
	ranks
    }
    data <- sort(data)
    n.excess <- unique(floor(length(data) - myrank(data)))
    points <- unique(data)
    nl <- length(points)
    n.excess <- n.excess[-nl]
    points <- points[-nl]
    excess <- cumsum(rev(data))[n.excess] - n.excess * points
    y <- excess/n.excess
    xx <- points[1:(nl-omit)] ; yy <- y[1:(nl-omit)]
    plot(xx, yy, xlab = "", ylab = "", ...)
    if(labels) title(xlab = "Threshold", ylab = "Mean Excess")
    invisible(list(x = xx, y = yy))
}

"qplot" <- 
function(data, xi = 0, trim = NA, threshold = NA, line = TRUE,
    labels = TRUE, ...)
{
    data <- as.numeric(data)
    if(!is.na(threshold)) data <- data[data >= threshold]
    if(!is.na(trim)) data <- data[data < trim]
    if(xi == 0) {
        add <- "Exponential Quantiles"
	y <- qexp(ppoints(data))
    }
    if(xi != 0) {
        add <- paste("GPD Quantiles; xi =", xi)
	y <- qgpd(ppoints(data), xi = xi)
    }
    plot(sort(data), y, xlab = "", ylab = "", ...)
    if(labels) title(xlab = "Ordered Data", ylab = add)
    if(line) abline(lsfit(sort(data), y))
    invisible(list(x = sort(data), y = y))
}

"records" <- 
function(data, do.plot = TRUE, conf.level = 0.95, ...)
{
    data <- as.numeric(data)
    record <- cummax(data)
    expected <- cumsum(1/(1:length(data)))
    se <- sqrt(expected - cumsum(1/((1:length(data))^2)))
    trial <- (1:length(data))[!duplicated(record)]
    record <- unique(record)
    number <- 1:length(record)
    expected <- expected[trial]
    se <- se[trial]
    if(do.plot) {
       	ci <- qnorm(0.5 + conf.level/2)
       	upper <- expected + ci * se
       	lower <- expected - ci * se
       	lower[lower < 1] <- 1
       	yr <- range(upper, lower, number)
       	plot(trial, number, log = "x", ylim = yr, xlab = "Trial",
             ylab = "Records", main = "Plot of Record Development", ...)
	lines(trial, expected)
	lines(trial, upper, lty = 2)
	lines(trial, lower, lty = 2)
    }
    data.frame(number, record, trial, expected, se)
}

"gpd" <- 
function(data, threshold = NA, nextremes = NA, method = c("ml","pwm"),
         information = c("observed","expected"), ...)
{
    data <- as.numeric(data)
    n <- length(data)
    if(is.na(nextremes) && is.na(threshold))
        stop("Enter either a threshold or the number of upper extremes")
    if(!is.na(nextremes) && !is.na(threshold))
        stop("Enter EITHER a threshold or the number of upper extremes")
    if(!is.na(nextremes))
        threshold <- findthresh(data, nextremes)
    exceedances <- data[data > threshold]
    excess <- exceedances - threshold
    Nu <- length(excess)
    xbar <- mean(excess)
    method <- match.arg(method)
    if(method == "ml") {
        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)
        if(fit$convergence)
            warning("optimization may not have 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)
	}
    }
    if(method == "pwm") {
        a0 <- xbar
        gamma <- -0.35
        delta <- 0
        pvec <- ((1:Nu) + delta)/(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)
        denom <- Nu * (1 - 2 * xi) * (3 - 2 * xi)
        if(xi > 0.5) {
            denom <- NA
      	    warning("Asymptotic standard errors not available for",
                    "PWM Method 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))
    p.less.thresh <- 1 - Nu/n
    out <- list(n = length(data), data = exceedances, threshold =
        threshold, p.less.thresh = p.less.thresh, n.exceed = Nu,
        method = method, par.ests = par.ests, par.ses = par.ses,
        varcov = varcov, information = information, converged =
        converged, nllh.final = nllh.final)
    names(out$par.ests) <- c("xi", "beta")
    names(out$par.ses) <- c("xi", "beta")
    class(out) <- "gpd"
    out
}

"gpd.q" <- 
function(x, pp, ci.type = c("likelihood","wald"), ci.p = 0.95,
         like.num = 50)
{
    if(x$dist != "gpd")
        stop("This function is used only with GPD curves")
    if(length(pp) > 1)
	stop("One probability at a time please")
    threshold <- x$lastfit$threshold
    par.ests <- x$lastfit$par.ests
    xihat <- par.ests["xi"]
    betahat <- par.ests["beta"]
    varcov <- x$lastfit$varcov
    p.less.thresh <- x$lastfit$p.less.thresh
    lambda <- 1
    if(x$type == "tail") lambda <- 1/(1 - p.less.thresh)
    a <- lambda * (1 - pp)
    gfunc <- function(a, xihat) (a^( - xihat) - 1) / xihat
    gfunc.deriv <- function(a, xihat)
        ( - (a^( - xihat) - 1)/xihat - a^( - xihat) * logb(a)) / xihat
    q <- threshold + betahat * gfunc(a, xihat)
    if(q < x$plotmax) abline(v = q, lty = 2)
    out <- as.numeric(q)
    ci.type <- match.arg(ci.type)
    if(ci.type == "wald") {
        if(class(x$lastfit) != "gpd")
	    stop("Wald method requires model be fitted with gpd (not pot)")
	scaling <- threshold
	betahat <- betahat/scaling
	xivar <- varcov[1, 1]
        betavar <- varcov[2, 2]/(scaling^2)
	covar <- varcov[1, 2]/scaling
	term1 <- betavar * (gfunc(a, xihat))^2
	term2 <- xivar * (betahat^2) * (gfunc.deriv(a, xihat))^2
	term3 <- 2 * covar * betavar * gfunc(a, xihat) * gfunc.deriv(a, xihat)
	qvar <- term1 + term2 + term3
	if(qvar < 0) stop("Negative estimate of quantile variance")
	qse <- scaling * sqrt(qvar)
	qq <- qnorm(1 - (1 - ci.p)/2)
	upper <- q + qse * qq
	lower <- q - qse * qq
        if(upper < x$plotmax) abline(v = upper, lty = 2, col = 2)
	if(lower < x$plotmax) abline(v = lower, lty = 2, col = 2)
	out <- as.numeric(c(lower, q, qse, upper))
	names(out) <- c("Lower CI", "Estimate", "Std.Err", "Upper CI")
    }
    if(ci.type == "likelihood") {
        parloglik <- function(theta, tmp, a, threshold, xpi)
	{
	    beta <- (theta * (xpi - threshold))/(a^( - theta) - 1)
	    if((beta <= 0) || ((theta <= 0) && (max(tmp) > ( - beta/theta))))
	        f <- 1e+06
	    else {
	        y <- logb(1 + (theta * tmp)/beta)
		y <- y/theta
		f <- length(tmp) * logb(beta) + (1 + theta) * sum(y)
	    }
	    f
	}
	theta <- xihat
	parmax <- NULL
	xp <- exp(seq(from = logb(threshold), to = logb(x$plotmax),
                      length = like.num))
        excess <- as.numeric(x$lastfit$data - threshold)
	for(i in 1:length(xp)) {
            fit2 <- optim(theta, parloglik, method = "BFGS", hessian = FALSE,
                tmp = excess, a = a, threshold = threshold, xpi = xp[i])
	    parmax <- rbind(parmax, fit2$value)
	}
	parmax <-  - parmax
	overallmax <-  - parloglik(xihat, excess, a, threshold, q)
	crit <- overallmax - qchisq(0.999, 1)/2
	cond <- parmax > crit
	xp <- xp[cond]
	parmax <- parmax[cond]
	par(new = T)
	dolog <- ""
        if(x$alog == "xy" || x$alog == "x") dolog <- "x"
	plot(xp, parmax, type = "n", xlab = "", ylab = "", axes = F,
	     xlim = range(x$plotmin, x$plotmax),
	     ylim = range(overallmax, crit), log = dolog)
	axis(4, at = overallmax - qchisq(c(0.95, 0.99), 1)/2,
             labels = c("95", "99"), tick = TRUE)
	aalpha <- qchisq(ci.p, 1)
	abline(h = overallmax - aalpha/2, lty = 2, col = 2)
	cond <- !is.na(xp) & !is.na(parmax)
	smth <- spline(xp[cond], parmax[cond], n = 200)
	lines(smth, lty = 2, col = 2)
	ci <- smth$x[smth$y > overallmax - aalpha/2]
	out <- c(min(ci), q, max(ci))
	names(out) <- c("Lower CI", "Estimate", "Upper CI")
    }
    out
}

"gpd.sfall" <- 
function(x, pp, ci.p = 0.95, like.num = 50)
{
    if(x$dist != "gpd")
       	stop("This function is used only with GPD curves")
    if(length(pp) > 1)
       	stop("One probability at a time please")
    threshold <- x$lastfit$threshold
    par.ests <- x$lastfit$par.ests
    xihat <- par.ests["xi"]
    betahat <- par.ests["beta"]
    varcov <- x$lastfit$varcov
    p.less.thresh <- x$lastfit$p.less.thresh
    lambda <- 1
    if(x$type == "tail") lambda <- 1/(1 - p.less.thresh)
    a <- lambda * (1 - pp)
    gfunc <- function(a, xihat) (a^( - xihat) - 1) / xihat
    q <- threshold + betahat * gfunc(a, xihat)
    s <- q + (betahat + xihat * (q - threshold))/(1 - xihat)
    if(s < x$plotmax) abline(v = s, lty = 2)
    out <- as.numeric(s)
    parloglik <- function(theta, tmp, a, threshold, xpi)
    {
        beta <- ((1 - theta) * (xpi - threshold)) /
          (((a^( - theta) - 1)/theta) + 1)
	if((beta <= 0) || ((theta <= 0) && (max(tmp) > ( - beta/theta))))
	    f <- 1e+06
	else {
	    y <- logb(1 + (theta * tmp)/beta)
	    y <- y/theta
	    f <- length(tmp) * logb(beta) + (1 + theta) * sum(y)
	}
	f
    }
    theta <- xihat
    parmax <- NULL
    xp <- exp(seq(from = logb(threshold), to = logb(x$plotmax),
                  length = like.num))
    excess <- as.numeric(x$lastfit$data - threshold)
    for(i in 1:length(xp)) {
        fit2 <- optim(theta, parloglik, method = "BFGS", hessian = FALSE,
                tmp = excess, a = a, threshold = threshold, xpi = xp[i])
        parmax <- rbind(parmax, fit2$value)
    }
    parmax <-  - parmax
    overallmax <-  - parloglik(xihat, excess, a, threshold, s)
    crit <- overallmax - qchisq(0.999, 1)/2
    cond <- parmax > crit
    xp <- xp[cond]
    parmax <- parmax[cond]
    par(new = T)
    dolog <- ""
    if(x$alog == "xy" || x$alog == "x") dolog <- "x"
    plot(xp, parmax, type = "n", xlab = "", ylab = "", axes = F, xlim = 
	 range(x$plotmin, x$plotmax), ylim =
         range(overallmax, crit), log = dolog)
    axis(4, at = overallmax - qchisq(c(0.95, 0.99), 1)/2,
         labels = c("95", "99"), tick = TRUE)
    aalpha <- qchisq(ci.p, 1)
    abline(h = overallmax - aalpha/2, lty = 2, col = 2)
    cond <- !is.na(xp) & !is.na(parmax)
    smth <- spline(xp[cond], parmax[cond], n = 200)
    lines(smth, lty = 2, col = 2)
    ci <- smth$x[smth$y > overallmax - aalpha/2]
    out <- c(min(ci), s, max(ci))
    names(out) <- c("Lower CI", "Estimate", "Upper CI")
    out
}

"plot.gpd" <- 
function(x, optlog = NA, extend = 1.5, labels = TRUE, ...)
{
    data <- as.numeric(x$data)
    threshold <- x$threshold
    xi <- x$par.ests["xi"]
    beta <- x$par.ests["beta"]
    choices <- c("Excess Distribution", "Tail of Underlying Distribution",
      	"Scatterplot of Residuals", "QQplot of Residuals")
    tmenu <- paste("plot:", choices)
    pick <- 1
    lastcurve <- NULL
    while(pick > 0) {
        pick <- menu(tmenu, title =
                     "\nMake a plot selection (or 0 to exit):")
        if(pick >= 3) {
            excess <- data - threshold
       	    res <- logb(1 + (xi * excess)/beta) / xi
            lastcurve <- NULL
	}
        if(pick == 3) {
      	    plot(res, ylab = "Residuals", xlab = "Ordering", ...)
	    lines(lowess(1:length(res), res))
	}
        if(pick == 4) qplot(res, ...)
        if(pick == 1 || pick == 2) {
            plotmin <- threshold
     	    if(extend <= 1) stop("extend must be > 1")
	    plotmax <- max(data) * extend
            xx <- seq(from = 0, to = 1, length = 1000)
      	    z <- qgpd(xx, xi, threshold, beta)
       	    z <- pmax(pmin(z, plotmax), plotmin)
       	    ypoints <- ppoints(sort(data))
       	    y <- pgpd(z, xi, threshold, beta)
	}
	if(pick == 1) {
	    type <- "eplot"
	    if(!is.na(optlog))
                alog <- optlog
       	    else alog <- "x"
	    if(alog == "xy")
	        stop("Double log plot of Fu(x-u) does\nnot make much sense")
	    yylab <- "Fu(x-u)"
       	    shape <- xi
	    scale <- beta
	    location <- threshold
	}
	if(pick == 2) {
	    type <- "tail"
	    if(!is.na(optlog))
	        alog <- optlog
	    else alog <- "xy"
	    prob <- x$p.less.thresh
	    ypoints <- (1 - prob) * (1 - ypoints)
	    y <- (1 - prob) * (1 - y)
	    yylab <- "1-F(x)"
	    shape <- xi
	    scale <- beta * (1 - prob)^xi
	    location <- threshold - (scale * ((1 - prob)^( - xi) - 1))/xi
	}
	if(pick == 1 || pick == 2) {
	    plot(sort(data), ypoints, xlim = range(plotmin, plotmax),
                 ylim = range(ypoints, y, na.rm = TRUE), xlab = "",
                 ylab = "", log = alog, axes = TRUE, ...)
	    lines(z[y >= 0], y[y >= 0])
	    if(labels) {
	        xxlab <- "x"
		if(alog == "x" || alog == "xy" || alog == "yx")
		    xxlab <- paste(xxlab, "(on log scale)")
	        if(alog == "xy" || alog == "yx" || alog == "y")
		    yylab <- paste(yylab, "(on log scale)")
		title(xlab = xxlab, ylab = yylab)
	    }
	    details <- paste("threshold = ", format(signif(threshold, 3)),
                             "   xi = ", format(signif(shape, 3)),
                             "   scale = ", format(signif(scale, 3)),
                             "   location = ", format(signif(location, 3)),
                             sep = "")
	    print(details)
	    lastcurve <- list(lastfit = x, type = type, dist = "gpd",
                plotmin = plotmin, plotmax = plotmax, alog = alog,
                location = as.numeric(location), shape = as.numeric(shape),
                scale = as.numeric(scale))
	}
    }
    invisible(lastcurve)
}

"quant" <- 
function(data, p = 0.99, models = 30, start = 15, end = 500,
	reverse = TRUE, ci = 0.95, auto.scale = TRUE, labels = TRUE,
	...)
{
    data <- as.numeric(data)
    n <- length(data)
    if(ci) qq <- qnorm(1 - (1 - ci)/2)
    exceed <- trunc(seq(from = min(end, n), to = start, length = models))
    if(p < 1 - min(exceed)/n) {
        cat("Graph may look strange !! \n\n")
	cat(paste("Suggestion 1: Increase `p' above",
            format(signif(1 - min(exceed)/n, 5)), "\n"))
	cat(paste("Suggestion 2: Increase `start' above ",
            ceiling(length(data) * (1 - p)), "\n"))
    }
    gpd.dummy <- function(nex, data)
    {
	out <- gpd(data = data, nex = nex, information = "expected")
	c(out$threshold, out$par.ests[1], out$par.ests[2],
          out$varcov[1, 1], out$varcov[2, 2], out$varcov[1, 2])
    }
    mat <- apply(as.matrix(exceed), 1, gpd.dummy, data = data)
    thresh <- mat[1,  ]
    xihat <- mat[2,  ]
    betahat <- mat[3,  ]
    lambda <- length(data)/exceed
    a <- lambda * (1 - p)
    gfunc <- function(a, xihat) (a^( - xihat) - 1) / xihat
    qest <- thresh + betahat * gfunc(a, xihat)
    l <- u <- qest
    yrange <- range(qest)
    if(ci) {
        xivar <- mat[4,  ]
	betavar <- mat[5,  ]
	covar <- mat[6,  ]
	scaling <- thresh
	betahat <- betahat/scaling
	betavar <- betavar/(scaling^2)
	covar <- covar/scaling
	gfunc.deriv <- function(a, xihat)
	    ( - (a^( - xihat) - 1)/xihat - a^( - xihat) * logb(a)) / xihat
	term1 <- betavar * (gfunc(a, xihat))^2
	term2 <- xivar * (betahat^2) * (gfunc.deriv(a, xihat))^2
	term3 <- 2 * covar * betavar * gfunc(a, xihat) * gfunc.deriv(a, xihat)
	qvar <- term1 + term2 + term3
	if(min(qvar) < 0)
	    stop(paste("Conditioning problems lead to estimated negative",
                       "quantile variance", sep = "\n"))
	qse <- scaling * sqrt(qvar)
	u <- qest + qse * qq
	l <- qest - qse * qq
	yrange <- range(qest, u, l)
    }
    mat <- rbind(thresh, qest, exceed, l, u)
    dimnames(mat) <- list(c("threshold", "qest", "exceedances", "lower",
        "upper"), NULL)
    index <- exceed
    if(reverse) index <-  - exceed
    if(auto.scale)
        plot(index, qest, ylim = yrange, type = "l", xlab = "", ylab = "",
             axes = FALSE, ...)
    else plot(index, qest, type = "l", xlab = "", ylab = "",
              axes = FALSE, ...)
    axis(1, at = index, lab = paste(exceed))
    axis(2)
    axis(3, at = index, lab = paste(format(signif(thresh, 3))))
    box()
    if(ci) {
       	lines(index, l, lty = 2, col = 2)
       	lines(index, u, lty = 2, col = 2)
    }
    if(labels) {
       	labely <- paste(p, "Quantile")
       	if(ci) labely <- paste(labely, " (CI, p = ", ci, ")", sep = "")
	title(xlab = "Exceedances", ylab = labely)
	mtext("Threshold", side = 3, line = 3)
    }
    invisible(mat)
}

"riskmeasures" <- 
function(x, p)
{
    u <- x$threshold
    par.ests <- x$par.ests
    xihat <- par.ests["xi"]
    betahat <- par.ests["beta"]
    p.less.thresh <- x$p.less.thresh
    lambda <- 1/(1 - p.less.thresh)
    quant <- function(pp, xi, beta, u, lambda)
    {
     	a <- lambda * (1 - pp)
       	u + (beta * (a^( - xi) - 1))/xi
    }
    short <- function(pp, xi, beta, u, lambda)
    {
      	a <- lambda * (1 - pp)
       	q <- u + (beta * (a^( - xi) - 1))/xi
       	(q * (1 + (beta - xi * u)/q)) / (1 - xi)
    }
    q <- quant(p, xihat, betahat, u, lambda)
    es <- short(p, xihat, betahat, u, lambda)
    rtn <- cbind(p, quantile = q, sfall = es)
    row.names(rtn) <- NULL
    rtn
}

"shape" <- 
function(data, models = 30, start = 15, end = 500, reverse = TRUE, ci = 
	0.95, auto.scale = TRUE, labels = TRUE, ...)
{
    data <- as.numeric(data)
    qq <- 0
    if(ci) qq <- qnorm(1 - (1 - ci)/2)
    x <- trunc(seq(from = min(end, length(data)), to = start, length = models))
    gpd.dummy <- function(nex, data)
    {
        out <- gpd(data = data, nex = nex, information = "expected")
	c(out$threshold, out$par.ests[1], out$par.ses[1])
    }
    mat <- apply(as.matrix(x), 1, gpd.dummy, data = data)
    mat <- rbind(mat, x)
    dimnames(mat) <- list(c("threshold", "shape", "se", "exceedances"), NULL)
    thresh <- mat[1,  ]
    y <- mat[2,  ]
    yrange <- range(y)
    if(ci) {
        u <- y + mat[3,  ] * qq
	l <- y - mat[3,  ] * qq
	yrange <- range(y, u, l)
    }
    index <- x
    if(reverse) index <-  - x
    if(auto.scale)
        plot(index, y, ylim = yrange, type = "l", xlab = "", ylab = "",
	     axes = FALSE, ...)
    else plot(index, y, type = "l", xlab = "", ylab = "", axes = FALSE, ...)
    axis(1, at = index, lab = paste(x), tick = FALSE)
    axis(2)
    axis(3, at = index, lab = paste(format(signif(thresh, 3))), tick = FALSE)
    box()
    if(ci) {
        lines(index, u, lty = 2, col = 2)
	lines(index, l, lty = 2, col = 2)
    }
    if(labels) {
        labely <- "Shape (xi)"
	if(ci) labely <- paste(labely, " (CI, p = ", ci, ")", sep = "")
	title(xlab = "Exceedances", ylab = labely)
	mtext("Threshold", side = 3, line = 3)
    }
    invisible(mat)
}

"tailplot" <- 
function(x, optlog = NA, extend = 1.5, labels = TRUE, ...)
{
    data <- as.numeric(x$data)
    threshold <- x$threshold
    xi <- x$par.ests["xi"]
    beta <- x$par.ests["beta"]
    plotmin <- threshold
    if(extend <= 1) stop("extend must be > 1")
    plotmax <- max(data) * extend
    xx <- seq(from = 0, to = 1, length = 1000)
    z <- qgpd(xx, xi, threshold, beta)
    z <- pmax(pmin(z, plotmax), plotmin)
    ypoints <- ppoints(sort(data))
    y <- pgpd(z, xi, threshold, beta)
    type <- "tail"
    if(!is.na(optlog))
	    alog <- optlog
    else alog <- "xy"
    prob <- x$p.less.thresh
    ypoints <- (1 - prob) * (1 - ypoints)
    y <- (1 - prob) * (1 - y)
    yylab <- "1-F(x)"
    shape <- xi
    scale <- beta * (1 - prob)^xi
    location <- threshold - (scale * ((1 - prob)^( - xi) - 1))/xi
    plot(sort(data), ypoints, xlim = range(plotmin, plotmax), ylim =
         range(ypoints, y, na.rm = TRUE), xlab = "", ylab = "", log = alog,
         axes = TRUE, ...)
    lines(z[y >= 0], y[y >= 0])
    if(labels) {
        xxlab <- "x"
        if(alog == "x" || alog == "xy" || alog == "yx")
	    xxlab <- paste(xxlab, "(on log scale)")
        if(alog == "xy" || alog == "yx" || alog == "y")
	    yylab <- paste(yylab, "(on log scale)")
        title(xlab = xxlab, ylab = yylab)
    }
    lastcurve <- list(lastfit = x, type = type, dist = "gpd",
        plotmin = plotmin, plotmax = plotmax, alog = alog, location = 
	as.numeric(location), shape = as.numeric(shape), scale = 
	as.numeric(scale))
    invisible(lastcurve)
}
"plot.pot" <- 
function(x, ...)
{
    rawdata <- x$data
    n <- length(as.numeric(rawdata))
    times <- attributes(rawdata)$times
    if(is.character(times) || inherits(times, "POSIXt") ||
       inherits(x, "date") || inherits(x, "dates")) {
        times <- as.POSIXlt(times)
        gaps <- as.numeric(difftime(times[2:n], times[1:(n-1)],
            units = "days")) * x$intensity
    }
    else gaps <- as.numeric(diff(times)) * x$intensity
    data <- as.numeric(rawdata)
    threshold <- x$threshold
    par.ests <- x$par.ests
    xi <- par.ests[1]
    beta <- par.ests[4]
    residuals <- logb(1 + (xi * (data - threshold))/beta)/xi
    choices <- c("Point Process of Exceedances", "Scatterplot of Gaps",
		 "Qplot of Gaps", "ACF of Gaps", "Scatterplot of Residuals",
		 "Qplot of Residuals", "ACF of Residuals", "Go to GPD Plots")
    tmenu <- paste("plot:", choices)
    pick <- 1
    lastcurve <- NULL
    while(pick > 0) {
        pick <- menu(tmenu, title = 
		     "\nMake a plot selection (or 0 to exit):")
        if(pick %in% c(4,7)) require("ts", quietly = TRUE)
        if(pick %in% 1:7) lastcurve <- NULL
	switch(pick,
            {
	       plot(times, rawdata, type = "h", sub = paste("Point process of",
	         length(as.numeric(rawdata)), "exceedances of threshold",
                 format(signif(threshold, 3))), ...)  
            },
	    {
	      plot(gaps, ylab = "Gaps", xlab = "Ordering", ...)
	      lines(lowess(1:length(gaps), gaps))
	    },
	      qplot(gaps, ...),
	      acf(gaps, lag.max = 20, ...),
	    {
	      plot(residuals, ylab = "Residuals", xlab = "Ordering", ...)
	      lines(lowess(1:length(residuals), residuals))
	    },
	      qplot(residuals, ...),
	      acf(residuals, lag.max = 20, ...),
	      lastcurve <- plot.gpd(x, ...))
    }
    invisible(lastcurve)
}

"pot" <- 
function(data, threshold = NA, nextremes = NA, run = NA,
    picture = TRUE, ...)
{
    n <- length(as.numeric(data))
    times <- attributes(data)$times
    if(is.null(times)) {
        times <- 1:n
        attributes(data)$times <- times
        start <- 1
        end <- n
        span <- end - start
    }
    else {
        start <- times[1]
        end <- times[n]
        span <- as.numeric(difftime(as.POSIXlt(times)[n],
            as.POSIXlt(times)[1], units = "days"))
    }
  
    if(is.na(nextremes) && is.na(threshold))
       	stop("Enter either a threshold or the number of upper extremes")
    if(!is.na(nextremes) && !is.na(threshold))
	stop("Enter EITHER a threshold or the number of upper extremes")
    if(!is.na(nextremes))
	threshold <- findthresh(as.numeric(data), nextremes)
    if(threshold > 10) {
	factor <- 10^(floor(log10(threshold)))
	cat(paste("If singularity problems occur divide data",
                  "by a factor, perhaps", factor, "\n"))
    }
    exceedances.its <- structure(data[data > threshold], times =
        times[data > threshold])
    n.exceed <- length(as.numeric(exceedances.its))
    p.less.thresh <- 1 - n.exceed/n
    if(!is.na(run)) {
       	exceedances.its <- decluster(exceedances.its, run, picture)
       	n.exceed <- length(exceedances.its)
    }
    intensity <- n.exceed/span
    exceedances <- as.numeric(exceedances.its)
    xbar <- mean(exceedances) - threshold
    s2 <- var(exceedances)
    shape0 <- -0.5 * (((xbar * xbar)/s2) - 1)
    extra <- ((length(exceedances)/span)^( - shape0) - 1)/shape0
    betahat <- 0.5 * xbar * (((xbar * xbar)/s2) + 1)
    scale0 <- betahat/(1 + shape0 * extra)
    loc0 <- 0
    theta <- c(shape0, scale0, loc0)
    negloglik <- function(theta, exceedances, threshold, span)
    {
        if((theta[2] <= 0) || (min(1 + (theta[1] * (exceedances -
            theta[3])) / theta[2]) <= 0))
	    f <- 1e+06
	else {
	    y <- logb(1 + (theta[1] * (exceedances - theta[3])) / theta[2])
	    term3 <- (1/theta[1] + 1) * sum(y)
	    term1 <- span * (1 + (theta[1] * (threshold - theta[3])) /
                         theta[2])^(-1/theta[1])
	    term2 <- length(y) * logb(theta[2])
	    f <- term1 + term2 + term3
	}
	f
    }
    fit <- optim(theta, negloglik, hessian = TRUE, ..., exceedances =
                 exceedances, threshold = threshold, span = span)
    if(fit$convergence)
        warning("optimization may not have succeeded")
    par.ests <- fit$par
    varcov <- solve(fit$hessian)
    par.ses <- sqrt(diag(varcov))   
    beta <- par.ests[2] + par.ests[1] * (threshold - par.ests[3])
    par.ests <- c(par.ests, beta)
    out <- list(n = length(data), period = c(start, end), data = 
        exceedances.its, span = span, threshold = threshold,
        p.less.thresh = p.less.thresh, n.exceed = n.exceed, run = run,
	par.ests = par.ests, par.ses = par.ses, varcov = varcov, 
	intensity = intensity, nllh.final = fit$value, converged
	= fit$convergence)
    names(out$par.ests) <- c("xi", "sigma", "mu", "beta")
    names(out$par.ses) <- c("xi", "sigma", "mu")
    class(out) <- "pot"
    out
}

"decluster" <- 
function(series, run = NA, picture = TRUE)
{
    n <- length(as.numeric(series))
    times <- attributes(series)$times
    if(is.null(times)) stop("`series' must have a `times' attribute")
    as.posix <- is.character(times) || inherits(times, "POSIXt") ||
      inherits(times, "date") || inherits(times, "dates")
    if(as.posix)
        gaps <- as.numeric(difftime(as.POSIXlt(times)[2:n],
            as.POSIXlt(times)[1:(n-1)], units = "days"))
    else gaps <- as.numeric(diff(times))
    longgaps <- gaps > run
    if(sum(longgaps) <= 1)
        stop("Decluster parameter too large")
    cluster <- c(0, cumsum(longgaps))
    cmax <- tapply(as.numeric(series), cluster, max)
    newtimes <- times[match(cmax, series)]
    newseries <- structure(series[match(cmax, series)], times = newtimes)
    n <- length(as.numeric(newseries))

    if(as.posix) {
        newgaps <- as.numeric(difftime(as.POSIXlt(newtimes)[2:n],
            as.POSIXlt(newtimes)[1:(n-1)], units = "days"))
        times <- as.POSIXlt(times)
        newtimes <- as.POSIXlt(newtimes)
    }
    else newgaps <- as.numeric(diff(newtimes))
    
    if(picture) {
      	cat("Declustering picture...\n")
       	cat(paste("Data reduced from", length(as.numeric(series)),
       		"to", length(as.numeric(newseries)), "\n"))
       	par(mfrow = c(2, 2))
        plot(times, series, type = "h")
	qplot(gaps)
        plot(newtimes, newseries, type = "h")
       	qplot(newgaps)
       	par(mfrow = c(1, 1))
    }
    newseries
}

"findthresh" <- 
function(data, ne)
{
	data <- rev(sort(as.numeric(data)))
	thresholds <- unique(data)
	indices <- match(data[ne], thresholds)
	indices <- pmin(indices + 1, length(thresholds))
	thresholds[indices]
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION			DESCRIPTION:
#  emdPlot			 Plots empirical distribution function
#  qqPlot			 Creates a quantile-quantile plot
#  qqbayesPlot		 Creates qq-Plot with 95 percent intervals
#  qPlot			 Creates exploratory QQ plot for EV analysis
#  mePlot			 Creates a sample mean excess plot
#   mxfPlot	          Plots the mean excess function
#   mrlPlot		      Returns a mean residual life plot with confidence levels
#  recordsPlot		 Plots records development
#   ssrecordsPlot	  Plots records development of data subsamples
#  msratioPlot		 Plots ratio of maximums and sums
#  xacfPlot			 Plots autocorrelations of exceedences
#  interactivePlot   Plots several graphs interactively
#  gridVector        Creates from two vectors rectangular grid points
################################################################################


emdPlot =
function(x, doplot = TRUE, plottype = c("", "x", "y", "xy"),
labels = TRUE, ...)
{	# A function imported from R-package evir

	# Description:
	#	Plots empirical distribution function
	
	# Arguments:
	#	plottype - which axes should be on a log scale: 
	#		"x" denotes x-axis only; "y" denotes y-axis only,
	#	    "xy" || "yx" both axes, "" denotes neither of the 
	#		axis

	# FUNCTION:
	
	# Settings:
	alog = plottype[1]
		
	# Convert x to a vector, if the input is a data.frame.
	if (is.data.frame(x)) x = x[,1] 
	xs = x = sort(as.numeric(x))
	ys = y = 1 - ppoints(x)
	
	if (plottype == "x") {
		xs = x[x > 0]
		ys = y[x > 0] }
	if (plottype == "y") {
		xs = x[y > 0]
		ys = y[y > 0] }
	if (plottype == "xy") {
		xs = x[x > 0 & y > 0]
		ys = y[x > 0 & y > 0] }
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "x"
		    ylab = "1-F(x)"
		    main = "Empirical Distribution" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		plot(xs, ys, log = alog, xlab = xlab, ylab = ylab, main = main, ...) 	
		if (labels) grid()
	}			
	
	# Result:
	result = data.frame(x, y)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)	
}


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


qqPlot =
function (x, doplot = TRUE, labels = TRUE, ...) 
{	# A function written by Diethelm Wuertz
	
	# Description:
	#	Creates Quantile-Quantile Plot
	
	# FUNCTION:
	
	# Convert x to a vector, if the input is a data.frame.
	if (is.data.frame(x)) x = x[,1] 
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Normal Quantiles"
		    ylab = "Empirical Quantiles"
		    main = "Normal QQ-Plot" 
		    print(main) }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		qqnorm(x, xlab = xlab, ylab = ylab, main = main, ...) 
		qqline(x) 
		if (labels) grid() 
	}
	
	# Return Value:
	if (doplot) return(invisible(x)) else return(x)
}


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


qqbayesPlot = 
function(x, doplot = TRUE, labels = TRUE, ...) 
{ 	# A function implemented by Diethelm Wuertz

	# Description:
	#	Example of a Normal quantile plot of data x to provide a visual
	# 	assessment of its conformity with a normal (data is standardised 	
	#	first).

	# Details:
	#	The ordered data values are posterior point estimates of the 
	#	underlying quantile function. So, if you plot the ordered data 
	#	values (y-axis) against the exact theoretical quantiles (x-axis), 	
	#	you get a scatter that should be close to a straight line if the 
	#	data look like a random sample from the theoretical distribution. 
	#	This function chooses the normal as the theory, to provide a 
	#	graphical/visual assessment of how normal the data appear.
	#	To help with assessing the relevance of sampling variability on 
	#	just "how close" to the normal the data appears, we add (very) 
	#	approximate posterior 95% intervals for the uncertain quantile 
	#	function at each point (Based on approximate theory) .

	# Author:
	#	Prof. Mike West, mw@stat.duke.edu 
	
	# Note:
	#	Source from
	#	http://www.stat.duke.edu/courses/Fall99/sta290/Notes/

	# FUNCTION:
	
	# Settings:
	mydata = x
   	n = length(mydata) 
   	p = (1:n)/(n+1)
   	x = (mydata-mean(mydata))/sqrt(var(mydata))
   	x = sort(x)
   	z = qnorm(p)
 
   	# Plot:
   	if (doplot) {
   		if (labels) {
		    xlab = "Standard Normal Quantiles"
   		    ylab = "Ordered Data"
   		    main = "Normal QQ-Plot with 95% Intervals" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
   		plot(z, x, xlab = xlab, ylab = ylab, main = main, ...)
   		abline(0, 1, col = "steelblue3")
   		if (labels) grid() 
   	}
  
   	# 95% Intervals:
   	s = 1.96*sqrt(p*(1-p)/n)
   	pl = p-s; i = pl<1&pl>0
   	lower = quantile(x, probs = pl[i])
   	if (doplot) lines(z[i], lower, col = 3)
   	pl = p+s; i = pl < 1 & pl > 0
   	upper = quantile(x, probs = pl[i])
   	if (doplot) lines(z[i], upper, col = 3)
   	
   	# Result:
	result = data.frame(lower, upper)
	
	# Return Value:
   	if (doplot) return(invisible(result)) else return(result)
}     


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


qPlot =
function(x, xi = 0, trim = NA, threshold = NA, doplot = TRUE, 
labels = TRUE, ...)
{	# A function imported from R-package evir
	
	# Description:
	#	Creates an exploratory QQplot for Extreme Value Analysis.

	# FUNCTION:
	
	# Settings:
	line = TRUE
	
	# Convert x to a vector, if the input is a data.frame.
	if(is.data.frame(x)) x = x[,1] 
	x = as.numeric(x)
	if (!is.na(threshold)) x = x[x >= threshold]
	if (!is.na(trim)) x = x[x < trim]
	if (xi == 0) {
		y = qexp(ppoints(x)) }
	if( xi != 0) {
		y = qgpd(ppoints(x), xi = xi) }
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Ordered Data"
			ylab = "Quantiles"
   		    if (xi == 0) {
				ylab = paste("Exponential", ylab) }
			if (xi != 0) {
				ylab = paste("GPD(xi=", xi, ") ",  ylab, sep = "") }
			main = "Exploratory QQ Plot" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		plot(sort(x), y, xlab = xlab, ylab = ylab, main = main, ...)
		if (line) abline(lsfit(sort(x), y)) 
		if (labels) grid()
	}
	
	# Result:
	result = data.frame(x = sort(x), y)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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


mxfPlot = 
function (x, tail = 0.05, doplot = TRUE, labels = TRUE, ...)     
{	# A function written by D. Wuertz
	
	# Description:
	#	Creates a simple mean excess function plot.
	
	# FUNCTION:
	
	# Convert x to a vector, if the input is a data.frame.
	if(is.data.frame(x)) x = x[,1] 
	u = rev(sort(x))
	n = length(x)
	u = u[1:floor(tail*n)]
	n = length(u)
	e = (cumsum(u)-(1:n)*u)/(1:n)
	
	# Plot
	if (doplot) {
		if (labels) {
		    xlab = "Threshold: u"
			ylab = "Mean Excess: e"
			main = "Mean Excess Function" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		plot (u, e, xlab = xlab, ylab = ylab, main = main, ...) 
		if (labels) grid()
	}
	
	# Result:
	result = data.frame(threshold = u, excess = e)
	
	# Return Values:
	if (doplot) return(invisible(result)) else return(result)
}


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


mrlPlot =
function(x, conf = 0.95, umin = NA, umax=NA, nint = 100, 
doplot = TRUE, plottype = c("autoscale", ""), labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz
	
	# Description:
	#	Create a mean residual life plot with
	#	confidence intervals.
	
	# Note:
	#	"autoscale" added by DW.
	
	# References:
	#	A function originally written by S. Coles
	
	# FUNCTION:
	
	# Settings:
	plottype = plottype[1]
	if (plottype == "autoscale") {
		autoscale = TRUE }
	else {
		autoscale = FALSE }
	
	# Convert x to a vector, if the input is a data.frame.
	if (is.data.frame(x)) x = x[,1] 
	if (is.na(umin)) umin = mean(x)
	if (is.na(umax)) umax = max(x)
	sx = xu = xl = rep(NA, nint)
	u = seq(umin, umax, length = nint)
	for(i in 1:nint) {
		x = x[x >= u[i]]
		sx[i] = mean(x - u[i])
		sdev = sqrt(var(x))
		n = length(x)
		xu[i] = sx[i] + (qnorm((1 + conf)/2) * sdev)/sqrt(n)
		xl[i] = sx[i] - (qnorm((1 + conf)/2) * sdev)/sqrt(n) }
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Threshold: u"
			ylab = "Mean Excess: e"
			main = "Mean Residual Live Plot" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		if (autoscale) {
			ylim = c(min(xl[!is.na(xl)]), max(xu[!is.na(xu)]))
			plot(u, sx, type = "l", lwd = 2, xlab = xlab, 
				ylab = ylab, ylim = ylim, main = main, ...) }
		else {
			plot(u[!is.na(xl)], sx[!is.na(xl)], type = "l", 
				lwd = 2, xlab = xlab, ylab = ylab, main = main, ...) } 
		lines(u[!is.na(xl)], xl[!is.na(xl)], col = "steelblue3")
		lines(u[!is.na(xu)], xu[!is.na(xu)], col = "steelblue3")
		if (labels) grid() 
	}
	
	# Result
	result = data.frame(threshold = u, mrl = sx)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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


mePlot =
function(x, doplot = TRUE, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz
	
	# Description:
	#	Create a Mean Excess Plot
	
	# Reference:
	# 	A function imported from R-package evir
	
	# FUNCTION:
	
	# Settings:
	omit = 0
	
	# Internal Function:
	myrank = function(x, na.last = TRUE){
		ranks = sort.list(sort.list(x, na.last = na.last))
		if(is.na(na.last))
			x = x[!is.na(x)]
		for(i in unique(x[duplicated(x)])) {
			which = x == i & !is.na(x)
			ranks[which] = max(ranks[which]) }
		ranks }
	
	# Convert x to a vector, if the input is a data.frame.
	if(is.data.frame(x)) x = x[,1] 
	x = as.numeric(x)
	n = length(x)
	x = sort(x)
	n.excess = unique(floor(length(x) - myrank(x)))
	points = unique(x)
	nl = length(points)
	n.excess = n.excess[-nl]
	points = points[-nl]
	excess = cumsum(rev(x))[n.excess] - n.excess * points
	y = excess/n.excess
	xx = points[1:(nl-omit)] 
	yy = y[1:(nl-omit)]
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Threshold: u"
			ylab = "Mean Excess: e"
			main = "Mean Excess Plot" }
		else {
			xlab = ""
			ylab = ""
			main = "" }
		plot(xx, yy, xlab = xlab, ylab = ylab, main = main, ...) 
		if (labels) grid()
	}
	
	# Results:
	result = data.frame(threshold = xx, me = yy)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
	
}


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


recordsPlot = 
function(x, conf = 0.95, doplot = TRUE, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz
	
	# Description:
	#	Creates a records plot.
	
	# Note:
	#	A function imported from R-package evir,
	#	original name in EVIR: records

	# FUNCTION:
	
	# Settings:
	conf.level = conf
	
	# Convert x to a vector, if the input is a data.frame.
	if (is.data.frame(x)) x = x[,1] 
	
	# Records:
	record = cummax(x)
	expected = cumsum(1/(1:length(x)))
	se = sqrt(expected - cumsum(1/((1:length(x))^2)))
	trial = (1:length(x))[!duplicated(record)]
	record = unique(record)
	number = 1:length(record)
	expected = expected[trial]
	se = se[trial]
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Trials"
       		ylab = "Records"
       		main = "Plot of Record Development" }
		else {
			xlab = ""
			ylab = ""
			main = "" }		
        ci = qnorm(0.5 + conf.level/2)
       	upper = expected + ci * se
       	lower = expected - ci * se
       	lower[lower < 1] = 1
       	yr = range(upper, lower, number)	
       	plot(trial, number, log = "x", ylim = yr, 
				xlab = xlab, ylab = ylab, main = main, ...) 
		lines(trial, expected)
		lines(trial, upper, lty = 2)
		lines(trial, lower, lty = 2) 
		if (labels) grid()
	}
		
	# Result:
	result = data.frame(number, record, trial, expected, se)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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


ssrecordsPlot = 
function (x, subsamples = 10, doplot = TRUE, plottype = c("lin", "log"),
labels = TRUE,  ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Creates a plot of records on subsamples.
	
	# note:
	#	Changes:
	#	2003/09/06 - argument list made consistent
	
	# FUNCTION:
	
	# Convert x to a vector, if the input is a data.frame.
	if(is.data.frame(x)) x = x[,1] 
	
	# Plot type:
	plottype = plottype[1]
	
	# Records:
	save = x 
	cluster = floor(length(save)/subsamples)
	records = c()
	for (i in 1:subsamples) {
		x = save[((i-1)*cluster+1):(i*cluster)]
		y = 1:length(x)
		u = x[1]
		v = x.records = 1
		while (!is.na(v)) {
			u = x[x > u][1]
			v = y[x > u][1]
			if(!is.na(v)) x.records = c(x.records, v) }	
		if (i == 1) {
			nc = 1:length(x)
			csmean = cumsum(1/nc)
			cssd = sqrt(cumsum(1/nc-1/(nc*nc)))
			ymax = csmean[length(x)]+2*cssd[length(x)]
			# Plot:
			if (doplot) {
				if (plottype == "log") nc = log(nc)
				if (labels) {
					if (plottype == "lin") xlab = "n"
					if (plottype == "log") xlab = "log(n)"
					ylab = "N(n)" }
					main = "Subsample Records Plot"
					plot (nc, csmean+cssd, type = "l", ylim = c(0, ymax),
						xlab = xlab, ylab = ylab, main = main, ...) }
				else {
					plot (nc, csmean+cssd, type = "l", ylim = c(0, ymax),
					 	...) } 
				lines(nc, csmean)  
				lines(nc, csmean-cssd) } 
		y.records = 1:length(x.records)
		x.records = x.records[y.records < ymax]
		if (doplot) {
			if (plottype == "log") x.records = log(x.records)
			points(x.records, y.records[y.records<ymax], col=i) }
		records[i] = y.records[length(y.records)]}
	
	# Result:
	subsample = 1:subsamples
	result = data.frame(subsample, records)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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


msratioPlot = 
function (x, p = 1:4, doplot = TRUE, plottype = c("autoscale", ""),
labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz
	
	# Description:
	#	Creates a Plot of maximum and sum ratio.
	
	# FUNCTION:
	
	# Settings:
	plottype = plottype[1]
	if (plottype == "autoscale") {
		autoscale = TRUE }
	else {
		autoscale = FALSE }
	if (autoscale) ylim = c(0,1)
	
	# Convert x to a vector, if the input is a data.frame.
	if(is.data.frame(x)) x = x[,1] 
	
	# Plot:
	if (doplot) {
		if (labels) {
		    xlab = "Trials"
       		ylab = "Records"
       		main = "Plot of Record Development" }
		else {
			xlab = ""
			ylab = ""
			main = "" }				
		if (autoscale) {
			plot(c(0, length(x)), y = ylim, xlab = xlab, 
					ylab = ylab, main = main, type = "n", ...) }
		else {
			plot(c(0, length(x)), xlab = xlab, 
					ylab = ylab, main = main, type = "n", ...) }
		if (labels) grid()
  	}
  	
	# Color numbering:
	i = 1
	
	# Suppress warnings for points outside the frame:
	ratios = matrix(rep(0, times=length(x)*length(p)), byrow=TRUE, 
		ncol=length(p))
	if (doplot) par(err=-1)
	
	# Loop over all exponents p:
  	for (q in p) {
    	rnp = cummax(abs(x)^q) / cumsum(abs(x)^q)
    	i = i + 1
		ratios[,q] = rnp
    	if (doplot) lines (rnp, col=i) }

	# Result:
	result = data.frame(ratios)
	
	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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


xacfPlot =
function(x, threshold = 0.95, lag.max = 15, doplot = TRUE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Creates plots of exceedences, one for the
	#	heights and one for the distances.
  	
	# FUNCTION:
	
	# Settings:
	# Sorry, user specified labels not yet implemented.
	labels = TRUE
	if (labels) {
		xlab = c("Index", "Lag")
		ylab = c("Heights", "Distances", "ACF")
		main = c("Heights over Threshold", "Distances between Heights", 
			"Series Heights", "Series Distances") }
			
	# Convert x to a vector, if the input is a data.frame.
	if (is.data.frame(x)) x = x[,1] 
	# Heights/Distances
    threshold = sort(x)[round(threshold*length(x))]
   	Heights = (x-threshold)[(x-threshold)>0]
   	Distances = diff((1:length(x))[(x-threshold)>0])
 	
	# Plot:
   	if (doplot) {
		plot (Heights, type="h", xlab = xlab[1], ylab = ylab[1], 
			main = main[1], ...)
   		plot (Distances,type="h", xlab = xlab[1], ylab = ylab[2], 
			main = main[2], ...) }
  	
	# Correlations:
	Heights = as.vector(acf(Heights, lag.max=lag.max, plot = doplot, 
		xlab = xlab[2], ylab = ylab[3], main = main[3], ...)$acf)
   	Distances = as.vector(acf(Distances, lag.max=lag.max, plot = doplot, 
   		xlab = xlab[2], ylab = ylab[3], main = main[4], ...)$acf)

	# Result:
	lag = as.vector(0:(lag.max))
	result = data.frame(lag, Heights, Distances)

	# Return Value:
	if (doplot) return(invisible(result)) else return(result)
}


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

	    
interactivePlot =
function(x, choices = paste("Plot", 1:9), 
plotFUN = paste("plot.", 1:9, sep = ""), which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for an object of class "template".
	
	# Arguments:
	#	x - an object to be plotted
	#	choices - the character string for the choice menu
	# 	plotFUN - the names of the plot functions
	#   which - plot selection, which graph should be 
	#	  displayed. If a character string named "ask" the 
	#     user is interactively asked which to plot, if
	#	  a logical vector of length N, those plots which
	#	  are set "TRUE" are displayed, if a character string
	#	  named "all" all plots are displayed.
	
	# Note:
	#	At maximum 9 plots are supported.

	# FUNCTION:
	
	# Some cecks:
	if (length(choices) != length(plotFUN)) 
		stop("Arguments choices and plotFUN must be of same length.")
	if (length(which) > length(choices)) 
		stop("Arguments which has incoorect length.")
	if (length(which) > length(plotFUN)) 
		stop("Arguments which has incorrect length.")
	if (length(choices) > 9)
		stop("Sorry, only 9 plots at max are supported.")
	
	# Internal "askPlot" Function:          	  
    multPlot = function (x, choices, ...) {
	    # Selective Plot:
		selectivePlot = function (x, choices, FUN, which){
			# Internal Function:
			askPlot = function (x, choices, FUN) {
				# Pick and Plot:
				pick = 1; n.plots = length(choices)
				while (pick > 0) { pick = menu (
					choices = paste("plot:", choices), 
					title = "\nMake a plot selection (or 0 to exit):")
				    if (pick > 0) match.fun(FUN[pick])(x) } }			        
			if (as.character(which[1]) == "ask") {
				askPlot(x, choices = choices, FUN = FUN, ...) }
			else { 
				for (i in 1:n.plots) if (which[i]) match.fun(FUN[i])(x) }
			invisible() }  
		# match Functions, up to nine ...
		if (length(plotFUN) < 9) plotFUN = 
			c(plotFUN, rep(plotFUN[1], times = 9 - length(plotFUN)))
		plot.1 = match.fun(plotFUN[1]); plot.2 = match.fun(plotFUN[2]) 
		plot.3 = match.fun(plotFUN[3]); plot.4 = match.fun(plotFUN[4]) 
		plot.5 = match.fun(plotFUN[5]); plot.6 = match.fun(plotFUN[6]) 
		plot.7 = match.fun(plotFUN[7]); plot.8 = match.fun(plotFUN[8]) 
		plot.9 = match.fun(plotFUN[9])   	
		pick = 1
	    while (pick > 0) { pick = menu (
	        ### choices = paste("plot:", choices),
	        choices = paste(" ", choices), 
	        title = "\nMake a plot selection (or 0 to exit):")
	        # up to 9 plot functions ...
	        switch (pick, plot.1(x), plot.2(x), plot.3(x), plot.4(x), 
	        	plot.5(x), plot.6(x), plot.7(x), plot.8(x), plot.9(x) ) } }
	        		          
	# Plot:
	if (which[1] == "all") which = rep(TRUE, times = length(choices))
	if (which[1] == "ask") {
		multPlot(x, choices, ...) }
	else {
		for ( i in 1:length(which) ) {
			FUN = match.fun(plotFUN[i])
			if (which[i]) FUN(x) } }
			
	# Return Value:
	invisible(x)
}


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


gridVector =
function(x, y)
{   # A function implemented by Diethelm Wuertz, GPL

    # Description:
    #   Creates from two vectors x and y all grid points
    
    # Details: 
    #   The two vectors x and y span a rectangular grid with nx=length(x) 
    #   times ny=length(y) points which are returned as a matrix of size
    #   (nx*ny) times 2.
    
    # Arguments:
    #   x, y - two numeric vectors of length m and n which span the 
    #   rectangular grid of size m times n.
    
    # Value:
    #   returns a list with two elements X and Y each of length m 
    #   times n
    
    # Example:
    #   > gridVector(1:3, 1:2)
    #             [,1] [,2]
    #       [1,]    1    1
    #       [2,]    2    1
    #       [3,]    3    1
    #       [4,]    1    2
    #       [5,]    2    2
    #       [6,]    3    2

    # FUNCTION: 
    
    # Prepare for Input:
    nx = length(x)
    ny = length(y)
    xoy = cbind(rep(x, ny), as.vector(matrix(y, nx, ny, byrow = TRUE)))
    X = matrix(xoy, nx * ny, 2, byrow = FALSE)
    
    # Return Value:
    list(X = X[,1], Y = X[,2])
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION			DESCRIPTION:
#  findThreshold     Finds threshold values
#  blocks            Creates data blocks on vectors and time series
#  blockMaxima		 Calculates block maxima on vectors and time series
#  deCluster 		 Declusters a point process
################################################################################


findThreshold =
function(x, n = NA)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Finds upper thresold for a given number of Extremes.
	
	# Arguments:
	# 	n 	- a numeric value or vector giving number of extremes 
	#         above the threshold. If "n" is not specified, "n"
	#		  is set to an integer representing 5% of the data 
	#         from the whole data set "x".
		
	# Note:
	#	Imported from R-package evir/EVIS.

	# FUNCTION:
	
	# Settings:
	if(is.na(n[1])) n = floor(0.05*length(x))
	
	# Continue:
	x = rev(sort(as.numeric(x)))
	thresholds = unique(x)
	indices = match(x[n], thresholds)
	indices = pmin(indices + 1, length(thresholds))	
	
	# Return Value:
	thresholds[indices]
}


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


blocks =
function(x, block = "month", FUN = max)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Creates data blocks on vectors and time series.

	# Note:
	#	Imported from R-package evir/EVIS.
	
	# FUNCTION:
	
    # Settings:
    data = x
    
    # Compute:
    n.all = length(data)
    if(is.character(block)) {
     	times = as.POSIXlt(attributes(data)$times) 
      	if(block %in% c("semester", "quarter")) {
      		sem = quart = times$mon
           	sem[sem %in% 0:5] = quart[quart %in% 0:2] = 0
           	sem[sem %in% 6:11] = quart[quart %in% 3:5] = 1
           	quart[quart %in% 6:8] = 2
          	quart[quart %in% 9:11] = 3 }
        grouping = switch(block,
          	semester = paste(times$year, sem),
           	quarter = paste(times$year, quart),
           	quarters = paste(times$year, quart),
          	month = paste(times$year, times$mon),
          	months = paste(times$year, times$mon),
          	year = times$year,
          	years = times$year,
          	stop("unknown time period"))
      	newdata = tapply(data, grouping, FUN=FUN) }
   	else {
       	data = as.numeric(data)
       	nblocks = (length(data) %/% block) + 1
       	grouping = rep(1:nblocks, rep(block, nblocks))[1:length(data)]
       	newdata = tapply(data, grouping, FUN=FUN)}
       
  	# Return Value: 
    result = newdata 
    result
}


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


blockMaxima = 
function(x, block = "month", details = FALSE, doplot = TRUE, ...) 
{  	# A function implemented by Diethelm Wuertz
	
	# Description:
	#	Calculates block maxima on vectors and time series.
	
	# Arguments:
	#	x 		- may be alternatively as.vector or as.ts
	#   block 	- as.numeric:   length of a block
	#			  as.character: year | semester | quarter | month
	
	# Note:
	#   Calls McNeils Splus function blocks()
	#	Output data as vector of transposed 
	#	result to get proper order of data!

	# FUNCTION:
	
	# Settings
	x = blocks(x, block)
  	
	# Plot:
	if (doplot) {
		plot(as.vector(x), type="h", ylab = "Block Maxima",  ...)
		title(main = paste(block, "- Block Maxima"))
		grid() }
  	
	# Details:
	# if details == FALSE a vector is returned, i.e details are removed!
	if (!details) x = as.vector(x[is.na(x) == FALSE])
	
	# Return Value:
	x
}


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


deCluster = 
function(x, run = NA, doplot = TRUE)
{	# A function implemented by Diethelm Wuertz 

	# Description:
	#	Declusters a point process
	
	# Note:
	#	Imported from R-package evir/EVIS.
	
	# FUNCTION:
	
	# Settings:
	labels = TRUE
	
	# Imported Function:
	series = x
	picture = doplot
	n = length(as.numeric(series))
    times = attributes(series)$times
    if (is.null(times)) 
        stop("`series' must have a `times' attribute")
    as.posix = is.character(times) || inherits(times, "POSIXt") || 
        inherits(times, "date") || inherits(times, "dates")
    if (as.posix) 
        gaps = as.numeric(difftime(as.POSIXlt(times)[2:n], 
        as.POSIXlt(times)[1:(n - 1)], units = "days"))
    else gaps = as.numeric(diff(times))
    longgaps = gaps > run
    if (sum(longgaps) <= 1) 
        stop("Decluster parameter too large")
    cluster = c(0, cumsum(longgaps))
    cmax = tapply(as.numeric(series), cluster, max)
    newtimes = times[match(cmax, series)]
    newseries = structure(series[match(cmax, series)], times = newtimes)
    n = length(as.numeric(newseries))
    if (as.posix) {
        newgaps = as.numeric(difftime(as.POSIXlt(newtimes)[2:n], 
            as.POSIXlt(newtimes)[1:(n - 1)], units = "days"))
        times = as.POSIXlt(times)
        newtimes = as.POSIXlt(newtimes) }
    else {
    	newgaps = as.numeric(diff(newtimes)) }
    
    # Plot:
    if (doplot) {
        # cat("Declustering picture...\n")
        # cat(paste("Data reduced from", length(as.numeric(series)), 
        #     "to", length(as.numeric(newseries)), "\n"))
        # par(mfrow = c(2, 2))
        if (labels) {
            main = "de-Clustering"
        	plot(times, series, type = "h", main = main)
        	qPlot(gaps)
        	plot(newtimes, newseries, type = "h", main = main)
        	qPlot(newgaps) }
    }
    
    # Result:
    ans = newseries

	# Return Value:
	ans
}   
  
 
# ******************************************************************************


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:				GEV DISTRIBUTION FAMILY: [USE FROM EVD]
#  devd					 Density for the GEV Distribution 
#   pevd				  Probability for the GEV Distribution
#   qevd				  Quantiles for the GEV Distribution
#   revd				  Random variates for the GEV Distribution
# FUNCTION:				GEV DISTRIBUTION FAMILY: [USE FROM EVIS]
#  dgev					 Density for the GEV Distribution 
#   pgev				  Probability for the GEV Distribution
#   qgev				  Quantiles for the GEV Distribution
#   rgev				  Random variates for the GEV Distribution
################################################################################


devd = 
function (x, loc = 0, scale = 1, shape = 0, log = FALSE) 
{
    # FUNCTION:
    
    if (min(scale) <= 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    x = (x - loc)/scale
    if (shape == 0) 
        d = log(1/scale) - x - exp(-x)
    else {
        nn = length(x)
        xx = 1 + shape * x
        xxpos = xx[xx > 0 | is.na(xx)]
        scale = rep(scale, length.out = nn)[xx > 0 | is.na(xx)]
        d = numeric(nn)
        d[xx > 0 | is.na(xx)] = log(1/scale) - xxpos^(-1/shape) - 
            (1/shape + 1) * log(xxpos)
        d[xx <= 0 & !is.na(xx)] = -Inf
    }
    if (!log) 
        d = exp(d)
    d
}


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


pevd =
function (q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE) 
{
    # FUNCTION:
    
    if (min(scale) <= 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    q = (q - loc)/scale
    if (shape == 0) 
        p = exp(-exp(-q))
    else p = exp(-pmax(1 + shape * q, 0)^(-1/shape))
    if (!lower.tail) 
        p = 1 - p
    p
}


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


qevd = 
function (p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE) 
{
    # FUNCTION:
    
    if (min(p, na.rm = TRUE) <= 0 || max(p, na.rm = TRUE) >= 
        1) 
        stop("`p' must contain probabilities in (0,1)")
    if (min(scale) < 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    if (!lower.tail) 
        p = 1 - p
    if (shape == 0) 
        return(loc - scale * log(-log(p)))
    else return(loc + scale * ((-log(p))^(-shape) - 1)/shape)
}


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


revd =
function (n, loc = 0, scale = 1, shape = 0) 
{
    # FUNCTION:
    
    if (min(scale) < 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    if (shape == 0) 
        return(loc - scale * log(rexp(n)))
    else return(loc + scale * (rexp(n)^(-shape) - 1)/shape)
}


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


dgev =
function(x, xi = 1, mu = 0, sigma = 1, log = FALSE)
{	# A function implemented from evd

  	# Description:
	#   GEV Density Function
	#   Note: 1 + xi*(x-mu)/sigma > 0
	#   xi > 0 Frechet
	#   xi = 0 Gumbel
	#   xi < 0 weibl

	# FUNCTION:
	
	# Settings:
	loc = mu
	scale = sigma
	shape = xi
	
	# Density function:
	if (min(scale) <= 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    x = (x - loc)/scale
    if (shape == 0) 
        d = log(1/scale) - x - exp(-x)
    else {
        nn = length(x)
        xx = 1 + shape * x
        xxpos = xx[xx > 0 | is.na(xx)]
        scale = rep(scale, length.out = nn)[xx > 0 | is.na(xx)]
        d = numeric(nn)
        d[xx > 0 | is.na(xx)] = log(1/scale) - xxpos^(-1/shape) - 
            (1/shape + 1) * log(xxpos)
        d[xx <= 0 & !is.na(xx)] = -Inf
    }
    if (!log) 
        d = exp(d)
	
	# Return Value:
	d
}


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


pgev =
function(q, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
{	# A function implemented from evd
 
  	# Description:
	#   GEV Probability Function
	#   Note: 1 + xi*(x-mu)/sigma > 0
	#   xi > 0 Frechet
	#   xi = 0 Gumbel
	#   xi < 0 Weibull

	# FUNCTION:
	
	# Settings:
	loc = mu
	scale = sigma
	shape = xi
	
	# Probability function:
	if (min(scale) <= 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    q = (q - loc)/scale
    if (shape == 0) 
        p = exp(-exp(-q))
    else p = exp(-pmax(1 + shape * q, 0)^(-1/shape))
    if (!lower.tail) 
        p = 1 - p
	
	# Return Value:
	p
}


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


qgev =
function (p, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
{	# A function implemented from evd

	# Description:
	#   GEV Quantile Function
	#   Note: 1 + xi*(x-mu)/sigma > 0
	#   xi > 0 Frechet
	#   xi = 0 Gumbel
	#   xi < 0 Weibull

	# FUNCTION:
	
	# Settings:
	loc = mu
	scale = sigma
	shape = xi
	
	# Return Value:
	if (min(p, na.rm = TRUE) < 0 || max(p, na.rm = TRUE) > 1) 
        stop("`p' must contain probabilities in (0,1)")
    if (min(scale) < 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    if (!lower.tail) 
        p = 1 - p
    if (shape == 0) 
        q = loc - scale * log(-log(p))
    else 
    	q = loc + scale * ((-log(p))^(-shape) - 1)/shape
    	
    # Return Value:
    q
}


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


rgev =
function (n, xi = 1, mu = 0, sigma = 1)
{	# A function implemented from evd

	# Description:
	#   GEV Random Variables
	#   Note: 1 + xi*(x-mu)/sigma > 0
	#   xi > 0 Frechet
	#   xi = 0 Gumbel
	#   xi < 0 Weibull

	# FUNCTION:
	
	# Settings:
	loc = mu
	scale = sigma
	shape = xi
	
	# Return Value:
	if (min(scale) < 0) 
        stop("invalid scale")
    if (length(shape) != 1) 
        stop("invalid shape")
    if (shape == 0) 
        r = loc - scale * log(rexp(n))
    else 
    	r = loc + scale * (rexp(n)^(-shape) - 1)/shape

	# Return Value:
	r
}
	
	
# ******************************************************************************



# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


# ##############################################################################
# FUNCTION:				   GEV MODELLING FROM EVIS:
#  gevSim					Simulates GEV including Gumbel rvs [EVIS/EVIR]
#  gevFit					Fits GEV Distribution
#   print.gevFit 		  	 Print Method for object of class "gevFit"
#   plot.gevFit 		  	 Plot Method for object of class "gevFit"
#   summary.gevFit       	 Summary Method for object of class "gevFit"
# FUNCTION:                ADDITIONAL PLOT:
#  gevrlevelPlot     		Calculates Return Levels Based on GEV Fit
################################################################################


gevSim = 
function(model = list(shape = 0.25, location = 0, scale = 1), n = 1000)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Generates random variates from a GEV distribution
	
	# FUNCTION:
	
	# Simulate:
	ans = rgev(n = n, xi = model$shape, mu = model$location, 
		sigma = model$scale)
		
	# Return Value:
	ans	
}


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


gevFit =
function(x, block = NA, type = c("mle", "pwm"), gumbel = FALSE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Fits parameters to a GEV distribution
	
	# Note:
	#	Argument named "method is already used for the selection
	#	of the MLE optimization algorithm, therfore we use here
	#	"type".
	
	# FUNCTION:
	
	# Settings:
	call = match.call()
	type = type[1]
	
	# Internal Function:
	gev.pwm = function(data, block = NA, ...) {
	# Probability Weighted Moment method.
    # Blocks and data:
    n.all = NA
    if (!is.na(block)) {
        n.all = length(data)
        if (is.character(block)) {
            times = as.POSIXlt(attributes(data)$times)
            if (block %in% c("semester", "quarter")) {
                sem = quart = times$mon
                sem[sem %in% 0:5] = quart[quart %in% 0:2] = 0
                sem[sem %in% 6:11] = quart[quart %in% 3:5] = 1
                quart[quart %in% 6:8] = 2
                quart[quart %in% 9:11] = 3 }
            grouping = switch(block, 
            	semester = paste(times$year, sem), 
            	quarter = paste(times$year, quart), 
            	month = paste(times$year, times$mon), 
            	year = times$year, 
            	stop("unknown time period"))
            data = tapply(data, grouping, max) }
        else {
            data = as.numeric(data)
            nblocks = (length(data)%/%block) + 1
            grouping = rep(1:nblocks, rep(block, nblocks))[1:length(data)]
            data = tapply(data, grouping, max) } }
    data = as.numeric(data)
    n = length(data)	
	
    # Internal Function - Sample Moments:
	sampwm = function (x, nmom) {
		# a = 0, b = 0, kind = 1
		x = rev(sort(x))
		moments = rep(0, nmom)
		moments[1] = mean(x)
		n = length(x)
		for (i in 1:n) {
			weight = 1/n
			for (j in 2:nmom) {
				weight = weight*(n-i-j+2)/(n-j+1)
				moments[j] = moments[j] + weight*x[i] } }
		return(moments) }
	
	# Internal Function:
	y = function(x, w0, w1, w2) { (3^x-1)/(2^x-1) - (3*w2 - w0)/(2*w1 - w0) }		
	# Calculate:
	w = sampwm(data, nmom = 3)
	w0 = w[1]
	w1 = w[2]
	w2 = w[3]	   
	xi = uniroot(f = y, interval = c(-5,+5), 
		w0 = w[1], w1 = w[2], w2 = w[3])$root
	sigma = beta = (2*w1-w0)*xi / gamma(1-xi) / (2^xi-1)
	mu = w0 + beta*(1-gamma(1-xi))/xi
	# Output:
	fit = list(n.all = n.all, n = n, data = data, bock = block, 
		par.ests = c(xi, sigma, mu), par.ses = rep(NA, 3),
		varcov = matrix(rep(NA, 9), 3, 3), converged = NA, 
		nllh.final = NA, call=match.call(), selected = "pwm")
	names(fit$par.ests) = c("xi", "sigma", "mu")
    names(fit$par.ses) = c("xi", "sigma", "mu")	
    # Return Value:
    class(fit) = "gev" 
	fit }
	
	# Internal Function:
	gumbel.pwm = function(data, block = NA, ...) {
	# "Probability Weighted Moment" method.
	# Blocks and data:
    n.all = NA
    if (!is.na(block)) {
        n.all = length(data)
        if (is.character(block)) {
            times = as.POSIXlt(attributes(data)$times)
            if (block %in% c("semester", "quarter")) {
                sem = quart = times$mon
                sem[sem %in% 0:5] = quart[quart %in% 0:2] = 0
                sem[sem %in% 6:11] = quart[quart %in% 3:5] = 1
                quart[quart %in% 6:8] = 2
                quart[quart %in% 9:11] = 3 }
            grouping = switch(block, 
            	semester = paste(times$year, sem), 
            	quarter = paste(times$year, quart), 
            	month = paste(times$year, times$mon), 
            	year = times$year, 
            	stop("unknown time period"))
            data = tapply(data, grouping, max) }
        else {
            data = as.numeric(data)
            nblocks = (length(data)%/%block) + 1
            grouping = rep(1:nblocks, rep(block, nblocks))[1:length(data)]
            data = tapply(data, grouping, max) } }
    data = as.numeric(data)
    n = length(data)
	# Sample Moments:
    x = rev(sort(data))
	lambda = c(mean(x), 0)
	for (i in 1:n) {
		weight = (n-i)/(n-1)/n
		lambda[2] = lambda[2] + weight*x[i] } 
	# Calculate Parameters:
    xi = 0
	sigma = beta = lambda[2]/log(2)
	mu = lambda[1] - 0.5772*beta
	# Output:
	fit = list(n.all = n.all, n = n, data = data, block = block, 
		par.ests = c(sigma, mu), par.ses = rep(NA, 2),
		varcov = matrix(rep(NA, 4), 2, 2), converged = NA, 
		nllh.final = NA, call = match.call(), selected = "pwm")
	names(fit$par.ests) = c("sigma", "mu")
    names(fit$par.ses) = c("sigma", "mu")
    # Return Value:
    class(fit) = "gev" # not gumbel!
	fit }
	
	# Estimate Parameters:
	if (gumbel) {	
		# Add Call and Type
		if (length(type) > 1) type = type[1]
		# Probability Weighted Moment Estimation
		if (type == "pwm") {
			fitted = gumbel.pwm(data = x, block = block, ...) }
		# Maximum Log Likelihood Estimation
		# Use Alexander McNeils EVIS:
		if (type == "mle") { 
			fitted = gumbel(data = x, block = block, ...) } }	
	else {
		# Add Call and Type
		if (length(type) > 1) type = type[1]
		# Probability Weighted Moment Estimation:
		if (type == "pwm") { 
			fitted = gev.pwm(data = x, block = block, ...) }
		# Maximum Log Likelihood Estimation
		# Use Alexander McNeils EVIS (renames as gev.mle)
		if (type == "mle") { 
			fitted = gev(data = x, block = block, ...) }	}
			
	# Compute Residuals:
	if (gumbel) {
		# GUMBEL:
		xi = 0
		sigma = fitted$par.ests[1]
		mu = fitted$par.ests[2] 
		fitted$residuals = exp( - exp( - (fitted$data - mu)/sigma)) }
	else {
		# GEV:
		xi = fitted$par.ests[1]
		sigma = fitted$par.ests[2]
		mu = fitted$par.ests[3]
		fitted$residuals = (1 + (xi * (fitted$data - mu))/sigma)^(-1/xi) }	
		
	# Make Unique:
	fit = list()
	fit$fit = fitted
	fit$call = call
	fit$type = c(if(gumbel) "gum" else "gev", type[1])
	fit$par.ests = fitted$par.ests
	fit$par.ses = fitted$par.ses
	fit$residuals = fitted$residuals
	fit$fitted.values = fitted$data - fitted$residuals
	fit$llh = fitted$nllh.final
	fit$converged = fitted$converged
		
	# Return Value:
	class(fit) = "gevFit"
	fit
}


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


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

	# Description:
	#	Print Method for an object of class "gevFit".
	
	# FUNCTION:
	
	# Function Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "")
	 		
	# Estimation Type:
	cat("\nEstimation Type:", x$type, "\n")	
	
	# Estimated Parameters:
	cat("\nEstimated Parameters:\n")
	print(x$par.ests)
	cat("\n")
	
	# Return Value:
	invisible(x)
}


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


plot.gevFit =
function(x, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for an object of class "gevFit".
	
	# Details:
	# 	plot.gev:
	#	Data are converted to unit exponentially distributed residuals
    #	under null hypothesis that GEV fits. Two diagnostics for iid
    #	exponential data are offered:
    #	"Scatterplot of Residuals" and "QQplot of Residuals"

	# FUNCTION:	
	
	# Internal Plot Functions:
	plot.1 <<- function(x) {
		plot(x$fit$data, type = "h", 
			xlab = "Index",
			ylab = "Data", 
			main = "Block Maxima") }
	plot.2 <<- function(x) {
	    plot(x$residuals, 
	    	xlab = "Ordering",
	    	ylab = "Residuals",  
	    	main = "Scatterplot of Residuals")
	   	lines(lowess(1:length(x$residuals), x$residuals)) }
	plot.3 <<- function(x) {            
		# evir::qplot
		qplot(x$residuals, 
			# xlab = "Ordered Data",
			# ylab = "Exponential Quantiles",
	 		main = "Quantile-Quantile Plot") }
	 		
	# Plot:
	interactivePlot(
		x = x,
		choices = c(
			"Block Maxima Plot", 
			"Scatterplot of Residuals", 
			"Quantile Quantile Plot"),
		plotFUN = c(
			"plot.1", 
			"plot.2", 
			"plot.3"),
		which = which)
	            	
	# Return Value:
	invisible(x)
}


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

summary.gevFit =
function(object, doplot = TRUE, which = "all", ...) 
{
	# A function implemented by Diethelm Wuertz

	# Description:
	#	Summary method for an object of class "gevFit".

    # FUNCTION:
    
    # Print:
    print(object, ...)
	
	# Summary:
	if (object$type[2] == "mle") {
		cat("\nStandard Deviations:\n"); print(object$par.ses)
		cat("\nLog-Likelihood Value: ", object$llh)
		cat("\nType of Convergence:  ", object$converged, "\n") } 
	cat("\n")
	
	# Plot:
	if (doplot) plot(object, which = which, ...)
	cat("\n")
	
	# Return Value:
	invisible(object)
}


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


gevrlevelPlot =
function(object, k.blocks = 20, add = FALSE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Calculates Return Levels Based on GEV Fit
	
	# FUNCTION:
	
	# Settings
	fit = object
	
	# Use "rlevel.gev":
	ans = rlevel.gev(out = fit$fit, k.blocks = k.blocks, add = add, ...)
	ans = c(min = ans[1], v = ans[2], max = ans[3])
	
	# Return Value:
	ans	
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:				   GEV MODELLING FROM ISMEV:
#  gevglmFit				Fits GEV Distribution
#   print.gevglmFit		     Print Method for object of class "gevglm"
#   plot.gevglmFit		     Plot Method for object of class "gevglm"
#   summary.gevglmFit		 Summary Method for object of class "gevglm"
# FUNCTION:                ADDITIONAL PLOTS:
#  gevglmprofPlot			Profile Log-likelihoods for Stationary GEV Models
#  gevglmprofxiPlot	        Profile Log-likelihoods for Stationary GEV Models
################################################################################


gevglmFit =
function(x, y = NULL, gumbel = FALSE, mul = NULL, sigl = NULL, shl = NULL, 
mulink = identity, siglink = identity, shlink = identity, show = FALSE, 
method = "Nelder-Mead", maxit = 10000, ...)
{	# A function written by Diethelm Wuertz

	# Description:
	#	Fits GEV Distribution
	
	# Note:
	#	This is a function wrapper to the functions 'gev.fit' and
	#	'gum.fit' which are part of the R package 'ismev'.
		
	# FUNCTION:
	
	# Fit - Use gev.fit() and gum.fit() from R's ismev Package:
	call = match.call()
	if (gumbel) {
		fitted = gum.fit(xdat = x, ydat = y, mul = mul, sigl = sigl,  
			mulink = mulink, siglink = siglink, show = show, 
			method = method, maxit = maxit, ...) }
	else {
		fitted = gev.fit(xdat = x, ydat = y, mul = mul, sigl = sigl, shl = shl, 
			mulink = mulink, siglink = siglink, shlink = shlink, show = show, 
			method = method, maxit = maxit, ...) }
	fitted$gumbel = gumbel
	
	# Standard Errors and Covariance Matrix:
	if (gumbel) {	
		# Parameters - We take the same order as in gevFit:
		mle = rev(fitted$mle)
		names(mle) = c("sigma", "mu")
		se = rev(fitted$se)
		names(se) = c("sigma", "mu")
		covar = fitted$cov
		covar[1,1] = fitted$cov[2,2]
		covar[2,2] = fitted$cov[1,1] }
	else {
		# Parameters - We take the same order as in gevFit:
		mle = rev(fitted$mle)
		names(mle) = c("xi", "sigma", "mu")
		se = rev(fitted$se)
		names(se) = c("xi", "sigma", "mu")
		covar = fitted$cov
		covar[1,1] = fitted$cov[3,3]
		covar[3,3] = fitted$cov[1,1]
		covar[1,2] = covar[2,1] = fitted$cov[2,3]
		covar[2,3] = covar[3,2] = fitted$cov[1,2]  }
	fitted$covar = covar
				
	# Calculate Residuals:
	if (gumbel) {
		# GUMBEL:
		xi = 0
		sigma = mle[1]
		mu = mle[2] 
		residuals = exp( - exp( - (fitted$data - mu)/sigma)) }
	else {
		# GEV:
		xi = fitted$mle[1]
		sigma = fitted$mle[2]
		mu = fitted$mle[3]
		residuals = (1 + (xi * (fitted$data - mu))/sigma)^(-1/xi) }		
	
	# Add:
	fit = list()
	fit$fit = fitted
	fit$call = match.call()
	fit$type = c(if(gumbel) "gumglm" else "gevglm", "mle")
	fit$par.ests = mle
	fit$par.ses = se
	fit$residuals = residuals
	fit$fitted.values = x - residuals
	fit$llh = fitted$nllh
	fit$converged = fitted$conv
	
	# Return Value:
	class(fit) = "gevglmFit"
	fit	
}


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


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

	# Description:
	#	Print method for objects of class 'gevglmFit'
	
	# FUNCTION:

	# Function Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "") 		
	
	# Estimation Type:
	cat("\nEstimation Type:", x$type, "\n")	
	
	# Fitted Parameters:
	cat("\nEstimated Parameters:\n")
	print(x$par.ests)
	cat("\n")
	
	# Return Value:
	invisible(x)
}


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


plot.gevglmFit = 
function(x, which = "ask", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for objects of class 'gevglmFit'
	
	# FUNCTION:
	
	# Settings:
	fit = x
	
	# Internal "plot.n" Function:
	plot.1 <<- function(x, ...) {
		if (x$gumbel) x$mle = c(x$mle, 0)
		gev.pp(x$mle, x$data) }
	plot.2 <<- function(x, ...) {
		if (x$gumbel) x$mle = c(x$mle, 0)
		gev.qq(x$mle, x$data) }
	plot.3 <<- function(x, ...) {
		if (x$gumbel) { 
			fit$mle = c(x$mle, 0)
			gum.rl(x$mle, x$cov, x$data) }
		else { 
			gev.rl(x$mle, x$cov, x$data) } }
	plot.4 <<- function(x, ...) {
		if (x$gumbel) x$mle = c(x$mle, 0)
		gev.his(x$mle, x$data) }
	plot.5 <<- function(x, ...) {
		n = length(x$data)
		z = (1:n)/(n + 1)
		plot(z, exp( - exp( - sort(x$data))), 
	    	xlab = "empirical", ylab = "model")
		abline(0, 1, col = 4)
		title("Residual Probability Plot") }
	plot.6 <<- function(x, ...) {
		n = length(x$data)
		z = (1:n)/(n + 1)
		plot( - log( - log(z)), sort(x$data), 
	   		xlab = "empirical", ylab = "model")
		abline(0, 1, col = 4)
		title("Residual Quantile Plot (Gumbel Scale)") }
				
	# Plot:
	if (fit$fit$trans) { 
		# Non-Stationary Plots: plot 11-12
		interactivePlot(
			x = x$fit,
			choices = c(
				"Residual Probability Plot",  
				"Residual Quantile Plot"),
			plotFUN = c(
				"plot.5", 
				"plot.6"),
			which = which) }	 	
	 else { 
		# Stationary Plots: plot 01-04
		interactivePlot(
			x = x$fit,
			choices = c(
				"Residual Probability Plot", 
		    	"Residual Quantile Plot",
				"Return Level Plot", 
				"Density Plot"),
			plotFUN = c(
				"plot.1", 
				"plot.2",
				"plot.3", 
				"plot.4"),
			which = which) }
	
	# Return Value:
	invisible(x)
}


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


summary.gevglmFit =
function(object, doplot = TRUE, which = "all", ...)
{   # A function implemented by Diethelm Wuertz

	# Description:
	#	Summary method for objects of class 'gevglmFit'
	
	# FUNCTION:
		    
	# Print:
	print(object, ...)
	
	# Summary:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$converged, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot(object, which = which, ...)
	cat("\n")
	
	# Return Result
	invisible(object)
}


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


gevglmprofPlot =
function(object, m, xlow, xup, conf = 0.95, nint = 100)
{	# A function implemented by Diethelm Wuertz

  	# Description:
  	#	Profile Log-likelihoods for Stationary GEV Models.
  	
  	# FUNCTION:
  	
  	# Compute:
  	if (object$fit$gumbel) {
		stop("Not for Gumbel type distributions") }
	else {
		gev.prof(z = object$fit, m = m, xlow = xlow, xup = xup , conf = conf, 
			nint = nint) }
}


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


gevglmprofxiPlot =
function(object, xlow, xup, conf = 0.95, nint = 100)
{	# A function implemented by Diethelm Wuertz

	# Description:
  	#	Profile Log-likelihoods for Stationary GEV Models.
  	
  	# FUNCTION:
	
	# Compute:
	if (object$fit$gumbel) {
		stop("Not for Gumbel type distributions") }
	else {
		gev.profxi(z = object$fit, xlow = xlow, xup = xup, conf = conf, 
			nint = nint) }
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


# ##############################################################################
# FUNCTION:				   MDA ESTIMATORS:
#  hillPlot					Plot Hill's estimator
#  shaparmPlot				Pickands, Hill & Decker-Einmahl-deHaan Estimator
#   shaparmPickands			 Auxiliary function called by shaparmPlot
#   shaparmHill				 ... called by shaparmPlot
#   shaparmDehaan			 ... called by shaparmPlot
################################################################################


hillPlot = 
function(x, option = c("alpha", "xi", "quantile"), start = 15, end = NA,
reverse = FALSE, p = NA, ci = 0.95, autoscale = TRUE, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz
 
	# Description:
	#	Plots the results from the Hill Estimator.
	
	# Note:
	# 	Imported from R-package evir

	# Settings:
	data = as.numeric(x)
	ordered = rev(sort(data))
	ordered = ordered[ordered > 0]
	n = length(ordered)
	option = match.arg(option)
	if((option == "quantile") && (is.na(p)))
		stop("Input a value for the probability p")
	if((option == "quantile") && (p < 1 - start/n)) {
		cat("Graph may look strange !! \n\n")
		cat(paste("Suggestion 1: Increase `p' above",
			format(signif(1 - start/n, 5)), "\n"))
		cat(paste("Suggestion 2: Increase `start' above ",
			ceiling(length(data) * (1 - p)), "\n"))
	}
	k = 1:n
	loggs = logb(ordered)
	avesumlog = cumsum(loggs)/(1:n)
	xihat = c(NA, (avesumlog - loggs)[2:n])
	alphahat = 1/xihat
	y = switch(option,
		alpha = alphahat,
		xi = xihat,
		quantile = ordered * ((n * (1 - p))/k)^(-1/alphahat))
	ses = y/sqrt(k)
	if(is.na(end)) end = n
	x = trunc(seq(from = min(end, length(data)), to = start))
	y = y[x]
	ylabel = option
	yrange = range(y)
	if(ci && (option != "quantile")) {
       		qq = qnorm(1 - (1 - ci)/2)
       		u = y + ses[x] * qq
       		l = y - ses[x] * qq
       		ylabel = paste(ylabel, " (CI, p =", ci, ")", sep = "")
       		yrange = range(u, l)
	}
	if(option == "quantile") ylabel = paste("Quantile, p =", p)
	index = x
	if(reverse) index =  - x
	if(autoscale)
		plot(index, y, ylim = yrange, type = "l", xlab = "", ylab = "",
			axes = FALSE, ...)
	else plot(index, y, type = "l", xlab = "", ylab = "", axes = FALSE, ...)
	axis(1, at = index, lab = paste(x), tick = FALSE)
	axis(2)
	threshold = findThreshold(data, x)
	axis(3, at = index, lab = paste(format(signif(threshold, 3))),
		tick = FALSE)
	box()
	if(ci && (option != "quantile")) {
       	lines(index, u, lty = 2, col = 2)
       	lines(index, l, lty = 2, col = 2)}
	if(labels) {
       	title(xlab = "Order Statistics", ylab = ylabel)
       	mtext("Threshold", side = 3, line = 3)}
	
	# Return Value:
	invisible(list(x = index, y = y))
}


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


shaparmPlot = 
function (x, revert = FALSE, standardize = FALSE, tails = 0.01*(1:10), 
doplot = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE), 
which = c(TRUE, TRUE, TRUE), doprint = TRUE, both.tails = TRUE, 
xi.range = c(-0.5, 1.5), alpha.range = c(0, 10))
{	# A function written by D. Wuertz 
	
	# Description:
	#	Displays Pickands, Einmal-Decker-deHaan, and Hill
    #	estimators together with several plot variants.
	
	# FUNCTION:
	
	# Settings:
	select.doplot = which
	if (revert) x = -x
	if (standardize) x = (x-mean(x))/sqrt(var(x))
	ylim1 = xi.range
	ylim2 = alpha.range
  	z = rep(mean(ylim2), length(tails))
	ylim1 = xi.range
	ylim2 = alpha.range
  	p1 = p2 = h1 = h2 = d1 = d2 = m1 = m2 = rep(0,length(tails))
  	for ( i in (1:length(tails)) ) {
    	tail = tails[i]
    	
	# Printing/Plotting Staff:
	if(doprint) cat("Taildepth: ", tail, "\n")
	if(select.doplot[1]) {
    		xi = shaparmPickands (x, tail, ylim1, doplot=doplot[i], 
			both.tails, ) 
      		p1[i] = xi$xi[1]; p2[i] = xi$xi[3] }
	if(select.doplot[2]) { 
    		xi = shaparmHill (x, tail, ylim1, doplot=doplot[i], 
			both.tails) 
      		h1[i] = xi$xi[1]; h2[i] = xi$xi[3] }
	if(select.doplot[3]) {
    		xi = shaparmDEHaan (x, tail, ylim1, doplot=doplot[i], 
			both.tails)
      		d1[i] = xi$xi[1]; d2[i] = xi$xi[3] }      
	if(doprint) {
		cat("Pickands - Hill - DeckerEinmaalDeHaan: \n")
		print(c(p1[i], h1[i], d1[i]))
   		if (both.tails) print(c(p2[i], h2[i], d2[i]))} 
		cat("\n") }

	
	# Plot Pickands' Summary:
  		if(select.doplot[1]) { 
			plot (tails, z, type="n", xlab="tail depth", ylab="alpha",
				ylim=ylim2, main="Pickands Summary")
				y1 = 1/p1
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=2); lines(x1, y1, col=2)
  			if (both.tails) { 
				y1 = 1/p2
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=3); lines(x1, y1, col=3)} }
	
	# Plot Hill Summary:
  		if(select.doplot[2]) { 
			plot (tails, z, type="n", xlab="tail depth", ylab="alpha", 
				ylim=ylim2, main="Hill Summary")
				y1 = 1/h1
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=2); lines(x1, y1, col=2)
			if (both.tails) { 
  				y1 = 1/h2
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=3); lines(x1, y1, col=3)} }
	
	# Plot Deckers-Einmahl-deHaan Summary
  		if(select.doplot[3]) { 
			plot (tails, z, type="n", xlab="tail depth", ylab="alpha", 
				ylim=ylim2, 
    				main="Deckers-Einmahl-deHaan Summary")
				y1 = 1/d1
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=2); lines(x1, y1, col=2) 
  			if (both.tails) { 
				y1 = 1/d2
  				x1 = tails [y1>ylim2[1] & y1<ylim2[2]]
  				y1 = y1 [y1>ylim2[1] & y1<ylim2[2]]
  				points (x1, y1, col=3); lines(x1, y1, col=3)} }
	
	# Return Value:
		lower = list(pickands=p1, hill=h1, dehaan=d1)
		if (both.tails) {
			upper = list(pickands=p2, hill=h2, dehaan=d2)
			result = list(tails=tails, lower=lower, upper=upper) }
		else {
			result = list(tails=tails, lower=lower) }
	result
}


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


shaparmPickands = 
function (x, tail, yrange, doplot = TRUE, both.tails = TRUE, ...)     
{	# A function written by D. Wuertz
  	
	# FUNCTION:
	
	# Order Residuals:
  	ordered1 = rev(sort(abs(x[x < 0])))
  	if (both.tails) ordered2 = rev(sort(abs(x[x > 0])))
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2) 
  	ordered1 = ordered1[1:floor(tail*n1)]    
  	if (both.tails) ordered2 = ordered2[1:floor(tail*n2)]
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2)
	
  	# Pickands Estimate:	
	k1 = 1:(n1%/%4)
	if (both.tails) k2 = 1:(n2%/%4) 
  	pickands1 = log ((c(ordered1[k1])-c(ordered1[2*k1])) /
   		(c(ordered1[2*k1])-c(ordered1[4*k1]))) / log(2)
  	if (both.tails) pickands2 = log ((c(ordered2[k2])-c(ordered2[2*k2])) /
    	(c(ordered2[2*k2])-c(ordered2[4*k2]))) / log(2)
  	
    # Prepare Plot:
  		y1 = pickands1[pickands1 > yrange[1] & pickands1 < yrange[2]]
  		x1 = log10(1:length(pickands1))[pickands1 > yrange[1] & 
			pickands1 < yrange[2]]
  	if (both.tails) {
		y2 = pickands2[pickands2 > yrange[1] & pickands2 < yrange[2]]
  		x2 = log10(1:length(pickands2))[pickands2 > yrange[1] & 
			pickands2 < yrange[2]] }
  	if (doplot) { 
    		par(err=-1)
		plot (x1, y1, xlab="log scale", ylab="xi", ylim=yrange, 
			main="Pickands Estimator", type="n")  
		title(sub=paste("tail depth:", as.character(tail)))
    		lines(x1, y1, type="p", pch=2, col=2)
    		if (both.tails) lines(x2, y2, type="p", pch=6, col=3) }
  	
    # Calculate invers "xi":
  	my1 = mean(y1, na.rm = TRUE)
	if (both.tails) my2 = mean(y2, na.rm = TRUE)
  	sy1 = sqrt(var(y1, na.rm = TRUE))
	if (both.tails) sy2 = sqrt(var(y2, na.rm = TRUE))
  	
	# Plot:
	if (doplot) {
    	par(err=-1)
		lines(c(x1[1], x1[length(x1)]), c(my1,my1), type="l", 
		lty=1, col=2)
    	if (both.tails) lines(c(x2[1], x2[length(x2)]), c(my2,my2), 
		type="l", lty=1, col=3) }
  	
	# Return Result: 
  	result = list(xi=c(my1, sy1))   
	if (both.tails) result = list(xi=c(my1, sy1, my2, sy2))
	result
}
 

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

   
shaparmHill = 
function (x, tail, yrange, doplot = TRUE, both.tails = TRUE, ...)     
{ 	# A Function written by D. Wuertz
   	
	# ORDER RESIDUALS:
  	ordered1 = rev(sort(abs(x[x < 0])))
  	if (both.tails) ordered2 = rev(sort(abs(x[x > 0])))
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2) 
  	ordered1 = ordered1[1:floor(tail*n1)]    
  	if (both.tails) ordered2 = ordered2[1:floor(tail*n2)]
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2)
  	# HILLS ESTIMATE:
  	hills1 = c((cumsum(log(ordered1))/(1:n1)-log(ordered1))[2:n1])     
  	if (both.tails) hills2 = c((cumsum(log(ordered2))/(1:n2) -
		log(ordered2))[2:n2])
  	# PREPARE PLOT:
  		y1 = hills1[hills1 > yrange[1] & hills1 < yrange[2]]
  		x1 = log10(1:length(hills1))[hills1 > yrange[1] & 
			hills1 < yrange[2]]
  	if (both.tails) {
		y2 = hills2[hills2 > yrange[1] & hills2 < yrange[2]]
  		x2 = log10(1:length(hills2))[hills2 > yrange[1] & 
			hills2 < yrange[2]]}
  	if (doplot) {
    		par(err=-1)
		plot (x1, y1, xlab="log scale", ylab="xi", ylim=yrange, 
			main="Hill Estimator", type="n")
		title(sub=paste("tail depth:", as.character(tail)))
    		lines(x1, y1, type="p", pch=2, col=2)
    		if (both.tails) lines(x2, y2, type="p", pch=6, col=3) }
  	# CALCULATE INVERSE XI:
  	my1 = mean(y1, na.rm = TRUE)
	if (both.tails) my2 = mean(y2, na.rm = TRUE)
  	sy1 = sqrt(var(y1, na.rm = TRUE))
	if (both.tails) sy2 = sqrt(var(y2, na.rm = TRUE))
  	if (doplot) {
    	par(err=-1)
		lines(c(x1[1], x1[length(x1)]), c(my1,my1), type="l",
		lty=1, col=2)
    	if (both.tails) lines(c(x2[1], x2[length(x2)]), c(my2,my2), 
		type="l",lty=1, col=3) }
  	# Return Result:
  	result = list(xi=c(my1, sy1))   
	if (both.tails) result = list(xi=c(my1, sy1, my2, sy2))
	result       
}
       

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


shaparmDEHaan = 
function (x, tail, yrange, doplot = TRUE, both.tails = TRUE, ...)     
{	# A function written by D. Wuertz
  	
	# ORDER RESIDUALS:
  	ordered1 = rev(sort(abs(x[x < 0])))
  	if (both.tails) ordered2 = rev(sort(abs(x[x > 0])))
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2) 
  	ordered1 = ordered1[1:floor(tail*n1)]    
  	if (both.tails) ordered2 = ordered2[1:floor(tail*n2)]
  	n1 = length(ordered1)
	if (both.tails) n2 = length(ordered2)
  	# DECKERS-EINMAHL-deHAAN ESTIMATE:
  	ns0 = 1
  	n1m = n1-1; ns1 = ns0; ns1p = ns1+1
  	bod1 = c( cumsum(log(ordered1))[ns1:n1m]/(ns1:n1m) -
    		log(ordered1)[ns1p:n1] ) 
  	bid1 = c( cumsum((log(ordered1))^2)[ns1:n1m]/(ns1:n1m) -
    		2*cumsum(log(ordered1))[ns1:n1m]*log(ordered1)[ns1p:n1]/(ns1:n1m) +
    		((log(ordered1))^2)[ns1p:n1] )
  	dehaan1 = ( 1.0 + bod1 + ( 0.5 / (  bod1^2/bid1 - 1 ) ))
  	if (both.tails) {
	n2m = n2-1; ns2 = ns0; ns2p = ns2+1
  	bod2 = c( cumsum(log(ordered2))[ns2:n2m]/(ns2:n2m) -
    		log(ordered2)[ns2p:n2] ) 
  	bid2 = c(  cumsum((log(ordered2))^2)[ns2:n2m]/(ns2:n2m) -
    		2*cumsum(log(ordered2))[ns2:n2m]*log(ordered2)[ns2p:n2]/(ns2:n2m) +
    		((log(ordered2))^2)[ns2p:n2] )
  	dehaan2 = ( 1.0 + bod2 + ( 0.5 / (  bod2^2/bid2 - 1 ) )) }
  	# PREPARE PLOT:
  		y1 = dehaan1[dehaan1 > yrange[1] & dehaan1 < yrange[2]]
  		x1 = log10(1:length(dehaan1))[dehaan1 > yrange[1] & 
			dehaan1 < yrange[2]]
  	if (both.tails) {
		y2 = dehaan2[dehaan2 > yrange[1] & dehaan2 < yrange[2]]
  		x2 = log10(1:length(dehaan2))[dehaan2 > yrange[1] & 
			dehaan2 < yrange[2]] }
  	if (doplot) { 
    		par(err=-1)
		plot (x1, y1, xlab="log scale", ylab="xi", ylim=yrange,
       			main="Deckers - Einmahl - de Haan Estimator", type="n")
		title(sub=paste("tail depth:", as.character(tail)))
    		lines(x1, y1, type="p", pch=2, col=2)
    		if (both.tails) lines(x2, y2, type="p", pch=6, col=3) }
  	# CALCULATE INVERSE XI:
  	my1 = mean(y1, na.rm = TRUE)
	if (both.tails) my2 = mean(y2, na.rm = TRUE)
  	sy1 = sqrt(var(y1, na.rm = TRUE))
	if (both.tails) sy2 = sqrt(var(y2, na.rm = TRUE))
  	if (doplot) {
    	par(err=-1)
		lines(c(x1[1], x1[length(x1)]), c(my1,my1), type="l", 
			lty=1, col=2)
    	if (both.tails) lines(c(x2[1], x2[length(x2)]), c(my2, my2), 
		type = "l", lty = 1, col = 3) }
  	# Return Result:
  	result = list(xi = c(my1, sy1))   
	if (both.tails) result = list(xi=c(my1, sy1, my2, sy2))
	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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


# ##############################################################################
# FUNCTION:             GPD DISTRIBUTION FAMILY:
# dgpd					 Density for the Generalized Pareto DF [USE FROM EVIS]
#  pgpd					  Probability for the Generalized Pareto DF
#  qgpd					  Quantiles for the Generalized Pareto DF
#  rgpd					  Random variates for the Generalized Pareto DF
################################################################################


dgpd = 
function(x, xi = 1, mu = 0, beta = 1)
{	# A function written by Diethelm Wuertz

	# FUNCTION:
	
	# Density:
	y = (x - mu)
	if (xi == 0) {
		d = (1-exp(-y/beta))/beta }
	else {
		d = 1/beta * (1 + (xi*y)/beta)^((-1/xi) - 1) }	
	
	d[y < 0] = 0
	if (xi < 0) d[y > (-1/xi)] = 0
	
	# Return Value:
	d
}


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


pgpd = 
function(q, xi = 1, mu = 0, beta = 1)
{	# A function written by Diethelm Wuertz

	# FUNCTION:
	
	# Probability:
	y = (q - mu)
	if (xi == 0) {
		p = y/beta + exp(-y/beta) -1 }
	else {
		p = (1 - (1 + (xi*y)/beta)^(-1/xi)) }	
	
	p[y < 0] = 0
	if (xi < 0) p[y > (-1/xi)] = 1
	
	# Return Value:
	p
}


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


qgpd = 
function(p, xi = 1, mu = 0, beta = 1)
{	# A function written by Diethelm Wuertz

	# FUNCTION:
	
	# Quantiles:
	if (xi == 0) 
		q = mu - beta*log(1-p)
	else
		q = mu + (beta/xi) * ((1 - p)^( - xi) - 1)
	
	# Return Value:
	q
}


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


rgpd = 
function(n, xi = 1, mu = 0, beta = 1)
{	# A function written by Diethelm Wuertz

	# FUNCTION:
	
	# Random variates:
	rvs = mu + (beta/xi) * ((1 - runif(n))^( - xi) - 1)
	
	# Return Value:
	rvs
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:					GPD MODELLING FROM EVIS:
#  gpdSim					 Simulates GPD rvs
#  gpdFit					 Fits GPD Distribution
#   print.gpd 		  		  Print Method for object of class "gpd"
#   plot.gpd 		  		  Plot Method for object of class "gpd"
#   summary.gpd       		  Summary Method for object of class "gpd"
# FUNCTION:                 ADDITIONAL PLOTS:
#  gpdtailPlot		         Plots Tail Estimate From GPD Model
#  gpdquantPlot	  		     Plots of GPD Tail Estimate of a High Quantile
#  gpdshapePlot      		 Plots for GPD Shape Parameter
#  gpdqPlot          		 Adds Quantile Estimates to plot.gpd
#  gpdsfallPlot      		 Adds Expected Shortfall Estimates to a GPD Plot
# FUNCTION:                 ADDITIONAL FUNCTION:
#   gpdriskmeasures   		 Calculates Quantiles and Expected Shortfalls
################################################################################


gpdSim = 
function(model = list(shape = 0.25, location = 0, scale = 1), n = 1000)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Generates random variates from a GPD distribution
	
	
	# FUNCTION:
	
	# Simulate:
	rgpd(n = n, xi = model$shape, mu = model$location, beta = model$scale)	
}


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


gpdFit =
function(x, threshold = NA, nextremes = NA, type = c("mle", "pwm"),
information = c("observed", "expected"), ...)
{   # A function implemented by Diethelm Wuertz

	# Description:
	#	Returns an object of class `"gpd"' 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:
	
    # Make the fit:
    call = match.call()
    type = type[1]
    # if (is.na(threshold) & is.na(nextremes)) threshold = min(x)
    if (type == "mle") type = "ml"
    fitted = gpd(data = x, threshold = threshold, nextremes = nextremes, 
    	method = type, information = information, ...) 
    	
	# Residuals:
    xi = fitted$par.ests["xi"]
    beta = fitted$par.ests["beta"]
	excess = as.numeric(fitted$data) - fitted$threshold
   	residuals = log(1 + (xi * excess)/beta)/xi
	
	# Make Unique:
	fit = list()
	fit$fit =  fitted
	fit$call = call
	fit$type = c("gpd", type[1])
	fit$par.ests = fitted$par.ests
	fit$par.ses = fitted$par.ses
	fit$residuals = residuals
	fit$fitted.values = fitted$data - residuals
	fit$llh = fitted$nllh.final
	fit$converged = fitted$converged
	
	# Return Value:
	class(fit) = "gpdFit"	
	fit
}


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


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

	# Description:
	#	Print Method for an object of class 'gpdFit'
	
	# FUNCTION:
	
	# Function Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "") 
			
	# Estimation Type:
	cat("\nEstimation Type:", x$type, "\n")	
	
	# Estimated Parameters:
	cat("\nEstimated Parameters:\n")
	print(x$par.ests)
	cat("\n")
	
	# Return Value:
	invisible(x)
}


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


plot.gpdFit =
function(x, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for objects of class 'gpdFit'

    # FUNCTION:
    
    # Plot Functions:
    plot.1 <<- function(x, ...) {
	    fit = x
    	data = fit$fit$data
    	xi = fit$par.ests[1]
    	beta = fit$par.est[2]
    	threshold = fit$fit$threshold
    	optlog = NA
    	extend = 1.5
    	labels = TRUE
    	# Start:
    	plotmin = threshold
    	if (extend <= 1) stop("extend must be > 1")
	  	plotmax = max(data) * extend
     	xx = seq(from = 0, to = 1, length = 1000)
     	z = qgpd(xx, xi, threshold, beta)
     	z = pmax(pmin(z, plotmax), plotmin)
       	ypoints = ppoints(sort(data))
       	y = pgpd(z, xi, threshold, beta) 
       	type = "eplot"
	    if (!is.na(optlog)) alog = optlog
       	else alog = "x"
	    if (alog == "xy") stop("Double log does not make much sense")
	    yylab = "Fu(x-u)"
       	shape = xi
	    scale = beta
	    location = threshold
	    plot(sort(data), ypoints, xlim = range(plotmin, plotmax),
        	ylim = range(ypoints, y, na.rm = TRUE), xlab = "",
     		ylab = "", log = alog, axes = TRUE, 
        	main = "Excess Distribution", ...)
		lines(z[y >= 0], y[y >= 0])    
	    xxlab = "x"
		if (alog == "x" || alog == "xy" || alog == "yx")
		    xxlab = paste(xxlab, "(on log scale)")
	    if (alog == "xy" || alog == "yx" || alog == "y")
		    yylab = paste(yylab, "(on log scale)")
		title(xlab = xxlab, ylab = yylab) } 
    plot.2 <<- function(x, ...) {
	    fit = x
    	data = fit$fit$data
    	xi = fit$par.ests[1]
    	beta = fit$par.est[2]
    	threshold = fit$fit$threshold
    	optlog = NA
    	extend = 1.5 #; if(extend <= 1) stop("extend must be > 1")
    	labels = TRUE
    	# Start:
    	plotmin = threshold
     	if (extend <= 1) stop("extend must be > 1")
		plotmax = max(data) * extend
      	xx = seq(from = 0, to = 1, length = 1000)
     	z = qgpd(xx, xi, threshold, beta)
     	z = pmax(pmin(z, plotmax), plotmin)
     	ypoints = ppoints(sort(data))
     	y = pgpd(z, xi, threshold, beta)
     	type = "tail"
		if (!is.na(optlog)) alog = optlog
	   	else alog = "xy"
	   	prob = fit$fit$p.less.thresh
	   	ypoints = (1 - prob) * (1 - ypoints)
	   	y = (1 - prob) * (1 - y)
	  	yylab = "1-F(x)"
	   	shape = xi
	  	scale = beta * (1 - prob)^xi
	  	location = threshold - (scale * ((1 - prob)^( - xi) - 1))/xi
	  	plot(sort(data), ypoints, xlim = range(plotmin, plotmax),
        	ylim = range(ypoints, y, na.rm = TRUE), xlab = "",
     		ylab = "", log = alog, axes = TRUE, 
        	main = "Tail of Underlying Distribution", ...)
		lines(z[y >= 0], y[y >= 0])    
	    xxlab = "x"
		if (alog == "x" || alog == "xy" || alog == "yx")
		    xxlab = paste(xxlab, "(on log scale)")
	    if (alog == "xy" || alog == "yx" || alog == "y")
		    yylab = paste(yylab, "(on log scale)")
		title(xlab = xxlab, ylab = yylab) }
    plot.3 <<- function(x, ...) {
	    res = x$residuals
    	plot(res, 
    		ylab = "Residuals", 
    		xlab = "Ordering", 
      	    main = "Scatterplot of Residuals", ...)
	    lines(lowess(1:length(res), res)) }      
	plot.4 <<- function(x, ...) {
	    qplot(x$residuals, 
	    	main = "QQ-Plot of Residuals", ...) }
	    	
	# Plot:
	interactivePlot(
		x = x,
		choices = c(
			"Excess Distribution",
			"Tail of Underlying Distribution",
			"Scatterplot of Residuals", 
			"QQ-Plot of Residuals"),
		plotFUN = c(
			"plot.1", 
			"plot.2", 
			"plot.3",
			"plot.4"),
		which = which)
	            	
	# Return Value:
	invisible(x)
}


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


summary.gpdFit = 
function(object, doplot = TRUE, which = "all", ...) 
{	# A function written by Diethelm Wuertz

	# Description:
	#	Summary method for objects of class "gpdFit"
	
	# FUNCTION:
    
    # Print:
    print(object, ...)
	
	# Summary:
	# For MLE print additionally:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$conv, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot(object, which = which, ...)
	cat("\n")
	
	# Return Value:
	invisible(object)
}


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


gpdtailPlot =
function(fit, optlog = NA, extend = 1.5, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plots Tail Estimate From GPD Model
	
	# FUNCTION:
	
	# Return Value:
	tailplot(x = fit, optlog = optlog, extend = extend, labels = labels, ...)
}


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


gpdquantPlot =
function(data, p = 0.99, models = 30, start = 15, end = 500,
reverse = TRUE, ci = 0.95, autoscale = TRUE, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plots of GPD Tail Estimate of a High Quantile
	
	# FUNCTION:

	# Return Value:
	quant(data = data, p = p, models = models, start = start, end = end,
		reverse = reverse, ci = ci, auto.scale = autoscale, labels = labels,
		...)
}


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


gpdshapePlot =
function(data, models = 30, start = 15, end = 500, reverse = TRUE, 
ci = 0.95, autoscale = TRUE, labels = TRUE, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plots for GPD Shape Parameter
	
	# FUNCTION:

	# Return Value:
	shape(data = data, models = models, start = start, end = end, 
		reverse = reverse, ci = ci, auto.scale = autoscale, 
		labels = labels, ...)
}


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


gpdqPlot =
function(x, pp, ci.type = c("likelihood", "wald"), ci.p = 0.95, like.num = 50)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Adds Quantile Estimates to plot.gpd
	
	# FUNCTION:

	# Return Value:
	gpd.q(x = x, pp = pp, ci.type = ci.type, ci.p = ci.p, like.num = like.num)
}


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


gpdsfallPlot =
function(x, pp, ci.p = 0.95, like.num = 50)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Adds Expected Shortfall Estimates to a GPD Plot	
	
	# FUNCTION:

	# Return Value:
	gpd.sfall(x = x, pp = pp, ci.p = ci.p, like.num = like.num)
}


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


gpdriskmeasures = 
function(x, plevels)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Calculates Quantiles and Expected Shortfalls	
	
	# FUNCTION:
	
	# Return Value:
	riskmeasures(x = x, p = plevels)
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:					GPD MODELLING FROM ISMEV:
#  gpdglmFit				 Fits GPD Distribution
#   print.gpdglmFit		      Print Method for object of class "gpdglm"
#   plot.gpdglmFit		      Plot Method for object of class "gpdglm"
#   summary.gpdglmFit		  Summary Method for object of class "gevglm"
# FUNCTION:                 ADDITIONAL PLOTS:
#  gpdglmprofPlot			 Profile Log-likelihoods for Stationary GPD Models
#  gpdglmprofxiPlot	         Profile Log-likelihoods for Stationary GPD Models
################################################################################


gpdglmFit =
function(x, threshold = min(x), npy = 365, y = NULL, sigl = NULL, shl = NULL, 
siglink = identity, shlink = identity, show = FALSE, method = "Nelder-Mead", 
maxit = 10000, ...)
{   # A function implemented by Diethelm Wuertz

	# Description:
	
	# FUNCTION:
    
    # Function Call:
    call = match.call()
    # Fit Parameters:
    fitted = gpd.fit(xdat = x, threshold = threshold, npy = npy, ydat = y, 
    	sigl = sigl, shl = shl, siglink = siglink, shlink = shlink, 
    	show = show, method = method, maxit = maxit, ...)
    # Add names attribute:
    names(fitted$se)  = names(fitted$mle) = c("sigma", "mle")
    	
    # Add:
    fit = list()
    fit$fit = fitted
	fit$call = call
	fit$type = c("gpdglm", "mle")
	fit$par.ests = fitted$mle
	fit$par.ses = fitted$se
	fit$residuals = fitted$residuals
	fit$fitted.values = x - fitted$residuals
	fit$llh = fitted$nllh
	fit$converged = fitted$conv
	
	# Return Value:
	class(fit) = "gpdglmFit"
    fit
}


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


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

	# Description:
	#	Print Method for an object of class 'gpdglmFit'
	
	# FUNCTION:

	# Print Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "") 	
	
	if (fit$fit$trans) {	
		# Still to do:
		print.default(x) }	
	else {
		# Estimation Type:
		cat("\nEstimation Type:", x$type, "\n")	
		# Fitted Parameters:
		cat("\nEstimated Parameters:\n")
		print(x$par.ests)
		cat("\n") }
		
	# Return Value:
	invisible(x)
}


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


plot.gpdglmFit =
function(x, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Print Method for an object of class 'gpdglmFit'
	
	# FUNCTION:
	
	# Plot Functions:
	if(x$fit$trans) {
       	plot.1 <<- function(x, ...) {
	       	n = length(x$data)
	       	plot( 
	       	 	x = (1:n)/(n + 1),
	       		y = 1 - exp( - sort(x$data)), 
	       		xlab = "Empirical", 
       			ylab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Probability Plot") }
       	plot.2 <<- function(x, ...) {
       		n = length(x$data)
	       	plot( 
       			x = - log( 1 - (1:n)/(n+1) ), 
       			y = sort(x$data), 
       			ylab = "Empirical", 
       			xlab = "Model")
       		abline(0, 1, col = 4)
       		title("Residual Quantile Plot (Exptl. Scale)") } }
    else {
       	plot.1 <<- function(x, ...) {
	       	gpd.pp(x$mle, x$threshold, x$data) }
       	plot.2 <<- function(x, ...) {
	       	gpd.qq(x$mle, x$threshold, x$data) }
       	plot.3 <<- function(x, ...) {
	       	gpd.rl(x$mle, x$threshold, x$rate, x$n, x$npy, 
       			x$cov, x$data, x$xdata) }
       	plot.4 <<- function(x, ...) {
	       	gpd.his(x$mle, x$threshold, x$data) } }
	       	       	
	# Plot:
	if (fit$fit$trans) {
		interactivePlot(
			x = x$fit,
			choices = c(
				"Excess Distribution",
				"QQ-Plot of Residuals"),
			plotFUN = c(
				"plot.1", 
				"plot.2"),
			which = which) }
	else {
		interactivePlot(
			x = x$fit,
			choices = c(
				"Probability Plot",
				"Quantile Plot",
				"Return Level Plot", 
				"Histogram Plot"),
			plotFUN = c(
				"plot.1", 
				"plot.2", 
				"plot.3",
				"plot.4"),
			which = which) }
			
	# Return Value:
    invisible(x)
}
	

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


summary.gpdglmFit =
function(object, doplot = TRUE, which = "all", ...) 
{	# A function written by Diethelm Wuertz

	# Description:
	#	Summary Method for an object of class 'gpdglmFit'
	
	# FUNCTION:
    
    # Print:
    print(object, ...)
	
	# Summary:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$converged, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot(object, which = which, ...)
	cat("\n")
	
	# Return Value:
	invisible(object)
}


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


gpdglmprofPlot =
function(fit, m, xlow, xup, conf = 0.95, nint = 100, ...)
{	# A function implemented by Diethelm Wuertz

  	# Description:
  	#	Profile Log-likelihoods for Stationary GPD Models
  	
  	# FUNCTION:
  	
  	# Compute:
	gpd.prof(z = fit$fit, m = m, xlow = xlow, xup = xup , conf = conf, 
		nint = nint) 
}


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


gpdglmprofxiPlot =
function(fit, xlow, xup, conf = 0.95, nint = 100, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
  	#	Profile Log-likelihoods for Stationary GPD Models

	# FUNCTION:
	
	# Compute:
	gpd.profxi(z = fit$fit, xlow = xlow, xup = xup, conf = conf, 
			nint = nint, ...) 
}


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


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:				   POT MODELLING FROM EVIS:
#  potSim   				Peaks over a threshold from arbitrary series
#  potFit					Fits with POT method
#   print.potFit 		  	 Print Method for object of class "potFit"
#   plot.potFit 		  	 Print Method for object of class "potFit"
#   summary.potFit           Summary Method for object of class "potFit"
# REQUIRES:
#  ts						Package ts (is preloaded)
################################################################################


potSim = 
function(x, threshold, nextremes = NA, run = NA)
{	# A function implemented by Diethelm Wuertz

	# Description:
	# 	Generates from an arbitray rvs sequence a series with
	# 	the peaks over a threshold
	
	# Settings:
	data = x
	n = length(as.numeric(data))
    times = attributes(data)$times
    if (is.null(times)) {
        times = 1:n
        attributes(data)$times = times
        start = 1
        end = n
        span = end - start}
    else {
        start = times[1]
        end = times[n]
        span = as.numeric(difftime(as.POSIXlt(times)[n], as.POSIXlt(times)[1], 
            units = "days"))}
            
    if (is.na(nextremes) && is.na(threshold)) 
        stop("Enter either a threshold or the number of upper extremes")
    if (!is.na(nextremes) && !is.na(threshold)) 
        stop("Enter EITHER a threshold or the number of upper extremes")
    if (!is.na(nextremes)) 
        threshold = findthresh(as.numeric(data), nextremes)
    if (threshold > 10) {
        factor = 10^(floor(log10(threshold)))
        cat(paste("If singularity problems occur divide data", 
            "by a factor, perhaps", factor, "\n")) }
            
    exceedances.its = structure(data[data > threshold], times = times[data > 
        threshold])
    n.exceed = length(as.numeric(exceedances.its))
    p.less.thresh = 1 - n.exceed/n
    if (!is.na(run)) {
        exceedances.its = decluster(exceedances.its, run, picture)
        n.exceed = length(exceedances.its) }
    intensity = n.exceed/span
    exceedances = as.numeric(exceedances.its)
    
    
    # Return Value:
    exceedances
}


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


potFit =
function(x, threshold = NA, nextremes = NA, run = NA, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Parameter Estimation for the POT model.

	# FUNCTION:
	
	# Call pot() from evir:	
	call = match.call()
	fitted = pot(data = x, threshold = threshold, nextremes = nextremes, 
		run = run, picture = FALSE, ...) 
		
	# Compute Residuals:
	xi = fitted$par.ests[1]
    beta = fitted$par.ests[4]
    threshold = fitted$threshold 
    fitted$residuals = 
    	as.vector(log(1 + (xi * (fitted$data - threshold))/beta)/xi)
    
    # Gaps:
	x = fitted
    x$rawdata = x$data
    n = length(as.numeric(x$rawdata))
    x$times = attributes(x$rawdata)$times
    if (is.character(x$times) || inherits(x$times, "POSIXt") || 
    	inherits(x$times, "date") || inherits(x$times, "dates")) {
        x$times = as.POSIXlt(x$times)
        x$gaps = as.numeric(difftime(x$times[2:n], x$times[1:(n - 1)], 
        	units = "days")) * x$intensity }
    else {
	    x$times = 1:n
    	x$gaps = as.numeric(diff(x$times)) * x$intensity }
    fitted$times = x$times
    fitted$rawdata = x$rawdata
    fitted$gaps = x$gaps
    
	# Add:
	fit = list()
	fit$fit = fitted
	fit$call = call
	fit$type = c("pot", "mle")
	fit$par.ests = fitted$par.ests
	fit$par.ses = fitted$par.ses
	fit$residuals = fitted$residuals
	fit$fitted.values = fitted$data - fitted$residuals
	fit$llh = fitted$nllh.final
	fit$converged = fitted$converged 		
	
	# Return Value:
	class(fit) = "potFit"
	fit
}


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


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

	# Description:
	#	Print Method for object of class "potFit"
	
	# FUNCTION:

	# Function Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "")
	 		
	# Estimation Type:
	cat("\nEstimation Type:", x$type, "\n")	
	
	# Estimated Parameters:
	cat("\nEstimated Parameters:\n")
	print(x$par.ests)
	cat("\n")
	
	# Decluster Run Length:
	if (!is.na(fit$fit$run))
	cat("\nDecluster Runlength:", x$fit$run, "\n")
	
	# Return Value:
	invisible(x)
}


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


plot.potFit =
function(x, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for objects of class "potFit".

	# FUNCTION:
             
    # Plot functions:   
    plot.1 <<- function(x, ...) {
        plot(x$times, x$rawdata , type = "h", 
           	main = "Point Process of Exceedances", ...) }
	plot.2 <<- function(x, ...) {
	    plot(x$gaps, ylab = "Gaps", xlab = "Ordering", 
        	main = "Scatterplot of Gaps", ...)
     	lines(lowess(1:length(x$gaps), x$gaps)) } 
	plot.3 <<- function(x, ...) {
	    qplot(x$gaps, 
        	main = "QQ-Plot of Gaps", ...) }           
   	plot.4 <<- function(x, ...) {
	 	acf(x$gaps, lag.max=20, 
         	main = "ACF of Gaps", ...) }
	plot.5 <<- function(x, ...) {
	  	plot(x$residuals, ylab = "Residuals", xlab = "Ordering", 
       		main = "Scatterplot of Residuals", ...)
    	lines(lowess(1:length(x$residuals), x$residuals)) }
	plot.6 <<- function (x, ...) {
	  	qplot(x$residuals, 
        	main = "QQ-Plot of Residuals", ...) }
	plot.7 <<- function (x, ...) { 
	  	acf(x$residuals, lag.max = 20, 
       		main = "ACF of Residuals", ...) }
    fit <<- fit; plot.8 <<- function (x, ...) { 
	    if (which == "ask") {
		    plot.gpd(x)
		    plot.potFit(fit, which = "ask") } }

    # Plot:
	interactivePlot(
		x = x$fit,
		choices = c(
			"Point Process of Exceedances", 
    		"Scatterplot of Gaps", 
        	"QQ-Plot of Gaps", 
        	"ACF of Gaps", 
        	"Scatterplot of Residuals", 
        	"QQ-Plot of Residuals", 
        	"ACF of Residuals",
        	"GOTO GPD Plots"),
		plotFUN = c(
			"plot.1", 
			"plot.2", 
			"plot.3",
			"plot.4", 
			"plot.5", 
			"plot.6",
			"plot.7",
			"plot.8"),
		which = which)
	            	
	# Return Value:
	invisible(x)
}


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


summary.potFit =
function(object, doplot = TRUE, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Summary Method for object of class "potFit"
	
	# FUNCTION:

	# Print:
	print(object, ...)
	
	# Summary:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$converged, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot.potFit(object, which = which, ...)
	cat("\n")
    
    # Return Value:
    invisible(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

# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:					POINT PROCESS MODELLING FROM ISMEV:
#  ppFit					 Fits Point Process Model
#   print.ppFit			      Print Method for object of class "ppFit"
#   plot.ppFit			      Plot Method for object of class "ppFit"
#   summary.ppFit			  Summary Method for object of class "ppFit"
################################################################################


ppFit =
function(x, threshold, npy = 365, y = NULL, mul = NULL, sigl = NULL, 
shl = NULL, mulink = identity, siglink = identity, shlink = identity, 
method = "Nelder-Mead", maxit = 10000, ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	
	# FUNCTION:
    
    # Function Call:
    call = match.call()
    # Fit Parameters:
	fitted = pp.fit(xdat = x, threshold = threshold, npy = npy, ydat = y, 
		mul = mul, sigl = sigl, shl = shl, mulink = mulink, siglink = siglink, 
		shlink = shlink, show = FALSE, method = method, maxit = maxit, ...)	
	names(fitted$mle) =names(fitted$se) = c("xi", "sigma", "mu")
	# Compute Residuals:
	residuals = NA
	# cat("\nVariance Covariance Matrix:\n")
	# covar = fit$cov
	# covar[1,1] = fit$cov[3,3]
	# covar[3,3] = fit$cov[1,1]
	# covar[1,2] = covar[2,1] = fit$cov[2,3]
	# covar[2,3] = covar[3,2] = fit$cov[1,2] 
	# print(covar)	
		
	# Add:
	fit= list()
	fit$fit = fitted
	fit$call = call
	fit$type = c("pp", "mle")
	fit$par.ests = fitted$mle
	fit$par.ses = fitted$se
	fit$residuals = residuals
	fit$fitted.values = x - residuals
	fit$llh = fitted$nllh
	fit$converged = fitted$conv	
	
	# Return Value:
	class(fit) = "ppFit"
    fit
}


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


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

	# Description:
	#	Print method for an object of class 'ppFit'
	
	# FUNCTION:
	
	# Print Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "") 	
	
	if (x$fit$trans) {	
		# Still to do:
		print.default(x) }	
	else {	
		# Parameters - We use the same order as in gevFit:
		cat("\nParameter Estimate:\n")
		print(x$par.ests) }
		
	# Return Value:
	invisible(x)
}
	

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


plot.ppFit =
function(x, which = "ask", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for an object of class 'ppFit'
	
	# FUNCTION:
	
	# Plot Functions:
    if (x$fit$trans) {   
	    plot.1 <<- function(x, ...) {
		    n <- length(x$data)
	    	xx <- (1:n)/(n + 1)
	    	plot(xx, sort(x$data), xlab = "empirical", ylab = "model")
	        abline(0, 1, col = 3)
	        title("Residual Probability Plot") }  
        plot.2 <<- function(x, ...) {
	        n <- length(x$data)
	    	xx <- (1:n)/(n + 1)
	    	plot(-log(1 - xx), -log(1 - sort(x$data)), ylab = "empirical", 
	            xlab = "model")
	        abline(0, 1, col = 3)
	        title("Residual quantile Plot (Exptl. Scale)") } }
    else {
	    plot.1 <<- function(x, ...) {
		    # Probability Plot:
        	pp.pp(x$mle, x$threshold, x$npy, x$data) }      
        plot.2 <<- function(x, ...) {
	        # Quantile Plot:
        	pp.qq(x$mle, x$threshold, x$npy, x$data) } }
	
	# Plot:
	if (x$fit$trans) {
		interactivePlot(
			x = x$fit,
			choices = c(
				"Residual Probability Plot",
				"Residual Quantile Plot"),
			plotFUN = c(
				"plot.1", 
				"plot.2"),
			which = which) }
	else {
		interactivePlot(
			x = x$fit,
			choices = c(
				"Probability Plot",
				"Quantile Plot"),
			plotFUN = c(
				"plot.1", 
				"plot.2"),
			which = which) }
			
    # Return Value:
    invisible(x)
}

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


summary.ppFit =
function(object, doplot = TRUE, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Summary method for an object of class 'ppFit'
	
	# FUNCTION:
	
	# Print:
	print(object, ...)
	
	# Summary:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$converged, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot.ppFit(object, which = which, ...)
	cat("\n")
	
	# Return Value:
	invisible(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

# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:					R-LARGEST ORDER MODELLING FROM ISMEV:
#  rlargFit				     Fits r-largest Order Statistic Model
#   print.rlargFit			  Print Method for object of class "rlargFit"
#   plot.rlargFit			  Plot Method for object of class "rlargFit"
#   summary.rlargFit		  Summary Method for object of class "rlargFit"
################################################################################


rlargFit =
function(x, r = dim(x)[2], y = NULL, mul = NULL, sigl = NULL, shl = NULL, 
mulink = identity, siglink = identity, shlink = identity, method = 
"Nelder-Mead", maxit = 10000, ...)
{   # A function implemented by Diethelm Wuertz

	# Description:
	#	Maximum-likelihood fitting for the order statistic model,
    #	including generalized linear modelling of each parameter.
    
	# FUNCTION:
   
    # Function Call:
    call = match.call()
    
    # Fit Parameters
    fitted = rlarg.fit(xdat = x, r = r, ydat = y, mul = mul, sigl = sigl, 
		shl = shl, mulink = mulink, siglink = siglink, shlink = shlink, 
		show = FALSE, method = method, maxit = maxit, ...)
    
	# Further Values:
	mle = rev(fitted$mle)
	se = rev(fitted$se)
	names(mle) = names(se) = c("xi", "sigma", "mu")
	covar = fitted$cov
	covar[1,1] = fitted$cov[3,3]
	covar[3,3] = fitted$cov[1,1]
	covar[1,2] = covar[2,1] = fitted$cov[2,3]
	covar[2,3] = covar[3,2] = fitted$cov[1,2]
	
	# Make Unique:
	fit = list()
	fit$fit = fitted	
	fit$call = call
	fit$type = c("mle", "rlarg")
	fit$par.ests = mle
	fit$par.ses = se
	fit$residuals = as.matrix(fitted$data)
	fit$fitted.values = as.matrix(x) - fit$residuals
	fit$cov = covar
	fit$llh = fitted$nllh 
	fit$converged = fitted$conv	
	
	# Return Value:
	class(fit) = "rlargFit"
    fit
}


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


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

	# Description:
	#	Print Method for object of class "rlargFit"
	
	# Notes:
	#	The ismev package has no print method. It uses the command
	#   > summary.rlargFit(fit = fit, details = FALSE, doplot = FALSE, ...) 

	# FUNCTION:
	
	# Function Call:
	cat("\nCall:\n")
	cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "") 
			
	# Estimation Type:
	cat("\nEstimation Type:", x$type, "\n")	
	
	# Estimated summaryParameters:
	cat("\nEstimated Parameters:\n")
	print(x$par.ests)
	cat("\n")
	
	# Return Value:
	invisible(x)
}


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


plot.rlargFit = 
function(x, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Plot method for objects of class "rlargFit".
	
	# FUNCTION:
	
	# Plot Functions:
	if (x$fit$trans) {
		# Non-Stationary:
		plot.1 <<-  function(x, ...) {
		    for (i in 1:z$r) {
				# Probability and Quantile Plots:
				rlarg.pp(c(0, 1, 0), x$data[, 1:x$r], i)
				rlarg.qq(c(0, 1, 0), x$data[, 1:x$r], i) } } }
	else {	
        # Stationary - GEV Plots:
        plot.1 <<- function(x, ...) {
	        gev.pp(x$mle, x$data[, 1]) }
        plot.2 <<- function(x, ...) {
	        gev.qq(x$mle, x$data[, 1]) }
        plot.3 <<- function(x, ...) {
	        gev.rl(x$mle, x$cov, x$data[, 1]) }
        plot.4 <<- function(x, ...) {
	        gev.his(x$mle, x$data[, 1]) }  
		fit <<- fit; plot.5 <<- function(x, ...) {
			par(ask = TRUE)
			for (i in 1:fit$fit$r) {
				# Probability and Quantile Plots:
				rlarg.pp(x$mle, x$data, i)
				rlarg.qq(x$mle, x$data, i) } 
			par(ask = FALSE) } } 
						
	# Plot:
	if (x$fit$trans) {
		interactivePlot(
			x = x$fit,
			choices = c(
				"Probability Plot",
				"Quantile Plot"),
			plotFUN = c(
				"plot.1", 
				"plot.2"),
			which = which) }
	else {
		interactivePlot(
			x = x$fit,
			choices = c(
				"GEV Probability Plot",
				"GEV Quantile Plot",
				"GEV Return Level Plot",
				"GEV Histogram Plot",
				"R-Largest PP and QQ Plots"),
			plotFUN = c(
				"plot.1", 
				"plot.2",
				"plot.3", 
				"plot.4",
				"plot.5"),
			which = which) }
		
	# Return Value:
	invisible(x)
}


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


summary.rlargFit =
function(object, doplot = TRUE, which = "all", ...)
{	# A function implemented by Diethelm Wuertz

	# Description:
	#	Summary Method for object of class "rlargFit".
	
	# FUNCTION:
	
	# Print:
	print(object, ...)
	
	# Summary:
	cat("\nStandard Deviations:\n"); print(object$par.ses)
	cat("\nLog-Likelihood Value: ", object$llh)
	cat("\nType of Convergence:  ", object$converged, "\n") 
	cat("\n")
	
	# Plot:
	if (doplot) plot(object, which = which, ...)
	cat("\n")

	# Return Value:
	invisible(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

# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


################################################################################
# FUNCTION:                 EXTREMAL INDEX:
#  exindexesPlot             Calculate and Plot Theta(1,2,3)
#  exindexPlot               Calculate Theta(1,2) and Plot Theta(1)
################################################################################


exindexesPlot = 
function (x, block=20, quantiles = seq(0.990,0.999,0.001), doplot = TRUE, ...)
{   # A function written by D. Wuertz
    
    # Description:
    # 	Calculates and Plots Theta(1,2,3)
    
    # FUNCTION:
    
    # Settings:
    main = "Extremal Index"
    doprint = FALSE
    
    # Block Size:
    blocklength = block # argument renamed
    
    # Note, in finance the x's should be residuals
    resid = x 
    
    # Extremal Index - Theta_1, Theta_2 and Theta_3
    k = floor(length(resid)/blocklength) # Number of blocks
    n = k*blocklength # Number of data points 
    
    # Now organize your residuels:
    # 1) truncate the rest of the time series,
    # 2) arrange them in matrix form,
    # 3) sort them in reverse order, ie. from high (pos) to low (neg)
    resid1 = resid[1:(k*blocklength)]
    resid1 = matrix(resid1, ncol = blocklength, byrow = TRUE)
    ordered1 = sort(resid1)
    
    # Threshold values associated to quantiles:
    z0 = ordered1[floor(quantiles*length(resid1))]
    
    # Printing:
    if (doprint) {print(z0); print(n); print(k) }
    
    # Presettings:
    theta1 = theta2 = theta3 = rep(0, times = length(quantiles))
    
    # Calculate Extremal Imdex:
    run = 0
    for ( z in z0 ) {
        run = run + 1
        # N - number of exceedences:
        N = length(resid1[resid1>z])
        # K - number of blocks with exceedences:
        K = sum(sign(apply(resid1,1,max)-z)+1)/2
        if (K/k < 1) theta1[run] = (k/n) * log(1-K/k) / log(1-N/n)
          else theta1[run] = NA 
        theta2[run] = K/N
        x = 1:n
        xx = diff(x[resid1 > z])
        xx = xx[xx>blocklength]
        theta3[run] = length(xx)/N
        # Printing: 
        if (doprint) {
            print(c(N, K, quantiles[run], z)) 
            print(c(theta1[run], theta2[run], theta3[run]))} }
    
    # Plotting:
    if (doplot) {
        plot(quantiles, theta1, 
            xlim = c(quantiles[1], quantiles[length(quantiles)]),
            ylim = c(0, 1.2), type = "b", pch = 1,
            ylab = " Theta 1,2,3", main = main, ...)
        points(quantiles, theta2, pch = 2, col = 3)
        points(quantiles, theta3, pch = 4, col = 4) }     
    
    # Return Value:
    data.frame(thresholds=z0, theta1=theta1, theta2=theta2, theta3=theta3)
}


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


exindexPlot = 
function(x, block = "month", start = 5, end = NA, 
plottype = c("thresh", "K"), labels = TRUE, autoscale = TRUE, ...)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #	Calculates Theta(1,2) and plots Theta(1)
 
    # Notes:
    #   Wraps "exindex" from Alexander McNeil's evir package

    # FUNCTION:
    
	# Wrapper:
	plottype = plottype[1]
	reverse = FALSE
	if (plottype == "K") reverse = TRUE
	ans = exindex(data = x, block = block , start = start, end = end, 
		reverse = reverse, auto.scale = autoscale, labels = labels, ...) 

    # Return Value:
    ans
}


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

# 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:
#  xmpfSeries            Popups the example menu
################################################################################

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

	# Description:
	#	Popups the example menu
	
	# FUNCTION:
	
	# Popup:	
	path = paste(.Library,"/fExtremes", 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: ", "fExtremes", "\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()
}
	
	
# ******************************************************************************


#*******************************************************************************
# fExtremes - A SOFTWARE COLLECTION FOR FINANCIAL ENGINEERS
# PART III: Beyond the Sample: Dealing with Extreme Values
#
# collected by Diethelm Wuertz
# Version 0.9
#*******************************************************************************


# 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
#   evir: original S functions (EVIS) by Alexander McNeil <mcneil@math.ethz.ch>
#     R port by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   ismev: Original S functions by Stuart Coles <Stuart.Coles@bristol.ac.uk>
#     R port/documentation by Alec Stephenson <a.stephenson@lancaster.ac.uk>
#   evd: Alec Stephenson <alec_stephenson@hotmail.com>


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


# Default Settings:
	
	xmpExtremes = function(prompt = "") {invisible(prompt)}


.First.lib = 
function(lib, pkg)
{ 	# A function implemted by D. Wuertz

	# Package:
	cat("\nfExtremes:  Beyond the Sample: Dealing with Extreme Values")
	
	# Requires:
	# minor = as.numeric(version$minor)
	# DEBUG = TRUE
	# sink("@sink@")           					# available from:
    # library(evd, warn.conflicts = DEBUG)		# CRANbin
	# library(ismev, warn.conflicts = DEBUG)	# CRANbin
    # library(evir, warn.conflicts = DEBUG)		# evir | Now on Cran
	# sink()
	# unlink("@sink@")
	

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


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

