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)