library(ismev) data(rain) ##install.packages("evir") ##library(evir) data(danish,package="evir") library(evd) data(portpirie) data(sealevel) ################################################################################## ### Functions to be given to the students. ################# Diagnostiques graphiques: Gumbel plot et Pareto Quantile plot. ### Gumbel plot : (-log(-log(i/n)), x_(i) ) gumPlot <- function(z, imin=1){ par(mfrow=c(1,1)) nz<-length(z) zord<-sort(z)[imin:nz] auxvec=imin:nz/(nz+1) gbquant<- -log(-log(auxvec)) plot(gbquant,zord, xlab="theoretical gumbel quantiles", ylab="order statistics", main="Adequacy of Gumbel model") reg.lin<-lm(zord~gbquant) coeffgb<- reg.lin$coefficients abline(coeffgb[1],coeffgb[2], col="red", lwd=2) print(summary( reg.lin)) invisible(reg.lin) } ## Pareto Quantile plot : (-log(1-i/(n+1));log(x_(i))) paretoPlot <- function(z,threshold=quantile(z,prob=0.1)) { y <- z[z>=threshold] n<-length(y) inds<-1:n yord<-sort(y) expquants <- -log(1-inds/(n+1)) plot(expquants,log(yord), main="Adequacy of heavy-tailed GPD model" , xlab="theoretical exponential quantiles", ylab="log(order statistics)") reg.lin<-lm( log(yord)~ expquants ) coeff<- reg.lin$coefficients abline(coeff[1],coeff[2], col="red", lwd=2) print(summary( reg.lin)) invisible(reg.lin) } ### extraction de maxima extractMax <- function(data, npp=365) { if(npp!=floor(npp)) stop("npp must be an integer valued") n <- length(data) nperiods <- n %/% npp nkept <- nperiods * npp Data <- matrix(data[1:nkept],ncol=nperiods) mdata <- apply(Data,2,max) return( mdata) } hill<- function(data, plot=FALSE, add=FALSE, starica=FALSE,...) { ## Computes the Hill estimates of gamma (Section 4.2) ## for a numeric vector of observations (data) and as a ## function of k ## ## If plot=TRUE then the estimates are plotted as a ## function of k ## ## If add=TRUE then the estimates are added to an existing plot. X <- sort(data) ##increasing order n <- length(X) Hill <- numeric(n) K <- 1:(n-1) ### Hill estimates for (k in 1:(n-1)) { ##k: number of observation kept to estimate the mean excess plot. Hill[k] <- (1/k)*sum(log(X[n:(n-k+1)])) - log(X[n-k]) ## print(k) } ### plots if TRUE if (plot || add){ if ( !add ) { ### plot estimates if(!starica) plot(K, Hill[K], type="l", ylab="gamma", xlab="k", ,main="Estimates of extreme value index", ...) else plot(log(K), Hill[K], type="l", ylab="gamma", xlab="log(k)", ,main="Estimates of extreme value index", ...) } else { ### adds estimates to existing plot if(!starica) lines(K, Hill[K], ...) else lines(log(K), Hill[K], ...) } } ### output list with values of k and ### corresponding Hill estimates list(k=K, gamma=Hill[K]) }