V=scan('C:/Users/rls/aug20/UNC/STOR834/Book/Chapter-1/venice.txt')
length(V)


V=matrix(V,ncol=10,byrow=F)
V=V[1:51,]
V=cbind(c(rep(10,4),6,rep(10,46)),V)

V=cbind(1931:1981,V)

V=data.frame(V)
names(V)=c('Year','Num','Y1','Y2','Y3','Y4','Y5','Y6','Y7','Y8','Y9','Y10')


write.csv(V,'C:/Users/rls/aug20/UNC/STOR834/Book/Chapter-1/Venice.csv',row.names=F)


# restart data analysis
V=read.csv('C:/Users/rls/aug20/UNC/STOR834/Book/Chapter-1/Venice.csv')

V=read.csv('https://rls.sites.oasis.unc.edu/s834-2023/Data/Venice.csv')

# likelihood function (compute NLLH - defaults to 10^10 if parameter values infeasible)
#
# data frame here is V
# par vector is (beta0, log psi, xi, beta1, beta2)
#
# r largest obs per year
#
# model=0 is no trend
# model=1 is linear trend
# model=2 is quadratic trend
# model=3 is linear trend + periodic of period P (years)
lh=function(par){
if(abs(par[2])>20){return(10^10)}
if(abs(par[3])>1){return(10^10)}
mu=par[1]
psi=exp(par[2])
xi=par[3]
f=0
for(i in 1:51){
ri=min(r,V$Num[i])
if(model==1){mu=par[1]+(i/10)*par[4]}
if(model==2){mu=par[1]+(i/10)*par[4]+(i/10)*(i/10)*par[5]}
if(model==3){mu=par[1]+(i/10)*par[4]+cos(2*pi*i/P)*par[5]+sin(2*pi*i/P)*par[6]}
f=f+ri*par[2]
s1=1+xi*(V[i,2+ri]-mu)/psi
if(s1<=0){return(10^10)}
s1=-log(s1)/xi
if(abs(s1)>50){return(10^10)}
f=f+exp(s1)
for(j in 1:ri){
s1=1+xi*(V[i,2+j]-mu)/psi
if(s1<=0){return(10^10)}
f=f+(1+1/xi)*log(s1)
}}
return(f)
}

model=0
r=10
par=c(mean(V$Y1),log(sd(V$Y1)),0.01)
lh(par)
opt1=optim(par,lh,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(opt1$val)
tab1=matrix(nrow=length(par),ncol=4)
tab1[,1]=opt1$par
tab1[,2]=sqrt(diag(solve(opt1$hessian)))
tab1[,3]=tab1[,1]/tab1[,2]
tab1[,4]=2*pnorm(-abs(tab1[,3]))
print(round(tab1,4))


model=1
r=10
par=c(mean(V$Y1),log(2*sd(V$Y1)),-0.001,0)
lh(par)
opt2=optim(par,lh,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(opt2$val)
tab1=matrix(nrow=length(par),ncol=4)
tab1[,1]=opt2$par
tab1[,2]=sqrt(diag(solve(opt2$hessian)))
tab1[,3]=tab1[,1]/tab1[,2]
tab1[,4]=2*pnorm(-abs(tab1[,3]))
print(round(tab1,4))

model=2
r=10
par=c(opt2$par,0)
lh(par)
opt3=optim(par,lh,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(opt3$val)
tab1=matrix(nrow=length(par),ncol=4)
tab1[,1]=opt3$par
tab1[,2]=sqrt(diag(solve(opt3$hessian)))
tab1[,3]=tab1[,1]/tab1[,2]
tab1[,4]=2*pnorm(-abs(tab1[,3]))
print(round(tab1,4))

model=3
r=10
P=18.62
par=c(opt2$par,0,0)
lh(par)
opt4=optim(par,lh,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(opt4$val)
tab1=matrix(nrow=length(par),ncol=4)
tab1[,1]=opt4$par
tab1[,2]=sqrt(diag(solve(opt4$hessian)))
tab1[,3]=tab1[,1]/tab1[,2]
tab1[,4]=2*pnorm(-abs(tab1[,3]))
print(round(tab1,4))

model=3
r=10
P=11
par=c(opt2$par,0,0)
lh(par)
opt5=optim(par,lh,method="BFGS",control=list(ndeps=rep(1e-6,length(par)),maxit=10000),hessian=T)
print(opt5$val)
tab1=matrix(nrow=length(par),ncol=4)
tab1[,1]=opt5$par
tab1[,2]=sqrt(diag(solve(opt5$hessian)))
tab1[,3]=tab1[,1]/tab1[,2]
tab1[,4]=2*pnorm(-abs(tab1[,3]))
print(round(tab1,4))

par(cex=1.3)
plot(V$Year,V$Y1,pch=20,xlab='Year',ylab='Largest Sea Level')

lines(V$Year,opt2$par[1]+(1:51)*0.1*opt2$par[4]+exp(opt2$par[2])*((log(2))^(-opt2$par[3])-1)/opt2$par[3])
lines(V$Year,opt3$par[1]+(1:51)*0.1*opt3$par[4]+(1:51)^2*0.01*opt3$par[5]+exp(opt3$par[2])*((log(2))^(-opt3$par[3])-1)/opt3$par[3],lty=2)

lines(V$Year,opt4$par[1]+(1:51)*0.1*opt4$par[4]+cos(2*pi*(1:51)/18.62)*opt4$par[5]+sin(2*pi*(1:51)/18.62)*opt4$par[6]
+exp(opt4$par[2])*((log(2))^(-opt4$par[3])-1)/opt4$par[3],lty=3)

lines(V$Year,opt5$par[1]+(1:51)*0.1*opt5$par[4]+cos(2*pi*(1:51)/11)*opt5$par[5]+sin(2*pi*(1:51)/11)*opt5$par[6]
+exp(opt5$par[2])*((log(2))^(-opt5$par[3])-1)/opt5$par[3],lty=4)

# analysis using ismev

V=read.csv('https://rls.sites.oasis.unc.edu/s834-2023/Data/Venice.csv')
V2=V[,3:12]
V2[5,7:10]=NA

library(ismev)
mod2=rlarg.fit(V2, r = 10, ydat = matrix((1:51),ncol=1), mul = 1, sigl = NULL,
shl = NULL, mulink = identity, siglink = identity, shlink = identity,
muinit = c(mean(V2[,1]),0), siginit = sqrt(var(V2[,1])), shinit = -0.001, show = F,
method = "Nelder-Mead", maxit = 10000)

rlarg.diag(mod2)