Baseline Methods for the Parameter Estimation of the Generalized Pareto Distribution

Here are functions for estimating the parameters of the Generalized Pareto distribution using Bayesian strategies. These methods are shown in Baseline Methods for the Parameter Estimation of the Generalized Pareto Distribution (https://doi.org/10.3390/e24020178).

Generalized Pareto distribution

Density (dGPD), distribution function (pGPD), quantile function (qGPD) and random generation (rGPD) for the Generalized Pareto distribution with shape parameter equal to xi and scale parameter equal to sigma.

dGPD <- function(x,xi=0,sigma=1){
  s <- x/sigma
  if(xi==0) out <- exp(-s)/sigma else 
    out <- (1+xi*s)^(-1/xi-1)/sigma
  if(xi>=0) out <- out*as.numeric(s>=0) else 
    out <- out*as.numeric((s>=0)&&(s<=-1/xi))
  out
}

pGPD <- function(q,xi=0,sigma=1,lower.tail=TRUE){
  s <- q/sigma
  if(xi==0) out <- 1-exp(-s) else 
    out <- 1-(1+xi*s)^(-1/xi)
  if(xi>=0) out <- out*as.numeric(s>=0) else 
    out <- out*as.numeric((s>=0)&&(s<=-1/xi))
  if(!lower.tail) out <- 1-out
  out
}

qGPD <- function(p,xi=0,sigma=1,lower.tail=TRUE){
  if(!lower.tail) p <- 1-p
  if(xi==0) out <- -sigma*log(1-p) else
    out <- -sigma/xi*(1+(1-p)^(-xi))
}

rGPD <- function(n=1,xi=0,sigma=1){
  ifelse(xi==0,-sigma*log(runif(n)),
         -sigma/xi*(1+runif(n)^(-xi)))
}

Arguments

x, q: vector of quantiles.

p: vector of probabilities.

n: number of observations. If length(n) > 1, the length is taken to be the number required.

xi: shape parameter.

sigma: scale parameter.

lower.tail: logical; if TRUE (default), probabilities are \(P(X\leq x)\) otherwise, \(P(X > x)\).

Bayesian estimation strategies

Here are functions to estimate the parameters of the Generalized Pareto distribution. Firstly, using the Metropolis-Hastings algorithm, and secondly, employing the strategy called Informative Priors Baseline Distribution for different base distributions.

Metropolis-Hastings algorithm

When the shape parameter is equalt to 0, the Generalized Pareto distribution becomes the exponential distribution. Therefore, Metropolis-Hastings algorithm is provided for both the exponential distribution (MHExp) and GPD with a non-zero shape parameter (MHGPDTail).

MHExp <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  u <- as.numeric(quantile(data,p))
  Xu <- data[data>u]-u
  m <- length(Xu)
  sigma <- 1/rgamma(ndraw,1+m,1+sum(Xu))
  print(paste("Scale = ",mean(sigma),sep=" "))
  list(threshold=u,
       excess=Xu+u,
       MCMC=list(scale=sigma))
}

MHGPDTail <- function(data,shape,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  
  u <- as.numeric(quantile(data,p))
  Xu <- data[data>u]-u
  m <- length(Xu)
  
  if(shape=="light"){
    MXu <- max(Xu)
    k0 <- 0.15; delta0 <- sd(data)/k0
    
    a1 <- 0.1; b1 <- 0.1; sdk <- 0.1
    a2 <- 0.1; b2 <- 0.1; sddelta <- 0.1
    
    ki <- k0; deltai <- delta0
    
    for(i in 1:nburn){
      
      kc <- rnorm(1,ki,sdk)
      while(kc<0) kc <- rnorm(1,ki,sdk)
      deltac <- rnorm(1,deltai,sddelta)
      while(deltac<MXu) deltac <- rnorm(1,deltai,sddelta)
      
      rk <- (a1-m-1)*log(kc/ki) + b1*(ki-kc) + (1/kc-1/ki)*sum(log(1-Xu/deltai))
      
      rdelta <- (a2-m-1)*log(deltac/deltai) + b2*(deltai-deltac) + 
        (1/ki-1)*sum(log(1-Xu/deltac)) - (1/ki-1)*sum(log(1-Xu/deltai))
      
      v <- log(runif(2))
      if(v[1]<rk) ki <- kc
      if(v[2]<rdelta) deltai <- deltac
      
    }
    draw <- matrix(c(-ki,ki*deltai),ncol=2,nrow=ndraw,byrow=TRUE)
    
    for(i in 2:ndraw){
      for(j in 1:nthin){
        
        kc <- rnorm(1,ki,sdk)
        while(kc<0) kc <- rnorm(1,ki,sdk)
        deltac <- rnorm(1,deltai,sddelta)
        while(deltac<MXu) deltac <- rnorm(1,deltai,sddelta)
        
        rk <- (a1-m-1)*log(kc/ki) + b1*(ki-kc) + (1/kc-1/ki)*sum(log(1-Xu/deltai))
        
        rdelta <- (a2-m-1)*log(deltac/deltai) + b2*(deltai-deltac) + 
          (1/ki-1)*sum(log(1-Xu/deltac)) - (1/ki-1)*sum(log(1-Xu/deltai))
        
        v <- log(runif(2))
        if(v[1]<rk) ki <- kc
        if(v[2]<rdelta) deltai <- deltac
        
      }
      draw[i,] <- c(-ki,ki*deltai)
    }
  } else {
    
    xi0 <- 1; sigma0 <- 100
    
    alpha <- 10^(-6); cc <- 10^(-3); sdxi <- 0.1
    a1 <- 0.1; b1 <- 0.1; sdsigma <- 5
    
    xii <- xi0; sigmai <- sigma0
    
    for(i in 1:nburn){
      
      xic <- rnorm(1,xii,sdxi)
      while(xic<cc) xic <- rnorm(1,xii,sdxi)
      sigmac <- rnorm(1,sigmai,sdsigma)
      while(sigmac<0) sigmac <- rnorm(1,sigmai,sdsigma)
      
      rxi <- (alpha+1)*log(xii/xic) + (1+1/xii)*sum(log(1+xii*Xu/sigmai)) - 
        (1+1/xic)*sum(log(1+xic*Xu/sigmai))
      
      rsigma <- (m+a1+1)*log(sigmai/sigmac) + b1*(1/sigmai-1/sigmac) + 
        (1+1/xii)*sum(log(1+xii*Xu/sigmai)) - (1+1/xii)*sum(log(1+xii*Xu/sigmac))
      
      
      v <- log(runif(2))
      if(v[1]<rxi) xii <- xic
      if(v[2]<rsigma) sigmai <- sigmac
      
    }
    draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
    
    for(i in 2:ndraw){
      for(j in 1:nthin){
        
        xic <- rnorm(1,xii,sdxi)
        while(xic<cc) xic <- rnorm(1,xii,sdxi)
        sigmac <- rnorm(1,sigmai,sdsigma)
        while(sigmac<0) sigmac <- rnorm(1,sigmai,sdsigma)
        
        rxi <- (alpha+1)*log(xii/xic) + (1+1/xii)*sum(log(1+xii*Xu/sigmai)) - 
          (1+1/xic)*sum(log(1+xic*Xu/sigmai))
        
        rsigma <- (m+a1+1)*log(sigmai/sigmac) + b1*(1/sigmai-1/sigmac) + 
          (1+1/xii)*sum(log(1+xii*Xu/sigmai)) - (1+1/xii)*sum(log(1+xii*Xu/sigmac))
        
        
        v <- log(runif(2))
        if(v[1]<rxi) xii <- xic
        if(v[2]<rsigma) sigmai <- sigmac
        
      }
      draw[i,] <- c(xii,sigmai)
    }
  }
  
  print(paste("Shape = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
  list(threshold=u,
       excess=Xu+u,
       MCMC=list(shape=draw[,1],scale=draw[,2]))
}

MHGPD <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  
  shape <- readline(prompt = "Indicate whether you want to estimate an exponential tail (exp), a light tail (light) or a heavy tail (heavy): ")
  
  if(shape=="exp"){
    MHExp(data,p,nthin,ndraw,nburn)
  } else{
    MHGPDTail(data,shape,p,nthin,ndraw,nburn)
  }
}

Arguments

data: a vector containing the observations.

p: order of the quantile representing the threshold. By default, it is 0.9

nthin: length of the chain for thinning.

ndraw: length of the chain for the parameters to be estimated.

nburn: length of the chain for training.

shape: character; if “exp”, the MHExp function is used to estimate the scale parameter of the exponential distribution, “light”, the MHGPDTail function is used to estimate the shape and scale parameters of the GPD where shape < 0 and “heavy” when shape > 0.

Value

Return a list object.

threshold: a number representing the p-th quantile.

excess: a vector containing the excess observations over the threshold.

MCMC: a list containing the MCMC chain for the estimated parameters.

It displays on the screen the mean of the chain for each parameter.

Informative Priors Baseline Distribution algorithm

Hence, Informative Priors Baseline Distribution algorithm is shown for different base distributions.

Stable distributions

IPBDMStable <- function(data,dist,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  
  n <- length(data)
  u <- as.numeric(quantile(data,p))
  Xu <- data[data>u] - u
  m <- length(Xu)
  
  dist <- readline(prompt = "Indicate whether you want to estimate an normal distribution (norm), Cauchy distribution (cauchy) or Lévy distribution (levy): ")

  switch(dist,
         "norm" = {
           MXu <- max(Xu)
           muxi0 <- -0.7 + 0.61*p; sdxi0 <- 0.03
           musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*sd(data) 
           c1 <- -46.24; c2 <- 83.55; c3 <- -41.58
         },
         "cauchy" = {
           muxi0 <- 1; sdxi0 <- 0.065
           c1 <- 323.57; c2 <- -588.51; c3 <- 266.13
           a <- 1; b <- 1
           scalei <- as.numeric((quantile(data,3/4)-quantile(data,1/4))/2)
           musigma0 <- 1/(pi*(1-p))*scalei
           sdscale <- 0.1
         },
         "levy" = {
           muxi0 <- 2; sdxi0 <- 0.1
           c1 <- 500.2; c2 <- -900.9; c3 <- 408.2
           scalei <- rgamma(1,1+n/2,1+sum(1/data)/2)
           musigma0 <- 4/(pi*(1-p)^2)*scalei
         }
  )
  
  sdsigma0 <- exp(c1*p^2+c2*p+c3)
  sdxi <- 0.1; sdsigma <- 0.1
  
  xii <- rnorm(1,muxi0,sdxi0)
  sigmai <- rnorm(1,musigma0,sdsigma0)
  if(dist=="norm"){
    while(MXu>=(-sigmai/xii)||sigmai<=0||xii>=0){
      xii <- rnorm(1,muxi0,sdxi0)
      sigmai <- rnorm(1,musigma0,sdsigma0)
    }
  } else {
    while(xii<0) xii <- rnorm(1,muxi0,sdxi0)
    while(sigmai<=0) sigmai <- rnorm(1,musigma0,sdsigma0)
  }
  
  for(i in 1:nburn){
    
    xic <- rnorm(1,xii,sdxi)
    sigmac <- rnorm(1,sigmai,sdsigma)
    
    if(dist=="norm"){
      while(MXu>=(-sigmac/xic)||sigmac<=0||xic>=0){
        xic <- rnorm(1,xii,sdxi)
        sigmac <- rnorm(1,sigmai,sdsigma)
      }
    } else {
      while(xic<0) xic <- rnorm(1,xii,sdxi)
      while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
    }
    
    switch (dist,
            "norm" = {
              scale <- sqrt(1/rgamma(1,1+n/2,1+sum(data^2)/2))
              musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*scale
            },
            "cauchy" = {
              scale <- rnorm(1,scalei,sdscale)
              while(scale<=0) scale <- rnorm(1,scalei,sdscale)
              rscale <- (a-n-1)*log(scale/scalei) - b*(scale-scalei) - 
                sum(log(1+data^2/scale^2))+sum(log(1+data^2/scalei^2))
              if(log(runif(1))<rscale) scalei <- scale
              musigma0 <- 1/(pi*(1-p))*scalei
            },
            "levy" = {
              scale <- rgamma(1,1+n/2,1+sum(1/data)/2)
              musigma0 <- 4/(pi*(1-p)^2)*scale
            }
    )
    
    if(dist=="norm"){
      
      r <- m*log(sigmai/sigmac) + 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 +
        0.5*((sigmai-musigma0)^2 - (sigmac-musigma0)^2)/sdsigma0^2 - 
        (1+1/xic)*sum(log(1+xic*Xu/sigmac)) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
      
      if(log(runif(1))<r){
        xii <- xic
        sigmai <- sigmac
      }
    } else {
      
      rxi <- 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 - 
        (1+1/xic)*sum(log(1+xic*Xu/sigmai)) + 
        (1+1/xii)*sum(log(1+xii*Xu/sigmai))
      
      rsigma <- m*log(sigmai/sigmac) + 0.5*(sigmai - musigma0)^2/sdsigma0^2 - 
        0.5*(sigmac-musigma0)^2/sdsigma0^2 - 
        (1+1/xii)*sum(log(1+xii*Xu/sigmac)) +
        (1+1/xii)*sum(log(1+xii*Xu/sigmai))
      
      v <- log(runif(2))
      if(v[1]<rxi) xii <- xic
      if(v[2]<rsigma) sigmai <- sigmac
      
    }
    
  }
  
  draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
  
  for(i in 2:ndraw){
    
    for(j in 1:nthin){
      
      xic <- rnorm(1,xii,sdxi)
      sigmac <- rnorm(1,sigmai,sdsigma)
      
      if(dist=="norm"){
        while(MXu>=(-sigmac/xic)||sigmac<=0||xic>=0){
          xic <- rnorm(1,xii,sdxi)
          sigmac <- rnorm(1,sigmai,sdsigma)
        }
      } else {
        while(xic<0) xic <- rnorm(1,xii,sdxi)
        while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
      }
      
      switch (dist,
              "norm" = {
                scale <- sqrt(1/rgamma(1,1+n/2,1+sum(data^2)/2))
                musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*scale
              },
              "cauchy" = {
                scale <- rnorm(1,scalei,sdscale)
                while(scale<=0) scale <- rnorm(1,scalei,sdscale)
                rscale <- (a-n-1)*log(scale/scalei) - b*(scale-scalei) - 
                  sum(log(1+data^2/scale^2))+sum(log(1+data^2/scalei^2))
                if(log(runif(1))<rscale) scalei <- scale
                musigma0 <- 1/(pi*(1-p))*scalei
              },
              "levy" = {
                scale <- rgamma(1,1+n/2,1+sum(1/data)/2)
                musigma0 <- 4/(pi*(1-p)^2)*scale
              }
      )
      
      if(dist=="norm"){
        
        r <- m*log(sigmai/sigmac) + 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 +
          0.5*((sigmai-musigma0)^2 - (sigmac-musigma0)^2)/sdsigma0^2 - 
          (1+1/xic)*sum(log(1+xic*Xu/sigmac)) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
        
        if(log(runif(1))<r){
          xii <- xic
          sigmai <- sigmac
        }
      } else {
        
        rxi <- 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 - 
          (1+1/xic)*sum(log(1+xic*Xu/sigmai)) + 
          (1+1/xii)*sum(log(1+xii*Xu/sigmai))
        
        rsigma <- m*log(sigmai/sigmac) + 0.5*(sigmai - musigma0)^2/sdsigma0^2 - 
          0.5*(sigmac-musigma0)^2/sdsigma0^2 - 
          (1+1/xii)*sum(log(1+xii*Xu/sigmac)) +
          (1+1/xii)*sum(log(1+xii*Xu/sigmai))
        
        v <- log(runif(2))
        if(v[1]<rxi) xii <- xic
        if(v[2]<rsigma) sigmai <- sigmac
        
      }
      
    }
    
    draw[i,] <- c(xii,sigmai)
  }
  
  print(paste("Shape = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
  list(threshold=u,
       excess=Xu+u,
       MCMC=list(shape=draw[,1],scale=draw[,2]))
}

Exponential distribution

IPBDMExp <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  
  n <- length(data)
  u <- as.numeric(quantile(data,p))
  Xu <- data[data>u] - u
  m <- length(Xu)
  
  sigmai <- mean(data)
  
  for(i in 1:nburn){
    lambdai <- rgamma(1,1+n,1+sum(data))
    sigmac <- 1/lambdai
    rsigma <- m*log(sigmai/sigmac) + sum(Xu)*(1/sigmai-1/sigmac)
    if(log(runif(1))<rsigma) sigmai <- sigmac
  }
  draw <- sigmai 
  
  for(i in 2:ndraw){
    for(j in 1:nthin){
      lambdai <- rgamma(1,1+n,1+sum(data))
      sigmac <- 1/lambdai
      rsigma <- m*log(sigmai/sigmac) + sum(Xu)*(1/sigmai-1/sigmac)
      if(log(runif(1))<rsigma) sigmai <- sigmac
    }
    draw <- c(draw,sigmai)
  }
  
  print(paste("Scale = ",mean(draw),sep=" "))
  list(threshold=u,
       excess=Xu+u,
       MCMC=list(scale=draw))
  
}

Gamma distribution

IPBDMGamma <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
  
  n <- length(data)
  u <- as.numeric(quantile(data,p))
  Xu <- data[data>u] - u
  m <- length(Xu)
  
  alphai <- mean(data)^2/var(data); sdalpha <- 0.5
  betai <- mean(data)/var(data); sdbeta <- 0.5
  xii <- -0.032+0.014/alphai
  sigmai <- 1/betai*(0.5+0.5*sqrt(alphai))
  
  for(i in 1:nburn){
    alphac <- rnorm(1,alphai,sdalpha)
    betac <- rnorm(1,betai,sdbeta)
    while(alphac<0||betac<0){
      alphac <- rnorm(1,alphai,sdalpha)
      betac <- rnorm(1,betai,sdbeta)
    }
    rB <- -n*log(gamma(alphac)) + (n*alphac-1)*log(betac) + 
      (alphac-1)*sum(log(data)) - betac*sum(data) + 
      0.5*log(alphac*trigamma(alphac)-1) + 
      n*log(gamma(alphai)) - (n*alphai-1)*log(betai) - 
      (alphai-1)*sum(log(data)) + betai*sum(data) -
      0.5*log(alphai*trigamma(alphai)-1)
    if(log(runif(1))<rB){
      alphai <- alphac
      betai <- betac
    }
    xic <- -0.032+0.014/alphai
    sigmac <- 1/betai*(0.5+0.5*sqrt(alphai))
    rIP <- -m*log(sigmac) - (1+1/xic)*sum(log(1+xic*Xu/sigmac)) + 
      m*log(sigmai) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
    if(log(runif(1))<rIP){
      xii <- xic
      sigmai <- sigmac
    }
  }
  draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
  
  for(i in 2:ndraw){
    for(j in 1:nthin){
      alphac <- rnorm(1,alphai,sdalpha)
      betac <- rnorm(1,betai,sdbeta)
      while(alphac<0||betac<0){
        alphac <- rnorm(1,alphai,sdalpha)
        betac <- rnorm(1,betai,sdbeta)
      }
      rB <- -n*log(gamma(alphac)) + (n*alphac-1)*log(betac) + 
        (alphac-1)*sum(log(data)) - betac*sum(data) + 
        0.5*log(alphac*trigamma(alphac)-1) + 
        n*log(gamma(alphai)) - (n*alphai-1)*log(betai) - 
        (alphai-1)*sum(log(data)) + betai*sum(data) -
        0.5*log(alphai*trigamma(alphai)-1)
      if(log(runif(1))<rB){
        alphai <- alphac
        betai <- betac
      }
      xic <- -0.032+0.014/alphai
      sigmac <- 1/betai*(0.5+0.5*sqrt(alphai))
      rIP <- -m*log(sigmac) - (1+1/xic)*sum(log(1+xic*Xu/sigmac)) + 
        m*log(sigmai) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
      if(log(runif(1))<rIP){
        xii <- xic
        sigmai <- sigmac
      }
    }
    draw[i,] <- c(xii,sigmai)
  }
  
  print(paste("Shape = ",mean(draw[1,])," Scale = ",mean(draw[,2]),sep=" "))
  list(threshold=u,
       excess=Xu+u,
       MCMC=list(shape=draw[,1],scale=draw[,2]))
  
}

Arguments

data: a vector containing the observations.

p: order of the quantile representing the threshold. By default, it is 0.9

nthin: length of the chain for thinning.

ndraw: length of the chain for the parameters to be estimated.

nburn: length of the chain for training.

dist: character; if “norm”, the IPBD algorithm is used to estimate the shape and scale parameter of the GPD when the base distribution is a normal distribution, “cauchy”, when the base distribution is a Cauchy distribution and “levy” when it is a Lévy distribution.

Value

Return a list object.

threshold: a number representing the p-th quantile.

excess: a vector containing the excess observations over the threshold.

MCMC: a list containing the MCMC chain for the estimated parameters.

It displays on the screen the mean of the chain for each parameter.