# Asymptotic normality of order statistic of heavy tailed distributions

Background: I have a sample which I want to model with a heavy tailed distribution. I have some extreme values, such that the spread of the observations are relatively large. My idea was to model this with a generalized Pareto distribution, and so I have done. Now, the 0.975 quantile of my empirical data (about 100 datapoints) is lower than the 0.975 quantile of the Generalized Pareto distribution that I fitted to my data. Now, I thought, is there some way to check if this difference is something to worry about?

We know that the asymptotic distribution of the quantiles are given as:

So I thought that it would be a good idea to entertain my curiosity by trying to plot the 95% confidence bands around the 0.975 quantile of a generalized Pareto distribution with the same parameters as I got from the fitting of my data.

As you see, we are working with some extreme values here. And since the spread is so enormous, the density function has extremely small values, making the confidence bands go to the order of $\pm 10^{12}$ using the variance of the asymptotic normality formula above:

$\pm 1.96\frac{0.975*0.025}{n({f_{GPD}(q_{0.975})})^2}$

So, this does not make any sense. I have a distribution with only positive outcomes, and the confidence intervals include negative values. So something is going on here. If I calculate the bands around the 0.5 quantile, the bands are not that huge, but still huge.

I proceed to see how this goes with another distribution, namely the $\mathcal{N}(1,1)$ distribution. Simulate $n=100$ observations from a $\mathcal{N}(1,1)$ distribution, and check if the quantiles are within the confidence bands. I do this 10000 times to see the proportions of the 0.975/0.5 quantiles of the simulated observations that are within the confidence bands.

    ################################################
# Test at the 0.975 quantile
################################################

#normal(1,1)

#find 0.975 quantile
q_norm<-qnorm(0.975, mean=1, sd=1)
#find density value at 97.5 quantile:
f_norm<-dnorm(q_norm, mean=1, sd=1)
#confidence bands absolute value:
band=1.96*sqrt((0.975*0.025)/(100*(f_norm)^2))
u=q_norm+band
l=q_norm-band

hit<-1:10000
for(i in 1:10000){
d<-rnorm(n=100, mean=1, sd=1)
dq<-quantile(d, probs=0.975)

if(dq[[1]]>=l & dq[[1]]<=u) {hit[i]=1} else {hit[i]=0}

}
sum(hit)/10000

#################################################################3
# Test at the 0.5 quantile
#################################################################
#using lower quantile:

#normal(1,1)

#find 0.7 quantile
q_norm<-qnorm(0.7, mean=1, sd=1)
#find density value at 0.7 quantile:
f_norm<-dnorm(q_norm, mean=1, sd=1)
#confidence bands absolute value:
band=1.96*sqrt((0.7*0.3)/(100*(f_norm)^2))
u=q_norm+band
l=q_norm-band

hit<-1:10000
for(i in 1:10000){
d<-rnorm(n=100, mean=1, sd=1)
dq<-quantile(d, probs=0.7)

if(dq[[1]]>=l & dq[[1]]<=u) {hit[i]=1} else {hit[i]=0}

}
sum(hit)/10000


EDIT: I fixed the code, and both quantiles gives approximately 95% hits with n=100 and with $\sigma=1$. If I crank up the standard deviation to $\sigma=2$, then very few hits are within the bands. So question still stands.

EDIT2: I retract what I claimed in the first EDIT above, as pointed out in the comments by a helpful gentleman. It actually looks like these CI's are good for the normal distribution.

Is this asymptotic normality of the order statistic just a very bad measure to use, if one wants to check if some observed quantile is probable given a certain candidate distribution?

Intuitively, it seems to me like there is a relationship between the variance of the distribution (which one thinks created the data, or in my R example, which we know created the data) and the number of observations. If you have 1000 observations and an enormous variance, these bands are bad. If one has 1000 observations and a small variance, these bands would maybe make sense.

Anybody care to clear this up for me?

• Your band is based on the variance of the asymptotic normal distribution, but should be based on the standard deviation of the asymptotic normal distribution (band=1.96*sqrt((0.975*0.025)/(100*(f_norm)^2)), and similarly for the generalized Pareto dist'n.) Try that instead and see what happens. Commented Dec 1, 2015 at 18:08
• @jbowman thank you for pointing that out! I'll fix it ! Commented Dec 1, 2015 at 20:49
• @jbowman that makes the band smaller and in the example with my R-code that actually gives a few less hits. It was another error as well, that made the computation wrong, but I fixed that now. You led me on to it, so I appreciate that very much! Smaller bands in the case of the GDP is very good news, but I'm afraid they are still so enormous they are impossible to use. I still cannot see any other takeaway than that the relationship sample size and variance is what should be large, not sample size alone. Commented Dec 2, 2015 at 8:10
• No worries! I note you've correctly got a $\sqrt(n)$ in front of your first formula; if you divide both sides by that, as in band = 1.96*sqrt((0.975*0.025)/(100*n*(f_norm)^2)), that may help. Sorry I missed that the first time through. (Maybe you fixed this too but haven't updated the relevant parts of the question.) Commented Dec 2, 2015 at 15:04
• Yes it does, I didn't pay attention. OTOH, when I run your code, changing sd=1 to sd=2 everywhere, I get almost exactly the same fraction of hits both times at the 0.975 quantile: 0.9683 and 0.9662 respectively. I wonder if you missed an sd=1 somewhere in the $\sigma = 2$ run? Commented Dec 3, 2015 at 17:57

I have a distribution with only positive outcomes, and the confidence intervals include negative values.

Well, given the normal approximation that makes sense. There is nothing stopping a normal approximation from giving you negative values, which is why it is a bad approximation for a bounded value when the sample size is small and/or the variance is large. If you crank up the sample size, then the intervals will shrink because the sample size is in the denominator of the expression for the width of the interval. The variance enters the problem through the density: for the same mean, a higher variance will have a different density, higher at the margins and lower near the center. A lower density means a wider confidence interval because the density is in the denominator of the expression. How the effects of changing sample size and variance together affect the width of the confidence interval and the quality of the approximation will depend on the distribution generating the data as well as the particular quantile.

A bit of googling found this page, among others, which uses the normal approximation to the binomial distribution to construct the confidence limits. The basic idea is that each observation falls below the quantile with probability q, so that the distribution is binomial. When the sample size is sufficiently large (that's important), the binomial distribution is well approximated by a normal distribution with mean $nq$ and variance $nq(1-q)$. So the lower confidence limit will have index $j = nq - 1.96 \sqrt{nq(1-q)}$, and the upper confidence limit will have index $k = nq - 1.96 \sqrt{nq(1-q)}$. There's a possibility that either $k > n$ or $j < 1$ when working with quantiles near the edge, and the reference I found is silent on that. I chose to just treat the maximum or minimum as the relevant value.

In the following re-write of your code I constructed the confidence limit on the empirical data and tested to see if the theoretical quantile falls inside of that. That makes more sense to me, because the quantile of the observed data set is the random variable. The coverage for n > 1000 is ~ 0.95. For n = 100 it is worse at 0.85, but that's to be expected for quantiles near the tails with small sample sizes.

#find 0.975 quantile
q <- 0.975
q_norm <- qnorm(q, mean=1, sd=1)

#confidence bands absolute value (note depends on sample size)
n <- 10000
band <- 1.96 * sqrt(n * q * (1 - q))

hit<-1:10000
for(i in 1:10000){
d<-sort(rnorm(n, mean=1, sd=1))
dq<-quantile(d, probs=q)
u <- ceiling(n * q + band)
l <- ceiling(n * q - band)
if (u > n) u = n
if (l < 1) l = 1
if(q_norm>=d[l] & q_norm<=d[u]) {hit[i]=1} else {hit[i]=0}

}
sum(hit)/10000


As far as determining what sample size is "big enough", well, bigger is better. Whether any particular sample is "big enough" depends strongly on the problem at hand, and how fussy you are about things like the coverage of your confidence limits.

• Thank you for contributing! I pointed out that I do not see how there exists any absolute "large" sample, and one has to account for the variance. Im curious to how this relates to my way of constructing the CI's, but also in general. As to the derivation, you can for example see here: math.mcgill.ca/~dstephens/OldCourses/556-2006/… The CI's I've constructed follows from the example in that link. You write that "I constructed the confidence limit on the empirical data..." and this makes more sense to you. Can you please elaborate a bit more on this irt my CI's? Commented Dec 1, 2015 at 14:34
• Ah, yes, you had the right derivation link. Sorry, my bad. Commented Dec 1, 2015 at 14:44
• OK, I edited it again to correctly describe how the variance of the distribution affects the approximation you're using, and a bit more discussion on what a "large" sample means. Your CI is centered on the theoretical value, while mine is centered on the empirical one. I think for comparing an empirical quantile with a theoretical one the intervals should be constructed on the empirical quantile. Also the approximation I used makes one less "normal" approximation because there is no appeal to the central limit theorem to start. Commented Dec 1, 2015 at 15:00
• I appreciate the effort, maybe my question could be clearer. I have already realized how the density and sample size affects the variance, that was my point in the first place. But, again, my bad, I could have been more clear. It is the "asymptotic" that I feel should be switched out with something that takes the variance into account. Well, you have also centered your CI's around the theoretical values. n*q is exactly your theoretical value. In constructing your bands, you have essentially done the same thing as I, only with a different method. Commented Dec 1, 2015 at 15:25