#Threshold exceedance for Normal variables
qu <- seq(1,5,by=0.02)
penult <- smith.penult(family = "norm", ddensF=function(x){-x*dnorm(x)},
method = 'pot', u = qu)
plot(qu, penult$shape, type='l', xlab='Quantile',
ylab='Penultimate shape', ylim=c(-0.5,0))
#Block maxima for Gamma variables -
#User must provide arguments for shape (or rate)
m <- seq(30, 3650, by=30)
penult <- smith.penult(family = 'gamma', method = 'bm', m=m, shape=0.1)
plot(m, penult$shape, type='l', xlab='Quantile', ylab='Penultimate shape')
#Comparing density of GEV approximation with true density of maxima
m <- 100 #block of size 100
p <- smith.penult(family='norm',
ddensF=function(x){-x*dnorm(x)}, method='bm', m=m, returnList=FALSE)
x <- seq(1, 5, by = 0.01)
plot(x, m*dnorm(x)*exp((m-1)*pnorm(x,log.p=TRUE)),type='l', ylab='Density',
main='Distribution of the maxima of\n 100 standard normal variates')
lines(x, mev::dgev(x,loc=p[1], scale=p[2], shape=0),col=2)
lines(x, mev::dgev(x,loc=p[1], scale=p[2], shape=p[3]),col=3)
legend(x = 'topright',lty = c(1,1,1,1), col = c(1,2,3,4),
legend = c('exact', 'ultimate', 'penultimate'), bty = 'n')
Run the code above in your browser using DataLab