### Suppose vector x contains the loss data
mu1 = 5        # the threshhold of 7B
 gpdmod = gpdFit(x,u=mu1)        #Fit GPD
 gpdpar = slot(gpdmod,"fit")    #get the parameters


 x1 = x[x >= mu1]    #get data exceeding the threshold


xeval = seq(mu1, (max(x1)+sd(x1)), length = 100)    #create a grid

#gpdeval = dgpd(xeval,xi=gpdpar$parfit$par[1],beta=gpdpar$fit$par[2],mu=mu1)
gpdeval = dgpd(xeval,xi=gpdpar$par.ests[1],beta=gpdpar$par.ests[2],mu=mu1)

hist(x1,probability=T, xlab="Loss data (B$)",ylab="PDF", ylim=range(gpdeval))
lines(xeval,gpdeval,col="red")

## Q-Q Plot

#empquant = c(1:length(x1))/(length(x1) + 1)
#gpdquant = qgpd(empquant,xi=gpdpar$fit$par[1],beta=gpdpar$fit$par[2],mu=mu1)
#gpdquant = qgpd(empquant,xi=gpdpar$par.ests[1],beta=gpdpar$par.ests[2],mu=mu1)
#x2=sort(x1)
#plot(x2,gpdquant,xlab="Observed Quantile", ylab="GPD Quantile",log="xy")	#log-log scale
#lines(gpdquant,gpdquant)
#grid()
### you can use qqplot to obtain the Q-Q plot of the GPD
qqplot(qgpd(ppoints(length(x1),a=0),xi=gpdpar$fit$par[1],beta=gpdpar$fit$par[2],mu=mu1), x1, xlab="Theoretical Quantile", ylab="Sample Quantile")
lines(sort(x1),sort(x1),col="red")
title(main="Q-Q plot from GPD")