*! version 1.4 2022-03-24

Title

bunching -- uses bunching to partially and point identify the elasticity of a response variable with respect to
changes in the slope of the budget set using different assumptions
on unobserved heterogeneity according to the procedures of Bertanha, McCallum, and Seegert (2022).
bunching is a wrapper function for three other commands: bunchbounds, bunchtobit, and bunchfilter.

Syntax

bunching
depvar [indepvars]
[if] [in]
[weight]
,
kink(#)
s0(#)
s1(#)
m(#)
[ nopic
savingbounds(filename[,replace])
binwidth(#)
grid(numlist)
numiter(#)
savingtobit(filename[,replace])
verbose
deltam(#)
deltap(#)
generate(newvar)
pctobs(#)
polorder(#)
]

Required
depvar must be one dependent variable (the response variable in logs in many applications)

kink(#) is the location of the kink point and must be a real number in the same units as the response variable

s0(#) is a real number; in many applications, it is the log of the slope before the kink point

s1(#) must be a real number that is strictly less than s0; in many applications, it is the log of the slope after the kink point

m(#) is the maximum slope magnitude of the probability density function (PDF) of the unobserved heterogeneity, a strictly positive scalar
(option of bunchbounds)

Optional
indepvars is a varlist of covariates; heterogeneity is a linear function of these covariates and an unobserved error that is normally distributed
conditional on these covariates (option of bunchtobit)

if | in like in any other Stata command, to restrict the working sample

weight follows Stata's weight syntax and only allows frequency weights, fweight

nopic suppresses displaying graphs; the default is to display graphs

savingbounds(filename[,replace]) saves filename.dta with coordinates of the partially-identified set as a function
of the slope magnitude of the heterogeneity PDF; use replace if filename.dta already exists in the working directory (option of bunchbounds)

binwidth(#) is the width of the bins for the histograms; it must be a strictly positive real number;
the default value is half of what is automatically produced by the command histogram

grid(numlist) is a numlist of integers from 1 to 99;
the values in the numlist correspond to percentages of the sample that define symmetric truncation windows around the kink point;
the truncated Tobit model is estimated on each of these samples
and also the full sample so that the number of estimates is always one more than the number of entries in numlist;
for example, if grid(15 82), then bunching estimates the Tobit model three times using 100, 82, and 15 percent of the data around the kink point;
the default value for the numlist is 10(10)90, which provides 10 estimates (option of bunchtobit)

numiter(#) is the maximum number of iterations allowed when maximizing the Tobit log likelihood;
it must be a positive integer and the default is 500 (option of bunchtobit)

savingtobit(filename[,replace]) saves filename.dta with Tobit estimates for each truncation window;
the filename.dta file contains eight variables corresponding to the matrices that the code stores in r();
see below for more details;
use replace if filename.dta already exists in the working directory (option of bunchtobit)

verbose displays detailed output from the Tobit estimation including iterations of maximizing the log likelihood;
non-verbose mode is the default (option of bunchtobit)

deltam(# real) is the distance between the kink point and the lower bound of the support of the friction error to be filtered;
it must be a real number in the same units as the response variable;
if this option is used,
then options deltap and generate must also be specified
(options of bunchfilter)

deltap(# real) is the distance between the kink point and the upper bound of the support of the friction error to be filtered;
it must be a real number in the same units as the response variable;
if this option is used,
then options deltam and generate must also be specified
(options of bunchfilter)

generate(varname) generates the filtered variable with a user-specified name of varname;
if this option is used,
then options deltam and deltap must also be specified
(options of bunchfilter)

pctobs(# real) for better fit, the polynomial regression uses observations in a symmetric window around the kink point that contains pctobs percent of the sample;
default value is 40 ( integer, min = 1, max = 99);
if this option is used,
then options deltam, deltap, and generate must also be specified
(options of bunchfilter)

polorder(# integer) order of polynomial for the filtering regression; default value is 7 (min = 2; max = 7);
if this option is used,
then options deltam, deltap, and generate must also be specified
(options of bunchfilter)

Description

The user enters the name of the response variable, the location of the kink point, the slopes before and after the kink point,
and the maximum slope magnitude of the heterogeneity PDF.
For example, in the income-tax application of Bertanha, McCallum, and Seegert (2022), dollars of taxable income and the dollar value of the kink point
are transformed by taking logs, and the slopes must be input as the log of one minus the marginal tax rates.
You need to have a dataset with the response variable drawn from a mixed continuous-discrete distribution.
The distribution is continuous except for the
kink value, which has a positive mass point.
The proportion of values at the kink value must be positive in your sample.
Check to see if count if y==kink (where y is the response variable) gives you the right number. In case it gives you zero when it should not, check the
value of the kink and if y is type double.
In case your dependent variable is such that the bunching mass is dispersed in a neighborhood around the kink point because of friction errors,
please refer to the last paragraph of this description.

First, bunching calls bunchbounds.
bunchbounds computes the maximum and minimum values of the elasticity that are consistent with the slope restriction on the PDF specified in m,
the observed distribution of the response variable, and values of the PDF of the response variable evaluated at the left and right limits approaching the kink.

Second, the code runs bunchtobit.
bunchtobit estimates multiple mid-censored Tobit regressions using specified sub-samples of the data.
It starts with the entire sample, then it truncates the sample to symmetric windows centered at the kink as specified by the user.
The elasticity estimate is plotted as a function of the percentage of data used in each truncation window.
The code also plots the histogram of the response variable along with the best-fit Tobit distribution for each truncation window.
The user has the option of entering covariates that help explain the unobserved heterogeneity.

In case there are friction errors in the dependent variable and the the bunching mass is dispersed in a neighborhood around the kink point,
the user must specify the following three options together: generate, deltam, and deltap.
bunching will then run bunchfilter before bunchbounds and bunchtobit and generate a variable without
friction errors.
The distribution of the data with error is continuous and its PDF typically exhibits a hump around the location of the mass point.
The distribution without error is
mixed continuous-discrete with one mass point at the kink.
The option generate specifies the name of the new variable without friction errors.
Options deltam and deltap specify the window around the mass point that contains the hump, that is, (kink - deltam, kink + deltap).

(Go up to Syntax)

Examples

The package includes the sample data file "bunching.dta". The income data in logs (i.e., "y")
were simulated using a middle-censored model with an elasticity of 0.5.
The budget constraint has a kink at ln(8) with log slopes of ln(1.3) and ln(0.9), respectively, to the left and right of the kink.
These correspond to tax rates to the left and right of the kink of -30 and 10 percent,
respectively.
The unobserved agent heterogeneity is a linear combination of three
binary covariates (x1, x2, and x3) plus a Gaussian error term.
The data also include randomly generated frequency weights "w".
We introduce nonsharp bunching by adding a random error to the income of bunching
individuals.
Such error term is a truncated normal with support [-ln(0.9), ln(1.1)],
which implies that deltam is 0.1054 and deltap is 0.0953.
The income variable with error is "yfric."

Load the test data file included in the package:

. use bunching

Run bunching with filter:

. bunching yfric x1 x2 x3, k(2.0794) s0(0.2624) s1(-0.1054) m(2) gen(yfilter) deltam(0.1054) deltap(0.0953) pctobs(30) polorder(7) savingtobit(dtobit.dta,replace) savingbounds(dbounds.dta, replace)

Run bunching without filter:

. bunching y x1 x2 x3, k(2.0794) s0(0.2624) s1(-0.1054) m(2) savingtobit(dtobit2.dta, replace) savingbounds(dbounds2, replace)

(Go up to Syntax)

Stored results

bunching stores the following in r():

Matrices

r(tobit_theta_l_hat) intercept of the left-hand side equation

r(tobit_theta_r_hat) intercept of the right-hand side equation

r(tobit_sigma_hat) estimated standard deviation of the error term

r(tobit_perc_obs) percentage of observations selected by truncation window

r(tobit_eps_hat) elasticity estimate

r(tobit_se_hat) standard error of the elasticity estimator

r(tobit_covcol) number of covariates whose coefficients were restricted because of collinearity

r(tobit_flag) dummy that equals one if the log likelihood optimization did not converge

Scalars

r(bounds_e_trap) estimated elasticity using trapezoidal approximation

r(bounds_emin_mmax) lower bound estimate for the elasticity using constant m entered by user

r(bounds_emax_mmax) upper bound estimate for the elasticity using constant m entered by user

r(bounds_M_hat) if the heterogeneity PDF has maximum slope magnitude greater than M_hat, then we cannot rule out PDFs that equal zero inside the bunching interval;
in this case, the upper bound is infinity

r(bounds_M_min_hat) if the heterogeneity PDF has maximum slope magnitude smaller than M_min_hat, then no heterogeneity PDF is consistent with
the observed distribution of income; bounds are not defined

r(bounds_M_data_min) minimum value of the slope magnitude in the estimated PDF of the observed distribution of income (continuous part of the PDF)

r(bounds_M_data_max) maximum value of the slope magnitude in the estimated PDF of the observed distribution of income (continuous part of the PDF)

r(bounds_emin_mhat) lower bound estimate for the elasticity if the choice of m were equal to M_hat;
available only when the choice of m entered by user is bigger than M_hat

r(bounds_emax_mhat) upper bound estimate for the elasticity if the choice of m were equal to M_hat;
available only when the choice of m entered by user is bigger than M_hat

r(tobit_bin_n) number of bins used in histograms of bunchtobit

r(filter_Bhat) estimated bunching mass

r(filter_R2) R-squared of polynomial regression

r(filter_vars_dropped) number of variables dropped out of the polynomial regression in case the initial set of explanatory variables had perfect collinearity

r(binwidth) value of bin width used in histograms of bunchfilter and bunchtobit;
this is not stored if option nopic is stated.

(Go up to Syntax)

Authors

Marinho Bertanha, mbertanha@nd.edu, University of Notre Dame

Andrew McCallum, andrew.h.mccallum@frb.gov, Federal Reserve Board

Nathan Seegert, nathan.seegert@eccles.utah.edu, University of Utah.

Reference

Bertanha, M., McCallum, A., Seegert, N. (2022), "Better Bunching, Nicer Notching". Working paper SSRN 3144539.

Bertanha, M., McCallum, A., Payne, A., Seegert, N. (2022), "Bunching Estimation of Elasticities Using Stata".
Stata Journal, forthcoming.