# revision of HartfordGEV.txt file

# incorporate ismev and use RLS gev.text function

# install this first
gev.test=function(xdat, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, 
    mulink = identity, siglink = identity, shlink = identity, 
    muinit = NULL, siginit = NULL, shinit = NULL,
    show = FALSE, method = "Nelder-Mead", maxit = 10000,nsim){
fit1=gev.fit(xdat, ydat = ydat, mul = mul, sigl = sigl, shl = shl, 
    mulink = mulink, siglink = siglink, shlink = shlink, 
    muinit = muinit, siginit = siginit, shinit = shinit,
    show = show, method = method, maxit = maxit)
u1=sort(exp(-(1+fit1$vals[,3]*(xdat-fit1$vals[,1])/fit1$vals[,2])^(-1/fit1$vals[,3])))
n=length(u1)
plotpos=qnorm((1:n-0.5)/n)
lg=cor(qnorm(u1),plotpos)
d1<-max(c((1:n)/n-u1,u1-(0:(n-1))/n))
w1<-sum((u1-((1:n)-0.5)/n)^2)+1/(12*n)
a1<--sum((2*(1:n)-1)*log(u1)+(2*n+1-2*(1:n))*log(1-u1))/n-n
count=rep(0,4)
for(isim in 1:nsim){
ysim=fit1$vals[,1]+fit1$vals[,2]*((-log(runif(n)))^(-fit1$vals[,3])-1)/fit1$vals[,3]
fit2=gev.fit(ysim, ydat = ydat, mul = mul, sigl = sigl, shl = shl, 
    mulink = mulink, siglink = siglink, shlink = shlink, 
    muinit = muinit, siginit = siginit, shinit = shinit,
    show = show, method = method, maxit = maxit)
u2=sort(exp(-(1+fit2$vals[,3]*(ysim-fit2$vals[,1])/fit2$vals[,2])^(-1/fit2$vals[,3])))
lg2=cor(qnorm(u2),plotpos)
d2<-max(c((1:n)/n-u2,u2-(0:(n-1))/n))
w2<-sum((u2-((1:n)-0.5)/n)^2)+1/(12*n)
a2<--sum((2*(1:n)-1)*log(u2)+(2*n+1-2*(1:n))*log(1-u2))/n-n
if(lg2<lg)count[1]<-count[1]+1
if(d2>d1)count[2]<-count[2]+1
if(w2>w1)count[3]<-count[3]+1
if(a2>a1)count[4]<-count[4]+1
}
ddf=matrix(c(lg,d1,w1,a1,count/nsim),byrow=T,ncol=4)
ddf=round(ddf,5)
ddf=cbind(c('Test','p-value'),ddf)
ddf=rbind(c('','LG','KS','CVM','AD'),ddf)
return(ddf)
}
# end of installation

# recreate Jenkinson dataset and run gev.fit
D=matrix(c(12,1,14,2,15,4,16,4,17,3,18,4,19,11,20,9,21,21,22,6,23,8,24,3,25,5,26,5,27,3,28,1,29,1,30,1),byrow=T,ncol=2)
y=0
for(i in 1:18){y=c(y,rep(D[i,1],D[i,2]))}
y=y[2:93]

library(ismev)
mod1=gev.fit(y, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL,
mulink = identity, siglink = identity, shlink = identity,
muinit = mean(y), siginit = sd(y), shinit = 0.001,
show = TRUE, method = "Nelder-Mead", maxit = 10000)

# output
# $mle
# [1] 19.6805415  3.4786618 -0.2574008
# 
# $se
# [1] 0.39674953 0.27357205 0.05978129
#
# Jenkinson p. 204 quoted 19.68, 3.48, -0.26


# ismev diagnostics
gev.diag(mod1)

# return value calculations by MLE

# endpoint
theta=mod1$mle[1]-mod1$mle[2]/mod1$mle[3]
g=c(1,-1/mod1$mle[3],mod1$mle[2]/mod1$mle[3]^2)
sqrt(g %*% mod1$cov %*% g)
# theta=33.19511 with SE=2.596093 - Jenkinson (p 212) gave 33.2, 2.39

# N-year return value
N=1000
theta=mod1$mle[1]+mod1$mle[2]*((-log(1-1/N))^(-mod1$mle[3])-1)/mod1$mle[3]
q=-log(1-1/N)
g=c(1,(q^(-mod1$mle[3])-1)/mod1$mle[3],-mod1$mle[2]*(q^(-mod1$mle[3])*log(q)/mod1$mle[3]+(q^(-mod1$mle[3])-1)/mod1$mle[3]^2))
sqrt(g %*% mod1$cov %*% g)
# 30.91132 with SE 1.303303 - Jenkinson gave 30.9, 1.21

# gof test (new calculation)

T1=Sys.time()
gev.test(y, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL,
mulink = identity, siglink = identity, shlink = identity,
muinit = mean(y), siginit = sd(y), shinit = 0.001,
show = FALSE, method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)

# [1,] ""        "LG"      "KS"      "CVM"     "AD"     
# [2,] "Test"    "0.98981" "0.13006" "0.19857" "0.95802"
# [3,] "p-value" "0.043"   "0"       "0.001"   "0.003"  

# not great result - try adding random noise to the data
set.seed(123)
y2=y-0.5+runif(92)

T1=Sys.time()
gev.test(y2, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL,
mulink = identity, siglink = identity, shlink = identity,
muinit = mean(y), siginit = sd(y), shinit = 0.001,
show = F, method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)

# [1,] ""        "LG"      "KS"      "CVM"    "AD"     
# [2,] "Test"    "0.99246" "0.08664" "0.1459" "0.69558"
# [3,] "p-value" "0.134"   "0.045"   "0.009"  "0.03"   

# better but still not great


data(venice)
y=venice[,2]

mod2=gev.fit(y, ydat = matrix((1:51),ncol=1), mul = 1, sigl = NULL,
shl = NULL, mulink = identity, siglink = identity, shlink = identity,
muinit = c(mean(y),0), siginit = sqrt(var(y)), shinit = -0.001, show = T,
method = "Nelder-Mead", maxit = 10000)

gev.diag(mod2)

T1=Sys.time()
gev.test(y, ydat = matrix((1:51),ncol=1), mul = 1, sigl = NULL,
shl = NULL, mulink = identity, siglink = identity, shlink = identity,
muinit = c(mean(y),0), siginit = sqrt(var(y)), shinit = -0.001, show = F,
method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)