#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")

#S&P 500 log returns
sp500All <- data.matrix(read.csv("Research/PredExtremes/Data/Financial/sp500.csv"))
sp500 <- rev(sp500All[,7])
l <- length(sp500)
logReturns <- log(sp500[2:l]) - log(sp500[1:(l-1)])
plot(logReturns, type = 'h', axes = F)
box()
axis(2)
axis(1, at = c(1, l - 15), labels = c("1980", "2009")) 
abline(v = l - 15, lty = 3)


#wave-surge bivariate data
library(ismev)
data(wavesurge)
plot(wavesurge)


#analysis of the Fort Collins data 1948-1990
dtotShort <- dtot[1:8541]
fcNonZero <- dtotShort[dtotShort > 0]
rate <- length(fcNonZero)/length(dtotShort)
rate
hist(fcNonZero, breaks = seq(0, 4.2, .1), freq = F)
#fit a gamma to non-zero observations
#first, an example
x <- seq(0, 4.2, .01)
y2 <- dgamma(x, shape = 1, rate = 2)
lines(x,y2, col = 2)
y3 <- dgamma(x, shape = 1, rate = 3)
lines(x, y3, col = 3)
#ml estimation
lhoodGamma <- function(par)
{
	#returns negative log-lhood:  lower is better
	out <- -sum(dgamma(fcNonZero, shape = par[1], 
				rate = par[2], log = T))
	return(out)
}
lhoodGamma(c(1,2))
lhoodGamma(c(1,3))

hist(fcNonZero, breaks = seq(0, 4.2, .05), freq = F)
opt <- optim(c(1,3), lhoodGamma)
opt
yOpt <- dgamma(x, shape = opt$par[1], rate = opt$par[2])
lines(x, yOpt, col = 4)

#what is probability an observation exceeds 6.18?
a <- 1 - pgamma(6.18, shape = opt$par[1], rate = opt$par[2]) 
a
b <- rate*a
b
pExc <- 1 - (1 - b)^214

#extract annual maxima
annMaxSeq <- numeric(42)
yearBase <- 1948
for(i in seq(1, 42))
{
	curYear <- yearBase + i
	dtotYr <- dtot[fc[, "year"] == yearBase + i]
	annMaxSeq[i] <- max(dtotYr)
}
gev.fit(annMaxSeq)
library(evd)
pExc2 <- 1 - pgev(6.18, 1.11, 0.46, 0.31)
pExc2
1/pExc2


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

#exploration of sample maxima


#m:  number of samples
#n:  sample size

#normals
n <- 100
m <- 1000

r100 <- matrix(rnorm(m*n), nrow = m, ncol = n)
sampleMax100 <- apply(r100, 1, max)
hist(sampleMax100, breaks = 
  seq(floor(min(sampleMax100)*4)/4, 
  	ceiling(max(sampleMax100)*4)/4, .25))

n <- 1000
m <- 1000

r1000 <- matrix(rnorm(m*n), nrow = m, ncol = n)
sampleMax1000 <- apply(r1000, 1, max)
op <- par(mfrow = c(1,2))
hist(sampleMax100, breaks = 
  seq(floor(min(sampleMax100)*4)/4, 
  	ceiling(max(sampleMax100)*4)/4, .25))
hist(sampleMax1000, breaks = 
  seq(floor(min(sampleMax1000)*4)/4, 
  	ceiling(max(sampleMax1000)*4)/4, .25))
par(op)


#exponentials
n <- 100
m <- 1000

r100 <- matrix(rexp(m*n), nrow = m, ncol = n)
sampleMax100 <- apply(r100, 1, max)
hist(sampleMax100, breaks = 
  seq(floor(min(sampleMax100)*4)/4, 
  	ceiling(max(sampleMax100)*4)/4, .25))

n <- 1000
m <- 1000

r1000 <- matrix(rexp(m*n), nrow = m, ncol = n)
sampleMax1000 <- apply(r1000, 1, max)
op <- par(mfrow = c(1,2))
hist(sampleMax100, breaks = 
  seq(floor(min(sampleMax100)*4)/4, 
  	ceiling(max(sampleMax100)*4)/4, .25))
hist(sampleMax1000, breaks = 
  seq(floor(min(sampleMax1000)*4)/4, 
  	ceiling(max(sampleMax1000)*4)/4, .25))
par(op)


#betas
n <- 100
m <- 1000

r100 <- matrix(rbeta(m*n, 3, 3), nrow = m, ncol = n)
sampleMax100 <- apply(r100, 1, max)
hist(sampleMax100, breaks = 10)

n <- 1000
m <- 1000

r1000 <- matrix(rbeta(m*n, 3, 3), nrow = m, ncol = n)
sampleMax1000 <- apply(r1000, 1, max)
op <- par(mfrow = c(1,2))
hist(sampleMax100, breaks = 10)
hist(sampleMax1000, breaks = 10)
par(op)


#t with 4 df
n <- 100
m <- 1000

r100 <- matrix(rt(m*n, 4), nrow = m, ncol = n)
sampleMax100 <- apply(r100, 1, max)
hist(sampleMax100, breaks = 10)

n <- 1000
m <- 1000

r1000 <- matrix(rt(m*n, 4), nrow = m, ncol = n)
sampleMax1000 <- apply(r1000, 1, max)
op <- par(mfrow = c(1,2))
hist(sampleMax100, breaks = 20)
hist(sampleMax1000, breaks = 20)
par(op)


##################################################################
#review of FTC rainfall
plot(annMaxSeq, type = 'h')
ftcGevFit <- gev.fit(annMaxSeq)
gev.diag(ftcGevFit)

#what is probability of ann max exceeding 6.18?
library(evd)
pExc <- 1 - pgev(6.18, 1.114, 0.4611, 0.31)
1/pExc


#hartford wind
data(wind)
wind[1:10,]
plot(wind[,c(1,2)], type = 'h')
hartfordWindFit <- gev.fit(wind[,2])
gev.diag(hartfordWindFit)
pExc <- 1 - pgev(70, 49.93, 5.02, 0.004)
pExc
1/pExc
#100-year RL?
qgev(.99, 49.93, 5.02, 0.004)

#method-of-moments?
library(lmom)
lmom <- samlmu(wind[,2])
lmom
out <- pelgev(lmom)

#profile likelihood
#compute CI's based on asymptotic normality
hartfordWindFit$cov
sqrt(diag(hartfordWindFit$cov))
ciNorm <- matrix( c(hartfordWindFit$mle + qnorm(.025)*hartfordWindFit$se, hartfordWindFit$mle + qnorm(.975)*hartfordWindFit$se), ncol = 2)
ciNorm

#profile likelihood
#xi
gev.profxi(hartfordWindFit, xlow = -.3, xup = .3)
abline(v = -.194, lty = 3)
abline(v = .202, lty = 3)

#RL
gev.prof(hartfordWindFit, m = 100, xlow = 60, xup = 110)