# If the best-fitting distribution has infinite variance, should low observed variance be troubling?

Suppose you have observations which, over the observed range of outcomes, are well-fitted by some distribution like the Pareto that, for certain parameter values, has a an infinite variance. For example, the famous (but not especially accurate) "80/20 law" implies a tail index of about 1.16, well into the infinite-variance range.

However, the sample variance that you observe is always finite. Moreover, if you observe such a process for a long period of time, with a large number of draws in each period, it will often happen that the sample variance, whether of the period means or the underlying observations, does not seem to vary wildly or to increase without limit as sample size increases. Tail means for income give an example where this is true.

Is the some level, either of the stability of the observed variance or of its value that implies you should reject your model? As compared to infinity, how low is to low?

I guess it depends on what you mean by "best fitting." I would point out two things. First, your sample variance will always be finite, regardless of how much sampling that you do. It is invalid as a tool if the shape parameter is less than two, but the sample variance will always be finite. Second, it may be better to think of variance as undefined rather than infinite. If you think of variance as a property that is either present or absent, rather than infinite or finite, it may make thinking about it simpler.

So, let us consider the simplest possible test, which would be a Bayesian test using a conjugate prior. This may not be the best test, but it is pretty simple. It also has nice properties if you use an honest, prior density.

Let's define our likelihood as $$L(\alpha|x_1,\dots,x_n)\propto\frac{\alpha^n\beta^{na}}{m^{\alpha+1}},$$ where $$\min(x_i)>\beta$$ and $$m=\prod_{i=1}^nx_i.$$

If you use a Gamma distribution as your prior distribution, then your prior density is:

$$\pi(\alpha|a,b)=\frac{\alpha^{a-1}\exp\left(\frac{-\alpha}{b}\right)}{\Gamma(a)b^a},\text{ where }\alpha>0$$ and $a$ and $b$ are hyperparameters.

The posterior density is also a gamma distribution with parameters $a'$ and $b'$ where $$a'=a+n$$ where n is the sample size and $$b'=\frac{1}{\frac{1}{b}+\ln(m)-n\ln(\beta)}$$ and $$b<(n\ln(\beta)-\ln(m))^{-1}.$$

From:

Arnold, Barry C., and Press, S. James. ”Bayesian Estimation and Prediction for Pareto Data”, Journal of the American Statistical Association 1989, 84, 1079-1084 quoted in https://www.johndcook.com/CompendiumOfConjugatePriors.pdf

I would use wolframalpha.com if you do not have access to a math package as the cumulative posterior density is a generalized, regularized, incomplete gamma function.

I would set $a=1$ and either set $b$ equal to the magnitude you think the observations have to the minimum value or at least small enough to meet the regularity condition required above. I would then integrate up to $\alpha<1$ to determine the probability a mean existed and $\alpha<2$ to determine if the variance is defined. You do have to choose $\beta$ yourself, but it should just be the smallest positive value that you believe is the true minimum. By setting $a=1$ you have only given the prior distribution the equivalent of one observation worth of weight.

EDIT

Is the some level, either of the stability of the observed variance or of its value that implies you should reject your model? As compared to infinity, how low is to low?

If you perform the above Bayesian method and end up with undefined variance, your question is "should I reject the model?" The answer is "no." In and of itself, the answer is no. The stability would trouble me on subsamples, but what I would do is test an alternate model. Bayesian methods allow you formally test which model has greater or lesser probability of being the near the data generating function.

I have had a similar experience with sample variance in another set that definitely has no defined variance. The variance was wild until an extreme value was observed and it was so extreme that it stabilized the sample variance from that point onward. The one observation was so extreme that every wild swing ignoring it was so small that it had the appearance of stabilizing the sample variance.

I had calculated the sample variance out of curiousity and I had a similar reaction that you had. I knew as a matter of mathematics a variance didn't exist, yet it hit a date and stabilized. Once I removed the offending observations, there were actually a few of them, then the sample variance was all over the place, but far smaller.

The offending observations, which I kept for parameter estimation, were over 1500 times the interquartile range. That meant that those handful of observations were over 1500 times the maximum value of the middle 25,000,000 observations in my data.

You could have things going on in your data, such as multiple Pareto distributions or some other effect that you have not accounted for. Also draw out the gamma distribution from the posterior density. If there is a significant reason to doubt the value of the parameter, it will show up in a very wide posterior density. The posterior accounts for all uncertainty in the data, presuming your likelihood function is correct. There may be a sizable mass for the parameter in the region that has a variance.

If you have a mixture distribution, or it does not follow a Pareto distribution, then there isn't a simple or easy solution. It requires a lot of domain knowledge to even begin talking about what things may be happening in the data. How you model it depends on the domain knowledge and not at all on the statistics.

• If I understand this correctly, & I am not sure I do since I don't have much Bayesian training, it doesn't address my question. Say I have annual income data for multiple years. I fit distributions to a sample, maybe using Ben's technique above. I get an answer like 1.16, so there is no variance, or infinite variance in the limit. But I still have observed variance within sample (finite, but very high} and between yearly sample means, which is not so high -- not so high as observe when I simulate of the mean of 126,000 independent draws from a 1.16 Pareto. So should I reject Pareto model? – andrewH Feb 5 '18 at 0:21
• (126,000 is roughly the top 0.1 percent of households. – andrewH Feb 5 '18 at 0:22
• @andrewH I posted a response. – Dave Harris Feb 5 '18 at 2:20
• @andrewH I just noticed that you can draw the Gamma and the cumulative Gamma in Microsoft Excel. Look to see how much mass is in the region with defined variance. Start alpha very small and run it to something around 3. – Dave Harris Feb 5 '18 at 2:26

It is possible to bypass any issue with the variance and look directly at the plausibility of the Pareto distribution by forming ancillary statistics and then making a QQ-plot of their distribution. If we have $X_1, ..., X_n \sim \text{IID Pareto} (x_{\text{min}}, \alpha)$ then $\ln X_1, ..., \ln X_n \sim \ln x_{\text{min}} + \text{IID Exp} (\alpha)$. Using the order statistics $X_{(1)} \leqslant ... \leqslant X_{(n)}$ we can form ancillary statistics using a variation on Renyi's reprentation. We have:

$$A_k \equiv (n-k) (\ln X_{(k)} - \ln X_{(k+1)}) \sim \text{IID Exp} (\alpha).$$

Taking $E_1, ..., E_{n-1} \sim \text{IID Exp} (1)$ we therefore have the distributional equivalence:

$$\frac{A_k}{\sum_{i=1}^{n-1} A_i} \sim \frac{E_k}{\sum_{i=1}^{n-1} E_i} \sim \text{Beta} (1, n-2).$$

This means that we can use our Pareto data to obtain a vector of ancillary statistics $A_{1}, ..., A_{n-1}$ that can be compared to the simulated proportions of independent exponential random variables with unit parameter. If they fit this distribution then the Pareto fits.

For small shape parameter values ($\alpha \leqslant 1$) the variance of a Pareto distribution is infinite and the expected value of the sample variance is also infinite. However, the quantiles of the distribution of the sample variance would still all be finite, and so it should be possible to compare the observed sample variance to its quantile function to see if it is "implausibly low". The distribution and quantile functions for the sample variance do not have a simple form and so you would need to simulate them, but this should be pretty simple.

An alternative approach would be to check the sample variance of some transform with a known distribution having finite variance. For example, if $X_1, ..., X_n \sim \text{IID Pareto} (1, \alpha)$ then we have $\ln X_1, ..., \ln X_n \sim \text{IID Exp} (\alpha)$. With this transformation the true variance is the square of the true mean, and so it should be possible to check to see if the sample variance is implausibly low. To get the distribution of the sample variance you could appeal to the CLT in this case, since you now have finite variance.

• Suppose the Pareto fits the annual data for some parameter value, but the variance of the annual means is more stable than you would expect from true independent draws from a Pareto distribution with that value? I think this is true of real income data. – andrewH Feb 4 '18 at 23:41
• It depends what you mean when you say it "fits the data". Presumably having a variance that differs from the theoretical variance constitutes a lack of fit, so you would need to be a little more specific on how you are assessing that the Pareto distribution "fits the data" despite this problem. – Ben - Reinstate Monica Feb 5 '18 at 0:18