Large Sample Properties of Partitioning-Based Series Estimators

Large Sample Properties of Partitioning-Based Series EstimatorsSupplemental Appendix

Setup and Assumptions

For a $d$-tuple $\mathbf{q}=(q_1,\ldots,q_d)\in\mathbb{Z}^d_+$, define $[\mathbf{q}]=\sum_{\ell=1}^{d}q_\ell$, $\mathbf{x}^\mathbf{q}=x_1^{q_1}x_2^{q_2}\cdots x_d^{q_d}$, $\mathbf{q}!=q_1!\cdots q_d!$ and $\partial^\mathbf{q} \mu(\mathbf{x})=\partial^{[\mathbf{q}]} \mu(\mathbf{x})/\partial x_1^{q_1}\ldots\partial x_d^{q_d}$. Unless explicitly stated otherwise, whenever $\mathbf{x}$ is a boundary point of some closed set, the partial derivative is understood as the limit with $\mathbf{x}$ ranging within it. Let $\mathbf{0}=(0,\cdots,0)'$ and $\mathbf{1}=(1,\cdots, 1)'$ be the length-$d$ vectors of zeros and ones respectively. The tensor product or Kronecker product operator is $\otimes$, and the entrywise division operator (Hadamard division) is $\oslash$. The smallest integer greater than or equal to $u$ is $\lceil u \rceil$.

We use several norms. For a column vector $\mathbf{v}=(v_1,\ldots,v_M)'\in\mathbb{R}^M$, we let $\|\mathbf{v}\|=(\sum_{i=1}^M v_i^2)^{1/2}$, $\|\mathbf{v}\|_\infty=\max_{1\leq i\leq M}|v_i|$ and $\dim(\mathbf{v})=M$. For a matrix $\mathbf{A}\in\mathbb{R}^{M\times N}$, $\|\mathbf{A}\|_1=\max_{1\leq j \leq N}\sum_{i=1}^{M}|a_{ij}|$, $\|\mathbf{A}\|=\max_i\sigma_i(\mathbf{A})$ and $\|\mathbf{A}\|_\infty=\max_{1\leq i\leq M}\sum_{j=1}^N|a_{ij}|$ for operator norms induced by $L_1$, $L_2$, and $L_\infty$ norms respectively, where $\sigma_i(\mathbf{A})$ is the $i$-th singular value of $\mathbf{A}$. $\lambda_{\max}(\mathbf{A})$ and $\lambda_{\min}(\mathbf{A})$ denote the maximum and minimum eigenvalues of $\mathbf{A}$. For a real-valued function $g(\mathbf{x})$, $\|g\|_{L_p(\mathcal{X})}=(\int_\mathcal{X}|g(x)|^p\,d\mathbf{x})^{1/p}$ and $\|g\|_{L_\infty(\mathcal{X})}=\operatorname*{ess\,sup}_{\mathbf{x}\in\mathcal{X}}|g(\mathbf{x})|$, and let $\mathcal{C}^s(\mathcal{X})$ denote the space of $s$-times continuously differentiable functions on $\mathcal{X}$.

We also use standard empirical process notation: $\mathbb{E}_n[g]=\mathbb{E}_n[g(\mathbf{x}_i)]=\frac{1}{n}\sum_{i=1}^ng(\mathbf{x}_i)$, $\mathbb{E}[g]=\mathbb{E}[g(\mathbf{x}_i)]$, and $\mathbb{G}_n[g]=\mathbb{G}_n[g(\mathbf{x}_i)]=\frac{1}{\sqrt{n}}\sum_{i=1}^n(g(\mathbf{x}_i)-\mathbb{E}[g(\mathbf{x}_i)])$. For sequences of numbers or random variables, we use $a_n \lesssim b_n$ to denote $\limsup_n |a_n/b_n|$ is finite, and $a_n \lesssim_\mathbb{P} b_n$ or $a_n=O_\mathbb{P}(b_n)$ to denote $\limsup_{\epsilon\rightarrow\infty}\limsup_n\mathbb{P}[|a_n/b_n|\geq\epsilon]=0$. $a_n=o(b_n)$ implies $a_n/b_n\rightarrow 0$, and $a_n=o_\mathbb{P}(b_n)$ implies that $a_n/b_n\rightarrow_\mathbb{P} 0$ where $\rightarrow_\mathbb{P}$ denotes convergence in probability. $a_n\asymp b_n$ implies that $a_n\lesssim b_n$ and $b_n\lesssim a_n$. $\rightsquigarrow$ denotes convergence in distribution, and for two random variables $X$ and $Y$, $X=_d Y$ implies that they have the same probability distribution.

We set $\mu(\mathbf{x}):=\partial^\mathbf{0}\mu(\mathbf{x})$ and likewise $\widehat{\mu}_j(\mathbf{x}):=\widehat{\partial^\mathbf{0}\mu}_j(\mathbf{x})$ for $j=0,1,2,3$. Let $\mathbf{X}=[\mathbf{x}_1, \ldots, \mathbf{x}_n]'$. Finally, let $C, C_1, C_2, \ldots$ denote universal constants which may be different in different uses.

Setup

As a complete reference, this section repeats the setup, assumptions, and notation used in the main paper.

We study the standard nonparametric regression setup, where $\{(y_i,\mathbf{x}_i'), i = 1, \ldots n\}$ is a random sample from the model

for a scalar response $y_i$ and a $d$-vector of continuously distributed covariates $\mathbf{x}_i = (x_{1,i}, \ldots, x_{d,i})'$ with compact support $\mathcal{X}$. The object of interest is the unknown regression function $\mu(\cdot)$ and its derivatives.

We focus on partitioning-based, or locally-supported, series (linear sieve) least squares regression estimators. Let $\Delta=\{\delta_l: 1\leq l \leq \bar{\kappa}\}$ be a collection of open and disjoint subsets of $\mathcal{X}$ such that the closure of their union is $\mathcal{X}$. Every $\delta_l\in\Delta$ is restricted to be polyhedral. An important special case is tensor-product partitions, in which case the support of the regressors is of tensor product form and each dimension of $\mathcal{X}$ is partitioned marginally into intervals.

Based on a particular partition, the dictionary of $K$ basis functions, each of order $m$ (e.g., $m=4$ for cubic splines) is denoted by \[\mathbf{x}_i \mapsto \mathbf{p}(\mathbf{x}_i) := \mathbf{p}(\mathbf{x}_i;\Delta,m)=(p_1(\mathbf{x}_i;\Delta,m),\cdots,p_K(\mathbf{x}_i;\Delta,m))'.\] For a point $\mathbf{x} \in \mathcal{X}$ and $\mathbf{q} = (q_1, \ldots, q_d)'\in\mathbb{Z}_+^d$, the partial derivative $\partial^\mathbf{q}\mu(\mathbf{x})$ is estimated by least squares regression \[\widehat{\partial^\mathbf{q}\mu}(\mathbf{x})=\partial^\mathbf{q} \mathbf{p}(\mathbf{x})'\widehat{\boldsymbol{\beta}},\qquad\qquad \widehat{\boldsymbol{\beta}}\in\underset{\mathbf{b}\in\mathbb{R}^K}{\operatorname*{arg\,min}}\,\sum_{i=1}^{n}(y_i - \mathbf{p}(\mathbf{x}_i)'\mathbf{b})^2,\] Our assumptions will guarantee that the sample matrix $\mathbb{E}_n[\mathbf{p}(\mathbf{x}_i) \mathbf{p}(\mathbf{x}_i)']$ is nonsingular with probability approaching one in large samples, and thus we write the estimator as

The order $m$ of the basis is usually fixed in practice, and thus the tuning parameter for this class of nonparametric estimators is $\Delta$. As $n\rightarrow\infty$, $\bar{\kappa}\rightarrow\infty$, and the volume of each $\delta_l$ shrinks proportionally to $h^d$, where $h=\max\{\operatorname*{diam}(\delta): \delta\in\Delta\}$.

Assumptions

We list our main assumptions for completeness. Detailed discussion can be found in the main paper. The first assumption concerns the data generating process.

The next assumption requires that the basis is locally supported. We employ the notion of active basis: a function $p(\cdot)$ on $\mathcal{X}$ is active on $\delta\in\Delta$ if it is not identically zero on $\delta$.

We remind readers that Assumption (ref) and (ref) implicitly relate the number of approximating series terms, the number of cells in $\Delta$, and the maximum mesh size: $K\asymp \bar{\kappa}\asymp h^{-d}$.

The next assumption gives an explicit high-level expression of the leading approximation error which is needed for bias correction and integrated mean squared error (IMSE) expansion. To simplify notation, for each $\mathbf{x}\in\mathcal{X}$ define $\delta_{\mathbf{x}}$ as the element of $\Delta$ whose closure contains $\mathbf{x}$ and $h_\mathbf{x}$ for the diameter of this $\delta_{\mathbf{x}}$.

Our last assumption concerns the basis used for bias correction. Specifically, for some $\tilde{m} > m$, let $\bm{\tilde{\mathbf{p}}}(\mathbf{x}) := \bm{\tilde{\mathbf{p}}}(\mathbf{x}; \tilde{\Delta}, \tilde{m})$ be a basis of order $\tilde{m}$ defined on a partition $\tilde{\Delta}$ with maximum mesh $\tilde{h}$. Objects accented with a tilde always pertain to this secondary basis and partition for bias correction. See the main text for complete discussion.

Technical Lemmas

We present a series of technical lemmas that will be used for pointwise and uniform analysis. To begin with, we introduce additional notation to simplify our expressions. The following are frequently used outer-product matrices:

In the main paper, we define four partitioning-based nonparametric estimators, numbered as $j=0,1,2,3$, which can be written in exactly the same form: \[ \widehat{\partial^\mathbf{q}\mu}_j(\mathbf{x}) := \widehat{\bm{\gamma}}_{\mathbf{q},j}(\mathbf{x})'\mathbb{E}_n[\bm{\bm{\Pi}}_{j}(\mathbf{x}_i) y_i ]. \] We summarize the definitions of these quantities here. For $j=0$ (classical estimator without bias correction),

For $j=1$ (high-order-basis bias correction),

For $j=2$ (least squares bias correction),

For $j=3$ (plug-in bias correction),

For $j=0, 1, 2, 3$, $\boldsymbol{\gamma}_{\mathbf{q}, j}(\mathbf{x})$ is defined as $\widehat{\bm{\gamma}}_{\mathbf{q},j}(\mathbf{x})$ in (ref), (ref), (ref) and (ref) but with sample averages replaced by their population counterparts. Finally, define

The next lemma establishes the rate of convergence and boundedness for $\widehat{\mathbf{Q}}_m$ and other relevant matrices, which will be crucial for both pointwise and uniform analysis. Since the orders of bases used to construct bias corrected estimators are fixed and mesh size ratio $\rho\rightarrow \rho_0\in(0,\infty)$ by Assumption (ref), the same conclusions also apply to $\widehat{\mathbf{Q}}_{\tilde{m}}$, $\mathbf{Q}_{\tilde{m}}$ and other related quantities, though we do not make it explicit in the statement.

The next lemma shows that the asymptotic error expansion specified in Assumption (ref) translates into the bias expansion for the classical point estimator $\widehat{\partial^\mathbf{q}\mu}_0(\cdot)$, which is (partly) reported in Lemma 3.1 of the main paper. Importantly, the following lemma also gives a uniform bound on the conditional bias which is crucial for least squares bias correction ($j=2$) and uniform inference.

Equation (ref) implies that the leading approximation error is approximately orthogonal to $\mathbf{p}(\cdot)$, in which case the $L_\infty$ approximation error coincides with the leading smoothing bias of $\widehat{\partial^\mathbf{q}\mu}_0(\cdot)$.

The next lemma shows that $\widehat{\bm{\gamma}}_{\mathbf{q},j}(\mathbf{x})$ in (ref), (ref), (ref), and (ref) converge to their population counterpart $\boldsymbol{\gamma}_{\mathbf{q}, j}(\mathbf{x})$ in a proper sense, and all the $\boldsymbol{\gamma}_{\mathbf{q}, j}(\mathbf{x})$ are sufficiently controlled in terms of both $L_\infty$- and $L_2$-operator norms.

The last lemma in this section proves that the asymptotic variance of classical and bias-corrected estimators is properly bounded. Importantly, we need not only an upper bound, but also a lower bound since the variance term appears in the denominator of $t$-statistics.

This lemma also implies that bias correction does not change the asymptotic order of the variance.

Integrated Mean Squared Error

We now state an integrated mean squared error (IMSE) result that specializes Theorem 4.1 in the main paper to the case of general tensor-product partitions but is less restricted than Theorem 4.2. In particular, we allow for any $\mathbf{q}$ and do not rely on Assumption 6(c). This result gives the key ingredient for characterizing the limit $\mathscr{V}_{\Delta,\mathbf{q}}$ and $\mathscr{B}_{\Delta,\mathbf{q}}$ in some generality.

To be specific, throughout this section, we assume the support of the regressors is of tensor-product form, that each dimension of $\mathcal{X}$ is partitioned marginally into intervals, and $\Delta$ is the tensor product of these intervals. Let $\mathcal{X}_\ell = [\underline{x}_\ell, \overline{x}_\ell]$ be the support of covariate $\ell = 1, 2, \ldots, d$, and partition this into $\kappa_\ell$ disjoint subintervals defined by $\{\underline{x}_\ell=t_{\ell,0}<t_{\ell,1}<\cdots<t_{\ell,\kappa_\ell-1}<t_{\ell,\kappa_\ell}=\overline{x}_\ell\}$. If this partition of $\mathcal{X}_\ell$ is $\Delta_\ell$, then a complete partition of $\mathcal{X}$ can be formed by tensor products of the one-dimensional partitions: $\Delta=\otimes_{\ell=1}^d\Delta_\ell$, with $\boldsymbol{\kappa}=(\kappa_1,\kappa_2,\dots,\kappa_d)'$ collecting the number of subintervals in each dimension of $\mathbf{x}_i$, and $\bar{\kappa} = \kappa_1\kappa_2\cdots\kappa_d$. A generic cell of this partition is the rectangle

Assumption (ref) can be verified by choosing the knot positions appropriately, often dividing $\mathcal{X}_\ell$ uniformly or by quantiles.

As before, $\delta_{\mathbf{x}}$ is the subrectangle in $\Delta$ containing $\mathbf{x}$, and we write $\mathbf{b}_{\mathbf{x}}$ for the vector collecting the interval lengths of $\delta_{\mathbf{x}}$ (as defined in Assumption 6). $\mathbf{t}_{\mathbf{x}}^L$ denotes the start point of $\delta_{\mathbf{x}}$. For $\ell=1, \ldots, d$, for a generic cell $\delta_{l_1\ldots l_d}$ as in (ref), we write $b_{\ell, l_\ell}=t_{\ell, l_\ell+1} - t_{\ell, l_\ell}$, $l_\ell=0, \ldots, \kappa_\ell-1$, $b_\ell=\max_{0\leq l_\ell\leq\kappa_\ell-1} b_{\ell,l_\ell}$ and $b=\max_{1\leq \ell \leq d} b_\ell$ ($b\asymp h$ by Assumption (ref)).

This theorem takes advantage of the assumed tensor-product structure of the partition $\Delta$ to express the leading bias and variance as proportional to the number of subintervals used for each regressor. The accompanying constants are expressed as sums over the local contributions of the basis function $\mathbf{p}(\cdot)$ used, and are easy to be shown bounded in general (they may still not converge to any well-defined limit at this level of generality).

Assumption 6(a) used in Theorem (ref) (and in Theorem 4.2) slightly strengthens the quasi-uniform condition imposed in Assumption (ref). It can be verified by choosing knot positions appropriately. Specifically, let $\Delta$ be a tensor-product partition (see Equation (ref)) with (marginal) knots chosen as \[t_{\ell,l} = G_\ell^{-1}\left(\frac{l}{\kappa_\ell}\right), \qquad l=0,1,\dots,\kappa_\ell,\qquad \ell=1,2,\dots,d,\] where $G_\ell(\cdot)$ is a univariate continuously differentiable distribution function and $G_\ell^{-1}(v)=\inf\{x\in\mathbb{R}:G_\ell(x)\geq v\}$. In this case the function $g_\ell(\cdot)$ in Assumption 6(a) is simply the density of $G_\ell(\cdot)$. Two examples commonly used in practice are: (i) evenly-spaced partitions, denoted by $\Delta_\mathtt{ES}$, where $G_\ell(x)=x$ and $g_\ell(x)=1$, and (ii) quantile-spaced partitions, denoted by $\Delta_\mathtt{QS}$, where $G_\ell(x)=\widehat{F}_\ell(x)$ with $\widehat{F}_\ell(x)$ the empirical distribution function for the $\ell$-th covariate, $\ell=1,2,\dots,d$. For the case of quantile-spaced partitioning, if $\widehat{F}_\ell(x)$ convergences to $F_\ell(x)=\mathbb{P}[x_\ell\leq x]$ in a suitable sense, $g_\ell(x)=dF_\ell(x)/dx$, i.e., the marginal density of $x_\ell$. See Agarwal-Studden_1980_AoS and Barrow-Smith_1978_QAM for a slightly more high-level condition in the context of univariate $B$-splines.

Assumption 6(b) in Theorem (ref) (and in Theorem 4.2) involves a scaling factor $\frac{h_\mathbf{x}^{2m-2[\mathbf{q}]}}{\mathbf{b}_\mathbf{x}^{\bm{u}_1+\bm{u}_2-2[\mathbf{q}]}}$. In Section (ref) we do not add specific restrictions on the shape of cells in $\Delta$ and thus the diameter of a cell is used to conveniently express the order of approximation error (denoted by $h$), but for tensor-product partitions the approximation error is usually characterized by the lengths of intervals on each axis (denoted by $\mathbf{b}_\mathbf{x}$), as illustrated in Section (ref). Therefore, the scaling factor makes the above theorem immediately apply to the three leading examples in Section (ref) based on tensor-product partitioning schemes (splines, wavelets, and piecewise polynomials).

Using Theorem (ref), the IMSE-optimal tuning parameter selector for partitioning-based series estimators becomes \[\boldsymbol{\kappa}_{\mathtt{IMSE}, \mathbf{q}} = \operatorname*{arg\,min}_{\boldsymbol{\kappa}\in\mathbb{Z}^d_{++}} \left\{ \frac{1}{n} \mathscr{V}_{\boldsymbol{\kappa},\mathbf{q}} + \mathscr{B}_{\boldsymbol{\kappa},\mathbf{q}} \right\}, \] which is still given in implicit form because the leading constants are not known to converge at this level of generality. It follows that $\boldsymbol{\kappa}_{\mathtt{IMSE},\mathbf{q}}=(\kappa_{\mathtt{IMSE},\mathbf{q},1},\dots,\kappa_{\mathtt{IMSE},\mathbf{q},d})'$ with $\kappa_{\mathtt{IMSE},\mathbf{q},\ell} \asymp n^{\frac{1}{2m+d}}$, $\ell=1,2,\cdots,d$.

Rates of Convergence

In this section, we discuss rates of convergence for the classical estimator and three bias-corrected estimators. Some results, such as the $L_2$ and uniform convergence rates for $\widehat{\partial^{\mathbf{q}}\mu}_0(\mathbf{x})$, are also reported in Theorem 4.3 of the main paper.

To begin with, we first report pointwise linearization which is an important intermediate step towards pointwise asymptotic normality.

We can use Lemma (ref) to establish the limiting distribution of the standardized $t$-statistics as given in Theorem 5.1 of the main paper. The proof is available in Section (ref).

For uniform convergence, we need to upgrade the above lemma to a uniform version at the cost of some stronger conditions:

Given Lemma (ref) and (ref), we are ready to derive the desired rates of convergence.

Relying on the uniform convergence results, we can further show the consistency of variance estimates $\widehat{\Omega}_j(\mathbf{x})$ defined in Section (ref). The results in the following theorem are partly reported in Theorem 5.2 in the main paper.

Uniform Inference

Strong Approximation

Now we move on to uniform inference. The following $t$-statistic processes are of interest: \[ \widehat{T}_j(\mathbf{x})=\frac{\widehat{\partial^\mathbf{q}\mu}_j(\mathbf{x})- \mu(\mathbf{x})}{\sqrt{\widehat{\Omega}_j(\mathbf{x})/n}}, \quad \mathbf{x}\in\mathcal{X},\quad j=0,1,2,3. \]

Let $r_n$ be a positive non-vanishing sequence, which will be used to denote the approximation error rate in the following analysis. As the first step, we employ Lemma (ref) and Theorem (ref) to show that the sampling and estimation uncertainty of $\widehat{\bm{\gamma}}_{\mathbf{q},j}(\mathbf{x})$ and $\widehat{\Omega}_j(\mathbf{x})$ are negligible uniformly over $\mathbf{x}\in\mathcal{X}$. Specifically, define \[t_j(\mathbf{x}) = \frac{\boldsymbol{\gamma}_{\mathbf{q},j}(\mathbf{x})'}{\sqrt{\Omega_j(\mathbf{x})}}\mathbb{G}_n[ \bm{\bm{\Pi}}_{j}(\mathbf{x}_i)\varepsilon_i], \quad \mathbf{x}\in\mathcal{X}, \quad j=0,1,2,3. \] The following lemma, which appears as Lemma 6.1 in the main paper, shows that $\widehat{T}_j(\cdot)$ can be approximated by $t_j(\cdot)$ uniformly.

The proof can be found in Section 8.2 of the main paper. Now our task reduces to construct valid distributional approximation to $t_j(\cdot)$ for each $j=0,1,2,3$, in a proper sense. Specifically, we want to show that on a sufficiently rich enough probability space there exists a copy $t'_j(\cdot)$ of $t_j(\cdot)$ and a Gaussian process $Z_j(\cdot)$ such that $\sup_{\mathbf{x}\in\mathcal{X}}|t'_j(\mathbf{x})-Z_j(\mathbf{x})|=o_\mathbb{P}(r_n^{-1})$. When such a construction is possible, we can use $Z_j(\cdot)$ to approximate the distribution of $t_j(\cdot)$, as well as $\widehat{T}_j(\cdot)$ in view of Lemma (ref). To save our notation, we denote this strong approximation by $\widehat{T}_j(\cdot)=_d Z_j(\cdot)+o_\mathbb{P}(r_n^{-1})$ in $\mathcal{L}^\infty(\mathcal{X})$ where $\mathcal{L}^\infty(\mathcal{X})$ refers to the set of all uniformly bounded real functions on $\mathcal{X}$ equipped with uniform norm.

We will show in the following that there exists an unconditional Gaussian approximating process \[Z_j(\mathbf{x})=\frac{\boldsymbol{\gamma}_{\mathbf{q},j}(\mathbf{x})'\bm{\bm{\Sigma}}_{j}^{1/2}} {\sqrt{\Omega_j(\mathbf{x})}}\mathbf{N}_{K_j}, \quad \mathbf{x}\in\mathcal{X}, \quad \text{for}\quad j=0,1,2,3, \] where $\mathbf{N}_{K_j}$ is a $K_j$-dimensional standard Normal vector with $K_j=\dim(\bm{\bm{\Pi}}_{j}(\cdot))$. As discussed in the main paper, we may have two strategies to construct strong approximations. The next theorem employs KMT coupling techniques.

Unfortunately, extending this to general multi-dimensional cases needs some non-trivial results which are not available in current strong approximation literature, to the best of our knowledge. Thus, when $d\geq 2$, we employ an improved version of Yurinskii's inequality due to Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, which allows us to replace the Euclidean norm in the original version of Yurinskii's inequality by a sup-norm. Since basis functions considered in this paper are locally supported, a metric of distance between random vectors in terms of sup-norm leads to a weaker rate restriction than that previously used in literature, though still stronger than what we obtained in Theorem (ref). The result is reported in Theorem 6.2 of the main paper. We restate it in the following for completeness. The proof can be found in Section 8.4 of the main paper.

Implementation

$\{Z_j(\mathbf{x}):\mathbf{x}\in\mathcal{X}\}$ is still an infeasible approximation of $\{\widehat{T}_j(\mathbf{x}):\mathbf{x}\in\mathcal{X}\}$. Our next objective is to construct practicably feasible distributional approximations. To be precise, we construct Gaussian processes $\{\widehat{Z}_j(\mathbf{x}):\mathbf{x}\in\mathcal{X}\}$, with distributions known conditional on the data $(\mathbf{y}, \mathbf{X})$, and there exists a copy $\widehat{Z}'_j(\cdot)$ of $\widehat{Z}_j(\cdot)$ in a sufficiently rich probability space such that (i) $\widehat{Z}'_j(\cdot)=_d \widehat Z_j(\cdot)$ conditional on the data $(\mathbf{y}, \mathbf{X})$ and (ii) for some positive non-vanishing sequence $r_n$ and for all $\eta>0$, \[\mathbb{P}^*\left[\sup_{\mathbf{x}\in\mathcal{X}} |\widehat{Z}'_j(\mathbf{x}) - Z_j(\mathbf{x})| \geq \eta r_n^{-1} \right] = o_\mathbb{P}(1),\] where $\mathbb{P}^*[\cdot]=\mathbb{P}[\cdot|\mathbf{y}, \mathbf{X}]$ denotes the probability operator conditional on the data. When such a feasible process exists, we write $\widehat{Z}_j(\cdot) =_{d^*} Z_j(\cdot) + o_{\mathbb{P}^*}(r_n^{-1})$ in $\mathcal{L}^\infty(\mathcal{X})$. From a practical perspective, sampling from $\widehat{Z}_j(\cdot)$, conditional on the data, is possible and provides a valid distributional approximation.

Our first construction is a direct plug-in approach using the conclusion of Theorem (ref) and (ref). All unknown objects are replaced by consistent estimators already used in the feasible $t$-statistics:

The validity of this method is shown in Theorem 6.3 of the main paper. We restate it in the following whose proof is available in Section 8.5 of the main paper.

Our second construction, which is not reported in the main paper, is inspired by the intermediate approximation in Lemma 8.2 of the main paper. Again, $\boldsymbol{\gamma}_{\mathbf{q},j}(\mathbf{x})$ and $\Omega_j(\mathbf{x})$ can be simply replaced by their sample analogues, and $\zeta_i$'s can be generated by sampling from an $n$-dimensional standard Gaussian vector independent of the data. The key unknown quantity is the conditional variance function $\sigma^2(\mathbf{x})=\mathbb{E}[\varepsilon_i^2|\mathbf{x}_i=\mathbf{x}]$. If the errors are homoskedastic, one can simply let $\sigma^2(\cdot)=1$. If not, we need an estimator $\widehat{\sigma}(\cdot)$ of the conditional variance satisfying a mild condition on uniform convergence, and then we can construct \[\widehat{z}_j(\mathbf{x})=\frac{\widehat{\bm{\gamma}}_{\mathbf{q},j}(\mathbf{x})'}{\sqrt{\widehat{\Omega}_{j}(\mathbf{x})}} \mathbb{G}_n[\bm{\bm{\Pi}}_{j}(\mathbf{x}_i)\zeta_i\widehat{\sigma}(\mathbf{x}_i)],\quad \mathbf{x}\in\mathcal{X}, \quad j=0,1,2,3. \] In practice one may simply construct a nonparametric estimator of $\sigma^2(\cdot)$ by using, for example, regression splines or smoothing splines. See Eggermont-LaRiccia_2009_Book for more details. We do not elaborate on this issue, but the next theorem shows the validity of this approach.

Our third approach to approximating the infeasible $Z_j(\cdot)$ employs an easy-to-implement wild bootstrap procedure. Specifically, for i.i.d. bounded random variables $\{\omega_i:1\leq i\leq n\}$ with $\mathbb{E}[\omega_i]=0$ and $\mathbb{E}[\omega_i^2]=1$ independent of the data, we construct bootstrapped $t$-statistics:

where $\widehat{\varepsilon}_{i,j}$ is defined in (ref), and the bootstrap Studentization $\widehat{\Omega}^*_{j}(\mathbf{x})$ is constructed using $\widehat{\bm{\Sigma}}_{j}^* = \mathbb{E}_n[\bm{\bm{\Pi}}_{j}(\mathbf{x}_i) \bm{\bm{\Pi}}_{j}(\mathbf{x}_i)' \omega_i^2 \widehat{\varepsilon}_{i,j}^2]$.

In comparison with previous plug-in approximations, the bootstrap method requires an additional Gaussian approximation step as in either Theorem (ref) or Theorem (ref) (conditional on the data). These results are stated in the following two theorems. $\mathbb{E}^*[\cdot]$ denotes the expectation conditional on the data.

The additional condition in (ii) of the above theorem is generally difficult to verify in practice unless some restrictive conditions, e.g. bounded support of $\varepsilon_i$, are imposed.

Application: Confidence Bands

An immediate application of Theorems (ref)--(ref) is to construct confidence bands for the regression function or its derivatives. Specifically, for $j=0,1,2,3$, and $\alpha\in(0,1)$, we seek a quantile $q_j(\alpha)$ such that \[\mathbb{P}\left[\sup_{\mathbf{x}\in\mathcal{X}} |\widehat{T}_j(\mathbf{x})| \leq q_j(\alpha) \right] =1-\alpha + o(1), \] which then can be used to construct uniform $100(1-\alpha)$-percent confidence bands for $\partial^\mathbf{q} \mu(\mathbf{x})$ of the form \[\left[\widehat{\partial^\mathbf{q} \mu}_j(\mathbf{x}) \pm q_j(\alpha)\sqrt{\widehat{\Omega}_{j}(\mathbf{x})/n} \;:\; \mathbf{x}\in\mathcal{X} \right]. \] The following theorem establishes a valid distributional approximation for the suprema of the $t$-statistic processes $\{\widehat{T}_j(\mathbf{x}):\mathbf{x}\in\mathcal{X}\}$ using Chernozhukov-Chetverikov-Kato_2014b_AoS to convert our strong approximation results into convergence of distribution functions in terms of Kolmogorov distance.

Examples

We now illustrate how Assumptions (ref), (ref), and (ref) are verified for popular basis choices, including $B$-splines, wavelets, and piecewise polynomials, and show when (ref) holds. Note well that even in the cases where (ref) holds, we will not make use of this property neither herein nor in our software implementation Cattaneo-Farrell-Feng_2019_lspartition, since both bias terms of Lemma (ref) may be important in finite samples.

The first three subsections illustrate each of these assuming rectangular support, $\mathcal{X}=\otimes_{\ell=1}^d\mathcal{X}_\ell\subset\mathbb{R}^d$, and correspondingly use tensor-product partitions. We note that these three in fact overlap in the special cases of $m=1$ on a tensor-product partition (see Equation (ref)), leading to so-called zero-degree splines, Haar wavelets, and regressograms (piecewise constant). The final subsection considers a general partition. Recall that whenever $\Delta$ covers only strict subset of $\mathcal{X}$, all our results hold on that subset; the first three subsections may also be interpreted in this light.

B-Splines on Tensor-Product Partitions

A univariate spline is a piecewise polynomial satisfying certain smoothness constraints. For some integer $m_\ell\geq 1$, let $\mathcal{S}_{\Delta_\ell,m_\ell}$ be the set of splines of order $m_\ell$ with univariate partition $\Delta_\ell$, i.e., \[\small \mathcal{S}_{\Delta_\ell,m_\ell}=\Big\{s\in \mathcal{C}^{m_\ell-2}(\mathcal{X}_\ell)\colon s(x)\text{ is a polynomial of degree}\; (m_\ell-1)\; \text{on each } [t_{\ell,l_\ell},t_{\ell,l_\ell+1}]\Big\}.\] If $m_\ell=1$, it is interpreted as the set of piecewise constant functions that are discontinuous on $\Delta_\ell$. $\mathcal{S}_{\Delta_\ell, m_\ell}$ is a vector space and can be spanned by many equivalent representing bases. $B$-splines as a local basis are well studied in literature and enjoy many nice properties Micula-Micula_1999_Handbook,Schumaker_2007_Book.

Define an extended knot sequence $\Delta_{\ell,e}$ such that \[t_{\ell,-m_\ell+1}=t_{\ell,-m_\ell+2}=\cdots=t_{\ell,0}<t_{\ell,1}<\cdots<t_{\ell,\kappa_{\ell}-1}<t_{\ell,\kappa_{\ell}}=\cdots=t_{\ell,\kappa_{\ell}+m_\ell-1},\] Then, the $m_\ell$-th order $B$-spline with knot sequence $\Delta_{\ell,e}$ is \[p_{l_\ell,m_\ell}(x_\ell)=(t_{\ell,l_\ell}-t_{\ell,l_\ell-m_\ell})[t_{\ell,l_\ell-m_\ell}, \cdots, t_{\ell,l_\ell}]_t(t-x_\ell)_+^{m_\ell-1},\; l_\ell=1,\ldots,K_\ell=\kappa_{\ell}+m_\ell-1,\] where $(a)_+=\max\{a,0\}$ and $[t_1,t_2,\cdots,t_l]_t\,\mu(t,x)$ denotes the divided difference of $\mu(t,x)$ with respect to $t$, given a sequence of knots $t_1\leq t_2\leq\cdots\leq t_l$.

When there is no chance of confusion, we shall write $p_{l_\ell}(x_\ell)$ instead of $p_{l_\ell, m_\ell}(x_\ell)$. Accordingly, the space of tensor-product polynomial splines of order $\mathbf{m}=(m_1,\ldots,m_d)$ with partition $\Delta$ is spanned by tensor products of univariate $B$-splines \[\mathcal{S}_{\Delta,\mathbf{m}}=\otimes_{\ell=1}^d \mathcal{S}_{\Delta_\ell,m_\ell}=\text{span}\{p_{l_1}(x_1)p_{l_2}(x_2)\cdots p_{l_d}(x_d)\}_{l_1=1,\ldots,l_d=1}^{K_1,\ldots,K_d}.\] We have a total of $K=\prod_{\ell=1}^{d}K_\ell$ basis functions. The order of univariate basis could vary across dimensions, but for simplicity we assume that $m_1=\cdots=m_d=m$ and hence write $\mathcal{S}_{\Delta, m}=\mathcal{S}_{\Delta, \mathbf{m}}$. Also, let $p_{l_1\ldots l_d}(\mathbf{x})=p_{l_1}(x_1)\cdots p_{l_d}(x_d)$ to further simplify notation.

Arrange this set of basis functions by first increasing $l_d$ from $1$ to $K_d$ with other $l_\ell$'s fixed at $1$ and then increasing $l_\ell$ sequentially. As a result, we construct a one-to-one correspondence $\varphi$ mapping from $\{(l_1,\ldots,l_d): 1\leq l_\ell\leq K_\ell, \ell=1,\ldots,d\}$ to $\{1,\ldots, K\}$. Then, we can write $p_k(\mathbf{x})=p_{\varphi^{-1}(k)}(\mathbf{x})$, $k=1,\ldots, K$, which is consistent with our notation in Section (ref).

The following lemma shows that Assumptions (ref) and (ref) hold for $B$-splines.

Equation (ref) gives a precise definition of the strong quasi-uniform condition on the partition scheme. Assumption (ref) requires that the volume of all cells vanish at the same rate but allows for any constant proportionality between neighboring cells. Presently, cells are further restricted to be of the same volume asymptotically, and further, a specific rate is required that is related to the smoothness of $\mu(\cdot)$. Note that, for example, equally spaced knots satisfy this conditional trivially. For other schemes, additional work may be needed. Under (ref), Barrow-Smith_1978_QAM obtained an expression of the leading asymptotic error of univariate splines, which was later used by Zhou-Shen-Wolfe_1998_AoS and Zhou-Wolfe_2000_SS, among others. Lemma (ref) extends previous results to the multi-dimensional case, in addition to showing that the high-level conditions in Assumption (ref) and (ref) hold for $B$-Splines.

When (ref) fails, it may be possible to obtain results, but additional cumbersome notation would be needed and the results may be less useful. However, it is straightforward to verify that without (ref), there still exists $s^*\in\mathcal{S}_{\Delta,m}$ such that for all $[\boldsymbol{\varsigma}]\leq \varsigma$,

where recall that $\mathcal{S}_{\Delta,m}$ denotes the linear span of $\mathbf{p}(\mathbf{x};\Delta,m)$. However, this result does not allow for bias correction or IMSE expansion.

Finally, $\Lambda_m$ contains only the multi-indices corresponding to $m$-th order partial derivatives of $\mu(\mathbf{x})$. This is due to the fact that, as a variant of polynomial approximation, the total order of tensor-product $B$-splines is not fixed at $m$, i.e.\ some higher-order components are retained that approximate terms with $m$-th order cross partial derivatives. This feature distinguishes tensor-product splines from multivariate splines which do control the total order of bases.

Wavelets on Tensor-Product Partitions

Our results apply to compactly supported wavelets, such as Cohen-Daubechies-Vial wavelets Cohen_1993_ACHA. For more background details see Meyer_1995_Book,Hardle-Kerkyacharian-Picard-Tsybakov_2012_Book,Chui_2016_Book, and references therein.

To describe these estimators, we first introduce the general definition of wavelets and then show that a large class of orthogonal wavelets satisfy our high-level assumptions. For some integer $m\geq 1$, we call $\phi$ a (univariate) scaling function or “father wavelet” of degree $m-1$ if (i) $\int_\mathbb{R}\phi(x_\ell)\,dx_\ell=1$, (ii) $\phi$ and all its derivatives up to order $m-1$ decrease rapidly as $|x_\ell|\rightarrow \infty$, and (iii) $\{\phi(x_\ell-l):l\in\mathbb{Z}\}$ forms a Riesz basis for a closed subspace of $L_2(\mathbb{R})$. A real-valued function $\psi$ is called a (univariate) “mother wavelet” of degree $m-1$ if (i) $\int_{\mathbb{R}}x_\ell^v\psi(x_\ell)dx_\ell=0$ for $0\leq v\leq m-1$, (ii) $\psi$ and all its derivatives decrease rapidly as $|x_\ell|\rightarrow\infty$, and (iii) $\{2^{s/2}\psi(2^sx_\ell-l):s,l\in\mathbb{Z}\}$ forms a Riesz basis of $L_2(\mathbb{R})$.

We restrict our attention to orthogonal wavelets with compact support. The father wavelet $\phi$ and the mother wavelet $\psi$ are both compactly supported, and for any integer $s_0\geq 0$, any function in $L_2(\mathbb{R})$ admits the following $(m-1)$-regular wavelet multiresolutional expansion: \[\mu(x_\ell)=\sum_{l=-\infty}^{\infty}a_{s_0l}\phi_{s_0l}(x_\ell)+\sum_{s=s_0}^{\infty}\sum_{l=-\infty}^{\infty}b_{sl}\psi_{sl}(x_\ell),\;x_\ell\in\mathbb{R},\quad \text{where}\]

and $\{\phi_{s_0l},l\in\mathbb{Z};\psi_{sl},s\geq s_0,l\in\mathbb{Z}\}$ is an orthonormal basis of $L_2(\mathbb{R})$. Accordingly, to construct an orthonormal basis of $L_2([0,1])$, one can pick those basis functions supported in the interior of $[0,1]$ and add some boundary correction functions. For details of construction of such a basis see, for example, Cohen_1993_ACHA. With a slight abuse of notation, in what follows we write $\{\phi_{s_0l},l\in\mathcal{L}_{s_0};\psi_{sl},s\geq s_0,l\in\mathcal{L}_s\}$ for an orthonormal basis of $L_2([0,1])$ rather than $L_2(\mathbb{R})$ where $\mathcal{L}_{s_0}$ and $\mathcal{L}_s$ are some proper index sets depending on $s_0$ and $s$ respectively. Again, we use tensor product to form a multidimensional wavelet basis, and then it fits into the context of our analysis.

The compactness of the father wavelet is needed to verify Assumption (ref). To see this, consider a standard multiresolutional analysis setting. A large space, say $L_2([0,1])$, is decomposed into a sequence of nested subspaces $\{0\}\subset\cdots\subset V_{-1}\subset V_0\subset V_1\cdots\subset L_2([0,1])$. Generally, $\{\phi_{sl},l\in\mathcal{L}_s\}$ constitutes a basis for $V_s$, and $\{\psi_{sl},l\in\mathcal{L}_s\}$ forms a basis for the orthogonal complement $W_s$ of $V_s$ in $V_{s+1}$. Since the support of $\phi$ is compact and $\{\phi_{sl}\}$ is generated simply by dilation and translation of $\phi$, $[0,1]$ can be viewed as implicitly partitioned. Specifically, by increasing the resolution level from $s$ to $s+1$, the length of the support of $\phi_{sl}$ is halved. Hence, it is equivalent to placing additional partitioning knots at the midpoint of each subinterval. In addition, as the sample size $n$ grows, we allow the resolution level $s$ to increase (thus written as $s_n$ when needed), but each basis function $\phi_{s_nl}$ is supported by only a finite number of subintervals since the generating scaling function is compactly supported.

The above discussion associates the generating process of wavelets at different levels with a partitioning scheme. It is easy to see that Assumption (ref) is automatically satisfied in this case, and given a resolution level $s_n$, the mesh width $b=2^{-s_n}$. On the other hand, given a starting level $s_0$, $\{\phi_{s_0l},l\in\mathcal{L}_{s_0};\psi_{sl},s_0\leq s\leq s_n,l\in\mathcal{L}_s\}$ does not satisfy Assumption (ref), since some functions in such a basis have increasing supports as the resolution level $s_n\to\infty$. Nevertheless, they span the same space as $\{\phi_{s_n,l},l\in\mathcal{L}_{s_n}\}$ and hence there exists a linear transformation that connects the two equivalent bases. Least squares estimators are invariant to such a transformation. Therefore, we will focus on the following tensor-product (father) wavelet basis

where $\boldsymbol{\phi}_{s_n}$ is a vector containing all functions in $\{\phi_{s_n,l},l\in\mathcal{L}_{s_n}\}$. By multiplying the basis by $2^{-s_n/2}$, we drop the normalizing constants that originally appear in the construction of orthonormal basis. As the next lemma shows, this large class of orthogonal wavelet bases satisfy our assumptions.

In addition to verifying that our high-level assumptions hold for wavelets, this result gives a novel asymptotic error expansion for multidimensional compactly supported wavelets. Our derivation employs the ideas in Sweldens-Piessens_1994_NM and exploits the tensor-product structure of the wavelet basis. The end result is similar to tensor-product splines, in the sense that the total order of the approximating basis is not fixed at $m$ and thus $\Lambda_m$ is the same as that for $B$-splines.

Generalized Regressograms on Tensor-Product Partitions

To construct the generalized regressogram, we first define the piecewise polynomials. What will distinguish these from splines is that (i) each polynomial is supported on exactly one cell, and relatedly (ii) no continuity is assumed between cells. First, for some fixed integer $m\in\mathbb{Z}_+$, let $\mathbf{r}(x_\ell)=(1,x_\ell,\ldots,x_\ell^{m-1})'$ denote a vector of powers up to degree $m-1$. To extend it to a multidimensional basis, we take the tensor product of $\mathbf{r}(x_\ell)$, denoted by a column vector $\mathbf{R}(\mathbf{x})$. The total order of such a basis is not fixed, and its behavior is more similar to tensor-product $B$-splines. Following Cattaneo-Farrell_2013_JoE, we exclude all terms with degree greater than $m-1$ in $\mathbf{R}(\mathbf{x})$. Hence the remaining elements in $\mathbf{R}(\mathbf{x})$ are given by $\mathbf{x}^{\boldsymbol{\alpha}}=x_1^{\alpha_1}\cdots x_d^{\alpha_d}$ for a unique $d$-tuple $\boldsymbol{\alpha}$ such that $[\boldsymbol{\alpha}]\leq m-1$. Then we “localize” this basis by restricting it to a particular subrectangle $\delta_{l_1\cdots l_d}$. Specifically, we write $\mathbf{p}_{l_1\ldots l_d}(\mathbf{x})=\mathds{1}_{\delta_{l_1\ldots l_d}}(\mathbf{x})\mathbf{R}(\mathbf{x})$, where $\mathds{1}_{\delta_{l_1\ldots l_d}}(\mathbf{x})$ is equal to 1 if $\mathbf{x}\in \delta_{l_1\ldots l_d}$ and $0$ otherwise. Finally, we rotate the basis by centering each basis function at the start point of the corresponding cell and scale it by interval lengths. We can arrange the basis functions according to a particular ordering $\varphi$: first order the subrectangles the same way as that for $B$-splines, and then within each subrectangle $\delta_{l_1\ldots l_d}$ the basis functions in $\mathbf{R}(\mathbf{x})$ are ordered ascendingly in $\boldsymbol{\alpha}$ and $\ell=1,\ldots, d$. Using the same notation as in Section (ref), we still write the basis as $\mathbf{p}(\mathbf{x})=(p_1(\mathbf{x}),\cdots, p_K(\mathbf{x}))'$.

The following lemma shows that Assumptions (ref) and (ref) hold for generalized regressograms.

The leading asymptotic error obtained in Lemma (ref) differs from the one in Cattaneo-Farrell_2013_JoE because it is expressed in terms of orthogonal polynomials. Here we employ Legendre polynomials $\bar{P}_m(x)$ orthogonal on $[-1,1]$ with respect to the Lebesgue measure, and then shift them to $P_m(x)=\bar{P}_m(2x-1)$. Thus the shifted Legendre polynomials are orthogonal on $[0,1]$. Further, we allow for more general partitioning schemes.

Generalized Regressograms on General Partitions

Moving away from tensor-product partitions will impede verification of Assumption (ref) in general. A more typical result, such as Equation (ref) can be more easily proven in generality. In practice, however, many data-driven partitioning selection procedures, such as those induced by regression trees, do not necessarily lead to a tensor-product partition. There is also a large literature in approximation theory discussing bases constructed on triangulations or other general partition schemes Micula-Micula_1999_Handbook. To demonstrate the power of our theory, we show presently that for the generalized regressogram we can obtain concrete results on general partitions.

The basis is as described above, but now allowing for cells of general shape. Let $\mathcal{X}$ be convex and $\mathbf{t}_\delta$ be an arbitrary point in each $\delta \in \Delta$. Hence, $\mathbf{t}_{\delta_\mathbf{x}}$ is simply a point in the cell containing $\mathbf{x}$. Construct the rotated piecewise polynomial basis as in Section (ref), but centering the basis at $\mathbf{t}_\delta$ and rescaling it by the diameter of $\delta$.

It is a challenging task to construct orthogonal polynomial bases on non-rectangular domains, which makes Equation (ref) hard to be satisfied when employing partitioning estimators on general partitions. Thus, our more general characterization (and correction) of the bias is quite useful in this case.

Discussion and Extensions

Connecting Splines and Piecewise Polynomials

There is an important connection between splines and piecewise polynomial bases. Essentially, the former can be viewed as a piecewise polynomial basis with certain continuity restrictions. Therefore, the estimators based on the two bases are linked by utilizing the well-known results about regressions with linear constraints. See Buse_1977_JASA for an illustration of cubic splines as a special case of restricted least squares.

Formally, let us start with the (rotated) piecewise polynomials $\mathbf{p}$ discussed in Section (ref). For expositional simplicity, we only discuss tensor-product partitions here. Clearly, $\mathbf{p}$ spans a vector space $\mathcal{P}$ containing all “piecewise” polynomials with degree no greater than $m-1$ on $\mathcal{X}$: \[\mathcal{P}:=\Big\{s(\cdot):\;s(\cdot)=\sum_{k=1}^{K}a_k p_k(\cdot), \; a_k\in\mathbb{R}\Big\}.\] Functions in this space are continuous within each subrectangle, but might have “jumps” along boundaries of cells. Then by imposing certain continuity restrictions on functions in $\mathcal{P}$, we can construct a subspace $\mathcal{S}\subset\mathcal{P}$ \[\mathcal{S}:=\Big\{s(\cdot):\; s(\cdot) = \sum_{k=1}^{K}a_k p_k(\cdot), \; a_k\in\mathbb{R}, \text{ and}\; \partial^{\boldsymbol{\varsigma}}s(\cdot)\in \mathcal{C}(\mathcal{X}),\; \forall [\boldsymbol{\varsigma}] \leq \iota\Big\}.\] where $\iota\leq m-2$ is a positive integer that controls the smoothness of the basis (and thus the smoothness of the estimated function). Since the derivatives of polynomials are linear in coefficients, the continuity constraints are linear.

Now consider the following restricted least squares:

where $\mathbf{P}=(\mathbf{p}(\mathbf{x}_1),\ldots, \mathbf{p}(\mathbf{x}_n))'$ and $\mathbf{R}$ is a $\vartheta_\iota\times K$ restriction matrix. $\vartheta_\iota$ denotes the number of restrictions depending on the required smoothness $\iota$. If there are no redundant constraints, $\mathbf{R}$ has full row rank. It is well known that given the restriction matrix $\mathbf{R}$, the least squares estimator can be written as \[\widehat{\boldsymbol{\beta}}_{\mathtt{r}}=[\mathbf{I}-(\mathbf{P}'\mathbf{P})^{-1}\mathbf{R}'(\mathbf{R}(\mathbf{P}'\mathbf{P})^{-1}\mathbf{R}')^{-1}\mathbf{R}](\mathbf{P}'\mathbf{P})^{-1}\mathbf{P}' \mathbf{y}.\] Since the unrestricted least squares estimator $\widehat{\boldsymbol{\beta}}_{\mathtt{ur}}=(\mathbf{P}'\mathbf{P})^{-1}\mathbf{P}'\mathbf{y}$, the above equation implies that $\widehat{\boldsymbol{\beta}}_{\mathtt{r}}=[\mathbf{I}-(\mathbf{P}'\mathbf{P})^{-1}\mathbf{R}'(\mathbf{R}(\mathbf{P}'\mathbf{P})^{-1}\mathbf{R}')^{-1}\mathbf{R}]\widehat{\boldsymbol{\beta}}_{\mathtt{ur}} =:(\mathbf{I}-\mathbf{U})\widehat{\boldsymbol{\beta}}_{\mathtt{ur}}$. Therefore,

where $\widehat{\mu}_{\mathtt{ur}}(\mathbf{x}):=\mathbf{p}(\mathbf{x})'\widehat{\boldsymbol{\beta}}_{\mathtt{ur}}$ is the unrestricted estimator. With this relation we can derive expressions of bias and variance for the restricted estimator. Clearly, the conditional variance $\mathbb{V}[\widehat{\mu}_{\mathtt{r}}(\mathbf{x})|\mathbf{X}]=\mathbf{p}(\mathbf{x})'(\mathbf{I}-\mathbf{U}) \mathbb{V}[\widehat{\mu}_{\mathtt{ur}}(\mathbf{x})|\mathbf{X}](\mathbf{I}-\mathbf{U})'\mathbf{p}(\mathbf{x})$. On the other hand, as shown in Lemma (ref), there exists $s^*(\mathbf{x}):=\mathbf{p}(\mathbf{x})'\boldsymbol{\beta}^*$ such that $\|\mu-s^*+\mathscr{B}_{m,\mathbf{0}}\|_{L_\infty(\mathcal{X})} =o(h^m)$. Hence

where $\pmb{\mathscr{B}}_{m, \mathbf{0}}=(\mathscr{B}_{m, \mathbf{0}}(\mathbf{x}_1), \ldots, \mathscr{B}_{m, \mathbf{0}}(\mathbf{x}_n))'$. $\mathbf{p}(\mathbf{x})'\mathbf{U}\boldsymbol{\beta}^*$ can be viewed as a measure of to what extent the continuity constraints are satisfied. When $s^*(\mathbf{x})$ is continuous up to order $\iota$, $\mathbf{R}\boldsymbol{\beta}^*$ is exactly $0$ and thus this term vanishes.

To repeat the previous analysis for such a restricted estimator, the main challenge is to analyze the asymptotic properties of the outer product of the restriction matrix. Once we know its limiting eigenvalue distributions (e.g., bounds on extreme eigenvalues), $\mathbf{U}$ can be properly bounded, and then the conclusions in previous sections may be established with similar proofs. Unfortunately, for general multidimensional cases, it is difficult and tedious to specify a non-redundant set of continuity constraints and analyze the eigenvalue distributions of $\mathbf{R}\mathbf{R}'$, thus damping the usefulness of such a method, whereas when $d=1$, the restriction matrix is quite straightforward and well bounded.

Formally, let $\mathbf{R}$ be a restriction matrix corresponding to $\bar{\iota}=\iota+1$ continuity constraints at each partitioning knot, implying that $s\in\mathcal{S}$ has $\iota$ continuous derivatives. We write the $\ell$th restriction at the $k$th knot as $\mathbf{r}_{k\ell}$, corresponding to the requirement that the $(\ell-1)$th derivatives of functions in $\mathcal{S}$ are continuous. Explicitly, the entire restriction matrix admits the following structure:

As $\kappa$ increases, the dimension of $\mathbf{R}$ also grows. Moreover, $\iota$ cannot exceed $m-2$ since when $m$ continuity constraints are imposed, $\mathcal{S}$ degenerates to the space of global polynomials of degree no greater than $m-1$.

The asymptotic behavior of the restricted estimator is closely related to the extreme eigenvalues of $\mathbf{R}\mathbf{R}'$. For a general restriction matrix $\mathbf{R}$ with fixed dimensions, if it contains non-redundant constraints only, the minimum eigenvalue of $\mathbf{R}\mathbf{R}'$ is nonzero. When the number of constraints $\vartheta_\iota\rightarrow\infty$, however, the limit of the minimum eigenvalue does not have to be nonzero, and its limiting behavior depends on the specific structure of constraints. The next lemma shows that for the particular restricted estimators considered here, the eigenvalues of $\mathbf{R}\mathbf{R}'$ are indeed bounded and bounded away from zero uniformly over the number of knots for $\iota\leq m-2$.

The proof of this lemma employs the specific structure of the restriction matrix. Generally, the outer product of $\mathbf{R}$ takes the following form:

where

Importantly, the form described in (ref) is usually referred to as a (tridiagonal) block Toeplitz matrix, meaning that it is a tridiagonal block matrix containing blocks repeated down the diagonals. It is well known that its asymptotic eigenvalue distribution is characterized by the Fourier transform $$\pmb{\mathscr{T}}_{\bar{\iota}}(\omega)=\mathbf{A}+(\mathbf{B}+\mathbf{B}')\cos\omega,\quad \omega\in[0, 2\pi].$$ As $\kappa\rightarrow\infty$, $\lambda_{\min}(\mathbf{R}\mathbf{R}')$ converges to the minimum attained by the minimum eigenvalue of $\pmb{\mathscr{T}}_{\bar{\iota}}(\omega)$ as a function of $\omega$ on $[0,2\pi]$. Similarly, the limit of $\lambda_{\max}(\mathbf{R}\mathbf{R}')$ is the maximum attained by the maximum eigenvalue of $\pmb{\mathscr{T}}_{\bar{\iota}}(\omega)$.

Comparison of Bias Correction Approaches

We make some comparison of the three bias correction approaches considered in this paper.

First, the higher-order correction ($j=1$) and least squares correction ($j=2$) are closely related. Let $\mathbf{P}=(\mathbf{p}(\mathbf{x}_1), \ldots, \mathbf{p}(\mathbf{x}_n))'$ and $\tilde{\mathbf{P}}=(\bm{\tilde{\mathbf{p}}}(\mathbf{x}_1), \ldots, \bm{\tilde{\mathbf{p}}}(\mathbf{x}_n))'$. Clearly,

where $\mathbf{M}_{\bm{\tilde{\mathbf{p}}}}:=\mathbf{I}-\tilde{\mathbf{P}}(\tilde{\mathbf{P}}'\tilde{\mathbf{P}})^{-1} \tilde{\mathbf{P}}'$. Importantly, when $\mathbf{p}$ and $\bm{\tilde{\mathbf{p}}}$ generate nested models, i.e., there exists a transformation matrix $\boldsymbol{\Upsilon}$ such that $\mathbf{p}(\cdot)=\boldsymbol{\Upsilon}\bm{\tilde{\mathbf{p}}}(\cdot)$, it is easy to see that $\mathbf{M}_{\bm{\tilde{\mathbf{p}}}}\mathbf{P}=\mathbf{0}$. Thus the higher-order and least squares bias correction approaches are equivalent. When $\bm{\tilde{\mathbf{p}}}$ and $\mathbf{p}$ are not nested bases, the two methods will typically differ in variance and bias.

To compare their variance, we generally have

When $\sigma^2(\mathbf{x})=\sigma^2$, the covariance term is $0$, and thus $\widehat{\mu}_2$ has variance no less than that of $\widehat{\mu}_1$. The same conclusion is true for asymptotic variance as shown in the proof of Lemma (ref).

Regarding their bias, let $\mathfrak{B}_{\tilde{m}, \mathbf{0}}(\mathbf{x}) =\mathbb{E}[\widehat{\mu}_1(\mathbf{x})|\mathbf{X}] - \mu(\mathbf{x})$ denote the conditional bias of $\widehat{\mu}_1(\mathbf{x})$, and $\boldsymbol{\mathfrak{B}}_{\tilde{m},\mathbf{0}}:= (\mathfrak{B}_{\tilde{m},\mathbf{0}}(\mathbf{x}_1),\ldots,\mathfrak{B}_{\tilde{m},\mathbf{0}}(\mathbf{x}_n))'$. Then $ \mathbb{E}[\widehat{\mu}_2(\mathbf{x})|\mathbf{X}]-\mu(\mathbf{x}) =\mathfrak{B}_{\tilde{m},\mathbf{0}}(\mathbf{x})- \mathbf{p}(\mathbf{x})'(\mathbf{P}'\mathbf{P})^{-1}\mathbf{P}'\boldsymbol{\mathfrak{B}}_{\tilde{m},\mathbf{0}} $. Clearly, the second term will asymptotically get close to the projection of $\mathfrak{B}_{\tilde{m},\mathbf{0}}(\mathbf{x})$ onto the space spanned by $\mathbf{p}$: $$\mathscr{L}_{\mathbf{p}}[\mathfrak{B}_{\tilde{m},\mathbf{0}}](\mathbf{x}) :=\mathbf{p}(\mathbf{x})'(\mathbb{E}[\mathbf{p}(\mathbf{x})\mathbf{p}(\mathbf{x})'])^{-1}\mathbb{E}[\mathbf{p}(\mathbf{x})\mathfrak{B}_{\tilde{m},\mathbf{0}}(\mathbf{x})]$$ where $\mathscr{L}_{\mathbf{p}}[\cdot]$ denotes the projection operator. Therefore, when $\mathbf{p}$ and $\bm{\tilde{\mathbf{p}}}$ are not nested bases, we typically have the bias of $\widehat{\mu}_2$ no greater than that of $\widehat{\mu}_1$ in terms of $\|\cdot\|_{F,L_2(\mathcal{X})}$, where for a real-valued function $g(\cdot)$ on $\mathcal{X}$, $\|g\|_{F,L_2(\mathcal{X})}=(\int_{\mathcal{X}}|g(\mathbf{x})|^2dF(\mathbf{x}))^{1/2}$.

According to the discussion above, the higher-order and least squares bias correction approaches do not dominate each other in general, and whether one is preferred to the other depends on the data generating process and the relation between the two approximation spaces (or more precisely, the approximation power of $\bm{\tilde{\mathbf{p}}}$ for functions in the linear span of $\mathbf{p}$). Then a natural question follows: is there an optimal weighting scheme when $\mathbf{p}$ and $\bm{\tilde{\mathbf{p}}}$ are not nested? Again, assume $\varepsilon_i$'s are homoskedastic for simplicity, and we take a weighted average of $\widehat{\mu}_1$ and $\widehat{\mu}_2$: $\widehat{\mu}_{w,\mathtt{bc}}:=w\widehat{\mu}_2+(1-w)\widehat{\mu}_1$ where $w\in[0,1]$. Then using a conclusion in the proof of Lemma (ref), the change in the integrated asymptotic variance (weighted by the design density $f(\mathbf{x})$) is \[w^2\sigma^2\int_{\mathcal{X}}\mathbf{p}(\mathbf{x})' \mathbf{Q}_m^{-1}(\mathbf{Q}_m-\mathbf{Q}_{m,\tilde{m}}\mathbf{Q}_{\tilde{m}}^{-1}\mathbf{Q}_{m, \tilde{m}}')\mathbf{Q}_m^{-1}\mathbf{p}(\mathbf{x})f(\mathbf{x})d\mathbf{x}=: w^2\bar{\mathscr{V}}. \] On the other hand, by the property of projection operator, the change in the integrated squared bias (weighted by the design density $f(\mathbf{x})$) is \[(w^2-2w)\int_{\mathcal{X}}\left(\mathscr{L}_{\mathbf{p}}[\mathfrak{B}_{\tilde{m},\mathbf{0}}](\mathbf{x})\right)^2f(\mathbf{x})d\mathbf{x}=:(w^2-2w)\bar{\mathscr{B}}. \] It is easy to see the optimal weight is $w^*=\bar{\mathscr{B}}/(\bar{\mathscr{B}}+\bar{\mathscr{V}})$. Clearly, when variance is very small (e.g., $\sigma^2$ is small), $w^*$ is close to $1$ and $\widehat{\mu}_2$ is preferred, whereas when bias is small, $w^*$ is close to $0$ and one may want to use $\widehat{\mu}_1$.

Next, the comparison of plug-in bias correction with the other two is more complicated since $\widehat{\mu}_3$ generally cannot be viewed as a regression estimator with additional covariates, and the covariance between $\widehat{\mu}_0$ and the estimated bias does not vanish. For piecewise polynomials, however, all three bias correction approaches are simply equivalent under certain conditions.

To see this, suppose $\mathbf{p}$ and $\bm{\tilde{\mathbf{p}}}$ are constructed on the same partitioning scheme $\Delta$, but the order of basis increases from $m$ to $m+1$. $\widehat{\mu}_0$ and $\widehat{\mu}_1$ are linked by (ref) since $\widehat{\mu}_0$ can be viewed as a restricted regression estimator compared with $\widehat{\mu}_1$. Specifically, one can construct a polynomial series of order $m+1$ (with degree no greater than $m$) within each cell $\delta\in\Delta$, and then implement a local regression restricting the coefficients of the polynomial terms of degree $m$ to be $0$. Then the restriction matrix in this case can be written as $\mathbf{R}=\left[\mathbf{0} \quad \mathbf{I}_\vartheta \right]$ where $\vartheta$ denotes the number of polynomial terms of degree $m$ and $\mathbf{R}$ is a $\vartheta\times\tilde{K}$ matrix. Plug it in (ref), and then use the formula for matrix inverse in block form to obtain $(\tilde{\mathbf{P}}'\tilde{\mathbf{P}})^{-1}$. It is easy to see that the second term on the RHS of (ref), $\mathbf{p}(\mathbf{x})'\mathbf{M}\hat{\boldsymbol{\beta}}_{\mathtt{ur}}$, is exactly the same as the leading bias derived in Cattaneo-Farrell_2013_JoE (see their proof of Theorem 3) with the $m$th derivative estimated by piecewise polynomials of order $m+1$. As explained in the proof of Lemma (ref), the leading approximation error can be alternatively expressed in terms of Legendre polynomials. Asymptotically, the two expressions are equivalent since the “locally” orthogonalized polynomials of degree $m$ will converge to the $m$th Legendre polynomials given by Lemma (ref).

For splines or wavelets, we do not have the above equivalence in general, since bases of different orders do not generate nested spaces, and the relative performance of the three approaches depends on the relation between these approximation spaces.

Implementation Details

In this section, we briefly discuss implementation details about choosing the IMSE-optimal tuning parameters. We restrict our attention to tensor-product partitions with the same number of knots used in every dimension. Thus the tuning parameter reduces to a scalar $\kappa$ which denotes the number of subintervals used in every dimension. Also, we let the weighting function $w(\mathbf{x})$ be the density of $\mathbf{x}_i$. We offer two approaches: rule-of-thumb (ROT) and direct plug-in (DPI).

Rule-of-Thumb Choice

The rule-of-thumb choice is based on the special case considered in Theorem 4.2 of the main paper. Specifically, assume $\mathbf{q}=\mathbf{0}$ and knots are evenly spaced. The implementation steps are summarized as follows.

Rather than assume a uniform design and evenly-spaced partition, one may use a trimmed-from-below Gaussian reference model to estimate the density of $\mathbf{x}_i$ and adjust the variance and bias constants based on Theorem 4.2. Clearly, this choice of $\kappa$ is derived based on simplifying assumptions, but it still has the correct rate ($\asymp n^{\frac{1}{2m+d}}$) even in other cases.

Direct Plug-in Choice

A direct plug-in (DPI) procedure is summarized in the following.

Simulations

In this section, we present detailed simulation results. We consider the following regression functions:

where $\operatorname*{sign}(x)=-1,0$, or $1$ if $x<0$, $x=0$, or $x>0$ respectively, $\phi(\cdot)$ is the standard Normal density function, and $\tau(x)=(x-0.5)+8(x-0.5)^2+6(x-0.5)^3-30(x-0.5)^4-30(x-0.5)^5$.

Given each regression function, we generate $(y_i, \mathbf{x}_i)_{i=1,\ldots, n}$ by $y_i=\mu(\mathbf{x}_i) + \varepsilon_i$ where $\mathbf{x}_i\sim i.i.d.\,\mathsf{U}[0,1]^d$ and $\varepsilon_i\sim i.i.d. \,\mathsf{N}(0,1)$. For each model, we generate $5,000$ simulated datasets with $n=1,000$.

We first present results for (tensor-product) spline regressions in Tables (ref)--(ref). For each model, we use linear splines ($m=2$) to form the classical point estimates. For robust inference, we use quadratic splines ($\tilde{m}=3$) to implement bias correction. Both evenly-spaced and quantile-spaced knot placements are considered, and for simplicity point estimators and bias correction employ the same knot placements ($\Delta=\tilde{\Delta}$).

For each model, we present three sets of simulation evidence. First, for three fixed points, we calculate the (simulated) root mean squared error (RMSE), coverage rate (CR) and average confidence interval length (AL). The nominal coverage is set to be $95\%$. Second, to evaluate the performance of our rule-of-thumb (ROT) and direct plug-in (DPI) knot selection procedures, we show some basic summary statistics (mean, median, standard deviation, etc.) of the selected number of knots ($\hat{\kappa}_{\mathtt{ROT}}$ and $\hat{\kappa}_{\mathtt{DPI}}$). Third, for uniform confidence bands, we calculate: the proportion of values covered with probability at least $95\%$ (CP), average coverage errors (ACE), and average width of confidence band (AW), and uniform coverage rate (UCR). The quantile estimation based on the plug-in and bootstrap methods uses $1000$ random draws conditional on the data for each simulated dataset.

Next, we evaluate the performance of (Daubechies) wavelet regressions. Several distinctive features of wavelets should be noted. First, the tuning parameter of wavelets ($s$) is referred to as “resolution level", which relates to the number of knots $\kappa$ by $\kappa=2^s$. Thus, as $s$ increases, the number of series terms grows very fast. This issue is even more severe for tensor-product wavelets. Given our relatively small sample size, we restrict the wavelet-based simulations to Model 1-5. Second, due to the lack of smoothness of low-order wavelet basis, the plug-in bias correction approach ($j=3$) is not feasible unless very high-order wavelet basis is used to estimate the derivatives in the leading bias, which may not be of practical interest. Thus, we only report results based on higher-order and least-squares bias correction (as well as classical estimators).

Results for wavelets are reported in Table (ref)--(ref), where Daubechies (father) wavelets of order $2$ ($m=2$) are used to form classical estimators, and Daubechies wavelets of order $3$ ($\tilde{m}=3$) for bias correction. In addition, as explained above, the IMSE-optimal resolution level is related to the number of knots by $s_{\mathtt{IMSE}}=\log_2\kappa_{\mathtt{IMSE}}$. Similarly, for the estimated resolution levels, $\hat{s}_{\mathtt{ROT}}=\log_2\hat\kappa_{\mathtt{ROT}}$, and $\hat{s}_{\mathtt{DPI}}=\log_2\hat\kappa_{\mathtt{DPI}}$.

Our simulation results show that the bias correction methods perform generally well in both pointwise and uniform inference. The following are some practical guidance and caveats. First, as discussed in Section (ref), no bias correction approach among the three dominates the others in general. When high-order bias is large, their differences may be more pronounced. Second, our simulation shows that when $d>1$, our estimators still perform relatively well if a small number of knots are used. If not, as expected, their performance could be poor since the number of regressors explodes in such cases. Third, the ROT selection procedure used in the simulation study employs a global polynomial regression of degree $m+2$ to form the estimates of derivatives in the leading bias, which may be too conservative for highly nonlinear models. In such cases, a higher-order global polynomial regression may offer a better initial choice of the tuning parameter and improve the performance of DPI selection.

In the end, we present a figure (see Figure (ref)) as a visual illustration, which compares confidence bands based on different estimators suggested in our paper. We use Model 1 to generate a simulated dataset, and both plug-in approximation and wild bootstrap are used in constructing bands.

Proofs

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Theorem 4.1

In this section, we provide the proof of Theorem 4.1 in the main paper.

Proof of Theorem (ref)

Proof of Theorem 4.2

Proof of Lemma (ref)

Proof of Theorem 5.1

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 Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

\newcounter{model} \setcounter{model}{1}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\forloop[1]{model}{\value{model}}{\value{model} < 8}{

}}}}}}}}}}}

\setcounter{model}{1}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\forloop[1]{model}{\value{model}}{\value{model} < 8}{

}}}}}}}}}}}

\setcounter{model}{1}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 8}{

\addtocounter{model}{1}\forloop[1]{model}{\value{model}}{\value{model} < 8}{

}}}}}}}}}}}

\setcounter{model}{1}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\forloop[1]{model}{\value{model}}{\value{model} < 6}{

}}}}}}}}}}}

\setcounter{model}{1}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\setcounter{model}{\value{model}}\ifthenelse{\value{model} < 6}{

\addtocounter{model}{1}\forloop[1]{model}{\value{model}}{\value{model} < 6}{

}}}}}}}}}}}