6
$\begingroup$

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)
$\endgroup$
1
  • $\begingroup$ The shape shape $\xi = 0.8$ corresponds to a very heavy-tailed distribution (with no variance), which is rarely used in practice. Also you can use the revdbayes package to efficiently compute/estimate the posterior, $\endgroup$
    – Yves
    Commented Nov 9, 2021 at 7:09

2 Answers 2

8
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ Thanks! That's a very good explanation. Perhaps I need not worry then about whether or not the GPD is applicable for data sets of less than 100? $\endgroup$
    – John Smith
    Commented Oct 15, 2021 at 23:52
  • $\begingroup$ What is the alternative? If you can collect the data from the body of the distribution, then you could jointly or piece-wise estimate the full distribution. In many cases EVT is used when it is not possible to collect all data but tails due to censoring. Consider the floods: it's impossible to collect data on every flood (your bathtub overflowing is a flood too), but the biggest ones. So, in many cases EVT or something that focuses on tails is the only way to go. $\endgroup$
    – Aksakal
    Commented Oct 16, 2021 at 15:10
  • $\begingroup$ your bathtub overflowing is a flood too :) $\endgroup$ Commented Oct 16, 2021 at 16:48
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.