################################################################
#### Section 8.2.2: The Peaks over Threshold Method
################################################################

# Parameter estimation of the GPD, Figure 8.13

# Functions

gpdfitls<-function(sorteddatavec,initguess)
{
n<-length(sorteddatavec)
optfun<-function(params){sum((params[2]/params[1]*(((1:n)/(n+1))^(-params[1])-1)-sorteddatavec)^2)}
optim(initguess,optfun)
}

gpdfitml<-function(datavec,initguess)
{
n<-length(datavec)
optfun<-function(params){n*log(params[2])+(1/params[1]+1)*sum(log(1+params[1]/params[2]*datavec))}
optim(initguess,optfun)
}

gpdparamspa3<-function(runs)
{
resmat<-matrix(0,runs,4)
tmpvals<-(1:1000)
for(i in (1:runs))
{
tmpvals<-sort(runif(1000)^(-1/3),dec=T)
resmat[i,1:2]<-gpdfitls(tmpvals[1:100]-tmpvals[101],c(1/3,1))$par
resmat[i,3:4]<-gpdfitml(tmpvals[1:100]-tmpvals[101],c(1/3,1))$par
}
resmat
}

# Program

gpdparmat<-gpdparamspa3(1000)
xmin<-min(gpdparmat[,1],gpdparmat[,3])
xmax<-max(gpdparmat[,1],gpdparmat[,3])
ymin<-min(gpdparmat[,2],gpdparmat[,4])
ymax<-max(gpdparmat[,2],gpdparmat[,4])

# Plot LS parameter estimates, Figure 8.13
plot(gpdparmat[,1],gpdparmat[,2],xlab="",ylab="",main="",xlim=c(xmin,xmax),ylim=c(ymin,ymax))
abline(h=(1-0.9)^(-1/3)/3)
abline(v=1/3)

# Plot ML parameter estimates, Figure 8.13
plot(gpdparmat[,3],gpdparmat[,4],xlab="",ylab="",main="",xlim=c(xmin,xmax),ylim=c(ymin,ymax))
abline(h=(1-0.9)^(-1/3)/3)
abline(v=1/3)