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


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