2.1 Partially Linear Regression Model

The partially linear regression model is of the form

$$ \begin{align*} y_i =\mathbf{x}_i'\boldsymbol{\theta}_0 + \gamma_0(\mathbf{w}_i) + \varepsilon_i, \end{align*} $$

where $\boldsymbol {\theta }_0$ is the parameter of interest, $\gamma _0(\cdot )$ is an unknown function, and where $\mathbb {E}[\varepsilon _i|\mathbf {x}_i,\mathbf {w}_i]=0$ . Defining $\dot {y}_{i,j}=y_i-y_j$ and $\dot {\mathbf {x}}_{i,j}=\mathbf {x}_i-\mathbf {x}_j$ , a pairwise difference estimator of $\boldsymbol {\theta }_0$ can be based on

$$ \begin{align*} m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{1}{2}(\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} )^2. \end{align*} $$

Setting $\Theta =\mathbb {R}^k$ , the minimization problem in (1.1) admits a closed-form solution (provided that a nonnegative kernel function is used), namely,

$$ \begin{align*} \widehat{\boldsymbol{\theta}}_n = \left( \sum_{i<j} \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}'K_{h_n}(\mathbf{w}_i-\mathbf{w}_j) \right)^{-1} \sum_{i<j} \dot{\mathbf{x}}_{i,j}\dot{y}_{i,j}K_{h_n}(\mathbf{w}_i-\mathbf{w}_j). \end{align*} $$

2.2 Partially Linear Logit Model

The partially linear logit model studied here is of the form

$$ \begin{align*} y_i = \mathbb{1} \{\mathbf{x}_i'\boldsymbol{\theta}_0 + \gamma_0(\mathbf{w}_i) + \varepsilon_i \geq 0 \}, \end{align*} $$

where $\boldsymbol {\theta }_0$ is the parameter of interest, $\gamma _0(\cdot )$ is an unknown function, and

$$ \begin{align*} \mathbb{P}\left[\varepsilon_i \leq u | \mathbf{x}_i,\mathbf{w}_i \right] = \Lambda(u), \qquad \Lambda(u) = \frac{\exp(u)}{1+\exp(u)}. \end{align*} $$

The parameter $\boldsymbol {\theta }_0$ can be estimated using a pairwise difference estimator with $\Theta =\mathbb {R}^k$ and

$$ \begin{align*} m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\mathbb{1}\{\dot{y}_{i,j}\neq 0\} \left(y_i\ln\Lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}) + y_j\ln\Lambda(-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta})\right). \end{align*} $$

The minimization problem in (1.1) does not admit a closed-form solution, but it is convex (provided that a nonnegative kernel function is used) because $u\mapsto -\ln \Lambda (u)$ is.

2.3 Partially Linear Tobit Model

The partially linear censored regression model studied here is of the form

$$ \begin{align*} y_i =\max\{\mathbf{x}_i'\boldsymbol{\theta}_0+\gamma_0(\mathbf{w}_i)+\varepsilon_i,0\}, \end{align*} $$

where $\boldsymbol {\theta }_0$ is the parameter of interest, $\gamma _0(\cdot )$ is an unknown function, $x_i \perp\!\!\perp \varepsilon_i \mid w_i$ , and the conditional distribution of $\varepsilon _i$ given $\mathbf {w}_i$ admits a Lebesgue density. A pairwise difference estimator of $\boldsymbol {\theta }_0$ can be obtained by setting $\Theta =\mathbb {R}^k$ and employing

$$ \begin{align*} m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})-\tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol0), \end{align*} $$

where

$$ \begin{align*} \tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \begin{cases} |y_i| - \left( \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} + y_j \right) \operatorname*{\mathrm{sgn}}(y_i) & \text{if } \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} \leq -y_j \\ \left|\dot{y}_{i,j} - \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} \right| & \text{if } -y_j < \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} < y_i \\ |y_j| + \left(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}-y_i \right) \operatorname*{\mathrm{sgn}}(y_j) & \text{if } y_i \leq \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}. \end{cases} \end{align*} $$

Because $\tilde {m}_{\mathtt {PLT}}(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol 0)$ does not depend on $\boldsymbol {\theta }$ , its presence in $m_{\mathtt {PLT}}(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol {\theta })$ does not affect the minimization problem defining the estimator. Nevertheless, it is theoretically attractive to work with $m_{\mathtt {PLT}}$ rather than $\tilde {m}_{\mathtt {PLT}}$ , as doing so allows for weaker regularity conditions for the existence of the expectation of the objective function.

For future reference, we note that $m_{\mathtt {PLT}}$ admits the alternative representation

$$ \begin{align*} m_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \begin{cases} \left|\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}\right|-\left|\dot{y}_{i,j}\right| &\text{ if } y_i>0,y_j>0 \\ \max\{ y_i-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}, 0\}-y_i &\text{ if } y_i>0,y_j=0 \\ \max\{ y_j+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}, 0\}-y_j &\text{ if } y_i=0,y_j>0 \\ 0 &\text{ if } y_i=0,y_j=0. \end{cases} \end{align*} $$

The function $\boldsymbol {\theta }\mapsto m_{\mathtt {PLT}}(\mathbf {z}_1,\mathbf {z}_2;\boldsymbol {\theta })$ is convex and therefore so is the minimization problem in (1.1) (provided that a nonnegative kernel function is used).

3.1 Regularity Conditions

The following assumption guarantees, among other things, that $\boldsymbol {\theta }_n$ is well defined for large n and that the bias component $\boldsymbol {\theta }_n-\boldsymbol {\theta }_0$ vanishes asymptotically; for details, see Lemma 1 of Section 4.1.

The next assumption enables us to analyze the asymptotic properties of the noise component $\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_n$ . To accommodate examples (such as the partially linear Tobit model) where $\boldsymbol {\theta } \mapsto m(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol {\theta })$ is not fully differentiable, we assume the existence of derivative-like functions $\mathbf {s}(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol {\theta })\in \mathbb {R}^k$ and $\mathbf {H}(\mathbf {w}_i,\mathbf {w}_j;\boldsymbol {\theta },\mathbf {t})\in \mathbb {R}^{k\times k}$ such that, for any direction $\mathbf {t} \in \mathbb {R}^k$ , the (remainder) terms

$$ \begin{align*} r_{\mathbf{t}}(\boldsymbol{\theta},\tau)=\frac{m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}+\mathbf{t} \tau)-m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})}{\tau} - \mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})' \mathbf{t} \end{align*} $$

and

$$ \begin{align*} R_{\mathbf{t}}(\boldsymbol{\theta},\tau)=\frac{\mathbb{E}[r_{\mathbf{t}}(\boldsymbol{\theta},\tau)|\mathbf{w}_1,\mathbf{w}_2]}{\tau} - \frac{1}{2}\mathbf{t}'\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t})\mathbf{t} \end{align*} $$

are suitably small for $\boldsymbol {\theta }$ near $\boldsymbol {\theta }_0$ , $\tau>0$ near zero, and $\mathbf {w}_1 \approx \mathbf {w}_2$ . As further discussed below, functions $\mathbf {s}$ and $\mathbf {H}$ satisfying the following assumption exist (and are relatively easy to find) in each of our motivating examples.

3.2 Small-Bandwidth Asymptotics

Defining

$$ \begin{align*} \mathbf{V}_n = \mathbf{V}_n(h_n) &=\boldsymbol{\Gamma}_0^{-1}\left[ n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(K)\right]\boldsymbol{\Gamma}_0^{-1}, \qquad \\ \boldsymbol{\Delta}_0(K) &= \mathbb{E}[ \boldsymbol{\Xi}_0(\mathbf{w})] \int_{\mathbb{R}^d} K^2(\mathbf{u})d\mathbf{u}, \end{align*} $$

and letting $\Phi _k$ denote the distribution function of a k-dimensional standard Gaussian random vector, we have the following result.

Theorem 1. Suppose Assumptions 1 and 2 hold. If $n^2h_n^d\to \infty $ and if $h_n\to 0$ , then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \mathbf{V}_n^{-1/2}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0. \end{align*} $$

Under the assumptions of Theorem 1, the convergence rate of $\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_n$ is given by the magnitude of $\mathbf {V}_n^{-1/2}$ , namely,

$$ \begin{align*} \rho_n = \sqrt{\min\left(n,\binom{n}{2}h_n^d\right)}. \end{align*} $$

Provided that the bias is “small” in the sense that $\rho _n\|\boldsymbol {\theta }_n-\boldsymbol {\theta }_0\|\to 0$ , Theorem 1 therefore encompasses the following three distinct large-sample regimes:

• • Asymptotic Linearity: If $nh_n^d\to \infty $ , then $\sqrt {n}(\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_0)$ satisfies (1.3) with $\boldsymbol {\psi }_0 = \boldsymbol {\Gamma }_0^{-1}\boldsymbol {\xi }_0$ . In particular, $\sqrt {n}(\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_0)$ converges in law to a centered Gaussian distribution with variance

  $$ \begin{align*} \lim_{n\to\infty} n \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1} \boldsymbol{\Sigma}_0 \boldsymbol{\Gamma}_0^{-1}. \end{align*} $$

• • Root-n Consistency without Asymptotic Linearity: If $nh_n^d\to 2c \in (0,\infty )$ , then $\widehat {\boldsymbol {\theta }}_n$ is not asymptotically linear, but it is $\sqrt {n}$ -consistent, $\sqrt {n}(\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_0)$ converging in law to a centered Gaussian distribution with variance

  $$ \begin{align*} \lim_{n\to\infty} n \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1}\left[ \boldsymbol{\Sigma}_0 + \frac{1}{c} \boldsymbol{\Delta}_0(K)\right]\boldsymbol{\Gamma}_0^{-1}. \end{align*} $$

• • Slower than Root-n Consistency: If $nh_n^d\to 0$ (but $n^2h_n^d\to \infty $ ), then $\widehat {\boldsymbol {\theta }}_n$ is neither asymptotically linear nor $\sqrt {n}$ -consistent, but $\sqrt {n^2h_n^d/2}(\widehat {\boldsymbol {\theta }}_n-\boldsymbol {\theta }_0)$ converges in law to a centered Gaussian distribution with variance

  $$ \begin{align*} \lim_{n\to\infty} \binom{n}{2}h_n^d \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1}\boldsymbol{\Delta}_0(K)\boldsymbol{\Gamma}_0^{-1}. \end{align*} $$

The small-bandwidth component (i.e., the term involving $\boldsymbol {\Delta }_0(K)$ ) in $\mathbf {V}_n$ captures the additional estimation uncertainty generated from increasing the localization of the observation pairs. Incorporating this component in the approximate variance is key to enabling us to replace the condition $nh_n^d\to \infty $ by the weaker condition $n^2h_n^d\to \infty $ when obtaining a Gaussian approximation. As demonstrated by Cattaneo et al. (Reference Cattaneo, Farrell, Jansson and Masini2025) in a related context, incorporating the small-bandwidth component can furthermore lead to a more accurate distributional approximation (in a higher-order asymptotic sense) even under asymptotic linearity.

3.3 Debiasing

In Theorem 1, we centered the estimator $\widehat {\boldsymbol {\theta }}_n=\widehat {\boldsymbol {\theta }}_n(h_n)$ at $\boldsymbol {\theta }_n=\boldsymbol {\theta }(h_n)$ to circumvent bias issues. This section focuses on the bias term $\boldsymbol {\theta }_n-\boldsymbol {\theta }_0$ and introduces an automatic debiasing approach under the assumption that $\boldsymbol {\theta }_n-\boldsymbol {\theta }_0$ can be expanded in even powers of $h_n$ . To be specific, we follow Honoré and Powell (Reference Honoré, Powell, Andrews and Stock2005, Section 3.3) and discuss debiasing under the following high-level condition.

Assumption 3. For some even $\mathtt {S} \geq 0$ , $\boldsymbol {\theta }(\cdot )$ admits $\mathbf {b}_{2l}\in \mathbb {R}^k$ (for $l=1,\dots , \mathtt {S}/2$ ) such that

$$ \begin{align*} \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{\mathtt{S}/2} \mathbf{b}_{2l} h^{2l} + o(h^{\mathtt{S}})\qquad \text{as } h\downarrow 0. \end{align*} $$

The ease with which Assumption 3 can be verified depends on the magnitude of $\mathtt {S}$ . For instance, Assumption 1 implies that Assumption 3 holds with $\mathtt {S}=0$ . Under additional smoothness conditions and using symmetry of K, the following result gives conditions under which Assumption 3 holds with $\mathtt {S}=2$ . When stating the result, we employ the following standard multi-index notation: for $\boldsymbol {\nu } = (\nu _1,\ldots ,\nu _d)'\in \mathbb {Z}_+^d$ , $\bar {\mathbf {w}} = (\bar {w}_1,\ldots ,\bar {w}_d)'\in \mathbb {R}^d$ , and a sufficiently smooth-in- $\bar {\mathbf {w}}$ function $f(\mathbf {w},\bar {\mathbf {w}}),$

$$ \begin{align*} \partial_{\bar{\mathbf{w}}}^{\boldsymbol{\nu}} f(\mathbf{w},\bar{\mathbf{w}}) = \frac{\partial^{|\boldsymbol{\nu}|}}{\partial \bar w_1^{\nu_1} \dots \partial \bar w_d^{\nu_d}} f(\mathbf{w},\bar{\mathbf{w}}), \qquad |\boldsymbol{\nu}| = \sum_{j=1}^d \nu_j. \end{align*} $$

Proposition 1. Suppose Assumptions 1 and 2 hold and that

1. (i) $\int _{\mathbb {R}^d} \|\mathbf {u}\|^2 K(\mathbf {u})d\mathbf {u}<\infty $ , and

2. (ii) with probability one, $\bar {\mathbf {w}} \mapsto \boldsymbol {\varphi }(\mathbf {w},\bar {\mathbf {w}}) = \mathbb {E}[\mathbf {s}(\mathbf {z}_1,\mathbf {z}_2;\boldsymbol {\theta }_0)|\mathbf {w}_1=\mathbf {w},\mathbf {w}_2=\bar {\mathbf {w}}]f_{\mathbf {w}}(\bar {\mathbf {w}})$ is twice continuously differentiable with $\mathbb {E}[\sup _{\bar {\mathbf {w}}\in \mathcal {W}}\|\partial _{\bar {\mathbf {w}}}^{\boldsymbol {\nu }} \boldsymbol {\varphi }(\mathbf {w},\bar {\mathbf {w}})\|]<\infty $ for all $\boldsymbol {\nu }\in \mathbb {Z}_+^d$ with $|\boldsymbol {\nu }|\leq 2$ .

Then $\boldsymbol {\theta }(\cdot )$ admits a $\mathbf {b}_2\in \mathbb {R}^k$ such that

$$ \begin{align*} \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 = \mathbf{b}_2 h^2 + o(h^2) \qquad \text{ as } h\downarrow 0. \end{align*} $$

The proof of Proposition 1 leverages convexity and may therefore be of independent interest. The convexity argument in question can furthermore be adapted to form the basis of a verification by induction of Assumption 3 with $\mathtt {S}>2$ . Details are provided in Section 4.4, which describes the induction step for general $\mathtt {S}$ and states explicit (smoothness) conditions under which Assumption 3 holds with $\mathtt {S} =4$ .

To describe the debiasing procedure based on generalized jackknifing, we maintain Assumption 3, define $c_0=1$ , choose a nonnegative integer L, and let $\mathbf {c}=(c_0,\dots , c_{L})'$ be a vector of (distinct) positive constants such that the following vector is well defined:

$$ \begin{align*} \begin{pmatrix} \lambda_0(\mathbf{c}) \\ \lambda_1(\mathbf{c}) \\ \vdots \\ \lambda_{L}(\mathbf{c}) \end{pmatrix} = \begin{pmatrix} 1 & 1 &\dots &1 \\ 1 & c_1^2 &\dots &c_{L}^2 \\ \vdots & &\ddots & \\ 1& c_1^{2L} & \dots & c_{L}^{2L} \end{pmatrix}^{-1} \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}. \end{align*} $$

The debiased estimator is

$$ \begin{align*} \widetilde{\boldsymbol{\theta}}_n = \widetilde{\boldsymbol{\theta}}_n(\mathbf{c},h_n) = \sum_{l=0}^{L} \lambda_{l}(\mathbf{c}) \widehat{\boldsymbol{\theta}}_n(h_{n,l}), \qquad h_{n,l}=c_lh_n, \end{align*} $$

the construction of which involves solving $L+1$ convex optimization problems. As defined, the debiased estimator is a generalization of the original pairwise difference estimator because if $L=0$ , then $\mathbf {c}=1=\lambda _0$ and therefore $\widetilde {\boldsymbol {\theta }}_n = \widehat {\boldsymbol {\theta }}_n$ .

The next theorem generalizes Theorem 1 by establishing the small-bandwidth Gaussian approximation for $\widetilde {\boldsymbol {\theta }}_n$ . To state the theorem, let

$$ \begin{align*} \bar{\boldsymbol{\theta}}_n = \bar{\boldsymbol{\theta}}_n(\mathbf{c},h_n) = \sum_{l=0}^{L} \lambda_l(\mathbf{c}) \boldsymbol{\theta}(h_{n,l}) \end{align*} $$

and

$$ \begin{align*} \bar{\mathbf{V}}_n = \bar{\mathbf{V}}_n(\mathbf{c},h_n) &= \boldsymbol{\Gamma}_0^{-1}\left[n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1} h_n^{-d}\boldsymbol{\Delta}_0(\bar{K}) \right]\boldsymbol{\Gamma}_0^{-1}, \qquad \\ \bar{K}(\mathbf{u}) &= \bar{K}(\mathbf{u};\mathbf{c}) = \sum_{l=0}^{L}\lambda_{l}(\mathbf{c})K_{c_l}(\mathbf{u}). \end{align*} $$

As they should, the expressions have the feature that if $L=0$ , then $\bar {\boldsymbol {\theta }}_n = \boldsymbol {\theta }_n$ and $\bar {\mathbf {V}}_n=\mathbf {V}_n$ . Another noteworthy feature of the expressions is that debiasing via generalized jackknifing affects the variance $\bar {\mathbf {V}}_n$ only through the kernel shape entering its small-bandwidth component.

Theorem 2. Suppose Assumptions 1 and 2 hold. If $n^2h_n^d\to \infty $ and if $h_n\to 0$ , then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \bar{\mathbf{V}}_n^{-1/2}(\widetilde{\boldsymbol{\theta}}_n-\bar{\boldsymbol{\theta}}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0. \end{align*} $$

As a consequence, if Assumption 3 also holds with $\mathtt {S} \geq 2L+2$ , and if $nh_n^{4(L+1)}\to 0$ , then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \bar{\mathbf{V}}_n^{-1/2}(\widetilde{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0 ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0. \end{align*} $$

The magnitude of $\bar {\mathbf {V}}_n^{-1/2}$ is the same as that of $\mathbf {V}_n^{-1/2}$ . With obvious modifications, the discussion of $\widehat {\boldsymbol {\theta }}_n$ following Theorem 1 therefore applies to $\widetilde {\boldsymbol {\theta }}_n$ , the only noteworthy difference being that (by design) the relevant “small bias” condition is different (and typically milder) in the case of $\widetilde {\boldsymbol {\theta }}_n$ .

When $L=0$ (i.e., one uses the original pairwise difference estimator), the second part of Theorem 2 imposes $nh_n^4\to 0$ , which coincides with the standard small bias condition in (1.2), while the first part of the theorem still accommodates the small-bandwidth distributional approximation.

It is worth noting that for $L\geq 1$ , the equivalent kernel $\bar {K}$ is of higher order, even though the debiased estimator $\widetilde {\boldsymbol {\theta }}_n$ only employs estimators constructed using second-order kernels, thereby retaining the desired convexity for implementation. To be specific, if $\int _{\mathbb {R}^d} \|\mathbf {u}\|^{2L+2} K(\mathbf {u})d\mathbf {u} < \infty $ , then $\bar {K}$ is of order $2L+2$ because

$$ \begin{align*} \int_{\mathbb{R}^d} \bar{K}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L}\lambda_l(\mathbf{c}) \int_{\mathbb{R}^d} K_{c_{l}}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L}\lambda_l(\mathbf{c}) = 1 \end{align*} $$

and, for $\boldsymbol {\nu }=(\nu _1,\dots ,\nu _d)'\in \mathbb {Z}^d_+$ with $0<|\boldsymbol {\nu }|\leq 2L+1$ ,

$$ \begin{align*} \int_{\mathbb{R}^d} \mathbf{u}^{\boldsymbol{\nu}} \bar{K}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L} \lambda_l(\mathbf{c}) \int_{\mathbb{R}^d} \mathbf{u}^{\boldsymbol{\nu}} K_{c_l}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L} \lambda_l(\mathbf{c}) c_l^{|\boldsymbol{\nu}|} \int_{\mathbb{R}^d} \mathbf{u}^{\boldsymbol{\nu}} K(\mathbf{u})d\mathbf{u} =0, \end{align*} $$

where the last equality uses the defining property of $\{\lambda _l(\mathbf {c})\}$ and symmetry of K, and $\mathbf {u}^{\boldsymbol {\nu }}$ denotes $\prod _{j=1}^d u_j^{\nu _j}$ for $\mathbf {u}=(u_1,\dots ,u_d)'\in \mathbb {R}^d$ .

3.4 Bootstrapping

To develop feasible inference procedures that do not require (explicit) estimation of $\bar {\mathbf {V}}_n$ , we consider nonparametric bootstrap-based approximations to the distribution of $\widetilde {\boldsymbol {\theta }}_n$ . (Since $\widetilde {\boldsymbol {\theta }}_n=\widehat {\boldsymbol {\theta }}_n$ when $L=0$ , results for $\widehat {\boldsymbol {\theta }}_n$ can be extracted by setting $L=0$ in what follows.)

Letting $\mathbf {z}^{\hspace{0.5pt}*}_{1,n},\dots ,\mathbf {z}^{\hspace{0.5pt}*}_{n,n}$ denote a random sample from the empirical distribution of $\mathbf {z}_1,\dots ,\mathbf {z}_n$ , the defining property of $\widehat {\boldsymbol {\theta }}_n^{\hspace{0.5pt}*}(h)$ , the nonparametric bootstrap analog of $\widehat {\boldsymbol {\theta }}_n(h)$ , is the following:

$$ \begin{align*} \widehat{M}_n^{\hspace{0.5pt}*}(\widehat{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}(h);h) \leq \inf_{\boldsymbol{\theta}\in\Theta} \widehat{M}_n^{\hspace{0.5pt}*}(\boldsymbol{\theta};h) + o_{\mathbb{P}}(n^{-1}), \end{align*} $$

where

$$ \begin{align*} \widehat{M}_n^{\hspace{0.5pt}*}(\boldsymbol{\theta};h) = \binom{n}{2}^{-1} \sum_{i<j}m(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n};\boldsymbol{\theta}) K_h(\mathbf{w}^{\hspace{0.5pt}*}_{i,n}-\mathbf{w}^{\hspace{0.5pt}*}_{j,n}). \end{align*} $$

Similarly, the nonparametric bootstrap analog of $\widetilde {\boldsymbol {\theta }}_n$ is

$$ \begin{align*} \widetilde{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*} = \widetilde{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}(\mathbf{c})= \sum_{l=0}^{L} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}^{\hspace{0.5pt}*}(h_{n,l}). \end{align*} $$

The following theorem characterizes the large sample properties of $\widetilde {\boldsymbol {\theta }}_n^{\hspace{0.5pt}*}-\widetilde {\boldsymbol {\theta }}_n$ , the bootstrap counterpart of $\widetilde {\boldsymbol {\theta }}_n-\bar {\boldsymbol {\theta }}_n$ . In perfect analogy with the results in Cattaneo et al. (Reference Cattaneo, Crump and Jansson2014b), we find that the bootstrap distribution estimator is consistent only when $nh_n^d\to \infty $ , but otherwise exhibits a variance inflation making the distributional approximation inconsistent. To state the result, let $\mathbb {P}_n^{\hspace{0.5pt}*}[\cdot ]$ denote $\mathbb {P}[\cdot | \mathbf {z}_1,\ldots ,\mathbf {z}_n]$ , let $\to _{\mathbb {P}}$ denote convergence in probability, and define

$$ \begin{align*} \bar{\mathbf{V}}_n^{\hspace{0.5pt}*} = \bar{\mathbf{V}}_n^{\hspace{0.5pt}*}(\mathbf{c},h_n) = \boldsymbol{\Gamma}_0^{-1}\left[n^{-1}\boldsymbol{\Sigma}_0 + 3\binom{n}{2}^{-1} h_n^{-d}\boldsymbol{\Delta}_0(\bar{K}) \right]\boldsymbol{\Gamma}_0^{-1}. \end{align*} $$

Theorem 3. Suppose Assumptions 1 and 2 hold and that, for $\boldsymbol {\theta }$ near $\boldsymbol {\theta }_0$ (and with probability one), $m(\mathbf {z},\mathbf {z};\boldsymbol {\theta })=0$ and $\mathbf {s}(\mathbf {z},\mathbf {z};\boldsymbol {\theta })=\mathbf {0}$ . If $n^2h_n^d\to \infty $ and if $h_n\to 0$ , then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^{\hspace{0.5pt}*}\left[ \bar{\mathbf{V}}_n^{*-1/2}(\widetilde{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}-\widetilde{\boldsymbol{\theta}}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to_{\mathbb{P}} 0. \end{align*} $$

Because $\bar {\mathbf {V}}_n^{-1}\bar {\mathbf {V}}_n^{\hspace{0.5pt}*}\to \mathbf {I}_k$ if and only if $nh_n^d\to \infty $ (where $\mathbf {I}_k$ denotes the k-dimensional identity matrix), under the assumptions of Theorem 3,

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^{\hspace{0.5pt}*}\left[ \widetilde{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}-\widetilde{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \bar{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] \right| \to_{\mathbb{P}} 0 \end{align*} $$

if and only if $nh_n^d\to \infty $ . In particular, if $\liminf _{n\to \infty } nh_n^d<\infty $ , then the nonparametric bootstrap is inconsistent, albeit conservative in the sense that the (approximate) variance under the bootstrap distribution is larger than the (approximate) variance of the asymptotic distribution: $\bar {\mathbf {V}}_n^{\hspace{0.5pt}*}> \bar {\mathbf {V}}_n$ in a positive-definite sense.

The variance inflation problem associated with the nonparametric bootstrap under the small-bandwidth regime can be easily fixed by appropriately rescaling the bandwidth used for the bootstrap implementation of the pairwise difference estimator. Employing

$$ \begin{align*} \breve{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*} = \breve{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}(\mathbf{c},h_n)= \sum_{l=0}^{L} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}^{\hspace{0.5pt}*}(3^{1/d}h_{n,l}) \end{align*} $$

and centering its distribution at

$$ \begin{align*} \breve{\boldsymbol{\theta}}_n = \breve{\boldsymbol{\theta}}_n(\mathbf{c},h_n)= \sum_{l=0}^{L} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}(3^{1/d}h_{n,l}) \end{align*} $$

automatically adjusts the bootstrap variance, leading to a consistent distributional approximation. Indeed, the following result is an immediate consequence of Theorems 2 and 3.

Corollary 1. If the assumptions of Theorem 3 hold, then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^{\hspace{0.5pt}*}\left[ \breve{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*} - \breve{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \bar{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] \right| \to_{\mathbb{P}} 0. \end{align*} $$

As a consequence, if Assumption 3 also holds with $\mathtt {S} \geq 2L+2$ , and if $nh_n^{4(L+1)}\to 0$ , then

$$ \begin{align*} \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^{\hspace{0.5pt}*}\left[ \breve{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*} - \breve{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \boldsymbol{\theta}_{0} \leq \mathbf{t} \right] \right| \to_{\mathbb{P}} 0. \end{align*} $$

The statement of Corollary 1 emphasizes the rate-adaptive nature of the consistency property enjoyed by the bootstrap distributional approximation.

3.5 Discussion

Our results provide a simple practical message: bootstrap-based inference for convex pairwise difference estimators is highly sensitive to bandwidth choice under conventional implementations, but this sensitivity can be substantially reduced by combining debiasing with a bandwidth rescaling. In particular, the classical bootstrap is valid only over a relatively narrow range of bandwidth sequences, whereas generalized jackknifing enlarges the set of admissible large bandwidths, small-bandwidth asymptotics justifies an implementation enlarging the set of admissible small bandwidths, and the combination of the two yields the broadest range of valid inference.

To summarize these conclusions, suppose $\mathtt {S}\geq 2L+2$ , let $\alpha \in (0,1)$ and $\mathbf {a}\in \mathbb {R}^k$ be fixed, and consider the following family of percentile bootstrap confidence intervals (in the terminology of van der Vaart, Reference van der Vaart1998) for $\mathbf {a}'\boldsymbol {\theta }_0$ :

$$ \begin{align*} \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(B,L) = \left[\mathbf{a}'\widetilde{\boldsymbol{\theta}}_n(\mathbf{c},h_n) - \widetilde{q}^{\hspace{0.5pt}*}_{1-\alpha/2,n}(Bh_n,L) ~ , ~ \mathbf{a}'\widetilde{\boldsymbol{\theta}}_n(\mathbf{c},h_n) - \widetilde{q}^{\hspace{0.5pt}*}_{\alpha/2,n}(Bh_n,L) \right], \end{align*} $$

where

$$ \begin{align*} \widetilde{q}^{\hspace{0.5pt}*}_{t,n}(h,L) = \inf\left\{ q \in \mathbb{R} : \mathbb{P}_n^{\hspace{0.5pt}*} [\mathbf{a}'\widetilde{\boldsymbol{\theta}}_n^{\hspace{0.5pt}*}(\mathbf{c},h) - \mathbf{a}'\widetilde{\boldsymbol{\theta}}_n(\mathbf{c},h)\leq q] \geq t \right\}. \end{align*} $$

The notation highlights the two tuning parameters that vary across the different inference procedures: given a choice of bandwidth h used in the point estimator, the scaling constant B used in the bootstrap approximation, and the debiasing order $L=\dim (\mathbf {c})-1$ , which together with $\mathbf {c}$ determines the generalized jackknife weights.

To sharpen the practical interpretation of these conclusions, it is useful to view each procedure through its coverage function

$$ \begin{align*} \mathsf{CF}_{n,1-\alpha}(B,L)=\mathbb{P}\!\left[\mathbf{a}'\boldsymbol{\theta}_0\in \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(B,L)\right]. \end{align*} $$

Our theory implies that the main differences across procedures can be understood in terms of how $\mathsf {CF}_{n,1-\alpha }(B,L)$ behaves as the bandwidth sequence moves from the small-bandwidth region to the large-bandwidth region. In particular, the classical bootstrap tends to become conservative when the bandwidth is too small, because the bootstrap variance is no longer correctly matched under small-bandwidth asymptotics, and it tends to undercover when the bandwidth is too large, because smoothing bias becomes non-negligible. Debiasing mitigates the latter problem by shifting the onset of bias-driven undercoverage toward larger bandwidths, whereas bootstrap rescaling mitigates the former problem by restoring correct variance matching in the small-bandwidth region. Combining the two therefore produces the flattest coverage profile, in the sense that the coverage function remains close to its nominal level over the widest range of bandwidth sequences.

We formally compare the four implementations, ordered from the conventional procedure to our recommended method.

• • Classical Method: This corresponds to $(B,L)=(1,0)$ . If Assumptions 1–3 hold and

  $$ \begin{align*} n h_n^d\to\infty \qquad\text{and}\qquad n h_n^{4} \to 0, \end{align*} $$
  then
  $$ \begin{align*} \lim_{n \to \infty} \mathbb{P}\left[\mathbf{a}'\boldsymbol{\theta}_0 \in \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(1,0) \right] = 1 - \alpha. \end{align*} $$

  The admissible range of bandwidth sequences is narrow and requires $d<4$ . In terms of the coverage function, the classical percentile bootstrap has two distinct failure modes. When the bandwidth is too small, so that $n h_n^d$ is no longer large while $n^2 h_n^d\to \infty $ still holds, the procedure becomes conservative: $\mathsf {CF}_{n,1-\alpha }(1,0)$ rises above its nominal level because the bootstrap variance is too large relative to the sampling variance. When the bandwidth is too large, so that the bias condition fails, $\mathsf {CF}_{n,1-\alpha }(1,0)$ falls below its nominal level because smoothing bias is no longer negligible. Thus, the classical method delivers accurate coverage only over a narrow range of bandwidth choices.

• • Classical Debiased Method: This corresponds to $(B,L)=(1,L)$ for some $L\geq 1$ . If Assumptions 1–3 hold and

  $$ \begin{align*} n h_n^d\to\infty \qquad\text{and}\qquad n h_n^{4(L+1)} \to 0, \end{align*} $$
  then
  $$ \begin{align*} \lim_{n \to \infty} \mathbb{P}\left[\mathbf{a}'\boldsymbol{\theta}_0 \in \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(1,L) \right] = 1 - \alpha. \end{align*} $$

  The range of allowable bandwidth sequences is wider on the large-bandwidth side than for the classical method and requires $d<4(L+1)$ . Relative to the classical method, debiasing leaves the small-bandwidth behavior essentially unchanged: if the bandwidth is too small, the procedure may still be conservative because the bootstrap variance mismatch remains. Its gain is instead on the large-bandwidth side, where the onset of bias-driven undercoverage is pushed outward from the $n h_n^4$ scale to the weaker $n h_n^{4(L+1)}$ scale. Consequently, $\mathsf {CF}_{n,1-\alpha }(1,L)$ stays close to nominal coverage over a wider range of large-bandwidth sequences than $\mathsf {CF}_{n,1-\alpha }(1,0)$ .

• • Small-Bandwidth Method: This corresponds to $(B,L)=(3^{1/d},0)$ . If Assumptions 1–3 hold and

  $$ \begin{align*} n^2 h_n^d\to\infty \qquad\text{and}\qquad n h_n^{4} \to0, \end{align*} $$
  then
  $$ \begin{align*} \lim_{n \to \infty} \mathbb{P}\left[\mathbf{a}'\boldsymbol{\theta}_0 \in \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(3^{1/d},0) \right] = 1 - \alpha. \end{align*} $$

  The range of allowable bandwidth sequences is wider on the small-bandwidth side than for the classical approach and requires $d<8$ . Relative to the classical method, bootstrap rescaling corrects the variance mismatch responsible for small-bandwidth conservativeness. As a result, $\mathsf {CF}_{n,1-\alpha }(3^{1/d},0)$ remains close to its nominal level throughout the admissible small-bandwidth region. Its limitation is on the large-bandwidth side: because the procedure is not debiased, coverage still deteriorates once the bias condition $n h_n^4\to 0$ fails.

• • Small-Bandwidth Debiased Method: This corresponds to $(B,L)=(3^{1/d},L)$ for some $L\geq 1$ . If Assumptions 1–3 hold and

  $$ \begin{align*} n^2 h_n^d\to\infty \qquad\text{and}\qquad n h_n^{4(L+1)} \to 0, \end{align*} $$
  then
  $$ \begin{align*} \lim_{n \to \infty} \mathbb{P}\left[\mathbf{a}'\boldsymbol{\theta}_0 \in \mathsf{CI}_{n,1-\alpha}^{\hspace{0.5pt}*}(3^{1/d},L) \right] = 1 - \alpha. \end{align*} $$

  The range of allowable bandwidth sequences is the widest among the four procedures and requires $d<8(L+1)$ . This method combines the two improvements described above. On the small-bandwidth side, bootstrap rescaling removes the variance mismatch that would otherwise induce conservativeness. On the large-bandwidth side, debiasing delays the onset of bias-driven undercoverage. Accordingly, $\mathsf {CF}_{n,1-\alpha }(3^{1/d},L)$ remains close to its nominal level over the broadest range of bandwidth sequences considered in this article.

Figure 1 presents a schematic plot summarizing the main theoretical conclusions. Overall, the comparison suggests a clear practical recommendation: bootstrap inference based on small bandwidth and debiasing corrections provides the most robust coverage accuracy as a function of bandwidth choice. Thus, for empirical settings where bandwidth choice is uncertain or data-driven, this method provides the most reliable default.

Figure 1

Schematic coverage functions for four bootstrap confidence interval procedures.

Figure 1 Long description

A line graph with a vertical y-axis labeled Coverage and a horizontal x-axis labeled h sub n. The y-axis has two marked points: 1 at the top and 1 minus alpha lower down. A horizontal dashed line extends from 1 minus alpha across the plot. Four vertical dashed lines divide the x-axis into regions, labeled from left to right as: n squared h sub n super d approaches infinity, n h sub n super d approaches infinity, n h sub n super 4 approaches 0, and n h sub n super 8 approaches 0.

Four data series are plotted:

* C F sub n, 1 minus alpha (1, 0): An orange dash-dot line that starts near 1 on the far left, curves down to meet the 1 minus alpha line at the second vertical threshold, and then follows it until dropping sharply after the third threshold.

* C F sub n, 1 minus alpha (1, 1): A teal dashed line that follows the same path as the orange line, starting near 1 and converging to 1 minus alpha.

* C F sub n, 1 minus alpha (3 super 1/d, 0): A red dotted line that stays consistently at the 1 minus alpha level from the start until it drops sharply after the third vertical threshold.

* C F sub n, 1 minus alpha (3 super 1/d, 1): A solid blue line that stays at the 1 minus alpha level until the fourth vertical threshold, where it begins a steep decline.

All four curves eventually trend downward toward zero as h sub n increases to the right.

4.1 A Useful Lemma

The following lemma is used in the proofs of Theorems 1–3.

Lemma 1. Suppose that Assumption 1 holds. Then $\operatorname *{\mbox {arg min}}_{\boldsymbol {\theta }\in \Theta }M(\boldsymbol {\theta };h)$ is non-empty for $h>0$ near zero and

$$ \begin{align*} \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 = o(1) \qquad \text{as } h \downarrow 0. \end{align*} $$

If Assumption 2 also holds, then $\mathbb {E}[\mathbf {s}(\mathbf {z}_1,\mathbf {z}_2;\boldsymbol {\theta }(h) ) K_h(\mathbf {w}_1-\mathbf {w}_2)]=\mathbf {0}$ for $h>0$ near zero.

Proof. For every $\boldsymbol {\theta } \in \Theta $ ,

$$ \begin{align*} & M(\boldsymbol{\theta};h)\\ & \quad = \int_{\mathcal{W}} \int_{\mathbb{R}^d} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}-\mathbf{u} h] f_{\mathbf{w}}(\mathbf{w}) f_{\mathbf{w}}(\mathbf{w} - \mathbf{u} h) K(\mathbf{u}) d\mathbf{u} d\mathbf{w} \\ &\quad \to \int_{\mathcal{W}} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}] f_{\mathbf{w}}(\mathbf{w})^2 d\mathbf{w} = M_0(\boldsymbol{\theta}) \qquad \text{as } h \downarrow 0, \end{align*} $$

the convergence being uniform on compact subsets of $\Theta $ because $\boldsymbol {\theta } \mapsto M(\boldsymbol {\theta };h)$ is convex (e.g., Hjort and Pollard, Reference Hjort and Pollard1993, Lemma 1).

Take any $\epsilon> 0$ with $\Theta ^\epsilon _0 = \{\boldsymbol {\theta } \in \mathbb {R}^k: \|\boldsymbol {\theta }-\boldsymbol {\theta }_0\| \leq \epsilon \} \subseteq \Theta $ . By the preceding paragraph,

$$ \begin{align*} \sup_{\boldsymbol{\theta} \in \Theta^\epsilon_0} |M(\boldsymbol{\theta};h) - M_0(\boldsymbol{\theta})| \leq \frac{1}{2} \left(\inf_{\boldsymbol{\theta}:\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|=\epsilon} M_0(\boldsymbol{\theta}) - M_0(\boldsymbol{\theta}_0)\right) \end{align*} $$

for $h>0$ near zero. For any such h and any $\boldsymbol {\theta } \in \Theta \setminus \Theta ^\epsilon _0$ , we have

$$ \begin{align*} \epsilon_{\boldsymbol{\theta}} = \frac{\epsilon}{\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|} \in (0,1), \end{align*} $$

and therefore, by convexity of $\boldsymbol {\theta } \mapsto M(\boldsymbol {\theta };h)$ ,

$$ \begin{align*} M(\epsilon_{\boldsymbol{\theta}} \boldsymbol{\theta} + (1-\epsilon_{\boldsymbol{\theta}}) \boldsymbol{\theta}_0;h) \leq \epsilon_{\boldsymbol{\theta}} M(\boldsymbol{\theta};h) + (1-\epsilon_{\boldsymbol{\theta}}) M(\boldsymbol{\theta}_0;h), \end{align*} $$

which rearranges as

$$ \begin{align*} M(\boldsymbol{\theta};h) - M(\boldsymbol{\theta}_0;h) \geq \frac{1}{\epsilon_{\boldsymbol{\theta}}} [M(\epsilon_{\boldsymbol{\theta}} \boldsymbol{\theta} + (1-\epsilon_{\boldsymbol{\theta}}) \boldsymbol{\theta}_0;h) - M(\boldsymbol{\theta}_0;h)] \geq 0. \end{align*} $$

As a consequence,

$$ \begin{align*} \inf_{\boldsymbol{\theta} \in \Theta}M(\boldsymbol{\theta};h) = \inf_{\boldsymbol{\theta} \in \Theta^\epsilon_0}M(\boldsymbol{\theta};h) = \min_{\boldsymbol{\theta} \in \Theta^\epsilon_0}M(\boldsymbol{\theta};h), \end{align*} $$

where the last equality uses continuity of $\boldsymbol {\theta } \mapsto M(\boldsymbol {\theta };h)$ and compactness of $\Theta ^\epsilon _0$ .

The above argument shows in particular that $\boldsymbol {\theta }(h) \in \Theta ^\epsilon _0$ for $h>0$ near zero.

If Assumption 2 also holds, then, for $\boldsymbol {\theta }$ near $\boldsymbol {\theta }_0$ , $h>0$ near zero, and for any $\mathbf {t} \in \mathbb {R}^k$ ,

$$ \begin{align*} \!\left| \!\frac{M(\boldsymbol{\theta}+\mathbf{t}\tau;h)-M(\boldsymbol{\theta};h)}{\tau} - \mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}) K_h(\mathbf{w}_1-\mathbf{w}_2)]' \mathbf{t} \right|\!\! &\leq\! \left| \mathbb{E}[r_{\mathbf{t}}(\boldsymbol{\theta},\tau) K_h(\mathbf{w}_1\! -\mathbf{w}_2)] \right| \\ &\to 0 \qquad \text{as } \tau \downarrow 0, \end{align*} $$

implying that for $h>0$ near zero, $\boldsymbol {\theta } \mapsto M(\boldsymbol {\theta };h)$ is (directionally) differentiable near $\boldsymbol {\theta }_0$ , the directional derivative $\mathbb {E}[\mathbf {s}(\mathbf {z}_1,\mathbf {z}_2;\boldsymbol {\theta }) K_h(\mathbf {w}_1-\mathbf {w}_2)]' \mathbf {t}$ being zero when $\boldsymbol {\theta } = \boldsymbol {\theta }(h)$ because $\boldsymbol {\theta }(h)$ minimizes $M(\boldsymbol {\theta };h).$

4.2 Proof of Theorems 1 and 2

Theorem 1 can be obtained from Theorem 2 by setting $L=0$ in the latter, so it suffices to prove Theorem 2. To do so, for $l \in \{0,\dots ,L\}$ , let $\widehat {\boldsymbol {\theta }}_{n,l}=\widehat {\boldsymbol {\theta }}_n(h_{n,l}),\boldsymbol {\theta }_{n,l}=\boldsymbol {\theta }(h_{n,l}),$ and

$$ \begin{align*} \widehat{\mathbf{U}}_{n,l} = \binom{n}{2}^{-1} \sum_{i<j} \mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j), \qquad \mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)= \mathbf{s}_{n,l}(\mathbf{z}_i,\mathbf{z}_j) - \mathbb{E}[\mathbf{s}_{n,l}(\mathbf{z}_1,\mathbf{z}_2)], \end{align*} $$

where

$$ \begin{align*} \mathbf{s}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)=\mathbf{s}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}_{n,l})K_{h_{n,l}}(\mathbf{w}_i-\mathbf{w}_j). \end{align*} $$

By Lemma 1, $\lim _{n \to \infty }\boldsymbol {\theta }_{n,l}=\boldsymbol {\theta }_0$ and $\mathbb {E}[\mathbf {s}_{n,l}(\mathbf {z}_i,\mathbf {z}_j)]=0$ for large n.

Suppose that

(4.1) $$ \begin{align} \widehat{\boldsymbol{\theta}}_{n,l}-\boldsymbol{\theta}_{n,l} = -\boldsymbol{\Gamma}_0^{-1}\widehat{\mathbf{U}}_{n,l} + o_{\mathbb{P}}(\rho_n^{-1}) \quad \text{for} \quad l \in\{0,\dots,L\}. \end{align} $$

Then

$$ \begin{align*} \widetilde{\boldsymbol{\theta}}_n-\bar{\boldsymbol{\theta}}_n = -\boldsymbol{\Gamma}_0^{-1}\widetilde{\mathbf{U}}_n + o_{\mathbb{P}}(\rho_n^{-1}), \end{align*} $$

where

$$ \begin{align*} \widetilde{\mathbf{U}}_n = \sum_{l=0}^{L} \lambda_l(\mathbf{c})\widehat{\mathbf{U}}_{n,l} \end{align*} $$

satisfies

(4.2) $$ \begin{align} \left[ n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(\bar{K})\right]^{-1/2} \widetilde{\mathbf{U}}_n \rightsquigarrow \mathsf{N} (\mathbf{0}_{k \times 1},\mathbf{I}_k) \end{align} $$

because, letting $\sum _i$ denote $\sum _{i=1}^n$ ,

$$ \begin{align*} \widetilde{\mathbf{U}}_n = \widetilde{\mathbf{L}}_n + \widetilde{\mathbf{W}}_n, \end{align*} $$

where

$$ \begin{align*} \widetilde{\mathbf{L}}_n = \frac{2}{n} \sum_i \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i), \quad \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i)=\sum_{l=0}^{L} \lambda_l(\mathbf{c}) \mathbb{E}[\mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)|\mathbf{z}_i] \quad (j \neq i), \end{align*} $$

and

$$ \begin{align*} \widetilde{\mathbf{W}}_n = \binom{n}{2}^{-1} \sum_{i<j} \widetilde{\boldsymbol{\omega}}_n(\mathbf{z}_i,\mathbf{z}_j), \quad \widetilde{\boldsymbol{\omega}}_n(\mathbf{z}_i,\mathbf{z}_j)=\sum_{l=0}^{L} \lambda_l(\mathbf{c}) \left[\mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j) - \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i) - \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_j)\right], \end{align*} $$

satisfy

$$ \begin{align*} \left(\begin{array}{c} \sqrt{n}\widetilde{\mathbf{L}}_n \\ \sqrt{\binom{n}{2}h_n^d} \widetilde{\mathbf{W}}_n \end{array}\right) \rightsquigarrow \mathsf{N}\left( \left[\begin{array}{c} \mathbf{0}_{k \times 1} \\ \mathbf{0}_{k \times 1} \end{array}\right], \left[\begin{array}{cc} \boldsymbol{\Sigma}_0 &\mathbf{0}_{k \times k} \\ \mathbf{0}_{k \times k} &\boldsymbol{\Delta}_0(\bar{K}) \end{array}\right] \right), \end{align*} $$

as can be shown by means of the Cramér–Wold device and the central limit theorem of Heyde and Brown (Reference Heyde and Brown1970), the latter being applicable because it follows from routine calculations that for every $\boldsymbol {\mu }_1,\boldsymbol {\mu }_2 \in \mathbb {R}^k$ , we have

$$ \begin{align*} \sigma_n^2 = \sum_i\mathbb{V}[g_{i,n}] = \boldsymbol{\mu}_1'\boldsymbol{\Sigma}_0\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2'\boldsymbol{\Delta}_0(\bar{K})\boldsymbol{\mu}_2 + o(1), \end{align*} $$

$$ \begin{align*} \sum_i \mathbb{E}[g_{i,n}^4] = o(1), \end{align*} $$

and

$$ \begin{align*} \mathbb{V}\left[\sum_i\sigma_{i,n}^2-\sigma_n^2\right] = o(1),\qquad \sigma_{i,n}^2=\mathbb{V}[g_{i,n}|\mathbf{z}_1,\dots,\mathbf{z}_{i-1}], \end{align*} $$

where

$$ \begin{align*} g_{i,n} = g_{i,n}(\boldsymbol{\mu}) = \frac{2}{\sqrt{n}} \boldsymbol{\mu}_1'\widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i) + \sqrt{\binom{n}{2}^{-1}h_n^{d}}\sum_{j=1}^{i-1}\boldsymbol{\mu}_2'\widetilde{\boldsymbol{\omega}}_{n}(\mathbf{z}_i,\mathbf{z}_j). \end{align*} $$

The proof of Theorem 2 can therefore be completed by verifying (4.1).

To do so, we leverage convexity. For any $l \in \{0,\dots ,L\}$ and any $\mathbf {t} \in \mathbb {R}^k$ , it can be shown that

$$ \begin{align*} \lim_{\tau \downarrow 0,h \downarrow 0,\boldsymbol{\theta} \to \boldsymbol{\theta}_0} \mathbb{E} \left[ \mathbb{E}[ r_{\mathbf{t}}(\boldsymbol{\theta},\tau) K_h(\mathbf{w}_1-\mathbf{w}_2) | \mathbf{z}_1 ] ^2 \right] = 0 \end{align*} $$

and

$$ \begin{align*} \lim_{\tau \downarrow 0,h \downarrow 0,\boldsymbol{\theta} \to \boldsymbol{\theta}_0} h^d \mathbb{E}[ r_{\mathbf{t}}(\boldsymbol{\theta},\tau)^2 K_h(\mathbf{w}_1-\mathbf{w}_2)^2 ] = 0, \end{align*} $$

and it therefore follows from a Hoeffding decomposition that

$$ \begin{align*} &\rho_n^2\left[\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \right] \\ &= \rho_n^2 [M(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - M(\boldsymbol{\theta}_{n,l};h_{n,l})] + \mathbf{t}'\rho_n\widehat{\mathbf{U}}_{n,l} + o_{\mathbb{P}}(1). \end{align*} $$

Moreover,

$$ \begin{align*} \rho_n^2 [M(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - M(\boldsymbol{\theta}_{n,l};h_{n,l}) ] \to \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t}, \end{align*} $$

and proceeding as in the proof of (4.2), it can be shown that $\rho _n\widehat {\mathbf {U}}_{n,l} = O_{\mathbb {P}}(1)$ . Because $\boldsymbol {\Gamma }_0$ is positive definite and because $\mathbf {t} \mapsto \widehat {M}_n(\boldsymbol {\theta }_{n,l}+\mathbf {t} \rho _n^{-1};h_{n,l})$ is convex (almost surely), the corollary following Hjort and Pollard (Reference Hjort and Pollard1993, Lemma 2) implies that (4.1) holds.

4.3 Proof of Theorem 3

The proof of Theorem 3 is a natural bootstrap analog of the proof of Theorem 2.

For $l \in \{0,\dots ,L\}$ , let $\widehat {\boldsymbol {\theta }}^{\hspace{0.5pt}*}_{n,l}=\widehat {\boldsymbol {\theta }}^{\hspace{0.5pt}*}_n(h_{n,l})$ and

$$ \begin{align*} \widehat{\mathbf{U}}^{\hspace{0.5pt}*}_{n,l} = \binom{n}{2}^{-1} \sum_{i<j} \mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}), \qquad \mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}) &= \mathbf{s}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}) \\ &\quad - \mathbb{E}^{\hspace{0.5pt}*}_n[\mathbf{s}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{1,n},\mathbf{z}^{\hspace{0.5pt}*}_{2,n})], \end{align*} $$

where $\mathbb {E}^{\hspace{0.5pt}*}_n[\cdot ]$ denotes $\mathbb {E}[\cdot |\mathbf {z}_1,\dots ,\mathbf {z}_n]$ .

It suffices to show that

(4.3) $$ \begin{align} \widehat{\boldsymbol{\theta}}^{\hspace{0.5pt}*}_{n,l}-\boldsymbol{\theta}_{n,l} = -\boldsymbol{\Gamma}_0^{-1}\left(\widehat{\mathbf{U}}^{\hspace{0.5pt}*}_{n,l} + \widehat{\mathbf{U}}_{n,l} \right) + o_{\mathbb{P}}(\rho_n^{-1}) \quad \text{for} \quad l \in\{0,\dots,L\} \end{align} $$

and that

(4.4) $$ \begin{align} \left[ n^{-1}\boldsymbol{\Sigma}_0 + 3 \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(\bar{K})\right]^{-1/2} \widetilde{\mathbf{U}}^{\hspace{0.5pt}*}_n \rightsquigarrow_{\mathbb{P}} \mathsf{N} (\mathbf{0}_{k \times 1},\mathbf{I}_k), \end{align} $$

where

$$ \begin{align*} \widetilde{\mathbf{U}}^{\hspace{0.5pt}*}_n = \sum_{l=0}^{L} \lambda_l(\mathbf{c})\widehat{\mathbf{U}}^{\hspace{0.5pt}*}_{n,l}. \end{align*} $$

For every $\mathbf {t} \in \mathbb {R}^k$ , using a Hoeffding decomposition and the fact that $m(\mathbf {z},\mathbf {z};\boldsymbol {\theta }_{n,l})=0$ and $\mathbf {s}(\mathbf {z},\mathbf {z};\boldsymbol {\theta }_{n,l})=\mathbf {0}$ for large n, we have

$$ \begin{align*} &\rho_n^2\left[\widehat{M}^{\hspace{0.5pt}*}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}^{\hspace{0.5pt}*}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \right] \\ &= \rho_n^2 (1+o(1)) [\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) ] + \mathbf{t}'\rho_n\widehat{\mathbf{U}}^{\hspace{0.5pt}*}_{n,l} + o_{\mathbb{P}}(1), \end{align*} $$

where it can be shown that $\rho _n\widehat {\mathbf {U}}^{\hspace{0.5pt}*}_{n,l}=O_{\mathbb {P}}(1)$ and where it follows from the proof of (4.1) that

$$ \begin{align*} \rho_n^2 \left[\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \right] = \frac{1}{2}\mathbf{t}'\boldsymbol{\Gamma}_0\mathbf{t} + \mathbf{t}'\rho_n\widehat{\mathbf{U}}_{n,l} + o_{\mathbb{P}}(1), \end{align*} $$

where $\rho _n\widehat {\mathbf {U}}_{n,l}=O_{\mathbb {P}}(1)$ . In other words,

$$ \begin{align*} &\rho_n^2\left[\widehat{M}^{\hspace{0.5pt}*}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}^{\hspace{0.5pt}*}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \right] \\ &= \frac{1}{2}\mathbf{t}'\boldsymbol{\Gamma}_0\mathbf{t} + \mathbf{t}'\rho_n\left(\widehat{\mathbf{U}}^{\hspace{0.5pt}*}_{n,l} + \widehat{\mathbf{U}}_{n,l}\right) + o_{\mathbb{P}}(1) \qquad \text{for every } \mathbf{t} \in \mathbb{R}^k. \end{align*} $$

Because $\boldsymbol {\Gamma }_0$ is positive definite and because $\mathbf {t} \mapsto \widehat {M}^{\hspace{0.5pt}*}_n(\boldsymbol {\theta }_{n,l}+\mathbf {t} \rho _n^{-1};h_{n,l})$ is convex (almost surely), the corollary following Hjort and Pollard (Reference Hjort and Pollard1993, Lemma 2) implies that (4.3) holds.

To prove (4.4), we begin by decomposing $\widetilde {\mathbf {U}}^{\hspace{0.5pt}*}_n$ as

$$ \begin{align*} \widetilde{\mathbf{U}}^{\hspace{0.5pt}*}_n = \widetilde{\mathbf{L}}^{\hspace{0.5pt}*}_n + \widetilde{\mathbf{W}}^{\hspace{0.5pt}*}_n, \end{align*} $$

where

$$ \begin{align*} \widetilde{\mathbf{L}}^{\hspace{0.5pt}*}_n = \frac{2}{n} \sum_i \widetilde{\boldsymbol{\ell}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n}), \qquad \widetilde{\boldsymbol{\ell}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n})=\sum_{l=0}^{L} \lambda_l(\mathbf{c}) \mathbb{E}^{\hspace{0.5pt}*}_n[\mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n})|\mathbf{z}^{\hspace{0.5pt}*}_{i,n}] \qquad (j \neq i), \end{align*} $$

and

$$ \begin{align*} &\widetilde{\mathbf{W}}^{\hspace{0.5pt}*}_n = \binom{n}{2}^{-1} \sum_{i<j} \widetilde{\boldsymbol{\omega}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}), \qquad \\ & \widetilde{\boldsymbol{\omega}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n})=\sum_{l=0}^{L} \lambda_l(\mathbf{c}) \left[\mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}) - \widetilde{\boldsymbol{\ell}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n}) - \widetilde{\boldsymbol{\ell}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{j,n})\right]. \end{align*} $$

Defining

$$ \begin{align*} \pi_n = \frac{\sqrt{nh_n^d}}{1+\sqrt{nh_n^d}}, \end{align*} $$

routine calculations can be used to show that for every $\boldsymbol {\mu }_1,\boldsymbol {\mu }_2 \in \mathbb {R}^k$ , we have

$$ \begin{align*} \widehat{\sigma}_n^2 = \sum_i\mathbb{V}^{\hspace{0.5pt}*}_n[g^{\hspace{0.5pt}*}_{i,n}] = \boldsymbol{\mu}_1'\left[\pi_n^2\boldsymbol{\Sigma}_0 + 4(1-\pi_n)^2\boldsymbol{\Delta}_0(\bar{K})\right]\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2'\boldsymbol{\Delta}_0(\bar{K})\boldsymbol{\mu}_2 + o_{\mathbb{P}}(1), \end{align*} $$

$$ \begin{align*} \sum_i \mathbb{E}^{\hspace{0.5pt}*}_n[g^{*4}_{i,n}] = o_{\mathbb{P}}(1), \end{align*} $$

and

$$ \begin{align*} \mathbb{V}^{\hspace{0.5pt}*}_n\left[\sum_i\widehat{\sigma}_{i,n}^2-\widehat{\sigma}_n^2\right] = o_{\mathbb{P}}(1), \qquad \widehat{\sigma}_{i,n}^2=\mathbb{V}^{\hspace{0.5pt}*}_n[g^{\hspace{0.5pt}*}_{i,n}|\mathbf{z}^{\hspace{0.5pt}*}_{1,n},\dots,\mathbf{z}^{\hspace{0.5pt}*}_{i-1,n}], \end{align*} $$

where $\mathbb {V}^{\hspace{0.5pt}*}_n[\cdot ]$ denotes $\mathbb {V}[\cdot |\mathbf {z}_1,\dots ,\mathbf {z}_n]$ and

$$ \begin{align*} g^{\hspace{0.5pt}*}_{i,n} = g^{\hspace{0.5pt}*}_{i,n}(\boldsymbol{\mu}) = \frac{\pi_n}{\sqrt{n}} 2\boldsymbol{\mu}_1'\widetilde{\boldsymbol{\ell}}^{\hspace{0.5pt}*}_n(\mathbf{z}^{\hspace{0.5pt}*}_{i,n}) + \sqrt{\binom{n}{2}^{-1}h_n^{d}}\sum_{j=1}^{i-1}\boldsymbol{\mu}_2'\widetilde{\boldsymbol{\omega}}^{\hspace{0.5pt}*}_{n}(\mathbf{z}^{\hspace{0.5pt}*}_{i,n},\mathbf{z}^{\hspace{0.5pt}*}_{j,n}). \end{align*} $$

The Cramér–Wold device and the central limit theorem of Heyde and Brown (Reference Heyde and Brown1970) therefore imply that if $\pi _n \to \pi _0 \in [0,1]$ , then

$$ \begin{align*} \left(\begin{array}{c} \sqrt{n}\pi_n\widetilde{\mathbf{L}}^{\hspace{0.5pt}*}_n \\ \sqrt{\binom{n}{2}h_n^d} \widetilde{\mathbf{W}}^{\hspace{0.5pt}*}_n \end{array}\right) \rightsquigarrow_{\mathbb{P}} \mathsf{N}\left( \left[\begin{array}{c} \mathbf{0}_{k \times 1} \\ \mathbf{0}_{k \times 1} \end{array}\right], \left[\begin{array}{cc} \pi_0^2\boldsymbol{\Sigma}_0 + 4(1-\pi_0)^2\boldsymbol{\Delta}_0(\bar{K}) &\mathbf{0}_{k \times k} \\ \mathbf{0}_{k \times k} &\boldsymbol{\Delta}_0(\bar{K}) \end{array}\right] \right). \end{align*} $$

Whether or not $\pi _n$ is convergent, the result (4.4) can be obtained from the preceding display by arguing along subsequences (if necessary).

4.4 Verifying Assumption 3

It follows from Lemma 1 that if Assumption 1 holds, then so does Assumption 3 with $\mathtt {S}=0$ . This observation provides the base case for an induction argument. To describe the induction step, suppose that for some even $S \geq 0$ , we have

$$ \begin{align*} \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{S/2} \mathbf{b}_{2l} h^{2l} + o(h^{S}) \qquad \text{ as } h\downarrow 0. \end{align*} $$

Suppose also that, for every $\mathbf {t} \in \mathbb {R}^k$ and some $\boldsymbol {\beta }_{S+2} \in \mathbb {R}^k$ , we have

$$ \begin{align*} &h^{-2(S+2)} \left[ M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{S/2} \mathbf{b}_{2l} h^{2l} + \mathbf{t} h^{S+2};h \right) - M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{S/2} \mathbf{b}_{2l} h^{2l}; h \right)\right] \\ &= \mathbf{t}' \boldsymbol{\beta}_{S+2} + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \qquad \text{as } h \downarrow 0. \end{align*} $$

Then, the corollary following Hjort and Pollard (Reference Hjort and Pollard1993, Lemma 2) implies that

$$ \begin{align*} h^{-(S+2)} \left( \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 - \sum_{l=1}^{S/2} \mathbf{b}_{2l} h^{2l} \right) &= \operatorname*{\mbox{arg min}}_{\mathbf{t} \in \mathbb{R}^k} M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{S/2} \mathbf{b}_{2l} h^{2l} + \mathbf{t} h^{S+2};h \right) \\ &= - \boldsymbol{\Gamma}_0^{-1} \boldsymbol{\beta}_{S+2} + o(1) \qquad \text{ as } h\downarrow 0; \end{align*} $$

that is, defining $\mathbf {b}_{S+2} = - \boldsymbol {\Gamma }_0^{-1} \boldsymbol {\beta }_{S+2}$ , we have

$$ \begin{align*} \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{(S+2)/2} \mathbf{b}_{2l} h^{2l} + o(h^{S+2}) \qquad \text{ as } h\downarrow 0. \end{align*} $$

To complete the proof of Proposition 1, it therefore suffices to note that (for every $\mathbf {t} \in \mathbb {R}^k$ and) under the assumptions of the proposition, we have

$$ \begin{align*} h^{-4} \left[ M (\boldsymbol{\theta}_0 + \mathbf{t} h^2;h ) - M (\boldsymbol{\theta}_0; h )\right] = \mathbf{t}' \boldsymbol{\beta}_2 + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \qquad \text{as } h \downarrow 0, \end{align*} $$

where

$$ \begin{align*} \boldsymbol{\beta}_2 = \frac{1}{2}\sum_{i=1}^d\sum_{j=1}^d \left(\left.\int_{\mathcal{W}} \frac{\partial^2 \boldsymbol{\varphi}(\mathbf{w},\bar{\mathbf{w}})}{\partial \bar w_i \partial \bar w_j}\right|_{\bar{\mathbf{w}}=\mathbf{w}} f_{\mathbf{w}}(\mathbf{w})d\mathbf{w}\right) \left(\int_{\mathbb{R}^d} u_i u_j K(\mathbf{u})d\mathbf{u}\right). \end{align*} $$

Similarly, if, in addition to the assumptions of Proposition 1, it is assumed that for every $\mathbf {t} \in \mathbb {R}^k$ and for some $\boldsymbol {\beta }_4 \in \mathbb {R}^k$ , we have

$$ \begin{align*} h^{-8} \left[ M (\boldsymbol{\theta}_0 + \mathbf{b}_2 h^2 + \mathbf{t} h^4;h ) - M (\boldsymbol{\theta}_0 + \mathbf{b}_2 h^2; h ) \!\right]\! = \mathbf{t}' \boldsymbol{\beta}_4 + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \quad \text{as } h \downarrow 0, \end{align*} $$

then Assumption 3 holds with $\mathtt {S}=4$ . One set of sufficient conditions for this to occur is that Assumptions 1 and 2 hold and that, for every $\mathbf {t} \in \mathbb {R}^k$ , the following are satisfied (with probability one):

1. (i) $\int _{\mathbb {R}^d} \|\mathbf {u}\|^4 K(\mathbf {u})d\mathbf {u}<\infty $ .

2. (ii) $\bar {\mathbf {w}} \mapsto \boldsymbol {\varphi }(\mathbf {w},\bar {\mathbf {w}}) = \mathbb {E}[\mathbf {s}(\mathbf {z}_1,\mathbf {z}_2;\boldsymbol {\theta }_0)|\mathbf {w}_1=\mathbf {w},\mathbf {w}_2=\bar {\mathbf {w}}]f_{\mathbf {w}}(\bar {\mathbf {w}})$ is four times continuously differentiable with $\mathbb {E}[\sup _{\bar {\mathbf {w}}\in \mathcal {W}}\|\partial _{\bar {\mathbf {w}}}^{\boldsymbol {\nu }} \boldsymbol {\varphi }(\mathbf {w},\bar {\mathbf {w}})\|]<\infty $ for all $\boldsymbol {\nu }\in \mathbb {Z}_+^d$ with $|\boldsymbol {\nu }|\leq 4$ .

3. (iii) $f_{\mathbf {w}}$ is twice continuously differentiable and $\bar {\mathbf {w}} \mapsto \mathbf {H}(\mathbf {w},\bar {\mathbf {w}};\boldsymbol {\theta }_0,\mathbf {t})$ is twice continuously differentiable with $\mathbb {E}[\sup _{\bar {\mathbf {w}}\in \mathcal {W}}\|\partial _{\bar {\mathbf {w}}}^{\boldsymbol {\nu }} \mathbf {H}(\mathbf {w},\bar {\mathbf {w}};\boldsymbol {\theta }_0,\mathbf {t})f_{\mathbf {w}}(\bar {\mathbf {w}})\|]<\infty $ for all $\boldsymbol {\nu }\in \mathbb {Z}_+^d$ with $|\boldsymbol {\nu }|\leq 2$ .

4. (iv) For some function $\dot {\mathbf {H}}(\mathbf {w},\bar {\mathbf {w}};\boldsymbol {\theta },\mathbf {t})\in \mathbb {R}^{k\times k}$ , $\bar {\mathbf {w}} \mapsto \dot {\mathbf {H}}(\mathbf {w},\bar {\mathbf {w}};\boldsymbol {\theta }_0,\mathbf {t})$ is continuous,

   $$ \begin{align*} \mathbb{E}\left[\sup_{\bar{\mathbf{w}}\in\mathcal{W}}\|\dot{\mathbf{H}}(\mathbf{w},\bar{\mathbf{w}};\boldsymbol{\theta}_0,\mathbf{t})f_{\mathbf{w}}(\bar{\mathbf{w}})\|\right]<\infty, \end{align*} $$
   $$ \begin{align*} \lim_{\tau\downarrow 0,(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})} \left\| \frac{\mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta}+\tau\mathbf{t},\mathbf{t})-\mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta},\mathbf{t})}{\tau}-\dot{\mathbf{H}}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta},\mathbf{t}) \right\| = 0, \end{align*} $$
   and, for some $\delta>0$ ,
   $$ \begin{align*} & \mathbb{E} \!\left[\! \sup_{ \tau\in (0,\delta),\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta, \mathbf{w}_2 \in\mathcal{W}} \!\left\| \!\frac{\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta}+\tau\mathbf{t},\mathbf{t})-\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t})}{\tau}-\dot{\mathbf{H}}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t}) \!\right\| \!\right]\\ & < \infty. \end{align*} $$

5.1 Partially Linear Regression Model

We take $\mathbf {s}=\mathbf {s}_{\mathtt {PLR}}$ and $\mathbf {H}=\mathbf {H}_{\mathtt {PLR}}$ , where

$$ \begin{align*} \mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{\partial}{\partial \boldsymbol{\theta}} m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\dot{\mathbf{x}}_{i,j} (\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}) \end{align*} $$

and

$$ \begin{align*} \mathbf{H}_{\mathtt{PLR}}(\mathbf{w}_i,\mathbf{w}_j) &= \frac{\partial^2}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}'} \mathbb{E}[ m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] = \frac{\partial}{\partial \boldsymbol{\theta}'} \mathbb{E}[ \mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] \\ &= \mathbb{E}[\dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' | \mathbf{w}_i,\mathbf{w}_j], \end{align*} $$

the latter depending on neither $\boldsymbol {\theta }$ nor $\mathbf {t}$ (because $m_{\mathtt {PLR}}(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol {\theta })$ is quadratic in $\boldsymbol {\theta }$ ).

Under mild conditions, Assumptions 1(iv), 1(v), and 2 hold with

$$ \begin{align*} \boldsymbol{\xi}_0(\mathbf{z}) &= -2\mathbb{E}[\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}) = 2 (\mathbf{x}-\mathbb{E}[\mathbf{x}|\mathbf{w}]) f_{\mathbf{w}}(\mathbf{w})\varepsilon, \end{align*} $$

$$ \begin{align*} \boldsymbol{\Xi}_0(\mathbf{w}) &= \mathbb{E}[\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}) \\ &= \mathbb{E}[ \dot{\mathbf{x}}_{1,2}\dot{\mathbf{x}}_{1,2}'(\varepsilon_1^2+\varepsilon_2^2)|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}), \end{align*} $$

and

$$ \begin{align*} \mathbf{G}_0(\mathbf{w}) &= \mathbf{H}_{\mathtt{PLR}}(\mathbf{w},\mathbf{w}) f_{\mathbf{w}}(\mathbf{w}) = 2 \mathbb{V}[\mathbf{x}|\mathbf{w}] f_{\mathbf{w}}(\mathbf{w}). \end{align*} $$

For instance, it suffices to set $b(\mathbf {z}) = (1+\|\mathbf {x}\|) (1+\|\mathbf {x}\|+|\gamma _0(\mathbf {w})|+|\varepsilon |)$ and to assume that

1. (i) The functions $\mathbf {w}\mapsto \gamma _0(\mathbf {w})$ , $\mathbf {w}\mapsto \mathbb {E}[\mathbf {x}|\mathbf {w}]$ , $\mathbf {w}\mapsto \mathbb {E}[\mathbf {x}\mathbf {x}'|\mathbf {w}]$ , $\mathbf {w}\mapsto \mathbb {E}[\varepsilon ^2|\mathbf {w}]$ , $\mathbf {w}\mapsto \mathbb {E}[\mathbf {x}\varepsilon ^2|\mathbf {w}]$ , and $\mathbf {w}\mapsto \mathbb {E}[\mathbf {x}\mathbf {x}'\varepsilon ^2|\mathbf {w}]$ are continuous on $\mathcal {W}$ .

2. (ii) $\mathbb {E}[(1+\|\mathbf {x}\|^4)\varepsilon ^4] + \sup _{\mathbf {w}\in \mathcal {W}}\mathbb {E}[(1+\|\mathbf {x}\|^4)\varepsilon ^4 |\mathbf {w} ] f_{\mathbf {w}}(\mathbf {w})<\infty $ and

   $$ \begin{align*} \mathbb{E}[(1 \!+\! \|\mathbf{x}\|^4)\gamma_0(\mathbf{w})^4 \!+\! \|\mathbf{x}\|^8] \!+\! \sup_{\mathbf{w}\in\mathcal{W}}\mathbb{E}[(1+\|\mathbf{x}\|^4)\gamma_0(\mathbf{w})^4 \!+\! \|\mathbf{x}\|^8|\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})\! < \!\infty. \end{align*} $$

3. (iii) With probability one, $\mathbb {V}[\mathbf {x}|\mathbf {w}]$ is positive definite and $\mathbb {V}[\varepsilon |\mathbf {x},\mathbf {w}]>0$ .

Under the additional assumptions of Theorem 1 and if $nh_n^d \to \infty $ , the pairwise difference estimator is asymptotically linear with influence function

$$ \begin{align*} \mathbf{z} \mapsto \mathbb{E}[\mathbb{V}[\mathbf{x}|\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})]^{-1} (\mathbf{x}-\mathbb{E}[\mathbf{x}|\mathbf{w}]) f_{\mathbf{w}}(\mathbf{w})\varepsilon. \end{align*} $$

Unless the distribution of $\mathbf {w}$ is uniform on $\mathcal {W}$ , the pairwise difference estimator is therefore asymptotically distinct from the estimators studied by Robinson (Reference Robinson1988) and Donald and Newey (Reference Donald and Newey1994), the influence function of these estimators being given by

$$ \begin{align*} \mathbf{z} \mapsto \mathbb{E}[\mathbb{V}[\mathbf{x}|\mathbf{w}]]^{-1} (\mathbf{x}-\mathbb{E}[\mathbf{x}|\mathbf{w}]) \varepsilon. \end{align*} $$

When $nh_n^d \to 2c<\infty $ , the pairwise difference estimator is still $\sqrt {n}$ -normal, but it ceases to be asymptotically linear. Similarly, the estimator of Donald and Newey (Reference Donald and Newey1994) is $\sqrt {n}$ -normal, but not asymptotically linear, when the associated tuning parameter is chosen appropriately (Cattaneo et al., Reference Cattaneo, Jansson and Newey2018), and in light of Linton (Reference Linton1995), it stands to reason that the same is true for the estimator of Robinson (Reference Robinson1988). Another $\sqrt {n}$ -normal estimator that is not asymptotically linear was proposed by Yatchew (Reference Yatchew1997). However, even under simplifying assumptions (e.g., uniformity of the distribution of $\mathbf {w}$ and/or conditional homoskedasticity of $\varepsilon $ ), it appears difficult to make insightful comparisons between the various estimators just mentioned.

5.2 Partially Linear Logit Model

We take $\mathbf {s}=\mathbf {s}_{\mathtt {PLL}}$ and $\mathbf {H}=\mathbf {H}_{\mathtt {PLL}}$ , where

$$ \begin{align*} \mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{\partial}{\partial \boldsymbol{\theta}} m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\dot{\mathbf{x}}_{i,j} (y_i - \Lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}))\mathbb{1}\{\dot{y}_{i,j}\neq0\} \end{align*} $$

and, defining $\lambda (u) = \partial \Lambda (u)/\partial u = \exp (u)/[1+\exp (u)]^2$ ,

$$ \begin{align*} \mathbf{H}_{\mathtt{PLL}}(\mathbf{w}_i,\mathbf{w}_j;\boldsymbol{\theta}) &= \frac{\partial^2}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}'} \mathbb{E}[ m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] = \frac{\partial}{\partial \boldsymbol{\theta}'} \mathbb{E}[ \mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] \\ &= \mathbb{E}[ \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta})\mathbb{1}\{\dot{y}_{i,j}\neq0\}| \mathbf{w}_i,\mathbf{w}_j], \end{align*} $$

where the latter does not depend on $\mathbf {t}$ (because $m_{\mathtt {PLL}}(\mathbf {z}_i,\mathbf {z}_j;\boldsymbol {\theta })$ is twice differentiable in $\boldsymbol {\theta }$ ).

Under mild conditions, Assumptions 1(iv), 1(v), and 2 hold with

$$ \begin{align*} \boldsymbol{\xi}_0(\mathbf{z}) &= -2\mathbb{E}[\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}), \end{align*} $$

$$ \begin{align*} \boldsymbol{\Xi}_0(\mathbf{w}) &= \mathbb{E}[\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}), \end{align*} $$

and

$$ \begin{align*} \mathbf{G}_0(\mathbf{w}) &= \mathbf{H}_{\mathtt{PLL}}(\mathbf{w},\mathbf{w};\boldsymbol{\theta}_0) f_{\mathbf{w}}(\mathbf{w}). \end{align*} $$

For instance, it suffices to set $b(\mathbf {z}) = 1+\|\mathbf {x}\|$ and to assume that, for some $\delta>0$ ,

1. (i) The function $\mathbf {w}\mapsto \gamma _0(\mathbf {w})$ is continuous on $\mathcal {W}$ . Also, the conditional distribution of $\mathbf {x}$ given $\mathbf {w}$ admits a density $f_{\mathbf {x}|\mathbf {w}}$ with respect to some measure $\rho $ such that $\mathbf {w}\mapsto f_{\mathbf {x}|\mathbf {w}}(\mathbf {x}|\mathbf {w})$ is continuous on $\mathcal {W}$ (with probability one) and

   $$ \begin{align*} \int_{\mathbb{R}^k} (1+ \|\mathbf{x}\|^2) \sup_{\|\mathbf{u}\|\leq\delta} f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}| \mathbf{w} + \mathbf{u})d\rho(\mathbf{x})<\infty \qquad \text{for every } \mathbf{w}\in\mathcal{W}. \end{align*} $$

2. (ii) $\mathbb {E}[\|\mathbf {x}\|^4] + \sup _{\mathbf {w}\in \mathcal {W}}\mathbb {E}[\|\mathbf {x}\|^4 |\mathbf {w}]f_{\mathbf {w}}(\mathbf {w})<\infty .$

3. (iii) With probability one, $\mathbb {V}[\mathbf {x}|\mathbf {w}]$ is positive definite.

5.3 Partially Linear Tobit Model

We take $\mathbf {s}=\mathbf {s}_{\mathtt {PLT}}$ and $\mathbf {H}=\mathbf {H}_{\mathtt {PLT}}$ , where

$$ \begin{align*} \mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\dot{\mathbf{x}}_{i,j} \left(\mathbb{1}\{y_i>\max(y_j+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta},0)\} - \mathbb{1}\{y_j>\max(y_i-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta},0)\} \right) \end{align*} $$

and

$$ \begin{align*} &\mathbf{H}_{\mathtt{PLT}}(\mathbf{w}_i,\mathbf{w}_j;\boldsymbol{\theta},\mathbf{t}) \\ &\quad = \mathbb{E}\left[\left. \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \left(\mathbb{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}>0\}+\mathbb{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}=0,\dot{\mathbf{x}}_{i,j}'\mathbf{t}\geq0\}\right)\eta_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})\right| \mathbf{w}_i,\mathbf{w}_j\right] \\ &\qquad + \mathbb{E}\left[\left. \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \left(\mathbb{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}<0\}+\mathbb{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}=0,\dot{\mathbf{x}}_{i,j}'\mathbf{t}<0\}\right)\eta_{\mathtt{PLT}}(\mathbf{z}_j,\mathbf{z}_i;\boldsymbol{\theta})\right| \mathbf{w}_i,\mathbf{w}_j\right], \end{align*} $$

with

$$ \begin{align*} &\eta_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) \\ &\quad = 2\int_0^{\infty} f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_i'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_i)+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}|\mathbf{w}_i) f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_j'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_j)|\mathbf{w}_j)d\varepsilon \\ &\qquad + f_{\varepsilon|\mathbf{w}}(-\mathbf{x}_i'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_i)+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}|\mathbf{w}_i) \int_{-\infty}^0 f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_j'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_j)|\mathbf{w}_j)d\varepsilon. \end{align*} $$

Under mild conditions, Assumptions 1(iv), 1(v), and 2 hold with

$$ \begin{align*} \boldsymbol{\xi}_0(\mathbf{z}) &= -2\mathbb{E}[\mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}), \end{align*} $$

$$ \begin{align*} \boldsymbol{\Xi}_0(\mathbf{w}) &= \mathbb{E}[\mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}), \end{align*} $$

and

$$ \begin{align*} \mathbf{G}_0(\mathbf{w}) &= \mathbf{H}_{\mathtt{PLT}}(\mathbf{w},\mathbf{w};\boldsymbol{\theta}_0,\boldsymbol{\theta}_0) f_{\mathbf{w}}(\mathbf{w}). \end{align*} $$

For instance, it suffices to set $b(\mathbf {z}) = 1+\|\mathbf {x}\|$ and to assume that, for some $\delta>0$ ,

1. (i) The function $\mathbf {w}\mapsto \gamma _0(\mathbf {w})$ is continuous on $\mathcal {W}$ . Also, the conditional distribution of $\mathbf {x}$ given $\mathbf {w}$ admits a density $f_{\mathbf {x}|\mathbf {w}}$ with respect to some measure $\rho $ such that $\mathbf {w}\mapsto f_{\mathbf {x}|\mathbf {w}}(\mathbf {x}|\mathbf {w})$ is continuous on $\mathcal {W}$ (with probability one) and

   $$ \begin{align*} \int_{\mathbb{R}^k} (1+ \|\mathbf{x}\|^2) \sup_{\|\mathbf{u}\|\leq\delta} f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}| \mathbf{w} + \mathbf{u})d\rho(\mathbf{x})<\infty \qquad \text{for every } \mathbf{w}\in\mathcal{W}. \end{align*} $$

   In addition, the function $(\varepsilon ,\mathbf {w})\mapsto f_{\varepsilon |\mathbf {w}}(\varepsilon |\mathbf {w})$ is continuous and bounded and the function

   $$ \begin{align*} (\mathbf{x},\mathbf{w})\mapsto \int_{\mathbb{R}} \sup_{|u|+\|\mathbf{u}\|\leq\delta}f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w})+u|\mathbf{w}+\mathbf{u})d\varepsilon \end{align*} $$
   is bounded.

2. (ii) $\mathbb {E}[\|\mathbf {x}\|^4] + \sup _{\mathbf {w}\in \mathcal {W}}(1+\mathbb {E}[\|\mathbf {x}\|^4 |\mathbf {w}])f_{\mathbf {w}}(\mathbf {w})<\infty .$

3. (iii) With probability one, $\mathbb {V}[\mathbf {x}|\mathbf {w}]$ is positive definite.

6.1 Simulation Design: Partially Linear Regression Model

We base the partially linear regression designs on equations (6.1) and (6.2) of Robinson (Reference Robinson1988). We consider

$$ \begin{align*} y_i = x_i + \gamma_0(\mathbf{w}_i) + \varepsilon_i. \end{align*} $$

The three designs are as follows. In Model 1, $x_i=2v_i$ , $v_i\sim \mathsf {N}(0,1)$ , $w_i=v_i+u_i$ , $u_i\sim \mathrm {N}(0,2)$ , $\gamma _0(w)=w^2+1$ , $\varepsilon _i\sim \mathsf {N}(0,1)$ , and $\{v_i,u_i,\varepsilon _i\}$ are jointly independent. In Model 2, $d=\dim (\mathbf {w})=2$ , $\gamma _0(\mathbf {w})=\mathbf {w}'\mathbf {w}+1$ , $(x_i,\mathbf {w}_i')'$ is normal with each component having mean $1$ and variance $3$ , and any pair in $(x_i,\mathbf {w}_i')'$ has correlation $2/3$ . Model 3 uses the same joint normal specification for $(x_i,\mathbf {w}_i')'$ as Model 2, but with $d=\dim (\mathbf {w})=3$ . In all models, $\theta _0=1$ . The regressor of interest $x_i$ is correlated with $\mathbf {w}_i$ , and thus, ignoring the $\gamma _0(\mathbf {w}_i)$ term will induce bias in the estimator of $\theta _0$ .

6.2 Simulation Design: Partially Linear Logit Model

We base the partially linear logit designs on those of Honoré and Powell (Reference Honoré, Powell, Andrews and Stock2005). As in the partially linear regression simulations, the three designs differ mainly in the dimension of $\mathbf {w}$ . We consider

$$ \begin{align*} y_i = \mathbb{1}\left\{ x_{1i} + x_{2i} + \gamma_0(\mathbf{w}_i) + \varepsilon_i\geq 0 \right\},\qquad \gamma_0(\mathbf{w}) = \mathbf{w}'\mathbf{w}-(1+\dim(\mathbf{w})), \end{align*} $$

where the CDF of $\varepsilon _i$ is $\Lambda , x_{2i}$ has a discrete distribution with $\mathbb {P}[x_{2i}=1]=1/2=\mathbb {P}[x_{2i}=-1]$ , $x_{1i}=v_i + \mathbf {w}_i'\mathbf {w}_i$ with $v_i\sim \mathsf {N}(0,1)$ , and $\{\varepsilon _i,\mathbf {w}_i,x_{2i},v_i\}$ are jointly independent. For Model d ( $d\in \{1,2,3\}$ ), $\dim (\mathbf {w})=d$ , and $\mathbf {w}_i$ has a normal distribution with each element having mean zero and variance one, equicorrelated with correlation $0.2$ . In all models, $\boldsymbol {\theta }_0=(1, 1)'$ .

As noted by Honoré and Powell (Reference Honoré, Powell, Andrews and Stock2005), ignoring the presence of $\gamma _0(\mathbf {w}_i)$ induces bias in the estimator of $\boldsymbol {\theta }_0$ , although the bias for the second element tends to be negligible relative to the standard error. We focus on constructing $95\%$ confidence intervals for the first element of $\boldsymbol {\theta }_0$ .

6.3 Simulation Results

Tables 1–3 display the coverage probabilities and average lengths of the four confidence intervals for the partially linear regression model, and Tables 4–6 for the partially linear logit model. Figure 2 complements the tables by plotting coverage probability as a function of bandwidth h for all six DGPs. For each data-generating process, the sample size was $n=2{,}000$ , we conducted $2{,}000$ simulation replications, and for each replication, we drew $2{,}000$ bootstrap samples to compute bootstrap quantiles.

Table 1

Bootstrap 95% confidence intervals for partially linear regression model: DGP 1 ( $d=1$ ).

Table 1 Long description

The table is structured with a primary header row dividing the data into two main groups: B equals 1 and B equals 3 super 1 all over d. Under each of these, a secondary header divides the data by L equals 0 and L equals 1. Each of these four sub-sections contains two columns: Coverage and Length.

The first column, h, lists bandwidth values from 0.35 to 1.50.

Key data trends:

* For B equals 1, Coverage remains stable near 0.95 across all h values for both L equals 0 and L equals 1, with Length staying constant at approximately 0.058 to 0.060.

* For B equals 3 super 1 all over d, Coverage and Length increase significantly as h increases. For L equals 0, Coverage starts at 0.956 at h equals 0.35 and reaches 1.000 by h equals 1.10, while Length grows from 0.058 to 0.119.

* For L equals 1 in the second group, Coverage starts at 0.956 and reaches 1.000 at h equals 1.30, with Length growing from 0.058 to 0.108.

Specific rows are highlighted: h equals 0.70 is highlighted in yellow, and h equals 1.00 is highlighted in orange, representing M S E-optimal bandwidths for L equals 0 and L equals 1 respectively.

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Table 2

Bootstrap 95% confidence intervals for partially linear regression model: DGP 2 ( $d=2$ ).

Table 2 Long description

The table is structured with a primary column for bandwidth h and two main super-columns for B = 1 and B = 3 super 1/d. Each super-column is subdivided into L = 0 and L = 1, which are further split into Coverage and Length metrics.

* Headers: The first column is h. The next four columns fall under B = 1 (subdivided into L = 0 and L = 1). The final four columns fall under B = 3 super 1/d (subdivided into L = 0 and L = 1).

* Data Trends: As h increases from 0.12 to 1.12:

* For B = 1, L = 0: Coverage drops significantly from 0.991 to 0.007, while Length decreases from 0.149 to 0.090.

* For B = 1, L = 1: Coverage remains relatively stable, starting at 0.994 and ending at 0.924, with Length decreasing from 0.171 to 0.088.

* For B = 3 super 1/d, L = 0: Coverage drops from 0.949 to 0.026, while Length stays relatively flat around 0.090 to 0.116.

* For B = 3 super 1/d, L = 1: Coverage remains very stable between 0.949 and 0.957, with Length decreasing from 0.122 to 0.098.

* Specific Values: At h = 0.30, Coverage for B = 1, L = 0 is 0.961 and for L = 1 is 0.978. For B = 3 super 1/d, L = 0 is 0.931 and L = 1 is 0.956.

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Table 3

Bootstrap 95% confidence intervals for partially linear regression model: DGP 3 ( $d=3$ ).

Table 3 Long description

The table displays Coverage and Length metrics for two primary groups: B = 1 and B = 3 super 1/d. Each group is further subdivided into L = 0 and L = 1.

* Header Structure: The first column is h (bandwidth). The subsequent columns are paired as Coverage and Length for the following four conditions:

1. B = 1, L = 0

2. B = 1, L = 1

3. B = 3 super 1/d, L = 0

4. B = 3 super 1/d, L = 1

* Data Trends:

As h increases from 0.20 to 1.50:

- For B = 1, L = 0: Coverage drops sharply from 0.998 to 0.000, while Length decreases from 0.325 to 0.118.

- For B = 1, L = 1: Coverage remains high (above 0.950) until h = 1.40, where it begins to decline to 0.887.

- For B = 3 super 1/d, L = 0: Coverage drops from 0.950 to 0.000, with Length decreasing from 0.210 to 0.156.

- For B = 3 super 1/d, L = 1: Coverage remains stable and high, ranging between 0.950 and 0.970 across all values of h, while Length decreases from 0.254 to 0.126.

Row h = 0.40 is highlighted in yellow, and row h = 1.00 is highlighted in orange, representing M S E-optimal bandwidths for L = 0 and L = 1 respectively.

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Table 4

Bootstrap 95% confidence intervals for partially linear logit model: DGP 1 ( $d=1$ ).

Table 4 Long description

The table presents Coverage and Length metrics for two main scenarios: B equals 1 and B equals 3 super 1/d. Each scenario is further divided into L equals 0 and L equals 1.

* Column Structure: The first column is h (bandwidth). The subsequent columns are paired as Coverage and Length for each combination of B and L.

* Data Trends: As h increases from 0.10 to 0.54, the Length of the confidence intervals generally decreases across all categories. For example, at B equals 1 and L equals 0, Length drops from 0.319 to 0.302. Coverage remains relatively stable near 0.95 for lower h values but shows a slight decline as h increases, particularly in the B equals 3 super 1/d and L equals 0 category, which drops to 0.916 at h equals 0.54.

* Specific Highlights:

- Row h equals 0.20 is highlighted in yellow, representing the M S E-optimal bandwidth for L equals 0, with coverage values of 0.950 and 0.943.

- Row h equals 0.36 is highlighted in orange, representing the M S E-optimal bandwidth for L equals 1, with coverage values of 0.951 and 0.946.

* Simulation Parameters: Based on 2,000 replications with a sample size n equals 2,000 and 2,000 bootstrap draws.

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Table 5

Bootstrap 95% confidence intervals for partially linear logit model: DGP 2 ( $d=2$ ).

Table 5 Long description

The table is structured with a primary header row dividing the data into two main groups: B = 1 and B = 3 super 1/d. Below these, a secondary header divides each group by L = 0 and L = 1. Each of these four sub-sections contains two columns: Coverage and Length.

The leftmost column is h, representing bandwidth, with 22 rows ranging from 0.14 to 0.68.

Key data trends:

* For B = 1 and L = 0: As h increases from 0.14 to 0.68, Coverage decreases from 0.982 to 0.858, and Length decreases from 0.507 to 0.338.

* For B = 1 and L = 1: Coverage decreases from 0.985 to 0.937, and Length decreases from 0.568 to 0.350.

* For B = 3 super 1/d and L = 0: Coverage decreases from 0.952 to 0.832, and Length decreases from 0.404 to 0.318.

* For B = 3 super 1/d and L = 1: Coverage remains relatively stable, starting at 0.952 and ending at 0.921, while Length decreases from 0.431 to 0.329.

Specific rows are highlighted: h = 0.22 and h = 0.36 are highlighted in yellow (corresponding to M S E-optimal bandwidths for L = 0), and h = 0.28 and h = 0.45 are highlighted in orange (corresponding to M S E-optimal bandwidths for L = 1).

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Table 6

Bootstrap 95% confidence intervals for partially linear logit model: DGP 3 ( $d=3$ ).

Table 6 Long description

The table is structured with a primary header row for B = 1 and B = 3 super 1 all over d. Under each B value, there are sub-headers for L = 0 and L = 1. Each of these four categories contains two columns: Coverage and Length. The leftmost column is the bandwidth h, ranging from 0.20 to 0.75.

Key data trends include:

* For B = 1, L = 0: Coverage starts at 0.996 at h = 0.20 and decreases to 0.720 at h = 0.75. Length decreases from 1.110 to 0.380.

* For B = 1, L = 1: Coverage starts at 0.995 and decreases to 0.902. Length decreases from 1.370 to 0.403.

* For B = 3 super 1 all over d, L = 0: Coverage starts at 0.953 and decreases to 0.663. Length decreases from 0.693 to 0.353.

* For B = 3 super 1 all over d, L = 1: Coverage starts at 0.944 and decreases to 0.855. Length decreases from 0.826 to 0.368.

Specific rows are highlighted: h = 0.39 is highlighted in yellow for L = 0 M S E-optimal bandwidths, and h = 0.50 is highlighted in orange for L = 1 M S E-optimal bandwidths.

Note: The table shows coverage probability and average length of confidence intervals over $2{,}000$ simulation replications. The sample size is $n=2{,}000$ and the number of bootstrap draws is $2{,}000$ . The rows highlighted in yellow and orange correspond to the estimated MSE-optimal bandwidths for the point estimators with $L=0$ and $L=1$ , respectively.

Figure 2

Coverage probabilities as a function of bandwidth h.

Note: Each panel plots the empirical coverage probability of $\mathsf {CI}_{n,0.95}^{\hspace{0.5pt}*}(B,L)$ against bandwidth h for the four procedures $(B,L)\in \{(1,0),(1,1),(3^{1/d},0),(3^{1/d},1)\}$ . “PLR” stands for the partially linear regression model, and “PLL” stands for the partially linear logit model. d denotes the dimension of $\mathbf {w}$ . The horizontal dotted line marks the $95\%$ nominal level. The two vertical lines mark the estimated MSE-optimal bandwidths $h_{L=0}$ and $h_{L=1}$ for the point estimators $\widehat {\boldsymbol {\theta }}_n$ ( $L=0$ ) and $\widetilde {\boldsymbol {\theta }}_n$ ( $L=1$ ), respectively. The sample size is $n=2{,}000$ and results are based on $2{,}000$ simulation replications with $2{,}000$ bootstrap samples.

Figure 2 Long description

A multi-panel figure with six line graphs labeled a through f. Each graph plots Coverage on the y-axis from 0.0 to 1.0 against bandwidth h on the x-axis. A horizontal dotted line at 0.95 represents the nominal level. Two vertical dashed lines mark M S E-optimal bandwidths h sub L equals 0 and h sub L equals 1.

* Panel a, P L R d equals 1: All four lines remain stable near the 0.95 level across h from 0.4 to 1.4.

* Panel b, P L L d equals 1: All four lines are tightly clustered at the 0.95 level across h from 0.1 to 0.5.

* Panel c, P L R d equals 2: The lines for C I sub n, 1 minus alpha (1, 0) and (3 to the 1 over d power, 0) show a sharp sigmoidal decline from 0.95 toward 0.0 as h increases beyond 0.4. The other two lines remain stable.

* Panel d, P L L d equals 2: Similar to panel c, but the decline is much more gradual, staying above 0.8.

* Panel e, P L R d equals 3: The decline for the two sensitive lines is even steeper than in d equals 2, reaching 0.0 by h equals 1.2.

* Panel f, P L L d equals 3: The decline for the two sensitive lines is more pronounced than in d equals 2 but remains significantly higher than the P L R counterpart, ending near 0.7.

In all panels, the legend identifies four procedures: C I sub n, 1 minus alpha (1, 0) as a brown dashed line, (1, 1) as a light blue dashed line, (3 to the 1 over d power, 0) as a brown dotted line, and (3 to the 1 over d power, 1) as a solid dark blue line.

The panel headed “ $B=1$ ” refers to $\mathsf {CI}_{n,0.95}^{\hspace{0.5pt}*}(1,L)$ , in which the bandwidth used to compute bootstrapped quantiles equals the bandwidth used for estimation. The panel headed “ $B=3^{1/d}$ ” refers to our proposed methods $\mathsf {CI}_{n,0.95}^{\hspace{0.5pt}*}(3^{1/d},L)$ , in which the bandwidth for bootstrap estimation is rescaled by the factor $3^{1/d}$ . Within each panel, $L=0$ is based on the estimator $\widehat {\boldsymbol {\theta }}_n$ and $L=1$ is based on the debiased estimator $\widetilde {\boldsymbol {\theta }}_n$ with $\mathbf {c}=(1, 2)'$ .

To form the bandwidth grid, we first estimate $h_{L=0}$ and $h_{L=1}$ by auxiliary simulations. These are the MSE-minimizing bandwidths for $\widehat {\boldsymbol {\theta }}_n$ and $\widetilde {\boldsymbol {\theta }}_n$ , respectively. The grid is then taken as the union of $h_{L=0}\times \{0.5,\dots ,1.5\}$ and $h_{L=1}\times \{0.5,\dots ,1.5\}$ (increments of $0.1$ ). Rows highlighted in yellow and orange in the tables mark $h_{L=0}$ and $h_{L=1}$ , respectively.

The results are broadly consistent with our theoretical predictions. We focus on two forces of particular interest: small-bandwidth regions, where the standard bootstrap variance may be inflated, and large-bandwidth regions, where smoothing bias may become non-negligible. The simulations also include one instructive qualification: in the partially linear regression model with $d=1$ , the conventional procedures are already very well calibrated, while our proposed procedure $\mathsf {CI}_{n,0.95}^{\hspace{0.5pt}*}(3^{1/d},1)$ is conservative. This is a specialized design in which the smoothing bias targeted by jackknife debiasing is negligible, so bandwidth rescaling and debiasing mostly widen the intervals without offsetting a meaningful bias. This qualification is specific to the partially linear regression $d=1$ design and does not describe the partially linear logit model with $d=1$ , where uncorrected intervals exhibit undercoverage as the bandwidth grows and bias correction is beneficial.