#
# Program to illustrate likelihood calculations for Jenkinson's Hartford data
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]
# first fit: GEV with no covariates
# neg log likelihood function
lh1=function(b){
if(abs(b[3])>1){return(1.0e10)}
if(abs(b[2])>10){return(1.0e10)}
mu=b[1]
sig=exp(b[2])
xi=b[3]
n=length(y)
lh1=0
for(i in 1:n){
aa=1+xi*(y[i]-mu)/sig
if(aa<0){return(1.0e10)}
aa=log(aa)
lh1=lh1+log(sig)+(1+1/xi)*aa
if(abs(aa/xi)>20){return(1.0e10)}
lh1=lh1+exp(-aa/xi)
}
return(lh1)
}
# use trial starting values (could vary these to confirm result)
par=c(mean(y),log(sqrt(var(y))),0.0011)
#
# initial log likelihood (before optimization)
#
print(lh1(par))
#[1] 265.085
#
# first optimization, note control parameters ndeps, maxit
#
m1=optim(par,lh1,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(m1$val)
# [1] 245.9712
V=solve(m1$hessian)
res=matrix(nrow=3,ncol=4)
for(i in 1:3){
res[i,1]=m1$par[i]
res[i,2]=sqrt(V[i,i])
res[i,3]=res[i,1]/res[i,2]
res[i,4]=2*pnorm(-abs(res[i,3]))
}
print(round(res,4))
#         [,1]   [,2]    [,3] [,4]
# [1,] 19.6809 0.3894 50.5438    0
# [2,]  1.2467 0.0786 15.8616    0
# [3,] -0.2575 0.0596 -4.3229    0
#
# second optimization using nlm procedure
#
m2=nlm(lh1,par,hessian=T)
print(m2$min)
# [1] 245.9712
V=solve(m2$hessian)
res=matrix(nrow=3,ncol=4)
for(i in 1:3){
res[i,1]=m2$est[i]
res[i,2]=sqrt(V[i,i])
res[i,3]=res[i,1]/res[i,2]
res[i,4]=2*pnorm(-abs(res[i,3]))
}
print(round(res,4))
#         [,1]   [,2]    [,3] [,4]
# [1,] 19.6809 0.3967 49.6083    0
# [2,]  1.2467 0.0786 15.8514    0
# [3,] -0.2575 0.0598 -4.3033    0
#
# Estimate endpoint of distribution and its SE
#
theta=res[1,1]-exp(res[2,1])/res[3,1]
grad=c(1,-exp(res[2,1])/res[3,1],exp(res[2,1])/res[3,1]^2)
se.theta=sqrt(grad %*% V %*% grad )
print(c(theta,se.theta))
# [1] 33.192655  2.597155
#
# estimate N-yr RV and its SE
#
N=100
theta=res[1,1]+exp(res[2,1])*((-log(1-1/N))^(-res[3,1])-1)/res[3,1]
grad=c(1,exp(res[2,1])*((-log(1-1/N))^(-res[3,1])-1)/res[3,1],
exp(res[2,1])*(-(-log(1-1/N))^(-res[3,1])*log(-log(1-1/N))/res[3,1]
-((-log(1-1/N))^(-res[3,1])-1)/res[3,1]^2))
se.theta=sqrt(grad %*% V %*% grad )
print(c(theta,se.theta))
# [1] 29.0589556  0.8315969
N=1000
theta=res[1,1]+exp(res[2,1])*((-log(1-1/N))^(-res[3,1])-1)/res[3,1]
grad=c(1,exp(res[2,1])*((-log(1-1/N))^(-res[3,1])-1)/res[3,1],
exp(res[2,1])*(-(-log(1-1/N))^(-res[3,1])*log(-log(1-1/N))/res[3,1]
-((-log(1-1/N))^(-res[3,1])-1)/res[3,1]^2))
se.theta=sqrt(grad %*% V %*% grad )
print(c(theta,se.theta))
# [1] 30.910389  1.304159
#
# Bayesian calculation: first load "AdaptMH" function
#
par=m1$par
npar=length(par)
C0=solve(m1$hessian)
nsim=200000
nn1=nsim/10
nn2=50
nn3=10
eps=0.001
scal=2.4
set.seed(123)
T1=Sys.time()
bay1=AdaptMH(lh1,par,npar,C0,nsim,nn1,nn2,nn3,eps,scal)
T2=Sys.time()
print(T2-T1)
print(bay1$acc0/nsim)
# discard first half of chain as warmup


bay1$parout=bay1$parout[10001:20000,]


# return value calculation
N=100
theta.out=rep(NA,nrow(bay1$parout))
for(i in 1:nrow(bay1$parout)){
theta.out[i]=bay1$parout[i,1]+exp(bay1$parout[i,2])*((-log(1-1/N))^(-bay1$parout[i,3])-1)/bay1$parout[i,3]
}
print(quantile(theta.out,c(0.025,0.5,0.975)))
# 28.09363 29.32322 32.02013 
N=1000
theta1.out=rep(NA,nrow(bay1$parout))
for(i in 1:nrow(bay1$parout)){
theta1.out[i]=bay1$parout[i,1]+exp(bay1$parout[i,2])*((-log(1-1/N))^(-bay1$parout[i,3])-1)/bay1$parout[i,3]
}
print(quantile(theta1.out,c(0.025,0.5,0.975)))
# 29.55644 31.27507 35.97726 - this is 95% credible interval

# endpoint calculation
theta2.out=rep(NA,nrow(bay1$parout))
for(i in 1:nrow(bay1$parout)){
theta2.out[i]=bay1$parout[i,1]-exp(bay1$parout[i,2])/bay1$parout[i,3]
}
print(quantile(theta2.out,c(0.025,0.5,0.975)))
# 30.75727 33.83849 48.64663 
#
# posterior density plots - we need to throw out extreme values of theta2.out to stabilize the density calculation
theta2.out=theta2.out[theta2.out>quantile(theta2.out,0.005)&theta2.out<quantile(theta2.out,0.995)]
par(cex=1.3)
plot(density(theta.out),xlim=c(28,50),xlab='Level',main='Posterior densities')
lines(density(theta1.out),lty=2)
lines(density(theta2.out),lty=3)


# alternative plots create smoother picture (double the bandwidth)

plot(density(theta.out,bw=2*density(theta.out)$bw),xlim=c(28,50),xlab='Level',main='Posterior densities')
lines(density(theta1.out,bw=2*density(theta1.out)$bw),lty=2)
lines(density(theta2.out,bw=2*density(theta2.out)$bw),lty=3)



### analysis using ismev

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)

gev.diag(mod1)