Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
51,774 characters · 16 sections · 59 citation commands
Tuning Parameter-Free Nonparametric Density Estimation from Tabulated Summary Data
Researchers can often access a tabulated summary of data more easily than the original data containing confidential information at the individual level. Examples include administrative tax data containing detailed information about individual-level records of income. Tax authorities often release summary statistics of the income distributions in a tabulated format, such as the number of taxpayers and their average income, grouped by bins of income levels. Despite their lack of details, such tabulated summary data are still useful for researchers in analyzing income distributions, especially for historically old times for which micro data are no longer available.
A typical econometric method for estimating the cross-sectional distribution of a continuous random variable, such as the kernel density estimator and the empirical cumulative distribution function, requires individual-level information, and thus is not suitable when researchers only have access to tabulated summaries. In this paper, we propose a novel method for nonparametrically estimating the probability density function of an absolutely continuous distribution from tabulated summary data based on the maximum entropy (ME) principle. This ME density estimator enjoys a number of desirable properties. First, it is piecewise exponential, which allows for analytical post-estimation integration to calculate cumulative distribution estimates such as the mean, variance, top income shares, Lorenz curve, and the Gini coefficient. Second, and more importantly, unlike alternative kernel-based estimators BlowerKelsall2002,Sun2014, the ME density estimator is free from tuning parameters such as the bandwidth. This feature is attractive in practice and provides a complete theoretical justification of the method, unlike the existing kernel-based methods for which the theory fails to formally account for the tuning parameter choice in practice. We establish the strong uniform consistency for the ME density and cumulative distribution estimators under the asymptotic framework in which the resolution of bins in the table becomes finer at certain rates as the underlying sample size increases.
To illustrate our proposed method, we consider two applications. The first is a simulation study of the calculation of top income shares. Compared to the existing methods, our proposed method generates smaller bias and root mean squared error (RMSE). The second is the estimation of the income distribution and top income shares from the tabulated summaries of the U.S. tax returns data. Unlike the popular Pareto interpolation method of Piketty2003, which relies on parametric assumptions in the upper tail, our method allows for the estimation in mid-sample as well as in the tails without parametric assumptions.
\paragraph{Related literature} Our paper is related to a large literature in economics on inequality measures as well as in statistics and econometrics of estimation with grouped data.
Applied researchers working with income inequality measures have long been studying the interpolation problem from grouped data; see for instance CowellMehta1982 for an early review. The parametric methods for estimating the Lorenz curve (and hence computing the Gini coefficient) of KakwaniPodder1976 and VillasenorArnold1989 have been used by the World Bank. FeenbergPoterba1993 and Piketty2003 interpolate the upper tail of the income distribution by the Pareto distribution with local Pareto exponents estimated from a tabulated summary. See PikettySaez2003 for an application to the U.S. income distribution. More recently, BlanchetFournierPiketty2022 apply spline interpolation to what they call the inverted Pareto coefficients. Compared to this applied literature, our approach enjoys several advantages such as that
In terms of estimating the income distribution with tabulated data, our proposed method differs from the existing methods, which can be categorized into parametric and nonparametric ones. For the former, Hajargasht2012 assume a parametric income distribution and develop a generalized method of moment (GMM) estimator for the unknown parameters. See also Chen2018 and HajargashtGriffiths2020 for other GMM methods relying on parametric assumptions. Despite their good fit in some data sets Jorda2021, parametric methods in general suffer from model misspecifications. Specifically, although the estimator of the coefficients in the parametric model might still converge to some pseudo true values, the implied density estimator is in general inconsistent. The theoretical effect of misspecification on further estimations of other features, say moments and top income shares, are unknown. They may exhibit large biases in finite samples, as we show by simulation studies in Section (ref). In contrast, our estimator is nonparametric and relatively robust to parametric assumptions.
For nonparametric estimation with binned/grouped data, ScottSheather1985 study the kernel density estimator when the observations are equally spaced. On the other hand, we do not require data to be equally spaced. Without the equal spacing restriction, BlowerKelsall2002 propose an alternative estimator with Gaussian kernel function, which Sun2014 extends by allowing for other kernel functions. These kernel approaches require a bandwidth as a tuning parameter, whose choice is challenging in practice especially under the nonstandard sampling setup of binned/grouped/tabulated data. In contrast, our proposed ME estimator is free from tuning parameters and therefore more attractive in practice. Furthermore, this tuning-parameter-free feature of our proposed method provides a complete theoretical justification under weak regularity condition. Reyes2016 derive the orders of magnitude of the bias and the variance of ScottSheather1985's estimator under very strong conditions:
The first condition implies that the group structure is negligible and hence the bias and the variance are asymptotically the same as in the case with individual observations. The second condition is strong and rules out some candidate distributions, such as double Pareto, which is used in empirical studies of income distribution. Unlike the existing kernel-based methods, our ME estimator only requires the underlying density to be Lipschitz continuous and is based on minimizing the Kullback-Leibler divergence, which lead to both numerical and theoretical advantages. See Section (ref) for details.
Finally, our paper is related to the large literature that applies the maximum entropy (ME) principle. Historically, the ME principle was developed in physics to infer the population distribution (\eg, energy distribution of gas molecules) from macroscopic variables (\eg, temperature); see jaynes1957a. In economics, applications of the ME method include general equilibrium theory foley1994,Toda2010ET,Toda2015ET, diagnosis of asset pricing models stutzer1995, derivative pricing stutzer1996, discretization of probability distributions and stochastic processes TanakaToda2013EL,TanakaToda2015SINUM,FarmerToda2017QE, among others. In this paper, we apply the ME principle as a tool to impose moment restrictions implied by the tabulated summary data. Applications in econometrics include kitamura-stutzer1997 and Wu2003, among others. To our knowledge, none of the existing ME methods work for tabulated summary data.
\paragraph{Organization of the paper} Section (ref) introduces the general data framework and previews the U.S. tax return data set used in the application. Section (ref) introduces the nonparametric density estimator and proves its strong uniform consistency. Section (ref) presents simulation studies. Section (ref) applies the proposed method to the U.S. tax return data set.
Consider the latent sample $\set{Y_i}_{i=1}^n$, which is not directly observed by the researcher. Denote the order statistics in descending order by $\set{Y_{(i)}}_{i=1}^n$ such that
For each $m \in \set{1,\dots, n}$, define the partial sum of top $m$ order statistics
Consider a positive number $K$ of bins, where $\set{t_k}_{k=0}^K$ denotes the sequence of bin threshold values such that
Let $n_k$ be the number of order statistics included in the top $k$ bins, that is, $n_k\coloneqq \#\set{i:Y_{(i)} \ge t_k}$. The tabulated data are summarized as $\set{(t_k,n_k,S_{n_k})}_{k=1}^K$, which is observed by the researcher.
As a concrete example of the data framework just described, consider the latent sample $\set{Y_i}_{i=1}^n$ of the values $Y_i$ of income, where $i$ indexes potential taxpayers and $n$ is the sample size. Due to confidentiality concerns, in general there is no public access to administrative data of income. Publicly available data on the income distribution released from tax authorities often take the form of the tabulated summary $\set{(t_k,n_k,S_{n_k})}_{k=1}^K$.
Table (ref) presents an example data set from the 2019 U.S.\ tax returns.\footnote{Table (ref) shows partial information from Internal Revenue Service, Statistics of Income (SOI) Individual Income Tax Returns Publication 1304 (https://www.irs.gov/statistics/soi-tax-stats-individual-income-tax-returns-complete-report-publication-1304), Table 1.4 under “Basic Tables”. Adjusted gross income (AGI) is AGI less deficit. We omit the row corresponding to negative income.} In this example, the number of income groups is $K=18$. Column (1) shows the lower threshold $t_k$ of adjusted gross income (AGI) for each income group $k \in \set{1,\dots,K}$. Column (2) shows the number of taxpayers within each income group $k$, which corresponds to $n_k-n_{k-1}$ in our notations, where we set $n_0=0$ by convention. Column (3) shows the total income (AGI) accruing to taxpayers in each income group $k$ in units of 1,000 U.S.\ dollars, which corresponds to $(S_{n_k}-S_{n_{k-1}})/1{,}000$ in our notations, where we set $S_0=0$ by convention. We thus observe the tabulated summary data $\set{(t_k,n_k,S_{n_k})}_{k=1}^K$ of AGI, where $K=18$.
This section investigates a method to characterize a well-behaved density function of the distribution of $Y$ from the tabulated summary data $\set{(t_k,n_k,S_{n_k})}_{k=1}^K$ introduced in Section (ref). We first consider the case when the sample size is infinite and there are no sampling errors in bin probabilities and conditional means. We next propose a feasible estimator and study its asymptotic properties as the sample size tends to infinity.
Let $F$ denote the true cumulative distribution function (CDF) of $Y$, which is assumed to be absolutely continuous with probability density function denoted by $f = F'$. Suppose that the thresholds satisfy
and let $I_k\coloneqq [t_k,t_{k-1})$ denote the interval for the top $k$-th bin with top fractile denoted by $p_k=\operatorname{P}(Y\ge t_k)=1-F(t_k)$. For each $k$, the bin probability and conditional mean are defined by
where $\mathbbm1_{I}(\cdot)$ denotes the indicator function indicating that the argument belongs to the set $I$.
Obviously, given only the finite tabulation $\set{(p_k,y_k)}_{k=1}^K$, we do not have sufficient moment restrictions to pin down the true density function $f$. The maximum entropy (ME) method is useful when only certain moment conditions are given. In our context of characterizing the distribution of $Y$ from a tabulation, we can proceed as follows.
Letting $g$ denote a generic density, the given moment conditions consistent with (ref) are
for each $k \in \set{1,\dots,K}$. The ME density $f^*$ is defined by the density $g$ on $I\coloneqq [\ubar{t},\infty)$ that minimizes the Kullback-Leibler divergence (with respect to the improper uniform density)
subject to the moment restrictions (ref). Below, we let $L_+^1(I)$ denote the equivalence class (identified by the $L^1$ norm) of nonnegative, measurable, and integrable functions $g:I\to \R$. The following proposition characterizes the solution to the ME problem.
This section proposes a feasible analog of the ME density characterized in Section (ref) for estimation of the true density function $f$ of $Y$. We then establish its strong uniform consistency.
We construct a feasible ME density estimator by replacing $q_k$ and $y_k$ with their empirical analogs $\widehat{q}_k=(n_k-n_{k-1})/n$ and $\widehat{y}_k=(S_{n_k}-S_{n_{k-1}})/(n_k-n_{k-1})$, respectively. Letting $t \coloneqq \set{t_k}^K_{k=1}$ denote the vector of thresholds, we thus define the sample-analog ME estimator $\widehat{f}$ of $f$ as the solution to the constrained optimization prblem of minimizing (ref) subject to (ref) with $\widehat{Q}_k$ and $\widehat y_k$ in place of $q_k$ and $y_k$, respectively.
We now establish the almost sure uniform consistency of $\widehat{f}$ for the true density function $f$ over any compact subset of the domain of $f$ as $n\to\infty$. To this end, consider the following conditions.
Condition (ref) assumes a random sample. Condition (ref) implies that $f$ is almost everywhere differentiable with a bounded derivative, which is much weaker than typical assumptions on kernel estimators that require high-order smoothness conditions. Condition (ref) requires that the length of any bin is neither too large nor too small, as well as that the interval $[\limsup_{n \to \infty}t_K, \liminf_{n \to \infty}t_1]$ covers the domain $D$ of $f$. With these conditions, the following theorem establishes the strong uniform consistency of the maximum entropy density estimator $\widehat{f}$ for the true density function $f$ over any compact subset of the domain and also bounds the convergence rate. The proof is non-trivial and deferred to Appendix (ref).
A few remarks are in order regarding this result on the convergence rate. First, while the rate is decomposed into the the deterministic part and the stochastic part, we do not have a control over the trade-off between these two components in the absence of a tuning parameter. This implies a drawback of the tuning parameter-free approach. Second, the rate depends on the parameters $r_1$ and $r_2$ of bin lengths. Slowly vanishing bin lengths (i.e., small $r_1$ and $r_2$) yield small variances at the expense of large biases. Quickly vanishing bin lengths (i.e., large $r_1$ and $r_2$) yield small biases at the expense of large variances. Third, suppose $r_1=r_2$ for simplicity. Then, $1/6 < r_1=r_2 < 1/4$ implies that the stochastic part dominates, while $r_1=r_2 < 1/6$ implies that the deterministic part dominates. In the latter case, the limit distribution has a biased center, and it is difficult to conduct statistical inference in general. This is another limitation of the tuning parameter-free approach.
Reyes2016 derive the convergence rate of the kernel estimator proposed by ScottSheather1985. Reyes2016 assume that the bin lengths are $o(h_n^2)$, where $h_n$ denotes the bandwidth satisfying $h_n\to 0$ and $nh_n \to\infty$ as $n\to\infty$. This condition implies that the group/bin structure is asymptotically negligible, and the resulting orders, $O(h_n^2)$ and $O_p(1/\sqrt{nh_n})$, of the deterministic and stochastic parts, respectively, are the same as those in the standard case with individual observations. On the one hand, choosing a certain bandwidth $h_n$ could lead to a smaller bias or variance than our ME estimator. On the other hand, the assumption that the bin lengths are $o(h_n^2)$ is very restrictive and could be violated in empirical studies where the bins are not too small. Furthermore, Reyes2016 assume that the underlying distribution function is seven-times differentiable with bounded derivatives. In contrast, our ME estimator only requires Lipschitz continuity for $f$, which is another advantage.
Having established the strong uniform consistency of the density, it is straightforward to establish the same for the cumulative distribution function (CDF) and quantiles. Define the estimator $\widehat{F}$ of $F$ by
The following corollary shows the strong uniform consistency of $\widehat{F}$.
Let $Q_\tau = \inf\set{ y : \tau \le F(y)}$ denote the $\tau$-th quantile of $F$. Given $\widehat{F}$, we can estimate $Q_\tau$ by the analog $\widehat{Q}_\tau = \inf\set{ y : \tau \le \widehat{F}(y)}$. The following corollary shows that this quantile estimator is also consistent.
In an early review of interpolation methods from grouped data of income, CowellMehta1982 list the following ten desirable properties (with slight rewording) that the hypothetical interpolated distribution should possess.
In addition to these properties, we would like to add:
The methods reviewed in CowellMehta1982 as well as those proposed thereafter satisfy only a few of these properties. For instance, the recent method of BlanchetFournierPiketty2022 does not satisfy (ref) and (ref) (because it uses polynomial interpolation), (ref) (because it interpolates the inverted Pareto coefficients, not the density), or (ref) (they do not provide formal theorems).
In contrast, our ME density estimator $\widehat{f}$ satisfies all properties except (ref) and (ref). To see this, property (ref) holds by construction and (ref) is established in Theorem (ref). All other properties (except (ref) and (ref)) hold because $\widehat{f}$ is piecewise exponential explicitly given by (ref) and hence is nonnegative, continuously differentiable, and monotonic on each interval. Regarding properties (ref) and (ref), they clearly hold except at bin thresholds.
To illustrate property (ref) further, we present some integral formulas that are useful when computing the CDF and top income shares when the density is piecewise exponential. Consider the piecewise exponential density (ref). To simplify the notation, let $q_k=q$, $\lambda_k^*=\lambda$, $t_k=a$, and $t_{k-1}=b$. Therefore, for $y\in [a,b)$, the density is
The counter CDF (tail probability) can be computed using
Applying integration by parts, the tail expectation can be computed using
Putting all the pieces together, we obtain the following closed-form expressions for the CDF and tail expectation. (We assume $\lambda_k\neq 0$ for simplicity, and we use the notation $y\vee t=\max\set{y,t}$.)
We conduct three simulation studies that examine the performance of our proposed ME method relative to existing methods.
We first present density estimates from one large simulation draw with the sample size comparable to those in our empirical data. We consider four typical models for income distribution, namely lognormal, gamma, Weibull, and double Pareto. In each case, we choose the scale parameter so that the population mean normalizes to 1. Table (ref) summarizes these distributions.
The simulation design is as follows. For each model, we generate a random sample $\set{Y_i}_{i=1}^n$ with size $n=10^7$. Such a large $n$ is coherent with the number of tax payers in our empirical data set; see Table (ref). We set the top fractiles to
(so $K=26$) and define the $k$-th threshold $t_k$ as the top $p_k$-th quantile of $\set{Y_i}_{i=1}^n$. We then estimate the ME density $\widehat{f}$ as in Section (ref).
Figure (ref) shows the population and estimated densities for each model. Because the population densities are skewed, for visibility we plot the density of $\log Y$. In each case, the two densities $f$ and $\widehat{f}$ are nearly identical.
Next, we estimate the top income shares, which correspond to the Lorenz curve flipped along the 45 degree line. There are many existing methods for estimating the Lorenz curve as discussed in Section (ref). We implement those proposed by KakwaniPodder1976 (henceforth KP) and VillasenorArnold1989 (henceforth VA), both of which have been used by the World Bank. In addition, we implement a more recently developed method by Hajargasht2012 (henceforth HGBRC), which has been further extended by Chen2018 and HajargashtGriffiths2020. The KP and VA methods impose some parametric assumptions on the Lorenz curve and essentially run linear regressions of the group mean ($\widehat{y}_k$ in our notation) on some transformation of the proportion of each group ($\widehat{q}_k$ in our notation). The HGBRC method imposes some parametric assumptions on the underlying density and constructs a generalized method of moments estimation. Regarding KP, we implement their Method III as described in their Section 4. Regarding VA, we implement their method with $a=1$ and $d=0$ as described in their Section 4. Regarding HGBRC, we adopt their assumption of the generalized beta distribution of the second kind (GB2, McDonald1984) and the diagonal weighting matrix as proposed by Chotikapanich2007.
Our data generating process is as follows. We suppose that the population distribution is double Pareto with parameters $\alpha=2.3$, $\beta=1.1$, and $M=1$ (normalization). There are two reasons for using the double Pareto distribution with these parameters. First, this distribution has been shown to fit the income distribution very well; see for instance Toda2012JEBO. Second, unlike other parametric distributions used in Figure (ref), the double Pareto distribution admits a closed-form CDF and Lorenz curve as discussed in Appendix (ref), which is convenient for numerical evaluation. Appendix (ref) considers other distributions.
We treat the population top $0.1, 1, 5,10,\dots,95,100$ percentiles as the observed thresholds ($K=22$) and compute the population top income shares. Next, we generate random samples with sizes $n=10^4, 10^5, 10^6$ from the population distribution\footnote{Since the logarithm of a double Pareto random variable is Laplace, which is double exponential, we can generate a double Pareto random variable using $\log (Y/M)=\frac{1}{\alpha}X_1-\frac{1}{\beta}X_2$, where $X_1,X_2$ are independent exponential random variables with parameter 1. Therefore $Y=MU_1^{-1/\alpha}U_2^{1/\beta}$, where $U_1,U_2$ are independent uniform random variables on $[0,1]$.} and record the proportions of observations and their average incomes within each group, which we treat as our data.
Implementing our proposed ME, the KP, the VA, and the HGBRC methods, we report their relative bias and relative root mean squared error (RMSE) for the income shares of the top $p_0$ fractile with $p_0 \in \set{0.001, 0.01, 0.05, 0.1, 0.2, \dots, 0.9}$. More specifically, let $s_0$ denote the true top income share and $\widehat{s}_m$ the estimator in the $m$-th simulation draw with $m\in \set{1,\dots,M}$. We define the relative bias and RMSE by
respectively. Table (ref) presents the results based on $M=1{,}000$ simulations.
The findings can be summarized as follows. First, our proposed ME method performs very well in terms of both bias and RMSE. They decrease as $n$ increases and are smaller than those of the other three methods for most of the $(n,p_0)$ combinations, especially the bias. Second, the KP and the VA methods both impose some parametric assumptions on the Lorenz curve and hence implicitly on the underlying density. In particular, the VA method imposes that the Lorenz curve is a part of an ellipse. This assumption implies that the underlying density $f(y)$ (after a location- and scale-transformation) is proportional to $(1+y^2/2)^{-3/2}$, which is the Student $t$ distribution with two degrees of freedom VillasenorArnold1989. The KP method introduces a new coordinate system and imposes another parametric form on the Lorenz curve. The implied density is still parametric but does not have a closed-form expression. The HGBRC method assumes the GB2 density, which has Pareo upper and lower tails. Therefore, its performance is substantially better than those of KP and VA. In summary, these and any other parametric assumptions could lead to large bias and RMSE caused by misspecification, which do not decrease with $n$.
Finally, we compare our proposed ME estimator with the nonparametric kernel estimator proposed by BlowerKelsall2002. Given a bandwidth $h$, define $K_h(u)$ as the PDF of the normal distribution with mean zero and variance $h^2$, that is,
We implement BlowerKelsall2002 by constructing the density estimator
where $\widehat{f}_{0}(s)$ is the histogram estimator
Using the fact that $\int_{I_{k}}\widehat{f}_{0}(s) \diff s=\widehat{q}_{k}$ and $K_{h}$ is the normal density, we can simplify $\widehat{f}_\mathrm{BK}$ as follows:
where $\Phi$ denotes the CDF of the standard normal distribution.
BlowerKelsall2002 do not derive any asymptotic properties of this estimator nor theoretical requirements on the choice of the bandwidth. We implement a variety of choices of $h$ to examine its finite sample performance. Specifically, we use the rule-of-thumb choice $h=c\widehat{\sigma}n^{-1/5}$, where $\widehat{\sigma}$ is the sample standard deviation based on individual observations (which is in principle infeasible given the tabulated data) and $c\in\set{0.1, 0.5, 1.0, 1.5}$ is some constant.
Figure (ref) presents the relative RMSE for the density estimators when the data generating process is double Pareto, lognormal, and gamma as in Section (ref) with sample size $n=10^4,10^5,10^6$. Although these figures are not necessarily easy to read, the RMSEs for the ME (BK) estimator are indicated with solid (dashed) lines. As is clear from this figure, the RMSEs for the ME estimator is generally closest to the horizontal axis uniformly across quantiles, so the performance of our proposed ME method is outstanding. In addition, we also implement the kernel estimator studied by Reyes2016. Its performance is substantially worse than that proposed by BlowerKelsall2002 and hence not reported.
We consider two empirical applications of our method. First, we estimate the distribution of U.S. income distribution for particular years. Second, we estimate the top income shares (including mid-sample) over the past century.
We estimate the distribution of U.S. adjusted gross income (AGI) in 1946 and 2019. We choose 2019 because it is the most recent year for which data is available. Before World War II, because only a small fraction of the population filed for taxes, the tax returns data is not representative for the population.\footnote{The fraction of tax filers among potential tax units has been stable at around 80--90% postwar but in the range of 1--20% before 1940; see the discussion in PikettySaez2003.} For this reason, we choose 1946 because it is one of the earliest years for which the tax returns data is representative for the population. Note that unlike in recent years, the tabulated summary data set is almost the only publicly available income data set in early years such as 1946.
To make the results comparable across years, we measure income in 2019 dollars by adjusting with the Consumer Price Index (CPI). The number of income groups is $K=48$ for 1946 and $K=18$ for 2019. Figure (ref) shows the ME density estimates $\widehat{f}(\e^x)\e^x$ of log income $x=\log y$ in a semi-log scale.
We can summarize the findings as follows. First, for each year the log income distribution is bell-shaped but slightly asymmetric. Second, observe that although $\widehat{f}$ is piecewise exponential, it is not necessarily continuous at the bin thresholds as can be seen from the spikes in the 1946 density. Third, the 2019 density is more spread-out than 1946, which suggests that income inequality has increased. Finally, Figure (ref) shows the tail probability $1-\widehat{F}(y)$ in a log-log scale, which is continuous. The fact that the 2019 tail probability is higher than 1946 implies that the 2019 (real) income distribution first-order stochastically dominates the 1946 one, possibly due to economic growth. The graphs also show a straight-line pattern for high incomes, which is consistent with a Pareto upper tail documented elsewhere; see for instance deVriesToda2022RIW and the references therein. Because the slope is steeper for 1946 than in 2019, the income Pareto exponent is smaller (top income inequality is higher) in 2019.
Because the ME density $\widehat{f}$ is piecewise exponential, which is analytically tractable, it is straightforward to compute statistics such as top income shares; see Section (ref). Figure (ref) shows the top income shares, both in original and log-log scales. Figure (ref) (original scale) is essentially the Lorenz curve flipped along the 45 degree line. The fact that the 2019 curve is above the 1946 one suggests that income inequality has increased. The straight-line pattern in log-log scale (Figure (ref)) is consistent with a Pareto upper tail.
Finally, we apply the proposed method to estimate the top $p$ fractile income share for various values of $p\in [0,1]$. To construct the top income shares, we use the following approach. First, we collect the tabulated summaries of income similar to Table (ref) for each year from the IRS Statistics of Income.\footnote{See the appendix in LeeSasakiTodaWangExponents for specific details.} These tables contain information on the number of tax units (an individual or a married couple with dependents if any) and their total income within each income group. As these tables contain only tax filers, we complement them with the total number of potential tax units and total income estimated by PikettySaez2003.\footnote{We obtain the total number of tax units from the spreadsheet https://eml.berkeley.edu/ saez/TabFig2018.xls, Table A0, Column B, and total income from Column I.} We suppose that non-filers are low income households and thus do not affect the calculation of the top $p$ fractile income share if $p$ is small enough. Because the fraction of tax filers among potential tax units exceeds 0.1 (0.8) since 1936 (1945), we construct the top $p$ fractile income share for $p\in \set{0.01,0.05,0.1}$ since 1936 and also for for $p\in \set{0.2, 0.4, 0.6, 0.8}$ since 1945. Figure (ref) shows the results.
The top 1%, 5%, 10% income shares exhibit an inverse U-shaped pattern, which is well known. To the best of our knowledge, the top income shares for mid-sample fractiles (\eg, $p=0.2, 0.4, 0.6, 0.8$) have not been reported in the previous literature. We find that the mid-sample top income shares exhibit an abrupt upward jump between 1986 and 1987. This could be due to the Tax Reform Act of 1986, which significantly altered the treatment of capital gains income.\footnote{According to PikettySaez2001WP, the fraction of capital gains income included in AGI was 100% until 1933, 70% in 1934--1937, 60% in 1938--1941, 50% in 1942--1978, 40% in 1979--1986, and 100% since 1987. Because high income earners tend to hold more financial asset (and generate more capital gains), the large exclusion of capital gains in 1934--1986 likely causes the top income shares to be biased downwards.}
We next compare our top income shares to those constructed by PikettySaez2003.\footnote{We obtain the top income shares from Table A3 in the spreadsheet https://eml.berkeley.edu/ saez/TabFig2018.xls. Piketty2001book provides the details of the method for constructing these top income shares. Since it is written in French, we describe the method in Appendix (ref) for the convenience of the readers.} Figure (ref) shows the top 1%, 5%, 10% income shares constructed in two ways. We find that post-1986, our top income shares are nearly identical to those from PikettySaez2003, even though our method is nonparametric while their method is parametric (assuming a Pareto upper tail). This is likely because the upper tail of the income distribution can be well approximated by the Pareto distribution. However, there are large discrepancies between the two series pre-1986 because we have used the raw AGI without adjusting for the excluded capital gains income discussed in Footnote (ref).
Existing standard nonparametric estimators of density and cumulative distribution functions require individual-level data. Even when individual-level information is difficult to access due to confidentiality concerns, tabulated summaries of such data are often publicly available. Administrative data of income are the leading examples. In this paper, we propose a novel method of maximum entropy density estimation from tabulated summary data and establish the strong uniform consistency of the density and cumulative distribution estimators. This method enjoys many desirable properties. First, the estimator is piecewise exponential, which is analytically tractable. Using its functional form, it is straightforward to compute statistics such as top income shares. Second, and more importantly, our estimator is free from tuning parameters unlike existing kernel-based methods, which is attractive in practice. This feature provides a complete theoretical justification that our proposed estimator works in practice, unlike the existing kernel-based estimators for which the theory does not formally account for the effects of bandwidth choice in practice.