# Mind the Income Gap: Behavior of Inequality Estimators from Complex Survey
Small Samples

Silvia De Nicolò  Department of Statistical Sciences, University of Padova
Maria Rosaria Ferrante  Department of Statistical Sciences ”P. Fortunati”,
University of Bologna  Silvia Pacei  Department of Statistical Sciences ”P.
Fortunati”, University of Bologna

###### Abstract

Income inequality measures are biased in small samples leading generally to an
underestimation. After investigating the nature of the bias, we propose a
bias-correction framework for a large class of inequality measures comprising
Gini Index, Generalized Entropy and Atkinson families by accounting for
complex survey designs. The proposed methodology is based on Taylor’s
expansions and generalized linearization method, and does not require any
parametric assumption on income distribution, being very flexible. Design-
based performance evaluation of the suggested correction has been carried out
using data taken from EU-SILC survey. Results show a noticeable bias reduction
for all measures. A bootstrap variance estimation proposal and a
distributional analysis follow in order to provide a comprehensive overview of
the behavior of inequality estimators in small samples. Results about
estimators distributions show increasing positive skewness and lepto-kurtosis
at decreasing sample sizes, confirming the non-applicability of the classical
asymptotic results in small samples and suggesting the development of
alternative methods of inference.

Keywords — Bias Correction, Complex Surveys, Income Inequality, Small Sample
Inference

##  1 Introduction

The estimation of income Inequality Measures (IM) in small samples can face
several issues related to bias and non-robustness to outliers. First of all,
those measure are known to be biased in small samples [Deltas(2003)] [Breunig
and Hutchinson(2008)], leading usually to an underestimation of inequality.
The amount of bias, generally negative, depends on the magnitude of the
coefficient of variation of the underlying distribution of the variable of
interest and, for some specific measures, it also depends on its skewness
[Breunig(2001)]. Consider that income is well-known to be positively skewed.
Moreover, the magnitude of the bias varies depending on the inequality
measure. This aspect deserves attention given that inequality measures values
are used for comparisons across time and location. Neglecting the bias may
bring out discrepancies in estimated quantities due to different sample sizes
or different amounts of distribution skewness rather than be true inequality
differences [Breunig and Hutchinson(2008)]. As regards Gini index bias, a rich
literature faces it such as [Jasso(1979)], [Lerman and Yitzhaki(1989)],
[Deltas(2003)], [Davidson(2009)], [Van Ourti and Clarke(2011)] in i.i.d.
samples and [Fabrizi and Trivisano(2016)] for the complex survey case.
However, concerning alternative measures such as Atkinson Indexes and the
Generalized Entropy (GE) measures, literature on bias is very scarce even in
i.i.d. case. Some contributions are provided by [Giles(2005)] and [Schluter
and van Garderen(2009)] for GE measures and [Breunig and Hutchinson(2008)] for
GE measures and Atkinson Indexes. The mentioned references adopt different
methodological approaches to correct or reduce bias in a i.i.d. context.
However, income data are usually collected via specific household survey with
a complex sampling design, based on stratification of the survey population
and on selection of sampling units in more than one stages. Thus, survey
sample selection process, together with ex-post treatment procedures such as
calibration and imputation, invariably introduces a complex correlation
structure in the data, which has to be taken into account and it makes the
development of a theoretically valid bias correction challenging.

Furthermore, the issues related to non-robustness to outliers of inequality
measures are widely explained by [Cowell and Victoria-Feser(1996)]. The non-
robustness issue is even exacerbated in income data applications, which are
traditionally affected by the extreme values problem, see [Van Kerm(2007)]. It
has been widely tested that, generally, those measures appear to be unrobust
even to an infinitesimal amount of data contamination, especially when dealing
with extreme values on the tails. This unrobustness depends clearly on the
type of measure we are dealing with and it becomes even more cumbersome to
handle in case of small samples. Moreover, the presence of extreme values can
increase the bias. The inequality measures considered here are Gini Index,
Atkinson Index and Generalized Entropy measures, together with the Coefficient
of Variation measure.

The aim of this paper is twofold. On one hand, we provide a comprehensive
discussion about the behavior of inequality estimators in small samples from
complex surveys, including issues generally addressed in a piecemeal manner,
such as bias evaluation, robustness, variance estimation and distributional
analysis. This analysis could be very useful for future developments of
inference strategies on such measures in the context of small samples
inference, where classical asymptotic results do not hold, above all in the
field of small domain estimation. On the other hand, the aim is to propose a
methodological framework for bias-correction in a finite population setting,
more specifically taking into account complex survey designs. The methodology
is based on Taylor’s expansions, where the estimation of design-based
estimator variances and covariances has been performed via the generalized
linearization method, [Deville(1999)] and [Demnati and Rao(2004)], relying on
the concept of influence functions. Any parametric assumptions on the income
distribution is not required, providing a very flexible framework applicable
to any income distributional assumption.

Our bias-correction proposal has been evaluated via simulations showing a
noticeable bias reduction for all the measures and leading in some cases to
approximate unbiased estimators. An in-depth analysis on the most sensitive
measures confirms the great impact outliers have on the magnitude of
estimators bias and variance.

A preliminary definition of the considered inequality measure can be found in
Section 2. The bias-correction strategy is set out in Section 3 and the bias-
correction estimation steps are detailed in Section 4. A design-based
simulation study involving European Survey of Income and Living Condition (EU-
SILC) income data [Clemenceau and Museux(2007)] is provided in Section 5, in
order to evaluate bias-correction and sensitivity to extreme values. A design-
aware bootstrap for bias-corrected estimators variance is proposed in Section
6, while a distributional analysis involving an additional model-based
simulation follows in Section 7. Conclusions are drawn in Section 8.

##  2 Inequality Measures

The most famous inequality measure is, indeed, the Gini concentration index,
employed in social sciences for measuring concentration in the distribution of
a positive random variable. There are several equivalent definitions of Gini
index [Ceriani and Verme(2015)], we will use the formulation of [Sen(1973)] as
follows. Suppose we have a finite population $\mathcal{U}$ of $N(<\infty)$
elements labeled as $\\{1,\dots,N\\}$. Let $y_{i}$ be a characteristic of
interest, in our case income, for the $i$-th unit of the finite population,
where $y_{i}\in R^{+}$, $\forall i=1,\dots,N$,

$\displaystyle G=\frac{1}{N^{2}\mu}\sum_{i}N_{i}y_{i}-\frac{N+1}{N},$

with $N_{i}$ denoting the rank of the $i$-th unit and $\mu$ the expected value
of income variable.

The estimation of alternative measures, as opposed to Gini index, may enable a
more meaningful assessments of different aspects of economic inequality. Gini
does not allow to decompose inequality into within groups and between groups
components, moreover it is positional (weakly) transfer sensitive, namely
variations in the index depends on the ranks of the donor and recipients. And
lastly Gini constitutes itself as a stochastic dominance measure, based on
partial ordering of probability distributions: two very different
distributions - one having more inequality amongst the poor, the other amongst
the rich can have exactly the same Gini Index value. When the distributional
dominance fails, welfare-based measures such as Atkinson Indexes, may provide
for a _complete_ ranking among alternative distributions, at the expense of
more stringent assumptions as to how to represent social welfare, leading to
weaker outcome robustness [Bellù and Liberati(2006)]. Furthermore, these
measures satisfy the additive decomposability axiom. Atkinson index has
support [0,1] and is defined as

$\displaystyle{\displaystyle
A(\varepsilon)={\begin{cases}1-{\frac{1}{\mu}}\left({\frac{1}{N}}\sum_{i}y_{i}^{1-\varepsilon}\right)^{1/(1-\varepsilon)}&{\mbox{for}}\
\varepsilon\neq 1\\\
1-{\frac{1}{\mu}}\left(\prod_{i}y_{i}\right)^{1/N}&{\mbox{for}}\
\varepsilon=1.\end{cases}}}$

The parameter $\varepsilon$ expresses the level of inequality aversion, as
$\varepsilon$ increases, the index becomes more sensitive to changes at the
lower end of the income distribution.

Another decomposable family of inequality measure is the Generalized Entropy
class. As opposed to the measures seen before, this class has the advantages
to be strongly transfer sensitive, meaning that it react to transfers
depending on donor and recipient income levels. It is based on the concept of
entropy which applied to income distributions has the meaning of deviations
from perfect equality:

$\displaystyle
GE(\alpha)={\begin{cases}{\frac{1}{N\alpha(\alpha-1)}}\sum_{i}\left[\left({\frac{y_{i}}{\mu}}\right)^{\alpha}-1\right]&\alpha\neq
0,1,\\\
{\frac{1}{N}}\sum_{i}{\frac{y_{{i}}}{\mu}}\ln{\frac{y_{{i}}}{\mu}}&\alpha\rightarrow
1,\\\ -{\frac{1}{N}}\sum_{i}\ln{\frac{y_{{i}}}{\mu}}&\alpha\rightarrow
0.\end{cases}}$

with $\alpha$ capturing the sensitivity of the index to a specific part of the
income distribution, a large $\alpha$ let the index be more sensitive to the
upper tail, vice versa a small $\alpha$ to the lower tail. $GE(0)$ is the Mean
Log Deviation, while $GE(1)$ is the more famous Theil index. Atkinson and
Generalized Entropy are two interrelated parametric families of measures, as a
transformation of the Atkinson Index is a member of the GE class:

$\displaystyle
A(\varepsilon)=1-[\varepsilon(\varepsilon-1)\cdot\text{GE}(1-\varepsilon)+1]^{1/(1-\varepsilon)}.$

In this paper, we will consider the estimation of both classes separately
since common parameter values used in one family does not correspond
deterministically to common parameter values used in the other one. Lastly, we
consider as inequality measure the Coefficient of Variation defined in
population as $CV=\sqrt{Var(y)}/\mu$. Its square has been used in some income
distribution analyses, including [OECD(2011)], but comparisons developed using
this measure seems to be very sensitive to top outliers [Atkinson(2015)].

##  3 Estimators Bias-Correction Proposal

The bias of the IM estimators in small samples drawn through a complex design
can be due to the structure of IM as non linear function of unbiased (or
sometimes even biased) estimators, leading them to be biased. The bias can be
either positive or negative, depending on the characteristics of the reference
variable distribution, except for the Mean Log Deviation ($GE(\alpha=0)$)
which has structurally negative bias; this aspect is made clearer in the
following. Among the measures with not predictable bias direction,
[Breunig(2001)] shows that $CV$ and $GE(2)$ bias direction is negatively
related with skewness. As regards the other measures, it depends on the
magnitude of the coefficient of variation related to the income variable and
to a monotonic transformation of it. This aspect could be analyzed in-depth by
imposing a distributional assumption on the income variable, however this is
out of scope.

######  Proposition 1.

As stated by [Breunig and Hutchinson(2008)], for a subset of the considered
measures, i.e. the ones pertaining to GE and Atkinson families, the
relationship between the expectation of the sample measure of inequality
$\hat{\theta}$ and its true population value $\theta$ is:

$\displaystyle\operatorname{\mathbb{E}}(\hat{\theta})=\theta+O_{p}\bigg{(}\frac{1}{n_{iid}}\bigg{)},$
(1)

with $n_{iid}$ denoting the sample size in the i.i.d. case.

###### Proof.

Let us consider a sample with i.i.d. elements $\\{y_{1},\dots,y_{n_{iid}}\\}$,
drawn from a population via simple random sampling, where $y_{i}$ is the
variable for the $i$-th unit with expected value $\mu$ and variance
$\sigma^{2}$. Let us consider also $\\{g(y_{1}),\dots,g(y_{n_{iid}})\\}$ with
$g(\cdot)$ defined a generic function $g:\mathbb{R}^{+}\rightarrow\mathbb{R}$
that changes for each measure, having expected value $\gamma$. Considering
equation (4) with $\hat{\mu}=\sum_{i=1}^{n_{iid}}y_{i}/n_{iid}$ and
$\hat{\gamma}=\sum_{i=1}^{n_{iid}}g(y_{i})/n_{iid}$, we can easily obtain
estimator moments as $\hat{\mu}\sim[\mu,\sigma^{2}/n_{iid}]$ and
$\hat{\gamma}\sim[\gamma,\phi^{2}/n_{iid}]$. Let us consider moreover that

$\displaystyle
Cov(\hat{\mu},\hat{\gamma})=\operatorname{\mathbb{E}}(\hat{\mu}\hat{\gamma})-\mu\gamma=\frac{1}{n_{iid}}(\operatorname{\mathbb{E}}(y\cdot
g(y))-\mu\gamma)=\frac{Cov(y,g(y))}{n_{iid}}.$ (2)

By expanding the inequality measure estimator $\hat{\theta}$ as a function of
$\hat{\mu}$ and $\hat{\gamma}$ via Taylor’s expansion around the population
values and considering its expected value:

$\displaystyle\operatorname{\mathbb{E}}(\hat{\theta})$
$\displaystyle=\theta+\frac{1}{2}f_{\gamma,\gamma}(\gamma,\mu)V(\hat{\gamma})+f_{\gamma,\mu}(\gamma,\mu)Cov(\hat{\gamma},\hat{\mu})+\frac{1}{2}f_{\mu,\mu}(\gamma,\mu)V(\hat{\mu})+O_{p}(n_{iid}^{-2})$
$\displaystyle=\theta+O_{p}(n_{iid}^{-1})+O_{p}(n_{iid}^{-1})+O_{p}(n_{iid}^{-1})+O_{p}(n_{iid}^{-2})$
$\displaystyle=\theta+O_{p}(n_{iid}^{-1}),$ (3)

where $f_{\gamma}=\frac{\partial f(\gamma,\mu)}{\partial\gamma}$ and
$f_{\gamma\mu}=\frac{\partial^{2}f(\gamma,\mu)}{\partial\gamma\partial\mu}$. ∎

Our bias-correction proposal constitutes as a generalization of the framework
of [Breunig and Hutchinson(2008)], developed for i.i.d. observations, to the
finite population and full design-based setting. At the same time, we extend
the proposal to a wider set of measures comprising Gini Index. We provide a
closed-form bias correction for complex designs which permits to avoid the use
of resampling techniques and it can be applied in a distribution-free setting
at once. This generalization has been developed considering Horwitz-Thompson
type design-based estimator, Ultimate Clusters technique and Influence
Function linearization technique for the estimation of design variances and
covariances.

We are interested in a variety of non linear functions of income values as IM
are. Let denote $\mathcal{S}$ a sample of size $n$, drawn using a complex
probability sampling design, with $p(s)=Pr(\mathcal{S}=s)$ the probability of
selecting the particular sample $s\subset\mathcal{U}$, thus $p(s)\geq 0$ and
$\sum_{s\in\mathcal{S}}p(s)=1$ . The inclusion probability of unit $i$ is
denoted with $\pi_{i}$ and defined such that $\pi_{i}=Pr(i\in\mathcal{S})$.

We consider the generic inequality measure written as a function of the mean
$\mu$ and the expected value of a generic function $g(\cdot)$ which we will
generically refer to as $\gamma=\operatorname{\mathbb{E}}(g(y))$. The
population value for the generic inequality measure is

$\displaystyle\theta=f(\mu,\gamma),$ (4)

with $f(\cdot)$ a twice-differentiable function. The related estimator in our
complex survey framework is $\hat{\theta}=f(\hat{\mu},\hat{\gamma})$ in which
Horvitz-Thompson complex survey estimators of the mean and the generic
component $\gamma$ are plugged in, i.e.

$\displaystyle\hat{\mu}=\frac{\sum_{i=1}^{n}w_{i}y_{i}}{\sum_{i=1}^{n}w_{i}}\quad\text{and}\quad\hat{\gamma}=\frac{\sum_{i=1}^{n}w_{i}g(y_{i},w_{i})}{\sum_{i=1}^{n}w_{i}}.$
(5)

where $w_{i}=1/\pi_{i}$ or a treated and calibrated version of it.
[Kakwani(1990)] uses a similar approach to express inequality indices in order
to derive their asymptotic standard error. By simply applying a second order
Taylor’s series expansion of the sample estimator around the population values
and evaluating its expected value, the bias is

$\displaystyle\operatorname{\mathbb{E}}(\hat{\theta}-\theta)=$ $\displaystyle
f_{\gamma}(\gamma,\mu)\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)+\frac{1}{2}f_{\gamma,\gamma}(\gamma,\mu)[V(\hat{\gamma})+\operatorname{\mathbb{E}}^{2}(\hat{\gamma}-\gamma)]+$
$\displaystyle+f_{\gamma,\mu}(\gamma,\mu)[Cov(\hat{\gamma},\hat{\mu})-\mu\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)]+\frac{1}{2}f_{\mu,\mu}(\gamma,\mu)V(\hat{\mu})+O_{p}(n^{-2}).$
(6)

Notice that $\hat{\mu}$ is unbiased and (6) appears as a generalization of
equation (3) to the finite population case and by considering $\hat{\gamma}$
estimator may be biased.

An alternative method, dealing with IM bias in small sample, is the
small-$\sigma$ approximation described in [Ullah(2004)]. However this
framework requires high order moments and cross-moments estimation when facing
non i.i.d. assumptions, which results quite challenging in cases of multi-
stage surveys and distribution-free setting.

So we propose approximate bias corrections for Gini Index, General Entropy,
Atkinson Indexes and Coefficient of Variation based on (6). In the following
the corrections for the whole set of inequality measures is made explicit.

###  3.1 Approximate Bias-Corrected Inequality Measures

  * •

Coefficient of Variation

$\displaystyle
CV=\frac{\sqrt{\gamma-\mu^{2}}}{\mu}\quad\text{with}\quad\gamma=\operatorname{\mathbb{E}}(Y^{2}).$
(7)

Considering that the complex survey estimator is

$\displaystyle\hat{CV}=\sqrt{\frac{n^{\prime}}{n^{\prime}-1}}\bigg{(}\frac{1}{\hat{N}}\sum_{i\in
S}w_{i}\bigg{(}\frac{y_{i}}{\hat{\mu}}\bigg{)}^{2}-1\bigg{)}^{1/2}=\sqrt{\frac{n^{\prime}}{n^{\prime}-1}}\frac{\sqrt{\hat{\gamma}-\hat{\mu}^{2}}}{\hat{\mu}},$
(8)

with $\hat{\gamma}=\sum_{i\in S}w_{i}y_{i}^{2}/\sum_{i\in S}w_{i}$ and
$\hat{\mu}$ as indicated in (5), $n^{\prime}$ the number of observations with
weights different than zero and $\sqrt{n^{\prime}/(n^{\prime}-1)}$ the
standard bias-correction adjustment for the weighted variance. Consider that
$\hat{\gamma}$ is unbiased.

The approximate bias in CV is:

$\displaystyle\operatorname{\mathbb{E}}[\hat{CV}]-CV\cong$
$\displaystyle\sqrt{\frac{n^{\prime}}{n^{\prime}-1}}\bigg{[}\bigg{(}-\frac{1}{8\mu}(\gamma-\mu^{2})^{-3/2}\bigg{)}V(\hat{\gamma})$
$\displaystyle-\frac{1}{2}\bigg{(}\mu^{-2}(\gamma-\mu^{2})^{-1/2}-(\gamma-\mu^{2})^{-3/2}\bigg{)}Cov(\hat{\gamma},\hat{\mu})$
(9)
$\displaystyle+\frac{1}{2}\bigg{(}2\mu^{-3}(\gamma-\mu^{2})^{1/2}+\mu^{-1}(\gamma-\mu^{2})^{-1/2}-\mu(\gamma-\mu^{2})^{-3/2}\bigg{)}V(\hat{\mu})\bigg{]}.$

  * •

Gini Index

By using [Sen(1973)] alternative formulation for Gini index:

$\displaystyle
G=\frac{2\gamma}{\mu}-1\quad\text{where}\quad\gamma=\operatorname{\mathbb{E}}(Y\cdot
F(Y)),$ (10)

considering $F(\cdot)$ the cumulative distribution function of the variable of
interest and the complex survey estimator proposed by [Langel and
Tillé(2013)]:

$\displaystyle\hat{G}$ $\displaystyle=\frac{2\sum_{i\in
S}w_{i}y_{i}(\hat{N_{i}}-w_{i}/2)}{\hat{N}^{2}\mu}-1=\frac{2\hat{\gamma}}{\hat{\mu}}-1,$
(11)

where $\hat{N}_{i}=\sum_{k\in S}w_{k}I(n_{k}\leq n_{i})$ with $n_{i}$ denoting
the rank of the unit $i$ in sample $\mathcal{S}$ and the biased estimator
defined as

$\displaystyle\hat{\gamma}=\frac{\sum_{i\in
S}w_{i}y_{i}(\hat{N_{i}}-\frac{w_{i}}{2})}{\hat{N}^{2}}.$ (12)

Thus, the approximate bias in small sample is

$\displaystyle\operatorname{\mathbb{E}}(\hat{G})-G$
$\displaystyle\cong\frac{2}{\mu}\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)+\frac{2\gamma}{\mu^{3}}V(\hat{\mu})-\frac{2}{\mu^{2}}[Cov(\hat{\mu},\hat{\gamma})-\mu\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)]$
(13)
$\displaystyle=\frac{4}{\mu}\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)+\frac{2\gamma}{\mu^{3}}V(\hat{\mu})-\frac{2}{\mu^{2}}Cov(\hat{\mu},\hat{\gamma}),$

The derivation of the approximate bias related to the weighted estimator
$\hat{\gamma}$ is not trivial. As explained in the context of Gini index
variance estimation by [Langel and Tillé(2013)], its numerator is not composed
of two simple sums. Indeed the quantity $\hat{N}_{k}$, an estimator of the
rank of unit $k$, is random since its value depends on the selected sample. An
heuristic solution is to consider the approximate bias of the corresponding
i.i.d. estimator,
$\operatorname{\mathbb{E}}(\hat{\gamma}-\gamma)=-\frac{1}{n_{iid}}\bigg{(}\gamma-\frac{\mu}{2}\bigg{)}$
as derived by [Davidson(2009)], with $n_{iid}$ the corresponding i.i.d. sample
size, so that:

$\displaystyle\operatorname{\mathbb{E}}(\hat{G})-G=\frac{-2G}{n_{iid}}+\frac{2\gamma}{\mu^{3}}V(\hat{\mu})-\frac{2}{\mu^{2}}Cov(\hat{\mu},\hat{\gamma}).$
(14)

This correction is in line with [Davidson(2009)] and [Fabrizi and
Trivisano(2016)] proposals, whereas these are based on a first order Taylor’s
expansions and thus limiting to the first term of the right hand side equation
(13), ours extends it to a second order expansion. This translates into the
fact that, while [Jasso(1979)] [Deltas(2003)], [Davidson(2009)] proposals
identify the adjusted Gini in i.i.d. contexts as
$n_{iid}(n_{iid}-1)^{-1}\hat{G}$ , our correction, with a further order of
approximation, reconsiders the shape of the adjusted estimator as

$\displaystyle\hat{G}_{adj}=\frac{n_{iid}}{n_{iid}-2}(\hat{G}-a),$ (15)

with $a$ equals to the sum of second and third term of (14).

  * •

Mean Log Deviation- GE($\alpha=0$)

$\displaystyle
GE(0)=\log(\mu)-\gamma\quad\text{with}\quad\gamma=\operatorname{\mathbb{E}}(\log
Y).$ (16)

Let us note that the survey estimator [Biewen and Jenkins(2006)] is

$\displaystyle\hat{GE}(0)=\frac{1}{\hat{N}}\sum_{i\in
S}w_{i}\log\frac{\hat{\mu}}{y_{i}}=\log(\hat{\mu})-\hat{\gamma},$ (17)

with $g(y_{i})=\log y_{i}\quad\forall i\in S$ and so $\hat{\gamma}=\sum_{i\in
S}w_{i}\log y_{i}/\sum_{i\in S}w_{i}$. The approximate bias is:

$\displaystyle\operatorname{\mathbb{E}}[\hat{GE}(0)]-GE(0)=\operatorname{\mathbb{E}}(\log\hat{\mu})-\log\mu\cong-\frac{1}{2\mu^{2}}V(\hat{\mu}),$
(18)

consistently with [Ferrante and Pacei(2019)]; note that it is structurally
negative.

  * •

Theil Index- GE($\alpha=1$)

$\displaystyle
GE(1)=\frac{\gamma}{\mu}-\log(\mu)\quad\text{with}\quad\gamma=\operatorname{\mathbb{E}}[Y(\log
Y)].$ (19)

Let us note that the survey estimator [Biewen and Jenkins(2006)] is

$\displaystyle\hat{GE}(1)=\frac{1}{\hat{N}}\sum_{i\in
S}w_{i}\frac{y_{i}}{\hat{\mu}}\log\frac{y_{i}}{\hat{\mu}}=\frac{\hat{\gamma}}{\hat{\mu}}-\log(\hat{\mu}),$
(20)

with $g(y_{i})=y_{i}\log y_{i}\quad\forall i\in S$ and so
$\hat{\gamma}=\sum_{i\in S}w_{i}y_{i}\log y_{i}/\sum_{i\in S}w_{i}$. The
approximate bias is:

$\displaystyle\operatorname{\mathbb{E}}[\hat{GE}(1)]-GE(1)\cong-\frac{Cov(\hat{\mu},\hat{\gamma})}{\mu^{2}}+\bigg{(}\frac{\gamma}{\mu^{3}}+\frac{1}{2\mu^{2}}\bigg{)}V(\hat{\mu}).$
(21)

  * •

GE($\alpha\neq 0,1$)

$\displaystyle GE(\alpha\neq
0,1)=\frac{1}{\alpha(\alpha-1)}\bigg{(}\frac{\gamma_{\alpha}}{\mu^{\alpha}}-1\bigg{)}\quad\text{with}\quad\gamma_{\alpha}=\operatorname{\mathbb{E}}(Y^{\alpha}).$
(22)

Considering that the survey estimator is

$\displaystyle\hat{GE}(\alpha\neq
0,1)=\frac{1}{\alpha(\alpha-1)}\frac{n^{\prime}}{n^{\prime}-1}\bigg{(}\frac{\hat{\gamma}_{\alpha}}{\hat{\mu}^{\alpha}}-1\bigg{)},$
(23)

[Biewen and Jenkins(2006)] with $\hat{\gamma}=\sum_{i\in
S}w_{i}y_{i}^{\alpha}/\sum_{i\in S}w_{i}$ and $\hat{\mu}$ as indicated in (5).
Since $\hat{\gamma}$ is unbiased as proved above, the approximate bias is

$\displaystyle\operatorname{\mathbb{E}}[\hat{GE}(\alpha)]-GE(\alpha)\cong\frac{n^{\prime}}{n^{\prime}-1}\bigg{[}-\frac{1}{(\alpha-1)\mu^{\alpha+1}}Cov(\hat{\gamma},\hat{\mu})+\frac{\alpha+1}{2(\alpha-1)}\frac{\gamma_{\alpha}}{\mu^{\alpha+2}}V(\hat{\mu})\bigg{]}.$
(24)

  * •

Atkinson Index with $\epsilon\neq 1$

$\displaystyle
A(\epsilon)=1-\frac{\gamma^{\frac{1}{1-\epsilon}}}{\mu}\quad\text{with}\quad\gamma=\operatorname{\mathbb{E}}(Y^{1-\varepsilon}).$
(25)

Let us note that the survey estimator [Biewen and Jenkins(2006)] is

$\displaystyle\hat{A}(\varepsilon)=1-\frac{1}{\hat{\mu}}\bigg{(}\frac{1}{\hat{N}}\sum_{i\in
S}w_{i}(y_{i})^{1-\varepsilon}\bigg{)}^{1/(1-\varepsilon)}=1-\frac{\hat{\gamma}^{\frac{1}{1-\varepsilon}}}{\hat{\mu}},$
(26)

with $g(y_{i})=y_{i}^{1-\varepsilon}\quad\forall i\in S$ and so
$\hat{\gamma}=\sum_{i\in S}w_{i}y_{i}^{1-\varepsilon}/\sum_{i\in S}w_{i}$. The
approximate bias is:

$\displaystyle\operatorname{\mathbb{E}}[\hat{A}(\varepsilon)]-A(\varepsilon)\cong-\frac{\varepsilon}{2(1-\varepsilon)^{2}}\mu^{-1}\gamma^{\frac{2\varepsilon-1}{1-\varepsilon}}V(\hat{\gamma})+\frac{1}{1-\varepsilon}\gamma^{\frac{\varepsilon}{1-\varepsilon}}\mu^{-2}Cov(\hat{\gamma},\hat{\mu})-\mu^{-3}\gamma^{\frac{1}{1-\varepsilon}}V(\hat{\mu}).$
(27)

  * •

Atkinson Index with $\varepsilon=1$

$\displaystyle
A(1)=1-\frac{e^{\gamma}}{\mu}\quad\text{with}\quad\gamma=\operatorname{\mathbb{E}}(\log
Y).$ (28)

Let us note that the survey estimator [Biewen and Jenkins(2006)] is

$\displaystyle\hat{A}(1)=1-\frac{1}{\hat{\mu}}\prod_{i\in
S}y_{i}^{w_{i}/\hat{N}}=1-\frac{1}{\hat{\mu}}\exp\bigg{\\{}\frac{\sum_{i\in
S}w_{i}\log y_{i}}{\hat{N}}\bigg{\\}}=1-\frac{e^{\hat{\gamma}}}{\hat{\mu}},$
(29)

with $g(y_{i})=\log y_{i}\quad\forall i\in S$ and so $\hat{\gamma}=\sum_{i\in
S}w_{i}\log y_{k}/\sum_{k\in S}w_{k}$ as in the case of the Mean Log
Deviation. The approximate bias is:

$\displaystyle\operatorname{\mathbb{E}}[\hat{A}(1)]-A(1)\cong-\frac{1}{2\mu}e^{\gamma}V(\hat{\gamma})+\frac{e^{\gamma}}{\mu^{2}}Cov(\hat{\gamma},\hat{\mu})-\frac{e^{\gamma}}{\mu^{3}}V(\hat{\mu}).$
(30)

Note that the approximate bias expressions for complex survey estimators in
some specific cases coincide with the ones for i.i.d. estimators made explicit
by [Breunig and Hutchinson(2008)], see formulas (18), (21), (27), (30), due to
the invariance properties related to expected values of sum of dependent
variables.

##  4 Bias Estimation

With the aim of estimating bias of estimators proposed in Section 3, we
replace $\mu$ and $\gamma$ with $\hat{\mu}$ and $\hat{\gamma}$ in formulas
(9), (14), (18), (21), (24), (27) and (30). Moreover, we have to consider that
those expressions (except (18)) depend on variances and covariances involving
$\hat{\gamma}$, generally non linear. Then a linearization of $\hat{\gamma}$
has been provided in order to make it tractable (section 4.1). Subsequently,
$V(\hat{\mu})$, $V(\hat{\gamma})$ and $Cov(\hat{\mu},\hat{\gamma})$ has been
estimated (section 4.2).

###  4.1 Linearization of Non-Linear Design Estimators

The linearization technique follows the intuition of approximating a nonlinear
statistic with a linear function of the observations. In doing so, the linear
approximation can be used to measure the precision and uncertainty associated
to the nonlinear statistic, using well-known linear estimators variances and
covariances. We apply the generalized linearization method [Deville(1999)],
[Demnati and Rao(2004)], and [Osier(2009)]. This method allows to encompass
more non linear statistics than the Taylor method, and at the same time it
does not involve more calculations, therefore being in general more flexible
[Osier(2009)]. Moreover, this methods works better when dealing with small
samples.

In particular, this procedure as stated by [Antal et al.(2011)Antal, Langel,
and Tillé], reconciles the two approaches introduced by [Deville(1999)] and
[Demnati and Rao(2004)], both relying on the concept of influence function, a
framework borrowed from robust statistics [Hampel(1974)]. The same method has
been used directly by [Graf and Tillé(2014)] to estimate the variance of some
inequality measures estimators via linearization. Following the theoretical
framework of [Antal et al.(2011)Antal, Langel, and Tillé], we could say that a
population parameter of interest $\theta$ in a population $\mathcal{U}$
undergoes the influence of a unit $k$, which depends on an infinitesimal
variation in the importance assigned to the unit. Let us express the parameter
as a functional $\theta=T(M)$ based on a measure $M(\cdot)$ such that

$\displaystyle\begin{cases}M(y)=1\quad y=y_{k}\quad\forall k\in\mathcal{U}\\\
M(y)=0\quad\text{otherwise},\end{cases}$

The general measure M turns to a discrete measure, leading T into a discrete
functional. The influence function of T, that in this case constitutes as the
linearized variable $z_{k}$, is defined as the functional derivative

$\displaystyle I[T(M)]_{k}=z_{k}=\lim_{t\rightarrow
0}\frac{T(M+t\delta_{k})-T(M)}{t}\quad\forall k\in\mathcal{U},$

where $\delta_{k}$ is a Dirac measure for unit $k$. However [Deville(1999)]
defines a linearized variable $\hat{z}_{k}$, or empirical influence function,
which takes into account that data comes from a limited sample $s$, replacing
the unknowns with the corresponding estimated sample quantities. [Antal et
al.(2011)Antal, Langel, and Tillé] present this approach with measure M
formalized taking into account sampling weights and $s$ such that

$\displaystyle\begin{cases}M(y)=\tilde{w}_{k}\quad y=y_{k}\quad\forall k\in
s\\\ M(y)=0\quad\text{otherwise},\end{cases}$

with $\tilde{w}_{k}$ generic weight associated to the observation $y_{k}$ that
in our case corresponds to $\tilde{w}_{k}=w_{k}/\sum_{i\in S}w_{i}$. This
alternative approach completely overcomes the fact that the starting point of
Deville’s approach is the population parameter rather than the sample
estimator. Since the functional is expressed as an explicit function of the
variables, i.e. the weights assigned by the measure M to the observations, the
linearized variable is merely a function of the partial derivatives with
respect to the weights:

$\displaystyle I[T(\hat{M})]_{k}=\hat{z}_{k}=\frac{\partial
T(\hat{M})}{\partial\tilde{w}_{k}}=\frac{\partial}{\partial\tilde{w}_{k}}\sum_{i\in
S}\tilde{w}_{i}g(y_{i})=g(y_{k}).$ (31)

Linearized variables for each inequality measure can be directly derived from
(31) by substituting the measure-specific function $g(\cdot)$ explicited in
Section 3.

###  4.2 Estimation of $V(\hat{\mu})$, $V(\hat{\gamma})$ and
$Cov(\hat{\mu},\hat{\gamma})$

As regards the estimation of the design variances and covariances of linear
and linearized estimators, we consider a complex survey design involving
stratification and multi-stage selection, with both Self-Representing (SR)
-included at the first stage with probability one - and Non Self-Representing
(NSR) strata. This design is consistent with the most common income survey
designs and, in general, with household surveys, whereas firm surveys are
based on a single stage stratified sample design, also a special case of the
two stage selection.

First of all we define the shape of the Horvitz Thompson variance estimator in
case of linear (or linearized) estimators, such as $\hat{\mu}=\sum_{i\in
s}\tilde{w}_{i}y_{i}$, when $\tilde{w}_{i}=1/\pi_{i}$, as:

$\displaystyle V(\hat{\mu})=\sum_{i\in
s}\frac{y_{i}^{2}}{\pi_{i}^{2}}(1-\pi_{i})+2\sum_{i\in s}\sum_{k\in s,i\neq
k}\frac{y_{i}}{\pi_{i}}\frac{y_{k}}{\pi_{k}}\frac{\pi_{ik}-\pi_{i}\pi_{k}}{\pi_{ik}},$
(32)

with $\pi_{ik}$, $\forall i,k\in\mathcal{U},i\neq k$ second order inclusion
probabilities i.e. the probability that the sample includes both $i$-th and
$k$-th units, $\pi_{ik}=\sum_{s\in\mathcal{S}_{i}\cap\mathcal{S}_{k}}p(s)$.
However generally (a) $\tilde{w}_{i}\neq 1/\pi_{i}$ and (b) $\pi_{ik}$,
$\forall i,k\in\mathcal{U},i\neq k$ are difficult to calculate under complex
sampling designs.

Therefore, the variance estimator to be considered constitutes as an
approximation relying on simplified assumptions. We assume Primary Sampling
Units (PSU) are sampled with replacement, and secondly multi-stage sampling is
reduced to a single stage process using the Ultimate Clusters technique
[Kalton(1979)]. Moreover we take into account the hybrid nature of the
probability scheme, blending a stratified design variance estimator for the SR
strata with a finite population correction factor and a typical Ultimate
Cluster approximation-based variance estimator for the multi-stage scheme in
NSR strata, widely used in official statistics. See [Osier et al.(2013)Osier,
Berger, and Goedeme] for Eurostat procedures. Therefore, considering without
loss of generality a two-stage scheme, let
$\hat{\mu}=\sum_{h}\sum_{d}\sum_{i}\tilde{w}_{hdi}y_{hdi}$ with $h$ stratum
indicator, $d$ Primary Sampling Unit (PSU) indicator and $i$ Secondary
Sampling Unit indicator (SSU), be a linear estimator for $\mu$, its standard
error estimate is as follows:

$\displaystyle\hat{V}(\hat{\mu})$
$\displaystyle=\sum_{h=1}^{H_{SR}}V(\hat{\mu}_{h})+\sum_{h=1}^{H_{NSR}}V(\hat{\mu}_{h})$
$\displaystyle=\sum_{h=1}^{H_{SR}}M_{h}^{2}(1-f_{h})\frac{s_{h}^{2}}{m_{h}}+\sum_{h=1}^{H_{NSR}}n_{h}s_{\hat{\mu}_{h}}^{2}$
(33)
$\displaystyle=\sum_{h=1}^{H_{SR}}M_{h}\frac{M_{h}-m_{h}}{m_{h}(m_{h}-1)}\sum_{i=1}^{m_{h}}(y_{hi}-\bar{y}_{h})^{2}+\sum_{h=1}^{H_{NSR}}\frac{n_{h}}{n_{h}-1}\sum_{d=1}^{n_{h}}(\hat{\mu}_{hd}-\bar{\mu}_{h})^{2},$

with $H_{SR}$ self-representative and $H_{NSR}$ non self-representative
strata, $M_{h}$ number of resident households in strata $h$, $m_{h}$ number of
sample households in strata $h$, $f_{h}=m_{h}/M_{h}$ finite population
correction factor, $n_{h}$ number of PSUs in strata $h$. Consider
$\bar{y}_{h}=\sum_{i=1}^{m_{h}}y_{hi}/m_{h}$,
$\hat{\mu}_{hd}=\sum_{i=1}^{m_{d}}\tilde{w}_{hdi}y_{hdi}$ with $i$ household
label of stratum $h$ and PSU $d$, with a total of $m_{i}$ per PSU and
$\bar{\mu}_{h}=\sum_{d=1}^{n_{h}}\hat{\mu}_{hd}/n_{h}$, lastly $n_{h}$ is the
number of PSU in stratum $h$. If however $n_{h}=1$ for some strata, the
estimator (33) cannot be used. A solution is to collapse strata to create
“pseudo-strata” so that each pseudo-stratum has at least two PSUs. Common
practice is to collapse a stratum with another one that is similar w.r.t. the
target variables of the survey [Rust and Kalton(1987)].

Secondly, an estimator of $V(\hat{\gamma})$ can be obtained after the
linearization of $\hat{\gamma}$, leading to $V(\hat{\gamma})\cong V(\sum_{i\in
S}\tilde{w}_{i}\hat{z}_{i})$, and by adopting the same strategy used for
$V(\hat{\mu})$ in (33).

Thirdly, as regards the estimation of the design-based estimators covariance,
$Cov(\hat{\gamma},\hat{\mu})$ let us consider that

$\displaystyle
Cov(\hat{\gamma},\hat{\mu})=\frac{1}{2}\bigg{[}V(\hat{\gamma}+\hat{\mu})-V(\hat{\gamma})-V(\hat{\mu})\bigg{]}.$
(34)

Thus, a possible estimator for the design covariance would be simply obtained
by plugging in the variance estimates previously mentioned, while the first
term of the left-hand side of the formula could be estimated by simply
considering $\hat{\gamma}+\hat{\mu}=\sum_{i\in
S}\tilde{w}_{i}(\hat{z}_{i}+y_{i})$ using (33). In fact, the estimation
procedure involves substituting the linearized variables $\hat{z}_{k}$,
different for each measure, in (33) and (34) in order to estimate
$V(\hat{\gamma})$ and $Cov(\hat{\mu},\hat{\gamma})$.

Lastly as regards Gini index, the $\hat{\gamma}$ bias (B) has already been
adapted to the complex survey case by [Fabrizi and Trivisano(2016)], with an
heuristic solution as follows

$\displaystyle\hat{B}(\hat{\gamma})=-\frac{\sum_{q\in
S}(\sum_{i=1}^{r_{q}}w_{qi})^{2}}{\hat{N}^{2}}\bigg{(}\hat{\gamma}-\frac{\hat{\mu}}{2}\bigg{)},$
(35)

with $\sum_{i=1}^{r_{q}}w_{qi}$ the sum of the weights associated with the
$r_{q}$ individuals living in household $q$. However, in our case we opt for
$-\frac{1}{n^{\prime}}\bigg{(}\hat{\gamma}-\frac{\hat{\mu}}{2}\bigg{)}$, since
involving survey weights could dramatically bias the estimator of the
effective sample size quantity. Thus our bias estimator is

$\displaystyle\hat{B}(\hat{G})$
$\displaystyle\cong-\frac{2G}{n^{\prime}}+\frac{2\hat{\gamma}}{\hat{\mu}^{3}}V(\hat{\mu})-\frac{2}{\hat{\mu}^{2}}Cov(\hat{\mu},\hat{\gamma}).$
(36)

##  5 Design-Based Simulation on Bias Correction

A first design-based simulation study has been carried out in order to
evaluate the bias correction framework proposed in the previous section. In
this simulation, the cross-section Italian EU-SILC sample (2017 wave) has been
assumed as synthetic population and the 21 NUTS-2 regions have been considered
as target domains. The study is based on real data, rather than use data
generated under some specific distribution models, in order to check whether
this specific framework can work with close-to-reality income data.

However, those data have been treated in order to circumvent non-robustness
problems. The issue of robust estimation of economic indicators based on a
semi-parametric Pareto upper tail model is well-established in literature see
[Brzezinski(2016)] for a review and [Alfons et al.(2013)Alfons, Templ, and
Filzmoser] for a specification suitable for survey data. On the contrary, the
issue of robust treatment of outlier in the lower tail of income distribution
appears less established, see [Van Kerm(2007)], [Masseran et
al.(2019)Masseran, Safari, and Ibrahim]. As regards the upper tail we operated
a semi-parametric Pareto-tail modeling procedure using the Probability
Integral Transform Statistic Estimator (PITSE) proposed by [Finkelstein et
al.(2006)Finkelstein, Tucker, and Alan Veeh], which blends very good
performances in small samples and a fast computational implementation, as
suggested by [Brzezinski(2016)]. As regards the lower tail extreme value
treatment, we used an inverse Pareto modification of PITSE estimator,
suggested by [Masseran et al.(2019)Masseran, Safari, and Ibrahim]. In our
simulations the treatment has been done at a regional level to the original
EU-SILC sample and the detection of outliers has been carried out following
[Safari and Ibrahim(2018)] by using a Generalized Boxplot outlier detection
procedure proposed by [Bruffaerts et al.(2014)Bruffaerts, Verardi, and
Vermandele]. In order to deal with skewed and heavy-tailed distribution,
[Bruffaerts et al.(2014)Bruffaerts, Verardi, and Vermandele] suggested
modifying the whiskers of the boxplot by a simple rank-preserving
transformation that allows for a so-called fitting for $Tukey$ $g-$ and $h-$
distribution. Based on the simulations performed by [Safari and
Ibrahim(2018)], the generalized boxplot gives a higher power value compared to
the adjusted and standard boxplot. Moreover, we compare the results obtained
after the treatment with those obtained by retained the outliers in the
synthetic population, in order to isolate the effect of outliers on the
performance of the corrections proposed (Table 1).

From the assumed population, we repeatedly select 1000 two-stage stratified
samples, mimicking the sampling strategy adopted in the survey itself: in the
first stage, Self-Representing (SR) strata are always included in the sample,
while a stratified sample of PSU in Non Self-Representing (NSR) strata is
selected; in the second stage, a systematic sample of households is selected
from each municipality included in the first stage. We repeated the drawing
for two different overall sampling rates, 1,5% and 3% respectively. Results
are set out in Table 1 with the Average Relative Bias (ARB) and the Average
Absolute Relative Error (AARE) in percentage, calculated for the 1000 samples
and the $21$ regions defined as:

$\displaystyle ARB$
$\displaystyle=\frac{1}{21}\sum_{r=1}^{21}\frac{1}{1000}\sum_{m=1}^{1000}\bigg{(}\frac{\hat{\theta}_{m,r}}{\theta_{r}}-1\bigg{)},$
$\displaystyle AARE$
$\displaystyle=\frac{1}{21}\sum_{r=1}^{21}\frac{1}{1000}\bigg{(}\sum_{m=1}^{1000}\bigg{|}\frac{\hat{\theta}_{m,r}}{\theta_{r}}-1\bigg{|}\bigg{)},$

where $\theta_{r}$ is the population value for region $r$ and
$\hat{\theta}_{m,r}$ is its estimate using sample $m$. In our simulation
setting the regional sample size ranges from 6 to 96 individuals (from 6 to 32
households) for the 1,5% sampling rate, and from 11 to 196 individuals (10 to
74 households) for the 3% sampling rate.

As clear from Table 1, bias reported on data without Extreme Value treatment
can be dramatically high due to the non-robustness properties of some measure
to extreme values, as the case of A($\varepsilon=2$), extremely sensible to
low income values (under $100$ euro per year) which is -48% biased on average
in case of 1,5% sample rate. Also GE with $\alpha=1,2$ values are highly
sensible to high income values being -18% and -23% biased. Moreover, the
negative correlation between sample size and bias is less marked in case of
non-treated data, since bias depends more on whether a sample is affected by
extreme values or not. On the contrary, the bias-correction seems to not
change in magnitude depending on the sample size and on the presence of
extreme values, showing robustness properties.

|  |  CV | GE(0) | GE(1) | GE(2) | A(0.5) | A(1) | A(2) | G
---|---|---|---|---|---|---|---|---|---
without Extreme Value treatment |  |  |  |  |  |
1.5% |  |  |  |  |  |  |  |  |
| ARB | -18.2 | -12.7 | -17.5 | -23.3 | -15.3 | -15.6 | -48.0 | -14.4
$\hat{\theta}\quad$ | AARE | 30.0 | 52.9 | 46.4 | 53.5 | 45.9 | 47.4 | 56.8 | 25.6
| ARB | -12.1 | -3.9 | -8.7 | -15.0 | -6.3 | -5.9 | -41.6 | 0.4
$\hat{\theta}_{corr}$ | AARE | 29.3 | 54.4 | 48.0 | 55.7 | 47.5 | 49.4 | 54.6 | 35.5
3.0% |  |  |  |  |  |  |  |  |
| ARB | -12.7 | -6.8 | -10.5 | -15.8 | -8.7 | -8.4 | -38.1 | -7.3
$\hat{\theta}\quad$ | AARE | 24.5 | 39.4 | 36.0 | 46.2 | 34.3 | 35.6 | 49.0 | 17.7
| ARB | -7.8 | -1.2 | -3.9 | -8.0 | -2.5 | -2.0 | -32.4 | 0.2
| ARB ($n\geq 20$) | -8.3 | -1.2 | -3.6 | -8.7 | -2.3 | -1.7 | -30.0 | -1.1
$\hat{\theta}_{corr}$ | AARE | 24.8 | 40.4 | 37.9 | 49.4 | 35.7 | 37.0 | 48.2 | 20.8
with Extreme Value treatment |  |  |  |  |  |
1.5% |  |  |  |  |  |  |  |  |
| ARB | -11.9 | -13.1 | -15.3 | -17.0 | -14.2 | -14.3 | -18.3 | -14.2
$\hat{\theta}\quad$ | AARE | 25.8 | 44.0 | 42.6 | 47.6 | 41.6 | 40.6 | 38.2 | 24.3
| ARB | -5.6 | -3.5 | -6.3 | -8.7 | -4.9 | -4.7 | -8.6 | 0.6
$\hat{\theta}_{corr}$ | AARE | 25.9 | 46.2 | 44.6 | 49.9 | 43.6 | 42.5 | 39.3 | 34.1
3.0% |  |  |  |  |  |  |  |  |
| ARB | -7.4 | -6.4 | -8.3 | -10.3 | -7.3 | -7.1 | -9.4 | -7.1
$\hat{\theta}\quad$ | AARE | 19.8 | 31.6 | 31.7 | 37.8 | 30.2 | 29.1 | 27.2 | 16.5
| ARB | -2.8 | -0.6 | -2.2 | -3.5 | -1.4 | -1.1 | -2.8 | 0.3
| ARB ($n\geq 20$) | -1.4 | -0.3 | -1.2 | -1.5 | -0.7 | -0.5 | -1.4 | -0.9
$\hat{\theta}_{corr}$ | AARE | 20.4 | 33.0 | 33.3 | 40.2 | 31.7 | 30.4 | 28.4 | 19.5
Table 1: Percentage ARB and AARE for the 21 target domains. Figure 1: Relative
Bias of non-corrected measures (grey line), and of corrected measures (blue
line) in 3% samples after extreme value treatment.

Concerning extreme value treated data results, Figure 1 clearly illustrates
the negative correlation between sample size and average relative bias in the
21 italian regions for both the design-based estimator $\hat{\theta}$ and the
bias corrected estimator $\hat{\theta}_{corr}$. The reduction of the bias
provided by the correction is noticeable for all measures, leading to slightly
biased (-1/-2%) or even approximately unbiased estimates when $n\geq 20$
depending on the measure. Notice that the bias correction works well for those
measures not particularly sensitive to extreme observations such as Gini
index, GE$(0)$, Atk$(0.5)$ and Atk$(1)$. In case of CV and GE$(2)$, the
correction provides good results, but it seems however to not capture totally
all the bias components. This confirms the results of [Breunig(2001)],
suggesting that the coefficient of variation squared and GE(2) bias depends on
the coefficient of skewness of the income distribution (remember that
GE$(2)=$CV${}^{2}/2$), not considered in our bias correction. Actually, a
reliable estimation of that quantity, while being straightforward in the
i.i.d. case, appears cumbersome in case of weighted data being basically
defined on a discrete grid of values. This leads to the non applicability of
the bias formula derived by [Breunig(2001)] in our case. Furthermore, the
bias-correction induces a slight but negligible error increase, except for
Gini index case that presents a relevant increase. This is due to the shape of
the unbiased estimators, as described by (15), where a sum of estimators is
multiplied by a factor $n_{iid}/(n_{iid}-2)$, substituted by
$n^{\prime}/(n^{\prime}-2)$ in the estimation, that inherently inflates the
variance by its square.

This results may constitute as a valuable reference guideline when measuring
inequality in small samples. When extreme values can be considered as a
consequence of data contamination, extreme value treatment and bias-correction
provides approximately unbiased or slightly bias estimates for a large class
of measures. On the other hand, when extreme values constitutes as
representative observations, it becomes necessary to restrict the attention to
the most robust measures such as GE with $\alpha=0$, Atkinson index with
$\varepsilon=1$ and Gini Index. Another important aspect to point out is that,
in certain countries, the EU-SILC is based on registers that better capture
top incomes, thus, a cross-country comparison of income inequality by effects
on a tail-sensitive measure must be another reason for caution
[Atkinson(2015)].

##  6 Bootstrap Variance of the Bias-Corrected Estimators

Concerning inequality design-based estimators, their variance estimation may
be easily carried out via linearization as seen in Section 3. Linearized
variables for each measure may be easily derived from (31) consistently with
[Langel and Tillé(2013)] results for Gini Index and [Biewen and Jenkins(2006)]
results for Generalized Entropy and Atkinson Indexes. On the other hand, the
estimation of bias-corrected estimators variance adds a new level of
complexity since the estimator formula is no longer the classical one. In fact
it comprises a bias-correction component that appears cumbersome to estimate
via linearization since it is inherently a results of several linearizations,
leading to an over-simplified and unreliable linearization loop.

We approach the variance estimation issue by taking into account strategies
based on resampling methods, that have been largely used in literature, since
our main aim is to provide a general ”turn-key” solution for variance
estimation. Specifically, we opt for a proper design-aware bootstrap procedure
as developed before by [Fabrizi et al.(2011)Fabrizi, Ferrante, Pacei, and
Trivisano] and [Fabrizi et al.(2020)Fabrizi, Ferrante, and Trivisano] for EU-
SILC data. A comprehensive review for bootstrap method in case of survey data
can be found in [Lahiri(2003)] and an interesting comparison between variance
estimation techniques for poverty and inequality measures in complex surveys
has been carried out by [De Santis et al.(2020)De Santis, Barabesi, and
Betti]. Generally, bootstrap algorithms for complex samples rely on the
assumption that the number of strata is large and on resampling of primary
units within strata: few primary units are sampled from each stratum so that
the sampling fraction at the first stage is negligible [Rao et al.(1999)Rao,
Chaudhuri, Eltinge, Fay, Ghosh, Ghosh, Lahiri, and Pfeffermann]. This
assumption is however not satisfied in small samples. An alternative solution
may be split up primary units into multiple parts to resample, to extremes
into secondary units. Thus, our strategy selects units with a stratified
single-stage design with replacement from the population of households
considering geographical macro-strata. After the drawing, a calibration
procedure has been put in place for each bootstrap sample, in order to adjust
weights to known gender and age classes totals similarly to the calibration
procedure applied to the original sample. The algorithm is similar to the ones
proposed by [Fabrizi et al.(2011)Fabrizi, Ferrante, Pacei, and Trivisano],
which provides estimates close to the ones obtained using linearization method
in case of simpler parameters.

Results on changes in variability depending on sample sizes are provided in
Figures 2 and 3. Notice that the analysis discriminates between the behavior
of the region whose presents a upper tail extreme values, called ”outlier
region”, see Figure 2, and the others. Specifically, bootstrap coefficient of
variation estimates on Italian NUTS-2 region are provided versus the sample
sizes in Figure 3. It is quite interesting to notice that a small-sized
extreme-valued-affected sample has always high variances estimates across all
measures. However the sensibility to small sizes changes depending on the
measure analyzed. In fact the more GE parameters $\alpha$ increases i.e. it
becomes more sensibile to upper tail values, the more the estimator variance
is affected only by the presence of upper extreme values rather than the size.
Vice versa, the more Atkinson $\varepsilon$ increases becoming more sensible
to lower tail, the more its variances becomes negatively correlated to sample
sizes.

Figure 2: Distribution of Equivalized Disposable Income for each NUTS-2
region, outlier region in orange.  Figure 3: Bootstrap Coefficient of
Variation estimates for each inequality measures.

##  7 Bias-Corrected Estimators Distribution

Literature on the probability distribution of inequality measures in small
sample is scarce. [Kakwani(1990)] and [Thistle(1990)] have studied the
asymptotic distribution of Generalized Entropy measures and Atkinson Index
proving their asymptotic normality. As regards Gini index, results on
asymptotic normality of the different estimators are well-established, see
[Giorgi and Gigliarano(2017)] for a comprehensive review. Interesting
performance comparisons between asymptotic and bootstrap inference for
inequality measures are provided by [Biewen(2002)] and [Davidson and
Flachaire(2007)].

Here, a distributional analysis has been developed in two parts. A first part
exploits empirical distributions from design-based simulation described in
Section 5, to evaluate for each measure the different rate of convergence to
well-behaved distributional shapes, considering varying sample sizes. In the
second part we perfomed a model-based simulation, using generated income data
and unequal probability sampling in a more controlled setting. The aim of the
latter one is to study the behavior of inequality measure at varying income
distributional assumptions and at the same time, to evaluate the fitting of
resulting IM distributions to assumed parametrical density with differing
responsiveness to outlier and robustness properties. We considered both mean
modelling and median modelling based distributions.

###  7.1 Design-based Simulation

Here is provided a brief analysis on the distribution of inequality measures
in samples of increasing size, to evaluate how quickly their distribution tend
to become symmetric. We considered regions as target domains and, using the
same simulation setting of Section 5, drawing $m=1000$ small samples by
mimicking the stratified multi-stage survey design with different sampling
rates, i.e. 10% (from 36 to 607 individuals, 28 to 254 households), 5% (from
16 to 337 individuals, 14 to 131 households) and 3% (from 11 to 196
individuals, 10 to 74 households). Coefficient of skewness $\eta_{3}$ and
excess kurtosis $\eta_{4}$ empirical values are set out in Table 2 for
different sampling rates. As its clear from the results, as the sample size
decreases, the empirical distributions tend to become more positively skewed
and leptokurtic. This is quite evident for the General Entropy measures,
similarly but to a lesser extent for the other measures. The empirical
distributions of the considered measures for each region are set out in Figure
4. The failing of asymptotic normality for these measures warns to not
consider for small samples parametric methods of inference based on such an
assumption.

|  | CV | GE(0) | GE(1) | GE(2) | A(0.5) | A(1) | A(2) | G
---|---|---|---|---|---|---|---|---|---
10% | $\hat{\eta_{3}}$ | 0.81 | 0.66 | 0.88 | 1.30 | 0.71 | 0.54 | 0.26 | 0.38
| $\hat{\eta_{4}}$ | 1.16 | 0.91 | 1.60 | 2.87 | 1.06 | 0.57 | -0.01 | 0.34
5% | $\hat{\eta_{3}}$ | 1.00 | 0.94 | 1.19 | 1.80 | 0.96 | 0.74 | 0.40 | 0.54
| $\hat{\eta_{4}}$ | 1.88 | 1.73 | 2.81 | 6.09 | 1.81 | 1.01 | 0.08 | 0.56
3% | $\hat{\eta_{3}}$ | 1.00 | 1.07 | 1.27 | 1.90 | 1.04 | 0.83 | 0.49 | 0.58
| $\hat{\eta_{4}}$ | 1.66 | 1.94 | 2.87 | 6.13 | 1.87 | 1.03 | 0.03 | 0.60
Table 2: Coefficients of Skewness and Excess Kurtosis for the Inequality
Measures Empirical Distributions from Design-Based simulation. Figure 4:
Inequality Measures Empirical Distributions by regions for the 3% samples from
Design-based Simulation.

###  7.2 Model-based Simulation

A model-based simulation has been also performed using the R package simFrame
[Alfons et al.(2010)Alfons, Templ, and Filzmoser], in order to provide a more
comprehensive study about the distributional form of inequality measures
estimators in small samples and complex survey context. In fact, the use of
simulated data allows us to carry on the analysis in a more controlled
setting, without potential confounding factors related to survey data
collection dynamics.

A first step involve a preliminary fitting on EU-SILC income data via weighted
pseudo-loglikelihood maximization of a plethora of suitable income
distribution such as Log-Normal, Dagum, Singh-Maddala and Generalized Beta of
the Second Kind ([Atkinson(2015)], pages 371-375). The pseudo-log-likelihood
constitutes as a weighted sum over the sample of the distribution log density,
where the weights are the sampling weights [Graf and Nedyalkova(2014)].
Maximizing it permits to take the sampling design into consideration. Results
related to some goodness-of-fit such as AIC and BIC measures are displayed in
Table 3.

The highly flexibility of the GB2 due to its four parameters allows the best
fitting in comparison to the other distributional assumptions. We decided
therefore to use the best two models (GB2 and log-normal ones) to generate
income data, in order to be able to analyse and capture differences in
estimators behavior to varying income distributional assumptions. We simulate
two different finite populations with size N=10.000, from Log-Normal
($\mu=9.64,\sigma=0.43$) and GB2 ($a=4.11,b=2.16\cdot 10^{4},p=0.47,q=0.92$)
with parameters resulting from maximum pseudo likelihood estimation.

Following the approach of [Alfons et al.(2013)Alfons, Templ, and Filzmoser],
in order to introduce a complex survey mimicking setting, we create an
auxiliary variable attached to each population unit $p=p_{1},\dots,p_{N}$
denoting probability weights. It has been constructed taking $s=100$ equally
spaced values between 1 and 10 as follows:

$\displaystyle
p_{i}=\begin{cases}1\quad\quad&x_{i}>F^{-1}_{\theta}(\frac{s-1}{s})\\\
10-\frac{9}{s-1}j\quad\quad&F^{-1}_{\theta}(\frac{j}{s})<x_{i}\leq
F^{-1}_{\theta}(\frac{j+1}{s})\quad\text{for any}\quad 1\leq j<s-1\\\
10\quad\quad&x_{i}\leq F^{-1}_{\theta}(\frac{1}{s}),\end{cases}$

with $F_{\theta}$ the cumulative distribution function of a Pareto
distribution. From each of the two populations, 3.000 samples are drawn with
three different sample sizes $n=(20,50,80)$.

The drawing has been performed by using Midzuno’s method for unequal
probability sampling [Midzuno(1952)] with inclusion probabilities proportional
to probability weights $p$. As a consequence, the observations with lower
incomes results to have higher inclusion probabilities and in turn lower
sample weights. This intuition comes from [Alfons et al.(2013)Alfons, Templ,
and Filzmoser], and it is motivated by the fact that in income survey,
generally higher inclusion probabilities are often assigned to larger
households, which typically have lower equivalized income than smaller
households. In addition, the above definition of the probability weights
induces large variation between weights. For each sample, bias-corrected
inequality estimators has been calculated using Horvitz-Thompson estimator
variances under Midzuno scheme [Narasimha Prasad and Srivenkataramana(1980)].

We transformed General Entropy measure with double-bounded support such as
Theil’s Index ($\alpha=1$) and GE($\alpha=2$) to relative Entropy measures
(RE), i.e. $\text{RE}(\alpha)=\text{GE}(\alpha)/\text{max}(\text{GE}(\alpha))$
in order to deal with measures defined on a unique support. Lastly, we
considered only measures defined on $[0,1]$ such as Atkinson Indexes, Gini
coefficient and RE, and we proceeded by fitting three distributions on
estimators samples differentiated for each fixed sample sizes $n$. Figure 5
displays how changes in family parameter values, underlying income
distribution and sample sizes, have an impact on bias-corrected estimators
distributions of the Atkinson measures. GE measures evolution depending on
varying parameter values does not show such relevant changes.

We considered the well-known Beta distribution, and some alternative
distributions defined on $(0,1)$ such as the Simplex distribution [Barndorff-
Nielsen and Jørgensen(1991)] and L-Logistic distribution. Simplex, pertaining
to the dispersion models family, is known to be an alternative to Beta
distribution in terms of over-dispersion control [Jørgensen(1997)], robustness
[Espinheira and de Oliveira Silva(2020)] and large left and right skewness
modeling [Carrasco and Reid(2019)]. The last two issues seems quite crucial
when dealing with IM estimators. It is a two parameters distribution, one
governing the mean, and the other the dispersion. On the other hand,
L-Logistic distribution, which was originally proposed by [Tadikamalla and
Johnson(1982)], through a transformation of the standard logistic
distribution, considers as parameters the median and a shape parameter. Its
choice is motivated by the fact that, if the data are skewed, since the median
is a natural robust measure of the center, the median modelling can be more
useful than mean modelling adopted in Beta and regression models, providing an
interesting alternative.

Results are set out in Table 4 with a goodness-of-fit evaluation via AIC and
BIC Information Criterions, consider that all the distributions have the same
number of parameters ($p=2$). Goodness-of-fit EDF statistics are not included
due to their unreliability in case of model parameters estimated from the
data. In all the cases robust alternatives distributions provide better
performance than Beta distribution. Whereas in case of log-normal income,
Simplex distribution provides greater fit, on the other hand in case of GB2
income, inequality estimators distributions presents large variability and
very high kurtosis values as displayed in Table 5 and thus a median-modeling
distribution such as L-Logistic seems to work better. Beta fitting is worse in
case of distributions too close to the zero-boundary, as in the case of
$RE(2)$ values with moderately large population sizes. Consider that
$\hat{RE}(2)=\hat{GE}(2)/\bigg{(}\frac{N-1}{2}\bigg{)}$, thus as long as the
true population size $N$ increases, $\hat{RE}(2)$ shrinks dramatically close
to zero-boundary, being difficult to be handled. This shrinks is less
pronounced in case of $\hat{RE}(1)=\hat{GE}(1)/\log N$.

| AIC | BIC
---|---|---
Log-Normal | 25.75 | 43.34
Singh-Maddala | 28.00 | 54.38
Dagum | 27.76 | 54.14
GB2 | -13.17 | 22.01
Table 3: Income distribution goodness-of-fit results from Model-based Simulation. Figure 5: Atkinson Indexes empirical distributions related to model-based simulation from Log-Normal income assumption (blue lines) and GB2 income assumption (orange lines), darker palette refers to decreasing sample sizes.  | Log-Normal Population | GB2 Population
---|---|---
| AIC | BIC | AIC | BIC
Beta |  |  |
RE(1) | -8719.05 | -8709.23 | -7369.32 | -7359.50
RE(2) | -10911.71 | -10901.90 | -9380.29 | -9370.47
A(0.5) | -5776.14 | -5766.33 | -4455.21 | -4445.40
A(1) | -4524.47 | -4514.66 | -3232.46 | -3222.65
A(2) | -3452.54 | -3442.72 | -1891.27 | -1881.45
G | -3487.33 | -3477.51 | -2772.32 | -2762.50
Simplex |  |  |
RE(1) | -8794.90 | -8785.08 | -7507.13 | -7497.31
RE(2) | -11019.40 | -11009.59 | -9601.15 | -9591.33
A(0.5) | -5840.08 | -5830.27 | -4550.89 | -4541.07
A(1) | -4576.11 | -4566.29 | -3280.62 | -3270.81
A(2) | -3482.09 | -3472.27 | -1901.25 | -1891.43
G | -3509.29 | -3499.48 | -2788.63 | -2778.82
L-Logistic |  |  |
RE(1) | -8783.11 | -8773.29 | -7589.79 | -7579.97
RE(2) | -11018.02 | -11008.20 | -9708.50 | -9698.68
A(0.5) | -5826.16 | -5816.34 | -4635.19 | -4625.38
A(1) | -4559.71 | -4549.90 | -3351.73 | -3341.92
A(2) | -3459.65 | -3449.83 | -1906.81 | -1896.99
G | -3496.99 | -3487.17 | -2880.50 | -2870.68
Table 4: Model-Based Simulation Results on bias-corrected estimator distributions fitting. |  | RE(1) | RE(2) | A(0.5) | A(1) | A(2) | G
---|---|---|---|---|---|---|---
| $\hat{\eta_{3}}$ | 1.68 | 2.19 | 1.47 | 1.26 | 0.88 | 0.80
logN | $\hat{\eta_{4}}$ | 4.50 | 8.02 | 3.37 | 2.38 | 1.03 | 0.77
| $\hat{\eta_{3}}$ | 3.45 | 4.63 | 2.83 | 1.98 | 0.71 | 1.50
GB2 | $\hat{\eta_{4}}$ | 19.87 | 32.92 | 14.39 | 7.61 | 0.46 | 4.07
Table 5: Mean values of Skewness and Kurtosis coefficients of the 21 regions
empirical distributions from Model-based Simulation.

##  8 Conclusions

A strategy to correct bias of inequality measures estimators in small samples
has been proposed and a sensitivity analysis has been conducted in order to
study the magnitude of the correction and its sensibility to extreme values
via a design-based simulation and EU-SILC data income for Italy. The
underlined heterogeneity of sensibilities and bias across measures can guide
analysts in choosing the most suitable inequality measure depending on the
context. Furthermore, this contribution may prove extremely useful in other
application contexts in which those measures are used. Overall, the most
famous Gini and Theil indexes are widely applied in several fields which make
use of survey data in the estimation of inequality and concentration of a
large plethora of socio-economic phenomena.

Generally, measures structurally more sensible to values on the tails appears
to be more biased, particularly GE($\alpha=2$) and Atkinson($\epsilon=2$),
reaching in some cases a bias of respectively more than -20% and more than
-45%. This problem can be circumvented, under the assumption of data
contamination via a semi-parametric Pareto-based treatment of the tails, and
bias can be corrected or reduced via bias-correction proposal set out. In
cases where extreme values cannot be considered as a consequence of data
contamination but rather constitutes as representative observations, we
suggest to use the most robust inequality measure such as GE with $\alpha=0$
and Atkinson index with $\varepsilon=1$, that, with our bias-correction, are
approximately unbiased when sample sizes $n\geq 20$.

Besides, a bootstrap variance estimation proposal and a distributional
evaluation of the corrected design-based estimators has been developed.
Results about the inequality estimators distributions show increasing positive
skewness and lepto-kurtosis at decreasing sample sizes, confirming the non-
applicability of gaussian assumption in small samples. As regards measures
defined on [0,1], model-based simulations shows that alternative robust
distributions rather than Beta have a better goodness-of-fit, paving the way
to the development of alternative parametric methods of inference when dealing
with inequality in small samples.

Further directions of research include, on one side, extending this framework
to other widely used inequality measures such as those based on quintiles, on
the other side exploring and applying variance reduction techniques to the
bias-corrected estimators in small samples in order to provide more reliable
inequality parameter estimates.

## References

  * [Alfons et al.(2010)Alfons, Templ, and Filzmoser] A. Alfons, M. Templ, and P. Filzmoser.  An object-oriented framework for statistical simulation: The R package simframe.  _Journal of Statistical Software_ , 37(3):1–36, 2010.  ISSN 15487660.  doi: 10.18637/jss.v037.i03.
  * [Alfons et al.(2013)Alfons, Templ, and Filzmoser] A. Alfons, M. Templ, and P. Filzmoser.  Robust estimation of economic indicators from survey samples based on pareto tail modelling.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 62(2):271–286, 2013.
  * [Antal et al.(2011)Antal, Langel, and Tillé] E. Antal, M. Langel, and Y. Tillé.  Variance estimation of inequality indices in complex sampling designs.  In _Proceedings 58th World Statistical Congress_ , 2011.
  * [Atkinson(2015)] A. Atkinson.  _Handbook of Income Distribution_ , volume 2.  2015\.  ISBN 9780444594303.  doi: 10.1016/B978-0-444-59428-0.09990-2.
  * [Barndorff-Nielsen and Jørgensen(1991)] O. E. Barndorff-Nielsen and B. Jørgensen.  Some parametric models on the simplex.  _Journal of multivariate analysis_ , 39(1):106–116, 1991.
  * [Bellù and Liberati(2006)] L. G. Bellù and P. Liberati.  Describing income inequality theil index and entropy class indexes.  _FAO EASYPol_ , 51, 2006.
  * [Biewen(2002)] M. Biewen.  Bootstrap inference for inequality, mobility and poverty measurement.  _Journal of Econometrics_ , 108(2):317–342, 2002\.
  * [Biewen and Jenkins(2006)] M. Biewen and S. P. Jenkins.  Variance estimation for generalized entropy and atkinson inequality indices: The complex survey data case.  _Oxford Bulletin of Economics and Statistics_ , 68(3):371–383, 2006.  ISSN 03059049.  doi: 10.1111/j.1468-0084.2006.00166.x.
  * [Breunig(2001)] R. Breunig.  An almost unbiased estimator of the coefficient of variation.  _Economics Letters_ , 70(1):15–19, 2001.
  * [Breunig and Hutchinson(2008)] R. Breunig and D. L.-A. Hutchinson.  Small sample bias corrections for inequality indices.  In W. N. Toggins, editor, _New Econometric Modelling Research_ , pages 61–83. Nova Science Publishers, 2008.
  * [Bruffaerts et al.(2014)Bruffaerts, Verardi, and Vermandele] C. Bruffaerts, V. Verardi, and C. Vermandele.  A generalized boxplot for skewed and heavy-tailed distributions.  _Statistics & probability letters_, 95:110–117, 2014\.
  * [Brzezinski(2016)] M. Brzezinski.  Robust estimation of the Pareto tail index: a Monte Carlo analysis.  _Empirical Economics_ , 51(1):1–30, 2016.  ISSN 03777332.  doi: 10.1007/s00181-015-0989-9.
  * [Carrasco and Reid(2019)] J. M. Carrasco and N. Reid.  Simplex regression models with measurement error.  _Communications in Statistics: Simulation and Computation_ , 0(0):1–16, 2019.  ISSN 15324141.  doi: 10.1080/03610918.2019.1626881.  URL https://doi.org/10.1080/03610918.2019.1626881.
  * [Ceriani and Verme(2015)] L. Ceriani and P. Verme.  Individual diversity and the gini decomposition.  _Social Indicators Research_ , 121(3):637–646, 2015.
  * [Clemenceau and Museux(2007)] A. Clemenceau and J.-M. Museux.  Eu-silc (community statistics on income and living conditions: general presentation of the instrument).  _Comparative EU statistics on Income and Living Conditions: Issues and Challenges_ , pages 11–36, 2007.
  * [Cowell and Victoria-Feser(1996)] F. A. Cowell and M.-P. Victoria-Feser.  Robustness properties of inequality measures.  _Econometrica: Journal of the Econometric Society_ , pages 77–101, 1996.
  * [Davidson(2009)] R. Davidson.  Reliable inference for the Gini index.  _Journal of Econometrics_ , 150(1):30–40, 2009\.  ISSN 03044076.  doi: 10.1016/j.jeconom.2008.11.004.  URL http://dx.doi.org/10.1016/j.jeconom.2008.11.004.
  * [Davidson and Flachaire(2007)] R. Davidson and E. Flachaire.  Asymptotic and bootstrap inference for inequality and poverty measures.  _Journal of Econometrics_ , 141(1):141–166, 2007\.  ISSN 03044076.  doi: 10.1016/j.jeconom.2007.01.009.
  * [De Santis et al.(2020)De Santis, Barabesi, and Betti] R. De Santis, L. Barabesi, and G. Betti.  Variance estimation techniques for poverty and inequality measures from complex surveys: a simulation study.  Technical report, Department of Economics, University of Siena, 2020.
  * [Deltas(2003)] G. Deltas.  The small-sample bias of the Gini coefficient: Results and implications for empirical research.  _Review of Economics and Statistics_ , 85(1):226–234, 2003.  ISSN 00346535.  doi: 10.1162/rest.2003.85.1.226.
  * [Demnati and Rao(2004)] A. Demnati and J. Rao.  Linearization variance estimators for survey data.  _Survey Methodology_ , 30(1):17–26, 2004.
  * [Deville(1999)] J. C. Deville.  Variance estimation for complex statistics and estimators: linearization and residual techniques.  _Survey methodology_ , 25(2):193–204, 1999.
  * [Espinheira and de Oliveira Silva(2020)] P. L. Espinheira and A. de Oliveira Silva.  Residual and influence analysis to a general class of simplex regression.  _Test_ , 29(2):523–552, 2020.  ISSN 18638260.  doi: 10.1007/s11749-019-00665-3.  URL https://doi.org/10.1007/s11749-019-00665-3.
  * [Fabrizi and Trivisano(2016)] E. Fabrizi and C. Trivisano.  Small area estimation of the gini concentration coefficient.  _Computational Statistics & Data Analysis_, 99:223–234, 2016.
  * [Fabrizi et al.(2011)Fabrizi, Ferrante, Pacei, and Trivisano] E. Fabrizi, M. R. Ferrante, S. Pacei, and C. Trivisano.  Hierarchical bayes multivariate estimation of poverty rates based on increasing thresholds for small domains.  _Computational statistics & data analysis_, 55(4):1736–1747, 2011.
  * [Fabrizi et al.(2020)Fabrizi, Ferrante, and Trivisano] E. Fabrizi, M. R. Ferrante, and C. Trivisano.  A functional approach to small area estimation of the relative median poverty gap.  _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ , 183(3):1273–1291, 2020.
  * [Ferrante and Pacei(2019)] M. R. Ferrante and S. Pacei.  Small area estimation of entropy inequality measures : a comparison between alternative distribution models.  In _Joint Statistical Meeting Proceedings, Survey Research Methods Section_ , Denver, Colorado, July 27 – August 1 2019. American Statistical Association.
  * [Finkelstein et al.(2006)Finkelstein, Tucker, and Alan Veeh] M. Finkelstein, H. G. Tucker, and J. Alan Veeh.  Pareto Tail Index Estimation Revisited.  _North American Actuarial Journal_ , 10(1):1–10, 2006.  ISSN 10920277.  doi: 10.1080/10920277.2006.10596236.
  * [Giles(2005)] D. E. Giles.  The bias of inequality measures in very small samples: some analytic results.  Technical report, Department of Economics, University of Victoria, 2005\.
  * [Giorgi and Gigliarano(2017)] G. M. Giorgi and C. Gigliarano.  The gini concentration index: a review of the inference literature.  _Journal of Economic Surveys_ , 31(4):1130–1148, 2017.
  * [Graf and Tillé(2014)] E. Graf and Y. Tillé.  Variance estimation using linearization for poverty and social exclusion indicators.  _Survey Methodology_ , 40(1):61–79, 2014.  ISSN 14920921.
  * [Graf and Nedyalkova(2014)] M. Graf and D. Nedyalkova.  Modeling of income and indicators of poverty and social exclusion using the generalized beta distribution of the second kind.  _Review of Income and Wealth_ , 60(4):821–842, 2014.  ISSN 14754991.  doi: 10.1111/roiw.12031.
  * [Hampel(1974)] F. R. Hampel.  The influence curve and its role in robust estimation.  _Journal of the american statistical association_ , 69(346):383–393, 1974.
  * [Jasso(1979)] G. Jasso.  On gini’s mean difference and gini’s index of concentration.  _American Sociological Review_ , 44(5):867–870, 1979.
  * [Jørgensen(1997)] B. Jørgensen.  _The Theory of Dispersion Models_ , volume 41.  1997\.  doi: 10.2307/1270747.
  * [Kakwani(1990)] N. Kakwani.  Large sample distribution of several inequality measures: With application to côte d’ivoire.  In _Contributions to econometric theory and application_ , pages 50–81. Springer, 1990.
  * [Kalton(1979)] G. Kalton.  Ultimate cluster sampling.  _Journal of the Royal Statistical Society: Series A (General)_ , 142(2):210–222, 1979.
  * [Lahiri(2003)] P. Lahiri.  On the impact of bootstrap in survey sampling and small-area estimation.  _Statistical Science_ , pages 199–210, 2003.
  * [Langel and Tillé(2013)] M. Langel and Y. Tillé.  Variance estimation of the Gini index: Revisiting a result several times published.  _Journal of the Royal Statistical Society. Series A: Statistics in Society_ , 176(2):521–540, 2013.  ISSN 09641998.  doi: 10.1111/j.1467-985X.2012.01048.x.
  * [Lerman and Yitzhaki(1989)] R. I. Lerman and S. Yitzhaki.  Improving the accuracy of estimates of gini coefficients.  _Journal of econometrics_ , 42(1):43–47, 1989\.
  * [Masseran et al.(2019)Masseran, Safari, and Ibrahim] Y. L. H. Masseran, Nurulkamal, M. A. M. Safari, and K. Ibrahim.  Power law behavior and tail modeling on low income distribution.  _Mathematics and Statistics_ , 7(3):70–77, 2019\.
  * [Midzuno(1952)] H. Midzuno.  On the sampling system with probability proportionate to sum of sizes.  _Annals of Institute of Statistical Mathematics_ , 3(2):99–108, 1952.
  * [Narasimha Prasad and Srivenkataramana(1980)] N. Narasimha Prasad and T. Srivenkataramana.  A modification to the horvitz-thompson estimator under the midzuno sampling scheme.  _Biometrika_ , 67(3):709–711, 1980.
  * [OECD(2011)] OECD.  _Divided we stand: Why inequality keeps rising_.  OECD Publishing Paris, 2011.
  * [Osier(2009)] G. Osier.  Variance estimation for complex indicators of poverty and inequality using linearization techniques.  _Survey Research Methods_ , 3(3):167–195, 2009\.  ISSN 18643361.  doi: 10.18148/srm/2009.v3i3.369.
  * [Osier et al.(2013)Osier, Berger, and Goedeme] G. Osier, Y. G. Berger, and T. Goedeme.  Standard error estimation for the eu-silc indicators of poverty and social exclusion.  _Eurostat Methodologies and Working Papers Series_ , 2013.
  * [Rao et al.(1999)Rao, Chaudhuri, Eltinge, Fay, Ghosh, Ghosh, Lahiri, and Pfeffermann] J. N. Rao, A. Chaudhuri, J. Eltinge, R. E. Fay, J. Ghosh, M. Ghosh, P. Lahiri, and D. Pfeffermann.  Some current trends in sample survey theory and methods (with discussion).  _Sankhyā: The Indian Journal of Statistics, Series B_ , pages 1–57, 1999.
  * [Rust and Kalton(1987)] K. Rust and G. Kalton.  Strategies for collapsing strata for variance estimation.  _Journal of Official Statistics_ , 3(1):69, 1987\.
  * [Safari and Ibrahim(2018)] M. N. Safari, Muhammad Aslam Mohd and K. Ibrahim.  Outliers detection for pareto distributed data.  In _American Institute of Physics Conference Series_ , volume 1940, 2018.
  * [Schluter and van Garderen(2009)] C. Schluter and K. J. van Garderen.  Edgeworth expansions and normalizing transforms for inequality measures.  _Journal of Econometrics_ , 150(1):16–29, 2009\.  ISSN 03044076.  doi: 10.1016/j.jeconom.2008.12.022.  URL http://dx.doi.org/10.1016/j.jeconom.2008.12.022.
  * [Sen(1973)] A. K. Sen.  _On economic inequality_.  Oxford: Clarendon Press, 1973.
  * [Tadikamalla and Johnson(1982)] P. R. Tadikamalla and N. L. Johnson.  Systems of frequency curves generated by transformations of logistic variables.  _Biometrika_ , 69(2):461–465, 1982.
  * [Thistle(1990)] P. D. Thistle.  Large sample properties of two inequality indices.  _Econometrica: Journal of the Econometric Society_ , pages 725–728, 1990.
  * [Ullah(2004)] A. Ullah.  _Finite Sample Econometrics_ , volume 53.  2004\.  ISBN 9788578110796.  doi: 10.1017/CBO9781107415324.004.
  * [Van Kerm(2007)] P. Van Kerm.  Extreme incomes and the estimation of poverty and inequality indicators from eu-silc.  Technical report, IRISS Working Paper Series 2007-01, CEPS/INSTEAD, 2007\.
  * [Van Ourti and Clarke(2011)] T. Van Ourti and P. Clarke.  A simple correction to remove the bias of the gini coefficient due to grouping.  _Review of Economics and Statistics_ , 93(3):982–994, 2011.