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.
124,191 characters · 21 sections · 102 citation commands
Identification and Inference in General Bunching Designs
The bunching method utilizes discontinuities in an incentive schedule to elicit behavioral responses from individuals. It is facilitated by a choice variable, also known as a running variable, manipulated by individuals to maximize their payoffs. Governments, or other authoritative entities, establish the incentive schedule applied to these individuals, which change discretely at specific cutoffs based on the chosen running variable. Such policies create varying incentives for the selection of the running variable across these cutoffs, leading to noticeable bunching -- manifesting as either excessive or missing mass in the distribution. Initiated by saezTaxpayersBunchKink2010, who observed bunching in the U.S. taxable income distribution and used it to measure the elasticity of taxable income, the bunching method has been further grounded in subsequent works (\citeauthor*{{chettyAdjustmentCostsFirm2011}}, {chettyAdjustmentCostsFirm2011}; \citeauthor*{{klevenUsingNotchesUncover2013}}, {klevenUsingNotchesUncover2013}, etc). Since then, it has found widespread applications across various economic studies: e.g., for the study of the effects of minimum wages on employment (\citeauthor*{{cengizEffectMinimumWages2019}}, {cengizEffectMinimumWages2019}; \citeauthor*{{derenoncourtMinimumWagesRacial2020}}, {derenoncourtMinimumWagesRacial2020}), the impact of tax base on corporate tax evasion in developing countries (\citeauthor*{{bestProductionRevenueEfficiency2015}}, {bestProductionRevenueEfficiency2015}), the impact of taxation and fiscal stimulus in the U.K. housing market (\citeauthor*{{bestHousingMarketResponses2018}}, {bestHousingMarketResponses2018}), the elasticity of intertemporal substitution for U.K. households (\citeauthor*{{bestEstimatingElasticityIntertemporal2020}}, {bestEstimatingElasticityIntertemporal2020}), regulatory costs faced by publicly listed firms (\citeauthor*{{ewensRegulatoryCostsBeing2024}}, {ewensRegulatoryCostsBeing2024}), and detection of score/data manipulation around cutoffs (\citeauthor*{{foremnyGhostCitizensUsing2017}}, {foremnyGhostCitizensUsing2017}; \citeauthor*{{deeCausesConsequencesTest2019}}, {deeCausesConsequencesTest2019}; \citeauthor*{{ghanemCensoredMaximumLikelihood2020}}, {ghanemCensoredMaximumLikelihood2020}) among many others.
This paper develops a formal econometric framework and tools for the identification and inference of a structural parameter in the bunching design. Our focus in this paper is on the kink design, in which a concave payoff function exhibits kinks at specific cutoffs. Another commonly used bunching design, known as the notch design, applies to cases where the incentive schedule features discontinuities, referred to as notches. Discussions on how the notch design could fare in our framework are deferred to the Appendix.
Building on a potential outcome framework similar to that of goffTreatmentEffectsBunching2024, we lay out an analytic framework for structural models in the bunching design. We analyze two dimensions of the counterfactual, where each counterfactual policy—referred to as the pre- and post-kink policies—has an incentive schedule that varies smoothly across the entire range of a choice variable. The counterfactual choices are defined as hypothetical choices an individual would make if subjected to either of these counterfactual policies. The prekink and postkink counterfactual choices are denoted by $Y^*(0)$ and $Y^*(1)$, respectively.
Identification in the bunching design hinges on two key ingredients: the counterfactual distribution of $Y^*(0)$ and the structural model of individuals' responses to policy changes. Our identification approach requires less stringent conditions on these ingredients than in existing methods. First, we explore the use of an analytic model for the counterfactual distribution for bunching identification. As noted by \citet*{blomquistBunchingIdentificationTaxable2021}, identification in bunching designs relies on implicit shape restrictions on the counterfactual distribution. Unlike many commonly used parametric approaches, our method is not restricted to an ex-ante finite-dimensional family of distributions, but it leverages a dense family of analytic functions. Under the analyticity condition, the counterfactual density and its derivatives can be stably approximated by polynomials of increasing orders. This allows us to achieve two seemingly conflicting goals: flexibly fitting the observed distribution while uniquely extrapolating the counterfactual into the unobserved region.
The idea of leveraging the analyticity assumption on the counterfactual distribution for bunching identification was explored in the earliest version of pollinger2024kinks, which developed identification and consistent estimation in an extended isoelastic model with intensive and extensive margin responses. Our approach is distinguished from it in two notable ways. First, we address both point and partial identification in a general class of structural models. Second, following the identification step, we focus on inference of a structural parameter, which requires a different set of tools and assumptions compared to estimation.
We allow for covariates included in the structural model to account for observable factors that influence individual decisions regarding the choice variable. By contrast, a common feature of structural models in bunching design is that the choice variable is determined by a scalar unobserved heterogeneity, whose distribution is directly linked to the structural parameter. As noted in \citet*{bestEstimatingElasticityIntertemporal2020}, which examines the elasticity of intertemporal substitution in dynamic consumption models, the inclusion of additional variables or parameters can significantly increase the complexity of the analysis. We address these challenges within a unified framework for structural modeling, identification, and subsequent procedures.
Previous approaches have permitted covariates only in limited ways. When certain covariates are continuously distributed, they need to be grouped into a small number of discrete categories. Moreover, such an approach requires partitioning the sample based on these categories, which could lead to aggregation bias or loss of power. This strategy has been employed in saezTaxpayersBunchKink2010 and bestEstimatingElasticityIntertemporal2020 among many others. Our approach integrates covariates into the model without requiring such coarseness.
We propose a partial identification approach that incorporates prior information on the counterfactual distribution. This method is characterized by encapsulating all potential instances of the counterfactual distribution function between the upper and lower envelope functions that are piecewise analytic, which leads to a set of inequality restrictions akin to first-order stochastic dominance. Under such a condition, the method can provide sharp identification bounds for a scalar structural parameter. This method applies to a wider range of counterfactual distributions than the first approach and accommodates various shape restrictions discussed in the literature. For instance, in the context of the isoelastic model, our result yields the elasticity bounds that are identical to those established in blomquistBunchingIdentificationTaxable2021 and \citet*{bertanhaBetterBunchingNicer2023} under the corresponding assumptions.
The polynomial strategy, initially proposed by chettyAdjustmentCostsFirm2011, is a widespread approach to estimating the unobserved counterfactual distribution. It proceeds by fitting a flexible polynomial to the observed distribution excluding a narrow window surrounding the cutoff. The fitted polynomial provides the initial estimate for the unobserved counterfactual density inside the window.\footnote{chettyAdjustmentCostsFirm2011 requires an iterative procedure to ensure that the integral constraint holds for the estimated counterfactual density. If excess bunching is small, this iterative procedure has no major impacts.} chettyAdjustmentCostsFirm2011 used the difference between the observed mass and the estimated counterfactual mass within the window as an estimate for the excess bunching fraction. The estimated excess bunching fraction $\hat{B}$, after being normalized by the counterfactual density $\hat{f}$ at the cutoff, identifies the average policy response via the small kink approximation $ {\hat{B}}/{\hat{f}} \approx E[Y^*(0)-Y^*(1)|Y^*(0)=K]. $ The structural model then identifies the average policy response as a function of the structural parameter $\theta$, represented by a function $h$ such that $E[Y^*(0)-Y^*(1)|Y^*(0)=K] = h(\theta)$. Such a restriction can justify using minimum distance estimation for $\theta$.
We highlight some weaknesses in the traditional polynomial estimation and propose an alternative scheme, termed the generalized polynomial strategy. First, the conventional polynomial estimation works under the assumption that $Y^*(1)$ and $Y^*(0)$ have proportional densities above the cutoff. However, this assumption is often incompatible with the commonly held structural assumption that both counterfactual choices are related but have different geneses. To address this issue, we propose a valid counterfactual adjustment procedure, referred to as the counterfactual correction. Consistency of the estimation and inference based on the proposed method is theoretically established and substantiated by numerical evidence.
The choice of the counterfactual adjustment procedure may have significant implications in practice. Figure (ref)(b) displays the counterfactual estimation results based on the U.S. tax records for married filers, used in saezTaxpayersBunchKink2010. The solid red curve in the right panel represents the estimated counterfactual density using the proportional adjustment by chettyAdjustmentCostsFirm2011, while the solid blue curve shows the estimated counterfactual density based on the counterfactual correction procedure. In this application, our method hinges on the structural assumption that $Y^*(0) = (0.84)^{-\theta} Y^*(1)$, where $\theta$ represents the taxable income elasticity and $0.84$ reflects the discrete change in the net-of-tax rate at the income cutoff. It proceeds by applying this transformation to all $Y^*(1)$ observed to the right of the window,\footnote{Our counterfactual correction applies to each hypothesized value of $\theta$. To draw Figure (ref)(b), we used the value of $\theta = 0.4$, which is drawn from a preliminary analysis as a plausible value.} in order to revert $Y^*(1)$ to $Y^*(0)$, followed by fitting a polynomial while excluding the data within the window. Applying this procedure produces a noticeably distinct counterfactual density compared to the one estimated using the proportional adjustment method. Numerically, $(\hat B, \hat f)$ are estimated as $(1.29\times 10^{-2},1.19\times 10^{-4})$ for the proportional adjustment method, and $(1.60\times 10^{-2},1.07\times 10^{-4})$ for the counterfactual correction method. Based on the small kink approximation, these estimates correspond to elasticities of $\hat\theta = 0.30$ and $\hat\theta = 0.41$, respectively.\footnote{Both estimates rely on the small kink approximation for fair comparison. Our estimates for $\theta$ remain stable across various polynomial degrees ranging from $5$ to $11$, up to the second decimal point. However, the estimates based on the proportional adjustment exhibit much higher sensitivity, decreasing from $0.44$ to $0.30$ as the degree ranges from $7$ to $10$.} This discrepancy could have substantial implications for policy making. According to saezUsingElasticitiesDerive2001, these elasticities translate into the optimal top rates of $0.63$ and $0.55$, assuming a common calibration of other parameters.\footnote{See footnote (ref) for details.}
Secondly, identification via the small kink approximation, which underpins many existing approaches, only reflects the first-order term in an infinite series expansion of the bunching moment. This potentially leads to a significant bias from omitting higher-order terms. For instance, incorporating the second-order term that accounts for the slope of the counterfactual density, as shown in Figure (ref)(b), corrects the overstatement of our estimate for $\theta$ by approximately 10%. For consistent estimation and correct inference, it is necessary to consider errors from finite-order approximation.
We close this gap by increasing the approximation order and incorporating omitted higher-order terms. Each of these terms can be consistently estimated from the sample based on their closed-form expression. We propose a sieve M-estimator based on a convex risk function, which can be solved easily using off-the-shelf numerical tools. We establish statistical guarantees of our estimator contingent on selecting appropriate tuning parameters, based on the theory of sieve derivative estimation and analytic function approximation.
Turning to the bunching inference problem, we focus on testing the hypothesis $H_0 : \theta_0 = \theta$. Based on our estimation strategy, we propose a t-test for a scalar parameter. To address multiple parameters, a Wald test can be developed using various bunching moments generated by choices of the weighting function. A valid confidence set can be constructed by inverting the proposed test. Our test adopts bias-aware critical values that account for potential approximation bias, whose upper bound can be calibrated using the data. The growth rates of the sieve and approximation orders are specified to ensure the asymptotic validity of counterfactual estimation and inference.
We assess the efficacy of the proposed inferential method compared to the polynomial estimator through a series of simulation studies. The findings suggest that the proposed test demonstrates good size and power properties and addresses bias in the polynomial estimator. The relevance of our method to real-world data is demonstrated through an empirical application revisiting the U.S. tax data originally used in saezTaxpayersBunchKink2010.
The paper is organized as follows. Section (ref) introduces the econometric framework. Section (ref) presents our main results on point and partial identification. Section (ref) develops a counterfactual estimation and inference method. Sections (ref) and (ref) illustrate the effectiveness of our method through Monte Carlo experiments and the empirical application. Finally, Section (ref) concludes the paper.
Let $\|A\| = \sup_{x\in\mathbb{R}^n,x\ne 0}\|Ax\|/\|x\|$ denote the spectral norm of a matrix $A \in \mathbb{R}^{m\times n}$. We denote by $A^-$ the pseudo-inverse of a square matrix $A \in \mathbb{R}^{m\times m}$. For a function $f:U\to\mathbb{R}$, where $U \subseteq \mathbb{R}^n$ is open, its mixed partial of order $\alpha = (\alpha_1, \alpha_2,\cdots, \alpha_n) \in \mathbb{Z}_+^n$ is denoted by $D_{x_1}^ {\alpha_1}D_{x_2}^{\alpha_2}\cdots D_{x_n}^ {\alpha_n} f = \frac{\partial^{|\alpha|} f}{\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\cdots \partial x_n ^{\alpha_n}}$ for $|\alpha| = \sum_{j=1}^n \alpha_j$. Here, $\mathbb{Z}_+ = \{ 0,1,2,\ldots \}$.
Denote $(\Omega, \mathcal{F})$ as the underlying measurable space, and $P$ as a probability measure on it. We denote $E_P[.]$ as the expectation computed with respect to $P$. For a random variable $X: \Omega \to \mathbb{R}$, let $F^{(P)}_X(x) = P(X \le x)$ denote the cumulative distribution function (CDF) of $X$. When $F^{(P)}_X$ is absolutely continuous, we denote $f^{(P)}_X(x) = dF^{(P)}_X(x)/dx$ as its probability density function (PDF). Likewise, we denote $F^{(P)}_{X|Y}(x, y)$ as the conditional CDF of $X$ given $Y=y$ and $f^{(P)}_{X|Y}(x , y)$ as the corresponding conditional PDF. For each $\tau \in [0,1]$, we denote $Q^{(P)}_{X|Y}(\tau, y) = \inf \{ x \in \mathbb{R}\cup\{ \pm \infty \}: F^{(P)}_{X|Y}(x | y) \ge \tau \}$ as the conditional $\tau$ quantile of $X$ given $Y = y$. We suppress $P$ in these functions when no confusion likely arises regarding the underlying measure. The indicator function, denoted by $ \mathbbm{1}\{ {.} \}$, takes the value $1$ if the event in the brackets holds and $0$ otherwise. The $(1-\alpha)$ quantile of the chi-squared distribution with $q$ degrees of freedom is denoted by $\chi^2_q(1-\alpha)$ for $\alpha \in (0,1)$.
For sequences of nonnegative real numbers $(a_n)_{n \in \mathbb{N}}$ and $(b_n)_{n \in \mathbb{N}}$, we write $a_n \lesssim b_n$ if there exists an absolute constant $C > 0$ such that $a_n \le C b_n$ for all $n \in \mathbb{N}$. For a sequence of random variables $(X_n)_{n \in \mathbb{N}}$ and a random variable $X$, we write $X_n \overset{P}{\to} X$ if $X_n$ converges to $X$ in probability, and $X_n \Rightarrow X$ if the law of $X_n$ converges weakly to that of $X$. If $(Y_n)_{n \in \mathbb{N}}$ is another sequence of random variables, we write $X_n = O_{P}(Y_n)$ or $X_n \lesssim_{P} Y_n$ if $(|X_n/Y_n|)_{n \in \mathbb{N}}$ is uniformly tight, and $X_n = o_{P}(Y_n)$ if $|X_n/Y_n| \overset{P}{\to} 0$. We adopt the convention $0/0 = 0$.
This section prepares an econometric framework for the analysis in the subsequent sections. We first present the isoelastic model as an illustrative example. We then expand on this model by formulating general conditions for the structural model to exhibit bunching under the compound policy to which individuals are actually subject. We discuss regularity conditions required for the structural model and the data, along with our treatment of optimization frictions. An extended framework for the notch design is relegated to the Appendix.
The isoelastic model refers to a model of consumption and income choices, where the individual preferences are represented by the utility function of the form
Here, $\eta_i > 0$ represents a unit-specific, unobserved preference parameter and $\theta_0 > 0$ is a structural parameter representing the taxable income elasticity to the net-of-tax rate.
In the isoelastic model, unit $i$ is assumed to select a bundle $(C_{{i}}^*, Y_{{i}}^*)$ that maximizes their utility subject to the income tax schedule specified by the government. Here, we focus on the piecewise linear tax schedule that applies to many countries including the U.S.
A piecewise linear tax system feature multiple tax brackets that are defined by values of the chosen $Y$. At the cutoffs between adjacent brackets, the marginal tax rates change discretely, while tax liabilities vary continuously. As demonstrated in bertanhaBetterBunchingNicer2023, we can focus on local analysis with only two brackets and a single cutoff. Specifically, in the global analysis, the local analysis can be applied sequentially to each pair of adjacent brackets, and the results can then be aggregated, leveraging the structure of the kink problem.
The payoff and the budget constraint faced by unit $i$ in this setup can be written as
Here, $A_0(y)$ and $A_1(y)$ represent the after-tax income defined by the linear tax schedules applied to the lower and upper tax brackets, respectively. We define $\tau_0 < \tau_1$ as the marginal tax rates for each bracket and $I$ as the intercept at the cutoff, denoted by $K$. Then, we can write $ A_d(y) = (1-\tau_d) (y - K) + I $ for each $d \in \{ 0,1 \}$. This implies that the actual after-tax income corresponding to the taxable income $y$ is given by $$ A(y) = (1-\tau_0)(y-K) \mathbbm{1}\{ {y \le K} \} + (1-\tau_1)(y-K) \mathbbm{1}\{ {y > K} \} + I, $$ which is a piecewise linear function of $y$ with a negative kink at $y = K$.
For each $d \in \{ 0,1 \}$, we define the counterfactual choice $Y_{{i}}^*(d)$ as the potential income that unit $i$ would have chosen if they had been subject to the linear tax schedule $c = A_d(y)$. I.e., $Y_{{i}}^*(d)$ is defined as the maximizer of the following counterfactual payoff function $U_d$:
Here, we have concentrated out $c$ from the original utility function by substituting in the budget constraint $c = A_d(y)$. This indirect utility representation allows us to focus on the choice of $Y$. By the first-order condition and the strict concavity of $U_d$, this implies the counterfactual choice
Using the counterfactual payoff functions defined above, we can formulate the actual income choice in (ref) as the solution to the following unconstrained problem:
where the actual payoff $U$ under the piecewise linear tax schedule is given by
This compound payoff function reflects the kink in the actual income tax schedule. Since $D_y U_0(y,\eta,\theta_0)|_{y=K} > D_y U_1(y,\eta,\theta_0)|_{y = K}$, the actual payoff function also exhibits a negative kink at $K$:
The actual payoff function is strictly concave in $y$ and smooth at all points except the cutoff.
As demonstrated in saezTaxpayersBunchKink2010, the actual income choice can be solved as follows:
Here, $[\underbar \eta ,\bar \eta]$ is given by $ \left[ K (1-\tau_0)^{-\theta_0}, K (1-\tau_1)^{-\theta_0}\right], $ which is a function of the unknown $\theta_0$ and other known parameters. Such an interval is referred to as the bunching interval, as it characterizes the range of the bunching subpopulation in the entire distribution.
We expand beyond the isoelastic model and explore bunching in a broader context. We consider a situation where two policies, indexed by $0$ and $1$, present smooth but different incentive schedules that are defined across the entire range of $y$. The distinct incentive schedules are embedded in the counterfactual payoff functions $U_d$, which represent the hypothetical utility that would result under policy $d \in \{ 0,1 \}$.
The counterfactual choice $Y_{{i}}^*(d)$ is defined as the maximizer of $U_d$ with respect to the argument $y$. For brevity, we consider a scalar running variable $Y \in \mathcal{Y}$, where $\mathcal{Y}$ is a compact choice interval in $\mathbb{R}$. We assume there is a single cutoff $K$ in the actual policy that lies in the interior of $\mathcal{Y}$. In the actual policy, whether $Y$ exceeds $K$ determines the treatment under each policy. The incentive schedule for policy $0$ is applied when $Y \le K$, and the incentive schedule for policy $1$ is applied otherwise. Polices $0$ and $1$ are also referred to as the prekink and postkink policies, respectively.
We allow for observed covariates in the counterfactual payoff function, denoted by a vector $X \in \mathcal{X} \subseteq \mathbb{R}^{d_x}$. This set of variables represents observably heterogeneous factors involved in the decision-making process of the running variable. In the tax context, $X$ may include demographic characteristics, occupation, additional income sources, etc., depending on the modeling details. As before, $\eta \in \mathcal{H} \subseteq \mathbb{R}$ represents a scalar unobserved preference heterogeneity. The counterfactual payoff is then summarized as a function $U_d(y,x,\eta)$.
Now, we make the following assumptions about the counterfactual payoff functions.
Condition (ref) requires that those choosing $Y < K$ find the postkink policy more beneficial than the prekink policy, and the converse is also true. The two policies offer the same payoff at the cutoff, which ensures the continuity of the compound payoff function. Such situations can be commonly found in many progressive tax systems and regulatory policies.
Condition (ref) introduces bunching in this general setup. Assuming differentiability, Condition (ref) implies that $ D_{y} U_0(y,x,\eta)|_{y = K} \ge D_{y} U_1(y,x,\eta)|_{y = K}. $ If the inequality is strict, the compound payoff function, which can be written as $$ U(y,x,\eta) = U_0(y,x,\eta) \mathbbm{1}\{ {y \le K} \} + U_1(y,x,\eta) \mathbbm{1}\{ {y > K} \}, $$ displays a negative kink at the cutoff. Consequently, the distribution of the actually chosen $Y^*$ exhibits a positive point mass at $K$ if there is a continuum of $\eta$ that satisfies
Assumption (ref) can be straightforwardly modified to allow for payoff maximization with equality constraints or multiple choice variables. Equality constraints can be integrated with the payoff function using the Lagrangian or direct substitution, as in the case of the isoelastic model.
The next proposition provides a simple connection between the actual and the counterfactual choices under Assumption (ref).
goffTreatmentEffectsBunching2024 (goffTreatmentEffectsBunching2024, {Lemma 1}) derived the same relationship in the payoff maximization problem under piecewise smooth budget constraints. Proposition (ref) offers two useful insights into the distribution of $Y^*$. First, the segments of $Y^*$ to the left and right of the cutoff reveal the distributions of the prekink and postkink counterfactual choices, respectively. Second, the bunching mass is the intersection of the right-censored portion of the prekink counterfactual distribution and the left-censored portion of the postkink counterfactual distribution.
Assume the structural equations for $Y^*(d)$ specified as follows:
The structural parameter is denoted by $\theta_0 \in \Theta \subseteq \mathbb{R}^{d_\theta}$. We assume that there exists a utility model satisfying Assumption (ref), from which the structural equations (ref) can be derived. For that utility model to satisfy Condition (ref), the structural equations must satisfy
in light of Proposition (ref). Conversely, if the structural equations meet Condition (ref), they are compatible with a utility model satisfying Condition (ref).
To proceed with our analysis, we require the following assumptions about the structural model.
Assumption (ref) imposes that a scalar unobservable heterogeneity $\eta$ determines the common quantiles of both counterfactual choices conditional on $X_i = x$. This assumption is standard in the literature on structural estimation ({matzkinNonparametricEstimationNonadditive2003}, {matzkinNonparametricEstimationNonadditive2003}; \citeauthor*{blundellIndividualCounterfactualsMultidimensional2017}, blundellIndividualCounterfactualsMultidimensional2017, etc). A similar rank invariance condition can be found in chernozhukovIVModelQuantile2005 in the context of IV quantile estimation. We may allow for an extension to cases where $y$ is a vector containing a running variable $y_1$ and $\eta$ is of the same dimension as $y$, although such an extension lies outside the scope of this paper.
A sufficient condition for the utility model to satisfy Assumption (ref)(i) is the following: for each $d \in \{ 0,1 \}$ and $x \in \mathcal{X}$, $(y,\eta)\mapsto U_d(y,x,\eta)$ has single crossing differences in $\eta$.\footnote{If $U_d(y,x,\eta)-U_d(y',x,\eta) \ge 0 $, then $U_d(y,x,\eta')-U_d(y',x,\eta') \ge 0$ for all $\eta' \ge \eta$. } This condition is implied by a stronger condition of increasing differences: ${\partial^2 U_d(y,x,\eta)}/{\partial y \partial \eta} \ge 0$ for each $d$ and $x$. Assumption (ref)(ii) is a standard regularity assumption ensuring continuous responses of individuals.
An important implication of Assumption (ref)(i) is the existence of the reversion, defined as
The reversion function couples one counterfactual choice with the other choice via the identities $Y_{{i}}^*(0) = R(Y_{{i}}^*(1),X_i,\theta)$ and $Y_{{i}}^*(1) = R^{-1}(Y_{{i}}^*(0),X_i,\theta)$. Here, $R^{-1}(y,x,\theta)$ represents the inverse of $R(y,x,\theta)$ with respect to the argument $y$, which is strictly increasing by Assumption (ref). Applying the reversion on both sides of $Y^*(1) \le K$, the bunching condition can be written equivalently as
Here, $R(K,x,\theta_0)-K \ge 0$ represents the propensity to bunch for those with $X_i = x$. This implies that the bunching moment is determined by the marginal distribution of $Y^*(0)$ and the conditional distribution of $R(K, X, \theta_0)$ given $Y^*(0)$, which highlights a reduction in dimensionality from $X$.
The spread of genuine bunchers around the cutoff is referred to as fuzzy bunching in this paper. The lack of tight bunching poses challenges to the bunching method, as the true bunching status is obscured by those incidentally located around the cutoff without the intention of bunching.
Fuzzy bunching has typically been attributed to optimization frictions in the literature. Various sources of these frictions, such as adjustment costs, inattention, and other factors, have been extensively discussed in the literature for their roles in dampening the observed elasticity (e.g., {chettyAdjustmentCostsFirm2011}, {chettyAdjustmentCostsFirm2011}; {klevenUsingNotchesUncover2013}, {klevenUsingNotchesUncover2013}). This paper chooses to treat these frictions as measurement errors without specifying their nature. A key assumption is the existence of a narrow window encapsulating the range of measurement errors. This is specified in the following assumptions along with those on the data-generating process.
We assume that the data are i.i.d. A form of dependence may also be accommodated, although such extensions are beyond the scope of this paper. The counterfactual distribution must be continuous according to Assumption (ref)(iii). This assumption is standard in the literature and implies no point mass observed in the absence of a kink.
Assumption (ref)(i) requires that observed choices outside the window equal the desired choices. This assumption is crucial for us to infer the counterfactual distribution from the observed distribution. This assumption could be economically justified, e.g., in a model where bounded adjustment costs are incurred conditional on one's decision to bunch, as in chettyAdjustmentCostsFirm2011.
blomquistBunchingIdentificationTaxable2021 have adopted the same assumption about the optimization error window. This assumption is also important in the polynomial strategy, which is a widespread practice in the bunching method. Practically, the choice of the window is often guided by visual inspection, which can be validated through some robustness checks.
There are alternative approaches to addressing optimization errors. In the notch literature, it is relatively common to model optimization frictions as a proportion of individuals who do not respond to the notches ({klevenUsingNotchesUncover2013}, {klevenUsingNotchesUncover2013}; {bestEstimatingElasticityIntertemporal2020}, {bestEstimatingElasticityIntertemporal2020}, etc). Recently, \citet*{anagolDiffuseBunchingFrictions2024} proposed a different model-based approach, where each unit $i$ receives a sparse set of offers $\{ Y_{{i}}^* + \varepsilon_1,\ldots,Y_{{i}}^* + \varepsilon_M \}$ and locates at $Y_{{i}}$ that maximizes their payoff among these offers. The offers are randomly drawn in the vicinity of the desired choice following a certain distribution, represented by $\varepsilon_1,\ldots,\varepsilon_M \overset{iid}\sim F_\varepsilon$. Their procedure is operationalized without specifying a priori range of optimization errors like ours. However, it requires researchers' insights into $M$ and the distribution of the offers. Their treatment of optimization frictions can seemingly be integrated with our method based on the analytic model, which could be a fruitful avenue for future research.
This section presents our bunching identification results. We first propose our semiparametric approach based on the analyticity condition of the counterfactual functions. It is demonstrated that the structural parameter can be point-identified as an implicit function of the bunching moment under the analyticity condition and mild regularity assumptions.
Our second approach relaxes the analyticity condition on the counterfactual density used in the first approach and considers instead (piecewise) analytic envelope functions bounding the counterfactual distribution. This approach accommodates a range of counterfactual scenarios, albeit at the expense of resulting in potentially multiple identified values of $\theta$.
blomquistBunchingIdentificationTaxable2021 demonstrated that without imposing further restrictions on the unobservable counterfactual density, the elasticity is unrestricted by the observables in the isoelastic model. Such identification failure occurs similarly in the general bunching design, intuitively because the bunching moment depends on latent component of the counterfactual distribution, which is independent of both the structural parameter and the observables.
This implies that the econometrician must narrow a model of the counterfactual density to reconstruct the unobservable counterfactual from the observables. We first discuss the required shape restrictions on the counterfactual density that allow point identification. Let $\mathcal{F}$ denote a model of $f_{Y^*(0)}: \mathcal{Y}_0 \to \mathbb{R}_+$. We assume that if $f\in \mathcal{F}$, $\tilde{f} \in \mathcal{F}$, and $f(y) = \tilde{f}(y)$ for all $y < K_0$, then $f \equiv \tilde{f}$ on $\mathcal{Y}_0$. If $\mathcal{F}$ satisfies this condition, we say that the model possesses the identity property. The identity property for the model of the counterfactual CDF was previously proposed by bertanhaBetterBunchingNicer2023 (bertanhaBetterBunchingNicer2023, {Assumption 1}).
Various models enjoy the identity property: many parametric classes such as exponential families with smooth sufficient statistics, their discrete and continuous mixtures, finite-degree polynomials, and most crucially for our analysis, the analytic class in light of the identity theorem from complex analysis (steinComplexAnalysis2010, steinComplexAnalysis2010). Generally speaking, there is a trade-off between a model's flexibility to fit the data and its capability of allowing unique extrapolation. For instance, many nonparametric classes, such as $C^k$ (where $k \in \mathbb{Z}_+$), Hölder spaces, or Sobolev spaces, offer nearly full flexibility in data fitting. However, at the expense of this freedom, they cannot provide a unique extrapolation into the unobservable region, highlighting the absence of the identity property. In this paper, we investigate the use of the analytic function class for bunching identification in a general structural model, as a semiparametric approach to balancing these conflicting goals. Similar ideas can be found in the econometric literature. Since its earliest version in 2020, pollinger2024kinks has explored the usefulness of the analyticity assumption for identifying the intensive and extensive margin elasticities in an extended isoelastic model. With similar motivation but in a different context, iariaRealAnalyticDiscrete2024 examined the role of real analyticity in identification, extrapolation, and numerical implementation in discrete choice models.
Now, we describe our first approach to point identification. For simplicity, we focus on the case where $\theta_0$ is a scalar parameter. However, it is possible to address a vector of structural parameters by developing multiple bunching moments. Throughout this paper, we follow the common convention of defining the status quo counterfactual choice as $Y^*(0)$ rather than $Y^*(1)$. Nevertheless, the roles of $Y^*(0)$ and $Y^*(1)$ can be interchanged by treating $\{ (-Y_{{i}}, X_i) \}_{i=1}^n$ as a new dataset with the corresponding adjustments to the cutoff and the optimization error window.
Let $T_i = T(X_i)$ be a bounded, nonnegative function of the covariates such that $E[T_i] > 0$. Given $T_i$, we define a new probability measure as ${P_{T}}(A) = E[\1_A T_i]/E[T_i]$. This offers a measure of bunching that may be different from the unweighted fraction, based on the weighting function $T(.)$. For various choices of $T(.)$, these measures can provide a set of moments that can be used as identifying restrictions for $\theta_0$. If $T_i$ is set to a constant 1, ${P_{T}}$ reduces to $P$, which is sufficient for identifying a scalar parameter.
For each $\tau \in (0,1)$ and $\theta \in \Theta$, we define
as the counterfactual density of $Y^*(0)$ and the conditional quantile function of $R(K_1,X,\theta)$ given $Y^*(0)$, both measured with respect to ${P_{T}}$. Here, $R(y,x,\theta)$ denotes the reversion as defined in (ref). We may interpret $R(y,x,\theta)$ as the predicted value of $Y^*(0)$ in response to the shift from policy $1$ to $0$, given that the unit was at $Y^*(1) = y$ and $X = x$. In this sense, $R(K,X,\theta) - K$ quantifies the predicted response to a potential repeal of the kink for individuals located at the right edge of the bunching population.
A function $f:[a,b] \to \mathbb{R}$ is said to be analytic if, for every $y \in [a, b]$, there is an open neighborhood around $y$ on which $f$ agrees with its Taylor expansion at $y$. Moreover, it is well-known that an analytic function $f:[a,b]\to\mathbb{R}$ can be extended to a complex-analytic function $\tilde{f}:U\to\mathbb{C}$ over an open planar set $U \subseteq \mathbb{C}$ containing $[a,b]$. Our identification theory leverages rich results in the theory of complex analysis applicable to such analytic continuations.
We now formally define the analyticity condition of the counterfactual density and conditional quantile functions.
The smoothness constants $(\rho,\beta,\delta)$ control the variability of the functions $f$ and $(g_{{\tau}})_{\tau \in (0,1)}$. The constant $\rho$ is at least as large as the radii of convergence of these functions. Given the value of $\rho$, the constants $\beta$ and $\delta$ more concretely govern variations in $f(y)$ and $g_{{\tau}}(y)$ within the $\rho$-neighborhood of $[a,b]$. In Assumption (ref), we require that $\delta < 1$ to ensure that bunching does not extend to the upper limit of the support. The intuition is that the higher $Y^*(0)$ is, the fewer incentives there should be for bunching, as will be clarified in the proof of Theorem (ref). Further technical remarks are relegated to the Appendix.
For each $\theta \in \Theta$, let $\bar{K}_1(\theta) \in [K, \bar {\mathcal{Y}}_0]$ be any constant satisfying ${P_{T}}( R(K_1, X_i,\theta) > \bar{K}_1(\theta), Y_{{i}} \in [K_0, K_1] ) = 0$. We impose the analyticity condition on the counterfactual density and the conditional quantile functions as follows.
Assumption (ref)(i) plays a central role in deriving our identification results. Loosely speaking, it requires the analyticity of three functions: the counterfactual density $f_{Y^*(0)}$, the conditional expectation function of $T_i$, and the conditional quantile functions of $R(K_1,X_i,\theta_0)$. Consequently, these functions must be uniformly approximated by a sequence of finite-degree polynomials with exponentially decreasing errors ({timanTheoryApproximationFunctions1994}, {timanTheoryApproximationFunctions1994}). Since the distribution of $Y_{{i}}^*(0)$ is observed only in a censored manner, the analyticity of these functions is untestable in the unobservable region. Moreover, there are no formal statistical procedures for testing the analyticity of densities or quantile functions, to the best of my knowledge, even when they are fully observable. Thus, we suggest visually supporting this assumption by examining various bin-based estimates of these functions in the observed region.
A sufficient condition for Assumption (ref)(i) is that the joint density of $(Y_{{i}}^*(0), R(K_1,X_i,\theta))$ is analytic and bounded away from $0$ in the support. This assumption can be further supported by the analyticity of the joint density of $(\eta_i, X_i)$. Such a condition would restrict $f_{Y^*(0),R(K_1,X,\theta)}$ to be a continuous function, as well as their derivatives.
Assumption (ref)(i) may appear to restrict $R(K_1,X_i,\theta)$ to be continuously distributed. However, it can be extended to cases where $R(K_1,X_i,\theta)$ has finitely many non-overlapping components in its support by conditioning on membership to each component. This situation arises, for example, when $X_i$ is a mixture of dummy variables and continuous covariates. In such cases, the requirement is that each conditional ${P_{T}}$-probability assigned to each component should be an analytic function of $y$.
Assumption (ref)(ii) stipulates that the bunching moment predicted by the model must be able to differentiate the true value from the others, according to the true distribution of $(Y^*(0), X)$. This condition ensures global identification through the comparison of moments in Theorem (ref). In the case where $\theta$ is a scalar parameter, the identification condition is fulfilled when the reversion $R$ is monotonic in $\theta$.
Now, we are prepared to present our point identification result.
Our point identification result crucially rests upon Assumption (ref)(i) to express a bunching moment as a series of functionals involving counterfactual quantities. If the analyticity condition is violated, our method would identify a pseudo-true value arising from the fitted distribution within the analytic model.
To see why $B(\theta)$ can be inferred from the observables, we use the fact that $Y_{{i}} = Y_{{i}}^*(0)$ for all units $i$ such that $Y_{{i}} < K_0$ by Proposition (ref). Then, we can express the $j$th summand on the right-hand side of (ref) as
which is determined by the distribution of $(Y,X)$. As a result, $\theta_0$ can be identified.
The observed vector $X$ affects the function $B(\theta)$ through the conditional moments of $R(K_1,X,\theta)$ given $Y^*(0)$. Notably, it is possible to accommodate a moderately large or high-dimensional $X$ without compromising efficiency, so long as $R(.,\theta)$ depends on a low-dimensional parameter $\theta$. Previous approaches have either partitioned the sample based on a limited number of $X$ values or simply consolidated $X$ into its empirical average in the data, which could serve as an approximation to the aggregation in (ref). Such methods may lead to concerns about aggregation bias or information loss, which can be addressed in our approach.
Consider the model $Y_{{i}}^*(0) = R(Y_{{i}}^*(1),X_i,\theta)$, where the observed bunching is assumed to be tight at $K$. Assume, however, that the correct specification is $Y_{{i}}^*(0) = R(Y_{{i}}^*(1),X_i,\theta_i)$, where $\theta_i \in \mathbb{R}$ represents a unit-specific parameter. Then, the pseudo-true value $\theta_*$ identified through the bunching method corresponds to the solution to the equation
Suppose that $f_{Y^*(0)}$ is flat around the bunching interval, that is, $D^{j} f_{Y^*(0)}(K) \approx 0$ for all $j \in \mathbb{N}$. By a slight extension of Theorem (ref), this implies that the equation above expands to
which shows that $\theta_*$ approximately satisfies $$ E[R(K, X_i, \theta_*)| Y_{{i}}^*(0) = K] \approx E[R(K, X_i, \theta_i)| Y_{{i}}^*(0) = K]. $$ In this sense, when $\theta_i$ is misspecified but $X_i$ is correctly specified, the pseudo-true value approximately identifies a weighted average of $\theta_i$: $$ \theta_* \approx \frac{E[ D_{\theta} R(K, X_i, \theta_*) \cdot \theta_i | Y_{{i}}^*(0) = K]}{E[ D_{\theta} R(K, X_i, \theta_*)| Y_{{i}}^*(0) = K]}, $$ where the weights are nonnegative if $R(K,x,\theta)$ is monotonic in $\theta$ for each $x \in \mathcal{X}$. This reaffirms a similar relationship derived in klevenUsingNotchesUncover2013 in the absence of $X$.
\paragraph*{Comparison to existing results} There are several existing approaches to bunching identification in the literature. The most closely related to our result is the small kink approximation, initially proposed by saezTaxpayersBunchKink2010 as a nonparametric identification method. It says that the normalized bunching fraction approximately identifies the local average of individual responses to policy shift at the cutoff:
provided that either the policy response $Y^*(0)-Y^*(1)$ is small or the counterfactual density is flat around the bunching interval. In the isoelastic model, after a further substitution $Y_{{i}}^*(0)-Y_{{i}}^*(1) = \theta_0 \log \left( \frac{1-\tau_0}{1-\tau_1}\right) K$ into (ref), this approximately identifies the true elasticity as
This identification method has been a predominant approach in the kink and notch design literature, with its theoretical foundation and applications further developed in subsequent works, such as chettyAdjustmentCostsFirm2011, klevenUsingNotchesUncover2013, blomquistBunchingIdentificationTaxable2021, and goffTreatmentEffectsBunching2024, among many others.
The small kink approximation relies on the first moment of the individual policy responses and the value of the counterfactual density at the cutoff to approximate the bunching fraction. Our result complements this approximate relationship with the omitted higher-order effects of the individual policy responses and the derivatives of the counterfactual density. Overlooking these terms may lead to significant bias in the counterfactual estimation, which can be addressed in our procedure.
bertanhaBetterBunchingNicer2023 proposed two semiparametric identification approaches to the elasticity in the isoelastic model. The first approach is based on the mid-censored Tobit regression facilitated by the observed covariates $X$. It is based on the quasi-maximum likelihood estimation (QMLE) of the following model:
where $\eta = X'\beta_0 + e$ and $e | X \sim N(0,\sigma_0^2)$ for some $\beta_0 \in \mathbb{R}^{d_x}$ and $\sigma_0^2 > 0$. The conditional normality assumption enables estimation of $(\theta_0, \beta_0, \sigma_0^2)$ through the QMLE on the observable $(y,X)$. bertanhaBetterBunchingNicer2023 emphasized that the identification of $\theta_0$ does not require the full power of the conditional normality assumption. Under a weaker condition that the unconditional distribution of $\eta$ coincides with a mixture of Gaussian distributions with mean $X'\beta_0$ and variance $\sigma_0^2$, they demonstrated that $\theta_0$ can be identified and consistently estimated.
Their identification hinges on the assumption $F_{\eta}(h) = E_X[\Phi( \frac{h - X'\beta_0}{\sigma_0})]$ for some pseudo-true values $(\beta_0,\sigma_0^2)$, where $\Phi$ denotes the normal CDF. They incorporated $X$ into the model primarily to improve its fit to the counterfactual distribution. In contrast, the our main motivation is to enhance the flexibility of the model. Notably, the Gaussian mixture CDF $F_\eta$ defined above satisfies Assumption (ref)(i).\footnote{We prove the global analyticity of $F_\eta$ and $f_\eta$ in the Appendix.} Therefore, $\theta_0$ can be identified through our means without requiring additional covariates. However, improving the estimation procedure by incorporating information provided by $X$ would be an interesting avenue for future research.
Their second approach is based on the restriction on a conditional quantile of $\eta$. They assume, for some $\tau \in (0,1)$ and $\beta_0$, the conditional $\tau$ quantile of $\eta$ given $X$ is given by $ Q_{\eta|X}(\tau, x) = x'\beta_0. $ This, in turn, implies that the conditional $\tau$ quantile of $y$ given $X$ assumes a parametric form similar to the mid-censored Tobit regression model. This can be estimated using a quantile regression under a suitable rank condition on $X$. This approach provides an alternative identification means through the conditional quantile restriction on the counterfactual distribution, instead of restricting the entire shape. On the other hand, it does not allow for semi- or nonparametric specification of the conditional quantile function due to the rank condition. Our approach would serve as a valuable complement in situations where researchers are uncertain about such restrictions.
We relax the analyticity condition on the counterfactual density by allowing the counterfactual cumulative distribution function to be bounded by a pair of analytic envelope functions. Specifically, we assume the following.
Condition (ref) imposes an inequality restriction similar to first-order stochastic dominance, except that $\bar{f}$ and $\underaccent{\bar}{f}$ are not required to be valid probability densities. Note that (ref) is implied by
However, (ref) provides a slightly more general restriction than the direct bounds on the counterfactual density.
Our approach to partial identification, which embraces a range of counterfactual scenarios, allows for identification under less stringent assumptions at the expense of potentially resulting in a continuum of identified values. It can accommodate several parametric and nonparametric shape restrictions in the literature, including the trapezoidal approximation ({saezTaxpayersBunchKink2010}, {saezTaxpayersBunchKink2010}), densities with bounded oscillations ({blomquistBunchingIdentificationTaxable2021}, {blomquistBunchingIdentificationTaxable2021}), bi-log-concave CDF ({goffTreatmentEffectsBunching2024}, {goffTreatmentEffectsBunching2024}), and Lipschitz continuous densities ({bertanhaBetterBunchingNicer2023}, {bertanhaBetterBunchingNicer2023}). We derive a set of sharp identification bounds in Theorem (ref) subject to the restriction (ref). That is, the lower bound is attained when $\underaccent{\bar}{f} \equiv f^{({P_{T}})}_{Y^*(0)}$, and the same applies to the upper bound. When applied to the isoelastic model, Theorem (ref) yields the same bounds as those established in blomquistBunchingIdentificationTaxable2021 (blomquistBunchingIdentificationTaxable2021, {Theorem 2}) and bertanhaBetterBunchingNicer2023 (bertanhaBetterBunchingNicer2023, {Theorem 2}) subject to the respective shape restrictions. This assertion is demonstrated in the Appendix. We now state the analyticity assumption on the envelope functions.
The following theorem presents our partial identification result.
It is possible to relax Assumption (ref) in Theorem (ref) to allow for piecewise analytic envelope functions. For this extension, we refer to Theorem (ref) in the Appendix. Such an extension is useful for handling piecewise linear envelope curves arising from the Lipschitz constraint ({bertanhaBetterBunchingNicer2023}, {bertanhaBetterBunchingNicer2023}).
The idea of bunching confidence region is to collect values of $\theta$ that are compatible with the observed pattern of bunching. More specifically, the identification results in Section (ref) enable us to draw empirically testable implications from the structural parameter, expressed as either moment equality or inequality restrictions. We construct a test for these restrictions under the hypothesis $H_0 : \theta_0 = \theta$.
In this section, we first present our counterfactual estimation scheme, termed the generalized polynomial strategy. We delineate our procedure based on the counterfactual correction and polynomial sieve estimation. An asymptotically valid test is proposed with the growth conditions to guide the proper selection of tuning parameters. Lastly, we propose a modified procedure for the inference on values in a partially-identified set derived from Theorem (ref).
This subsection details our counterfactual estimation scheme, referred to as the generalized polynomial strategy. We assume conditions in Theorem (ref) for the point identification of $\theta_0$. Our aim is to construct a confidence set $C_n$ satisfying $ \liminf_{n\to\infty} P(\theta_0 \in C_n) \ge 1- \alpha, $ where $\alpha \in (0,1)$ is a prescribed coverage level. It is well-known that such a confidence region can be constructed by inverting a test for the hypothesis $H_0: \theta_0 = \theta$, which will be our focus.
Let $\theta \in \Theta$ be a hypothesized value and assume $H_0:\theta_0 = \theta$ is correct. Using Theorem (ref), we can express the observed bunching moment as the following infinite sum:
We test this relationship by estimating the functions $f_j(y) := E[T_i w_i^j |Y_{{i}}^*(0) = y] f_{Y^*(0)}(y)$. Here, we denote $w_i = R(K_1,X_i,\theta) - K_0 \ge 0$, where the nonnegativity follows from the restriction (ref) and the monotonicity of the reversion.
Our procedure comprises two steps: (i) the construction of the estimation sample and (ii) the sieve estimation of the counterfactual functions. We present a pseudo-code in Algorithm (ref) to provide an overview of the procedure. The construction of the estimation sample necessitates an appropriate counterfactual adjustment procedure. We name the proposed method as the counterfactual correction. The method crucially hinges on the existence of the reversion, $R(y, x, \theta)$, which is directly applied to transform the data.
Since $f_j$ is analytic on $\mathcal{Y}_0$, as established in Lemma (ref), these functions as well as their derivatives can be globally approximated by a sequence of polynomials of increasing degrees. This underpins our estimation scheme using the polynomial sieve. For each $j \in \mathbb{N}$, note that $f_j$ is proportional to the density $D_y {P}_{j}(Y_{{i}}^*(0) \le y)$, where ${d{P}_{j}}/{d P} = T_i w_i^j/E[T_i w_i^j]$. Motivated by $$ f_j \in \operatornamewithlimits{arg\hspace{0.1em} max}_{{f > 0}}\hspace{0.1em} E[T_i w_i^j \log (f(Y_{{i}}^*(0))) \mathbbm{1}\{ {Y_{{i}}^*(0) \in \mathcal{S}} \}] - \int_{\mathcal{S}} f(y)dy $$ for any region $\mathcal{S} \subseteq \mathcal{Y}_0$, we propose a sieve M-estimator for $f_j$ with mid-censored data. By recasting the problem as the estimation of a reweighted density, we can estimate $f_j$ without binning the observations.
\SetKwInOut{KwFinOut}{Final output}
This subsection presents our proposal for counterfactual adjustment that allows consistent estimation of the counterfactual quantities. For the estimation of the counterfactual related to $Y^*(0)$, we need inputs that correctly reflect the distribution of $Y^*(0)$ despite censoring. On the one hand, one can observe $f_{Y^*(0)}$ to the left of the excluded window. This enables us to extrapolate $f_{Y^*(0)}$ into the unobservable region. However, using only the data in the left segment is not the most efficient idea, as it renders counterfactual estimation highly noisy. There are significant benefits in taking advantage of the data to the right of the window for improving counterfactual estimation. To this end, an appropriate counterfactual adjustment procedure is necessitated since the data on the right represent $f_{Y^*(1)}$, not $f_{Y^*(0)}$.
The proportional adjustment is a common counterfactual adjustment scheme in the polynomial strategy in the presence of an excess bunching fraction. In the proportional adjustment procedure, the adjustment of the form $\hat f_{Y^*(0)} = (1+c_j) \hat f_{Y^*(1)}$ is sequentially applied to $\hat f_{Y^*(1)}$, where the factor $c_j>0$ is updated so as to ensure the resulting estimate of $\hat f_{Y^*(0)}$ integrates to $1$. However, the proportional adjustment may not align with the structural assumptions, leading to incorrect analysis. For instance, assume that $Y^*(0) = R(Y^*(1))$. Then, the change of variables implies that $f_{Y^*(0)}(y) = f_{Y^*(1)}(R^{-1}(y)) \frac{1}{R'(R^{-1}(y))}$, which in general differs from $f_{Y^*(0)}(y) = (1+c) f_{Y^*(1)}(y)$.
Our counterfactual correction procedure leverages the truth of the hypothesis $H_0:\theta_0 = \theta$. The procedure can be described as follows. Let $(Y_{{i}}, X_i)$ be a data point. If $Y_{{i}} > K_1$, we impute the missing value of $Y_{{i}}^*(0)$ as $Y_{{i}}(0) = R(Y_{{i}},X_i,\theta)$ by applying the reversion. If $Y_{{i}} < K_0$, we define $Y_{{i}}(0) = Y_{{i}}$. If $Y_{{i}} \in [K_0, K_1]$, unit $i$ is excluded from the estimation sample. Iterate this over all units and collect the values of $(Y_{{i}}(0),X_i)$ for the retained units. The initial estimation sample is then defined as $\{ (Y_{{i}}(0),X_i):Y_{{i}}<K_0 \text{ or } Y_{{i}} > K_1 \}$, where
Assuming that $\theta_0 = \theta$ and the structural equations are correctly specified, we have $Y_{{i}}(0) = Y_{{i}}^*(0)$ for all units in the estimation sample. However, this does not immediately imply that valid counterfactual estimation is possible based on this estimation sample. To see this point, we note that
This equation accounts for the fact that units engaged in bunching are excluded from the estimation sample not based on their $Y^*(0)$ values, but on a random selection criterion. Since the proportion of bunchers may be nonzero for some values of $Y_{{i}}(0)$ in the estimation sample, a further truncation is necessitated.
We resolve this issue by truncating $Y_{{i}}(0)$ above a threshold beyond which the bunching fraction becomes $0$. To implement this, we recall the bunching condition $\{K_0 \le Y_{{i}}^*(0) \le R(K_1, X_i, \theta)\}$. This implies that the bunching will not occur beyond
where $\mathcal{B}_i = \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1],T_i > 0} \}$ indicates the observed bunching status. By construction, we have $\bar{K}_1(\theta) \ge \max_{i: \mathcal{B}_i=1}Y_{{i}}^*(0)$. Under mild conditions, $\bar{K}_1(\theta)$ is a consistent estimator for an upper bound of $Y^*(0)$ among bunchers. This allows us to treat $\bar{K}_1(\theta)$ as if it were a deterministic constant for the sake of counterfactual estimation.
This leads to the final estimation sample, defined as $\{ (Y_{{i}}(0),X_i): Y_{{i}}(0) \in \mathcal S(\theta)\}$, where $\mathcal S(\theta) := [\underaccent{\bar}{\mathcal{Y}}_0,K_0]\cup(\bar{K}_1(\theta),\bar {\mathcal{Y}}_0]$ denotes the support of $Y_{{i}}(0)$ in the estimation sample. This sample will be used for the counterfactual estimation related to $Y^*(0)$. For notational convenience, the estimation sample will be denoted as $\{ (Y_{{i}}(0), X_i)\}_{i \in N(\theta)}$. Here, $ \mathbbm{1}\{ {i \in N(\theta)} \} = \mathbbm{1}\{ {Y_{{i}}(0) \in \mathcal S(\theta)} \} = \mathbbm{1}\{ {Y_{{i}}^*(0) \in \mathcal S(\theta)} \}$ indicates unit $i$'s membership in the estimation sample. Hereafter, we omit the dependence on $\theta$ where it is unlikely to cause confusion.
This section describes our method to estimate $f_j$ for each $j \in \mathbb{N}$. Let $\mathcal{P}_k$ denote the linear span of all polynomials of degree less than $k \in \mathbb{N}$. The polynomial basis of order $k\in\mathbb{N}$ is defined as
We denote the coefficient as $\gamma = \gamma_k \in \mathbb{R}^{k}$ with $k$ suppressed for simplicity. We define our sieve space of order $k$ for $f_j$ as
where $c_{3}>0$ is chosen such that $f_{k j}$, the best approximation of $f_j$ in $\mathcal{F}_{k j}$ in a suitable sense, lies in the interior of $\mathcal{F}_{k j}$. The linear growth of $|\log f_j|$ as $j$ increases is ensured by Lemma (ref). This manifests that our search space for $f_{kj}$ will be limited to functions comparable in size to $f_j$. We will not be specific about this constant as it mainly serves as a technical and numerical device. With some abuse of notation, we denote $\gamma \in\mathcal{F}_{kj}$ if $z_k'\gamma \in \mathcal{F}_{kj}$.
For illustration, we focus on a single bunching moment generated by a weighting function $T_i$. If there are multiple bunching moments, they can be addressed individually in a similar fashion. For each $k \in \mathbb{N}$ and $j \in \mathbb{N}$, define $\hat \gamma_{k j} \in \mathbb{R}^k$ as
Using the Lagrangian method, $\hat{\gamma}_{\kappa j}$ is the solution to the following convex optimization problem:
Assuming an interior solution, the first-order condition for $\hat\gamma_{k j}$ is given by
This leads to our estimator for $f_j$ using the order $k$ polynomial sieve:\footnote{In practice, it is recommended to use an orthogonal polynomial basis for numerical stability, which can be obtained by applying a linear basis change.}
where $\hat \gamma_{k j,j}$ represents the $j$-th component of $\hat \gamma_{k j}=(\hat \gamma_{k j,1},\ldots,\hat \gamma_{k j,k})'$, associated with the term $(y-K_0)^{j-1}$.
For each $k \in \mathbb{N}$ and $1 \le l \le k$, let
where $\hat{E}[T_i \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1]} \}] : = \frac{1}{n} \sum_{i=1}^n T_i \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1]} \}$. Here, $k$ and $l$ represent the sieve and the approximation orders, respectively. These tuning parameters reduce the infinite-dimensional array of coefficients to a finite-dimensional one. They are chosen to be $k = \kappa_n$ and $l = \ell_n$ based on a pre-specified schedule determined by the sample size. We allow $\kappa_n$ and $\ell_n$ to increase with $n$, where the subscript $n$ will be suppressed at times. The growth conditions for $\kappa_n$ and $\ell_n$ will be specified later. Under suitable regularity conditions, we establish that $\hat \mu_n := \hat \mu_{\kappa \ell} = b_{\ell} + o_{P}(1)$ under the null hypothesis (Lemma (ref)), where $b_{l} := \sum_{j=l+1}^\infty \frac{1}{j!} D_{y}^{j-1}[ E[T_i (R(K_1,X_i,\theta)-K_0)^j |Y_{{i}}^*(0) = y] f_{Y^*(0)}(y)]_{y = K_0}$ represents the approximation error at truncation order $l$.
It remains to quantify the uncertainty associated with $\hat \mu_n$. Let $s_{k l i}$ be an influence function in the asymptotic linear expansion of $\hat \mu_{k l}$, i.e., for each $k \in \mathbb{N}$ and $1 \le l \le k$, $$ \sqrt{n}(\hat\mu_{k l} - b_{l}) = \frac{1}{\sqrt n}\sum_{i=1}^n s_{k l i} + o_{P}(\sigma_{k l}). $$ Then, we estimate $\sigma_{k l}^2 = \operatorname{var}(s_{k l i})$ by $ \hat \sigma_{k l}^2 = \frac{1}{n} \sum_{i=1}^n \hat s_{k l i}^2, $ where $\hat s_{k l i}$ is a suitable estimate of the individual influence function that satisfies $$ \frac{1}{n} \sum_{i=1}^n |\hat{s}_{k l i}^2 - s_{k l i}^2| = o_{P}( \sigma^2_{k l}). $$ For example, we can obtain $\hat s_{k l i}$ by collecting the influence functions derived from the score functions of $\hat{E}[T_i \mathbbm{1}\{ {Y_{{i}} \in [K_0,K_1]} \}]$ and $\frac{1}{j}\hat \gamma_{k j, j}$'s. (See Algorithm (ref).) This method is adopted across all our applications as well as in the proof of Theorem (ref).
Let $\hat \sigma_n^2 := \hat \sigma^2_{\kappa \ell}$ denote the estimated standard error of $\hat \mu_n$ and $\mathcal T_n(\theta) := \sqrt{n}|\hat\mu_{n}|/\hat \sigma_{n}$ be the absolute t-statistic. To derive the null distribution of $\mathcal T_n$, we should make additional regularity assumptions on the data-generating process and the choice of the sieve and approximation orders. We first assume that the asymptotic variances of $\hat \mu_{k l}$ are uniformly bounded away from $0$, which is innocuous in most situations.
We assume the existence of the interior solution to the sample and the corresponding population maximization program in the following. We emphasize that this assumption is not too restrictive since $c_{3}$ can be set large enough to ensure an interior solution.
To state the conditions for an appropriate polynomial order, we introduce a measure of the quality of the polynomial extrapolation. For each $k \in \mathbb{N}$, let $H_k \in \mathbb{R}^{k \times k}$ and $o_k(.)$ be defined as $H_{k} = \int_{\mathcal{Y}_0} z_k(y) z_k(y)' dy$ and $o_k(y) = H_k^{-1/2} z_k(y)$, respectively. Then, $o_k(.)$ represents the orthonormal polynomial basis of order $k$ with respect to the Lebesgue measure on $\mathcal{Y}_0$. Let $0 < \chi_k \le 1$ be defined as the smallest eigenvalue of the matrix $\int_{\mathcal S} o_k(y) o_k(y)' dy$. The number $1/\chi_k$ represents the worst-case ratio of the extrapolation errors within $\mathcal S^c$ to the fitted errors in $\mathcal S$, quantifying the stochastic variability of the counterfactual estimation relative to the observed fit. We refer to this constant as the extrapolation norm. The sequence $(1/\chi_k)_{k\in\mathbb{N}}$ increases unboundedly, at most at an exponential rate of $C^{k}$ for some $C > 1$ by Lemma (ref)(ii). The magnitude of this sequence depends positively on the size of the extrapolated region relative to that of the entire support, reflecting the difficulty of extrapolation.
We need to ensure $\chi_{\kappa_n}^{-1}$ does not diverge too quickly relative to the sample size, as it would result in a highly noisy extrapolation into the unobserved region despite the goodness of the observed fit. This underscores the need to gradually expand the sieve space due to the ill-posed nature of the extrapolation problem. The concrete conditions are presented in Assumption (ref), which guarantees consistent counterfactual estimation and correct inference.
Notice that $n^{c_{5}/2} \lesssim \chi_\kappa^{-1} \lesssim C^\kappa$ requires $\kappa \gtrsim \log n$. Thus, the sieve dimension should increase at a logarithmic rate in $n$, where the proportional factor depends on the size of the extrapolated region.
Next, we make regularity assumptions about the smoothness of the counterfactual functions, represented by $\rho$, and the lower bound of the estimated function. The assumed conditions ensure that the errors from polynomial approximation decay sufficiently fast.
Given these additional assumptions, it is possible to derive an asymptotic distribution of $\mathcal T_n$, as presented in Theorem (ref) below. It allows us to construct an $\alpha$-sized test for $H_0 : \theta_0 = \theta$ as
Here, $\operatorname{cv}_{{\alpha}}(b)$ represents the $(1-\alpha)$ quantile of $|N(b,1)|$ and $\bar{b}_\ell = \beta \delta^{\ell} E[T_i \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1]} \}]$ serves as an upper bound for the approximation bias derived in Theorem (ref). We may estimate $\bar{b}_\ell$ or another upper bound based on the sample (see the remarks following Theorem (ref)). Consequently, a confidence set with an asymptotic coverage of $(1-\alpha)$ can be constructed as
Some remarks on the choice of $\ell$ and the calibration of approximation bias are in order. First, it could be informative to report confidence sets using critical values with $\bar{b}_\ell$ set to $0$. These confidence sets should cover the pseudo-true values that are identified despite the approximation errors. If $\ell$ is appropriately chosen and the approximation errors are small in absolute terms, these pseudo-true values can still provide meaningful information about the true parameter. We adopt this convention in reporting our results in the Monte Carlo experiments and the empirical application.
Second, to obtain an upper bound for the approximation bias, we can estimate the smoothness constants from the data. To implement this, we first set a grid $\mathcal{G}$ for values of $\rho$. Next, we estimate $f^{({P_{T}})}_{Y^*(0)}(y)$ and $g^{({P_{T}})}(\tau, \theta, y)$ along with their derivatives for $\tau$ picked from a dense grid in $(0,1)$. The conditional quantile function can be estimated using the weighted polynomial sieve quantile regression akin to the proposed estimation procedure. Given estimates for these functions, for each $\rho \in \mathcal{G}$, one can compute $(\beta_\rho,\delta_\rho)$ according to Definition (ref), and bound $|b_\ell| \le \min_{\rho \in \mathcal{G}} \beta_\rho \delta_\rho^{\ell} E[T_i \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1]} \}]$. The bias can be controlled more efficiently using the upper bound provided in Lemma (ref)(iii), which is expressed as a Fourier coefficient related to these functions. Such methods could further support an appropriate choice of the approximation order by evaluating the approximation quality.
This subsection outlines our proposal for inference on the values in the identified set
derived from Theorem (ref). As before, our aim is to construct a confidence region $C_n$ with the correct coverage $ \liminf_{n\to\infty} \inf_{\theta \in \Theta_0}P(\theta \in C_n) \ge 1- \alpha. $ This can be accomplished by inverting a test for $H_0 : \theta_0 = \theta$ that satisfies $\limsup_{n\to\infty} \sup_{\theta \in \Theta_0} E[\phi_n(\theta)] \le \alpha$. We consider a single bunching moment weighted by $T_i$.
We assume the existence of preliminary estimators for the envelope functions and their derivatives, which permit an asymptotic linear expansion around $y = K_0$. This allows for flexible incorporation of prior knowledge about $\underaccent{\bar}{f}$ and $\bar{f}$.
Let $\theta\in \Theta_0$ be a hypothesized value. We then estimate the conditional moment functions, defined as $m_j(y) = E_{P_{T}}[w_i^j | Y_{{i}}^*(0) = y]$. These functions can be estimated in various ways using the weighted sieve regression. For instance, we can define a least-squares estimator as
or an alternative estimator as
Choose $\kappa \in \mathbb{N}$ and $\ell \in [1, \kappa]$ so that Assumption (ref) is satisfied. Let
The following moment inequality restrictions are tested under $H_0 : \theta_0 = \theta$.
where $A$ is a diagonal matrix with diagonal entries $(1,-1)$, $\mu = (\mu_1,\mu_2)' = \operatornamewithlimits{plim}_{n\to\infty} (\hat \mu_{1n}, \hat \mu_{2n})'$, and $\bar{b}=(\bar{b}_{1\ell},\bar{b}_{2\ell})'$. Here, $\bar{b}_{1 \ell}$ and $\bar{b}_{2 \ell}$ represent upper bounds for the finite approximation errors of $\underaccent{\bar}{B}(\theta_0)$ and $\bar{B}(\theta_0)$, respectively.
There is a large body of literature on inference for parameters identified by moment inequality restrictions ({romanoPracticalTwoStepMethod2014}, {romanoPracticalTwoStepMethod2014}; {andrewsInferenceLinearConditional2023}, {andrewsInferenceLinearConditional2023}; {coxSimpleAdaptiveSizeExact2023}, {coxSimpleAdaptiveSizeExact2023}, etc). Among these methods, coxSimpleAdaptiveSizeExact2023 proposed to use a chi-square test with data-adaptive critical values. Their procedure does not require simulation of data or bootstrap, making it well-suited for our application. Following their procedure, we define a quasi-likelihood-ratio test statistic as
where the asymptotic variance, $\hat V_n(\theta)$, can be derived from the asymptotic linear expansion of $\sqrt{n}(\hat \mu_{1n}(\theta),\hat \mu_{2n}(\theta))'$. The critical value is set to $\chi^2_{\hat {df}(\theta)}(1-\alpha)$, where $\hat {df}(\theta)$ denotes the rank of $[A \hat V(\theta)^-]_{\hat J}$. Here, we denote by $[A\hat V(\theta)^-]_{\hat J}$ the submatrix of $A \hat V(\theta)^-$ formed by the rows corresponding to the indices in $\hat J$, and by $\hat J$ the set of indices for the binding inequalities at the solution to (ref). For instance, suppose $\hat V_n(\theta)^{-1}$ exists and only one of the inequalities is binding across all values in $\Theta$. Then, we reject $H_0 : \theta \in \Theta_0$ when $\mathcal{T}_n(\theta) > \chi_1^2(1-\alpha)$. The corresponding confidence region can be constructed as $ C_n = \{ \theta \in\Theta : \mathcal{T}_n(\theta) \le \chi^2_1(1-\alpha) \}. $
To demonstrate the efficacy of the proposed method, we carry out a series of Monte Carlo experiments. We are particularly interested in the size and power properties of our test proposed in Section (ref). We will compare the performance of our method with that of a version of the polynomial estimator to be detailed in this section.
The design of our experiments is described as follows. We generate a random sample $\{ Y_i, X_i \}_{i=1}^n \subseteq \mathbb{R}\times \mathbb{R}$ of size $n = 100,000$ from the augmented isoelastic model. It consists of the equations
where $\eta$ possesses a positive density $f_\eta$ and $X$ is a uniform random variable on $[-1,1]$ independent of $\eta$. The optimization errors are represented by $\varepsilon$, which assumes a triangular density in the window $[K_0, K_1]$. The structural parameters refer to the average income elasticity $\theta_0 = E[\theta(X_i)]$ and the average differential $\omega_0 = E[{\partial \theta(X_i)}/{\partial X_i}]$. We set $\theta_0 = 0.5$ and $\omega_0$ to either $0$ or $0.25$ in the experiment. The case where $\omega_0 = 0$ corresponds to the standard isoelastic model, while $\omega_0 \ne 0$ suggests heterogeneous income responses across different $X$ groups. The policy parameters and the excluded window are set to $(\tau_0, \tau_1, K, K_0, K_1) = (0.00,0.20,2.0,1.7,2.3)$.
To impose the analyticity condition on $f_{\eta} = f_{Y^*(0)}$, we consider two types of distributions: a polynomial density and a Gaussian mixture density, which are referred to as DGP 1 and DGP 2, respectively.\footnote{Details about the construction of these DGPs are relegated to the Appendix.} Figure (ref) exhibits each DGP through the histograms of $Y$ and $Y^*(0)$, where $\theta_0 = 0.5$ and $\omega_0=0$. The left panel for DGP 1 is based on a seventh-degree polynomial fitted to the U.S. taxable income data from 1960 to 1969, which was used in saezTaxpayersBunchKink2010. This corresponds to the counterfactual density estimated using the polynomial strategy with the iterative proportional adjustment. DGP 2 is obtained by fitting a flexible location-scale Gaussian mixture distribution to a mixed skewed generalized error distribution. DGP 2 exemplifies an analytic density not covered by well-known parametric families.
We compare the performance of our inferential method using the generalized polynomial strategy with that of the polynomial estimator, which are referred to as GPS and PE, respectively, hereafter. The approximation order is set to $\ell = 5$ in the GPS method. This choice is seemingly reasonable as it is capable of controlling approximation biases. We adopted critical values disregarding potential approximation bias, which effectively performs inference on the pseudo-true value for $\ell=5$. The selected polynomial degrees range from $7$ to $11$. The size of the extrapolation norm $\chi_{\kappa}^{-1}$ is monitored and bounded as $\chi_{\kappa}^{-1} \approx 10 = n^{1/5}$ for all selected $\kappa$. We trimmed the lower 1% and upper 5% tails of the distribution for each DGP. Such trimming may be necessitated in practice for various reasons: to prevent the counterfactual density from hitting the zero bound, to restrict the support of the data, or to exclude regions where singularities are suspected.
For the inference of $\theta_0 > 0$, where $\omega_0 = 0$ is treated as known, all tests are based on the unweighted bunching moment, $P(Y_{{i}} \in [K_0, K_1])$. Figure (ref) presents the results of the size and power comparison. To make this comparison, we used a version of the polynomial estimator that closely follows the procedure specified in chettyAdjustmentCostsFirm2011. In the Appendix, we show that the iteratively adjusted polynomial estimator can be formulated using the IV (instrumental variable) estimator in the following regression:
where $\mathcal{J} = \{ 1, \ldots, |\mathcal{J}| \}$ represents the set of equispaced bins partitioning the data support, sorted from left to right. The dependent variable is the observed fraction of the data, $f_j = \frac{1}{n}\sum_{i=1}^n \mathbbm{1}\{ {Y_{{i}} \in j} \}$, at bin $j$. The first regressor is an order $k$ polynomial in the center of bin $j$, denoted by $c_j$. This term represents the counterfactual fraction of bin $j$ before the introduction of the kink.
Let ${j}_0+1, \ldots, {j}_1$ be the set of bins located within the window. The second set of regressors is coupled with the coefficients $\beta_l$, ${j}_0 < l\le {j}_1$, which represent the amount of excess fraction relative to the counterfactual at bin $l$. According to the given specification, the observed fraction exceeds the counterfactual fraction by $\beta_j$ if bin $j$ is located within the window. If bin $j$ falls to the right of the window, $f_j$ falls short of the counterfactual fraction by $(\sum_{l={j}_0+1}^{{j}_1} \beta_l)/(n^{-1}\sum_{i} \mathbbm{1}\{ {Y_{{i}} > K_1} \}) f_j$, in proportion to the counterfactual fraction. This assumption ensures that the excess and the missing fractions net to $0 = \sum_{{j}_0 < l\le {j}_1} \beta_l - (\sum_{{j}_0 < l\le {j}_1} \beta_l)/(n^{-1}\sum_i \mathbbm{1}\{ {Y_{{i}} > K_1} \}) \sum_{j > {j}_1} f_j$, thereby the resulting estimate of the counterfactual density satisfies the integral constraint $\sum_{j \in \mathcal{J}} z_k(c_j)'\gamma = 1$ without adjustments.
The endogeneity introduced by the adjustment terms, $- \mathbbm{1}\{ {j > {j}_1} \}/(\sum_{j'> {j}_1}f_{j'})f_j$, which are correlated with the dependent variable, is addressed using an equal number of IVs $$ \left[ \mathbbm{1}\{ {j = {j}_0+1} \},\ldots, \mathbbm{1}\{ {j = {j}_1} \} \right]', $$ which consist of the indicators for the bins within the window. Let $(\hat \beta',\hat \gamma')'$ be the coefficient estimated by the IV estimator using these instruments. We construct the counterfactual density and the fraction of excess bunching using the formulas
where $j^*$ denotes the bin containing the cutoff and $\operatornamewithlimits{mesh}(\mathcal{J})$ represents the width of each bin.\footnote{The original polynomial estimator by chettyAdjustmentCostsFirm2011 used $\hat f = \frac{1}{{j}_1 - {j}_0}\sum_{{j}_0 + 1}^{{j}_1} z_k(c_j)' \gamma / \operatornamewithlimits{mesh}(\mathcal{J})$. We modified it to more closely follow the identification strategy based on the small kink approximation. Under both of our specifications, there were no significant differences between the original and modified estimators and the resulting CIs.} This leads to the final estimate $ \hat \theta = \frac{\hat B / \hat f }{K \log \left( \frac{1-\tau_0}{1-\tau_1}\right)}. $ The standard errors are computed based on the delta method\footnote{When $|\mathcal{J}|$ is sufficiently large, $\varepsilon_j$, $j=1,\ldots, |\mathcal{J}|$, can be treated as approximately independent errors.} and heteroscedasticity-robust standard errors, which allow us to construct an $\alpha$-level confidence intervals about $\theta_0$ as $[\hat \theta \pm z_{1-\alpha/2} \operatornamewithlimits{se}(\hat \theta)]$.
Figure (ref) presents our findings in the isoelastic model. The top panel shows that our tests maintain correct size control across all selected degrees, whereas the power deteriorates as the order increases. In terms of the average CI lengths, degree 11 is approximately 30% wider than degree 7. In contrast, the tests based on the polynomial estimator suffer from negative biases, leading to over-rejection of the correct null hypothesis. A significant source of such biases is suspected to be the lack of an appropriate counterfactual adjustment. Figure (ref)(a) suggests that the PE method likely overestimates the counterfactual density, as the polynomial is fitted to the observed distribution rather than $f_{Y^*(0)}$. Since bunching fundamentally arises from the counterfactual distribution compressing near the cutoff, a proper adjustment should reconstruct it before this compression.
Figure (ref)(b) shows more variability compared to Figure (ref)(a), which is seemingly due to the fact that the counterfactual distribution in DGP 2 exhibits more hunches than in DGP 1. This poses greater difficulty for precise counterfactual estimation.\footnote{ Some preliminary analysis shows that polynomials of degree 7 or 9 are quite poor at fitting the density in DGP 2 globally, especially near the hunch at the right peak. The polynomial of degree 11 exhibits a noticeably improved fit, and the 14th-degree polynomial achieves an almost perfect fit.} When we restrict the data to the 1st through the 60th percentiles, issues with the fitness of the estimated density are mostly resolved. Based on 59% of the sample, we obtain Figure (ref)(c) that exhibits a pattern similar to Figure (ref)(a). Together with Figure (ref)(b), this illustrates that our method has a reduced sieve-approximation bias as the degree increases or the estimated region shrinks. At degree $11$, our method achieves correct size control when 94% of the sample is used. In contrast, the PE methods are consistently oversized under all specifications. Regarding the power properties, the average widths of the CIs, measured by the area between the power curve and the $y = 1$ line, increase with the selected polynomial degree. This implies the bias-variance trade-off in selecting the polynomial degree.
In the augmented model, we test the joint hypothesis $H_0 : \theta_0=0.5, \omega_0 = \omega$ for various values of $\omega$. The test requires at least two pieces of moments. This is facilitated by an additional bunching moment, $E[\exp(X_i) \mathbbm{1}\{ {Y_{{i}} \in [K_0, K_1]} \}]$, which assigns higher weights to those with higher $\theta_i$. The Wald statistic is constructed based on the weighted and the unweighted bunching moments. The critical value is set to $\chi^2_{2}(0.95)$ at the nominal level of 5%, disregarding approximation bias.
Figure (ref) presents the power curves for the joint test under DGPs 1 and 2. In both panels, U-shaped power curves are observed, indicating nearly correct size control and the nontrivial power of the test. The actual size was approximately 7% under DGP 2 when degree $11$ was selected, which appears to be a moderate degree of size distortion. In Figure (ref)(b), we omitted the power curve for degree $7$, as the over-rejection at the true value is already demonstrated in Figure (ref)(b). In light of near size control when degrees $9$ and $11$ are chosen, we reaffirm a trade-off between the size and the power of the test when selecting the polynomial degree.
As an empirical illustration, we revisit saezTaxpayersBunchKink2010, who studied the elasticity of taxable income to marginal tax rate changes assuming the isoelastic model. Specifically, we applied our method to the U.S. tax records data for joint married filers from 1960 to 1969, used for Figure 6A in saezTaxpayersBunchKink2010. Due to limited accessibility to the original data, we instead used the aggregate data attached to the replication files of blomquistBunchingIdentificationTaxable2021. For our purposes, these aggregate data are allegedly sufficient. For details on the institutional setting, we refer the readers to Table 3 in saezTaxpayersBunchKink2010.
By fitting a seventh-degree polynomial to the data outside the range $K \pm \$4000$, where $K = \$20000$, blomquistBunchingIdentificationTaxable2021 found a point estimate of 0.454 using the polynomial estimator by chettyAdjustmentCostsFirm2011. We conduct an empirical analysis of the same dataset based on our method and the polynomial estimator by chettyAdjustmentCostsFirm2011 along with the comprehensive robustness checks across choices of the polynomial degree, the excluded window, and the fitted range of the sample. The approximation order remains $\ell = 5$, as in the previous section, which is seemingly reasonable given the valid size control in the data-based Monte Carlo studies. We apply our testing method and the Wald test based on the polynomial estimator introduced in Section (ref). The degrees are chosen from $\{ 7,9,11 \}$. For the polynomial estimator, we used an equispaced grid with a width of $\$500$ following the analyses in blomquistBunchingIdentificationTaxable2021 and saezTaxpayersBunchKink2010. The excluded window is set either to $K \pm \$4000$ or $K \pm \$5000$, each serving as the main specification and the robustness check. These windows are selected by examining the observed income distribution and identifying the points where the distribution begins to rise and returns to a normal level (see Figure (ref)). The sample is truncated from the 1st to the $u$th percentile, where $u$ ranges over $\{ 60,80,90,95 \}$. The upper boundary $u$ represents the degree of locality in the counterfactual estimation.
The CIs from our analysis are presented in Table (ref). Several crucial observations can be made from this table. First, as the polynomial degree increases, the GPS CIs become less sensitive to sample truncation. This suggests potential issues with the model fitness of the seventh and ninth-degree polynomials over the entire data range. In contrast, the eleventh-degree polynomial provides relatively stable CIs as more data are used, indicating its fitness for the empirical distribution. The GPS CIs appear to converge on the same range as the estimation becomes more local to the window. Consequently, $[0.250, 0.450]$ appears to be the most credible GPS CI robust to specifications.
On the other hand, the PE CIs show greater variability compared to the GPS CIs across the chosen degrees and the truncation levels. Given a fixed level of truncation, the selected polynomial degree affects the final result quite significantly. This makes it more involved to choose an appropriate polynomial degree in the PE method than in the GPS method. The degrees 7 and 9 appear to be poor at fitting the entire range of distribution in light of the sensitivity to the truncated region. Consequently, it seems reasonable to select one of the CIs with degree 11 among the PE CIs. However, the resulting PE CI would significantly differ from the GPS CI.
In Panel B, we repeated the same exercises as above except with the excluded window set to $K \pm \$ 5000$. We observe that the GPS CIs are mildly widened as the excluded window is expanded. This may be attributed to increases in the standard errors. Parallel to Panel A, there is little ambiguity in the resulting GPS CIs when 60% of the sample is used. This substantiates our analysis using the main specification of the window, which is unlikely driven by an inappropriate guess on the error range. However, this observation is unclear for the PE method.
Summing up, we found a CI of $[0.25, 0.45]$ for the taxable income elasticity, with $0.350$ being the most likely value. This result differs from the previous estimates of 0.390 and 0.454, reported in blomquistBunchingIdentificationTaxable2021 using the trapezoidal approximation and polynomial estimator of degree 7, respectively. These differences may well matter for the empirical implementation of the optimal tax rates. Calculations based on saezUsingElasticitiesDerive2001 reveal that these differences in the compensated income elasticity translate to 2.7%p and 6.4%p increases in the revised optimal top tax rates, respectively.\footnote{ According to saezUsingElasticitiesDerive2001, the optimal top rate in the absence of income effect is given by $\tau^* = (1-\bar{g})/(1-\bar{g} + a \theta)$, where $a>0$ denotes the Pareto parameter associated with the tail income distribution and $\bar{g} \ge 0$ represents the social welfare weight assigned to the top bracket tax-payers. In all calculations in this paper, these parameters are calibrated as $a = 2$ and $\bar{g} = 0$, respectively.} If these elasticities apply only to a small neighborhood around the first kink, the 11% and 30% decreases in bunching estimates translate to approximately the same percentage increases in the odds of optimal local tax rates, assuming all else remains equal.
Our analysis reveals that the polynomial estimator may be sensitive to the choice of various user-chosen parameters, such as the polynomial degree, excluded window, or the fitted region. It suggests that the counterfactual correction method reduces the variability of estimation results across polynomial degrees, enhancing their reliability. Our findings demonstrate that the researcher can obtain credible confidence regions for the structural parameter by combining bunching evidence with principled methodological practices.
Bunching in the choice data offers a valuable opportunity to infer a structural parameter by leveraging it as a collective response to quasi-experimental variation in the incentive schedule. Since its initial proposal by saezTaxpayersBunchKink2010 and subsequent development into a formal framework ({chettyAdjustmentCostsFirm2011}, {chettyAdjustmentCostsFirm2011}; {klevenUsingNotchesUncover2013}, {klevenUsingNotchesUncover2013}), it has found a broad range of applications in economic studies. However, despite its growing popularity, the econometric foundations of its theoretical and practical aspects have been lacking until recently.
This paper develops an econometric framework and tools for identification and inference in bunching designs. Our approach can incorporate observably heterogeneous responses into the model, a capability that was only partially permitted in previous methods. Our identification scheme hinges on the analyticity condition of the counterfactual, which does not necessitate an ex-ante parametric family of distributions. We develop a suite of tools for counterfactual estimation and inference, termed the generalized polynomial strategy. The proposed method restores the merits of the traditional polynomial strategy while addressing weaknesses in the widespread practice. In an alternative approach to partial identification, the assumption of analytic counterfactual is relaxed, which could enhance credibility of the procedure.
In this paper, our focus has primarily been on the setting where the counterfactual can be deduced from the observed distribution in the data. In other important applications, such an assumption may not be realistic, for example, due to responses in participation margin, limited data availability, etc. Exploring viable methodologies in these contexts could be a fruitful avenue for future research.