R/uvgev-gen-retlev.R

In POT: Generalized Pareto Distribution and Peaks Over Threshold

#############################################################################
#   Copyright (c) 2014 Mathieu Ribatet              
#   Copyright (c) 2022 Christophe Dutang => replace fitted to object 
#                                                                                                                                                                        
#   This program is free software; you can redistribute it and/or modify                                               
#   it under the terms of the GNU General Public License as published by                                         
#   the Free Software Foundation; either version 2 of the License, or                                                   
#   (at your option) any later version.                                                                                                            
#                                                                                                                                                                         
#   This program is distributed in the hope that it will be useful,                                                             
#   but WITHOUT ANY WARRANTY; without even the implied warranty of                                          
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                                 
#   GNU General Public License for more details.                                                                                    
#                                                                                                                                                                         
#   You should have received a copy of the GNU General Public License                                           
#   along with this program; if not, write to the                                                                                           
#   Free Software Foundation, Inc.,                                                                                                              
#   59 Temple Place, Suite 330, Boston, MA 02111-1307, USA                                                             
#                                                                                                                                                                         
#############################################################################


retlev.uvevd <- function(object, npy=1, main, xlab,
                         ylab, xlimsup, ci = TRUE, points = TRUE,
                         ...)
{
  if(!inherits(object, "uvevd"))
    stop("Use only with 'uvevd' objects")
  data <- object$data
  loc <- object$estimate["loc"]
  scale <- object$estimate["scale"]
  shape <- object$estimate["shape"]
  
  n <- length(data)
  
  gev.fun <- function(T)
  {
    p <- rp2prob(T, npy)[, "prob"]
    return(qgev(p, loc, scale, shape))
  }
  
  eps <- 10^(-3)
  
  if (!is.null(object$noy))
    npy <- n / object$noy
  else if (missing(npy)){
    warning("Argument ``npy'' is missing. Setting it to 1.")
    npy <- 1
  }
  if (missing(main)) main <- 'Return Level Plot'
  if (missing(xlab)) xlab <- 'Return Period (Years)'
  if (missing(ylab)) ylab <- 'Return Level'
  if (missing(xlimsup)) xlimsup <- prob2rp((n - .35)/n, npy)[,"retper"] 
  
  plot(gev.fun, from= 1 / npy + eps, to = xlimsup, log='x',
       xlab = xlab, ylab = ylab, main = main, ...)
  
  if (points){
    p_emp <- (1:n -.35) / n
    points(1 / ( npy * (1 - p_emp) ), sort( data ), pch = 1)
  }
  
  if (ci){
    p_emp <- (1:n - .35 ) / n
    samp <- rgev(1000*n, loc, scale, shape)
    samp <- matrix(samp, n, 1000)
    samp <- apply(samp, 2, sort)
    samp <- apply(samp, 1, sort)
    ci_inf <- samp[25,]
    ci_sup <- samp[975,]
    lines( 1 / ( npy * (1 - p_emp) ), ci_inf, lty = 2)
    lines( 1 / ( npy * (1 - p_emp) ), ci_sup, lty = 2)
  }
  
  invisible(gev.fun)
}