Mathias Silva
email: mathias.silva-vazquez@univ-amu.fr.

Michel Lubrano

Stephen Bazen

Emmanuel Flachaire

Philippe Van Kerm

Markus Jäntti

Jonathan Goupille-Lebret

Parametric models of income distributions integrating misreporting and non-response mechanisms

Keywords: 'Missing rich', GB2, Bayesian Inference. JEL Code: D31, C18, C11

Several representativeness issues affect the available data sources in studying populations' income distributions. High-income under-reporting and non-response issues have been evidenced to be particularly significant in the literature, due to their consequence in under-estimating income growth and inequality. This paper bridges several past parametric modelling attempts to account for high-income data issues in making parametric inference on income distributions at the population level. A unified parametric framework integrating parametric income distribution models and popular data replacing and reweighting corrections is developped. To exploit this framework for empirical analysis, an Approximate Bayesian Computation approach is developped. This approach updates prior beliefs on the population income distribution and the high-income data issues pressumably affecting the available data by attempting to reproduce the observed income distribution under simulations from the parametric model. Applications on simulated and EU-SILC data illustrate the performance of the approach in studying population-level mean incomes and inequality from data potentially affected by these high-income issues.

Introduction

The recent literature on income inequality has paid increasing attention to the dynamics and the measurement of top incomes (e.g., [START_REF] Atkinson | Top incomes over the twentieth century: a contrast between continental european and english-speaking countries[END_REF][START_REF] Leigh | The Oxford handbook of economic inequality[END_REF][START_REF] Atkinson | Top incomes in the long run of history[END_REF][START_REF] Burkhauser | Top incomes and inequality in the UK: reconciling estimates from household survey and tax return data[END_REF]. The slowly-rising availability of tax data for research purposes along with findings concerning the recent rises in the share of incomes accumulated at the very top quantiles of the distribution (e.g., [START_REF] Lakner | Global income distribution: From the fall of the Berlin wall to the great recession[END_REF]Milanovic 2016, Alvaredo et al. 2018) have jointly brought forward the multiple deficiencies affecting the typical methods and data sources used to study income distributions.

A robustly evidenced shortcoming of these conventional approaches involves the limited quality of publicly-available household survey data, the most commonly used data source on the subject, in capturing the magnitude and trends of the income shares of the highest incomes in their population (e.g., [START_REF] Deaton | Measuring poverty in a growing world (or measuring growth in a poor world)[END_REF][START_REF] Burdín | Inequality and top incomes in Uruguay: a comparison between household surveys and income tax micro-data[END_REF][START_REF] Jenkins | Pareto models, top incomes and recent trends in UK income inequality[END_REF][START_REF] Higgins | The rich underreport their income: Assessing bias in inequality estimates and correction methods using linked survey and tax data[END_REF][START_REF] Lustig | Measuring the distribution of household income, consumption and wealth: State of play and measurement challenges[END_REF].

Typically these measurement and coverage issues around the upper tail of a population's income distribution imply non-random missing information in the data (i.e., the errors are more likely or larger in magnitude for higher incomes) and therefore induce bias into any resulting distributional estimate. When ignored, this can have many clear policy implications as it leads to underestimation of income growth and inequality at the population level, along with a biased reading of their relationship and dynamics. This has motivated a vast literature on correction and estimation methods to overcome this data issue for the study of income distributions.

An important implication of this almost universal problem of missing or misreported high incomes is that any empirical strategy seeking to overcome it requires a decision on the magnitude and distribution of such errors affecting the data. As put forward by [START_REF] Bourguignon | Simple adjustments of observed distributions for missing income and missing people[END_REF], adjusting for measurement errors on high incomes requires a value for some or all of three key parameters: the income level beyond which measurement errors are to be corrected, the true population share of incomes above this level, and the share of undercovered population incomes.

Although external data sources can be instrumentally used to formulate informative choices for these parameters (e.g., [START_REF] Atkinson | Top incomes over the twentieth century: a contrast between continental european and english-speaking countries[END_REF], Chapter 2, Bustos 2015[START_REF] Institut | The weight of the rich: Improving surveys using tax data[END_REF][START_REF] Jorda | Global inequality: How large is the effect of top incomes?[END_REF], Flachaire et al. 2022), correcting for measurement or coverage errors on high incomes is conditioned by the uncertainty around them. Broadly speaking, the precision with which inference can be made on a population's income distribution depends on the uncertainty around the form and magnitude of measurement or coverage errors affecting the available data. This paper proposes a new empirical strategy bridging several previous results in the income inequality literature. Firstly, a parametric modelling approach is developped in the interest of integrating within a single framework all assumptions about the form of the population's income distribution and the form of the measurement or coverage issues affecting the available data. This parametric framework allows for exploiting several previously explored parametric corrections for high incomes data issues in making inference at the population level.

Secondly, a Bayesian estimation strategy allows for inference on the population's income distribution through data presumably affected by representativeness issues on the upper tail. This strategy extends the Approximate Bayesian Computation approach recently explored in [START_REF] Kobayashi | Approximate Bayesian computation for Lorenz curves from grouped data[END_REF] and [START_REF] Silva | Parametric estimation of income distributions using grouped data: an Approximate Bayesian computation approach[END_REF] in the context of income distributions. In exploiting this approach to estimate income distributions under the proposed parametric framework, the magnitudes and forms of the representativeness issues may be uncertain. Past knowledge on the possible nature of these under similar settings poses information that may be used in dealing with this uncertainty through the use of informative prior beliefs.

Finally, several applications over simulated and household survey data from the European Union's Statistics on Income and Living Conditions (EU-SILC) illustrate the performance of the proposed approach in controlled and observational settings. These applications evidence the several biases that hinder making inference on a population's income distribution if high-income representativeness issues affecting the available data are ignored. Additionally, the presented estimates suggest the presence of both high-income under-reporting and high-income non-response issues in selected EU-SILC samples. This results in population-level estimates of average incomes and inequality that are at higher levels and with higher uncertainty than their sample counterparts.

The following section presents an overview on the common causes and corrections for data errors on high incomes explored in the previous literature. Section 3 develops a parametric framework integrating popular forms of such data errors to parametric income distributions. The fourth section introduces an Approximate Bayesian Computation routine for inference on a population's income distribution through the proposed parametric framework and under magnitudes and forms for high-income data issues that might be uncertain. Section 5 presents simulated and EU-SILC data applications of the method under typical parametric forms. The sixth and final section of the paper presents concluding remarks with proposals for future work in studying high-income data issues through the proposed approach.

Dealing with 'missing rich' issues

In modelling income distributions, the use of parametric models is a standard. Some work has fruitfully explored the use of non-or semi-parametric methods for income distribution analysis (e.g., [START_REF] Jenkins | Did the middle class shrink during the 1980s? UK evidence from kernel density estimates[END_REF]), yet there is vast evidence of parametric models fitting real data on incomes better than these alternatives in many different settings (e.g., [START_REF] Darvas | Global interpersonal income inequality decline: The role of China and India[END_REF][START_REF] Jorda | Inequality measurement with grouped data: parametric and non-parametric methods[END_REF].

The usual modelling step involves assuming that individuals' incomes y i are distributed across its population following some three-(i.e., Θ ⊆ R 3 ) or four-parameter (i.e., Θ ⊆ R 4 ) distribution y i ∼ f y (.; θ) , θ ∈ Θ. Popular choices for f y (.; θ) include the Generalized Beta family of distributions (e.g., see [START_REF] Mcdonald | Some generalized functions for the size distribution of income[END_REF][START_REF] Jenkins | Distributionally-sensitive inequality indices and the GB2 income distribution[END_REF][START_REF] Graf | Modeling of income and indicators of poverty and social exclusion using the generalized beta distribution of the second kind[END_REF][START_REF] Chotikapanich | Using the GB2 income distribution[END_REF][START_REF] Jorda | Global inequality: How large is the effect of top incomes?[END_REF], in particular the fourparameter Generalized Beta distribution of the second kind (GB2, which is taken as illustratory reference in what follows) and the three-parameter Singh-Maddala (Burr XII) distribution, and the Double Pareto Log-Normal distribution.

There are several virtues to the parametric approach aside from its generally good fit to real data on incomes. Of particular relevance is its flexibility with respect to the format of data available. Several estimators following parametric expressions for microdata or bracketed/grouped data from incomes following y i ∼ f y (.; θ) are available for most distributions such as Generalized Method-of-Moments (GMM), Maximum Likelihood Estimation (MLE), or Bayesian inference methods.

A central consideration required in analysing data prone to high-income representativeness issues is that the distribution of observed incomes y Obs i will very unlikely follow the form of the population's income distribution y i ∼ f y (.; θ). Within a parametric approach, however, the observed distribution can be derived under assumed parametric forms for the errors affecting the data. Jointly modelling the population income distribution component f y (.; θ) and the high-income issues is an attempt at separating which aspects of the data reflect those of the population income distribution and which aspects are due to the high-income problems considered.

In describing the nature of the 'missing rich' (MR) problem, [START_REF] Lustig | The missing rich in household surveys: Causes and correction approaches[END_REF] points at the many different issues affecting the upper tail of the observed income distribution in usual data sources. In the context of survey data, the main focus of this paper, one first source of observed MR may arise from noncoverage errors in the sampling design itself as a consequence of the sparseness and irreplaceability of high income households. High-income households are generally so few and so dissimilar between themselves that households on any part of the upper tail of the distribution may have a zero probability of inclusion in the achieved survey sample.

A second possible source for MR in survey data involves reporting issues either in the form of unit or item nonresponse (i.e., high-income households refusing to respond to the survey or particularly to the items concerning their income level, respectively) or in the form of under-reporting of income levels when responding to the survey. Even if the sampling scheme is designed to be representative of the income distribution of the entire population of interest, unit or item non-response may yield an achieved sample which is not and particularly so when this nonresponse occurs more significantly for households on the upper tail of the income distribution. In a similar way, with income under-reporting the achieved survey sample may yield an income distribution which is not representative of the population's true income distribution if under-reporting is particularly present for high-income households.

Finally, a third possible source can be found within the data provision procedures commonly used by the institutions in charge of distributing publicly-available household survey datasets. In the interest of statistical disclosure control, it is common for such publicly-available datasets to contain a measure of household incomes which is top-coded (i.e., right-censored) meaning that reported incomes above a certain threshold cannot be observed as measured and only an indicator of being above this income threshold is presented.

Deriving parametric distributions for the sample distribution of incomes y Obs i observed under simple forms of measurement errors is the focus of early literature in the field. Models obtained from simple two-parameter distributions 'distorted' through classical measurement errors (i.e., independent of incomes) brought forward implications that would be in strong contrast with recent empirical observations: classical measurement errors can yield sample inequality estimates that overestimate inequality at the population level (e.g., see [START_REF] Krishnaji | Characterization of the Pareto distribution through a model of underreported incomes[END_REF][START_REF] Hartley | On the estimation of the Pareto law from under-reported data[END_REF][START_REF] Hinkley | Estimation of the Pareto law from underreported data: A further analysis[END_REF][START_REF] Van Praag | The influence of classification and observation errors on the measurement of income inequality[END_REF][START_REF] Ransom | Income distribution functions with disturbances[END_REF][START_REF] Chesher | Welfare measurement and measurement error[END_REF].

More recent literature, in change, has focused in characterizing under-reporting phenomena affecting income data. The robustly evidenced progresiveness of under-reporting with respect to income levels has resulted in more appropriate nonclassical parametric expressions for these measurement errors (e.g., see [START_REF] Gottschalk | Are earnings inequality and mobility overstated? the impact of nonclassical measurement error[END_REF][START_REF] Bourguignon | Simple adjustments of observed distributions for missing income and missing people[END_REF][START_REF] Blanchet | The weight of the rich: improving surveys using tax data[END_REF], and Flachaire et al. 2022) and has consistently found that high-income under-reporting yields sample inequality measures that underestimate inequality at the population level. This recent exploration of high-income under-reporting has given way to what are known as replacing corrections: incomes pressumed to be under-reported in the data are replaced by imputations from external data sources such as administrative tax data or by imputations from a model for the under-reporting mechanism. The 'corrected' data is then treated as a representative sample of the population's incomes following f y (.; θ).

Another source of MR issues, that of missing observations, has also been treated under parametric approaches. The case of item non-response has received significantly more treatment than the more complex case of unit non-response (e.g., see [START_REF] Brunori | Estimating inequality with missing incomes[END_REF] for a recent survey). The main aspect determining how to proceed concerns the distinction between observations Missing Completely at Random (MCAR), Missing at Random (MAR), or Missing Not at Random (MNAR), following the works of [START_REF] Rubin | Inference and missing data[END_REF], [START_REF] Rubin | Formalizing subjective notions about the effect of nonrespondents in sample surveys[END_REF], and [START_REF] Greenlees | Imputation of missing values when the probability of response depends on the variable being imputed[END_REF]. In the MCAR case the probability of a unit/item being missing in the data is independent of any characteristics of the unit and is constant across all units, inducing no particular biases to any distributional estimates from the data. The MAR case allows this probability to change with the characteristics of the unit but requires it to be independent of the unit's income level. Finally, the more complex MNAR case allows this probability to also change with the unit's income level and is therefore the only mechanism capable of representing the empirically evidenced negative relationship between response probabilities and income levels in survey data (e.g., [START_REF] Bollinger | Trouble in the tails? what we know about earnings nonresponse 30 years after Lillard, Smith, and Welch[END_REF][START_REF] Hlasny | Nonresponse bias in inequality measurement: Cross-country analysis using Luxembourg Income Study surveys[END_REF].

The biases introduced by MAR or MNAR missing data mechanisms in distributional analysis have mostly been treated under the assumption that unit/item missingness is due to non-response. Namely, the assumption that conditional on being sampled high-income units are less likely to report their incomes (in the case of item non-response) or any information at all (in the case of unit non-response) than other units. This approach has motivated the use of reweighting corrections: the empirical distribution of incomes in the sample is reweighted by the related distribution of (imputed) response probabilities. Like with replacing, the reweighted data is then treated as a representative sample of the population's incomes following f y (.; θ).

A particularly lacking aspect of the recent replacing/reweighting approaches in deal-ing with MR is the lack of unified parametric frameworks integrating the modelling assumptions on the income distribution f y (.; θ) and those on the under-reporting and/or non-response mechanisms. This has several consequences on the applicability and generalizability of these methods.

As a model for the data directly as it is observed, a unified parametric approach can allow for deriving expressions and estimation strategies suitable for microdata but also for other formats such as bracketed or grouped data (under known grouping mechanisms). Additionally, this may prove useful in dealing with the challenge that recent approaches in the literature face concerning the choice of correction quantities (i.e., the share of missing and/or under-reported incomes and their distribution). In general these quantities are hand-set by the analyst or are set to match quantities given from more reliable external data.

While setting correction quantities ad hoc relies entirely on the analyst's knowledge about the population's true income distribution, setting these quantities taking external data as reference poses several issues of its own. Firstly, it is not always the case that more reliable external data sources on incomes are available for research purposes as there may be access restrictions to such data or the data may suffer from MR issues of their own such as those induced by tax evasion and tax avoidance on administrative tax data. Secondly, even when external data is available it is generally the case that the population coverage and income components covered differ from those in the primary data available for the analysis and this implies that several compatibilizations must be made in transfering quantities from the former to complement the latter. This compatibilizations often come at the cost of forcing different income concepts to represent the same and of a loss in being able to quantify the statistical uncertainty around the resulting distributional estimates and in particular how these are affected by the inherent uncertainty concerning the correction quantities.

The parametric framework proposed in what follows builds on the recent literature exploring replacing and reweighting corrections for MR by integrating in a single distribution function both the population income distribution model f y (.; θ) and any assumed form for measurement or missing data issues affecting incomes. Several formats of data and inference may be analyzed through the scope of this framework including that of learning about plausible values for the different MR correction quantities from the data itself and of integrating the uncertainty around these quantities into distributional estimates such as the Gini coefficient.

A parametric replacing and reweighting framework

Let individual i's true income be denoted by y i ∼ f y (.; θ), with probability density function (pdf) f y (.; θ) and cumulative distribution function (CDF) F y (.; θ), and consider a sample of observed individual incomes y Obs = {y Obs i } N i=1 . This sample may be affected by two types of 'missing rich' issues: high-income under-reporting, in which case income y Obs i is observed but differs from y i following under-reporting of high incomes, and nonresponse, in which case no income is observed for the individual i.

To introduce a parametric model for data under both of these possible issues, let ϕ(y, X; ν) a response probability function defining the probability for an individual to report her income after being sampled from the population. In its most general formulation this probability, parametrized by the vector ν, may depend on the individual's income y i and/or other characteristics X i,. but also on others' incomes y and/or characteristics X more generally. Additionally, denote by m(y, X; η) an income reporting function, defining the link between i's income y i and her income reported in the data, if any, y Obs i ≡ m(y, X; η). In particular, consider this latter function to be invertible, such that y i ≡ m -1 (y Obs , X; η) is a replacing function.

Within this framework, we can relate i's income to her observed income y Obs i , if reported, and to some unobservable income y N Obs i , in case of non-response, following 1 :

y i = m -1 (y Obs i , X i ; η) , with probability ϕ(y i , X i ; ν) y N Obs i , with probability 1 -ϕ(y i , X i ; ν)
If no measurement or non-response issues are believed to affect the data, then this amounts to setting (m -1 (y Obs i , X i ; η), ϕ(y i , X i ; ν)) ≡ (y Obs i , 1).

Whenever some form of measurement error is assumed to affect incomes in the data, then this may be introduced through a specific choice for the replacing function m -1 (y Obs i , X i ; η). This function serves the purpose of introducing any replacing or imputation step where i's income is set as a function of her observed income.

Any m -1 (y Obs i , X i ; η) representing progressive under-reporting of high incomes should imply an increasing and convex quantile ratio r(i; η) defined as:

r(i; η) = m -1 (y Obs (i) ; η) y Obs (i) , ∂r(i; η) ∂y Obs (i) ≥ 0 , ∂ 2 r(i; η) ∂ 2 y Obs (i)
≥ 0 with y Obs (i) denoting the i-th quantile of y Obs . This restricts relative discrepancies between observed y Obs i and replaced m -1 (y Obs i ; η) incomes to be non-decreasing with incomes.

Recently popular replacing approaches can easily be expressed as deterministic forms for m -1 (y Obs i , X i ; η) including:

Piecewise linear quantile-ratio replacing (e.g., [START_REF] Flachaire | Underreporting of top incomes and inequality: A comparison of correction methods using simulations and linked survey and tax data[END_REF]):

m -1 (y Obs (i) ; {p k } K-1 k=1 , {β k } K-1 k=1 , {δ k } K-1 k=1 ) ≡          y Obs (i) , if p (i) ≤ p1 y Obs (i) × K-1 k=1 1(p k < p (i) ≤ pk+1 ) × (β k + δ k p (i) )
Linear replacing weights for incomes in the k-th segment.

with δ j ≤ δ j+1 < ∞ , j = 1, ..., K -2, pK = 1, and with p (i) = F y (y (i) ; θ) denoting ordered-incomes individual y Obs (i) 's percentile in the population's income distribution 2 . In absence of missing data, the sample percentile p Obs (i) ≡ (i) N , (i) = 1, ..., N is 1 For simplification reasons, all derivations in what follows are under the assumption that i's both response probabilities and reported income depend only on i's income and characteristics: m -1 (y Obs , X; ν) ≡ m -1 (y Obs i , X i ; ν), ϕ(y, X; ν) ≡ ϕ(y i , X i ; ν) 2 In what follows 1(.) represents the identity function, taking value 1 whenever the condition it takes as argument holds true and 0 otherwise. also a valid estimate of F y (y (i) ; θ). The central assumption under this approach is that the quantile ratio r(

i; η) = m -1 (y Obs (i) ;{p k } K-1 k=1 ,{β k } K-1 k=1 ,{δ k } K-1 k=1 ) y Obs (i)
can be represented as a continuous piecewise linear function. This piecewise representation allows for progressive under-reporting of high incomes across segments (p k ; pk+1 ] of the income distribution at the cost of introducing 3 additional parameters (p k ; β k ; δ k ) per segment.

Linear progressive under-reporting (LPU, Bourguignon 2018):

m -1 (y Obs i ; p, δ) ≡ y Obs i ×         1 + 1(y Obs i > F -1 y (p; θ))
       
with δ ∈ [0, 1) and with F -1 y (p; θ) denoting the p-th population income quantile. This replacing scheme assumes that all individuals with incomes above the p-th percentile under-report their incomes in the observed sample and do so in a linearly progressive manner with under-reporting increasing by δ with every additional unit of income. The incomes quantile ratio implied under LPU r(

i; η) = m -1 (y Obs (i) ;p,δ) y Obs (i)
m -1 (y Obs i ; µ, σ, ζ) ≡ y Obs i ×            1 + 1(y Obs i > µ) Indiv. with observed incomes above µ under-report ×       1 -p i -p 1-p -ζ -1 ζ    × σ -(y Obs i -µ)   
Observed incomes are replaced by corresponding percentile under a Generalized Pareto dist.

          
where (µ, σ, ζ) are respectively the location, scale, and shape parameters of a Generalized Pareto distribution GPD(µ, σ, ζ) with CDF given by:

F y (y i ; µ, σ, ζ) =    1 -1 + ζ(y i -µ) σ -1 ζ , if ζ = 0 1 -e -( y i -µ σ ) , if ζ = 0 , y i > µ
Under this replacing scheme, any income above a level µ is assumed to be underreported. True incomes are assumed to follow a Generalized Pareto distribution (or some specific case such as the Pareto I (i.e., GPD( σ ζ , σ, ζ) , ζ > 0) or Pareto II (i.e., GPD(µ, σ, ζ) , ζ > 0)) above this income level. In absence of missing data issues an individual in the p i -th sample percentile with p i > p (equivalently, y Obs i > µ) has true income corresponding to the p i -p 1-p -th quantile on this Pareto distribution. Similarly to LPU, the incomes quantile ratio implied under Generalized

Pareto replacing r(i; η) = m -1 (y Obs (i) ;µ,σ,ζ) y Obs (i)
is strictly convex for income levels above µ, representing progressiveness of the under-reporting.

These common replacing schemes all exploit the assumption that under-reporting is a deterministic function of individual incomes (or their sample percentile/rank equivalently), and that individuals have the same rank in the population's income distribution as in the observed sample. It's also important to note that each specific replacing scheme implies within it specific assumptions on under-reporting behaviour at the individual level.

Figure 1 illustrates a comparative example of the quantile ratios r(i, η) under these three common forms for m -1 (.; η). The respective parameter values η are set to represent a same progressive under-reporting pattern: LPU affecting observed incomes from the .75th percentile of the income distribution upwards. A first observation is that a piecewise linear approximation to this under-reporting pattern introducing six parameters in η in total (i.e., a linear approximation with two segments) is not flexible enough to correctly represent it. Secondly, replacing under a Generalized Pareto tail all incomes above the .75-th percentile can represent the reference LPU pattern accurately, except for the top of the income distribution. Finally, this similarity across GPD and LPU schemes suggests the latter as the more stable and parsimonious alternative of the two. Note: Three common replacing schemes: LPU, piecewise linear quantile ratio, and Generalized Pareto replacing (GP in legend), as applied to a same income distribution following y i ∼ GB2(2.257, 17393, 1, 1.033) (see following sections for details on the GB2 distribution) and affected by LPU with p = .75 (represented by the dashed vertical line) and δ = .67. The piecewise linear approximation was calibrated to fit this LPU pattern at the p1 = .75 and p2 = .9375 sample percentiles. The GPD coefficients ζ and σ were estimated conditional on µ being the .995-th sample quantile as a typical empirical practice (e.g., see Jenkins ( 2017)) and imposing finite variance (ζ < 1 2 ).

Further choices for m -1 (y Obs i , X i ; η) may exploit other individual characteristics X i to represent larger heterogeneities in under-reporting patterns. A popular choice when data on individual consumption is available without reporting errors of its own is to define the income reporting function as an Engel curve (e.g., see [START_REF] Pissarides | An expenditure-based estimate of Britain's black economy[END_REF][START_REF] Lyssiotou | Estimates of the black economy based on consumer demand approaches[END_REF][START_REF] Hurst | Are household surveys like tax forms? evidence from income underreporting of the self-employed[END_REF]. Additionally, a stochastic component may be introduced in the definition of m -1 (y Obs i , X i ; η) to allow for heterogeneity in underreporting behaviour across individuals with same level of incomes (e.g., [START_REF] Flachaire | Underreporting of top incomes and inequality: A comparison of correction methods using simulations and linked survey and tax data[END_REF].

Concerning the modelling assumptions for the response probabilities ϕ(y i , X i ; ν), the several possible types of missing data mechanisms may be considered, following [START_REF] Rubin | Inference and missing data[END_REF]. If income non-response follows a random process which is unrelated to incomes y and other characteristics X, then the mechanism corresponds to a MCAR process. In the MCAR case, we observe incomes y Obs = {y Obs i } N i=1 which are a random sample of the population's incomes and therefore no particular bias is induced by the missing data. A simple MCAR mechanism is such that ϕ(y i , X i ; ν) ≡ ϕ(y i , X i ; p) ≡ p , p ∈ (0, 1], where all individuals are just as likely to report incomes after they have been sampled.

A second potential mechanism concerns the case where non-response in incomes are not random but where the misingness can be fully explained by other non-missing characteristics of the individuals and/or by the observed incomes, i.e., ϕ(y i , X i ; ν) ≡ ϕ(y Obs i , X i ; ν). This mechanism represents an MAR process and is an appropriate representation for scenarios of item non-response, where sampled individuals report information about their characteristics X i but not about their income, as long as their unobserved income y N obs i is unnecessary to account for the non-random non-response probabilities. MAR mechanisms are usually dealt with in analysis through multiple imputations of incomes for those individuals in the data with missing incomes but observed characteristics.

Finally, it may be the case that response probabilities may not be fully accounted for from observed information. For instance, it may be the case that the reason why individuals do not report their incomes in the data has everything to do with their unobserved level of incomes y N Obs i . This corresponds to the MNAR scenario and is particularly complex to deal with, as it may include non-random unit non-response mechanisms, where sampled individuals do not report neither incomes nor characteristics and where their unobserved incomes y N Obs i are a determinant of this. Forms for ϕ(y i , X i ; η) suitable for MAR mechanisms have been the focus of the recent survey in [START_REF] Brunori | Estimating inequality with missing incomes[END_REF]. Recently popular reweighting approaches allowing for dealing also with MNAR mechanisms can easily be expressed as deterministic forms for ϕ(y i , X i ; η) including:

Right-truncation (e.g., Alvaredo 2011, Jorda and Niño-Zarazúa 2019):

ϕ(y (i) ; t, α) ≡ α , if p (i) ≤ t 0 , if p (i) > t
which amounts to assuming that any and all individuals above the t-th percentile on the population income distribution will not report incomes in the data, while anyone below this threshold will report an income with probability α. The limiting case α → 1 corresponds to assuming that any unit with income below the t-th percentiles will always report an income when sampled.

Regional non-response reweighting (e.g., [START_REF] Korinek | An econometric method of correcting for unit nonresponse bias in surveys[END_REF], Hlasny and Verme 2018):

ϕ(y i , X i ; β) ≡ e g(y i ,X i ;β) 1 + e g(y i ,X i ;β) with g(y i , X i ; β) being a twice continuously differentiable function of observed unit i's characteristics parametrized by the vector β. The comparative analysis in [START_REF] Hlasny | Top incomes and the measurement of inequality: A comparative analysis of correction methods using Egyptian, EU, and US survey data[END_REF] suggests a simple logarithmic form for g taking as input a linear combination of income y i and region indicator variables to be equally efficient as more complex specifications in many scenarios. This approach infers response probabilities for units from modelling the relationship between non-response rates and units' characteristics at aggregate (i.e., regional) levels, when this information is available. The key assumption is that individual characteristics relate to individual response probabilities in the same way that they do at the aggregate level (i.e., that ecological inference is feasible), such that individual response probabilities ϕ(y i , X i ; β) may be properly estimated and used for the purpose of reweighting the observed data.

Income-proportional reweighting (e.g., [START_REF] Blanchet | The weight of the rich: improving surveys using tax data[END_REF]:

ϕ(y i ; γ 0 , γ 1 , t, α) = e γ 0 (y i ) -γ 1 , γ 1 > 0 , if p (i) > t α , if p (i) ≤ t
This scheme corresponds to a non-response mechanism where individuals with true incomes above the t-th percentile have increasingly lower response probabilities, with the parameter γ 1 representing the income elasticity of non-response (i.e., how much response probabilities decrease with an increase in incomes of 1%). For a given value of such elasticity, γ 0 serves as an intercept to assure the continuity of ϕ(y i ; γ 0 , γ 1 , t, α) at t. Similarly to right-truncation, the parameter α represents the (constant) response probability for units with incomes below the t-th percentile. Moreover, this reweighting scheme includes the right-truncation ϕ(y (i) ; t, α) as the limiting case γ 1 → ∞. Figure 2 provides an illustratory example of how these two cases relate and their resulting contrasts with the population's income distribution.

Both replacing and reweighting corrections may interact within this framework in a straightforward manner. Response probabilities are modelled as a function of true incomes y, yet these may differ from observed incomes y Obs following an assumed form for m -1 (y Obs i , X i ; η). The interaction allowing for modelling MNAR non-response mechanisms through observed incomes directly simply amounts to the composite function ϕ(m -1 (y Obs i , X i ; η), X i ; ν).

A crucial question this parametric framework allows to answer is: given i) an assumed form for the population income distribution f y (.; θ), ii) an assumed income reporting form m(y i , X i ; η), and iii) an assumed response probability function ϕ(y i , X i ; ν), then what distribution f y Obs (y Obs i ; θ, η, ν) will observed incomes under this framework follow? This distribution can be obtained applying the deterministic transformation m -1 (.; η) to f y (.; θ) and reweighting the resulting density by the response probabilities ϕ(.; ν), yielding the relationship3 :

f y Obs (y Obs i ; θ, η, ν) = Reporting function: Replacing transformation of y f y (m -1 (y Obs i ; η); θ) × ∂m -1 (y Obs i ; η) ∂y Obs i × Non-response: Reweighting of fy ϕ(m -1 (y Obs i ; η); ν) f y (m -1 (y Obs i ; η); θ) × ϕ(m -1 (y Obs i ; η); ν) × ∂m -1 (y Obs i ; η) ∂y Obs i dy Obs Normalizing constant (1)
The main application of the result in (1) is that of parametrically integrating all assumptions about the population income distribution and the 'missing rich' issues affecting the data in a model for the observed data itself. This model for the data constitutes therefore a standard case of continuous model expansion to accomodate for non-response or measurement errors (e.g., see [START_REF] Nandram | Hierarchical Bayesian nonresponse models for binary data from small areas with uncertainty about ignorability[END_REF][START_REF] Gustafson | On model expansion, model contraction, identifiability and prior information: Two illustrative scenarios involving mismeasured variables[END_REF], Gelman et al. 2013, Chapter 7,). Fitting such a model to data is an attempt at identifying separately features characteristic of the population income distribution, captured by the θ vector, and features representing the 'missing rich' forms pressumed to affect the data, captured by the η and ν vectors.

There are several virtues to integrating the replacing and/or reweighting corrections considered relevant into a model to be taken to the data as-is. Firstly, because the correction quantities are completely defined through η and/or ν this approach guarantees that all corrections are done on the income concept and population being analyzed. This avoids the issue of manipulating these concepts to be compatible with correction quantities defined in terms of different income concepts or population. On a related note, if external data informative on the forms and magnitudes of 'missing rich' are available then these should be introduced by specifying adequate representantions m -1 (y Obs i ; η) and ϕ(y Obs i ; ν) and setting η and ν to quantify these magnitudes.

A second virtue of this integrated approach is that it makes it feasible to compute measures of uncertainty, such as standard errors, for the estimated parameter values. Importantly, these measures of uncertainty can allow for probabilistic assessment of the relevance of any corrections considered, as well as allowing for considering the uncertainty in the estimated features of the population's income distribution that also considers the uncertainty surrounding the corrections quantities.

Thirdly, the 'building blocks' nature of the framework allows for exploring several candidate forms for replacing and/or reweighting corrections leaving other components unchanged in a straightforward manner. In particular, this allows for studying the robustness of the estimated θ to different assumptions on the form of 'missing rich' affecting the data.

Finally, stating the model as a properly defined parametric distribution implies that all observed units are re-weighted under any assumed form for φ(.; ν), either directly through the reweighting of units prone to non-response in the numerator of (1) (i.e., through downweighting the density at their respective level of income with respect to the population density.) or indirectly through the correction for missing observations in the normalization constant of (1). This 'indirect' reweighting accomodates for the fact that if some units are under-represented in the data due to higher non-response probabilities then necessarily the rest of units are over-represented and therefore need to be reweighted under any correction for these non-response probabilities.

Making inference about the features of the population income distribution and the 'missing rich' aspects of the data simultaneously poses several challenges. Issues of identifiability, in particular, require attention as a given model specified following (1) might fit equally well a sample of observed incomes for very different values of the θ , η , and ν parameters, making inference on them invalid. The type of continuous model expansion underlying (1) to introduce uncertainty about the specific form and magnitudes of the 'missing rich' issues affecting the data falls in line with previous empirical strategies within Bayesian inference (e.g., see [START_REF] Nandram | Hierarchical Bayesian nonresponse models for binary data from small areas with uncertainty about ignorability[END_REF]. The use of prior probabilities on parameter values under a Bayesian approach can overcome some identifiability issues. The following section details a Bayesian inference approach for this purpose which can exploit external information on the 'missing rich' correction quantities in dealing with this.

Parameter inference under 'missing rich'

Under the framework developed in the previous section, inference on the population's income distribution f y (.; θ) is made through inference on the values of the θ vector given the sample of observed incomes y Obs . This task is considerably less complex whenever the correction quantities η and ν are given fixed values. However, it is rarely the case that sound candidate values for these quantities are available. The central challenge in making inference on θ is therefore to exploit the framework under an empirical strategy that can properly estimate these parameters but also η and ν at the same time.

The goal is to learn about which values for the parameters η) , and ν ∈ Θ ν ⊆ R dim(ν) are more likely to have generated the observed data y Obs than others within some region of possible values

θ ∈ Θ θ ⊆ R dim(θ) , η ∈ Θ η ⊆ R dim(
Θ ≡ Θ θ × Θ η × Θ ν .
In the Bayesian framework, this information takes the form of a posterior probability distribution π(θ, η, ν|y Obs ) defined by two main components under Bayes' theorem. Firstly, all prior beliefs about the values of the (θ, η, ν) parameters must be ellicited through a prior probability distribution p(θ, η, ν) over Θ. Secondly, for any fixed value for the parameters ( θ, η, ν) the model's likelihood L(y Obs | θ, η, ν) quantifies how likely the observed data y Obs is to have been generated from f y Obs (.; θ, η, ν) 4 . π(θ, η, ν|y Obs ) is then a probability distribution proportional to the prior probability distribution updated (or reweighted, equivalently) by the likelihood function:

π(θ, η, ν|y Obs ) ∝ L(y Obs |θ, η, ν) × p(θ, η, ν) (2)
As an evidence-weighted conversion of prior beliefs, the information contained in the π(θ, η, ν|y Obs ) posterior distribution can be interpreted as all remaining uncertainty on the values of (θ, η, ν) after having 'learnt' from the data through the likelihood L(y Obs |θ, η, ν). Whenever the data are informative about these parameters, the posterior distribution reflects less uncertainty around their values than that in p(θ, η, ν).

Estimating a posterior distribution π(θ, η, ν|y Obs ) for the model parameters faces several complexities. As is usual in most Bayesian inference settings, it is rarely the case that π(θ, η, ν|y Obs ) admits a known form given a model L(y Obs |θ, η, ν) and a prior p(θ, η, ν). This is typically circumvented by studying the posterior distribution through samples generated to converge to π(θ, η, ν|y Obs ) under the Monte Carlo principle 5 or the Markov Chain Monte Carlo (MCMC) extension of this principle (e.g., see Gelman et al. 2013, Chapter 11).

A second complexity in estimating π(θ, η, ν|y Obs ) concerns the possible 'nonidentifiability' of at least some of the parameters in (θ, η, ν). As an illustratory example of this issue, consider a model specified following (1) with a parameter λ θ ∈ θ ruling the right tail of the f y (.; θ) income distribution and a replacing correction m -1 (.; η) with parameter λ η ∈ η also affecting only the right tail. It can be the case that a same 4 For example, in the case of y Obs = {y Obs i } N i=1 being N independent observations their joint likelihood follows L(y Obs |θ, η, ν) = N i=1 f y Obs (y Obs i ; θ, η, ν) 5 The Monte Carlo principle states that any quantity of π(θ, η, ν|y Obs ) which can be expressed as an expectation can be studied through a sufficiently large sample of J independent draws {( θ(j) , η(j) , ν(j) )} J j=1 from this distribution ( θ(j) , η(j) , ν(j) ) ∼ π(θ, η, ν|y Obs ) sample of incomes y Obs may be equally well fit under two different parameter values ( θ, η, ν) ∈ Θ and ( θ , η , ν ) ∈ Θ including ( λθ , λη ) and ( λ θ , λ η ) respectively. This can render the model incapable of separately identifying variations in high incomes in y Obs that would occur with changes in λ θ and those due to λ η . If y Obs is not informative about differences in the respective likelihoods L(y Obs | θ, η, ν) and L(y Obs | θ , η , ν ), then prior beliefs on these values will not be updated. The respective posterior probabilities π( θ, η, ν|y Obs ) and π( θ , η , ν |y Obs ) will therefore be dominated entirely by differences in prior beliefs p( θ, η, ν) and p( θ , η , ν ). If available, external information about the plausibility of ( θ, η, ν) and ( θ , η , ν ) may be exploited to set an informative prior distribution p(θ, η, ν) giving a lower prior probability to the one set of parameter values less compatible with this external data amongst the two. Informative priors are a way of exploiting prior knowledge to justify differences in posterior densities for parameter values where y Obs is uninformative through the model L(y Obs |θ, η, ν). Returning to the illustratory example above, consider an application of (1) as a model for a survey's sample on incomes y Obs = {y Obs i } N i=1 integrating a GB2 income distribution f GB2 y (.; θ) , θ = (α, β, p, q) with an LPU form for m -1 (.; η) , η = (p, δ). LPU affects only the tail above the p-th percentile of the income distribution, while the p and q parameters of the GB2 distribution rule its right tail. This allows for identifiability issues as described above, as there might be configurations 'trading' values of p and q with values of p and δ while representing two observably identical income distributions. In this example, external information might be introduced in the form of prior probabilities by setting the marginal prior distributions for p and δ around previous empirical findings on 'missing rich' issues in similar settings6

Several sampling algorithms can be devised to obtain samples {( θ(j) , η(j) , ν(j) )} J j=1 from π(θ, η, ν|y Obs ) under a model following (1) and an informative prior p(θ, η, ν).

The Metropolis-Hastings (MH) algorithm defines a type of MCMC sampler suitable for estimating parametric income distribution models in several contexts (e.g., see [START_REF] Chotikapanich | Applications: posterior distributions for the Gini coefficient using grouped data[END_REF][START_REF] Peters | Bayesian inference, monte carlo sampling and operational risk[END_REF][START_REF] Chotikapanich | Estimating income distributions using a mixture of Gamma densities[END_REF].

A standard MH sampler for the joint parameter vector φ = (θ, η, ν) is possible following algorithm (1) below. Such an MH algorithm yields as output a sample { θ(j) , η(j) , ν(j) } J j=1 resulting from a global exploration of the support of π(θ, η, ν|y Obs ) through local accept-reject steps. Any j-th, j = 1, ..., J, local accept-reject step is defined by the MH acceptance probability:

ρ (j) = min 1, π( φ(j) |y Obs ) × g( φ(j-1) , φ(j) ) π( φ(j-1) |y Obs ) × g( φ(j) , φ(j-1) ) , φ = (θ, η, ν)
with g( φ(j) , φ(j-1) ) denoting a candidate function from which the j-th candidate value φ(j) is sampled, given the previously retained value φ(j-1) .

Algorithm 1: A Metropolis-Hastings algorithm (MH ). 1: Initialization: 2: Until L(y Obs |φ (0) ) > 0:

3:

Sample φ(0) from p(φ)

4: Sampling: for j = 1, ..., J do 5: Sample φ(j) ∼ g(φ, φ(j-1) ) from the candidate g 6: Accept and store φ(j) with probability: j) where u (j) is a draw from a U nif orm(0, 1) distribution otherwise store φ(j) = φ(j-1) end A common choice of candidate function is that of the Adaptive Random-Walk Metropolis (AM ) algorithm [START_REF] Haario | An adaptive Metropolis algorithm[END_REF]. In this case the proposal g ≡ g Σ is defined by the following adaptive random walk process:

ρ (j) = min              1, ∝π( φ(j) |y Obs ) under (2) L(y Obs | φ(j) ) × p( φ(j) ) × g( φ(j-1) , φ(j) ) L(y Obs | φ(j-1) ) × p( φ(j-1) ) × g( φ(j) , φ(j-1) )              e.g., if u (j) ≤ ρ (
(θ (j) , η (j) , ν (j) ) ≡ φ (j) ∼ g Σ (j-1) (φ, φ(j-1) ) ⇒ φ(j) = φ(j-1) + ˜ (j) (j) ∼ N d (0, Σ (j-1) ) Σ (j-1) = Σ (0) , if j ≤ J 0 s d × 1 (j-1) (j-1) i=1 φ(i) φ(i) -i × φ φ + s d × χ × I d , if j > J 0 , 0 < χ 1
with φ denoting the mean value of all draws up to and including the (j -1)-th and with s d suggested, following [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF], to be set to s d = 2.4 2 d where d is the number of parameters in φ7 . Under this proposal distribution the j-th candidate value φ(j) is obtained by sampling from a multivariate Gaussian distribution centered at the previously retained draw φ(j-1) and with covariance matrix Σ (j-1) . Being initially set to a given matrix Σ (0) , this covariance matrix starts adapting exploiting all past draws after a sufficiently large initial period J 0 following the sample covariance matrix. An (AM ) algorithm can thus focus on sampling more densely in regions near values φ with high posterior density and less densely in regions of low posterior density. It is also possible to extend the scope of the local accept-reject exploration by sampling M candidates at once from g Σ (j-1) (φ, φ(j-1) ) in the spirit of the multiple-try Metropolis sampler of [START_REF] Liu | The multiple-try method and local optimization in Metropolis sampling[END_REF].

Implementing the (AM ) algorithm requires being able to compute the likelihood function L(y Obs |θ, η, ν).

For a model following (1), the joint likelihood for a sample of independent microdata y Obs = {y Obs i } N i=1 can be computed as L(y Obs |θ, η, ν) = N i=1 f y Obs (y Obs i ; θ, η, ν). However, a computable likelihood function may not be available in several contexts, such as when data on incomes is only available at group level (e.g., see [START_REF] Kobayashi | Approximate Bayesian computation for Lorenz curves from grouped data[END_REF]Kakamu 2019, Eckernkemper and[START_REF] Eckernkemper | Classical and bayesian inference for income distributions using grouped data[END_REF]. A more general MCMC sampling algorithm, requiring only being able to simulate data from a model following ( 1) is available through the Approximate Bayesian Computation (ABC) approach (e.g., see [START_REF] Kobayashi | Approximate Bayesian computation for Lorenz curves from grouped data[END_REF]Kakamu 2019, Silva 2023).

ABC constitutes a class of simulation-based Bayesian inference methods which approximate the unavailable likelihood L(y Obs |θ, η, ν) through simulating data. This approximation requires a way to assess how closely the observed data y Obs resembles data simulated from the model ỹObs ∼ f y Obs (.; θ, η, ν) for any given parameter values ( θ, η, ν) ∈ Θ8 . Following [START_REF] Silva | Parametric estimation of income distributions using grouped data: an Approximate Bayesian computation approach[END_REF], this comparison may be done by representing y Obs and ỹObs through their respective empirical Generalized Lorenz curves (GLC, Shorrocks 1983) defined as9 Given the observed-data GLC, denoted {GLC Obs k } N k=1 , and the analogue GLC from a simulated sample, denoted {GLC Obs k ( θ, η, ν)} N k=1 , the overall degree of discrepancy between both empirical income distributions may be summarized by the following unidimensional metric:

GLC Obs k = k i=1 y Obs (i) N i=1 y Obs (i) s Obs k × 1 N N i=1 y Obs (i) µ Obs = k i=1 y Obs
d(y Obs , ỹObs ) = d({GLC Obs k } N k=1 , {GLC Obs k ( θ, η, ν)} N k=1 ) = N k=1 |(GLC Obs k -GLC Obs k-1 ) -(GLC Obs k ( θ, η, ν) -GLC Obs k-1 ( θ, η, ν)|
which corresponds to the empirical Wasserstein-1 distance [START_REF] Kantorovich | The mathematical method of production planning and organization[END_REF] in the case of microdata10 . Explored in the context of ABC by [START_REF] Bernton | Approximate Bayesian computation with the Wasserstein distance[END_REF], this distance summarizes the absolute discrepancies between all order statistics across observed and simulated data |y Obs (i) -ỹObs (i) | , i = 1, ..., N .

In approximating L(y Obs |θ, η, ν) for an ABC implementation of the (AM ) algorithm, parameter values ( θ, η, ν) ∈ Θ yielding simulated data ỹObs resembling y Obs more closely under d(., .) than others should be given larger importance. This is commonly introduced exploiting a kernel function K τ giving increasingly larger weight to parameter values with a lower discrepancy ε ≡ d(y Obs , ỹObs ). A common 'smooth' kernel for this purpose is the Gaussian kernel (e.g., see Ratmann 2010):

K gauss τ (ε) = 1 τ × 1 √ 2π × exp - 1 2 ε τ 2 , ε ≡ d(y Obs , ỹObs )
. Under this kernel the ABC discrepancies are weighted following a Normal distribution centered at zero (i.e., highest weight is given to values ( θ, η, ν) ∈ Θ exactly reproducing y Obs ), and with a standard deviation of τ .

The bandwidth parameter τ rules the strictness of the ABC approximation to L(y Obs |θ, η, ν) by defining how the weights K gauss τ (ε) decrease with an increase in the discrepancy ε. The approximation is exact when τ → 0, as only parameter values which exactly reproduce the observed income distribution are given a non-zero weight, and τ → ∞ amounts to considering any and all parameter values in Θ equally likely to have generated the observed data (i.e., the likelihood is approximated as a flat function).

Following [START_REF] Silva | Parametric estimation of income distributions using grouped data: an Approximate Bayesian computation approach[END_REF], an ABC (AM ) algorithm with these settings can be devised extending the [START_REF] Marjoram | Markov chain monte carlo without likelihoods[END_REF] ABC implementation of the (MH ) algorithm. Algorithm (2) below presents a possible implementation, denoted (ABC-AM ) in what follows.

At any j-th step, the (ABC-AM ) algorithm draws M candidate parameter values { φ(m) } M m=1 from the adaptive proposal g Σ (j-1) , simulates a single income distribution from the model for each such candidate, and computes their respective discrepancies with respect to the observed income distribution. The candidate with the lowest discrepancy is then taken as the j-th candidate φ(j) , along with its associated ABC discrepancy ε(j) , in the same spirit as [START_REF] Clarté | Componentwise approximate Bayesian computation via Gibbs-like steps[END_REF]. Finally, the MH accept-reject rule is computed with respect to the ABC approximation K τ (ε (j) ) of the likelihood.

Introducing the ABC approximation K τ (ε (j) ) to the likelihood L(y Obs | θ(j) , η(j) , ν(j) ) implies sampling from an ABC approximation to the target posterior distribution π τ (θ, η, ν|{GLC Obs k } N k=1 ). This posterior distribution might differ from that in (2) for several reasons. One first source of differences lies on the quality of the approximation. The main determinant of this is the choice for the bandwidth parameter τ . In practice, the choice for this bandwidth results from calibrating the sampling algorithm through several initial runs balancing strictness of the approximation and computational cost.

The second main source for differences between π τ (θ, η, ν|{GLC Obs k } N k=1 ) and π(θ, η, ν|y Obs ) concerns the possible loss of information due to summarizing the data through the GLC and not through the microdata directly. If what can be learnt about (θ, η, ν) from the data represented through the GLC is less than what can be learnt from microdata then their respective estimated posterior distributions will differ even when the ABC approximation to the likelihood is exact (i.e., when τ → 0). Together, the parametric framework for income distributions under 'missing rich' issues developed in the previous section along with the Bayesian empirical strategy Algorithm 2: An AM ABC (ABC-AM ) algorithm.

1: Initialization:

2: Set Σ (0) , J 0 , M , τ 3: Until K gauss τ (ε (0) ) > 0: 4: Sample ( θ(0) , η(0) , ν(0) ) ≡ φ(0) from p( φ) 5:
Generate {GLC Obs k ( φ(0) )} N k=1 by simulating from f Obs y (.; φ(0) )

6: Generate ε(0) = d({GLC Obs k } N k=1 , {GLC Obs k ( φ(0) )} N k=1 ) 7: Sampling: for j = 1, ..., J do 8: Sample { φ(m) } M m=1 ∼ g Σ (j-1) (φ, φ(j-1) ) from the candidate g Σ (j-1) 9: Generate {GLC Obs k ( φ(m) )} N k=1 by simulating from f Obs y (.; φ(m) ) , m = 1, ..., M 10: Generate ε(j) = min m∈{1,...,M } d({GLC Obs k } N k=1 , {GLC Obs k ( φ(m) )} N k=1 ) and candidate φ(j) = arg min φ(m) {d({GLC Obs k } N k=1 , {GLC Obs k ( φ(m) )} N k=1 )} M m=1
11: Accept and store ( φ(j) , ε(j) ) with probability:

ρ (j) = min 1, K gauss τ (ε (j) ) × p( φ(j) ) × g Σ (j-1) ( φ(j-1) , φ(j) ) K gauss τ (ε (j-1) ) × p( φ(j-1) ) × g Σ (j-1) ( φ(j) , φ(j-1) )
e.g., if u (j) ≤ ρ (j) where u (j) is a draw from a U nif orm(0, 1) distribution otherwise store ( φ(j) , ε(j) ) = ( φ(j-1) , ε(j-1) ) if j > J 0 then 12: Update Σ (j) end end presented in this section allow for a broad range of applications. The following section illustrates some of the main income distribution analysis possible under this approach.

Applications and examples

Applications on simulated data

Simulated-data applications can give insight on the performance of the ABC approach in making inference on (θ, η, ν) in a controlled setting exploiting a model following (1). Consider a hypothetical population's income distribution following a GB2 distribution y i ∼ f GB2 y (.; θ) ≡ GB2(α, β, p, q), with parameters α, p, and q ruling the shape of the distribution and β ruling the scale 11 . Typically, these parameters are the focus of the analysis of the income distribution. However, if the available data y Obs is presumably affected by any of the 'missing rich' forms considered in the previous sections, additional parameters ruling assumed parametric forms for these issues must also be introduced into the analysis.

Assume that microdata samples from this population's income distribution may be jointly affected by high-income under-reporting following an LPU scheme with parameters (p, δ) and high-income non-response following a right-truncation scheme with parameter t where t p. Noting that under this joint scheme 12 : , a model for the observable data y Obs can be obtained applying (1) 13 :

∂m -1 (y Obs i ; p, δ) ∂y Obs i = 1 1 -δ × 1(y Obs i > F -1;GB2 y (p; α, β, p, q))
f y Obs (y Obs i ; θ, p, δ, t) = f GB2 y m -1 (y Obs i ; p, δ); θ × ∂m -1 (y Obs i ;p,δ) ∂y Obs i × ϕ(m -1 (y Obs i ; p, δ); t) f GB2 y (m -1 (y Obs i ; p, δ); θ) × ∂m -1 (y Obs i ;p,δ) ∂y Obs i × ϕ(m -1 (y Obs i ; p, δ); t)dy Obs = f GB2 y y Obs i + 1(y Obs i > F -1;GB2 y (p; θ)) × δ(y Obs i -F -1;GB2 y (p;θ)) 1-δ ; θ × 1(y Obs i ≤ (1 -δ)F -1;GB2 y (t; θ) + δF -1;GB2 y (p; θ)) t × (1 -δ × 1(y Obs i > F -1;GB2 y (p; θ)))
(3) Equation ( 3) expands the GB2 distribution to allow for LPU (whenever p t and δ > 0) and for non-response in the form of a right-truncation (whenever p t < 1). For

11 The GB2 distribution GB2(α, β, p, q) has pdf:

y i ∼ f GB2 y (y i |α, β, p, q) = αy αp-1 i β αp B(p, q) 1 + yi β α p+q , (y i , β, α, p, q) ∈ R 5 +
where B(p, q) denotes the Beta function, defined as B(p, q) = 1 0 t p-1 (1 -t) q-1 dt. See [START_REF] Chotikapanich | Using the GB2 income distribution[END_REF] for other parametric expressions under this distribution such as the Gini coefficient.

12 In what follows, F -1;GB2 (.; α, β, p, q) denotes the quantile function of the GB2 distribution.

13 For simplicity of notation, θ = (α, β, p, q) in what follows.

illustratory purposes, a first experiment of interest consists in estimating the posterior distribution π(θ, p, δ, t|y Obs ) using the (ABC-AM ) algorithm through this model over a sample of N simulated incomes y Obs = {y Obs i } N i=1 . In particular, this exercise is most interesting when the simulated data is effectively affected by LPU and right-truncation forms of 'missing rich' jointly.

Benchmark parameter values can be set to α, β 1000 , p, q = (2.3, 10, 1.75, 1.25) and (p, δ, t) = (.5, .15, .99) in this interest14 . These correspond to a population income distribution with an average income of 15054 and a Gini coefficient of 0.348. Data simulated under this setting corresponds to a sample from a GB2 distribution which starts being affected by LPU above the median with a slope of δ = .15 and which contains no observations for units above the .99-th population's income distribution percentile.

Data can be simulated in this specific case by sampling N 1+t incomes from the GB2 distribution and applying the LPU and right-truncation transformations under the benchmark values. This yields a single random sample of N observed incomes. The sample used in this exercise was generated in this way, for a hypothetical population of 10000 units (i.e., N = 9900). Figure 3 below illustrates how a sample generated under this setting relates to the theoretical observed incomes' distribution f y Obs under (3) and to the respective complete population's f GB2 y income distribution.

Several conditions must be considered in elliciting a joint prior probability distribution for the model parameters p(θ, p, δ, t). Firstly, this joint prior distribution can be set as the product of several marginal prior distributions:

p(θ, p, δ, t) = p(α) × p β 1000 × p(p) × p(q) × p(p) × p(δ) × p(t)
Secondly, given the high flexibility of the GB2 distribution it is possible to represent virtually any specific case of this distribution in a constrained range of parameter values. In this sense, the marginal prior distributions for the GB2 parameters were set as follows:

α ∼ p(α) ≡ Gamma(1, 1) β 1000 ∼ p β 1000 ≡ Gamma(5, 2) p ∼ p(p) ≡ Gamma(1, 1) q ∼ p(q) ≡ Gamma(1, 1)
. This amounts to prior beliefs on the shape parameters α, p, and q following a right-skewed Gamma distribution with mode at the value 1 and to prior beliefs on β 1000 following another right-skewed Gamma distribution with mode approximately at the value 8. Thirdly, reflecting a strong prior belief on the presence 'missing rich' issues in the data, the (p, δ, t) parameters were given the following prior distributions:

p ∼ p(p) ≡ Beta(8, 5) δ ∼ p(δ) ≡ Beta(1, 5) (1 -t) ∼ p(1 -t) ≡ Beta(1, 25)
These reflect empirically-relevant values for the literature using right-truncation forms for non-response (e.g., see Jorda and Niño-Zarazúa 2019) and that exploring high-income under-reporting in survey data (e.g., see [START_REF] Flachaire | Underreporting of top incomes and inequality: A comparison of correction methods using simulations and linked survey and tax data[END_REF]. Importantly, these prior beliefs also give considerable probability to the 'complete data' scenario where no under-reporting or non-response issues affect the sample. This is, it is also made likely a priori that the observed income distribution may be correctly represented by a single GB2 distribution without introducing 'missing rich' phenomena.

Finally, several constraints may be imposed on the ellicited joint prior distribution to further constrain the parameter space. Imposing restrictions for finite variance on the GB2 income distribution amounts to giving 0 prior probability to parameter values with α < 2 q and α < -1 p . Additionally, because under-reported incomes have no relevance if they correspond to a true income above the truncation point, the restriction t > p is imposed 15 . Figure 4 below summarizes these ellicited prior distributions for each of the model's parameters. Three central scenarios are explored in applying the (ABC-AM ) algorithm. Firstly, to evidence the biases that these forms of 'missing rich' induce if the issue is not taken into consideration, a simple GB2 distribution is fit to the data. A second scenario consists of estimating the income distribution parameters (α, β, p, q) under (3) conditional on fixing the correction quantities (η, ν) = (p, δ, t) at their true values. Finally, a third scenario consists of estimating all parameters in (3) elliciting prior uncertainty in (η, ν). In all cases, the algorithm is set with parameters τ = 25, M = 10, J 0 = 15000, and Σ (0) = diag(.01, .1, .01, .01, .01, .01, .01), and J = 250000 MCMC samples are obtained, discarding the initial 50000 considered as the burn-in period where the algorithm's adaptive terms are being calibrated. For computational ease, the simulated data taken as estimating sample was summarized by its GLC computed at sample centiles {GLC Obs k } 100 k=1 .

Figure 5 below illustrates the goodness-of-fit of the resulting estimates for all three scenarios explored, along with the estimating sample {GLC Obs k } 100 k=1 . The posterior predictive distributions computed following {GLC Obs k ( φ(j) )} 100 k=1 , j = 1, ..., 200000 graphically match the estimating data very closely for all scenarios, with very narrow 95% highest posterior density intervals 16 reflecting little uncertainty around the predicted values. This close fit to the data even holds at the top 10% of the income distribution, where the 'missing rich' issues are progressively more present. An additional implication of these results is that both the sample Gini coefficient and the sample mean income 15 Formally, these restrictions impose the following joint prior distribution:

p(θ, p, δ, t) = p(α) × p β 1000 × p(p) × p(q) × p(p) × p(δ) × p(t) × 3 i=1 C (i) (θ, p, δ, t) with        C (1) (θ, p, δ, t) = 1 α > 2 q C (2) (θ, p, δ, t) = 1 α > -1 p C (3) (θ, p, δ, t) = 1 (t > p)
16 Highest posterior density intervals are computed in what follows as the narrowest interval within the estimated posterior distribution accumulating 95% of the mass of the distribution.

are accurately fit in all three scenarios, as both are quantities of the GLC. Because the GLC is determined jointly by all parameters in the model, however, this overall good fit may hide important biases or differences in the estimated posterior distributions for each parameter individually. k=1 . Estimates separately obtained without 'missing rich' corrections (GB2, in legend), conditional on the true (p, δ, t) = (.5, .15, .99) correction parameters, and with prior uncertainty on these. 95% highest posterior density interval bounds in respective dashed lines. Bottom: Same as top, focusing only on top 10% of the income distribution.

Estimated ABC marginal posterior distributions for each of the income distribution parameters in (3) are summarized in figure 6 below. As a first observation, all scenarios yield posterior distribution estimates which significantly differ from the ellicited prior distributions, effectively updating these prior beliefs. The most relevant result is the strong bias for all parameters affecting the estimates obtained without considering 'missing rich' issues. In contrast, both scenarios correcting for these issues yield estimated posterior distributions centered at their true value. Additionally, introducing uncertainty on the (p, δ, t) parameters yields posterior distributions for the income distribution parameters which reflect higher uncertainty than the respective estimates obtained under the known true values for these. Note: Kernel density estimates for ABC marginal posterior distributions over J = 150000 samples { α (j) , β (j) 1000 , p (j) , q (j) } 150000 j=1 from the (ABC-AM ) algorithm after 50000 burn-in samples. Estimates separately obtained without 'missing rich' corrections (GB2, in legend), conditional on the true (p, δ, t) = (.5, .15, .99) correction parameters, and with prior uncertainty on these. Traceplots of the underlying MCMC samples for estimates with uncertainty on (η, ν) below, with burn-in period in gray. True parameter values in dashed black lines.

For the scenario where (p, δ, t) are also to be inferred from the data, figure 7 below summarizes the estimated posterior distributions for these parameters. The estimates showcase a significant update of the ellicited prior distributions, with posterior distributions centered at the true values for these correcting quantities. This result provides support of the ABC approach as a fruitful empirical strategy for parametric inference on income distributions through data affected by 'missing rich' issues of uncertain magnitude. Note: Kernel density estimates for ABC marginal posterior distributions over J = 150000 samples { p(j) , δ (j) , t (j) } 150000 j=1 from the (ABC-AM ) algorithm after 50000 burn-in samples. Traceplots of the underlying MCMC samples for estimates with uncertainty on (η, ν) below, with burn-in period in gray. True parameter values in dashed black lines.

One additional illustration that this exercise provides concerns the impact of correcting for 'missing rich' issues in making inference on income growth and inequality at the population level. Figure 8 presents the posterior predictive distributions of the population's mean income and Gini coefficient, both determined by the (α, β, p, q) coefficients alone. These distributions evidence a significant under-estimation of both growth (through the mean income) and inequality (through the Gini coefficient) when the 'missing rich' issues are neglected. Only the scenarios estimated under corrections through (p, δ, t) achieve estimates closely reproducing the true value at the population level. Despite all three estimates closely reproducing the estimating sample GLC, the bias affecting parameter estimates when 'missing rich' issues are not considered make inference on the population's income distribution invalid.

Gini

Note: Kernel density estimates for ABC posterior predictive distributions of population's mean income (µ) and Gini coefficient over J = 150000 samples { p(j) , δ (j) , t (j) } 150000 j=1 from the (ABC-AM ) algorithm after 50000 burn-in samples. Estimates separately obtained without 'missing rich' corrections (GB2, in legend), conditional on the true (p, δ, t) = (.5, .15, .99) correction parameters, and with prior uncertainty on these. Traceplots of the underlying MCMC samples for estimates with uncertainty on (η, ν) below, with burn-in period in gray. True values in dashed black lines.

Real data applications:

As an illustratory example on real data, the European Union's Statistics on Income and Living Conditions (EU-SILC) provide an interesting household survey setting. EU-SILC data provides information on people and households within the EU representative at the country level, covering most countries in the EU yearly since 2005, and under a common framework defining the exact definitions of incomes and populations to be surveyed.

Although several calibrations are done over EU-SILC samples and sample weights for enforcing population representativeness on several dimensions, these are not done on income variables. Recent analysis have explored 'missing rich' phenomena on EU-SILC data (e.g., see Hlasny and Verme 2018, [START_REF] Bartels | An integrated approach for a top-corrected income distribution[END_REF][START_REF] Angel | What did you really earn last year?: explaining measurement error in survey income data[END_REF][START_REF] Carranza | Top income adjustments and inequality: An investigation of the EU-SILC[END_REF][START_REF] Ederer | Distributional national accounts (DINA) with household survey data: Methodology and results for European countries[END_REF], suggesting this issue to be present with different magnitudes in all countries and periods.

A key income variable in income distribution analysis on EU-SILC data is household disposable income under the OECD-modified equivalence scale17 (HX090). This considers all gross household incomes in the data net of regular taxes on wealth, regular interhousehold transfers paid, and regular taxes on income and social insurance contributions.

Although this aggregate variable includes definitions of income variables that are common to all countries covered by EU-SILC, the sources from which the data is obtained differ across cross-sectional waves of data and countries. In particular, while some countries rely entirely on survey responses to measure these income variables, other countries source these variable entirely or partially from administrative registers. These differences in sources across waves and countries introduce large heterogeneities in the quality of the data in terms of under-reporting as registers are considered more reliable than survey responses. Additionally, the rising use of register sources determines that for some countries some of the waves of data are sourced from surveys and other waves are sourced from register. This can produce trends in the observed income distributions across waves without necessarily reflecting trends of the true population's income distribution.

Information on non-response rates for each country and wave of EU-SILC is publicly available through the corresponding quality reports published by the European Commission. Household non-response rates, in particular, can be informative about the overall degree of non-response affecting an observed distribution of household incomes. These rates are computed from a country-level household reponse rate, which is the product of address contact rates (i.e., the share of households in the sampling frame that were successfully contacted) and household response rates (i.e., the share of households in the sampling frame that completed their survey after being successfully contacted).

Table 1 below summarizes EU-SILC samples for five selected countries (Austria, Germany, France, Spain, and Italy) and for the 2005, 2007, 2011, and 2016 waves. The mean and Gini coefficient for household disposable income distributions summarize the observable trends in growth and inequality across countries and waves. Because these distributions are presumable affected by 'missing rich' issues, these values serve as a lower bound estimate for the corresponding population's mean and Gini coefficient. With many heterogeneities, all countries experienced mean income growth and, with the exception of Italy, income inequality increased from 2005 to 2016. Finally, overall household non-response rates show large disparities across countries and years in terms of levels and trends, illustrating possible heterogeneities in the incidence of this issue on the respective observed income distributions. Under the same settings as in the simulated data application, model (3) can be fit to the EU-SILC samples through the (ABC-AM ) algorithm. Figure 9 below summarizes the resulting estimates in terms of goodness-of-fit to the weighted-sample mean income and Gini coefficient, as well as the estimated population values of these. As a first observation, the obtained results reproduce very accurately the levels and trends of mean incomes and Gini coefficients for all waves and countries considered.

Concerning the estimates of the population income distribution, a first observed result is these do not match their sample counterparts for any case. This suggests that 'missing rich' corrections are required in order to reproduce the observed distributions so accurately. Consequently, the estimated GB2 population income distribution parameters imply levels of income growth and inequality above those observed in the data. These population level estimates reproduce similar changes in mean incomes across waves as those in the observed samples, with some increase in the uncertainty around these quantities for the latter years in the case of France and Italy. Population estimates that reproduce the same dynamic of mean incomes as in the observed sample can be indicative of the total mass of incomes affecting the data due to 'missing rich' issues changing very little across waves.

As summarized by the Gini coefficient, the income inequality dynamics implied by these estimates at the population level pose some contrasts with their sample counterpart. In all cases inequality is estimated to be higher at the population level than what sample estimates suggest. The non-overlap of the computed credible intervals across populationand sample-level estimates of the Gini provide strong evidence of this. Another observation is that the uncertainty around the population Gini is relatively stable across EU-SILC waves for all countries, with strong heterogeneities across countries. This uncertainty is in some cases high enough to make inequality increases and decreases across time equally likely, as is the case for the Austrian Gini between the 2011 and 2016 waves. In other cases, however, the estimates provide clearer evidence of significant increases in income inequality across periods, as can be seen for France and Germany across 2005 and 2016 waves.

Conditional on the assumed forms for 'missing rich' issues in these applications, the estimated parameters can suggest margins where the representativeness of the data changes across waves. The detailed estimates in Appendix B suggest that the right-truncation parameter t introducing high-income non-response in the model ( 3) is estimated to be within the top .1% of the population income distribution in all waves and countries. This can suggest that significant non-response issues are mostly concentrated on households within this income group.

Concerning under-reporting issues, these estimates suggest strong heterogeneities in the share of the population affected by progressive under-reporting and the progressiveness of under-reporting across waves and countries. As quantified by the estimates for the p and δ parameters these are estimated to range from .5401 to .8525, and from .0779 to .3997 respectively across the selected EU-SILC samples. Taken together with the estimates for t, these results illustrate that reproducing the observed income distributions accurately under model (3) always requires jointly correcting for progressive under-reporting of incomes above the median and for right-truncation non-responses somewhere within the top .1% of the income distribution. , η (j)

, ν (j)

} 150000 j=1 from the (ABC-AM ) algorithm after 50000 burn-in samples.

Concluding remarks

Building on previous 'missing rich' correction methods explored in the literature, a new framework for parametric income distributions is proposed. This framework directly deals with these corrections based on expanding pre-existing parametric distributions with functional forms for both reporting and non-response mechanisms. As a model for data on incomes presumably affected by 'missing rich' this framework allows for devising empirical strategies to infer jointly features of the associated population's income distribution and features of the 'missing rich' issues in the data.

In dealing with the several constraints that must be faced by such an empirical strategy, the ABC approach is proposed as a fruitful method. This Bayesian estimator allows for updating prior uncertainty on the true population income distribution and the often uncertain 'missing rich' quantities affecting the observed data. This is achieved by attempting to reproduce the observed GLC with simulated incomes from the specified model.

The illustrative applications presented in this paper evidence some of the virtues of the ABC approach for inference on the parameters of GB2 income distributions under data affected by 'missing rich'. In a simulated-data setting, the analysis illustrates the several biases affecting inference on a population's income growth and inequality when 'missing rich' issues affecting the estimating data are neglected. This experiment also suggests the ABC approach to be fruitful in learning about uncertain 'missing rich' quantities given informative prior beliefs about these.

Applications to cross-sectional EU-SILC data on incomes give insight on the suitability of the framework in a typical household survey setting. The resulting estimates imply that reproducing the observed income distributions accurately requires considering both high-income under-reporting and non-response in all settings analysed. The analysis also illustrates how inference on population income distributions can be made under uncertain 'missing rich' quantities, uncovering contrasts between some of the observed trends at the sample level and those inferred for the respective population.

Further work seeking to understand the pitfalls of the proposed approach could explore further setups in terms of the specified 'missing rich' parametric forms. In principle, if these forms are capable of representing similar patterns of under-reporting and non-response, then estimates obtained under them should yield very similar results at the population level.

Another possible line of analysis for future work concerns exploring this empirical strategy in a linked-data setting. If comparing survey-sourced incomes and registersourced incomes evidences some form of progressive under-reporting and high-income non-response, then estimates obtained under this framework using the survey data alone should be found to reproduce 'missing rich' patterns consistent with these.

In the specific case of the EU-SILC, an additional direction for future work concerns integrating available external information on the representativeness of the observed income distribution into the prior beliefs for the 'missing rich' quantities. In particular, household non-response rates and information about the specific sources for the observed incomes for a given wave of data and country can be exploited in setting up informative prior probabilities. This could help in accounting for possible artificial trends arising from changes in the sampling or income sources and not from actual changes in the population's income distribution.

Finally, a possible extension of this framework involves making inference on income distributions of populations defined at aggregate levels such as regions or the globe. Taking the mixture of all countries' income distributions, estimates obtained accounting for 'missing rich' issues at the country level can be used to study patterns of income growth and distribution on aggregate levels. This result supports the use of the Wasserstein-1 distance as a common unidimensional discrepancy allowing for ABC inference either with microdata or grouped data summarized through the GLC.

To compute the Wasserstein-1 distance on grouped data in the form of K groups' GLC coordinates (i.e., observed incomes are split into K segments with bounds z (k) , k = 1, ..., K) we could have:

× K k=1 ȳ(k) × n (k) K - ȳ(k) × n (k) K k=1 ȳ(k) × n (k) × K k=1 ȳ(k) × n (k) K = 1 K K k=1 ȳ(k) × n (k) -ȳ(k) × n (k) = 1 K K k=1 (ȳ (k) -ȳ(k) ) × n (k) = 1 K K k=1 (ȳ (k) -ȳ(k) ) × n (k) = 1 K K k=1 n i=1 (y (i) -ỹ(i) ) × I(z (k) ≥ y (i) ≥ z (k-1) ) × I(z (k) ≥ ỹ(i) ≥ z (k-1) )
which, in the trivial case of having K = n groups (i.e., one observation per group) corresponds to the expression for this distance on microdata.

These results suggest that in the case of grouped data we can exploit the discrepancies between GLC curves through their first-order difference (i.e., through the approximation to the Wasserstein-1 distance).

Geometrically, the Wasserstein-1 distance computes the average absolute difference between the quantile functions of two distributions. When only grouped data is available, this average distance could be approximated by first computing the area between both empirical quantile curves within each interval of the grouped data, sum these areas across all intervals and divide by the number of intervals. The approximation comes the fact that in computing these areas the curves might cross within an interval and so we would have no way of accounting for those differences which counteract within the interval (i.e., the absolute value is applied at the interval level in the grouped data expression above). For a same population size n, however, the quality of the approximation always increases with K.

Having access to microdata allows a computationally-cheap alternative in which userspecified groups or bins can be defined for exploiting the grouped-data approximation to the Wasserstein-1 distance. For instance, instead of grouping the data on sample deciles, one could define broader groups for lower incomes and finer groups for higher incomes

  is strictly convex for income levels above F -1 y (p; θ). Generalized Pareto replacing: (e.g., Atkinson and Piketty 2007, Chapter 2,, Jenkins 2017, Hlasny and Verme 2022, Charpentier and Flachaire 2022):

Figure 1 :

 1 Figure 1: Quantile ratios under common replacing schemes

Figure

  Figure 2: Income-proportional reweighting schemes

  income share up to the k-th observation in the ordered sample and µ Obs denoting the sample average income.

  δ); α, β, p, q) × ϕ(m -1 (y Obs i ; p, δ); t) × ∂m -

Figure 3 :

 3 Figure 3: Population, theoretical, and sample densities for model (3)

Figure

  Figure 4: Prior distributions

Figure 5 :

 5 Figure 5: Posterior predictive estimates of the GLC

Figure 6 :

 6 Figure 6: Estimated ABC marginal posterior distributions for income distribution parameters

Figure 7 :

 7 Figure 7: Estimated ABC marginal posterior distributions for 'missing rich' parameters

Figure 8 :

 8 Figure 8: Posterior predictive estimates of population mean income and Gini coefficient

Figure 9 :

 9 Figure 9: Posterior predictive estimates of mean income and Gini coefficients -EU-SILC countries

  ρ( fy , fy (.; θ))= K k=1 |(GLC(y (k) ) -GLC(y (k-1) )) -(GLC(ỹ (k) ) -GLC(ỹ (k-1)

Table 1 :

 1 EU-SILC sample descriptives for selected countries Note: EU-SILC samples for Austria (AT), Germany (DE), France (FR), Spain (ES), and Italy (IT) from2005, 2007, 2011, and 2016 waves. Only considers households with reported household disposable income (HX090) of at least 1 euro. Household non-reponse rates as reported in the publicly-available quality reports for each wave. Weighted-sample estimates of mean incomes (µ Obs ) and Gini coefficients.

	Country	Wave	N	Household non-response rate	µ Obs	Gini
		2005	5146	0.38	20212.24 0.27
	Austria (AT)	2007 2011	6805 6182	0.22 0.23	20405.17 0.28 23948.16 0.29
		2016	5992	0.27	26274.72 0.28
		2005 13078	0.35	18078.73 0.27
	Germany (DE)	2007 14047 2011 13473	-0.21	20084.84 0.31 21047.33 0.30
		2016 13260	0.23	23424.24 0.31
		2005	9745	0.16	18237.49 0.29
	France (FR)	2007 10485 2011 11348	0.14 0.18	18423.25 0.27 23934.22 0.31
		2016 11446	0.17	25788.05 0.30
		2005 12865	0.28	12289.05 0.33
	Spain (ES)	2007 12234 2011 12993	0.23 0.22	13520.56 0.32 16535.78 0.33
		2016 14168	0.20	16151.14 0.34
		2005 21874	0.15	16648.63 0.33
	Italy (IT)	2007 20809 2011 19234	0.14 0.25	17422.81 0.32 18491.99 0.32
		2016 20966	0.21	18839.71 0.32

Source: Own calculations from EU-SILC.

For simplicity and without loss of generality, only forms of the type m(.; η) ≡ m(y i ; η) and ϕ(.; η) ≡ ϕ(y i ; η) are considered in what follows.

For example, in their study comparing household survey incomes to linked tax return data for Uruguay[START_REF] Flachaire | Underreporting of top incomes and inequality: A comparison of correction methods using simulations and linked survey and tax data[END_REF] find evidence of progressive under-reporting potentially affecting the survey data above p = .50. In studying similar linked data for the Austrian case,[START_REF] Angel | What did you really earn last year?: explaining measurement error in survey income data[END_REF] find evidence of progressive under-reporting of wages potentially affecting their survey above the p = .50 percentile. The degree of progresiveness of under-reporting can be quantified in terms of δ under a linear approximation to the observed under-reporting patterns.

The addition of the diagonal matrix χ × I d is needed with an insignificantly small but non-zero χ to assure the non-singularity of Σ(j-1) and assure the convergence of the MCMC sampling distribution to π(θ, η, ν|y Obs ).

For simplicity, it is assumed in what follows that simulated data are in the form of independent microdata ỹ = {ỹ Obs i } N i=1 .

For a grouped-data formulation see[START_REF] Silva | Parametric estimation of income distributions using grouped data: an Approximate Bayesian computation approach[END_REF].

See the derivations in Appendix A for a grouped-data implementation of this discrepancy.

The benchmark value for β being 10000, it is here scaled by 1000 in the interest of numerical stability when applying (ABC-AM ).

The OECD-modified equivalence scale computes a household's size HX050 as: HX050 = 1+0.5×(#household members aged 14 and over -1)+0.3×(#household members aged 13 or less)

The project leading to this publication has received funding from the French government under the "France 2030" investment plan managed by the French National Research Agency (reference: ANR-17-EURE-0020) and from Excellence Initiative of Aix-Marseille University -A*MIDEX. Declaration of interest: None. All R replication code for the results shown in this paper are available upon request and will soon be made available on the author's Github site.

Appendix A

Comparing income distributions through the Wasserstein distance.

When microdata on a sample of incomes {y (i) } n i=1 is available we can estimate empirical quantiles of the distribution and base the inference on a parametric model's parameters for the population income distribution around fitting the model at these quantiles. This is so because, quantiles are informative on both shape and scale of the distribution. This allows for an ABC approach without the need of summarizing the data through a small set of summary statistics.

The Wasserstein distance, originally developed in [START_REF] Kantorovich | The mathematical method of production planning and organization[END_REF], was recently explored for the purpose of summary-free ABC inference in [START_REF] Bernton | Approximate Bayesian computation with the Wasserstein distance[END_REF] and [START_REF] Drovandi | A comparison of likelihood-free methods with and without summary statistics[END_REF].

The Wasserstein distance between an income distribution f y with quantile function F -1 y and a parametric distribution model for this distribution f y (.; θ) with quantile function F -1 y (.; θ) follows:

. In the case of p = 1 and d(x, y) = |x -y| this can be consistently estimated from the sample of incomes {y (i) } n i=1 with empirical distribution fy and a simulated sample of equal size from the model with empirical distribution fy (.; θ) as (e.g., [START_REF] Berthet | A central limit theorem for Wasserstein type distances between two distinct univariate distributions[END_REF]:

This latter formulation can be interpreted as a metric comparing all sample order statistics (i.e., all sample quantiles). In essence, this metric estimates the average absolute difference between quantiles of the two distributions.

A metric ρ( fy , fy (.; θ)) may be specified under a similar logic comparing the first-order differences of all empirical GLC coordinates (i.e., estimates of the quantiles by definition of the GLC) between these microdata samples instead of order statistics directly then this amounts simply to taking the Wasserstein distance as defined above:

allowing for a particularly stricter fit on the upper tail of the distribution.

Appendix B

EU-SILC application.