How to know if my data fits Pareto distribution?

Question

I have a sample which is a vector with 220 numbers. Here is a link to a histogram of my data.. And I wish to check if my data fits a Pareto distribution, but I don't want to see QQ plots with that distribution, but I need an exact answer with p-value in R, such as Anderson-Darling test for normality (ad.test). How could I do that? Please be as specific as you can.

Answer 1

(PS) First of all I think Glen_b is right in his above comments on the usefulness of such a test: real data are surely not exactly Pareto distributed, and for most practical applications the question would be "how good is the Pareto approximation?" – and the QQ plot is a good way to show the quality of such an approximation.

Any way you can do your test with the Kolmogorov-Smirnov statistic, after estimating the parameters by maximum likelihood. This parameter estimation prevents to use the $p$-value from ks.test, so you can do parametric bootstrap to estimate it. As Glen_b tells in the comment, this can be connected to Lilliefors test.

Here are a few lines of R code.

First define the basic functions to deal with Pareto distributions.

# distribution, cdf, quantile and random functions for Pareto distributions
dpareto <- function(x, xm, alpha) ifelse(x > xm , alpha*xm**alpha/(x**(alpha+1)), 0)
ppareto <- function(q, xm, alpha) ifelse(q > xm , 1 - (xm/q)**alpha, 0 )
qpareto <- function(p, xm, alpha) ifelse(p < 0 | p > 1, NaN, xm*(1-p)**(-1/alpha))
rpareto <- function(n, xm, alpha) qpareto(runif(n), xm, alpha)

The following function computes the MLE of the parameters (justifications in Wikipedia).

pareto.mle <- function(x)
{
  xm <- min(x)
  alpha <- length(x)/(sum(log(x))-length(x)*log(xm))
  return( list(xm = xm, alpha = alpha))
}

And this functions compute the KS statistic, and uses parametric bootstrap to estimate the $p$-value.

pareto.test <- function(x, B = 1e3)
{
  a <- pareto.mle(x)

  # KS statistic
  D <- ks.test(x, function(q) ppareto(q, a$xm, a$alpha))$statistic

  # estimating p value with parametric bootstrap
  B <- 1e5
  n <- length(x)
  emp.D <- numeric(B)
  for(b in 1:B)
  {
    xx <- rpareto(n, a$xm, a$alpha);
    aa <- pareto.mle(xx)
    emp.D[b] <- ks.test(xx, function(q) ppareto(q, aa$xm, aa$alpha))$statistic
  }

  return(list(xm = a$xm, alpha = a$alpha, D = D, p = sum(emp.D > D)/B))
}

Now, for example, a sample coming from a Pareto distribution:

> # generating 100 values from Pareto distribution
> x <- rpareto(100, 0.5, 2)
> pareto.test(x)
$xm
[1] 0.5007593

$alpha
[1] 2.080203

$D
         D 
0.06020594 

$p
[1] 0.69787

...and from a $\chi^2(2)$:

> # generating 100 values from chi square distribution
> x <- rchisq(100, df=2)
> pareto.test(x)
$xm
[1] 0.01015107

$alpha
[1] 0.2116619

$D
        D 
0.4002694 

$p
[1] 0

Note that I do not claim that this test is unbiased: when the sample is small, some bias can exist. The parametric bootstrap doesn’t take well into account the uncertainty on the parameter estimation (think to what would happen when using this strategy to test naively if the mean of some normal variable with unknown variance is zero).

PS Wikipedia says a few words about this. Here are two other questions for which a similar strategy was suggested: Goodness of fit test for a mixture, goodness of fit test for a gamma distribution.