What likelihood to use to model sample means from a Pareto-like distribution?

Suppose there is a random variable with Lomax (Pareto Type II) probability density

$$P(x; c) = \frac{c}{(1 + x )^{c + 1}}, \quad x \ge 0, c > 0.$$

Let's draw n_samples=30000 samples of length sample_len=500 and compute mean in each sample. For a distribution with finite mean and variance, a distribution of sample means should be close to normal due to a central limit theorem. But for Lomax with $$c \le 2$$ variance is not finite, and sample means distribution is not normal (see a plot below):

Are there any likelihoods that can be used to model such distribution of central means?

Looks like an inverse-gamma model can provide a decent fit (see the plot), but is there any theoretical justification for this? I've seen there is a generalized central limit theorem that applies to a sum of Pareto random variables, but have not seen any simple closed-form expressions for a limiting distribution.

Thanks!

Code to reproduce the plot:

import numpy as np
import scipy.stats as stats
import plotly.graph_objects as go

np.random.seed(7)

c = 1.7
sample_len = 500
n_samples = 30000

exact_dist = stats.lomax(c=c)
samp = exact_dist.rvs(size=(n_samples, sample_len))
means = np.array([x.mean() for x in samp])

clt_like_mean = samp.mean()
clt_like_stdev = means.std()

invgamma_alpha = (means.mean()**2 / means.std()**2 + 2)
invgamma_beta = means.mean() * (means.mean()**2 / means.std()**2 + 1)
#arbitrary corrections for better fit
invgamma_beta = invgamma_beta / 2
invgamma_loc = 0.7

x = np.linspace(0.001, 30, 10000)

fig = go.Figure()
mode='lines', line_dash='dash', name='Original Distribution'))
mode='lines', name='CLT-like Normal'))
mode='lines', name='InvGamma'))
fig.update_layout(title='Sample Means Distribution',
xaxis_title='x',
yaxis_title='Prob Density',
hovermode="x",
height=550)
fig.update_layout(xaxis_range=[0, 5])
fig.show()

• The density* of sums of IID Pareto random variates is discussed in Ramsay, C. M. (2008), The Distribution of Sums of I.I.D. Pareto Random Variables with Arbitrary Shape Parameter, Communications in Statistics - Theory and Methods, Volume 37, Issue 14, p 2177-2184 (from which an answer to your question could be obtained) ... but a more useful question would be the frame-challenge why take the average at all?; that's not a good way to summarize the information in a Pareto. . . . * (NB do not call it the likelihood) Commented Oct 29, 2022 at 8:31
• I've tried to derive something similar (Deriving the limiting distribution of a sum of Pareto distributed variables) and for $\alpha = 0.5$ it seems to approach a Levy distribution. I guess that for other values $\alpha < 2$ the average of a sample will also approach an alpha-stable function. I stopped at that point, but I guess you could use a series representation of those distributions. Commented Oct 29, 2022 at 9:11
• Another related question where a sum of Pareto variables occurs: Expectation of the ratio of sum (XY) and sum(X) In that question a sum of the inverse of uniform variables occurs, which relates to Pareto variables with $\alpha = 1$ and the related limiting distribution is the Landau distribution. Commented Oct 29, 2022 at 9:37
• Yes, "why take the average at all?" is a valid concern :) Indeed, I have omitted most of the context from the question. Probably I'll try to post a separate question with a more detailed description of my problem. Thanks! @Glen_b Commented Oct 30, 2022 at 3:05
• @SextusEmpiricus, Glen_b, thanks for the links! I'll take a look! Commented Oct 30, 2022 at 3:07