library(ismev)

#point process simulation for extremes
rpareto <- function(n = 1, alpha = 1)
{
	#simulates a draw from a pareto distribution
	#uses the inverse law
	#pareto:  F(x) = 1 - x^(-alpha)
	
	u <- runif(n)
	rp <- (1-u)^(-1/alpha)
	
	return(rp)
}

alpha <- 3

n <- 10000

x <- rpareto(n, alpha)
bn <- n^(1/alpha)
an <- 1/alpha*n^(1/alpha)
pts <- cbind(seq(1:n)/n, (x - bn)/an)
plot(pts, ylim = c(min(pts[,2]), max(max(pts[,2]), 2)))
abline(h = -1.5, lty = 2)



#Fort Collins Precip
fc <- data.matrix(read.table("Research/RetiredProjects/FrontRange/Data/OrigUliData/fortCollins.dat"))
fc[c(1:10, 10870:10876),]
dtot <- fc[,8]/100
plot(dtot, axes = F, type = 'h', ylab = "daily precip (in)", xlab = "date")
box()
axis(2)
axis(1, at = c(306, 2434, 4525, 6649, 8541, 10478), labels = c("1950", "1960", "1970", "1980", "1990", "2000"))
title("Ft. Collins Summer Precip")
abline(h = 0.5, lty = 3)

#fit a gpd
ftcGPDFit <- gpd.fit(dtot, threshold = 0.5)
gpd.diag(ftcGPDFit)

#fit a PP representation
ftcPPFit <- pp.fit(dtot, threshold = 0.5)


###################
#threshold selection Fort Collins Data
ftcGPDFit2 <- gpd.fit(dtot, threshold = 0.05)
gpd.diag(ftcGPDFit2)

ftcGPDFit3 <- gpd.fit(dtot, threshold = 2)
gpd.diag(ftcGPDFit3)


###################
# a better threshold selection example.
set.seed(39)
randomT <- rt(10000, = 4)
plot(randomT)

#an okay choice
gpd.fit(randomT, 3)

#????
gpd.fit(randomT, 2)

#too low
gpd.fit(randomT, 0.5)

#even worse
gpd.fit(randomT, 0)

#too high
gpd.fit(randomT, 5)


#################
#mrl.plot
mrl.plot(dtot)
mrl.plot(dtot, umin = 0, umax = 3)

mrl.plot(randomT)
mrl.plot(randomT, umin = 0, umax = 7)


###################
gpd.fitrange(dtot, umin = 0, umax = 4)
gpd.fitrange(dtot, umin = 0, umax = 2)

gpd.fitrange(randomT, umin = 0, umax = 7)