Arthur Charpentier
email: arthur.charpentier@uqam.ca

Emmanuel Flachaire
email: emmanuel.flachaire@univ-amu.fr

Pareto Models for Top Incomes and Wealth by

Keywords: C46, D31 Pareto distribution, top incomes, inequality measures

Top incomes are often related to Pareto distribution. To date, economists have mostly used Pareto Type I distribution to model the upper tail of income and wealth distribution. It is a parametric distribution, with interesting properties, that can be easily linked to economic theory. In this paper, we first show that modeling top incomes with Pareto Type I distribution can lead to biased estimation of inequality, even with millions of observations. Then, we show that the Generalized Pareto distribution and, even more, the Extended Pareto distribution, are much less sensitive to the choice of the threshold. Thus, they can provide more reliable results. We discuss different types of bias that could be encountered in empirical studies and, we provide some guidance for practice. To illustrate, two applications are investigated, on the distribution of income in South Africa in 2012 and on the distribution of wealth in the United States in 2013.

Introduction

Income and wealth distributions are skewed to the right, with thick upper tails. Thus, microdata samples often exhibit outliers, from which sample inequality indices can be severly distorted, even in large sample (Cowell and Victoria-Feser 1996). The upper tail of income distribution is then often modelled with a parametric distribution, such as the Pareto distribution.

In his initial work on income and wealth distributions, Vilfredo Pareto suggested several parametric models (Pareto 1895, 1896). The term "Pareto distribution" refers to both Pareto I and Generalized Pareto distributions. Rytgaard (1990) wrote that the Pareto I distribution is the commonly used definition of the Pareto distribution in Europe, and the Generalized Pareto distribution in America. To some extent, this distinction could also be made between economists and statisticians.

To date, economists have mainly used the Pareto I distribution to fit the upper tail of income distribution with microdata (e.g. Davidson and Flachaire 2007, Cowell and Flachaire 2007, Burkhauser et al. 2011, Alfons et al. 2013, Hlasny and Verme 2016, Vermeulen 2018, Higgins et al. 2018, Lustig 2018) and with tabulated data (e.g. Atkinson, 2017, Atkinson et al. 2011, Lakner and Milanovic 2015). 1 The density, the cumulative distribution function and the quantile function are simple power functions, with a single parameter to estimate, which can be related to theoretical models that can explain the generation of his thick upper tail (Gabaix 2009, Jones 2015, Benhabib and Bisin 2018).

Recently, Jenkins (2017) addressed the question: What model should be fitted to top incomes? In an empirical application, based on tax data in U. K. over 1995-2011, he finds that the preferred specification is the Generalized Pareto distribution (GPD). In this paper, we provide theoretical foundation to this result, and we go further.

Our contribution is manifold. We first show that the Pareto I distribution is very sensitive to the choice of the threshold. In particular, a threshold too low can lead to over-estimate the heaviness of the distribution and, thus, to over-estimate inequality. This bias comes from a misspecification, it does not disappear as the sample size increases. Next, we show that the GPD is less sensitive to the threshold, but its estimation is less accurate. We also show that the Pareto I behaves like the GPD only at (much) higher threshold. Then, we introduce the Extended Pareto distribution (EPD), which is even less sensitive to the threshold and which provides more reliable results. Finally, we discuss different types of bias that could lead to under-or over-estimate inequality in practice, and we illustrate our findings through two applications, on income distribution in South-Africa in 2012 and on wealth distribution in the U.S. in 2013.

The importance of choosing the right threshold for fitting a model is widely acknowledged in Extreme Value Theory (Embrechts et al. 1997, Beirlant et al. 2004). Two approaches are often used to address this problem: by improving the estimation or by using more flexible model. Related to the first approach, the literature has accumulated many different estimation methods for the tail index. Fedotenkov (2018) provides a comprehensive review of more than one hundred Pareto-tail index estimators, from which there is no single best choice. The current paper is related to the second approach. Obviously, using a more flexible model is an easy way to reduce bias, but it comes at the cost of higher variance. Our results show that using the GPD may provide more variance, but much less with the EPD, which is also more flexible.

An important implication of our results, discussed in the conclusion, is to reconsider two common beliefs in empirical studies of income inequality. Indeed, it is widely believed that inequality measures are: (1) under-estimated from surveys, because these data are subject to topcoding, censoring and underreporting of the rich ; (2) more reliable from tax data, because there is no topcoding or censoring and these data are much less sensitive to misreporting. Our results show that inequality can be strongly over-estimated from surveys and from tax data, if the selected threshold is too low.

In section 2, we present Pareto I and Generalized Pareto distributions. We discuss how sensitive they are to the choice of the threshold, and we introduce the Extended Pareto distribution. In section 3, we conduct several simulation experiments. In section 4, we discuss different type of biases. In section 5, two applications are investigated. Section 6 concludes.

Pareto models

Pareto models have been often used for modelling the upper tail of distributions in economic inequality and economic losses, in finance and insurance (Albrecher et al. 2017, Beirlant et al. 2004, Kleiber and Kotz 2003). In this section, we first present strict Pareto distributions, that is, Pareto Type I and Generalized Pareto distributions. We discuss the close relationship between them and their sensitive to the choice of the threshold. Then, we present a Pareto-type distribution, namely the Extended Pareto distribution, from which Pareto Type I and GPD are two special cases. In the following, we consider distributions with finite mean, that is, with Pareto tail parameter greater than one (α > 1). For more details on properties of Pareto distributions, see Arnold (2015).

Pareto I distribution

For Pareto Type I distribution bounded from below by u > 0, with tail parameter α, the probability density function and the cumulative density function (CDF) are, respectively, equal to

f (x) = αu α x α+1 and F (x) = 1 - x u -α , for x ≥ u (1)
If a random variable X has (1) as its CDF, we will write X ∼ P 1 (u, α).

In the economic inequality literature, Pareto Type I distribution has been always used to fit the upper tail of income and wealth distributions. This distribution has an attractive property: the average above a threshold is proportional to the threshold, and it does not depend on the scale parameter u,

E(X|X > u ) = αu α -1 , α > 1 (2)
where u ≥ u. For instance, if the inverted Pareto coefficient α/(α -1) = 2, the average income of individuals with income above $100,000 is $200,000 and, the average income of individuals with income above $1 million is $2 million (Piketty 2007, Atkinson et al. 2011). This nice property makes easy calculation of the mean and leads to simple formulas of inequality measures (see appendix B). Moreover, simple theoretical economic models can be used to explain the basic mechanism that gives rise to Pareto distributions (Jones 2015).

Generalized Pareto distribution

For Generalized Pareto Distribution (GPD) with a tail parameter α > 0, a convenient parametrization of the cumulative density function (CDF) is provided by

F (x) = 1 -1 + x -u σ -α for x ≥ u (3)
where σ, the scale parameter, is positive, and u is a lower-bound.2 It is the form also known as the Pareto Type II distribution in the terminology of Arnold (2015). If a random variable X has (3) as its CDF, we will write X ∼ GPD(u, σ, α). The GPD(0, σ, α) distribution is also called "Lomax" in the literature (Lomax 1954).

From ( 1) and (3), we can see that Pareto I is used to model the distribution of relative excesses X/u given X > u, while GPD is used to model the distribution of absolute excesses X -u given X > u. And because of the additional parameter, there is an affine transformation that links these relative and absolute excesses. The GPD is "generalized" in the sense that Pareto I distribution is a special case, when σ = u:

GPD(u, u, α) = P 1 (u, α) (4)
GPD lets the model decide if the upper tail is better modeled as a power law from relative excesses (σ = u) or from absolute excesses (σ = u). Overall, GPD is more flexible than Pareto I distribution.

The average above a higher threshold, u ≥ u, depends on all parameters of the distribution,

E(X|X > u ) = σ -u α -1 + α α -1 u , (5)
The linearity of this function characterizes the GPD class (see Guess and Proschan 1988 and Ghosh and Resnick 2010).

In the statistical literature, GPD is often used to fit the upper tail of heavy-tailed distributions (Beirlant et al. 2004). Indeed, one of the most important result in the extreme value theory states that, for most heavytailed distributions, the conditional excess distribution function, above a threshold u, converges towards a GPD distribution as u goes to infinity (Pickands 1975, Balkema and de Haan 1974), for some parameters α and σ,

F u (x) -→ GPD (or Pareto II) as u → +∞ (6)
where

F u (x) = P (X -u ≤ x|X > u)
. This result is known as the Pickands-Balkema-de Haan theorem, also called the second theorem in extreme value theory (as discussed in footnote 33, page 28). It provides strong theoretical support for modelling the upper tail of heavy-tailed distributions with GPD, or Pareto Type II distribution. In a very general setting, it means that there are α and σ such that F u can be approximated by the CDF of a GPD(u, σ, α), see Embrechts et al. (1997, Theorem 3.4.13, p.165).

Threshold selection

Estimation of the upper tail of a distribution with a Pareto distribution proceeds with a preliminary choice of the threshold u by the researcher. In statistics, the choice of the threshold is notoriously difficult, and the results are often quite sensitive to this choice in empirical studies.

The choice of the threshold means that only the k largest observations are considered, and the distribution of those k observations is supposed to be Pareto, with tail index α. There is a bias-variance trade-off: if the threshold is too high, k is small and estimators can be very volatile, while if the threshold is too small, there is less variability, but the bias can be large. In practice, an optimal threshold would then be the lowest threshold that does not generate significant bias in parameter estimates.

Whether to use a Pareto I or GPD model to fit the upper tail is related to the choice of the threshold. From ( 1) and (3), we have seen that Pareto I is a special case of GPD, when σ = u. They differ by an affine transformation when σ = u. A key property of Pareto distributions is that, if a distribution is Pareto for a fixed threshold u, it is also Pareto with the same tail parameter α for a higher threshold u ≥ u. For GPD, we have

F u = GPD(u, σ, α) ⇒ F u = GPD(u , σ + u -u, α). ( 7
)
where F u is the survival excess function above u. 3 Thus, Pareto I and GPD are the same for all u ≥ u, if σ = u. Otherwise, we have σ + u -u ≈ u for very large values of u only. It follows that a GPD above a threshold will behave approximately as a Pareto I above a higher threshold, much higher as σ differs from u.

To illustrate, Figure 1 shows boxplots of the maximum likelihood estimator (MLE) of the tail index α, from Pareto I (left) and GPD (right) models, as the threshold increases. The values are obtained from 1,000 samples of 1,000 observations drawn from a GPD where σ differs from u, GPD(0.5, 1.5, 2).

For Pareto I model (left), the MLE of the tail index is the Hill estimator,

α k = 1 k n i=n-k+1 log x (i) -log x (n-k+1) -1 (8)
3 F u is a truncated Pareto distribution, with density equals to f (x)/(1 -F (u )). Note that this property is quite intuitive, since the GPD distribution appears as a limit for exceeding distributions, and limit in asymptotic results are always fixed points: the Gaussian family is stable by addition (and appears in the Central Limit Theorem) while Fréchet distribution is max-stable (and appears in the first theorem in extreme value theory).

from a sample {x 1 , • • • , x n }, with the ordered version x (1) ≤ x (2) ≤ • • • x (n) ,
where only the k largest values are considered. We can see that the boxes are always below the true value (dashed line), and the dispersion increases as the threshold increases. It shows that, if the threshold is not very high, Pareto I model provides severe bias in the estimation of the tail parameter. For instance, with a threshold u = 0.5, we have α ≈ 1, that is, the estimated tail parameter is half as small as the true value, α = 2.

For GPD model (right), we can see that there is no bias, but the dispersion is much greater. Moreover, several values are extremes, it can be more than twice the true value (circles in the plot). Here, the MLE of the tail index has no analytical solution, and numerical methods are required for the estimation. It is known that the MLE of the GPD can be very unstable and, in some cases, has no solution (Castillo and Hadi 1997, Hosking and Wallis 1987, del Castillo and Daoudi 2009).

[Figure 1 about here.] Overall, this Figure illustrates that Pareto Type I distribution behaves like a GPD only at much higher threshold, when σ differs from u. From his analysis on UK income data, Jenkins (2017) concludes that "the Pareto I model is as good as GPD only at extremely high incomes, beyond a range of thresholds usually considered". It is consistent with the values of σ -u that he found, always very different from zero.4

Extended Pareto distribution

We have seen that lower threshold can be used with GPD model, compared to Pareto I. Nevertheless, a (very) high threshold may still be required to have (nearly) unbiased estimation of the tail parameter with GPD model. To select not too high a threshold, we need to consider Pareto-type models, rather than strictly Pareto models (Pareto I and GPD). It leads us to consider higher order of regular variation, from which an Extended Pareto distribution (EPD) can be derived (see Appendix A for more details). The EPD provides better fit of the upper tail when the distribution becomes strictly Pareto in the very top of the distribution only.

The Extended Pareto Distribution (EPD), proposed by Beirlant et al. (2009), is defined as follows:

F (x) = 1 - x u 1 + δ -δ x u τ -α for x ≥ u (9)
where τ ≤ 0 and δ > max(-1, 1/τ ). If a random variable X has (9) as its CDF, we will write X ∼ EPD(u, δ, τ, α).

Pareto I is a special case when δ = 0 and GPD is a special case when τ = -1:

EPD(u, 0, τ, α) = P 1 (u, α) (10) EPD(u, δ, -1, α) = GPD(1, u/(1 + δ), α) (11)
The mean over a threshold for the EPD distribution has no closed form expression.5 Numerical methods can be used to calculate it. Since u ≥ u > 0, X given X > u is a positive random variable and

E[X|X > u ] = ∞ 0 F u (x)dx ( 12
)
where

F u (x) = P[X > x|X > u ] for x > u, i.e. F u (x) = F (x) F (u )
where F is the s.d.f. of X

Thus

E[X|X > u ] = u + 1 F (u ) ∞ u F (x)dx (14)
The integral in ( 14) can be computed numerically. Since numerical integration over a finite segment could be easier, we can make a change of variable (1/x) to obtain an integral over a finite interval:

E u = ∞ u F (x)dx = 1/u 0 1 x 2 F 1 x dx (15)
The Extended Pareto distribution has a stable tail property: if a distribution is EPD for a fixed threshold u, it is also EPD for a higher threshold u ≥ u, with the same tail parameter α. Indeed, deriving a truncated EPD distribution, we find

F u = EPD(u, δ, τ, α) ⇒ F u = EPD(u , δ , τ, α). ( 16
)
where

δ = δ(u /u) τ /[1 + δ -δ(u /u) τ ].

Simulations

In this section, we use Monte Carlo experiments to study the sensitivity of Pareto models to the choice of the threshold and to the degree of deviation from a strict-Pareto distribution in the upper tail.

Sensitivity to the threshold

Figure 2 shows boxplots of values of the tail parameter α, estimated by maximum likelihood from Pareto I (left), GPD (middle) and EPD (right) models, as the threshold increases. The x-axis is the percentage of the largest observations used to fit Pareto distributions. The values are obtained from 1,000 samples of 50,000 observations drawn from a Singh-Maddala distribution that closely mimics the US 2013 income distribution, SM (2.07, 1.14, 1.75).6

[Figure 2 about here.]

For Pareto I model (left), we can see that the income distribution is estimated too heavy, whatever the threshold is (blue boxes are always below the true value, given by the dashed line). Even when the Pareto I distribution is estimated with the top 1% observations only (highest threshold, k = 500), the box is clearly below the true value. In other words, the Hill estimator is biased downward: the lower the threshold, the more biased it is. 7 For GPD model (middle), we can see that the income distribution is estimated less heavy than it is in reality (green boxes are always above the true value). It is only for very high thresholds that the estimates/boxes are around the true value. Moreover, we can see that the dispersion is much higher than for Pareto I model.

For Extended Pareto model (right), we can see that the estimation of the tail parameter outperforms Pareto I and GPD models. Indeed, red boxes are always much closer to the true value, compared to blue boxes (Pareto I) and green boxes (GPD), for given threshold. The dispersion is slightly larger than Pareto I, and much smaller than GPD.8

[Figure 3 about here.]

It is important to note that the bias doesn't disappear as the sample size increases, since it is a misspecification problem. To illustrate, Figure 3 shows similar results from huge samples, that is, from samples of 1 million of observations. Compared to Figure 2, Figure 3 shows clearly that the dispersion decreases as the sample size increases, not the bias. Thus, having several millions of observations does not reduce the bias.

Sensitivity to deviations from a strict-Pareto

Using a Singh-Maddala (Burr) distribution allows us to control the deviation from a strict-Pareto, with the first-and second-order of the regular variation (RV), see appendix A for more details. Moreover, this distribution is often used to estimate density of income distribution in parametric approaches, see Cowell and Flachaire (2015).

In our experiments, we generate 10,000 samples of 50,000 observations from Singh-Maddala distributions, SM (-ρ, 0.5, -α/ρ), with an upper tail that deviates more and more from a strict-Pareto, that is, with an increasing second-order RV parameter ρ = -2, -1, -0.75, -0.5. For each sample, we estimate the tail index α by maximum likelihood based on Pareto I, GPD and EPD models, fitted on, respectively, the top 10%, 5% and 1% largest observations:

-Pareto 1 models fitted on, respectively, the top 10%, 5%, 1% largest observations are, respectively, denoted P1, P1', P1", -GPD models fitted on, respectively, the top 10%, 5%, 1% largest observations are, respectively, denoted GPD, GPD', GPD", -EPD models fitted on, respectively, the top 10%, 5%, 1% largest observations are, respectively, denoted EPD, EPD', EPD". Atkinson et al. (2011, p.15) observe that, in practice, the tail index α of income distributions typically vary between 1.5 and 3. When α < 2, the variance is infinite. In the following, we first consider the case of finite variance (α = 3) and the case of infinite variance (α = 1.5).

Finite variance

[Figure 4 about here.]

Figure 4 shows boxplots of the tail index estimates α, when the true value is α = 3. The case ρ = -2 (top-left) is quite similar to the case used in the previous subsection and the results have similar patterns as those in Figure 2. The other cases, where ρ increases, are used to analyze what happens when the deviation from a strict Pareto in the upper tail increases. We can see that:

• Pareto I model (Hill estimator) shows severe under-estimation of the tail index α, when ρ increases (boxes are increasingly below the horizontal dashed line). As expected by Schluter ( 2018), the distribution is estimated too heavy with a Pareto I model, see ( 33).

• GPD model exhibits less bias, but much more dispersion, when ρ increases (boxes are closer to the dashed line, but wider). Several outliers are far from the true value. It suggests that GPD can perform better than Pareto I, but it can also provide very poor results.9

• EPD model outperforms Pareto I model. The dispersion is slightly increased, but the bias is largely reduced, when ρ increases. However, the EPD model does not correct completely the bias when ρ is large (boxes are still below the horizontal dashed line).10

Infinite variance

[Figure 5 about here.]

Figure 5 shows boxplots of the tail index estimates α, when the true value is α = 1.5. The results are quite similar to the case of finite variance. It provides poor results with Pareto I model. Fitted on the top 10% of the sample (P1), with a large deviation to a strict-Pareto (ρ = -0.5), the tail index is estimated smaller than 1, which suggests that the mean doesn't exist and Pareto models cannot be used. Note that similar results have been obtained in empirical studies on wealth distribution. For instance, Cowell (2013) finds α = 0.48, 0.52, 0.73 from a Pareto I model fitted on, respectively, top 10%, 5% and 1% observations of the U.S. wealth distribution.

Overall, the simulation results show that Pareto I model can lead to severe over-estimation of the heaviness of the distribution, when the threshold is not large enough and when the distribution deviates from a strict-Pareto. At given threshold, the EPD model can provide significant improvement over Pareto I model.

From theory to practice

In empirical studies, the analysis of Pareto-type upper tail is often inferred from some plots. For Pareto Type I distribution, with tail parameter α, rearranging (1) and applying a logarithmic transformation, we have

log(1 -F (x)) = c -α log x ( 17
)
where c = α log u. Given a sample {x 1 , • • • , x n }, let us consider the plot of the log of incomes in the x-axis and the log of the survival function in the y-axis:

{log x, log (1 -F (x))} (18)
or its empirical version

log x i , log 1 -F (x i ) (19)
This plot, known as the Pareto diagram (Cowell 2011) or Zipf plot (Cirillo 2013), shows the proportion of the population with x or more against x itself on a double-logarithmic diagram. If the distribution is strictly Pareto, it should exhibit a linear function, where the slope coefficient (taken positively) is equal to the tail index α. If the distribution is not strictly Pareto, the curve obtained from a Pareto diagram is not linear.11 This interpretation can be carried over to the general case of Pareto-type distributions since then ultimately for x → ∞ the Pareto diagram is still linear with slope α at some set of largest values (Albrecher et al. 2017, p.70).

In practice, tail index estimation is subject to three potential sources of bias: misspecification bias, estimation bias and sampling bias.

Misspecification bias

Let us assume that the CDF is known, F (x).

[Figure 6 about here.]

Figure 6 shows a Pareto diagram, based on a Singh-Maddala distribution that closely mimics the US 2013 income distribution, SM (2.07, 1.14, 1.75). 6 This plot exhibits a concave curve, which becomes linear in the right part, when log x > 1. As expected, it suggests that the distribution behaves like a Pareto distribution in the upper tail. It also exhibits two linear functions, as defined in ( 17), obtained when we consider two different threshods, log u = -2, 1. Since α is both the slope of the linear function (taken positively) and the tail index of the Pareto distribution, we can see that:

• α increases as the threshold u increases, until it is constant,

• the distribution is then estimated too heavy if u is not large enough. Most Pareto diagrams based on the CDF exhibit typically a concave-like curvature, at least for large x, when the distribution is heavy-tailed.12 Thus, we may expect that distributions are often estimated too heavy with Pareto I model, if the threshold is too low. 13 This problem occurs because the linear approximation is not appropriate. It is important to note that it is a misspecification bias, which does not disappear as the sample size increases to infinity (see section 3.1).

The simulation results in

Estimation bias

In empirical studies, the CDF is in general unknown. We replace it by a consistent estimate, F (x). When we consider the estimation of α, the replacement of the CDF by an estimate introduce another type of bias: an estimation bias.

A standard choice is the EDF, which is a consistent estimator of the CDF. If we sort a sample in ascending order,

x (1) < x (2) < • • • < x (n)
, the Pareto diagram is then the log of the ordered observations in the x-axis and the log of the empirical survival function in the y-axis:14

Pareto diagram: log x (i) , log 1 - i n + 1 i = 1, . . . , n (20)
The Pareto diagram often exhibits an erratic behavior in the extreme right part. It is because the EDF provides a poor estimation of the upper tail of the CDF, especially when the distribution is heavy-tailed.15 Indeed, sample quantiles poorly estimate higher population quantiles when the distribution is heavy-tailed. 16 To illustrate the consequences of this erratic behavior, let us consider the value of the maximum income in the sample, x (n) . The Pareto diagram tells us that the proportion of population with the maximum income or more is equal to 1/(n + 1). It follows that,

• if the maximum income is too large (Bill Gates' income belongs to the sample),17 the proportion of the population with income more or equal to x (n) is over-estimated: 1/(n + 1) > 1 -F (x (n) ). Thus, the heaviness of the upper tail is over-estimated.

• if the maximum income is too small (topcoding, censoring, truncation or underreporting),18 the proportion of the population with income more or equal to

x (n) is under-estimated: 1/(n + 1) < 1 -F (x (n) ).
Thus, the heaviness of the upper tail is under-estimated.

Figure 7 illustrates these features. A Pareto diagram is obtained from a sample of 20,000 observations drawn from the Singh-Maddala distribution SM (2.07, 1.14, 1.75). Firstly, we gradually increase the largest incomes, that is, those such that log(x) > 2. The Figure shows that the presence of outliers pushes the extreme part of the curve to the right, which becomes convex before becoming concave again. Secondly, we topcode the largest incomes, that is, every log-income greater than 2 is set equal to 2. The Figure shows that the topcoding pushes the extreme part of the curve to the left, with a more pronounced concavity.

[Figure 7 about here.]

Overall, the estimation bias can go in opposite directions. It will depend on the exact nature of the problem. The presence of outliers in the sample is known to biased downward classical estimators of the tail index (Hubert et al. 2013). Thus, it tends to estimate the upper tail too heavy. To the opposite, top-coding or underreporting from the rich tend to biased upward estimators of the tail index.19 Thus, it tends to under-estimate the heaviness of the upper tail.

In applications, the main question is which type of bias is the most pronounced and to what extent? Examining the Pareto diagrams can help to answer this question, as we will show in the application section.

Sampling bias

It is often argued that surveys do not capture well the top of income distributions, because the rich may be harder to reach or more likely to refuse to participate (Atkinson 2007, Atkinson et al. 2011). If some members of the population are more or less likely to be included than others, a sampling bias is introduced in the estimation. To correct this bias, surveys often include weights to make the sample representative of the overall population. 20 It is particularly important to use the weights when provided by data producers, otherwise the estimated distribution maybe quite different from the population distribution.

To illustrate, recall that a sample {x 1 , • • • , x n } is based on individual observations x i of individuals who agreed to respond. Let X denote the (true) variable of interest (individual income, or wealth), and D denote the response (to the survey) variable (1 if the individual responds). Let Z denote some possible covariates (age, gender, etc) and X the income. Since all computations should be conditional on Z, we will skip it to avoid too heavy notations. 21 Following Horvitz and Thompson (1952), assume that the variable of interest X has distribution F θ . The Horvitz-Tompson estimator of θ should be based on weights W = P[D = 1|X] -1 . Here, our sample of x i 's can be seen as realizations of variables X|D = 1. Using Bayes formula, and the weights, it will be possible to link the true distribution of income, and the one of our sample. Hence,

P[X|D = 1] = P[X] • P[D = 1|X] P[D = 1] ∝ P[X] =F θ • P[D = 1|X] =W -1 (21)
Consider here some sort of proportional hazard model for the weight function, with either

P[D = 1|X = x] = S x ( 22
)
where S is some baseline with value in [0, 1], or

P[D = 1|X = x] = x -a (23)
If X is Pareto distributed with tail index α, with density proportional to x -(1+α) , then, with model ( 22), density of X|D = 1 should be proportional to

x → x -(1+α) • S x (24)
and the survey distribution has a Weibull distribution. With model ( 23), density of X|D = 1 should be proportional to

x → x -(1+α) • x -a (25)
which is another Pareto distribution.

In practical applications, even if the true distribution of X is Pareto I, if the response variable is (strongly) correlated with the income, either the sample is also Pareto I, with possibly another tail index, or there might be some second order effect. Thus, to prevent this possible bias, weights can be used. With a sample {(x 1 , w 1 ), • • • , (x n , w n )}, where w i 's are given by statistical institutes providing the datasets, our estimators should take them into account. With maximum likelihood estimation, we need to rewrite the likelihood function as the product of each individual contribution multiplied by its weight. For instance, the weighted version of the Hill estimator would be

α k,w = n i=n-k+1 w (i) wk log y (i) -log y (n-k+1) -1 (26)
where wk = n i=n-k+1 w (i) . The special case w i = c, where c is a constant value, corresponds to the standard Hill estimator, as defined in (8). The estimation of GPD and EPD distributions with weights, as well as plotting Pareto diagram with weights, are not implemented in standard softwares. We develop functions in R to do this, which we make available on GitHub. 22

Applications

In this section, we consider two applications: the income distribution in South-Africa in 2012 and the wealth distribution in the United-States in 2013. We consider top share indices, which are derived for the three Pareto models in appendix B.

In practice, a plot of the estimates of the tail index αk with respect to the number of upper ordered observations k used for the estimation, also known as Hill plot, is often used to select a threshold. If the Pareto distribution is appropriate a stable horizontal straight line should be plotted. An appropriate threshold is then selected from a stable region of the graph. However, Embrechts et al. (1997, p.194) and Resnick (2007, p.87) illustrate that the Hill estimator can perform very poorly if the distribution is not strictly Pareto, i.e. not Pareto I or GPD. 23 It can lead to very volatile Hill plots, also known as Hill horror plots, from which a threshold is difficult to select. In such cases, we could expect to obtain more stable plots with an EPD-based estimator and, thus, a threshold could be more easily selected.24

Income distribution in South-Africa

The income distribution in South-Africa in 2012 is obtained from the Luxembourg Income Study (LIS) Cross-National Data Center in Luxembourg. We use similar incomes as those used to generate the Key Figures, except that we do not bottom-and top-code incomes. 25 Income is equal to disposable household income divided by the square root of the number of household members (x = DHI/ √ NHHMEM). All households where disposable income is missing or exactly equal to zero are excluded, and we use person-level ad-justed weights (w = HWGT * NHHMEM). The number of observations is equal to n = 7,990.

Figure 8 shows values of the tail index estimated by (weighted) maximum likelihood, α, as k the number of the largest observations used for Pareto estimation increases. The vertical dashed lines marked q 90 , q 95 and q 99 correspond to, respectively, the 90%, 95% and 99% weighted quantiles. 26 The Pareto I curve shows α obtained from Pareto I model, also known as Hill plot. The two other curves show α obtained from GPD and EPD models. We can see that the Pareto I curve decreases as k increases. It is never stable and horizontal, the Hill plot is then not very revealing. The GPD curve is unstable for k < 200, otherwise it increases as k increases. To the opposite, the EPD curve looks stable when more than 200 observations are used to fit the Pareto distribution, suggesting that the tail index is slightly greater than 2. These results are consistent with those obtained by simulations, in Figure 2 and3.

[Figure 8 about here.]

Figure 9 shows values of the top 1% share, as k increases. The horizontal line is the value of the sample top 1% share, given in (36). We can see that the top 1% share obtained from Pareto I model increases as k increases. Since more inequality is expected when the distribution is heavier (α is smaller), these results are consistent with those obtained in Figure 8. To the opposite, the top 1% obtained from GPD and EPD models are much more stable. The values given by the EPD model are slightly higher than those given by the GPD model.

[Figure 9 about here.]

The following Table shows values from Figures 8 and9 for three thresholds, q 90 , q 95 and q 99 , that is, when the Pareto distribution is fitted on, respectively, the top 10%, 5% and 1% observations. tail index top 1% threshold q 90 q 95 q 99 q 90 q 95 q 99 Pareto I

1 8: it does not capture decreases of the tail index from GPD model, as the threshold increases, from q 90 to q 94 .27

Figure 10 shows a Pareto diagram, with the Pareto I, GPD and EPD distributions fitted on the top 10% observations. We can see that the straight lines of the Pareto I model decreases more slowly than the curves of the GPD and EPD models, above the top 1% income (when log(x) ≥ q 99 ). The EPD model captures deviations from a strict Pareto distribution, it seems to better fit the Pareto diagram. The extreme right part exhibits an erratic behavior, but we should not pay attention to this part (see section 4.2).

[Figure 10 about here.] Overall, we can see that a threshold too low leads to under-estimate the tail index and to over-estimate the top share, with a Pareto I model. It illustrates the misspecification bias, as explained in section 4.1, which is not expected to disappear as the sample size increases. Figures ( 8), ( 9) and ( 10) suggest that the EPD model provides more reliable results. 28 The top 1% is equal to 14.8%, with an EPD distribution fitted on the top 10% observations (q 90 ), which is slightly more than the sample top 1% share equals to 13.3%.

Wealth distribution in the United-States

The wealth distribution in the United-States in 2013 is obtained from the Luxembourg Wealth Study (LWS) database. We use disposable household net worth divided by the square root of the number of household members (x = DNW/ √ NHHMEM). All households where disposable net worth is missing are excluded. We do not bottom-and top-code wealth, and we use personlevel adjusted weights (w = HWGT * NHHMEM). The number of observations is equal to n = 6,015. 29 Figure 11 shows values of the tail index estimated by (weighted) maximum likelihood, α, as k the number of the largest observations used for Pareto estimation increases. The vertical dashed lines marked q 90 , q 95 and q 99 correspond to, respectively, the 90%, 95% and 99% weighted quantiles. The Pareto I curve shows α obtained from Pareto I model, also known as Hill plot. The two other curves show α obtained from GPD and EPD models. We can see that the Pareto I curve decreases as k increases. It becomes quite stable when the threshold is below q 99 (above the top 1% observations), with a tail index slightly greater than 1.5. The EPD curve also decreases as k increases. It is not very stable when the threshold is between q 95 and q 99 , but it isn't still decreasing, suggesting a tail index slightly greater than 1.5. The GPD curve is not very different from the EPD curve, but it is slightly less stable.

[Figure 11 about here.]

Figure 12 shows values of the top 1% share, as k increases. We can see that the top 1% shares obtained from Pareto I model with a threshold greater than q 99 , and from GPD and EPD models with a threshold greater than q 95 , are very stable and very similar to the sample top 1% share (given by the horizontal line). Otherwise, the top 1% share increases as k increases, faster with Pareto I model.

[Figure 12 about here.]

The following Table shows values from Figures 11 and 12 for three thresholds, q 90 , q 95 and q 99 , that is, when the Pareto distribution is fitted on, respectively, the top 10%, 5% and 1% observations. tail index top 1% threshold q 90 q 95 q 99 q 90 q 95 q 99 Pareto I 0. We can see that the Pareto I model is very sensitive to the choice of the threshold, compared to the others. For instance, the Pareto I model with q 95 and q 99 provides very different estimates of the top 1% share (76.7% and 42.1%), while EPD model gives quite similar values (46.6% and 43.6%). Moreover, Pareto I model with q 90 gives inconsistent results: α < 1, meaning that the mean doesn't exist and, thus, the top share is undefined. If we assume that the top 1% of the population has less than 100% of income or wealth, the tail index has to be greater than one, α > 1.

Figure 13 shows a Pareto diagram, with the Pareto I, GPD and EPD models fitted on the top 3% observations (q 97 ). We can see that the GPD and EPD models provide a much better fit of the Pareto diagram, compared to the Pareto I model.

[Figure 13 about here.] Overall, we can see that the EPD model is much less sensitive to the choice of the threshold.

Finally, it is important to notice that the measurement of inequality of wealth is somewhat more challenging than the analysis of income or consumption, because most sample data include a substantial fraction of negative net worth, and the variance may not exist. 30 Negative observations make the interpretation of standard measures of inequality not so easy, see Cowell and Van Kerm (2015). For instance, the top share is not defined over the interval [0, 1], but over ] -∞, +∞[. 31 It can take any negative or positive value, even greater than one.

Conclusion

Since income and wealth distributions are often assumed to be heavy-tailed, economists have mainly used the Pareto Type I distribution, or Pareto I model, to estimate the upper tail in empirical studies.

In this paper, we first show that the Pareto I model (Hill estimator) is very sensitive to the choice of the threshold. When the threshold is too low, it can lead to severe under-estimation of the tail index and, thus, to severe over-estimation of inequality.

Then, we provide evidence that the GPD model, based on the Generalized Pareto Distribution, is less sensitive to the choice of the threshold. We show that Pareto I model is as good as the GPD model but only for higher threshold, much higher as the scale parameter differs from the threshold (σ = u). However, the estimation of the GPD distribution turns out to be quite inaccurate in our simulation results.

Next, we introduce the EPD model, based on the Extended Pareto Distribution proposed by Beirlant et al. (2009), which allows to capture deviations from a strict Pareto distribution, and from which Pareto I and GPD are two special cases. Our simulation results and two applications show that the EPD model is less sensitive to the choice of the threshold, and its estimation is more accurate than the GPD model. It suggests that the EPD model should be preferred to the GPD model in practice.

Finally, we discuss the fact that the tail index estimation is sensitive to several biases, which can go in opposite directions. To summarize, the tail index estimation of a Pareto model is mostly: (1) downward biased, when the threshold is not high enough and/or in the presence of outliers ; (2) upward biased with topcoding, censored data and when the rich underreport their income. It leads us to provide new highlights on several common beliefs.

It is widely believed that the tail index estimation of Pareto I model based on surveys is biased upward, because these data are subject to topcoding, censoring and underreporting of the rich. In a seminal paper, Atkinson et al. (2011, footnote 8, p.11) write that: "The Pareto parameter is estimated using the ratio of the top 5 percent income share to the top decile income share (. . . ). Because those top income shares are often based on survey data (and not tax data), they likely underestimate the magnitude of the changes at the very top."32 It is true, under the (strong) assumption that the distribution is strictly Pareto I above the top decile, and if there is no outlier in the sample. Otherwise, the tail index can be biased downward, and the upper tail might just as easily be over-estimated.

It is also widely believed that the tail index estimation is much more reliable with tax data, because there is no topcoding or censoring and, above all, these data are much less sensitive to misreporting. Indeed, tax data are much less sensitive, if not at all, to estimation bias. However, they are still sensitive to misspecification bias. Thus, a threshold too low may lead to substantial downward bias of the tail index and to severe over-estimation of inequality, even with millions of observations. From tax data in the U.K., for several years 1995-2010, Jenkins (2017, p.279) finds that the optimal thresholds for Pareto I model are at around the 99.5-percentile or higher, that is, well above the thresholds commonly used. This suggests that fitting Pareto models to tax data should be done with more caution than has been done so far. A First-and second-order regular variations

Pareto-type models are more flexible than strictly Pareto models. It leads us to consider higher order of regular variation, from which an Extended Pareto distribution (EPD) can be derived.

The tail index α is related to the max-domain of attraction of the underlying distribution, while parameter σ is simply a scaling parameter. 33 The shape of the conditional excess cumulative distribution function is a power function (the Pareto distribution) if the threshold is large enough. Tails are then said to be Pareto-type, and can be described using so called regularly varying functions (see Bingham et al. 1987).

First and second order regular variation were originally used in extreme value theory, to study respectively the tail behavior of a distribution and the speed of convergence of the extreme value condition (see Bingham et al. 1987, de Haan and Stadtmüller 1996, Peng and Qi 2004, or

A function regularly varying with index γ = 0 is said to be slowly varying.

Observe that any regularly varying function of index -γ can be written

H(x) = x -γ (x)
where is some slowly varying function.

Consider a random variable X, its distribution is regularly-varying with index -γ if, up-to some affine transformation, its survival function is regularly varying. Hence,

lim t→∞ x -γ F (tx) F (t) = 1 (28) or F (x) = x -γ (x), (29)
where F (x) = 1 -F (x). A regularly varying survival function is then a function that behaves like a power law function near infinity. Distributions with survival function as defined in ( 29) are called Pareto-type distributions.

It means that the survival function tends to zero at polynomial (or power) speed as x → ∞, that is, as x -γ . For instance, a Pareto I distribution, with survival function F (x) = x -α u α , is regularly varying with index -α, and the associated slowly varying function is the constant u α . And a GPD or Pareto II distribution, with survival function

F (x) = 1 + σ -1 x -α , is
regularly varying with index -α, for some slowly varying function. But in a general setting, if the distribution is not strictly Pareto, will not be constant, and it will impact the speed of convergence. In de Haan and Stadtmüller (1996), a concept of second-order regular variation function is introduced, that can be used to derive a probabilistic property using the quantile function,34 as in Beirlant et al. (2004). Following Beirlant et al. (2009), we will consider distributions such that an extended version of equation ( 28) is satisfied,

lim t→∞ x -γ F (tx) F (t) = 1 + x ρ -1 ρ , for some ρ ≤ 0, ( 30
)
that we can write, up to some affine transformation,

F (x) = cx -γ [1 -x ρ (x)], (31)
for some slowly varying function and some second-order tail coefficient ρ ≤ 0. The corresponding class of Pareto-type distributions defined in ( 31) is often named the Hall class of distributions, referring to Hall (1982). It includes the Singh-Maddala (Burr), Student, Fréchet and Cauchy distributions. A mixture of two strict Pareto-I distributions will also belong to this class. 35 Since ρ ≤ 0, x → 1 -x ρ (x) is slowly varying, and therefore, a distribution F that satisfies (31) also satisfies (29). More specifically, in (31), the parameter ρ captures the rate of convergence of the survival function to a strict Pareto distribution. Smaller is ρ, faster the upper tail behaves like a Pareto, as x increases. Overall, we can see that

• γ is the first-order of the regular variation, it measures the tail parameter of the Pareto distribution,

• ρ is the second-order of the regular variation, it measures how much the upper tail deviates from a strictly Pareto distribution.

In the following, we will write RV(-γ, ρ). There are connections between tail properties of the survival function F and the density f (see also Karamata theory for first order regular variation). More specifically, if F is RV(-γ, ρ), with γ > 1, then f is RV(-γ -1, ρ).

For instance, consider a Singh-Maddala (Burr) distribution, with survival distribution F (x) = [1 + x a ] -q , then a second order expansion yields

F (x) = x -aq [1 -qx -a + o(x -a )] as x → ∞ (32)
which is regularly varying of order -aq and with second order regular variation -a, that is RV(-aq, -a). Schluter (2018) shows that the bias, induced by higher order regular variation in the Burr case (32), for the estimator of the inverted tail parameter ξ = 1/α in Pareto I distribution, obtained from an OLS regression of log sizes on log ranks, is equal to

b k,n = 1 2 ξ 2 -ρ (1 -ρ) 2 n k ρ , ( 33
)
where the k largest observations are considered to fit a Pareto I and n is the sample size. Since ρ ≤ 0, the bias is positive, it increases as ρ and k increases. Thus, the distribution will be estimated heavier than it is really, as ρ tends to zero, and the bias can be large if the threshold is not high enough. For values of ρ close to zero, the bias of the tail parameter of Pareto I model is then expected to be negligible for extremely high threshold only. Beirlant et al. (2009) show that (31) can be approximated by

F (x) = [x (1 + δ -δx τ )] -α for x ≥ 1 (34)
where τ ≤ 0 and δ > max(-1, 1/τ ). 36 The main feature of this function is that it captures the second-order regular variation of the Hall class of distributions, that is, deviation to a strictly Pareto tail. For more details, see also Albrecher et al. (2017, section 4.2.1). From ( 34), we can define the Extended Pareto Distribution (EPD) in ( 9), as proposed by Beirlant et al. (2009).

B Lorenz curve and top shares

In this appendix, we derive inequality measures based on distributions being Pareto in the upper tail (with finite mean, α > 1).

The top p100% income share can be defined as follows:

TS p = p E(X|X > Q(1 -p)) E(X) (35)
where Q(t) is the quantile function, Q(t) = inf{x ∈ R : t ≤ F (x)}. A Lorenz curve is defined by (p, L(p)), where L(p) = 1 -TS 1-p and p ∈ (0, 1). A top income share is then related to a single value from the Lorenz curve.

From a sample of n incomes and weights, {(x 1 , w 1 ), • • • , (x n , w n )}, the sample top income share, based on the Empirical Distribution Function (EDF), is computed as:

TS (edf) p = n i=1 w i x i 1(x i > Q(1 -p)) n i=1 w i x i (36)
where 1(•) is an indicator function and Q(1 -p) is a sample weighted quantile. 37 Without weights, the sample top income share can be computed from (36), with w i = 1.

Let us consider a distribution where the upper tail, above a threshold or cut-point Q(1 -q) = u, is modelled by a Pareto distribution P, and the remaining distribution is modelled by another distribution F. For a two-component sliced distribution, with a Pareto distribution P for the top q100% and a distribution F for the bottom (1 -q)100%, we need to consider two cases. When the top p100% share is located in the upper tail modelled by the Pareto distribution, we have

TS p≤q = p E P (X|X > Q P (1 -p/q)) (1 -q)E F (X) + qE P (X) if p ≤ q (37)
When the top p100% share is located below the upper tail modelled by the Pareto distribution, we have

TS p≥q = 1 - (1 -p) E F (X|X < Q F (1 -p)) (1 -q)E F (X) + qE P (X) if p > q (38)
If F is the empirical distribution function (EDF) and if the sample is given with weights, the threshold u is a weighted quantile. 38 Pareto I and GPD models: we consider that the upper tail, above a threshold Q(1 -q) = u, is modelled by a GPD distribution, and the remaining distribution is modelled by the empirical distribution function (EDF). The quantile function of a GPD is equal to

Q GPD (t) = σ(1 -t) -1/α -σ + u (39)
Using ( 5) and ( 39) in ( 37) and ( 38), we obtain

TS (GPD) p,q =          p [α/(α -1)]σ(p/q) -1/α + u -σ (1 -q)x q + qσ/(α -1) + qu if p ≤ q 1 - (1 -p) xp (1 -q)x q + qσ/(α -1) + qu if p > q (40)
where xq (x p ) is the weighted mean of the (1-q)100% ((1-p)100%) smallest ordered observations. Top shares from Pareto I model for top incomes are obtained from ( 40), with σ = u.

EPD model: we consider that the upper tail, above a threshold Q(1q) = u, is modelled by an Extended Pareto Distribution (EPD), and the remaining distribution is modelled by the empirical distribution function 38 With n i=1 wi = 1, the cut-point u is the ordered observation x (l) satisfying l-1 i=1 w (i) ≤ 1 -q and n i=l+1 w (i) ≤ q.

(EDF). There are no closed form expressions for the quantile function and for the mean of EPD distribution. However, we can use numerical methods. The quantile function is defined as

Q EPD (t) = {x ≥ u : (x/u)[1 -δ + δ(x/u) τ ] = (1 -t) -1/α } (41)
The equation in parenthesis can be solved numerically, by findings root for the difference between the two terms, or by minimizing the square of the difference.39 From ( 14), ( 15), ( 37) and ( 38), top p100% shares from EPD model are defined as follows:40

TS (EPD) p,q =        pu + qE u (1 -q)x q + q(u + E u ) if p ≤ q 1 - (1 -p) xp (1 -q)x q + q(u + E u ) if p > q (42)
where E u and E u are obtained by numerical integration from (15), and u is obtained from (41) with t = 1 -p/q. Figure 2: Boxplots of maximum likelihood estimators of the tail index: 1,000 samples of 50,000 observations drawn from a Singh-Maddala distribution, SM (2.07, 1.14, 1.75). From the left to the right, the x-axis is the threshold (percentile) used to fit Pareto I (blue), GPD (green) and EPD (red) models. Figure 3: Boxplots of maximum likelihood estimators of the tail index: 1,000 samples of 1,000,000 observations drawn from a Singh-Maddala distribution, SM (2.07, 1.14, 1.75). From the left to the right, the x-axis is the threshold (percentile) used to fit Pareto I (blue), GPD (green) and EPD (red) models.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' ρ=-2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' ρ=-1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' ρ=-0.75 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' Figure 4: Finite variance (α = 3) -Boxplots of tail index estimates α based on 10,000 samples of 50,000 observations drawn from Singh-Maddala distributions with a upper tail that deviates more and more from a strict-Pareto (as ρ increases). Pareto I (blue), GPD (green) and EPD (red) models, fitted on, respectively, the top 10%, 5% and 1% of the sample. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD''

1.0 1.5 2.0 2.5 ρ=-2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD''

1.0 1.5 2.0 2.5 ρ=-1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' 1.0 1.5 2.0 2.5 ρ=-0.75 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q P1 P1' P1'' GPD GPD' GPD'' EPD EPD' EPD'' q 90 q 95 q 99 log(1-F(x))

Pareto 1 GPD EPD q 90 q 95 q 99 log(1-F(x))

Pareto 1 GPD EPD q 90 q 95 q 99 q 97

  Figures 2 and 3 are obtained with samples drawn from the Singh-Maddala distribution used to draw the Pareto diagram in Figure 6. The concavity of the Pareto diagram explains why the Pareto I model (Hill estimator) under-estimates the tail index and, thus, overestimates the heaviness of the distribution.

  Rytgaard, M. (1990). Estimation in the Pareto distribution. ASTIN Bulletin 20, 201-216. Schluter, C. (2018). Top incomes, heavy tails, and rank-size regressions. Econometrics -MDPI 6 (1), 1-16. Schluter, C. (2019). On Zipf's law and the bias of Zipf regression. Toda, A. A. (2012). The double power law in income distribution: Explanations and evidence. Journal of Economic Behavior & Organization 84 (1), 364 -381. Vermeulen, P. (2018). How fat is the top tail of the wealth distribution? Review of Income and Wealth 64, 357-387.

Figure 1 :

 1 Figure 1: Boxplots of maximum likelihood estimators of the tail index from Pareto I (left) and GPD (right) models, as the threshold increases: 1,000 samples of 1,000 observations drawn from a GPD distribution, GPD(0.5, 1.5, 2).

Figure 5 :

 5 Figure5: Infinite variance (α = 1.5) -Boxplots of tail index estimates α based on 10,000 samples of 50,000 observations drawn from Singh-Maddala distributions with a upper tail that deviates more and more from a strict-Pareto (as ρ increases). Pareto I model (blue), GPD model (green) and EPD (red), fitted on, respectively, the top 10%, 5% and 1% of the sample.

Figure 6 :

 6 Figure 6: Pareto diagram based on the true CDF, with two linear approximations based on log x ≥ -2 and log x ≥ 1

Figure 7 :

 7 Figure 7: Pareto diagram based on a sample (black), with an artificial increase of the largest observations (blue) and with topcoding (green).

Figure 8 :

 8 Figure 8: Incomes in South-Africa 2012: MLE estimates of the tail index, α, from Pareto I (blue), GPD (green) and EPD (red) models, as the number of the k-largest observations used for Pareto estimation increases.

Figure 9 :

 9 Figure 9: Incomes in South-Africa 2012: Top 1% income shares obtained from Pareto I (blue), GPD (green) and EPD (red) models, as the number of the k-largest observations used for Pareto estimation increases.

Figure 10 :

 10 Figure 10: Incomes in South-Africa 2012: Pareto diagram, with Pareto I (blue), GPD (green) and EPD (red) models fitted on the top 10% incomes.

Figure 11 :

 11 Figure 11: Wealth in the United States 2013: MLE estimates of the tail index, α, from Pareto I (blue), GPD (green) and EPD (red) models, as the number of the k-largest observations used for Pareto estimation increases.

Figure 12 :

 12 Figure 12: Wealth in the United States 2013: Pareto diagram, with Pareto I (blue) and EPD (red) models fitted on the top 5% incomes.

Figure 13 :

 13 Figure 13: Wealth in the United States 2013: Pareto diagram, with Pareto I (blue) and EPD (red) models fitted on the top 5% incomes.

  We can see that the GPD and EPD models provides very stable values for the top 1% shares, while the Pareto I is quite sensitive to the choice of the threshold. The tail index values for q 99 may appear odd. They are not reliable, since they use k = 26 observations only. Finally, this Table is not as informative as Figure

		.742 1.881 2.492	0.192 0.171 0.146
	GPD	2.689 2.935 19.249	0.142 0.141 0.139
	EPD	2.236 2.255 4.198	0.148 0.149 0.139

A few exceptions include Jenkins (2017) with microdata and[START_REF] Blanchet | Generalized Pareto curves: Theory and applications[END_REF] with tabulated data.

The limiting case α → ∞ corresponds to an exponential distribution.

The plots of u -σ can be found in bottom Figures in Appendix H, in the online supplementary material ofJenkins' paper,

Albrecher et al. (2017, section 4.6) give an approximation, based on (1 + δ -δy τ ) -α ≈ 1 -αδ + αδy τ , which can be very poor. Thus, we do not recommend to use it.

These parameters are obtained from the estimation of a Singh-Maddala distribution with the US incomes in 2013 provided by LIS cross-national center in Luxemburg (a = 2.07 and q = 1.75 are shape parameters, b = 1.14 is a scale parameter)

This results is expected from Schluter's result in (33).

We have also made similar experiments for two extensions of the GPD model, proposed in the literature: Pareto IV[START_REF] Arnold | Pareto Distributions[END_REF] and Extended Generalized Pareto, or EGPD3[START_REF] Papastathopoulos | Extended generalised pareto models for tail estimation[END_REF]. We do not report the results here, since the dispersion is always much larger and they are outperformed by Extended Pareto model.

See the discussion in section 2.3. Note that GPD seems to have no bias when ρ = -1 (boxes are around the dashed line). It is not surprising, since the Singh-Maddala distribution is a GPD when the shape parameter a = 1 in (32), which corresponds to ρ = -1.

This model is based on the approximation (1 + δ -δy ρ ) -α ≈ 1 -αδ + αδy ρ , which is better as ρ is more and more negative when y > 1.

Another popular plot is the Pareto QQ-plot, {-log(1 -F (xi)), log xi}. It is obtained from the tail quantile function of Pareto I distribution, Q(p) = (1 -p) -1/α . Hence, log Q(p) = -1 α log(1 -p).If the distribution is Pareto I, we expect to see a linear function with a slope coefficient equal to 1/α. Note that these plots may not be helpful to distinguish between a Lognormal and a Pareto distributions[START_REF] Cirillo | Are your data really Pareto distributed?[END_REF].

Schluter (2019, appendix A) shows formal derivations on the concavity of Pareto QQplots, with the Hall class of distributions. Using many different parameters and parametric distributions, we always found concave Pareto diagrams, based on the CDF.

 13 Beirlant et al. (2004, p.113) argue that, in many cases, the Hill estimator overestimates the population value of 1/α.

With R, we can plot a Pareto diagram with: plot(log(sort(x)),log((n:1)/(n+1))), and a Pareto QQ-plot with plot(-log((n:1)/(n+1)),log(sort(x))).

This is why[START_REF] Davidson | Asymptotic and bootstrap inference for inequality and poverty measures[END_REF] and[START_REF] Cowell | Income distribution and inequality measurement: The problem of extreme values[END_REF] proposed to fit the upper tail with a parametric Pareto distribution to improve inference.

The (asymptotic) variance of sample quantiles is equal to n-1 F (x)(1 -F (x))/f (x) 2 .For Pareto I distribution, F (x) = 1-x -α , the variance is then equal to n -1 (x α+2 -x 2 )/α 2 , which explodes as x → +∞.

More precisely, if x (n) is greater than the n/(n + 1)-percentile of the population.

More precisely, if x (n) is smaller than the n/(n + 1)-percentile of the population.

Linking a restrictive subsample of Uruguay's official household survey to tax data,[START_REF] Higgins | The rich underreport their income: Assessing bias in inequality estimates and correction methods using linked survey and tax data[END_REF] show that the rich tend to underreport their income.

To correct the sampling weights for unit nonresponse, see[START_REF] Korinek | Survey nonresponse and the distribution of income[END_REF][START_REF] Korinek | An econometric method of correcting for unit nonresponse bias in surveys[END_REF]

For instance, P[X|D = 1, Z] could be written P Z [X|D = 1], but instead of using P Z in all computations, we will simply use P

All R programs developed in this paper and a R package can be downloaded from https://github.com/freakonometrics/TopIncomes

More precisely, when the slowly varying function is not constant in (29)

The threshold selection from a stable range should not have a significant impact on the uncertainty of inequality estimates, since tail parameter estimates are quite similar over this range. But, it should have an impact if the plot does not exhibit a stable range.

see http://www.lisdatacenter.org/data-access/key-figures/methods/

Note that the 200 largest observations corresponds to the top 5% incomes (x ≥ q95), but not to 5% of observations in the sample (n = 7,990) because of the weights.

That is, the GPD curve that increases as k increases, for k > 200, in Figure8.

Figure8suggests that similar results could be obtained with Pareto I and GPD models, with a higher threshold. However, the choice of the threshold is difficult based on Pareto I and GPD models only, since the curves are never stable and horizontal.

The data provider conducted a multiple imputation procedure to impute missing values. The number of observations in the original file, equals to

30,075, is five times the actual number of households, because the imputations are stored as five successive implicates of each record (see the LWS database user guide). In our sample, we use the first imputation of each record.

In our sample, 747 observations are negative (12.4%). The variance does not exist when α < 2, which can lead to unreliable inequality estimates[START_REF] Fontanari | Gini estimation under infinite variance[END_REF]).

From (35), top share can take large positive (negative) values when the overall mean is close to zero and positive (negative). It is equal to infinity when the overall mean is 0.

For more details on the Pareto parameter calculation, seeAtkinson (2007, p.24).

Historically, extremes were studied through block-maximum -yearly maximum, or maximum of a subgroup of observations. Following[START_REF] Fisher | Limiting forms of the frequency distribution of the largest or smallest member of a sample[END_REF], up to some affine transformation, the limiting distribution of the maximum over n i.i.d observations is either Weibull (observations with a bounded support), Gumbel (infinite support, but light tails, like the exponential distribution) or Fréchet (unbounded, with heavy tails, like Pareto distribution).[START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] and[START_REF] Balkema | Residual life time at great age[END_REF] obtained further that not only the only possible limiting conditional excess distribution is GPD, but also that the distribution of the maximum on subsamples (of same size) should be Fréchet distributed, with the same tail index γ, if γ > 0. For instance in the U.S., if the distribution of maximum income per county is Fréchet with parameter γ (and if county had identical sizes), then the conditional excess distribution function of incomes above a high threshold is a GPD distribution with the same tail index γ.

The quantile function U is defined as U (x) = F -1 (1 -1/x).

Note that this second order Pareto distribution is not the double power law discussed in Toda (2012), where a power law is consider for both tails, the lower and the upper ones.

Using the expansion (1 + y a ) b ≈ 1 + by a , for small y a , in (34) yields (31).

n i=1 wi = 1, the weighted quantile Q(1 -p) is the ordered observation x (k) satisfying k-1 i=1 w (i) ≤ 1 -p and n i=k+1 w (i) ≤ p.

The function qepd of the ReIns R package minimizes the square of the difference, with u = 1.

Since the distribution is Pareto in the upper tail only, we have F (u ) = p/q.

Acknowledgments: Arthur Charpentier received financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC-2019-07077) and the AXA Research Fund, and Emmanuel Flachaire from the the French National Research Agency (ANR-17-EURE-0020, ANR-17-CE41-0007-02, ANR-19-FRAL-0006), and by the Excellence Initiative of Aix-Marseille University -A*MIDEX.