Nonlinear Binscatter Methods

Nonlinear Binscatter Methods Supplemental Appendix

\thispagestyle{empty}

\pagestyle{plain}

\doublespacing \setcounter{page}{1}

{\onehalfspacing \setcounter{tocdepth}{2} \addtocontents{toc}{\thispagestyle{plain}} }

Introduction

This supplemental appendix is a comprehensive collection of all our new theoretical results for nonlinear binscatter estimators with semi-linear covariate-adjustment and random partitioning. Many of our results contribute to the broader literature on nonparametric estimation and inference, particularly when using series estimators, and are thus of independent interest outside of binned scatter plots. To help place our results in the literature, we include a remark labelled “Improvements over literature” at the end of each technical subsection that discusses in detail the technical improvements presented in that subsection and gives related references.

Here we give a brief summary of this appendix, including pointers to some of the major new results. We proceed as follows. The next subsection lists notation used throughout; further notation is defined throughout Section (ref) and at the outset of Section (ref). Then, Section (ref) describes the setup for nonlinear binscatter methods, including the statistical model, parameters of interest, and assumptions, as well as the (random) partitioning and estimation details. Specifically, Assumption (ref) imposes some basic conditions on the data generating process. Assumption (ref) imposes some technical conditions that characterizes and restricts the statistical model of interest. The loss function specified there is general enough to cover many practically important examples such as mean regression, quantile regression, logit/probit estimation, and Huber regression. Assumption (ref) imposes some mild high-level conditions on the estimation procedure. Assumption (ref) summarizes the key conditions on the partitioning scheme used in our theory. We allow for a large class of random partitions. Importantly, the “convergence” of the random partition (Assumption (ref)(ii)) is not necessary for most of our main theoretical results, thereby allowing for flexible data-driven partitioning methods, including certain recursive adaptive partitioning methods: see Devroye-etal2013_book, zhang2010recursive, and references therein.

Section (ref) presents some preliminary technical lemmas for analyzing nonlinear binscatter (and thus also partitioning-based estimators more broadly). New results include precise non-asymptotic concentration results related to Gram matrices (Lemma (ref)), asymptotic variances (Lemmas (ref) and (ref)), approximation errors (Lemma (ref)), and uniform convergence (Lemma (ref)). Sharp control of these objects is a crucial ingredient for obtaining results under weak conditions, as below.

Section (ref) presents a tight uniform (in $x$) Bahadur representation for nonlinear binscatter (Theorem (ref)). This is our first main result. We allow for random partitions and much weaker rate restrictions (on the tuning parameter $J$) than previously imposed in the literature, in addition to additional controls. The data-dependent partitioning means our series estimator uses random basis functions, and this is entirely new. In terms of tuning parameter rate restrictions, previous results required $J^4/n\to 0$ (up to $\log(n)$ terms) or something stricter, while our restriction is that $J^{\frac{2\nu}{\nu-1}}/n\to 0$ (up to $\log(n)$ terms), with $\nu>2$ denoting the number of finite moments of the “score”, and thus may be substantially weaker. Note that our class of models is often broader than prior work also. Importantly, our results can now allow for piecewise constant binscatters, i.e., with degree $p=0$, which is excluded by prior results in the literature (i.e., for previous technical results there was no sequence $J \to\infty$ such that the bias and variance are simultaneously controlled). In addition, employing our novel uniform Bahadur representation, we can establish the uniform convergence rates of nonlinear binscatter (Corollary (ref)) and variance estimators (Theorem (ref)) under similarly weak restrictions.

Section (ref) studies the pointwise distributional approximation for nonlinear binscatter estimators. These results are omitted from the main paper to save space, but are standard properties of interest in the nonparametrics literature and thus are included for completeness. The main result is Theorem (ref), which establishes pointwise asymptotic Normality for our point estimators, again allowing for random (and possibly “non-convergent”) partitions, and under mild rate restrictions similar to those for the (uniform) Bahadur representation.

Section (ref) presents a new Nagar-type approximate IMSE expansion for nonlinear binscatter estimators with semi-linear covariate-adjustment and random partitions (Theorem (ref)), which has no antecedent in the literature. Our results can be used to design data-driven procedures for selecting IMSE-optimal choices of tuning parameters for nonlinear binscatter. Again, these results are novel in their breadth, the weakness of the assumptions, and the conditions on the partitioning. Here we do require an extra assumption on the partitioning in order to characterize the leading terms in the expansion: intuitively, the random partitioning must “settle” to a population partition so that the leading constants of the expansion can be expressed. For example, sample quantiles converge to population quantiles, so this assumption is satisfied.

Uniform inference is dealt with in the next several sections of this appendix. First, Section (ref) establishes a uniform (in $x$) distributional approximation for nonlinear binscatter estimators. The two main results, which are combined into one in the main text (Theorem 2), are the (conditional) strong approximation in Theorem (ref) and the feasible implementation thereof in Theorem (ref). Again, We allow for a large class of random partitions, a broad class of (possibly) nonlinear and nonsmooth models, and additional controls. Here the partitions do not need to be “convergent” in any sense. These results are obtained under weak assumptions, including in particular mild rate restrictions, as in the case of the uniform Bahadur representation, all of which improves on the literature in various directions as explained in the text below. Finally, Theorem (ref) shows a distributional approximation for the suprema of the $t$-statistic processes in the case of the convergent partition (as in the previous paragraph).

Sections (ref)--(ref) employ the strong approximation results to study uniform inference for various parameters of interest in the specific context of nonlinear binscatter. These results rely on, and inherit the novelty of, Theorems (ref) and (ref). New results include valid uniform confidence bands (Theorem (ref)), consistent hypothesis tests about parametric specification (Theorem (ref)) and tests for shape restrictions (Theorem (ref)). All these results explicitly account for the possibly random partitioning scheme and semi-linear covariate-adjustment with random evaluation points.

Section (ref) discusses implementation details for nonlinear binscatter, including standard error computation, feasible data-driven number of bins selector, and choices of polynomial orders given a fixed number of bins. For a more explicit treatment of the package binsreg per se, see Cattaneo-Crump-Farrell-Feng_2024_Stata and \url{https://nppackages.github.io/binsreg/}.

Finally, Section (ref) contains the proofs for all the technical results in Section (ref).

Notation

See vandevarrt-Wellner_1996_book, Bhatia_2013_book, Gine-Nickl_2016_book, and references therein, for background definitions.

Matrices and Norms. For (column) vectors, $\|\cdot\|$ denotes the Euclidean norm, $\|\cdot\|_1$ denotes the $L_1$ norm, $\|\cdot\|_\infty$ denotes the sup-norm, and $\|\cdot\|_0$ denotes the number of nonzeros. For matrices, $\|\cdot\|$ is the operator matrix norm induced by the $L_2$ norm, and $\|\cdot\|_\infty$ is the matrix norm induced by the supremum norm, i.e., the maximum absolute row sum of a matrix. For a square matrix $\mathbf{A}$, $\lambda_{\max}(\mathbf{A})$ and $\lambda_{\min}(\mathbf{A})$ are the maximum and minimum eigenvalues of $\mathbf{A}$, respectively. $[\mathbf{A}]_{ij}$ denotes the $(i,j)$th entry of a generic matrix $\mathbf{A}$. We will use $\mathcal{S}^{L}$ to denote the unit circle in $\mathbb{R}^L$, i.e., $\|\mathbf{a}\|=1$ for any $\mathbf{a}\in \mathcal{S}^L$. For a real-valued function $g(\cdot)$ defined on a measure space $\mathcal{Z}$, let $\|g\|_{\mathbb{Q}, 2}:=(\int_\mathcal{Z} |g|^2d\mathbb{Q})^{1/2}$ be its $L_2$-norm with respect to the measure $\mathbb{Q}$. In addition, let $\|g\|_\infty=\sup_{z\in\mathcal{\mathcal{Z}}}|g(z)|$ be $L_\infty$-norm of $g(\cdot)$, and if $g$ is a univariate function, let $g^{(v)}(z)=d^vg(z)/dz^v$ be the $v$th derivative for $v\geq 0$.

Asymptotics. For sequences of numbers or random variables, we use $l_n\lesssim m_n$ to denote that $\limsup_n|l_n/m_n|$ is finite, $l_n\lesssim_\mathbb{P} m_n$ or $l_n=O_\mathbb{P}(m_n)$ to denote $\limsup_{\varepsilon\to\infty}\limsup_n\mathbb{P}[|l_n/m_n|\geq\varepsilon]=0$, $l_n=o(m_n)$ implies $l_n/m_n\to 0$, and $l_n=o_\mathbb{P}(m_n)$ implies that $l_n/m_n\to_\mathbb{P} 0$, where $\to_\mathbb{P}$ denotes convergence in probability. Accordingly, we write $l_n\gtrsim m_n$ if $m_n\lesssim l_n$, and $l_n\gtrsim_\mathbb{P} m_n$ if $m_n\lesssim_\mathbb{P} l_n$. $l_n\asymp m_n$ implies that $l_n\lesssim m_n$ and $m_n\lesssim l_n$.

Empirical Process. We employ standard empirical process notation: $\mathbb{E}_n[g(\mathbf{v}_i)]=\frac{1}{n}\sum_{i=1}^ng(\mathbf{v}_i)$, and $\mathbb{G}_n[g(\mathbf{v}_i)]=\frac{1}{\sqrt{n}}\sum_{i=1}^n(g(\mathbf{v}_i)-\mathbb{E}[g(\mathbf{v}_i)])$ for a sequence of random variables $\{\mathbf{v}_i\}_{i=1}^n$. In addition, we employ the notion of covering number extensively in the proofs. Specifically, given a measurable space $(A, \mathcal{A})$ and a suitably measurable class of functions $\mathcal{G}$ mapping $A$ to $\mathbb{R}$ equipped with a measurable envelop function $\bar{G}(z)\geq \sup_{g\in\mathcal{G}}|g(z)|$, the covering number of $N(\mathcal{G}, L_2(\mathbb{Q}), \varepsilon)$ is the minimal number of $L_2(\mathbb{Q})$-balls of radius $\varepsilon$ needed to cover $\mathcal{G}$ for a measure $\mathbb{Q}$. The covering number of $\mathcal{G}$ relative to the envelope is denoted as $N(\mathcal{G}, L_2(\mathbb{Q}), \varepsilon\|\bar{G}\|_{\mathbb{Q},2})$.

Other. $\lceil z \rceil$ outputs the smallest integer no less than $z$ and $a\wedge b=\min\{a,b\}$. “w.p.a. $1$” means “with probability approaching one”.

Setup

Suppose that $(y_i,x_i, \mathbf{w}_i')$, $1\leq i\leq n$, is a random sample where $y_i\in\mathcal{Y}$ is a scalar response variable, $x_i\in\mathcal{X}$ is a scalar covariate, and $\mathbf{w}_i\in\mathcal{W}$ is a vector of additional control variables of dimension $d$. Let $\mathbf{D}= [(y_i, x_i, \mathbf{w}_i')' : i=1,2,\dots,n]$.

For a loss function $\rho(\cdot; \cdot)$ and a strictly monotonic transformation function $\eta(\cdot)$, define

where $\mathcal{M}$ is a space of functions satisfying certain smoothness conditions to be specified later.

This setup is general. For example, consider $\boldsymbol{\gamma}_0=\mathbf{0}$. If $\rho(\cdot;\cdot)$ is a squared loss and $\eta(\cdot)$ is the identity function, $\mu_0(x)$ is the conditional expectation of $y_i$ given $x_i=x$. Let $\mathds{1}(\cdot)$ denote the indicator function. If $\rho(y;\eta)=(q-\mathds{1}(y<\eta))(y-\eta)$ for some $0<q<1$ and $\eta(\cdot)$ is an identity function, then $\mu_0(x)$ is the $q$th conditional quantile of $y_i$ given $x_i=x$. Introducing a transformation function $\eta(\cdot)$ is useful. For instance, it may accommodate logistic regression for binary responses. When $\boldsymbol{\gamma}_0\neq \mathbf{0}$, the parametric and the nonparametric components are additively separable, and thus (ref) becomes a generalized partially linear model.

Binscatter estimators are typically constructed based on a (possibly random) partition of the support of the covariate $x_i$. Specifically, the relevant support of $x_i$ is partitioned into $J$ disjoint intervals, leading to the partitioning scheme $\widehat{\Delta} = \{\widehat{\mathcal{B}}_1, \widehat{\mathcal{B}}_2, \dots, \widehat{\mathcal{B}}_J\}$, where \[ \widehat{\mathcal{B}}_j =

, \] One popular choice in binscatter applications is the quantile-based partition: $\widehat\tau_j=\widehat{F}_X^{-1}(j/J)$ with $\widehat{F}_X(u)=n^{-1}\sum_{i=1}^n\mathds{1}(x_i\leq u)$ the empirical cumulative distribution function and $\widehat{F}_X^{-1}$ its generalized inverse. Our theory is general enough to cover other partitioning schemes satisfying certain regularity conditions specified below. An innovation herein is accounting for the additional randomness from the partition $\widehat\Delta$. The number of bins $J$ plays the role of the tuning parameter for the binscatter method, and is assumed to diverge: $J\to\infty$ as $n\to\infty$ throughout the supplement, unless explicitly stated otherwise.

The piecewise polynomial basis of degree $p$, for some choice of $p=0,1,2,\dots$, is defined as \[\Big[

\Big]' \otimes \Big[

\Big]', \] where $\mathds{1}_\mathcal{A}(x)=\mathds{1}(x\in \mathcal{A})$ and $\otimes$ is the Kronecker product operator. For convenience of later analysis, we use $\widehat{\mathbf{b}}_{p,0}(x)$ to denote a standardized rotated basis, the $j$th element of which is given by \[ \sqrt{J}\times\mathds{1}_{\widehat{\mathcal{B}}_{\bar{j}}}(x)\times\Big(\frac{x-\widehat{\tau}_{\bar{j}-1}}{\hat{h}_{\bar{j}}} \Big)^{j-1-(\bar{j}-1)(p+1)}, \quad j=1, \cdots, (p+1)J, \] where $\bar{j}=\lceil j/(p+1)\rceil$, $\lceil\cdot\rceil$ is the ceiling operator, and $\hat{h}_{\bar{j}}=\widehat\tau_{\bar{j}}-\widehat\tau_{\bar{j}-1}$. Thus, each local polynomial is centered at the start of each bin and scaled by the length of the bin. $\sqrt{J}$ is an additional scaling factor which helps simplify some expressions of our results. The standardized rotated basis $\widehat{\mathbf{b}}_{p,0}(x)$ is equivalent to the original piecewise polynomial basis in the sense that they represent the same (linear) function space.

To impose the restriction that the estimated function is $(s-1)$-times continuously differentiable for $1\leq s \leq p$, we introduce the following basis \[ \widehat{\mathbf{b}}_{p,s}(x)=\Big(\widehat{b}_{p,s,1}(x),\ldots, \widehat{b}_{p,s,K_{p,s}}(x)\Big)' =\widehat{\mathbf{T}}_s\widehat{\mathbf{b}}_{p,0}(x), \qquad K_{p,s}=(p+1)J-s(J-1), \] where $\widehat{\mathbf{T}}_s:=\widehat{\mathbf{T}}_s(\widehat{\Delta})$ is a $K_{p,s}\times(p+1)J$ matrix depending on $\widehat{\Delta}$, which transforms a piecewise polynomial basis into a smoothed binscatter basis. Some useful properties of $\widehat\mathbf{T}_s$ are given in Lemma (ref) in Section (ref), and the explicit representation of $\widehat\mathbf{T}_s$ is available in the proof of Lemma SA-3.2 in \citet*{Cattaneo-Crump-Farrell-Feng_2024_AER}. When $s=0$, we let $\widehat{\mathbf{T}}_0=\mathbf{I}_{(p+1)J}$, the identity matrix of dimension $(p+1)J$. When $s=p$, $\widehat{\mathbf{b}}_{p,s}(x)$ is the well-known $B$-spline basis of order $p+1$ with simple knots, which is $(p-1)$-times continuously differentiable. When $0<s<p$, they can be defined similarly as $B$-splines with knots of certain multiplicities. See Definition 4.1 in Section 4 of Schumaker_2007_book for more details about spline functions. We require $s\leq p$, since if $s=p+1$, $\widehat{\mathbf{b}}_{p,s}(x)$ reduces to a global polynomial basis of degree $p$.

Given a choice of basis, we consider the following generalized binscatter estimator:

where $\widehat{\mathbf{b}}_{p,s}^{(v)}(x)=\frac{d^v}{dx^v}\widehat{\mathbf{b}}_{p,s}(x)$ for some $v\in\mathbb{Z}_+$ such that $v\leq p$. This estimator can be written as:

The representation (ref) allows us to be more general and agnostic about the estimation of $\boldsymbol{\gamma}_0$, and also simplifies some of the proofs. More specifically, our theory requires only a sufficiently fast convergence rate of $\widehat{\boldsymbol{\gamma}}$ (see Assumption (ref) below), which in nonlinear estimation models cases can be justified in different ways, e.g., joint estimation, backfitting, profiling, and split-sampling, among other possibilities. Our software implementation \citep*{Cattaneo-Crump-Farrell-Feng_2024_Stata} relies on joint estimation, as done by the default base estimation packages in Python, R, and Stata.

In this supplement, we focus on estimation and inference of the following three parameters:

where $\mathsf{w}$ is a user-chosen evaluation point of the control variables, and thus these parameters are viewed as functions of $x$ only in our theory. Nevertheless, all our results are readily applied to other linear or nonlinear transformations of $\mu_0(x)$, such as the higher-order derivatives $\frac{\partial^v}{\partial x^v}\eta(\mu_0(x)+\mathsf{w}'\boldsymbol{\gamma}_0)$. Given the binscatter estimates $\widehat{\mu}_{p,s}(x)$ and $\widehat\boldsymbol{\gamma}$ in (ref), the estimators of the three parameters defined above are given by

respectively, for some consistent estimate $\widehat\mathsf{w}$ (non-random or generated based on $\{\mathbf{w}_i\}_{i=1}^n$) of the evaluation point $\mathsf{w}$. As a reminder, we need to require $p\geq v$ to get $\widehat\mu_{p,s}^{(v)}(x)$, $p\geq 0$ to get $\widehat\vartheta_{p,s}(x,\widehat\mathsf{w})$, and $p\geq 1$ to get $\widehat\zeta_{p,s}(x,\widehat\mathsf{w})$.

Recall that in the main text we always set $s=p$ and omit the dependence of estimators on $s$. Thus, $\widehat{\mu}_p^{(v)}(x)=\widehat{\mu}_{p,p}^{(v)}(x)$, $\widehat{\vartheta}_p(x,\widehat\mathsf{w})= \widehat{\vartheta}_{p,p}(x,\widehat\mathsf{w})$, and $\widehat{\zeta}_p(x,\widehat\mathsf{w})= \widehat{\zeta}_{p,p}(x,\widehat\mathsf{w})$. In this supplement, however, all our results hold for a general choice of the degree and the smoothness of the basis. For ease of notation, the subscripts $p$ and $s$ of the above estimators are dropped hereafter: \[ \widehat{\mu}^{(v)}(x):=\widehat{\mu}_{p,s}^{(v)}(x), \quad \widehat{\vartheta}(x,\widehat\mathsf{w}):=\widehat{\vartheta}_{p,s}(x,\widehat\mathsf{w}),\quad \text{and}\quad \widehat{\zeta}(x,\widehat\mathsf{w}):=\widehat{\zeta}_{p,s}(x,\widehat\mathsf{w}). \]

We first assume the following basic conditions on the data generating process.

Next, we impose several technical conditions related to the statistical model of interest.

Our next assumption imposes mild high-level conditions on the estimator $\widehat{\boldsymbol{\gamma}}$ of the coefficient vector $\boldsymbol{\gamma}_0$, the estimator $\widehat\mathsf{w}$ of the evaluation point $\mathsf{w}$ for control variables, and the estimator of the function $\Psi_1$ defined previously in Assumption (ref)(iv).

Note that $\Upsilon(x,\mathbf{w})=\Psi_1(x, \mathbf{w};\eta(\mu_0(x)+\mathbf{w}'\boldsymbol{\gamma}_0))$ in the main paper to streamline the presentation. Part (i) is a mild condition on the convergence of $\widehat{\boldsymbol{\gamma}}$ and $\widehat\mathsf{w}$. Part (ii) is a high-level condition that ensures we have a valid feasible estimator of the Gram matrix ($\bar\mathbf{Q}$ or $\mathbf{Q}_0$ defined at the outset of Section (ref) below). Note that the convergence rate of $\eta^{(1)}(\widehat{\mu}(x_i)+\mathbf{w}_i'\widehat{\boldsymbol{\gamma}})$ can be deduced from Corollary (ref) below. Thus, part (ii) can be largely viewed as a restriction on $\widehat{\Psi}_1$ only. Note that $\widehat{\Psi}_1$ does not have to be consistent for $\Psi_1$ in any sense; it suffices that the estimator $\mathbb{E}_n[\widehat{\mathbf{b}}_{p,s}(x_i)\widehat{\mathbf{b}}_{p,s}(x_i)' \widehat{\varkappa}(x_i, \mathbf{w}_i)]$ based on $\widehat{\Psi}_1$ as a whole is consistent. See Section (ref) for several examples of the estimator $\widehat{\Psi}_1$.

Partitions

We need some regularity conditions on the partitioning scheme, which can be verified in a case-by-case basis. We first define a family of “quasi-uniform” partitions for some absolute constant $C>0$:

where $h_j(\Delta)$ denotes the length of the $j$th bin in the partition $\Delta$. Roughly speaking, (ref) says that the bins in any $\Delta\in\Pi_C$ do not differ too much in length. Also, let $\mathbf{X} = [x_1, \ldots, x_n]'$, $\mathbf{W}=[\mathbf{w}_1,\cdots, \mathbf{w}_n]'$ and $\mathbf{Y}=[y_1, \cdots, y_n]'$.

Part (i) is the key condition for our main results and will be imposed throughout. First, it requires that the possibly random partition $\widehat\Delta$ be independent of the outcome $\mathbf{Y}$ given the covariates $(\mathbf{X},\mathbf{W})$. This conditional independence assumption is trivially satisfied if $\widehat\Delta$ is deterministic (e.g., equally-spaced partition) or depends on $\mathbf{X}$ and $\mathbf{W}$ only (e.g., quantile-spaced partition based on $\mathbf{X}$). It also holds if a sample splitting scheme is used: a subsample (including the information about the outcome) is used for constructing the partition, and the other is employed to construct the binscatter estimator (so that $\widehat\Delta$ is independent of the data $(\mathbf{X}, \mathbf{W}, \mathbf{Y})$). Second, $\widehat\Delta$ is required to be “quasi-uniform” with large probability. It is trivially true for equally-spaced partitions and can be verified for quantile-spaced partitions under the mild conditions on the covariates density imposed before (see Lemma (ref)). However, this condition may be too restrictive for other modern machine-learning-based partitioning methods, in which case some additional regularization may be necessary to recover the quasi-uniformity property.

Part (ii) requires that the random partition $\widehat\Delta$ “stabilizes” to a fixed one in large samples. This is true if the partition is non-deterministic or generated by sample quantiles (since sample quantiles converge to population quantiles), but more generally, it is not always possible. Fortunately, this “convergence” requirement is not necessary for most of our main results (except Theorem (ref) and Theorem (ref)). Thus, we will always make it clear if part (ii) of Assumption (ref) is imposed.

Given the random partition $\widehat{\Delta}$, we use the notation $\mathbb{E}_{\widehat{\Delta}}[\cdot]$ to denote the expectation operator with the partition $\widehat{\Delta}$ viewed as fixed. To further simplify notation, let $\hat{h}_j=\hat{\tau}_j-\hat{\tau}_{j-1}$ be the width of the $j$th bin $\widehat{\mathcal{B}}_j$, and when the “limiting” partition $\Delta_0=\{\mathcal{B}_1, \cdots, \mathcal{B}_J\}$ is defined (Assumption (ref)(ii) holds), let $h_j$ be the width of $\mathcal{B}_j$. Analogously to $\widehat{\mathbf{b}}_{p,s}(x)$, $\mathbf{b}_{p,s}(x)$ denotes the binscatter basis of degree $p$ that is $(s-1)$-times continuously differentiable and is constructed based on the nonrandom partition $\Delta_0$. We sometimes write $\mathbf{b}_{p,s}(x;\Delta)=(b_{p,s,1}(x;\Delta), \ldots, b_{p,s,K_{p,s}}(x;\Delta))'$ to emphasize a binscatter basis is constructed based on a particular partition $\Delta$. Therefore, $\widehat{\mathbf{b}}_{p,s}(x)=\mathbf{b}_{p,s}(x;\widehat{\Delta})$ and $\mathbf{b}_{p,s}(x)=\mathbf{b}_{p,s}(x;\Delta_0)$. Accordingly, we use $\mathbf{T}_s$ to denote the transformation matrix based on the non-random partition $\Delta_0$ (which transforms $\mathbf{b}_{p,0}(x)$ to $\mathbf{b}_{p,s}(x)$).

Main Results

We introduce the following quantities that will be extensively used throughout the supplement:

Recall that in the main text we always set $s=p$ and omit the dependence on $s$ whenever there is no confusion. Thus,

In this supplement, however, all our results hold for a general choice of the degree and the smoothness of the basis. For ease of notation, the subscripts $p$ and $s$ of the above quantities are dropped hereafter:

In addition, given the family $\Pi_C$ of the quasi-uniform partitions defined in (ref), for any $\Delta\in \Pi$, we let $\boldsymbol{\beta}_0(\Delta)\in\mathbb{R}^{K_{p,s}}$ be any vector such that for every $v\leq p$,

Let $r_{0,v}(x;\Delta)=\mu_0^{(v)}(x)-\mathbf{b}_{p,s}^{(v)}(x;\Delta)'\boldsymbol{\beta}_0(\Delta)$ denote the corresponding approximation error. Accordingly, given the random partition $\widehat\Delta$, we let $\widehat\boldsymbol{\beta}_0:=\boldsymbol{\beta}_0(\widehat\Delta)$, and $\widehat{r}_{0,v}(x)=\mu_0^{(v)}(x)-\widehat\mathbf{b}_{p,s}^{(v)}(x)'\widehat\boldsymbol{\beta}_0$ denote the corresponding approximation error. The existence of such vectors is guaranteed by Assumptions (ref) and (ref)(v), and is verified in Lemma (ref) in Section (ref).

Preliminary Lemmas

The next lemma shows that the limiting variance is bounded from above and below.

The next lemma gives a bound on the variance component of the nonlinear binscatter estimator.

Our result holds without imposing any smoothness restrictions on the estimation space. Specifically, the estimation procedure (ref) searches for solutions in $\mathbb{R}^{K_{p,s}}$, leading to an estimation space $\{\widehat{\mathbf{b}}_{p,s}(x)'\boldsymbol{\beta}: \boldsymbol{\beta}\in\mathbb{R}^{K_{p,s}}\}$. In contrast, many studies of series (or sieve) methods restrict the functions in the estimation space to satisfy certain smoothness conditions, e.g., Lipschitz continuity, to derive the uniform consistency. See, for example, Chernozhukov-Imbens-Newey_2007_JoE and references therein.

Bahadur Representation

The following corollary is an immediate result of Lemma (ref) and Theorem (ref), and hence its proof is omitted.

The next theorem shows that the proposed variance estimator is consistent.

Pointwise Inference

Starting from this section, we consider statistical inference on $\mu_0^{(v)}(x)$, $\vartheta_0(x,\mathsf{w})$ and $\zeta_0(x,\mathsf{w})$ based on the following Studentized $t$-statistics:

The next theorem shows the pointwise asymptotic normality of the binscatter estimators.

Integrated Mean Squared Error

In this section we give a Nagar-type approximate IMSE expansion for each of the three estimators $\widehat\mu^{(v)}(x)$, $\widehat\vartheta(x,\widehat\mathsf{w})$ and $\widehat\zeta(x,\widehat\mathsf{w})$, with explicit characterization of the leading constants. Define

where $\mathscr{E}_m(\cdot)$ is the $m$th Bernoulli polynomial for each $m\in\mathbb{Z}_+$, $\tau_x^{\mathtt{L}}$ is the start of the interval in the non-random partition $\Delta_0$ containing $x$ and $h_x$ denotes its length.

In general, $\mathscr{B}_n(p,s,v)\gtrsim 1$ (see Remark SA-3.7 in Cattaneo-Crump-Farrell-Feng_2024_AER), and thus the above theorem implies that the (approximate) IMSE-optimal number of bins satisfies that $J_{\mathtt{AIMSE}}\asymp n^{\frac{1}{2p+3}}$. Relying on the IMSE expansion in Theorem (ref), one may design a data-driven procedure to select the IMSE-optimal number of bins for nonlinear binscatter-based M-estimators.

Uniform Inference

Recall that $(a_n: n\geq 1)$ is a sequence of non-vanishing constants. We will first show that the (feasible) Studentized $t$-statistic processes $T_{\mu^{(v)},p}(\cdot)$, $T_{\vartheta,p}(\cdot)$ and $T_{\zeta,p}(\cdot)$ can be approximated by Gaussian processes in a proper sense at certain rate.

The approximating processes $\bar{Z}_{\mu^{(v)},p}(\cdot)$, $\bar{Z}_{\vartheta,p}(\cdot)$ and $\bar{Z}_{\zeta,p}(\cdot)$ are Gaussian processes conditional on $\mathbf{X}$, $\mathbf{W}$ and $\widehat{\Delta}$, and $Z_{\mu^{(v)},p}(\cdot)$, $Z_{\vartheta,p}(\cdot)$ and $Z_{\zeta,p}(\cdot)$ are Gaussian processes conditional on $\widehat{\Delta}$ by construction. In practice, one can replace all unknowns in $\bar{Z}_{\mu^{(v)},p}(\cdot)$, $\bar{Z}_{\vartheta, p}(\cdot)$ and $\bar{Z}_{\zeta, p}(\cdot)$ (or $Z_{\mu^{(v)},p}(\cdot)$, $Z_{\vartheta, p}(\cdot)$ and $Z_{\zeta, p}(\cdot)$) by their sample analogues, and then construct the following feasible (conditional) Gaussian processes:

where $\mathbf{N}_{K_{p,s}}^\star$ denotes a $K_{p,s}$-dimensional standard normal vector independent of the data $\mathbf{D}$ and the partition $\widehat{\Delta}$.

For ease of presentation, we will always require a fast convergence rate of $\widehat\mathsf{w}$ hereafter: $\|\widehat\mathsf{w}-\mathsf{w}\|=o_\mathbb{P}(a_n^{-1}\sqrt{J/n})$. Nevertheless, note that as shown in Theorem (ref), such a rate restriction on $\widehat\mathsf{w}$ can be different for inference of $\vartheta_0(x,\mathsf{w})$ and $\zeta_0(x,\mathsf{w})$ and are unnecessary for inference of $\mu_0^{(v)}(x)$.

Theorems (ref) and (ref) offer a way to approximate the distribution of the whole $t$-statistic process based on $\widehat\mu^{(v)}(\cdot)$, $\widehat\vartheta(\cdot, \widehat\mathsf{w})$ or $\widehat{\zeta}(\cdot,\widehat\mathsf{w})$. A direct application of these results is the distributional approximations to the suprema of these $t$-statistic processes.

Confidence Bands

Let

be confidence bands for $\mu_0^{(v)}(\cdot)$, $\vartheta_0(\cdot,\mathsf{w})$ and $\zeta_0(\cdot,\mathsf{w})$ respectively, where $\mathfrak{c}_{\mu^{(v)}}$, $\mathfrak{c}_\vartheta$ and $\mathfrak{c}_\zeta$ are corresponding critical values to be specified. Recall that $\mathsf{w}$ here is taken as a fixed evaluation point for the control variables, and these bands are constructed based on a certain choice of $J$ and the $p$th-order binscatter basis. Using the previous results, we have the following theorem.

Parametric Specification Tests

As another application, we can test parametric specifications of $\mu_0^{(v)}(x)$, $\vartheta_0(x,\mathsf{w})$ and $\zeta_0(x,\mathsf{w})$. We introduce the following tests:

where $m(x;\boldsymbol{\theta})$ is some known function depending on some finite dimensional parameter $\boldsymbol{\theta}$. This testing problem can be viewed as a two-sided test where the equality between two functions holds uniformly over $x\in\mathcal{X}$. In this case, we introduce $\widetilde{\boldsymbol{\theta}}$ and $\widetilde{\boldsymbol{\gamma}}$ as consistent estimators of $\boldsymbol{\theta}$ and $\boldsymbol{\gamma}_0$ under $\dot{\mathsf{H}}_0^{\mu^{(v)}}$. Then we rely on the following test statistic: \[ \dot{T}_{\mu^{(v)},p}(x) :=\frac{\widehat{\mu}^{(v)}(x)-m^{(v)}(x;\widetilde{\boldsymbol{\theta}})} {\sqrt{\widehat{\Omega}_{\mu^{(v)}}(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}}|\dot{T}_{\mu^{(v)},p}(x)|>\mathfrak{c}_{\mu^{(v)}}$ for some critical value $\mathfrak{c}_{\mu^{(v)}}$.

Similarly, to test the specification of $\vartheta_0(x,\mathsf{w})$, we introduce

where $M(x,\mathsf{w};\boldsymbol{\theta},\boldsymbol{\gamma}_0)=\eta(m(x;\boldsymbol{\theta})+\mathsf{w}'\boldsymbol{\gamma}_0)$. We rely on the following test statistic: \[ \dot{T}_{\vartheta,p}(x) :=\frac{\widehat{\vartheta}(x,\widehat\mathsf{w})-M(x,\widehat\mathsf{w};\widetilde{\boldsymbol{\theta}},\widetilde\boldsymbol{\gamma})} {\sqrt{\widehat{\Omega}_\vartheta(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}}|\dot{T}_{\vartheta,p}(x)|>\mathfrak{c}_\vartheta$ for some critical value $\mathfrak{c}_\vartheta$.

To test the specification of $\zeta_0(x,\mathsf{w})$, we introduce

where $M^{(1)}(x,\mathsf{w};\boldsymbol{\theta},\boldsymbol{\gamma}_0):=\eta^{(1)}(m(x;\boldsymbol{\theta})+\mathsf{w}'\boldsymbol{\gamma}_0)m^{(1)}(x;\boldsymbol{\theta})$. We rely on the following test statistic: \[ \dot{T}_{\zeta,p}(x) :=\frac{\widehat{\zeta}(x,\widehat\mathsf{w})-M^{(1)}(x,\widehat\mathsf{w};\widetilde{\boldsymbol{\theta}},\widetilde\boldsymbol{\gamma})} {\sqrt{\widehat{\Omega}_\zeta(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}}|\dot{T}_{\zeta,p}(x)|>\mathfrak{c}_\zeta$ for some critical value $\mathfrak{c}_\zeta$.

Shape Restriction Tests

The third application of our results is to test certain shape restrictions on $\mu_0^{(v)}(x)$, $\vartheta_0(x,\mathsf{w})$ and $\zeta_0(x,\mathsf{w})$. To be specific, consider the following problem: \[

\]

This testing problem can be viewed as a one-sided test where the inequality holds uniformly over $x\in\mathcal{X}$. Importantly, it should be noted that under both $\ddot{\mathsf{H}}_0^{\mu^{(v)}}$ and $\ddot{\mathsf{H}}_\text{A}^{\mu^{(v)}}$, we fix $\bar{\boldsymbol{\theta}}$ and $\bar{\boldsymbol{\gamma}}$ to be the same values in the parameter space. In such a case, we introduce $\widetilde{\boldsymbol{\theta}}$ and $\widetilde{\boldsymbol{\gamma}}$ as consistent estimators of $\bar{\boldsymbol{\theta}}$ and $\bar{\boldsymbol{\gamma}}$ under both $\ddot{\mathsf{H}}_0^{\mu^{(v)}}$ and $\ddot{\mathsf{H}}_\text{A}^{\mu^{(v)}}$. Then we will rely on the following test statistic: \[ \ddot{T}_{\mu^{(v)},p}(x):= \frac{\widehat{\mu}^{(v)}(x)-m^{(v)}(x;\widetilde{\boldsymbol{\theta}})} {\sqrt{\widehat{\Omega}_{\mu^{(v)}}(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}} \ddot{T}_{\mu^{(v)},p}(x)>\mathfrak{c}_{\mu^{(v)}}$ for some critical value $\mathfrak{c}_{\mu^{(v)}}$.

Similarly, define the test for the shape of $\vartheta_0(x,\mathsf{w})$: \[

\] We will rely on the following test statistic: \[ \ddot{T}_{\vartheta,p}(x):= \frac{\widehat{\vartheta}(x,\widehat\mathsf{w})-M(x,\widehat\mathsf{w};\widetilde{\boldsymbol{\theta}},\widetilde{\boldsymbol{\gamma}})} {\sqrt{\widehat{\Omega}_\vartheta(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}} \ddot{T}_{\vartheta,p}(x)>\mathfrak{c}_\vartheta$ for some critical value $\mathfrak{c}_\vartheta$.

Also, define the test for the shape of $\zeta_0(x,\mathsf{w})$: \[

\] We will rely on the following test statistic: \[ \ddot{T}_{\zeta,p}(x):= \frac{\widehat{\zeta}(x,\widehat\mathsf{w})-M^{(1)}(x,\widehat\mathsf{w};\widetilde{\boldsymbol{\theta}},\widetilde{\boldsymbol{\gamma}})} {\sqrt{\widehat{\Omega}_\zeta(x)/n}}. \] The null hypothesis is rejected if $\sup_{x\in\mathcal{X}} \ddot{T}_{\zeta,p}(x)>\mathfrak{c}_\zeta$ for some critical value $\mathfrak{c}_\zeta$.

The following theorem characterizes the size and power of such tests.

Implementation Details

Standard Error Computation

In Section (ref), we have given the variance formulas $\widehat\Omega_{\mu^{(v)}}(x)$, $\widehat\Omega_\vartheta(x)$ and $\widehat\Omega_\zeta(x)$ that can be used to obtain the standard errors of $\widehat{\mu}^{(v)}(x)$, $\widehat{\vartheta}(x,\widehat\mathsf{w})$ and $\widehat{\zeta}(x,\widehat\mathsf{w})$. Recall that the formula for the estimator $\widehat{\boldsymbol{\Sigma}}$ of $\boldsymbol{\Sigma}_0$ is $$ \widehat{\boldsymbol{\Sigma}}= \mathbb{E}_n\Big[\widehat{\mathbf{b}}_{p,s}(x_i)\widehat{\mathbf{b}}_{p,s}(x_i)' \psi(y_i, \widehat{\eta}_i)^2\eta^{(1)}(\widehat{\mu}(x_i)+\mathbf{w}_i'\widehat{\boldsymbol{\gamma}})^2\Big]. $$ It only relies on known or estimable quantities such as the derivative of the loss function $\psi(\cdot)$, the derivative of the inverse link function $\eta^{(1)}(\cdot)$, the residual $\widehat{\epsilon}_i$ and the binscatter estimates $\widehat{\mu}(\cdot)$ and $\widehat\boldsymbol{\gamma}$. Thus, $\widehat{\boldsymbol{\Sigma}}$ and other types of heteroskedasticity-robust “meat” matrix estimators can be easily constructed using the data. Then, it remains to obtain an estimator $\widehat{\mathbf{Q}}$ of $\bar\mathbf{Q}$ (or $\mathbf{Q}_0$), which in general relies on an estimator $\widehat{\Psi}_1(\cdot)$ of $\Psi_1(\cdot)$ and can be constructed in a case-by-case basis. In the following we discuss several examples.

Example 1 (Least Squares Regression). For least squares regression, the loss function $\rho(y; \eta)=\frac{1}{2}(y-\eta)^2$ and the (inverse) link function $\eta(\theta)=\theta$. Therefore, $\psi(y_i,\eta_i)=-\epsilon_i$ and $\eta_{i,1}=1$. Thus, the formula for $\widehat{\mathbf{Q}}$ given in Section (ref) reduces to $\mathbb{E}_n[\widehat\mathbf{b}_{p,s}(x_i)\widehat\mathbf{b}_{p,s}(x_i)']$, which is immediately feasible in practice.

Example 2 (Logistic Regression). For logistic regression, the loss function is given by the corresponding likelihood function, i.e., $-\rho(y;\eta)=y\log \eta+(1-y)\log (1-\eta)$, and the inverse link is given by the logistic function $\eta(\theta)=\frac{e^\theta}{1+e^\theta}$. Accordingly, an estimator of $\mathbf{Q}_0$ is given by \[ \widehat{\mathbf{Q}}=\mathbb{E}_n\Big[\widehat{\mathbf{b}}_{p,s}(x_i)\widehat{\mathbf{b}}_{p,s}(x_i)'\widehat{\eta}_i(1-\widehat{\eta}_i)\Big], \quad \widehat{\eta}_i=\eta(\widehat{\mu}(x_i)+\mathbf{w}_i'\widehat{\boldsymbol{\gamma}}). \]

Example 4 (Quantile Regression). For quantile regression, $\rho(y;\eta)=(q-\mathds{1}(y<\eta))(y-\eta)$ for some $q\in(0,1)$ and $\eta(\theta)=\theta$. Accordingly, $\psi(y_i,\eta_i)=\mathds{1}(\epsilon_i<0)-q$, and one needs to estimate \[ \mathbf{Q}_0=\mathbb{E}\Big[\mathbf{b}_{p,s}(x_i)\mathbf{b}_{p,s}(x_i)'f_{Y|XW}(\mu_0(x_i)+\mathbf{w}_i'\boldsymbol{\gamma}_0|x_i,\mathbf{w}_i)\Big]. \] The key is to estimate the conditional density $f_{Y|XW}(\cdot|x_i,\mathbf{w}_i)$ evaluated at the conditional quantile of interest $(\mu_0(x_i)+\mathbf{w}_i'\boldsymbol{\gamma}_0)$, whose reciprocal is termed “sparsity function” in the literature. Many different methods have been proposed. For example, the sparsity function is simply the derivative of the conditional quantile function with respect to the quantile, which can be estimated by using the difference quotient of the estimated conditional quantile function. Alternatively, $\mathbf{Q}_0$ can be viewed as a matrix-weighted density function, and one can construct a corresponding estimator based on kernel density estimation ideas. In addition, one can use bootstrapping methods to estimate the variance, avoiding the technical difficulty of estimating the sparsity function. See Section 3.4 and Section 3.9 of Koenker_2005_book for more discussion of variance estimation for quantile regression.

Number of Bins Selector

We discuss the implementation details for data-driven selection of the number of bins, based on the approximate integrated mean squared error expansion in Theorem (ref).

We offer two procedures for estimating the bias and variance constants, and once these estimates ($\widehat{\mathscr{B}}_n(p,s,v)$ and $\widehat{\mathscr{V}}_n(p,s,v)$) are available, the estimated optimal $J$ is \[ \widehat{J}_{\mathtt{IMSE}}=\bigg\lceil \bigg(\frac{2(p-v+1)\widehat{\mathscr{B}}_n(p,s,v)} {(1+2v)\widehat{\mathscr{V}}_n(p,s,v)}\bigg)^{\frac{1}{2p+3}} n^{\frac{1}{2p+3}} \bigg\rceil. \] We always let $\omega(x)=f_X(x)$ as weighting function for concreteness.

Rule-of-thumb Selector

A rule-of-thumb choice of $J$ can be obtained based on Corollary SA-3.2 in Cattaneo-Crump-Farrell-Feng_2024_AER, which gives an explicit characterization of the variance and bias constants for least squares binscatter using piecewise polynomials ($s=0$).

Specifically, the variance constant $\mathscr{V}(p,0,v)$ is estimated by \[ \widehat{\mathscr{V}}(p,0,v)=\operatorname*{trace}\Big\{\Big(\int_{0}^{1}\bm{\varphi}(z)\bm{\varphi}(z)'dz\Big)^{-1}\int_0^1\bm{\varphi}^{(v)}(z)\bm{\varphi}^{(v)}(z)'dz\Big\}\times \frac{1}{n}\sum_{i=1}^{n}\widehat{\sigma}^2(x_i,\mathbf{w}_i)\widehat{f}_X(x_i)^{2v} \] where $\bm{\varphi}(z)=(1, z, \ldots, z^p)'$, $\widehat{\sigma}^2(x_i,\mathbf{w}_i)$ is some estimate of the conditional variance $\mathbb{V}[y_i|x_i,\mathbf{w}_i]$ and $\widehat{f}_X(x_i)$ is some estimate of the density $f_X(x_i)$. On the other hand, the bias constant $\mathscr{B}(p,0,v)$ is estimated by \[ \widehat{\mathscr{B}}(p,0,v)=\frac{\int_{0}^{1}[\mathscr{B}_{p+1-v}(z)]^2 dz}{((p+1-v)!)^2}\times \frac{1}{n}\sum_{i=1}^{n}\frac{[\widehat{\mu}^{(p+1)}(x_i)]^2}{\widehat{f}_X(x_i)^{2p+2-2v}}. \] where $\mathscr{B}_p(z)=(-1)^p\sum_{k=0}^p\binom{p}{k}\binom{p+k}{k}(-z)^k/\binom{2p}{p}$ for each $p\in\mathbb{Z}_+$ and $\widehat\mu^{(p+1)}(x_i)$ is some preliminary estimate of $\mu_0^{(p+1)}(x_i)$. The details about getting the estimates $\widehat\sigma^2(x_i,\mathbf{w}_i)$, $\widehat{f}_X(x_i)$ and $\widehat\mu^{(p+1)}(x_i)$ can be found in Section SA-4.1 in Cattaneo-Crump-Farrell-Feng_2024_AER.

This procedure still yields a choice of $J$ with the correct rate, though the constant approximations are inconsistent for general loss.

Direct-plug-in Selector

The direct-plug-in selector is implemented based on nonlinear binscatter estimators, which applies to any user-specified $p$, $s$ and $v$. It requires a preliminary choice of $J$, for which the rule-of-thumb selector previously described can be used.

More generally, suppose that a preliminary choice $J_{\mathtt{pre}}$ is given, and then a binscatter basis $\widehat{\mathbf{b}}_{p,s}(x)$ (of order $p$) can be constructed immediately on the preliminary partition. Implementing a nonlinear binscatter estimation using this basis and partitioning, we can obtain the variance constant estimate using the variance matrix estimators discussed in Section (ref).

Regarding the bias constant, the key unknown in the expression of the leading approximation error $r_{0,v}^\star(x)$ in Theorem (ref) is $\mu_0^{(p+1)}(x)$, which can be estimated by implementing a nonlinear binscatter estimation of order $p+1$ (with the preliminary partition unchanged). Also, an estimate of $f_X(x_i)^{-1}$ in $r_{0,v}^\star(x)$ is $J\hat{h}_{x_i}$ where $\hat{h}_{x_i}$ denotes the length of the interval in $\widehat\Delta$ containing $x_i$. All other quantities in the expression of $\mathscr{B}(p,s,v)$ can be replaced by their sample analogues. Then, a bias constant estimate is available.

By this construction, the direct-plug-in selector employs the correct rate and consistent constant approximations for any nonlinear binscatter with any choice of $p$, $s$ and $v$.

Fixed $J$ and choice of polynomial order

Our main theory treats $J$ as diverging with the sample size. This reflects the standard approach wherein a researcher selects $p$ and $s$ in advance (often $s=p=0$ or $s=p=3$) and then selects $J$ given the data. The partition must get finer to remove the nonparametric smoothing bias in estimating the function $\mu_0(x)$ (and along with it, $\vartheta_0(x,\mathsf{w})$ or $\zeta_0(x,\mathsf{w})$). Correct recovery (either for estimation or visualization) of these functions is the primary use of binscatter. However, researchers sometimes prefer to pre-specify a fixed $J = {\tt J}$, and we also discuss implementation and interpretation of binscatter in this case.

Instead of modeling $J$ as diverging and searching for the optimal choice, a researcher may desire a fixed (often small and round) number of $J$, which we denote by ${\tt J}$. This is done either to make the estimate more visually appealing or because the results can be directly interpreted. In this case, instead of recovering the functions $\mu_0^{(v)}(x)$, $\vartheta_0(x,\mathsf{w})$, and $\zeta_0(x,\mathsf{w})$, the binscatter is interpreted as estimating their coarsened versions: the distribution of $y_i$ conditional on $x_i$ lying in a (fixed) bin, rather than at a single point. For some ${\tt J}$, this remains interpretable and all our inference results apply to this case. For example, in our application we can take ${\tt J} = 10$ and study the distribution of uninsured rate for each decile of income. The confidence bands then become pointwise confidence intervals that are simultaneously valid. For example, this could be used to examine inequality in health care access by asking if median uninsured rates are statistically different between the top and bottom decile.

A fixed ${\tt J}$ is also interpretable, and applicable, if $x_i$ is discrete. Then each mass point is given its own bin and the results apply to the conditional distribution of $y_i$ given $x_i = x$. Again, our theoretical results apply directly to this case and one obtains simultaneous inference over the set of points. Cattaneo-Crump-Farrell-Feng_2024_AER give further discussion and examples.

As a practical compromise between the visual appeal and interpretation of a small, fixed ${\tt J}$ and the demand for consistent estimation, we propose a novel, albeit ad-hoc, adjustment to the estimator aimed at addressing the smoothing bias left by fixing ${\tt J}$ by adjusting the choice of polynomial order $p$. The standard approach fixes $p$ in advance and selects $J$ based on the data, but we can invert this and search for the value of $p$ for which the fixed ${\tt J}$ would be IMSE-optimal. That is, we solve for

where in principle the set $\mathcal{P}$ is all nonnegative integers, but in practice $\mathcal{P} = \{p_\mathtt{min},p_\mathtt{min}+1,\dots,p_\mathtt{max}-1,p_\mathtt{max}\}$, for some integers $0 \leq p_\mathtt{min} \leq p_\mathtt{max}$. The (in)flexibility of fixed $J = {\tt J}$ is offset by changing the polynomial approximation. This may have some practical appeal, but our theoretical results in the next section continue in the standard case of fixed $p$ and diverging $J$.

To implement the data-driven choice $p_{\mathtt{IMSE}}({\tt J},v)$, users needs to specify the desired (often small) number of bins $\mathtt{J}$, the derivative order $v$ of interest, and a (finite) set $\mathcal{P}$ of acceptable polynomial orders. The size of $\mathcal{P}$ is usually small since in practice $p=3$ or $4$ often suffices to yield a small IMSE-optimal number of bins. Then, for each value of $p$ in $\mathcal{P}$, we can implement the rule-of-thumb or direct plug-in procedure as described in Section (ref) to obtain $J_{\mathtt{IMSE}}(p,v)$. The “optimal” choice $p_{\mathtt{IMSE}}({\tt J},v)$ is the value of $p$ with the resulting $J_{\mathtt{IMSE}}(p,v)$ closest to $\mathtt{J}$.

Proofs

We begin with a subsection collecting some technical lemmas used in the proofs of our main results. We then collect all the proof of the results presented in this supplemental appendix, which are in several cases more general than those discussed in the main text. Some of our technical results may be of more broad independent interest in the nonlinear series estimation literature.

Technical Lemmas

We first give several simple facts about $\widehat\Delta$ in the following lemma, which are immediate from Assumption (ref)(ii).

The next lemma then verifies Assumption (ref)(ii) for the special case of quantile-spaced partitions. The proof is available in the supplemental appendix of Cattaneo-Crump-Farrell-Feng_2024_AER (see Section SA-3.1 therein) and thus omitted here.

The next three lemmas (ref)--(ref) concern the properties of binscatter basis functions. Their proofs are the same as those for quantile-based partitions that are available in the supplemental appendix of Cattaneo-Crump-Farrell-Feng_2024_AER (see Section SA-3.1 therein) and are omitted here to conserve space.

The following lemma provides a particular way to define $\boldsymbol{\beta}_0(\Delta)$ and $\widehat\boldsymbol{\beta}_0$ so that the required approximation rate is achieved. We define $$ \boldsymbol{\beta}_0^{\mathtt{LS}}(\Delta):=\operatorname*{arg\,min}_{\boldsymbol{\beta}\in\mathbb{R}^{K_{p,s}}}\; \mathbb{E}[(\mu_0(x_i)-\mathbf{b}_{p,s}(x_i;\Delta)'\boldsymbol{\beta})^2], \quad \widehat\boldsymbol{\beta}_0^{\mathtt{LS}}=\boldsymbol{\beta}_0^{\mathtt{LS}}(\widehat\Delta). $$

Next, the following maximal inequality is useful in our analysis. Its proof is available in Cattaneo-Feng-Underwood_2024_jasa and thus omitted here.

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

References

\begingroup \endgroup