# Deriving the limiting distribution of a sum of Pareto distributed variables

For a series of independent and identical Pareto distributed variables $$X_i$$ with $$\alpha > 2$$, their sum $$S_n = \sum_{i=1}^{n} X_i$$ has a normal distribution as limiting distribution for $$n\to \infty$$

$$\frac{1}{\sqrt{n}\sigma_X} S_n -\mu_X \quad \sim \quad N(0,1)$$

But what is the situation when we have a shape parameter $$\alpha \leq 2$$? Are there some series of constants $$a_n$$ and $$b_n$$ such that the following scaled and translated sum approaches a distribution?

$$a_n S_n + b_n \quad \sim \quad ?$$

Currently I am thinking about trying to derive that it must be a stable distribution by using the characteristic function for $$a_n S_n + b_n$$ (for simplicity I set the scale parameter $$x_m =1$$).

$$\begin{array}{} \varphi_{a_n S_n + b_n}(t) & =& e^{it\,b_n} \alpha^n (-it \, a_n)^{n\alpha} \Gamma(-\alpha, -it \, a_n)^n \end{array}$$

For the $$\alpha > 2$$ case we would scale by $$a_n = \sigma_X \,n^{-0.5}$$, and for the $$\alpha \leq 2$$ case we will, I guess/suspect, need something like $$a_n \propto n^{-1/\alpha}$$.

• In order to get a more intuitive feel for this, I started to look for a derivation of that characteristic function myselve. But it seems not so easy. I have posted a question about it on math.stackexchange.com/questions/4309671 Nov 18 '21 at 14:16

### Partial results

• Below is a trial by comparing the sum of Pareto variables (with $$\alpha = 0.5$$) with a Levy distribution. The shifting and scaling are done based on the median and interquartile range. The convergence is not very fast but it does seem to work.

### function to get a scaled sample mean
require(actuar)
sum_sample = function(n,s) {
a_n = n^(-1/s)
b_n = 1*(0.5)^(-1/s)
sum( rpareto1(n,s,1)-b_n ) * a_n
}

### get a sample
set.seed(1)
x <- replicate(10^4,sum_sample(10^5,0.5))

### scale the sample according to median and interquartile range
x = x-median(x)
iqr = diff(quantile(x, probs = c(0.25,0.75)))
x = x*(rmutil::qlevy(0.75)-rmutil::qlevy(0.25))/iqr

### histogram
hist(x, freq = 0, breaks = c(seq(-5,5, 0.1),max(x)),
xlim = c(-4,4), main = "histogram of shifted and scaled sample means \n compared with levy distribution")
### levy curve
xs = seq(0.001,10,0.01)
lines(xs-rmutil::qlevy(0.5),rmutil::dlevy(xs),col=2)