Do you need large amounts of data to estimate parameters in extreme value distributions?

Question

There is probably not a hard answer for this, but I am wondering if you need to collect more data when trying to estimate the parameters of generalized pareto distribution well?

The reason I ask is because I am trying to estimate the parameters of a generalized pareto distribution using Bayesian estimation, and my parameter estimates seem to be very good when I have lots of data (say, 1000+ data points), but when I drop the data size to say 100 then the estimates can be very poor.

For example, if I have a generalized pareto distribution with true parameters $\mu=0$, $\sigma=1.2$, and $\xi=0.8$, and I sample $N=1000$ observations then (running my Bayesian algorithm) I get estimates of $\hat\sigma=1.27$ 95% CI: $(1.12, 1.46)$ and $\hat\xi=0.83$ 95% CI: $(0.72, 0.98)$.

However, if I drop the same size to $N=100$ I get $\hat\sigma=0.87$ 95% CI: $(0.55, 1.24)$ and $\hat\xi=0.94$ 95% CI: $(0.61, 1.38)$.

If I keep decreasing $N$, the point estimates of $\sigma$ and $\xi$ only get worst. Is there a rule of thumb about how much data is need for extreme value distributions? In most cases, 100 data points would be sufficient for modeling a distribution that doesn't have too extreme of values (say, normal, exponential, gamma, etc.). In my application I will always be dealing with less than 100 data points and so is it a bad idea to use the generalized pareto distribution?

Here is an example of code doing what I am trying to explain:

# log-likelihood
likelihood <- function(x, xi, sigma){
  
  llik <- -log(sigma) - (1 / xi + 1) * log(1 + xi * x / sigma)
  lik  <- sum(llik)
  
  return(lik)
}

# log(prior)
prior <- function(xi, sigma){
  
  prior1 <- dgamma(xi, .01, .01, log = TRUE)
  prior2 <- dgamma(sigma, .01, .01, log = TRUE)
  prior <- (prior1 + prior2)
  
  return(prior)
}

# log(posterior)
posterior <- function(x, xi, sigma){
  
  post <- likelihood(x, xi, sigma) + prior(xi, sigma)
  
  return(post)
}


##############################################################
### Function to simulate data from GPD
##############################################################

gpd <- function(n, mu, sigma, xi){
  
  u <- runif(n)
  x = mu + sigma * (u^-xi - 1) / xi
  
  return(x)
}

set.seed(4)
N = 1000 # Number of data points
x = gpd(N, 0, 1.2, .8)  # Here mu = 0, sigma = 1.2, and xi = 0.8



S <- 10000
xi <- rep(NA, S)
sigma <- rep(NA, S)
xi[1] <- 1
sigma[1] <- 1

for(i in 2:S){
  
  
  # MCMC for xi
  xi.star = xi[i-1] + rnorm(1,0)
  
  if(xi.star < 0){
    alpha = 0
  }else{
    ratio <- exp(posterior(x, xi.star, sigma[i-1]) - posterior(x, xi[i-1], sigma[i-1]))
    alpha <- min(1, ratio)
  }

  if(runif(1) < alpha){
    xi[i] <- xi.star
  }else{
    xi[i] <- xi[i - 1]
  }
  
  
  # MCMC for sigma
  sigma.star = xi[i-1] + rnorm(1,0)
  
  if(sigma.star < 0){
    alpha = 0
  }else{
    ratio <- exp(posterior(x, xi[i-1], sigma.star) - posterior(x, xi[i-1], sigma[i-1]))
    alpha <- min(1, ratio)
  }
  
  if(runif(1) < alpha){
    sigma[i] <- sigma.star
  }else{
    sigma[i] <- sigma[i - 1]
  }  
}


sigma <- sigma[5000:S]
xi <- xi[5000:S]

mean(sigma)
mean(xi)

Answer 1

It's good to have more data , always :) However, consider why we have EVT: to work with less data! Why would you need EVT if you could collect infinite amount of data? You'd simply fit the underlying distribution and calculate any metrics on it. Because only a fraction of data goes to tails, we'd need to collect enormous amount of data before we get something going in the tails. That's where EVT comes handy: it focuses on the tails. So it allows us to study tails with much smaller data sets than otherwise would be required

Answer 2

The Fisher information matrix tells you how much information there is in each observed value about your parameters. If your observations are independent, then the information in $n$ samples is $n$ times the Fisher information matrix. The inverse of the Fisher information matrix is a lower bound on the covariance of your (unbiased) estimate (Cramer-Rao bound). So if you know how accurately you want to measure your parameters, you can invert that and divide by the elements of the Fisher Information to get a rough estimate for $n$. If your estimators are not efficient you may need more.

There's an R package mle.tools for calculating Fisher information - I've not looked to see if it handles the generalised Pareto distribution, but if not, it should at least give you a starting point for some references. Or if you have the log-likelihood, the hessian() function in package numDeriv may help.

As a rule, extreme value distributions are not necessarily any harder to estimate. It depends instead on how varying the parameters changes the shape of the distribution. If varying the parameters affects the tails but leave the central part virtually unchanged, then you need data from the tails to get a good estimate. But if the parameters change the shape of the central part as well, the bulk of the information comes from here. If you're interested, you can investigate that by considering $f(x,p) \log(f(x,p+\epsilon)/f(x,p))$ where $f(x,p)$ is the pdf at $x$ for parameter $p$, subject to a small perturbation $\epsilon$. The term $\log(f(x,p+\epsilon)/f(x,p))$ is proportional to the information in the observation of an outcome $x$ about the parameter $p$, and multiplying by $f(x,p)$ gives you the average per observation. It tells you where in the distribution the information about that parameter is typically coming from.