A Localized Neural Network with Dependent Data:
Estimation and Inference

Jiti Gao†, Bin Peng† and Yanrong Yang∗

†Monash University and ∗The Australian National University

In this paper, we propose a localized neural network (LNN) model and then
develop the LNN based estimation and inferential procedures for dependent data
in both cases with quantitative/qualitative outcomes. We explore the use of
identification restrictions from a nonparametric regression perspective, and
establish an estimation theory for the LNN setting under a set of mild
conditions. The asymptotic distributions are derived accordingly, and we show
that LNN automatically eliminates the dependence of data when calculating the
asymptotic variances. The finding is important, as one can easily use
different types of wild bootstrap methods to obtain valid inference
practically. In particular, for quantitative outcomes, the proposed LNN
approach yields closed-form expressions for the estimates of some key
estimators of interest. Last but not least, we examine our theoretical
findings through extensive numerical studies.

Keywords: Activation Function; Binary Structure; Deep Learning; Kernel
Regression; Time Series

JEL classification: C14, C22, C45

##  1 Introduction

Neural network (NN) architecture has received increasing attention over the
last several decades. On relevant topics, a large number of papers have been
published in different journals such as Econometrica, The Annals of
Statistics, Journal of Machine Learning Research, Neural Networks,
Neurocomputing, etc. by experts from different disciplines. Apparently, we
cannot exhaust the literature, but refer the interested reader to Bartlett et
al. (2021) and Fan et al. (2021) for extensive reviews from a methodological
point of view.

NN usually includes three ingredients: input layer, hidden layer(s), and
output layer. We now briefly comment on them one by one. The input layer is
possibly the easiest one to understand, as it includes regressors only.
Because NN usually concerns about a generic unknown function of multiple
regressors of interest, it is unnecessary to involve any interaction terms
among these regressors in the input layer. The hidden layer(s) include the
parameters which need to be determined by training data. Once the parameters
are estimated, one can load the test dataset in order to evaluate the
performance of the newly constructed NN. Finally, the output layer receives
the estimation/decision/prediction which varies under different context.
Usually, there are two types of data for the output layer: quantitative
outcome (e.g., $y$ is continuous), and qualitative outcome (e.g., $y$ is
binary).

That said, in this paper we account for the features of input and output
layers, and propose to consider two specific models as follows:

$\displaystyle y=g(\mathbf{x})+\varepsilon\quad\text{for continuous outcome},$
(1.1) $\displaystyle
y=\bigg{\\{}\begin{array}[]{ll}1&g(\mathbf{x})-\varepsilon\geq 0\\\
0&\text{otherwise}\end{array}\quad\text{for binary outcome},$ (1.4)

in which $\mathbf{x}=(x_{1},\ldots,x_{d})^{\top}$ includes $d$ regressors with
$d$ being finite, and $\varepsilon$ is the idiosyncratic error component.
Further, we suppose that $g\,:\,[-a,a]^{d}\to\mathbb{R}$, where $a$ is
fixed111We work with fixed $a$ in the main text in order to be consistent with
the majority of the existing NN literature, and explain how to allow
$g(\cdot)$ to be defined on $\mathbb{R}^{d}$ in Appendix A.3.. The main goal
of this paper is to infer $g(\cdot)$ through a neural network architecture for
both models. In addition, for Model (1.4) only, we suppose that the respective
probability density function (PDF) and the cumulative distribution function
(CDF) of $\varepsilon$ are known, and denote the PDF and CDF by
$\phi_{\varepsilon}(\cdot)$ and $\Phi_{\varepsilon}(\cdot)$ respectively.
Here, the information about $\phi_{\varepsilon}(\cdot)$ and
$\Phi_{\varepsilon}(\cdot)$ is necessary for carrying on likelihood
estimation. Possibly, one may further merge both models into one framework
(say, $y=g(\mathbf{x},\varepsilon)$), and consider a generic loss function,
such as Chen and White (1999) and Section 1 of Bartlett et al. (2021), which
however will involve some high level conditions to facilitate the development.
As we will derive some results which are for Model (1.4) in particular and
will also show that Model (1.1) may enjoy some closed-form estimation
procedures, we do not further pursue a more general setup in this paper, but
would like to investigate the general setup in future research.

Until very recently, the investigation on NN architecture mainly focuses on
some fixed design (e.g., Cybenko, 1989 and many follow-up studies since then),
or uses independent and identically distributed (i.i.d.) data (e.g., Kohler
and Krzyzák (2017), Bauer and Kohler, 2019, Schmidt-Hieber, 2020, and many
references therein). There are only limited studies available for us to
understand NN with dependent data from theoretical perspective (see, Chen and
Shen, 1998, Chen, 2007, for example), although NN based methods have been
widely used to study different types of time series data in practice (e.g.,
Hill et al., 1996, Chen et al., 2001, Gu et al., 2021, Gu et al., 2020, just
to name a few). We would like to contribute along this line of research, so
throughout we assume the following time series data are observable:

$\displaystyle\\{(y_{t},\mathbf{x}_{t})\,|\,t\in[T]\\},$ (1.5)

where, for a positive integer $T$, $[T]$ stands for $\\{1,\ldots,T\\}$. The
values of $\\{y_{t}\\}$ vary with respect to the model being studied. We shall
not mention it again unless misunderstanding may arise.

Up to this point, it is worth further commenting on the hidden layer, which
usually involves lots of neurons. Here, a neuron stands for an activation
function, which maps a linear combination of the regressors from the input
layer to a certain range of the real line. As a consequence, lots of
parameters are actually involved in the hidden layer. The literature seems to
agree that in order to ensure the number of effective parameters within a
reasonable range, sparsity should be imposed (e.g., Schmidt-Hieber, 2020; Wang
and Lin, 2021; Fan and Gu, 2022; Zhong et al., 2022, just to name a few), and
in particular Schmidt-Hieber (2020) writes “the network sparsity assumes that
there are only few non-zero/active network parameters”, which essentially is
reflected by the number of active neurons. In this paper we consider Sigmoidal
function, $\sigma(\cdot):\mathbb{R}\to[0,1]$, as the activation function, and
acknowledge the growing literature on the non-sigmoidal rectified linear unit
(ReLU) function. We conjecture using ReLU function one may derive results
similar to what we are about to present below, but it requires some careful
investigation in future research.

All things considered, we draw Figure 1 for the purpose of illustration, which
presents a neural network architecture with only a small number of active
neurons. Specifically, the elements of $\mathbf{x}$ form the input layer,
while, in the hidden layer, only neurons in the dark area are activated.
Eventually, a final result is passed on to the output layer. To understand
Figure 1 fully, a few questions arise naturally:

  1. 1.

Provided a set of dependent time series data, how do we train a neural
network, such as that in Figure 1 ?

  2. 2.

Why are there only a small of number of neurons getting activated ?

  3. 3.

Given a set of activated neurons, can we further reduce the number of
parameters ? In other words, how does sparsity come to play just among the
activated neurons ?

  4. 4.

When the least absolute shrinkage and selection operator (i.e., LASSO) is
employed, how do we define the set of true parameters ?

  5. 5.

Can any inference (such as a confidence interval) be established ? and so
forth.

$x_{1}$$\vdots$$x_{d}$$\begin{array}[]{c}\sigma:\mathbf{x}\to[0,1]\\\
\vdots\\\
\sigma:\mathbf{x}\to[0,1]\end{array}$$\begin{array}[]{c}\sigma:\mathbf{x}\to[0,1]\\\
\vdots\\\
\sigma:\mathbf{x}\to[0,1]\end{array}$$\begin{array}[]{c}\sigma:\mathbf{x}\to[0,1]\\\
\vdots\\\ \sigma:\mathbf{x}\to[0,1]\end{array}$$\begin{array}[]{ll}\text{Model
\ref{Eq1}}&\widehat{g}(\mathbf{x})\\\ \text{Model
\ref{Eq2}}&\Phi_{\varepsilon}(\widehat{g}(\mathbf{x}))\end{array}$ Hidden
layer Input layer Output layer Figure 1: A Neural Network with Sparsity

In Figure 1, $\widehat{g}(\cdot)$ is an estimated version of $g(\cdot)$.

With respect to the existing literature, it seems that so far the entire
literature has been focusing on prediction errors, and barely talks about how
to build feasible inferential procedures, such as constructing confidence
intervals. A few exceptions known to us are Du et al. (2021), Farrell et al.
(2021) and Chen et al. (2022) for example on estimation and inference for the
average treatment effect rather than on $g(\cdot)$ we are interested in this
paper. In other words, to the best of our knowledge, there have been neither
asymptotic distribution results available for $\widehat{g}(\cdot)$ nor
inferential results available for $g(\cdot)$. Therefore, the main objective
and contribution of this paper is that we develop the LNN based approach to
addressing important estimation and inferential issues for $g(\cdot)$
associated with time series data.

Another challenge which arises with the complexity of the NN architecture is
the transparency of algorithms. It is actually associated with our discussion
about lacking of understanding on how to establish a feasible version
practically. Although a variety of software packages (e.g., Günther and
Fritsch, 2010 and references therein) have been well adopted, the detailed
implementation is still very vague in our view. When writing up this paper, we
cannot find any existing paper using NN approach on social science topics that
provides description about numerical implementation. While we agree with Athey
(2019) that “machine learning approaches will have dramatic impacts on
different fields of social science within a short time frame”, a transparent
algorithm is much needed to ensure that findings are relevant and useful in
practice.

Below, we provide some examples to illustrate our concern about the
transparency of algorithms further. In the literature, the “neuralnet” R
package (see Günther and Fritsch, 2010 for detailed illustration) has been
well adopted. When training a NN, a key parameter is called “hidden” that is a
vector of integers specifying the number of hidden neurons in each layer by R
document, and the package refers to Murata et al. (1994) regarding the choice
of the number of neurons. However, it is worth pointing out that Murata et al.
(1994) use a modified AIC criterion to investigate the case with the number of
neurons (as well as the number of parameters) being finite which is reflected
in their asymptotic development given in the appendix. As a consequence, the
arguments made in Murata et al. (1994) no longer hold when the number of
neurons is diverging. As we are about to show below, having a diverging number
of neurons is the minimum requirement to achieve asymptotic consistency, and
the rate of divergence is also associated with the sample size under a set of
minor conditions. Similar issues also apply to “deepnet” package of R,
“torch.nn.Linear” and “tfl.layers.Linear” of Python, “feedforwardnet” of
Matlab, etc. All in all, the fundamental question behind these is how to
establish NN based estimation and inference for practical implementation.

Having those said, our contributions are four-fold. First, we explore the use
of identification restrictions from a nonparametric regression perspective,
and establish the LNN based estimation method for the effective parameters
under a set of mild conditions. Second, asymptotic distributions are derived
accordingly for inferential purposes, and we show that the LNN method
automatically eliminates the dependence of data when calculating the
asymptotic variances. The finding is important, as one can easily use
different types of wild bootstrap methods to obtain valid inference
practically. Also, it is worth emphasizing that for the first time in the
relevant literature, we are able to rigorously establish asymptotic
distributions for the proposed estimators of the unknown functions using the
LNN approach. Third, for quantitative outcomes, the proposed LNN approach
yields closed-form expressions for the estimators of the parameters of
interest. Last but not least, we examine our theoretical findings through
extensive numerical studies.

The rest of paper is organised as follows. Section 2 provides some basic
assumptions with justifications. Section 3 first introduces the LNN
architecture from a nonparametric regression perspective in Section 3.1, and
then establishes asymptomatic results associated with Models (1.1) and (1.4)
in Section 3.2 and Section 3.3 respectively. Section 3.4 summarizes some key
points outlined in this paper, and further discusses some important issues
associated with LNN. We provide extensive simulation studies in Section 4 to
examine the theoretical findings. Section 5 includes two empirical studies
using Model (1.1) and Model (1.4) respectively. Section 6 concludes. We
provide some additional information about Sigmoidal squasher in Appendix A.1,
some extra simulation results in Appendix A.2, an important extension from
$\mathbf{x}\in[0,1]^{d}$ to $\mathbf{x}\in\mathbb{R}^{d}$ in Appendix A.3,
important lemmas in Appendix A.4, and all the necessary proofs in Appendix
A.5.

##  2 Notation and Assumptions

In this section, we introduce some notation that will be repeatedly used
throughout the rest of this paper, and lay out some necessary assumptions.

We start from mathematical symbols. Vectors and matrices are always expressed
in bold font. Further, $I(\cdot)$ defines the indicator function; $\|\cdot\|$
denotes the Euclidean norm of a vector or the Frobenius norm of a matrix;
$\mathbf{0}_{a}$ and $\mathbf{1}_{a}$ are respectively $a\times 1$ vectors of
zeros and ones for $a\in\mathbb{N}$; for a vector of nonnegative integers
$\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{d})^{\top}\in\mathbb{N}_{0}^{d}$
in which $\mathbb{N}_{0}=0\cup\mathbb{N}$, let
$\boldsymbol{\alpha}!=\alpha_{1}!\cdots\alpha_{d}!$; $\mathtt{c}$,
$\mathtt{C}$ and $O(1)$ always stand for fixed constants, and may be different
at each appearance; $\to_{P}$ and $\to_{D}$ stand for convergence in
probability and convergence in distribution, respectively. For $g(\mathbf{x})$
defined in Model (1.1) and Model (1.4), if the partial derivative of
$g(\mathbf{x})$ exists, we write
$\frac{\partial^{|\boldsymbol{\alpha}|}g(\mathbf{x})}{\partial\mathbf{x}^{\boldsymbol{\alpha}}}=\frac{\partial^{|\boldsymbol{\alpha}|}g(\mathbf{x})}{\partial
x_{1}^{\alpha_{1}}\cdots\partial x_{d}^{\alpha_{d}}}$ for short, where
$|\boldsymbol{\alpha}|=\sum_{j=1}^{d}\alpha_{j}$. We let
$\|g\|_{\infty}=\sup_{\mathbf{x}\in[-a,a]^{d}}|g(\mathbf{x})|$.

Having these symbols in hand, we are ready to proceed. First, we regulate the
function of interest and the activation function. Below, we present a
definition, and then formally state the first assumption of the paper.

######  Definition 2.1 (Continuity).

Let $p=q+s$ for some $q\in\mathbb{N}$ and $0<s\leq 1$. A function
$m\,:\,[-a,a]^{d}\mapsto\mathbb{R}$ is called $(p,\mathscr{C})$-smooth, if for
every $\boldsymbol{\alpha}\in\mathbb{N}_{0}^{d}$ with
$|\boldsymbol{\alpha}|=q$ the partial derivative
$\frac{\partial^{|\boldsymbol{\alpha}|}m(\mathbf{x})}{\partial\mathbf{x}^{\boldsymbol{\alpha}}}$
exists and satisfies that for all $\mathbf{x},\mathbf{z}\in[-a,a]^{d}$,

$\displaystyle\left|\frac{\partial^{|\boldsymbol{\alpha}|}m(\mathbf{x})}{\partial\mathbf{x}^{\boldsymbol{\alpha}}}-\frac{\partial^{|\boldsymbol{\alpha}|}m(\mathbf{z})}{\partial\mathbf{z}^{\boldsymbol{\alpha}}}\right|\leq\mathscr{C}\cdot\|\mathbf{x}-\mathbf{z}\|^{s}.$

Definition 2.1 basically defines a family of sufficiently smooth functions.

######  Assumption 1.

  1. 1.

Let $g(\mathbf{x})$ be $(p,\mathscr{C})$-smooth as in Definition 2.1, and
$\max_{|\boldsymbol{\alpha}|\leq
q}\|\frac{\partial^{|\boldsymbol{\alpha}|}g(\mathbf{x})}{\partial\mathbf{x}^{\boldsymbol{\alpha}}}\|_{\infty}\leq\mathtt{c}.$

  2. 2.

Let Sigmoidal function $\sigma(\cdot)\,:\,\mathbb{R}\to[0,1]$ satisfy that

    1. (a)

$\sigma(\cdot)$ is at least $q+1$ times continuously differentiable with
bounded derivatives.

    2. (b)

A point $u_{\sigma}\in\mathbb{R}$ exists, where all derivatives up to the
order $q$ of $\sigma(\cdot)$ are different from zero.

Assumption 1.1 is widely adopted in the literature of nonparametric regression
(e.g., Li and Racine, 2007). The main point is that each component in Taylor
expansion of $g(\cdot)$ is bounded and also sufficiently smooth.

Assumption 1.2 imposes very limited conditions on the activation function, and
nests a wide class of activation functions commonly used in the relevant
literature as special cases, e.g., Sigmoidal squasher (i.e.,
$\sigma(x)=\frac{1}{1+\exp(-x)}$), Error function (i.e.,
$\sigma(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-w^{2})dw$), etc. We refer the
interested reader to Dubey et al. (2022) for a comprehensive review on
different activation functions. In addition, we note that the activation
function is usually chosen by users, so in general it can be much smoother
than the function of interest (i.e., $g(\mathbf{x})$). A detailed example is
provided in the online supplementary Appendix A.1 to show Assumption 1.2 can
be easily fulfilled in practice. Obviously, Assumption 1.2 rules out ReLU
function, which will be left for future research as explained in the
introduction.

We then regulate the behaviour of the time series involved in (1.5).

######  Assumption 2.

  1. 1.

$\\{(\mathbf{x}_{t},\varepsilon_{t})\,|\,t\in[T]\\}$ are strictly stationary
and $\alpha$–mixing with mixing coefficient

$\displaystyle\alpha(t)=\sup_{A\in\mathcal{F}_{-\infty}^{0},B\in\mathcal{F}_{t}^{\infty}}|P(A)P(B)-P(AB)|$

satisfies $\sum_{t=1}^{\infty}\alpha^{\nu/(2+\nu)}(t)<\infty$ for some
$\nu>0$, which is the same as involved in Assumption 2.3 below, where
$\mathcal{F}_{-\infty}^{0}$ and $\mathcal{F}_{t}^{\infty}$ are the
$\sigma$-algebras generated by $\\{(\mathbf{x}_{s},\varepsilon_{s}):s\leq
0\\}$ and $\\{(\mathbf{x}_{s},\varepsilon_{s}):s\geq t\\}$, respectively.

  2. 2.

The density function $f_{\mathbf{x}}(\cdot)$ of $\mathbf{x}_{1}$ is Lipschitz
continuous on $[-a,a]^{d}$, and is also bounded away from 0 below and bounded
away from $\infty$ above on $[-a,a]^{d}$.

  3. 3.

Model (1.1): $E[\varepsilon_{1}\,|\,\mathbf{x}_{1}]=0$,
$E[\varepsilon_{1}^{2}\,|\,\mathbf{x}_{1}]=\sigma_{\varepsilon}^{2}$ almost
surely (a.s.), and
$E[|\varepsilon_{1}|^{2+\nu}\,|\,\mathbf{x}_{1}]\leq\mathtt{c}$ a.s.

Model (1.4): $\phi_{\varepsilon}(\cdot)$ and $\Phi_{\varepsilon}(\cdot)$ are
known, and
$[1-\Phi_{\varepsilon}(g(\mathbf{x}))]\Phi_{\varepsilon}(g(\mathbf{x}))\neq 0$
uniformly on $[-a,a]^{d}$.

Assumption 2 is standard in the literature of time series regression (Fan and
Yao, 2003; Gao, 2007). As we shall embed the regressors $\\{\mathbf{x}_{t}\\}$
in the activation function under consideration that is bounded from both below
and above, we do not impose many assumptions on the behaviour of
$\\{\mathbf{x}_{t}\\}$. In a sense, LNN automatically rescales the regressors
by design, which somewhat simplifies our asymptotic analysis. Also, we note
that the uncorrelated structure imposed between $\\{\mathbf{x}_{t}\\}$ and
$\\{\varepsilon_{t}\\}$ is not necessary in this paper. For Model (1.1), we in
fact can allow for heterogeneity, such as

$\displaystyle
E[\varepsilon_{1}^{2}\,|\,\mathbf{x}_{1}]=\psi(\mathbf{x}_{1}),$

which however makes notation even more complicated than what we need to
involve. Similarly, for Model (1.4), we may also introduce heterogeneity in a
more general setup such as

$\displaystyle\Pr(y_{1}=1\,|\,\mathbf{x}_{1})=\Phi_{\varepsilon}(g(\mathbf{x}_{1})\,|\,\widetilde{\psi}),$

in which $\widetilde{\psi}$ stands for an unknown variance of
$\varepsilon_{1}$ to be determined by data. Again, it will complex our
discussion and deviate from our main goal.

Finally, we would like to point out that $\\{\mathbf{x}_{t}\\}$ do not have to
be strictly stationary, and can have other more complex structures such as
linear processes, locally stationarity, heterogeneity, deterministic trends,
etc. Surely, the corresponding development and asymptotic results will need to
be modified accordingly, but it will not add too many extra credits to the
original idea of this paper. Therefore, we adopt the current setting
throughout the rest of this paper.

We now move on to construct the so–called LNN architecture for both models,
and will provide some necessary comments along the way.

##  3 Methodology and Theory

In this section, we first introduce the LNN architecture in Section 3.1, and
then establish asymptotic results for both models in Sections 3.2 and 3.3
respectively. Section 3.4 further discusses some issues which may guide our
future research.

One important goal of this paper is to understand the logic of NN. Therefore,
it is worth to be more specific on how NN architecture works conceptually. In
the literature of machine learning, one prefers to target the entire area that
the unknown function is defined on, and normally uses training set to pre-
specify a lot of parameters, which do not change with respect to different
observations in the test set. By doing so, one just needs to pay some price
when calculating the parameters in the first time, and no longer needs to
update them with respect to different values in the test set. Having this
mind, we are ready to consider the estimation of both Models (1.1) and (1.4).

###  3.1 The Localized Neural Network Architecture

In contrast to the current literature that usually allows for all parameters
to be estimated from data, we start by presenting some identification
conditions, which can help reduce the number of effective parameters
significantly. In a sense, we bring in sparsity from a different perspective.

In view of the literature (e.g., Bauer and Kohler, 2019; Schmidt-Hieber,
2020), one key idea behind NN architecture is that it can approximate
polynomial terms, of which as well understood the linear combination can
further approximate unknown functions by standard nonparametric analysis. In
this paper, we do the same, but the difference is that we introduce a
bandwidth parameter $h$ below, which serves the same purpose as that employed
in nonparametric kernel estimation. This is basically why we define the neural
network under consideration as a localized neural network (LNN). More
importantly, the length of $h$ has a direct control on activated and non-
activated neurons involved in the hidden layer. We shall be clear on this
point when presenting the estimation method in Sections 3.2 and 3.3.

To proceed further, we develop several new and useful results in Lemmas
3.1-3.3 to show how to approximate $g(\cdot)$ by a localized neural network
system.

######  Lemma 3.1 (Identification).

Suppose that Assumption 1.2 holds. For $\forall x_{0}\in\mathbb{R}^{d}$, we
can find $\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{q+1})^{\top}$ and
$\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{q+1})^{\top}$ such that

$\displaystyle\sup_{|x-x_{0}|\leq
h}\left|\sum_{k=1}^{q+1}\gamma_{k}\sigma(\beta_{k}\cdot(x-x_{0})+u_{\sigma})-(x-x_{0})^{q}\right|=O\left(h^{q+1}\right)$

holds, where $\mathtt{C}$ is a constant, $h$ is a bandwidth, and
$\gamma_{k}=\frac{(-1)^{q+k-1}{\mathtt{C}}^{q}}{\sigma^{(q)}(u_{\sigma})}\binom{q}{k-1}$
and $\beta_{k}=\frac{k-1}{\mathtt{C}}$ for $k\in[q+1]$.

######  Remark 3.1.

  1. 1.

Note that Lemma 3.1 is independent of data. Only the number of polynomial
terms to be approximated depends on the smoothness of $g(\cdot)$. In fact,
even the smoothness is predetermined by users, which is equivalent to choosing
among local constant, local linear, and local higher order polynomials in the
literature of kernel regression (Fan and Gijbels, 1996). That said,
$\boldsymbol{\gamma}$, $\boldsymbol{\beta}$ and $u_{\sigma}$ are fully decided
by the activation function and the polynomial terms, so they are known prior
to regression.

The constant $\mathtt{C}$ raises an issue of identifiability, so we simply let
$\mathtt{C}=1$ throughout. Also, as required by Assumption 1.2,
$\sigma^{(q)}(u_{\sigma})\neq 0$, which can be easily realized in view of
Figure A.1 of the online supplementary appendix.

  2. 2.

$h$ is equivalent to the bandwidth of kernel regression, so it will be decided
by the sample size eventually. Later on, we show that $h$ reflects the
relationship between the sample size and the total number of neurons, and also
controls how many neurons are activated.

  3. 3.

$k-1$ of $\beta_{k}$ is to invoke the property of the Stirling number of the
second kind in the development. As a result, the LNN architecture with a
suitable Sigmoidal function can approximate the polynomial terms, of which the
linear combination further ensures an unknown smooth function can be
approximated.

  4. 4.

Lemma 3.1 in fact yields a recursive relationship, as one can repeatedly
invoke Lemma 3.1 to replace $(x-x_{0})$ inside the activation. It will then
yield a LNN architecture with multiple hidden layers. However, we do not see
any benefit of doing so unless $g(\cdot)$ in both models has certain specific
structure. That said, we focus on LNN with one hidden layer in this study.

To carry on, we further introduce a few more symbols. Denote that

$\displaystyle\mathscr{P}_{q}=\left\\{\text{Linear span of the monomials
}\prod_{k=1}^{d}x_{k}^{n_{k}}\text{ with }0\leq|\mathbf{n}|\leq q\right\\},$
(3.1)

where $\mathbf{n}=(n_{1},\ldots,n_{d})^{\top}\in\mathbb{N}_{0}^{d}$ and
$|\mathbf{n}|=\sum_{i=1}^{d}n_{i}$. The dimension of $\mathscr{P}_{q}$ is
equivalent to the number of distinctive terms in the expansion of
$(1+x_{1}+\cdots+x_{d})^{q}$, so

$\displaystyle\text{dim}\mathscr{P}_{q}=\binom{d+q}{d}:=d_{q},$

where $d_{q}$ is introduced for notational simplicity, and is fixed due to the
fact that both $d$ and $q$ are finite. Furthermore, we define

$\displaystyle\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})=(m_{1}(\mathbf{x}\,|\,\mathbf{x}_{0}),\ldots,m_{d_{q}}(\mathbf{x}\,|\,\mathbf{x}_{0}))^{\top},$
(3.2)

where $m_{j}(\mathbf{x}\,|\,\mathbf{x}_{0})$’s are the basis monomials
(centred at $\mathbf{x}_{0}$) in $\mathscr{P}_{q}$. Also, we denote a set
$C_{\mathbf{x}_{0},h}$:

$C_{\mathbf{x}_{0},h}=\\{\mathbf{x}\,|\,|x_{j}-x_{0,j}|\leq h\text{ for
}j\in[d]\\},$ (3.3)

where $x_{j}$ and $x_{0,j}$ are the $j^{th}$ elements of $\mathbf{x}$ and
$\mathbf{x}_{0}$ respectively.

With these symbols in hand, we present the next lemma, which demonstrates the
feasibility of the LNN architecture.

######  Lemma 3.2 (Feasibility).

Suppose that Assumption 1.2 holds. For
$\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$, we define
$p(\mathbf{x}\,|\,\mathbf{x}_{0},\boldsymbol{\lambda})=\boldsymbol{\lambda}^{\top}\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0}),$
where $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{d_{q}})^{\top}$ with
$\|\boldsymbol{\lambda}\|\leq\mathtt{c}$. Then a LNN of the type

$\displaystyle
s(\mathbf{x}\,|\,\mathbf{x}_{0},\widetilde{\boldsymbol{\lambda}})=(\widetilde{\boldsymbol{\lambda}}\otimes\boldsymbol{\gamma})^{\top}\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})\quad\text{with}\quad\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})=\\{\sigma\left([1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top}]\boldsymbol{\pi}_{j}\right)\\}_{d_{q}(q+1)\times
1}$

exists such that

$\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}|s(\mathbf{x}\,|\,\mathbf{x}_{0},\widetilde{\boldsymbol{\lambda}})-p(\mathbf{x}\,|\,\mathbf{x}_{0},\boldsymbol{\lambda})|=O(h^{q+1}),$

where $\widetilde{\boldsymbol{\lambda}}=\mathbf{D}^{\top}\boldsymbol{\lambda}$
with $\mathbf{D}$ being a rotation matrix, and

$\displaystyle(\boldsymbol{\pi}_{1},\ldots,\boldsymbol{\pi}_{d_{q}(q+1)})$
$\displaystyle=$
$\displaystyle\frac{\operatorname*{\textrm{diag}}\\{h,\mathbf{I}_{d}\\}\mathbf{W}(\mathbf{I}_{d_{q}}\otimes\boldsymbol{\beta}^{\top})}{d+1}+\bigg{[}\begin{matrix}u_{\sigma}\\\
\mathbf{0}_{d}\end{matrix}\bigg{]}\otimes\mathbf{1}_{d_{q}(q+1)}^{\top},$

in which $\mathbf{W}$ is a user chosen matrix satisfying that
$\max_{j}\|\mathbf{w}_{j}\|\leq\sqrt{d+1}$ and
$\mathbf{w}_{j_{1}}\neq\mathbf{w}_{j_{2}}$ for any given
$j_{1},j_{2}\in[d_{q}]$, and $\mathbf{w}_{j}$ stands for the $j^{th}$ column
of $\mathbf{W}$.

######  Remark 3.2.

  1. 1.

Note that we use vector form to rewrite each element of
$\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})$, i.e.,

$\displaystyle\sigma\left([1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top}]\boldsymbol{\pi}_{j}\right)$

which is exactly what Sigmoidal activation function does (i.e., mapping a
liner combination of regressors plus a location parameter to $[0,1]$). By
design, $\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})$ naturally and
automatically explores different interaction terms of the regressors.

  2. 2.

By choosing $\boldsymbol{\lambda}$ and $\mathbf{x}_{0}$,
$p(\mathbf{x}\,|\,\mathbf{x}_{0},\boldsymbol{\lambda})$ can approximate
$g(\mathbf{x})$ reasonably well in a small neighbourhood of $\mathbf{x}_{0}$.
This is not hard to see, as the leading terms of the $q^{th}$ order Taylor
expansion of $g(\mathbf{x})$ at $\mathbf{x}_{0}$ can be written as follows:

$\displaystyle g(\mathbf{x})\simeq\sum_{0\leq|\mathbf{J}|\leq
q}\frac{1}{\mathbf{J}!}\cdot\frac{\partial^{|\mathbf{J}|}g(\mathbf{x}_{0})}{\partial\mathbf{x}^{\mathbf{J}}}(\mathbf{x}-\mathbf{x}_{0})^{\mathbf{J}}:=\boldsymbol{\lambda}^{\top}\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0}):=p(\mathbf{x}\,|\,\mathbf{x}_{0},\boldsymbol{\lambda}),$

where the definition of $\boldsymbol{\lambda}$ should be obvious in view of
the definition of $\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})$ according to
(3.2). As a result,
$s(\mathbf{x}\,|\,\mathbf{x}_{0},\widetilde{\boldsymbol{\lambda}})$
essentially approximates $g(\mathbf{x})$ but in a small neighbourhood only,
and the rate produced in Lemma 3.2 is proportional to $h^{q+1}$ as a bias
term.

  3. 3.

Lemma 3.2 infers that the LNN architecture rotates the parameters of interest
(i.e., $\boldsymbol{\lambda}$) with a pre-determined $d_{q}\times d_{q}$ full
rank matrix $\mathbf{D}$, so the new parameters of interest become
$\widetilde{\boldsymbol{\lambda}}$. Therefore, accounting for identification
conditions properly ensures the number of effective parameters to be estimated
is far less than it looks like. Without addressing identification issues, one
has to determine $d_{q}(q+1)(d+2)$ parameters. In other words, utilizing
identification conditions automatically reduces the number of effective
parameters which need to be estimated.

  4. 4.

Without loss of generality, we can let
$\mathbf{w}_{j}=\frac{\sqrt{d+1}}{q}\mathbf{w}_{j}^{*}$, where
$\\{\mathbf{w}_{j}^{*}\\}$ are the vectors corresponding to the powers of the
distinctive terms in the expansion of $(1+x_{1}+\cdots+x_{d})^{q}$. Thus,
$\|\mathbf{w}_{j}\|\leq\frac{\sqrt{d+1}}{q}|\mathbf{w}_{j}^{*}|=\sqrt{d+1},$
which further ensures Lemma 3.1 can be invoked.

Using Lemma 3.2, the final form of LNN requires us to focus on a compact
subspace of $\mathbb{R}^{d}$ (i.e., $[-a,a]^{d}$), as we need to partition it
into lots of small cubes. As discussed in Appendix A.3 below, the main results
remain valid when allowing $a\rightarrow\infty$ along with the sample size.
These small cubes may have different names in each appearance. For example,
they are referred to as (hyper-)cubes in Bauer and Kohler (2019) and Schmidt-
Hieber (2020), and are referred to as localization in Farrell et al. (2021).
In our view, each cube is corresponding to an effective sample set, so it is
very much similar to the effective sample range in nonparametric kernel
regression, which is jointly determined by the point of interest and the
bandwidth.

That said, for a given integer $M\geq 1$, we subdivide $[-a,a]^{d}$ into
$M^{d}$ cubes of side length as follows:

$\displaystyle h=\frac{a}{M},$

and for comprehensibility, we number these cubes by
$C_{\mathbf{x}_{0\mathbf{i}},h}$ with $\mathbf{i}\in[M]^{d}$, where
$\mathbf{x}_{0\mathbf{i}}$ represents the center of
$C_{\mathbf{x}_{0\mathbf{i}},h}$. Mathematically, we can write

$\displaystyle
C_{\mathbf{x}_{0\mathbf{i}},h}=\\{\mathbf{x}\,|\,|x_{j}-x_{0\mathbf{i},j}|\leq
h\text{ for }j\in[d]\\},$ (3.4)

where $x_{j}$ and $x_{0\mathbf{i},j}$ are the $j^{th}$ elements of
$\mathbf{x}$ and $\mathbf{x}_{0\mathbf{i}}$ respectively. Then the following
lemma holds.

######  Lemma 3.3 (NN Architecture).

Suppose that Assumption 1 holds. There exists a LNN approximate of
$g(\mathbf{x})$ as follows:

$\displaystyle\|g(\mathbf{x})-\widetilde{s}(\mathbf{x}\,|\,\widetilde{\boldsymbol{\Lambda}})\|_{\infty}=O(h^{p}),$

where
$\widetilde{\boldsymbol{\Lambda}}=\\{\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}\,|\,\mathbf{i}\in[M]^{d}\\}$,
and

$\displaystyle\widetilde{s}(\mathbf{x}\,|\,\widetilde{\boldsymbol{\Lambda}})=\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
s(\mathbf{x}\,|\,\mathbf{x}_{0\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})\quad\text{with}\quad
I_{\mathbf{i},h}(\mathbf{x})=I(\mathbf{x}\in C_{\mathbf{x}_{0\mathbf{i}},h}).$

######  Remark 3.3.

  1. 1.

LNN considers the estimation of
$\\{\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}\,|\,\mathbf{i}\in[M]^{d}\\}$
simultaneously, where $\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$
corresponds to $\boldsymbol{\lambda}_{\mathbf{i}}$ up to a rotation matrix
$\mathbf{D}$, and $\boldsymbol{\lambda}_{\mathbf{i}}$ is decided by the Taylor
expansion of $g(\mathbf{x})$ at the point $\mathbf{x}_{0\mathbf{i}}$.

  2. 2.

The use of the indicator function is consistent with the sparsity setting of
Schmidt-Hieber (2020), in which the author argues that “the network sparsity
assumes that there are only few non-zero/active network parameters”. In fact,
our study further clarifies the definition of “the network sparsity” by
providing two categories: (1) defining non-active neurons (such as those in
Figure 1) which is realised through the use of indicator function
$I_{\mathbf{i},h}(\mathbf{x})$, and (2) pointing out the number of effective
parameters. In particular, the second point is important in our view, as one
can now define a true set of parameters when using thresholding techniques as
in Wang and Lin (2021), and Fan and Gu (2022).

  3. 3.

The bandwidth $h$ (i.e., the number of cubes) is data driven. This point will
become obvious when we move on to Sections 3.2 and 3.3 below.

###  3.2 On Model (1.1)

In this subsection, we consider Model (1.1) only. As explained previously, LNN
targets the entire area that $g(\mathbf{x})$ is defined on, and uses all
training data to pre-specify a large number of parameters which do not change
with respect to the test set.

In view of Lemma 3.3, we define the following set of LNN candidates:

$\displaystyle\mathcal{S}=\left\\{\widetilde{s}(\mathbf{x}\,|\,\boldsymbol{\Theta})=\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
s(\mathbf{x}\,|\,\mathbf{x}_{0\mathbf{i}},\boldsymbol{\theta}_{\mathbf{i}})\,\Big{|}\,\|\boldsymbol{\theta}_{\mathbf{i}}\|<\infty\right\\},$
(3.5)

in which
$\boldsymbol{\Theta}=\\{\boldsymbol{\theta}_{\mathbf{i}}\,|\,\mathbf{i}\in[M]^{d}\\}$
with $\boldsymbol{\theta}_{\mathbf{i}}$’s being $d_{q}\times 1$ vectors. The
objective function is defined as

$Q(\boldsymbol{\Theta})=\sum_{t=1}^{T}[y_{t}-\widetilde{s}(\mathbf{x}_{t}\,|\,\boldsymbol{\Theta})]^{2}.$
(3.6)

Accordingly the OLS estimator of $\widetilde{\boldsymbol{\Lambda}}$ defined in
Lemma 3.3 is obtained by

$\displaystyle\widehat{\boldsymbol{\Theta}}=\operatorname*{\arg\\!\min}_{\widetilde{s}\in\mathcal{S}}Q(\boldsymbol{\Theta})\quad\text{with}\quad\widehat{\boldsymbol{\Theta}}=\\{\widehat{\boldsymbol{\theta}}_{\mathbf{i}}\,|\,\mathbf{i}\in[M]^{d}\\}$
(3.7)

and, for $\forall\mathbf{x}_{0}\in[-a,a]^{d}$, the estimator of
$g(\mathbf{x}_{0})$ is then defined by

$\displaystyle\widehat{g}(\mathbf{x}_{0})=\widetilde{s}(\mathbf{x}_{0}\,|\,\widehat{\boldsymbol{\Theta}}).$
(3.8)

By design, LNN accounts for the points on the boundary of $[-a,a]^{d}$
automatically.

######  Remark 3.4.

Our investigation assumes that the number of regressors is correctly
specified. When they are over-specified, sparsity naturally arises. It is now
even clearer why thresholding methods (such as those in Wang and Lin, 2021;
Fan and Gu, 2022) should be adopted, and what they are really penalizing. In a
sense identification restrictions bridge LNN and thresholding approaches in a
better fashion. In this work, we do not further explore regression with
penalty terms in order not to deviate from our main goal.

Equation (3.7) admits a closed–form estimator for each
$\widehat{\boldsymbol{\theta}}_{\mathbf{i}}$. To see this, we write

$\displaystyle\frac{\partial
Q(\boldsymbol{\Theta})}{\partial\boldsymbol{\theta}_{\mathbf{i}}}$
$\displaystyle=$
$\displaystyle-2\sum_{t=1}^{T}[y_{t}-\widetilde{s}(\mathbf{x}_{t}\,|\,\boldsymbol{\Theta})]\cdot
I_{\mathbf{i},h}(\mathbf{x}_{t})\cdot\frac{\partial
s(\mathbf{x}_{t}\,|\,\mathbf{x}_{0\mathbf{i}},\boldsymbol{\theta}_{\mathbf{i}})}{\partial\boldsymbol{\theta}_{\mathbf{i}}}$
$\displaystyle=$
$\displaystyle-2\sum_{t=1}^{T}[y_{t}-\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\boldsymbol{\theta}_{\mathbf{i}}]\cdot\widetilde{\mathbf{x}}_{\mathbf{i},t},$

where
$\widetilde{\mathbf{x}}_{\mathbf{i},t}=I_{\mathbf{i},h}(\mathbf{x}_{t})(\mathbf{I}_{d_{q}}\otimes\boldsymbol{\gamma}^{\top})\boldsymbol{\sigma}(\mathbf{x}_{t}\,|\,\mathbf{x}_{0\mathbf{i}})$,
and the second equality follows from the fact that
$I_{\mathbf{i},h}(\mathbf{x}_{t})I_{\mathbf{j},h}(\mathbf{x}_{t})=0$ for
$\mathbf{i}\neq\mathbf{j}$. Thus, for $\forall\mathbf{i}$, the first order
condition yields

$\displaystyle\widehat{\boldsymbol{\theta}}_{\mathbf{i}}$ $\displaystyle=$
$\displaystyle\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}y_{t}.$
(3.9)

######  Remark 3.5.

We note that by (A.11) and (A.12) of the online supplementary appendix:

$\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0\mathbf{i}},h}}\left\|\mathbf{D}^{-1}\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0\mathbf{i}})-(\mathbf{I}_{d_{q}}\otimes\boldsymbol{\gamma}^{\top})\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0\mathbf{i}})\right\|=O(h^{q+1}),$

where $\mathbf{D}$ is a full rank predetermined rotation matrix, and has been
explained in Remark 3.3. In connection with the definition of
$\widetilde{\mathbf{x}}_{\mathbf{i},t}$, we immediately obtain

$\displaystyle\left\|\mathbf{D}^{-1}I_{\mathbf{i},h}(\mathbf{x}_{t})\cdot\mathbf{m}(\mathbf{x}_{t}\,|\,\mathbf{x}_{0\mathbf{i}})-\widetilde{\mathbf{x}}_{\mathbf{i},t}\right\|=O_{P}(h^{q+1}).$

Thus, $\widetilde{\mathbf{x}}_{\mathbf{i},t}$ is equivalent to
$I_{\mathbf{i},h}(\mathbf{x}_{t})\cdot\mathbf{m}(\mathbf{x}_{t}\,|\,\mathbf{x}_{0\mathbf{i}})$
up to a predetermined rotation matrix $\mathbf{D}$. The finding is consistent
with Remark 3.2, and ensures the invertibility of
$\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}$,
at least asymptotically.

After carefully studying (3.9) for each $\mathbf{i}$ and repeatedly invoking
$I_{\mathbf{i},h}(\mathbf{x}_{t})I_{\mathbf{j},h}(\mathbf{x}_{t})=0$ for
$\mathbf{i}\neq\mathbf{j}$, the first theorem of this paper is stated as
follows.

######  Theorem 3.1.

Suppose that Assumptions 1 and 2 hold. For
$\forall\mathbf{x}_{0}\in[-a,a]^{d}$,

$\displaystyle\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{x}_{0}}^{-1}(\widehat{g}(\mathbf{x}_{0})-g(\mathbf{x}_{0})+O_{P}(h^{p}))\to_{D}N(0,1),$

where
$\widehat{\sigma}_{\mathbf{x}_{0}}^{2}=\sigma_{\varepsilon}^{2}\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})^{\top}\mathbf{H}\boldsymbol{\Sigma}_{\mathbf{x}_{0}}^{-1}\mathbf{H}\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})$,
$\boldsymbol{\Sigma}_{\mathbf{x}_{0}}=f_{\mathbf{x}}(\mathbf{x}_{0})\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}$,
$\mathbf{H}=\operatorname*{\textrm{diag}}\\{H_{1},\ldots,H_{d_{q}}\\}$ with
$H_{j}=\prod_{k=1}^{d}h^{-n_{j,k}}$, and
$\mathbf{n}_{j}=(n_{j,1},\ldots,n_{j,d})^{\top}$ includes the corresponding
power terms of $m_{j}(\mathbf{x}\,|\,\mathbf{x}_{0})$ defined in (3.2).

It is noted that our LNN based estimation method is simple and easy to
implement. As pointed out in the introduction, moreover, we are probably among
the first being able to rigorously establish the asymptotic normality of the
LNN based estimator given by $\widehat{g}(\mathbf{x}_{0})$ in Theorem 3.1.

It is also noted out that our simulation studies in Appendix A.2 show that
$\widehat{g}(\mathbf{x}_{0})$ has better finite–sample properties than those
of the conventional local–constant kernel estimation method. The reason to
explain the finite–sample advantages is probably that the LNN based estimation
method naturally makes the best use of all possible linear combinations of
$\mathbf{x}_{t}$ as explained in Remark 3.2.1 above. By contrast, the
conventional local–constant kernel estimation method only evaluates the
estimates at the single data points: $\mathbf{x}_{t}$.

The same comments apply to Theorem 3.3 below.

######  Remark 3.6.

  1. 1.

When establishing the asymptotic distribution, the terms
$E[\varepsilon_{1}\varepsilon_{1+t}]$ for $t\geq 1$ all vanish in the
asymptotic covariance matrix due to the partition of (3.4) and the use of the
indicator function. Here, the indicator function is a special nonparametric
kernel function, so it is not very surprising that the finding is consistent
with Chapter 5 of Fan and Yao (2003), for example, in the literature of
nonparametric kernel regression with dependent data. As a result, when
constructing confidence interval, one can turn to wild bootstrap procedure
directly, which also bypasses the rotation matrix $\mathbf{D}$ (see Theorem
3.2 below for details).

  2. 2.

Ideally, we would like to replace the indicator function by a general kernel
function such as $K(\cdot)$ in (3.5). However, by doing so,
$s(\mathbf{x}\,|\,\mathbf{x}_{0},\widetilde{\boldsymbol{\lambda}})$ will not
be able to approximate $g(\mathbf{x}_{0})$ using Lemma 3.2.

  3. 3.

We do require $g(\mathbf{x})$ to be defined on a compact set by the design of
LNN, but do not impose restriction on the range of $\\{\mathbf{x}_{t}\\}$. In
Appendix A.3, we further explain how to get a valid result for
$\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$ from a nonparametric regression
perspective.

To close this subsection, we propose the following bootstrap procedure to
establish inference in practice.

  1. 1.

By (3.8), we calculate
$\widehat{\varepsilon}_{t}=y_{t}-\widehat{g}(\mathbf{x}_{t})$ for $t\in[T]$.

  2. 2.

Collect i.i.d. draws of $\\{\eta_{t}\,|\,t\in[T]\\}$ from $N(0,1)$, and
construct the bootstrap version dependent variables as follows:
$y_{t}^{*}=\widehat{g}(\mathbf{x}_{t})+\widehat{\varepsilon}_{t}\eta_{t}$. We
re-estimate $g(\cdot)$ using $\\{(y_{t}^{*},\mathbf{x}_{T})\,|\,t\in[T]\\}$ as
in (3.8), and denote the estimate as $\widehat{g}^{*}(\cdot)$.

  3. 3.

Repeat Step 2 $R$ times, where $R$ is sufficiently large.

For the bootstrap procedure, the following result holds immediately.

######  Theorem 3.2.

Let the conditions of Theorem 3.1 hold. Suppose further that $Th^{d+2p}\to 0$.
For $\forall\mathbf{x}_{0}\in[-a,a]^{d}$, we have

$\displaystyle\sup_{w}|\text{Pr}^{*}(\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{x}_{0}}^{-1}[\widehat{g}^{*}(\mathbf{x}_{0})-\widehat{g}(\mathbf{x}_{0})]\leq
w)-\Pr(\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{x}_{0}}^{-1}[\widehat{g}(\mathbf{x}_{0})-g(\mathbf{x}_{0})]\leq
w)|=o_{P}(1).$

where $\text{Pr}^{*}$ is the probability measure induced by the bootstrap
procedure.

The condition $Th^{d+2p}\to 0$ is required to eliminate the bias terms. For
the LNN estimators, the bias terms arise from two resources: (1) LNN
approximates each polynomial term, and (2) the polynomial terms approximate
$g(\cdot)$. As a result, we could not derive the detailed form of the biases
but providing its order only. That said, to ensure the bootstrap procedure can
provide a valid confidence interval in practice, we require $Th^{d+2p}\to 0$.

###  3.3 On Model (1.4)

In this subsection, we investigate Model (1.4), which has been widely adopted
for a variety of decision making processes. See Athey (2019) for discussions
on different examples.

Direct calculation shows that

$\displaystyle\Pr(y=1\,|\,\mathbf{x})=\Phi_{\varepsilon}(g(\mathbf{x})),$
$\displaystyle\Pr(y=0\,|\,\mathbf{x})=1-\Phi_{\varepsilon}(g(\mathbf{x})),$

which yields $E[y\,|\,\mathbf{x}]=\Phi_{\varepsilon}(g(\mathbf{x}))$. Thus,
provided a set of time series observations given in (1.5), the likelihood
function is specified as follows:

$\displaystyle
L(g)=\prod_{t=1}^{T}[1-\Phi_{\varepsilon}(g(\mathbf{x}_{t}))]^{1-y_{t}}\cdot\Phi_{\varepsilon}(g(\mathbf{x}_{t}))^{y_{t}}.$

Accordingly, the log-likelihood function is defined below:

$\displaystyle\log L(g)$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{T}l_{t}(g(\mathbf{x}_{t}))$ (3.10) $\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\left\\{(1-y_{t})\cdot\log[1-\Phi_{\varepsilon}(g(\mathbf{x}_{t}))]+y_{t}\cdot\log\Phi_{\varepsilon}(g(\mathbf{x}_{t}))\right\\},$

where the definition of $l_{t}(\cdot)$ is obvious.

Our interest is still inferring $g(\cdot)$. We construct the neural network as
in (3.5), and consider the following objective function:

$\displaystyle\log L(\widetilde{s}(\cdot\,|\,\boldsymbol{\Theta}))$
$\displaystyle=$ $\displaystyle\sum_{t=1}^{T}\log
l_{t}(\widetilde{s}(\mathbf{x}_{t}\,|\,\boldsymbol{\Theta}))$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\left\\{(1-y_{t})\cdot\log[1-\Phi_{\varepsilon}(\widetilde{s}(\mathbf{x}_{t}\,|\,\boldsymbol{\Theta}))]+y_{t}\cdot\log\Phi_{\varepsilon}(\widetilde{s}(\mathbf{x}_{t}\,|\,\boldsymbol{\Theta}))\right\\},$

which yields the following maximum likelihood estimator:

$\displaystyle\widehat{\boldsymbol{\Theta}}=\operatorname*{\arg\\!\max}_{\widetilde{s}\in\mathcal{S}}\log
L(\widetilde{s}(\cdot\,|\,\boldsymbol{\Theta})).$ (3.11)

Note that our LNN method involves a general unknown function form. As a
result, the classical results of likelihood estimation, such as those in Newey
and McFadden (1994), no longer hold. Therefore, before establishing an
asymptotic distribution using $\widehat{\boldsymbol{\Theta}}$, we state a
lemma to show the overall feasibility of the LNN architecture when modelling
binary outcomes.

######  Lemma 3.4.

Under Assumptions 1 and 2,

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}[\Phi_{\varepsilon}(g(\mathbf{x}_{t}))-\Phi_{\varepsilon}(\widetilde{s}(\mathbf{x}_{t}\,|\,\widehat{\boldsymbol{\Theta}}))]^{2}=o_{P}(1).$

Lemma 3.4 provides the consistency, and also bridges the likelihood estimation
and the nonlinear least squares approach to some extent. To be precise, Lemma
3.4 does not provide any specific consistency for
$\forall\,\widehat{\boldsymbol{\theta}}_{\mathbf{i}}$. Instead, it evaluates
the overall performance of LNN. More importantly, it says when modelling a
binary outcome, the likelihood estimation using the LNN architecture is
approximately equivalent to implementing a nonlinear least squares method
provided the distribution of $\varepsilon_{t}$ is correctly specified. In
addition, Lemma 3.4 further infers that

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widehat{\boldsymbol{\Theta}})]^{2}=o_{P}(1),$

which is nice, as many remarks made in Section 3.2 can be directly applied.
Last but not least, Lemma 3.4 facilitates numerical implementation in
practice, which will be further discussed in Section 4.2.

Below, we establish the following asymptotic distribution for Model (1.4).

######  Theorem 3.3.

Suppose that Assumptions 1 and 2 hold. For
$\forall\mathbf{x}_{0}\in[-a,a]^{d}$,

$\displaystyle\sqrt{Th^{d}}\widetilde{\sigma}_{\mathbf{x}_{0}}^{-1}(\widehat{g}(\mathbf{x}_{0})-g(\mathbf{x}_{0})+O_{P}(h^{p}))\to_{D}N(0,1),$

where $\widehat{g}(\mathbf{x})$ is defined in the same form as (3.8),
$\widetilde{\sigma}_{\mathbf{x}_{0}}^{2}=\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})^{\top}\mathbf{H}\widetilde{\boldsymbol{\Sigma}}_{\mathbf{x}_{0}}^{-1}\mathbf{H}\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})$,
and

$\displaystyle\widetilde{\boldsymbol{\Sigma}}_{\mathbf{x}_{0}}=\frac{f_{\mathbf{x}}(\mathbf{x}_{0})\phi_{\varepsilon}(g(\mathbf{x}_{0}))^{2}}{[1-\Phi_{\varepsilon}(g(\mathbf{x}_{0}))]\Phi_{\varepsilon}(g(\mathbf{x}_{0}))}\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}.$

Remark 3.6 can be applied to Theorem 3.3 with some obvious modifications.

In light of Patrick and Andres (2012), we then propose a score based wild
bootstrap approach for inferential purposes as follows.

  1. 1.

For each bootstrap replication, we collect i.i.d. draws of
$\\{\eta_{t}\,|\,t\in[T]\\}$ from $N(0,1)$, and calculate

$\displaystyle\widehat{\boldsymbol{\theta}}_{\mathbf{i}}^{*}=\widehat{\boldsymbol{\theta}}_{\mathbf{i}}+\left(\sum_{t=1}^{T}\frac{\partial^{2}\log
l_{t}(\widetilde{s}(\mathbf{x}_{t}\,|\,\widehat{\boldsymbol{\Theta}}))}{\partial\boldsymbol{\theta}_{\mathbf{i}}\partial\boldsymbol{\theta}_{\mathbf{i}}^{\top}}\right)^{-1}\sum_{t=1}^{T}\frac{\partial\log
l_{t}(\widetilde{s}(\mathbf{x}_{t}\,|\,\widehat{\boldsymbol{\Theta}}))}{\partial\boldsymbol{\theta}_{\mathbf{i}}}\eta_{t},$
(3.12)

where $l_{t}(\cdot)$ is defined in (3.10).

  2. 2.

Repeat Step 1 $R$ times, where $R$ is sufficiently large.

It is worth pointing out that the above procedure is computationally efficient
in the sense that the right hand side of (3.12) enjoys a closed–form
expression, which is given in (A.29) and (A.30) specifically for the sake of
space. Practically, the bootstrap procedure may require much less time
compared with the estimation of (3.11) itself. We will further comment on the
computational issues associated with Model (1.4) in Section 4.2.

The following theorem holds for the above bootstrap procedure.

######  Theorem 3.4.

Let the conditions of Theorem 3.3 hold. Suppose further that $Th^{d+2p}\to 0$.
For $\forall\mathbf{x}_{0}\in[-a,a]^{d}$, we have

$\displaystyle\sup_{w}|\text{Pr}^{*}(\sqrt{Th^{d}}\widetilde{\sigma}_{\mathbf{x}_{0}}^{-1}[\widehat{g}^{*}(\mathbf{x}_{0})-\widehat{g}(\mathbf{x}_{0})]\leq
w)-\Pr(\sqrt{Th^{d}}\widetilde{\sigma}_{\mathbf{x}_{0}}^{-1}[\widehat{g}(\mathbf{x}_{0})-g(\mathbf{x}_{0})]\leq
w)|=o_{P}(1),$

where $\text{Pr}^{*}$ is the probability measure induced by the bootstrap
procedure, and $\widehat{g}^{*}(\cdot)$ is yielded by the bootstrap draws in
an obvious manner.

The comments and discussion below Theorem 3.2 apply to Theorem 3.4 with some
minor modifications, so we do not repeat them here.

###  3.4 Further Discussion

Up to this point, we would like to point out that those questions raised in
Section 1 have all been answered, so Figure 1 can be understood better. We
summarize some key points which may have been discussed previously here and
there, and further discuss some issues left behind.

Neurons — Having established the results in Sections 3.2 and 3.3, it is now
clear that for LNN, the total number of neurons is

$\displaystyle M^{d}\cdot d_{q}(q+1)=\left(\frac{2a}{h}\right)^{d}\cdot
d_{q}(q+1),$

of which only $d_{q}(q+1)$ neurons are activated when loading test data. Among
the activated neurons, the number of effective parameters is only $d_{q}$,
while the rest of the parameters are predetermined.

Multiple Hidden Layers — Lemma 3.1 yields a recursive relationship, as one can
repeatedly invoke Lemma 3.1 to replace $(x-x_{0})$ inside the activation. It
will then yield a LNN architecture with multiple hidden layers. Although
having multiple hidden layers is achievable, at this stage it is not clear to
us why we should do so. As discussed in Remark 3.1, this step is completely
independent of data, so we do not see any benefit of doing so unless
$g(\cdot)$ in both models has certain specific structure. Under some extra
structure on $g(\cdot)$, however, the necessity of developing a LNN
architecture with multiple hidden layers deserves extra attention in future
research.

Dependence — As shown for both models, LNN automatically eliminates the
correlation of observations from different time periods when establishing the
asymptotic distribution. This finding is important, as one can easily use
different types of wild bootstrap methods to obtain valid inference.

ReLU Function — The current study focuses on Sigmoidal activation function. We
conjecture for ReLU function a corresponding LNN approach can be constructed
similarly, which, however, requires some careful consideration.

Thresholding — Our investigation assumes that the number of regressors is
correctly specified. When they are over-specified, sparsity naturally arises.
In this case, it is clear that why a thresholding method (such as those in
Wang and Lin, 2021; Fan and Gu, 2022) should be adopted, and what it is really
penalizing. Utilizing identification conditions allows one to define the true
set of parameters, so dimension reduction techniques such as LASSO
(Tibshirani, 1996) and bridge estimation (Huang et al., 2008) techniques can
be adopted straight away.

Trending — In this paper, we assume that the regressors
$\\{\mathbf{x}_{t}\,|\,t\in[T]\\}$ are strictly stationary and mixing. In
fact, they can have other more complex structures, such as linear processes,
locally stationarity, heterogeneity, deterministic trends, etc. As a result,
many climate models (such as those in Mudelsee, 2019) may be better captured.
For such cases, one may need to revise the assumptions and proofs accordingly
depending on detailed research questions.

##  4 Simulation

In this section we conduct simulations to examine the theoretical findings. We
consider Models (1.1) and (1.4) in Sections 4.1 and 4.2 respectively.

Before proceeding further, we would like to point out the following simulation
design focuses on coverage rates for the entire test set. Also, no statistical
software package is required, as the algorithm becomes obvious after
presenting Section 3. We now provide some details of the numerical
implementation. As shown in Section 3, many parameters are involved in the LNN
architecture. It would be extremely difficult to systematically check every
single one in one paper, so we have to be selective. That said, the following
quantities are pre-fixed without loss of generality.

  * •

Throughout, we use the Sigmoidal squasher, $\sigma(w)=1/(1+\exp(-w))$, as the
activation function. The derivatives of Sigmoidal squasher are provided in
Appendix A.1 of the online supplementary file for the sake of space.

  * •

$\boldsymbol{\pi}_{j}$’s are generated in exactly the same way as mentioned in
Remark 3.2.

  * •

Let $R=200$ for the bootstrap procedure.

  * •

Let $g(\mathbf{x})=1+\sin(\mathbf{x}^{\top}\mathbf{1}_{d}/d)$.

Surely, we can explore different settings, but these are less important based
on the development of Section 3. Also, we are constrained by computing power.

###  4.1 Model (1.1)

Consider the following regression model:

$y_{t}=g(\mathbf{x}_{t})+\varepsilon_{t},$ (4.1)

where $\varepsilon_{t}=0.5\varepsilon_{t-1}+N(0,0.75)$, and the $j^{th}$
element of $\mathbf{x}_{t}$ is generated as $x_{t,j}\sim U(-a,a)$. The
bandwidth $h$ is set as $h=a/M$, where $M$ is the integer closest to $a/h_{1}$
with $h_{1}=2.5\cdot T^{-1/(d+2p-0.5)}$. Here, $-0.5$ is to ensure
$\sqrt{T}h^{p+d/2}\to 0$ holds. In fact, $h$ is very close to $h_{1}$, and the
current setup is simply to guarantee $M$ is a large positive integer. When
designing the simulations, our impression is that the results are not
sensitive to the choices of $g(\mathbf{x})$ and the bandwidth. Therefore, in
what follows, we only vary the values of

$\displaystyle(q,d,u_{\sigma}).$ (4.2)

To measure the finite sample performance, we select $L^{d}$ test points from
$[-a,a]^{d}$ as follows:

$\displaystyle\mathbf{x}_{L,\mathbf{j}}=\left(-a+\frac{2a}{L-1}(j_{1}-1),\ldots,-a+\frac{2a}{L-1}(j_{d}-1)\right)^{\top},$
(4.3)

where $\mathbf{j}\in[L]^{d}$. With each dataset, we first estimate all
$g(\mathbf{x}_{L,\mathbf{j}})$ using the approach of Section 3.2, and then
construct the corresponding 95% confidence interval using the bootstrap
procedure documented in Theorem 3.2 for each point. After $n$ replications, we
report

$\text{RMSE}_{g}=\left\\{\frac{1}{nL^{d}}\sum_{i=1}^{n}\sum_{\mathbf{j}\in[L]^{d}}[\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})-g(\mathbf{x}_{L,\mathbf{j}})]^{2}\right\\}^{1/2}$

and

$\text{CR}_{g}=\frac{1}{nL^{d}}\sum_{i=1}^{n}\sum_{\mathbf{j}\in[L]^{d}}I(\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})-g(\mathbf{x}_{L,\mathbf{j}})\in\text{CI}_{i,\mathbf{j}}),$

where $\widehat{g}_{i}(\cdot)$ stands for the estimate of $g(\cdot)$ at the
$i^{th}$ replication, and $\text{CI}_{i,\mathbf{j}}$ is the 95% confidence
interval of
$\widehat{g}_{i}^{*}(\mathbf{x}_{L,\mathbf{j}})-\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})$
based on the bootstrap draws from the $i^{th}$ replication. We set
$T\in\\{800,1600,2400\\}$ to examine the convergence of NN approach, and
let222We provide some extra simulation results with larger $d$ value in the
online supplementary appendices of the paper, and comment on some issues
related to computational overhead. $u_{\sigma}\in\\{-0.5,0.5\\}$ and
$d\in\\{2,3\\}$. We further set $a=3$, $L=20$ and $n=200$ without loss of
generality.

We start presenting the results below. First, we draw some plots for the case
with $d=2$. In both Figures 2 and 3, the first sub-plot is always the true
$g(\mathbf{x})$. For the rest of sub-plots, each has three layers. The middle
one is the average of estimates over $n$ replications. The top and bottom
layers are the averages of the bootstrap draws corresponding to the 97.5% and
2.5% quantiles respectively over $n$ replications. A few facts emerge.
Overall, the LNN approach can recover the unknown function reasonably well.
When $q=4$, both figures have smoother plots, which should be expected
provided $q$ is considered as fixed. Also, both figures are very similar, so
the results are not sensitive to the choice of $u_{\sigma}$ as explained in
Remark 3.1.

More detailed numbers are summarized Table 1. As expected, when $T$ goes up,
RMSEg converges to 0 and CRg converges to 0.95. In particular, for CRg, when
the sample size is not large enough, we tend to get a narrower CI rather than
wider CI. Also, we note that when $d$ increases, RMSEg increases rather
significantly, which is understood as the curse of dimensionality. The curse
is also reflected in CRg, but the results are still acceptable in our view.
Recall that Figures 2 and 3 suggest that when $q=4$ the plots are much
smoothers, but Table 1 actually says larger $q$ is corresponding to larger
RMSEg. However, the change in $q$ does not alter CRg much, which shows certain
robustness of the bootstrap procedure suggested in Section 3.2. Finally, we
note that the time consuming part of this simulation design comes from
generating inference for a large number of points defined in (4.3) using the
proposed bootstrap procedure.

Figure 2: Simulation Results of Example 1.1 ($u_{\sigma}=-0.5,d=2$) Figure 3: Simulation Results of Example 1.1 ($u_{\sigma}=0.5,d=2$) Table 1: Simulation Results of Example 1.1 |  |  | RMSEg |  | CRg
---|---|---|---|---|---
|  | $T\setminus d$ | 2 | 3 |  | 2 | 3
$u_{\sigma}=-0.5$ | $q=3$ | 800 | 0.282 | 0.637 |  | 0.894 | 0.871
|  | 1600 | 0.261 | 0.543 |  | 0.925 | 0.884
|  | 2400 | 0.211 | 0.356 |  | 0.929 | 0.916
| $q=4$ | 800 | 0.324 | 0.932 |  | 0.917 | 0.837
|  | 1600 | 0.222 | 0.529 |  | 0.927 | 0.896
|  | 2400 | 0.180 | 0.408 |  | 0.934 | 0.916
$u_{\sigma}=0.5$ | $q=3$ | 800 | 0.291 | 0.608 |  | 0.912 | 0.873
|  | 1600 | 0.261 | 0.388 |  | 0.917 | 0.907
|  | 2400 | 0.236 | 0.293 |  | 0.918 | 0.910
| $q=4$ | 800 | 0.327 | 0.126 |  | 0.915 | 0.837
|  | 1600 | 0.226 | 0.579 |  | 0.922 | 0.891
|  | 2400 | 0.184 | 0.433 |  | 0.927 | 0.911

###  4.2 Model (1.4)

Next, we consider the following data generating process:

$\displaystyle
y_{t}=\left\\{\begin{array}[]{ll}1&g(\mathbf{x}_{t})-\varepsilon_{t}\geq 0\\\
0&\text{otherwise}\end{array}\right.,$

in which $\\{\mathbf{x}_{t}\\}$, $\\{\varepsilon_{t}\\}$ and $g(\cdot)$ are
exactly the same as those in Section 4.1.

We first note a computational issue. We reply on “fminunc” function of Matlab
to find the solution of

$\displaystyle\operatorname*{\arg\\!\min}_{\widetilde{s}\in\mathcal{S}}[-\log
L(\widetilde{s}(\cdot\,|\,\boldsymbol{\Theta}))],$

which is the same as that in (3.11). In order to invoke the minimization
process in any statistical software (including R, Matlab, etc.), one needs to
provide initial values to the parameters under estimation. As a consequence,
the numbers reported below are affected by the initial values more or less.
Although it is not our intention to tackle this complicated computational
issue in this paper, Lemma 3.4 does become useful in this case. Recall that
Lemma 3.4 bridges the log likelihood estimation and the nonlinear least
squares estimation. Therefore, we first conduct an OLS estimation using the
approach of Section 4.1 as the initial value of $\boldsymbol{\Theta}$ for each
generated $\\{(y_{t},\mathbf{x}_{t})\,|\,t\in[T]\\}$. We then invoke log
likelihood estimation as our final estimate of $\boldsymbol{\Theta}$ for each
dataset. Even in this case, the computation is rather slow, and the
computational time increases dramatically when the number of parameters goes
up.

That said, due to the limitation of computing power, we consider $d=2$ only in
this section. The rest settings are the same as those in Section 4.1. On top
of the criteria $\text{RMSE}_{g}$ and $\text{CR}_{g}$ introduced above, we
further introduce

$\displaystyle\text{RMSE}_{g}^{*}=\left\\{\frac{1}{\sharp\text{Q}_{L,\mathbf{j}}\cdot
L^{d}}\sum_{\mathbf{j}\in[L]^{d}}\sum_{\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})\in\text{Q}_{L,\mathbf{j}}}[\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})-g(\mathbf{x}_{L,\mathbf{j}})]^{2}\right\\}^{1/2},$

where $\text{Q}_{L,\mathbf{j}}$ defines the values between 1st and 3rd
quartiles of $\\{\widehat{g}_{i}(\mathbf{x}_{L,\mathbf{j}})\,|\,i\in[n]\\}$.
By doing so, we can remove some impacts of the initial value issues associated
with “fminunc”. A similar treatment can also be found in the simulation study
of Bauer and Kohler (2019), in which some quantiles of the simulation results
are reported in order to eliminate the impacts of outliers.

Again, we draw a few plots as in Figures 2 and 3. In this case, to minimize
the impacts of initial value issues caused by the program, we take the median
based on the simulation results when plotting each layer. As shown in Figures
4 and 5 below, the results are clearly not as good as those for Model 1.1 due
to lack of closed-form solutions. Still, we are able to see the consistency as
the sample size increases. Both figures show that the patterns are very much
similar to those in Figures 2 and 3, so we omit the discussion here.

We further summarize the detailed numbers in Table 2. A few facts should be
mentioned. First, the coverage rates are reasonably well. The values of
RMSE${}_{g}^{*}$ are comparable to those reported in Table 1, in which the
numbers are generated from closed-form estimates. Third, the results are not
changing much with respect to the value of $u_{\sigma}$, which again confirms
our argument in Remark 3.1.

Figure 4: Simulation Results of Example 1.4 ($u_{\sigma}=-0.5,d=2$) Figure 5: Simulation Results of Example 1.4 ($u_{\sigma}=0.5,d=2$) Table 2: Simulation Results of Example 1.4 |  | $T$ | RMSEg | RMSE${}_{g}^{*}$ | CRg
---|---|---|---|---|---
$u_{\sigma}=-0.5$ | $q=3$ | 800 | 1.212 | 0.790 | 0.904
|  | 1600 | 1.042 | 0.537 | 0.912
|  | 2400 | 0.967 | 0.428 | 0.917
| $q=4$ | 800 | 1.378 | 0.813 | 0.915
|  | 1600 | 0.797 | 0.405 | 0.923
|  | 2400 | 0.525 | 0.282 | 0.934
$u_{\sigma}=0.5$ | $q=3$ | 800 | 1.106 | 0.600 | 0.898
|  | 1600 | 1.023 | 0.577 | 0.905
|  | 2400 | 0.795 | 0.449 | 0.907
| $q=4$ | 800 | 1.110 | 0.616 | 0.906
|  | 1600 | 0.687 | 0.384 | 0.919
|  | 2400 | 0.480 | 0.278 | 0.930

##  5 Empirical Study

In this section, we consider two empirical examples to illustrate the two
models that have been investigated in Section 3.

###  5.1 On Climate Data

First, we would like to recall that as mentioned in Section 2, we assume that
$\\{\mathbf{x}_{t}\,|\,t\in[T]\\}$ are strictly stationary for the purpose of
simplicity. In practice, it is not surprising that data may present certain
time trend or seasonality as we shall show below. Such features do not alter
our theoretical results fundamentally. For example, we can suppose that

$\displaystyle\mathbf{x}_{t}=\mathbf{u}(\tau_{t})+\mathbf{s}_{t}+\mathbf{v}_{t},$
(5.1)

where $\tau_{t}=\frac{t}{T}$, $\mathbf{u}(\cdot)$ is a vector of deterministic
trending functions, $\mathbf{s}_{t}$ captures some seasonal effects, and
$\mathbf{v}_{t}$ mimics randomness. With some obvious modifications on
Assumption 2, all the results established previously still hold.

The data to be investigated are collected from Weather Undergroud API, and has
also been extensively studied by researchers from a very diversified
background on Kaggle333See https://www.kaggle.com/datasets/sumanthvrao/daily-
climate-time-series-data?resource=download for details.. This dataset includes
four time series (temperature, humidity, wind speed, atmospheric pressure)
from $1^{st}$ January 2013 to $24^{th}$ April 2017 for the city of Delhi,
India. As a routine exercise, we normalize each time series to get sample mean
0 and sample standard deviation 1 respectively. Typically, researchers use the
time period from $1^{st}$ January 2017 to $24^{th}$ April 2017 as the test
data, and use the rest time period as the training data. We do the same, so we
have $T=1461$ observations to train our LNN, and the size of test data is
$T^{*}=114$.

We plot the normalized time series in Figure 6 below. In each sub-plot, the
observations on the left hand of the vertical line are the training data,
while the observations on the right hand side are the test data. Temperature
and pressure move in opposite direction, while humidity and wind speed look
like white noises more or less. Although wind speed has some spikes from time
to time, it is rather stationary overall. When training our LNN, we consider
Model 1.1 and fit temperature at period $t$ (i.e., $y_{t}$) on humidity, wind
speed and atmospheric pressure at period $t-1$ (i.e., $\mathbf{x}_{t}$).

Figure 6: Climate Data of Delhi

In view of the data range, we specifically consider the following settings
when running the regression:

$\displaystyle a=2,\quad q\in\\{3,4,5,6\\},\quad u_{\sigma}\in\\{-0.5,0.5\\}$
(5.2)

The rest settings are the same as those in simulation studies. First, we
report RMSE using test set. Specifically, after training our LNN, we calculate

$\displaystyle\text{RMSE}=\left\\{\frac{1}{T^{*}}\sum_{t=1}^{T^{*}}[y_{t}-\widehat{g}(\mathbf{x}_{t})]^{2}\right\\}^{1/2},$
(5.3)

where $(y_{t},\mathbf{x}_{t})$ are from the test set, and $\widehat{g}(\cdot)$
is obtained using the training set.

We summarize RMSE in Table 3 below, in which the results are very close.
Regardless the value of $u_{\sigma}$, when we increase the value of $q$, RMSE
always shows a U shape, which can be used as a criterion to pick the “optimal”
$q$ in practice. In what follows, we focus on the cases
$(u_{\sigma},q)=(0.5,5)$ and $(u_{\sigma},q)=(-0.5,4)$, as both combinations
are corresponding to the bottom points of each U–curve.

In both Figures 7 and 8, the black dotted line stands for the true data, the
red dashdotted line stands for the estimated values, and the blue solid lines
present the 95% CI based on the bootstrap draws. Figure 7 shows that the LNN
approach can capture the trend of the entire test period reasonably well over
all, and Figure 8 provides a zoomed-in version by focusing on the test period
only. It is clear that both subplots are almost identical, although the values
of $(u_{\sigma},q)$ are different. Therefore, we conclude that the numerical
results show certain robustness of NN architecture.

Table 3: RMSE for the Test Set $u_{\sigma}\setminus q$ | 3 | 4 | 5 | 6
---|---|---|---|---
0 | 0.5447 | 0.5227 | 0.5287 | 0.5354
0.5 | 0.5719 | 0.5267 | 0.5238 | 0.5368
Figure 7: Plot of the Fitted Curve on the Entire Period Figure 8: Plot of the
Fitted Curve of the Test Period

###  5.2 On Clean Energy Index

As pointed by Christoffersen and Diebold (2006) and many follow-up studies,
the sign of stock market returns may be predictable even if the returns
themselves are not predictable. In this section, we specifically use Model 1.4
to investigate the clean energy index (ECO) assembled by WiderShares, LLC,
which is a modified equal weighted index comprised of companies that are
publicly traded in the United States and that are engaged in the business of
the advancement of cleaner energy and conservation. A very detailed data
description about ECO can be found at https://wildershares.com/about.php, so
we omit them here for simplicity. In connection with the recent trend of
investigating climate change, the importance of understanding clean energy
sector is more important than ever. Not to mention that for finance
discipline, studying energy sector is always a topical interest (e.g., Jin and
Jorion, 1996; El-Sharif et al., 2005; and many follow-up studies).

In what follows, we focus on the next model

$\displaystyle
y_{t}=\left\\{\begin{array}[]{ll}1&g(\mathbf{x}_{t})-\varepsilon_{t}\geq 0\\\
0&\text{otherwise}\end{array}\right.,$

where $y_{t}=1$ if ECO index return is positive at time $t$, otherwise
$y_{t}=0$; and $\mathbf{x}_{t}$ includes oil price change rate at time $t-1$,
gas price change rate at time $t-1$, and CBOE volatility index (VIX) change
rate at time $t-1$. The choices of regressors are based on Section III of Jin
and Jorion (1996) and Section 7 of Christoffersen and Diebold (2006). One
definitely can consider more choices and different forms of regressors, which
may lead to another research paper with comprehensive empirical study. Again,
we normalize the observations of each regressor to ensure mean 0 and standard
deviation 1. Finally, we present the data in Figure 9.

Figure 9: Clean Energy Index Data

The entire time period we consider covers 4th January 2010 to 31st December
2019, in which we use the period from 1st July 2019 to 31st December 2019 as
the test period. As a result, we have 2369 observations as training data
(i.e., $T=2369$), and 146 observations as test data (i.e., $T^{*}=146$).

After training our LNN, we calculate the mean absolute error (MAE) as follows:

$\displaystyle\text{MAE}=\frac{1}{T^{*}}\sum_{t=1}^{T^{*}}|y_{t}-\widehat{y}_{t}|,$

where $\widehat{y}_{t}=1$ if $\widehat{g}(\mathbf{x}_{t})\geq 0.5$, and
$\widehat{y}_{t}=0$ otherwise. Here, we arbitrarily use 0.5 as the threshold.
One may further explore other options in practice. We consider different
choices of $u_{\sigma}$ and $q$, and summarize results in Table 4 below.
Regardless the value of $u_{\sigma}$, we obtain the minimum MAE with $q=4$. In
the best case scenario, the successful rate is around 59% which might be
further improved by accounting for volatility as well as other possible
explanatory variables through a comprehensive empirical investigation.

Table 4: MAE for the Test Set $u_{\sigma}\setminus q$ | 3 | 4 | 5
---|---|---|---
-0.5 | 0.5068 | 0.4384 | 0.5137
0.5 | 0.5205 | 0.4110 | 0.4589

In what follows, we focus on the results by setting $q=4$. In Figure 10, using
test data we respectively plot the estimated $\widehat{g}(\cdot)$, and the
corresponding probabilities (i.e., $\Phi_{\varepsilon}(\widehat{g}(\cdot))$).
In each sub-figure, the dash-dotted line is corresponding to the choice of
$(u_{\sigma}=-0.5,q=4)$, while the dotted line is corresponding to the choice
of $(u_{\sigma}=0.5,q=4)$. The results are almost identical when $u_{\sigma}$
varies, so it again justifies our argument on the choice of $u_{\sigma}$ and
the robustness of the approach investigated above.

Figure 10: Plot of $\widehat{g}(\mathbf{x})$ and associated Probabilities of
the Test Period

Finally, we point out that Christoffersen and Diebold (2006) write “As
volatility moves, so too does the probability of a positive return: the higher
the volatility, the lower the probability of a positive return”. With good
sign prediction, we can then remove the “bad volatility” and keep the “good
volatility” (e.g., Patton and Sheppard, 2015). This is important, as
traditionally a primary goal of portfolio analysis is to minimize the
volatility (e.g., Engle et al., 2019), so many good volatilities are wiped out
by classic approaches. Model (1.4) associated with the LNN approach naturally
marries the above studies and improves the prediction of positive returns, so
a better portfolio strategy can be implemented in practice. We will no longer
proceed further along this line, and would like to leave this in another
research paper by systematically comparing different approaches adopted when
constructing portfolio.

##  6 Conclusion

NN has gained considerable attentions over the past a few decades. Yet,
questions as those stated in Section 1 have not been answered well in the
literature. In this paper, we bring in identification restrictions to the LNN
framework from a nonparametric regression perspective, and consider the LNN
based estimation with dependent data for both cases with
quantitative/qualitative outcomes. We then establish the LNN based estimation
theory under a set of minor conditions. The asymptotic distributions are
derived accordingly, and we show that LNN automatically eliminates the
dependence of data when calculating the asymptotic variances. The finding is
important, as one can easily use different types of wild bootstrap methods to
obtain valid inference for practical implementation. In particular, for
quantitative outcomes, the LNN approach yields closed–form expressions for the
estimates of some key parameters of interest. Last but not least, we examine
our theoretical findings through extensive numerical studies.

Several major comments have been made here and there in Section 3, and some
future research directions have been acknowledged along the way. Finally, we
hope the current article will shed light on how to produce transparent
algorithms to ensure that our research findings are relevant and useful for
practical implementations and applications.

##  7 Acknowledgements

Gao, Peng and Yang acknowledge financial support from the Australian Research
Council Discovery Grants Program under Grant Numbers: DP200102769, DP210100476
and DP230102250, respectively.

## References

  * (1)
  * Abramovitz and Stegun (1972) Abramovitz, M. and Stegun, I. A. (1972), Handbook of mathematical functions, Dover Publications, New York, U.S.
  * Athey (2019) Athey, S. (2019), The impact of machine learning on economics, in J. G. Ajay Agrawal and A. Goldfarb, eds, ‘The Economics of Artificial Intelligence: An Agenda’, pp. 507–547.
  * Bartlett et al. (2021) Bartlett, P. L., Montanari, A. and Rakhlin, A. (2021), ‘Deep learning: A statistical viewpoint’, Acta Numerica 30, 87201.
  * Bauer and Kohler (2019) Bauer, B. and Kohler, M. (2019), ‘On deep learning as a remedy for the curse of dimensionality in nonparametric regression’, The Annals of Statistics 47(4), 2261–2285.
  * Chen (2007) Chen, X. (2007), ‘Large Sample Sieve Estimation of Semi–Nonparametric Models, Chapter 76 edited by James J. Heckman and Edward E. Leamer’, Handbook of Econometrics 6B, 5549–5632.
  * Chen et al. (2022) Chen, X., Liu, Y., Ma, S. and Zhang, Z. (2022), Casual inference of general treatment effects using neural networks with a diverging number of confounders.  Available at https://arxiv.org/abs/2009.07055v5.
  * Chen et al. (2001) Chen, X., Racine, J. and Swanson, N. (2001), ‘Semiparametric ARX neural network models with an application to forecasting inflation’, IEEE Transactions on Neural Networks 12(6), 674–683.
  * Chen and Shen (1998) Chen, X. and Shen, X. (1998), ‘Sieve extremum estimates for weakly dependent data’, Econometrica 66(2), 298–314.
  * Chen and White (1999) Chen, X. and White, H. (1999), ‘Improved rates and asymptotic normality for nonparametric neural network estimators’, IEEE Transactions on Information Theory 45(6), 682–691.
  * Christoffersen and Diebold (2006) Christoffersen, P. F. and Diebold, F. X. (2006), ‘Financial asset returns, direction-of-change forecasting, and volatility dynamics’, Management Science 52(8), 1273–1287.
  * Cybenko (1989) Cybenko, G. (1989), ‘Approximation by superpositions of a sigmoidal function’, Mathematics of Control, Signals and Systems 2, 303–314.
  * Du et al. (2021) Du, X., Fan, Y., Lv, J., Sun, T. and Vossler, P. (2021), Dimension-free average treatment effect inference with deep neural networks.  Available at https://doi.org/10.48550/arXiv.2112.01574.
  * Dubey et al. (2022) Dubey, S. R., Singh, S. K. and Chaudhuri, B. B. (2022), ‘Activation functions in deep learning: A comprehensive survey and benchmark’, Neurocomputing 503, 92–108.
  * El-Sharif et al. (2005) El-Sharif, I., Brown, D., Burton, B., Nixon, B. and Russell, A. (2005), ‘Evidence on the nature and extent of the relationship between oil prices and equity values in the UK’, Energy Economics 27(6), 819–830.
  * Engle et al. (2019) Engle, R. F., Ledoit, O. and Wolf, M. (2019), ‘Large dynamic covariance matrices’, Journal of Business & Economic Statistics 37(2), 363–375.
  * Fan and Gijbels (1996) Fan, J. and Gijbels, I. (1996), Local Polynomial Modelling and its Applications, Chapman & Hall/CRC.
  * Fan and Gu (2022) Fan, J. and Gu, Y. (2022), Factor augmented sparse throughput deep relu neural networks for high dimensional regression.  Available at https://doi.org/10.48550/arXiv.2210.02002.
  * Fan et al. (2021) Fan, J., Ma, C. and Zhong, Y. (2021), ‘A selective overview of deep learning’, Statistical Science 36(2), 264–290.
  * Fan and Yao (2003) Fan, J. and Yao, Q. (2003), Nonlinear Time Series: Nonparametric and Parametric Methods, Springer-Verlag.
  * Farrell et al. (2021) Farrell, M. H., Liang, T. and Misra, S. (2021), ‘Deep neural networks for estimation and inference’, Econometrica 89(1), 181–213.
  * Gao (2007) Gao, J. (2007), Nonlinear Time Series: Semiparametric and Nonparametric Methods, Chapman and Hall.
  * Gu et al. (2020) Gu, S., Kelly, B. and Xiu, D. (2020), ‘Empirical asset pricing via machine learning’, The Review of Financial Studies 33(5), 2223–2273.
  * Gu et al. (2021) Gu, S., Kelly, B. and Xiu, D. (2021), ‘Autoencoder asset pricing models’, Journal of Econometrics 222(1, Part B), 429–450.
  * Günther and Fritsch (2010) Günther, F. and Fritsch, S. (2010), ‘Neuralnet: Training of neural networks’, R Journal 2, 30–38.
  * Hill et al. (1996) Hill, T., O’Connor, M. and Remus, W. (1996), ‘Neural network models for time series forecasts’, Management Science 42(7), 1082–1092.
  * Huang et al. (2008) Huang, J., Horowitz, J. L. and Ma, S. (2008), ‘Asymptotic properties of bridge estimators in sparse high-dimensional regression models’, The Annals of Statistics 36(2), 587–613.
  * Jin and Jorion (1996) Jin, Y. and Jorion, P. (1996), ‘Firm value and hedging: Evidence from U.S. oil and gas producers’, The Journal of Finance 61(2), 893–919.
  * Kohler and Krzyzák (2017) Kohler, M. and Krzyzák, A. (2017), ‘Nonparametric regression based on hierarchical interaction models’, IEEE Transactions on Information Theory 63(3), 1620–1630.
  * Li et al. (2016) Li, D., Tjøstheim, D. and Gao, J. (2016), ‘Estimation in nonlinear regression with Harris recurrent Markov chains’, The Annals of Statistics 44(5), 1957–1987.
  * Li and Racine (2007) Li, Q. and Racine, J. (2007), Nonparametric Econometrics Theory and Practice, Princeton University Press, New Jersey.
  * Magnus and Neudecker (2007) Magnus, J. R. and Neudecker, H. (2007), Matrix Differential Calculus with Applications in Statistics and Econometrics, third edn, John Wiley & Sons Ltd.
  * Minai and Williams (1993) Minai, A. A. and Williams, R. D. (1993), ‘On the derivatives of the sigmoid’, Neural Networks 6(6), 845–853.
  * Mudelsee (2019) Mudelsee, M. (2019), ‘Trend analysis of climate time series: A review of methods’, Earth-Science Reviews 190, 310–322.
  * Murata et al. (1994) Murata, N., Yoshizawa, S. and Amari, S. (1994), ‘Network information criterion-determining the number of hidden units for an artificial neural network model’, IEEE Transactions on Neural Networks 5(6), 865–872.
  * Newey and McFadden (1994) Newey, W. K. and McFadden, D. (1994), Chapter 36 large sample estimation and hypothesis testing, Vol. 4 of Handbook of Econometrics, Elsevier, pp. 2111–2245.
  * Newey and Powell (2003) Newey, W. K. and Powell, J. L. (2003), ‘Instrumental variable estimation of nonparametric models’, Econometrica 71(5), 1565–1578.
  * Patrick and Andres (2012) Patrick, K. and Andres, S. (2012), ‘A score based approach to wild bootstrap inference’, Journal of Econometric Methods 1(1), 23–41.
  * Patton and Sheppard (2015) Patton, A. J. and Sheppard, K. (2015), ‘Good volatility, bad volatility: Signed jumps and the persistence of volatility’, The Review of Economics and Statistics 97(3), 683–697.
  * Sauer (2006) Sauer, T. (2006), Polynomial interpolation in several variables: Lattices, differences, and ideals, in K. Jetter, M. D. Buhmann, W. Haussmann, R. Schaback and J. Stöckler, eds, ‘Topics in Multivariate Approximation and Interpolation’, Vol. 12 of Studies in Computational Mathematics, pp. 191–230.
  * Schmidt-Hieber (2020) Schmidt-Hieber, J. (2020), ‘Nonparametric regression using deep neural networks with ReLU activation function’, The Annals of Statistics 48(4), 1875–1897.
  * Tibshirani (1996) Tibshirani, R. (1996), ‘Regression shrinkage and selection via the Lasso’, Journal of the Royal Statistical Society. Series B (Methodological) 58(1), 267–288.
  * Wang and Lin (2021) Wang, H. and Lin, W. (2021), Harmless overparametrization in two-layer neural networks.  Available at https://doi.org/10.48550/arXiv.2106.04795.
  * Zhong et al. (2022) Zhong, Q., Mueller, J. and Wang, J.-L. (2022), ‘Deep learning for the partially linear Cox model’, The Annals of Statistics 50(3), 1348–1375.

Online Supplementary Appendix to
“A Localized Neural Network with Dependent Data:
Estimation and Inference”

Jiti Gao†, Bin Peng† and Yanrong Yang∗

†Monash University and ∗Australian National University

In this file, we first provide extra information about Sigmoidal squasher in
Appendix A.1. Appendix A.2 provides some extra simulation results. In Appendix
A.3, we explain how to relax the restriction about the compact set to allow
for $\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$. We then introduce some notations
and preliminary lemmas in Appendix A.4; their proofs are given in Appendix
A.5. It is noted that Lemmas A.3 and A.4 below are of general and independent
interest, as they establish useful asymptotic properties for the
point–estimators of the unknown parameters of interest.

Appendix A

A.1 Derivatives of the Sigmoidal Squasher

As we adopt Sigmoidal squasher (i.e., $\sigma(x)=\frac{1}{1+\exp(-x)}$)
throughout the numerical studies, we provide some of its properties below,
which are useful for numerical implementation. By Minai and Williams (1993),
we have

$\displaystyle\sigma^{(n)}(x)=\sum_{k=1}^{n+1}(-1)^{k+1}(k-1)!S(n+1,k)\cdot\sigma(x)^{k},$

in which

$\displaystyle S(n+2,k)=kS(n+1,k)+S(n+1,k-1).$ (A.1)

is the recursion relation for Stirling numbers of the second kind.

Sigmoidal squasher is easy to use in the sense that we can arbitrarily choose
$u_{\sigma}$ of Lemma 3.1. Without loss of generality, we let
$u_{\sigma}\in\\{-0.5,0.5\\}$ in the numerical studies. In Figure A.1, we plot
$\sigma^{(k)}(x)$ for $k=0,\ldots,5$ for the purpose of demonstration:

Figure A.1: Plots of Sigmoidal Squasher with Its Derivatives

A.2 Extra Simulation Results

In this section, we increase the value of $d$ to see whether the LNN based
estimation method still has good finite sample performance. The following
exercises are in line with the same spirit of those provided in the
supplementary file of Farrell et al. (2021), where they show, even for the
i.i.d. data, one needs around 10,000 observations to generate some meaningful
results with large $d$. We focus on Model (1.1) only. Note that when $d$
increases, we can no longer evaluate the pre-specified points by following
(4.3), as the computational overhead increases exponentially. In fact, even
for the case $d=5$, we cannot finish the computation within a few weeks using
the cluster facility of Monash University if we stick to (4.3). Therefore, we
consider the following points:

$\displaystyle\mathbf{x}_{L,i}=-2.5+(i-1)\cdot
0.2\mathbf{1}_{d}\quad\text{with}\quad i=1,\ldots,26.$ (A.2)

As explained in Section 4, our focus is still the coverage rate. Also, for the
purpose of comparison, we also consider the kernel regression with Uniform and
Epanechnikov kernels respectively, and the confidence intervals are
constructed using a bootstrap procedure similar to that in Section 3.2, so we
omit the details. The rest settings are the same as in Section 4.

As shown in Table A.1, the LNN approach shows much better performance when
$d=5$, and the coverage rates move towards to 95% much faster. For the kernel
methods, the coverage rates increase extremely slow for Epanechnikov kernel,
and does not seem to increase at all for Uniform kernel. Note that the points
selected in (A.2) are not on the boundary, so the results should not suffer
boundary effects. It is not our intention to criticise the traditional kernel
method, but LNN does show some superiority in terms of finite sample study.

Due to limits of our computational power, we no longer further increase the
value of $d$, and also acknowledge the fact our code may not be as efficient
as in some well established statistical packages.

Table A.1: Coverage rates ($d=5$) $T$ |  NN ($u_{\sigma}=-0.5$) |  NN ($u_{\sigma}=0.5$) | Epanechnikov | Uniform
---|---|---|---|---
2400 | 0.658 | 0.672 | 0.641 | 0.726
4800 | 0.861 | 0.848 | 0.682 | 0.722
9600 | 0.899 | 0.889 | 0.707 | 0.718

A.3 From $[-a,a]^{d}$ to $\mathbb{R}^{d}$

Consider Model (1.1) only for simplicity. We provide two treatments to relax
the restriction on $a$.

Treatment 1: Note that in Section 3, we require $g(\mathbf{x})$ to be defined
on a compact set, but do not impose restriction on the range of
$\\{\mathbf{x}_{t}\\}$. In fact, for time series data, it may make more sense
to assume that $a$ is diverging, which is indeed achievable. Suppose that
$\mathbf{x}_{t}$ follows a sub-Gaussian distribution. In this case, we still
need to require $\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{i}}^{-1}\to\infty$ in
view of Lemma A.3 below, which after some algebra is equivalent to asking for

$\displaystyle\sqrt{\log(Th^{d})}\cdot a\to\infty.$

Apparently, there is a price that we have to pay, which is the slow rate of
convergence. A similar treatment has also been discussed in Li et al. (2016)
for example, so we do not further elaborate it here.

Treatment 2: Alternatively, we can modify the construction of LNN from a
nonparametric viewpoint. In the literature of kernel regression, one normally
pre-specifies a point of interest (e.g., $\mathbf{x}_{0}$), and investigates a
small area nearby only which is usually decided by some bandwidth(s)
converging to 0. As a result, the parameters obtained from the estimation
procedure usually vary with respect to the point of interest. In other words,
when evaluating different points from the test set, the number of parameters
to be estimated will be proportional to the cardinality of the training set
(although estimation is always carried on using the same training set).
Provided a large test set, it may create lots of overhead from a computational
viewpoint, but the advantage is that in theory it allows us to consider
$\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$.

That said, consider Model (1.1) as an example. For
$\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$, we consider the following objective
function:

$Q_{L}(\boldsymbol{\theta})=\sum_{t=1}^{T}[y_{t}-s(\mathbf{x}_{t}\,|\,\mathbf{x}_{0},\boldsymbol{\theta})]^{2}\cdot
I_{0,h}(\mathbf{x}_{t}),$ (A.3)

where the subscript $L$ infers the local version, $\boldsymbol{\theta}$ is
$d_{q}\times 1$ vector satisfying $\|\boldsymbol{\theta}\|<\infty$, we let
$I_{0,h}(\mathbf{x}_{t}):=I(\mathbf{x}_{t}\in C_{\mathbf{x}_{0},h})$ for
short, and $C_{\mathbf{x}_{0},h}$ is defined in (3.3) already. Then the OLS
estimate of $\widetilde{\boldsymbol{\lambda}}$ defined in Lemma 3.2 is
obtained by

$\displaystyle\widehat{\boldsymbol{\theta}}=\operatorname*{\arg\\!\min}_{\boldsymbol{\theta}}Q_{L}(\boldsymbol{\theta})$

and, accordingly, the estimate of $g(\mathbf{x}_{0})$ is defined by

$\displaystyle\widehat{g}_{L}(\mathbf{x}_{0})=s(\mathbf{x}_{0}\,|\,\mathbf{x}_{0},\widehat{\boldsymbol{\theta}}).$

Note that

$\displaystyle\widehat{g}_{L}(\mathbf{x}_{0})-g(\mathbf{x}_{0})$
$\displaystyle\simeq$ $\displaystyle
s(\mathbf{x}_{0}\,|\,\mathbf{x}_{0},\widehat{\boldsymbol{\theta}})-s(\mathbf{x}_{0}\,|\,\mathbf{x}_{0},\widetilde{\boldsymbol{\lambda}})$
(A.4) $\displaystyle\simeq$
$\displaystyle\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})^{\top}\mathbf{H}\mathbf{H}^{-1}\mathbf{D}^{\top,-1}(\widehat{\boldsymbol{\theta}}-\widetilde{\boldsymbol{\lambda}}),$

where the first step is due to Lemma A.1, and the last step follows from
(A.15) of the online supplementary appendix. Here, it should be understood
that $\widetilde{\boldsymbol{\lambda}}$ is decided by $\mathbf{x}_{0}$ only.

We summarize the results in the following corollary.

######  Corollary A.1.

Suppose that Assumptions 1 and 2 hold. As $(1/h,Th^{d})\to(\infty,\infty)$,
for $\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$,

  1. 1.

$\sqrt{Th^{d}}\,\mathbf{H}^{-1}\left(\mathbf{D}^{\top}\right)^{-1}(\widehat{\boldsymbol{\theta}}-\widetilde{\boldsymbol{\lambda}}+O_{P}(h^{p}))\to_{D}N\left(\mathbf{0},\sigma_{\varepsilon}^{2}\boldsymbol{\Sigma}_{\mathbf{x}_{0}}^{-1}\right),$

  2. 2.

$\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{x}_{0}}^{-1}(\widehat{g}_{L}(\mathbf{x}_{0})-g(\mathbf{x}_{0})+O_{P}(h^{p}))\to_{D}N(0,1),$

where $\widetilde{\boldsymbol{\lambda}}$ is uniquely determined by
$g(\mathbf{x}_{0})$.

The first result of Corollary A.1 in fact provides CLT for the LNN parameters.
While some CLTs are established previously, the focus has been on the unknown
function without much inference on the effective parameters (e.g., Gu et al.,
2021; Farrell et al., 2021; Chen et al., 2022).

A.4 Preliminary Lemmas

For $\mathbf{a}=(a_{0},a_{1},\ldots,a_{d})^{\top}\in\mathbb{R}^{d+1}$,
$\mathbf{r}=(r_{0},r_{1},\ldots,r_{d})^{\top}\in\mathbb{N}_{0}^{d+1}$, and
$\mathbf{x}\in\mathbb{R}^{d}$, we let

$\displaystyle
f_{\mathbf{a}}(\mathbf{x})=[(1,\mathbf{x}^{\top})\mathbf{a}]^{q}=\sum_{|\mathbf{r}|=q}\binom{q}{\mathbf{r}}\cdot
a_{0}^{r_{0}}\prod_{k=1}^{d}a_{k}^{r_{k}}x_{k}^{r_{k}}.$

Obviously, we have $f_{\mathbf{a}}\in\mathscr{P}_{q}$, where $\mathscr{P}_{q}$
is defined in (3.1).

######  Lemma A.1.

Let $f\,:\,\mathbb{R}^{d}\to\mathbb{R}$ be a $(p,\mathscr{C})$-smooth
function. For $\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$, let

$\displaystyle
p_{q}(\mathbf{x}\,|\,\mathbf{x}_{0})=\sum_{0\leq|\mathbf{J}|\leq
q}\frac{1}{\mathbf{J}!}\cdot\frac{\partial^{|\mathbf{J}|}f(\mathbf{x}_{0})}{\partial\mathbf{x}^{\mathbf{J}}}(\mathbf{x}-\mathbf{x}_{0})^{\mathbf{J}}.$

Then

$\displaystyle\|f(\mathbf{x})-p_{q}(\mathbf{x}\,|\,\mathbf{x}_{0})\|_{\infty}\leq
O(1)\|\mathbf{x}-\mathbf{x}_{0}\|^{p},$

where $p=q+s$, $O(1)$ depends on $d$ and $q$ only.

######  Lemma A.2.

For almost all $\mathbf{a}_{1},\ldots,\mathbf{a}_{d_{q}}\in\mathbb{R}^{d+1}$
(with respect to the Lebesgue measure in $\mathbb{R}^{(d+1)\times d_{q}}$), we
have that
$\\{f_{\mathbf{a}_{1}}(\mathbf{x}),\ldots,f_{\mathbf{a}_{d_{q}}}(\mathbf{x})\\}$
is a basis of the linear vector space $\mathscr{P}_{q}$.

######  Lemma A.3 (On Model (1.1)).

Suppose Assumptions 1 and 2 hold. As $(1/h,Th^{d})\to(\infty,\infty)$, for
each $\mathbf{i}\in[M]^{d}$,

$\displaystyle\sigma_{\varepsilon}^{-1}\boldsymbol{\Sigma}_{\mathbf{i}}^{1/2}\sqrt{Th^{d}}\mathbf{H}^{-1}\mathbf{D}^{\top,-1}(\widehat{\boldsymbol{\theta}}_{\mathbf{i}}-\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}+O_{P}(h^{p}))\to_{D}N\left(\mathbf{0},\mathbf{I}\right),$

where
$\boldsymbol{\Sigma}_{\mathbf{i}}=f_{\mathbf{x}}(\widetilde{\mathbf{x}}_{\mathbf{i}})\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}$.

######  Lemma A.4 (On Model (1.4)).

Under Assumptions 1 and 2, for each $\mathbf{i}\in[M]^{d}$,

  1. 1.

$\frac{1}{T}\sum_{t=1}^{T}I_{\mathbf{i},h}(\mathbf{x}_{t})[g(\mathbf{x}_{t})-s(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widehat{\boldsymbol{\theta}}_{\mathbf{i}})]^{2}=o_{P}(1)$,

  2. 2.

$\left\|\frac{1}{T}\mathbf{H}\mathbf{D}\frac{\partial^{2}\log
L(\widetilde{s}(\cdot\,|\,\widehat{\boldsymbol{\Theta}}))}{\partial\boldsymbol{\theta}_{\mathbf{i}}\partial\boldsymbol{\theta}_{\mathbf{i}}^{\top}}\mathbf{D}^{\top}\mathbf{H}-\widetilde{\boldsymbol{\Sigma}}_{\mathbf{i}}\right\|=o_{P}(1)$,

  3. 3.

$\left\|\frac{1}{T}\mathbf{H}\mathbf{D}\frac{\partial^{2}\log
L(s(\cdot\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\theta}_{\mathbf{i}}^{*}))}{\partial\boldsymbol{\theta}_{\mathbf{i}}\partial\boldsymbol{\theta}_{\mathbf{i}}^{\top}}\mathbf{D}^{\top}\mathbf{H}-\widetilde{\boldsymbol{\Sigma}}_{\mathbf{i}}\right\|=o_{P}(1)$,

  4. 4.

$\widetilde{\boldsymbol{\Sigma}}_{\mathbf{i}}^{-1/2}\frac{1}{\sqrt{Th^{d}}}\mathbf{H}\mathbf{D}\frac{\partial\log
L(g(\cdot))}{\partial\boldsymbol{\theta}_{\mathbf{i}}}\to_{D}N(\mathbf{0},\mathbf{I})$,

where $\boldsymbol{\theta}_{\mathbf{i}}^{*}$ lies between
$\widehat{\boldsymbol{\theta}}_{\mathbf{i}}$ and
$\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$, and

$\displaystyle\widetilde{\boldsymbol{\Sigma}}_{\mathbf{i}}=\frac{f_{\mathbf{x}}(\widetilde{\mathbf{x}}_{\mathbf{i}})\phi_{\varepsilon}(g(\widetilde{\mathbf{x}}_{\mathbf{i}}))^{2}}{[1-\Phi_{\varepsilon}(g(\widetilde{\mathbf{x}}_{\mathbf{i}}))]\Phi_{\varepsilon}(g(\widetilde{\mathbf{x}}_{\mathbf{i}}))}\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}.$

A.5 Proofs

Proof of Lemma A.1:

This is Lemma 8 of Bauer and Kohler (2019), so the derivation is omitted.
$\blacksquare$

Proof of Lemma A.2:

In what follows, let

$\displaystyle\mathbf{a}_{k}=(a_{k,0},a_{k,1},\ldots,a_{k,d})^{\top}\
\text{for}\
k\in[d_{q}],\quad\text{and}\quad\widetilde{\mathbf{r}}=\\{\mathbf{r}\in\mathbb{N}_{0}^{d+1}\,|\,|\mathbf{r}|=q\\}.$

It suffices to show that
$f_{\mathbf{a}_{1}}(\mathbf{x}),\ldots,f_{\mathbf{a}_{d_{q}}}(\mathbf{x})$ are
linearly independent. To do this, let $b_{1},\ldots,b_{d_{q}}\in\mathbb{R}$ be
such that

$\displaystyle\sum_{k=1}^{d_{q}}b_{k}f_{\mathbf{a}_{k}}(\mathbf{x})=0.$ (A.5)

As we explained under (3.1), the monomials involved in (3.1) are linearly
independent. Thus, (A.5) implies that

$\displaystyle\sum_{k=1}^{d_{q}}b_{k}\prod_{j=0}^{d}a_{k,j}^{r_{j}}=0\ \text{
for }\ \forall\mathbf{r}\in\widetilde{\mathbf{r}}.$

Note that $\sharp\widetilde{\mathbf{r}}=d_{q}$ by design, so we can construct
a one-to-one relationship between $k\in[d_{q}]$ and
$\mathbf{r}\in\widetilde{\mathbf{r}}$. Using this relationship, we can
construct $p(\mathbf{x})\in\mathscr{P}_{q}$ as follows:

$\displaystyle
p(\mathbf{x})=\sum_{k=1}^{d_{q}}b_{k}\prod_{j=1}^{d}x_{j}^{r_{j}}=\sum_{\mathbf{r}\in\widetilde{\mathbf{r}}}b_{\mathbf{r}}\prod_{j=1}^{d}x_{j}^{r_{j}},$

which satisfies

$\displaystyle p(\mathbf{a}_{k})=0\quad\text{for}\quad\forall k\in[d_{q}].$
(A.6)

Position 4 in Sauer (2006) implies that (A.6) has the only solution
$p(\cdot)=0$ in $\mathscr{P}_{q}$ for Lebesgue almost all
$\mathbf{a}_{1},\ldots,\mathbf{a}_{d_{q}}\in\mathbb{R}^{d+1}$, which in turn
implies $b_{1}=\cdots=b_{d_{q}}=0$. The proof is now completed. $\blacksquare$

Proof of Lemma 3.1:

Before proceeding, we would like to point out that in what follows,
$\mathtt{C}$ is a constant and is for the purpose of rescaling only.

By Assumption 1.2, there is a point $u_{\sigma}\in\mathbb{R}$ such that none
of the derivatives up to the order $q$ is 0 at $u_{\sigma}$. Thus, we
construct the following one-layer NN:

$\displaystyle\sum_{k=1}^{q+1}(-1)^{q+k-1}\cdot\frac{\mathtt{C}^{q}}{\sigma^{(q)}(u_{\sigma})}\binom{q}{k-1}\cdot\sigma\left(\frac{k-1}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right)$
(A.7) $\displaystyle=$
$\displaystyle\sum_{k=0}^{q}(-1)^{q+k}\cdot\frac{\mathtt{C}^{q}}{\sigma^{(q)}(u_{\sigma})}\binom{q}{k}\cdot\sigma\left(\frac{k}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right)$
$\displaystyle=$
$\displaystyle(-1)^{q}\frac{\mathtt{C}^{q}}{\sigma^{(q)}(u_{\sigma})}\sum_{k=0}^{q}(-1)^{k}\cdot\binom{q}{k}\cdot\sigma\left(\frac{k}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right),$

in which the definitions of $\gamma_{k}$’s and $\beta_{k}$’s are obvious.

By Assumption 1.2 again, $\sigma(\cdot)$ is $q+1$ times continuously
differentiable. Thus, it can be expanded in a Taylor series with Lagrange
remainder around $u_{\sigma}$ up to order $q$:

$\displaystyle\sum_{k=0}^{q}(-1)^{k}\cdot\binom{q}{k}\cdot\sigma\left(\frac{k}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right)$
(A.8) $\displaystyle=$
$\displaystyle\sum_{k=0}^{q}(-1)^{k}\cdot\binom{q}{k}\cdot\left(\sum_{j=0}^{q}\frac{\sigma^{(j)}(u_{\sigma})\cdot((x-x_{0})k)^{j}}{\mathtt{C}^{j}\cdot
j!}+\frac{\sigma^{(q+1)}(\xi_{k})\cdot((x-x_{0})k)^{q+1}}{\mathtt{C}^{q+1}\cdot(q+1)!}\right)$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{q}\frac{\sigma^{(j)}(u_{\sigma})\cdot(x-x_{0})^{j}}{\mathtt{C}^{j}\cdot
j!}\sum_{k=0}^{q}(-1)^{k}\cdot k^{j}\cdot\binom{q}{k}$
$\displaystyle+\frac{(x-x_{0})^{q+1}}{\mathtt{C}^{q+1}\cdot(q+1)!}\sum_{k=0}^{q}(-1)^{k}\cdot
k^{q+1}\cdot\sigma^{(q+1)}(\xi_{k})\cdot\binom{q}{k},$

where
$\xi_{k}\in[u_{\sigma}-\frac{k}{\mathtt{C}}\cdot|x-x_{0}|,u_{\sigma}+\frac{k}{\mathtt{C}}\cdot|x-x_{0}|]$
for all $0\leq k\leq q$.

Note that

$\displaystyle\sum_{k=0}^{q}(-1)^{k}\cdot k^{j}\cdot\binom{q}{k}$
$\displaystyle=$ $\displaystyle
q!(-1)^{q}\cdot\frac{1}{q!}\sum_{k=0}^{q}(-1)^{k-q}\cdot(q-(q-k))^{j}\cdot\binom{q}{q-k}$
$\displaystyle=$ $\displaystyle
q!(-1)^{q}\cdot\frac{1}{q!}\sum_{k=0}^{q}(-1)^{q-k}\cdot(q-(q-k))^{j}\cdot\binom{q}{q-k}$
$\displaystyle=$ $\displaystyle q!(-1)^{q}\cdot\Big{\\{}\begin{array}[]{c}j\\\
q\end{array}\Big{\\}},$

where $\Big{\\{}\begin{array}[]{c}j\\\ q\end{array}\Big{\\}}$ is the Stirling
number of the second kind. The Stirling number of the second kind describes
the number of options to split a set of $j$ elements into $n$ non-empty
subsets, which is equal to 0 for $0\leq j<n$, and is equal to 1 for $j=n$. The
result holds true for all $j,n\in\mathbb{N}$ (Abramovitz and Stegun, 1972, p.
825).

Thus, we can further simplify the right hand side of (A.8), and write

$\displaystyle\sum_{k=0}^{q}(-1)^{k}\cdot\binom{q}{k}\cdot\sigma\left(\frac{k}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right)$
$\displaystyle=$
$\displaystyle\frac{\sigma^{(q)}(u_{\sigma})\cdot(x-x_{0})^{q}}{\mathtt{C}^{q}}\cdot(-1)^{q}$
$\displaystyle+\frac{(x-x_{0})^{q+1}}{\mathtt{C}^{q+1}\cdot(q+1)!}\sum_{k=0}^{q}(-1)^{k}\cdot
k^{q+1}\cdot\sigma^{(q+1)}(\xi_{k})\cdot\binom{q}{k},$

which in connection with (A.7) yields that

$\displaystyle(-1)^{q}\frac{\mathtt{C}^{q}}{\sigma^{(q)}(u_{\sigma})}\sum_{k=0}^{q}(-1)^{k}\cdot\binom{q}{k}\cdot\sigma\left(\frac{k}{\mathtt{C}}\cdot(x-x_{0})+u_{\sigma}\right)$
$\displaystyle=$
$\displaystyle(-1)^{q}\frac{\mathtt{C}^{q}}{\sigma^{(q)}(u_{\sigma})}\Big{\\{}\frac{\sigma^{(q)}(u_{\sigma})\cdot(x-x_{0})^{q}}{\mathtt{C}^{q}}\cdot(-1)^{q}$
$\displaystyle+\frac{(x-x_{0})^{q+1}}{\mathtt{C}^{q+1}\cdot(q+1)!}\sum_{k=0}^{q}(-1)^{k}\cdot
k^{q+1}\cdot\sigma^{(q+1)}(\xi_{k})\cdot\binom{q}{k}\Big{\\}}$
$\displaystyle=$
$\displaystyle(x-x_{0})^{q}+\frac{(-1)^{q}(x-x_{0})^{q+1}}{\mathtt{C}\cdot\sigma^{(q)}(u_{\sigma})\cdot(q+1)!}\sum_{k=0}^{q}(-1)^{k}\cdot
k^{q+1}\cdot\sigma^{(q+1)}(\xi_{k})\cdot\binom{q}{k}.$

In view of Assumption 1.2 and $\mathtt{C}$ being a fixed value, the proof is
now completed. $\blacksquare$

Proof of Lemma 3.2:

By Lemma A.2, we can reconstruct all of
$\\{m_{i}(\mathbf{x}\,|\,\mathbf{x}_{0})\\}$ as follows:

$\displaystyle\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})=\mathbf{D}\mathbf{A}(\mathbf{x}\,|\,\mathbf{x}_{0})$

where

$\displaystyle\mathbf{D}$ $\displaystyle=$
$\displaystyle\\{d_{ij}\\}_{d_{q}\times
d_{q}}=(\mathbf{d}_{1},\ldots,\mathbf{d}_{d_{q}}),$
$\displaystyle\mathbf{A}(\mathbf{x}\,|\,\mathbf{x}_{0})$ $\displaystyle=$
$\displaystyle\left([(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{1}]^{q},\ldots,[(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{d_{q}}]^{q}\right)^{\top}.$

Note that the rotation matrix $\mathbf{D}$ is determined by
$\boldsymbol{\alpha}_{j}$’s only, so they are fixed.

Apparently, we have an issue of identification here, because for example we
can arbitrarily rescale $\boldsymbol{\alpha}_{j}$’s, and modify $\mathbf{D}$
accordingly without changing $m_{i}(\mathbf{x}\,|\,\mathbf{x}_{0})$ as
follows:

$\displaystyle\mathbf{D}\mathbf{A}(\mathbf{x}\,|\,\mathbf{x}_{0})=\mathbf{D}\mathbf{B}\mathbf{B}^{-1}\mathbf{A}(\mathbf{x}\,|\,\mathbf{x}_{0}),$

in which $\mathbf{B}$ is full rank. Therefore, for the purpose of
identification, we regulate $\boldsymbol{\alpha}_{j}$’s as follows:

$\displaystyle[\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{d_{q}}]=\frac{1}{d+1}\cdot\mathbf{I}_{h}\mathbf{W}$
(A.10)

in which $\mathbf{I}_{h}=\operatorname*{\textrm{diag}}\\{h,\mathbf{I}_{d}\\}$,
and $\mathbf{W}$ is defined in the body of this lemma. As a result, for
$\forall j\in[d_{q}]$,

$\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}|(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}|$
$\displaystyle=$ $\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\frac{1}{d+1}|(h,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\mathbf{w}_{j}|$
$\displaystyle\leq$
$\displaystyle\frac{h}{d+1}\sqrt{d+1}\cdot\max_{j}\|\mathbf{w}_{j}\|=h,$

so we can invoke Lemma 3.1 later on.

Treating $(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}$
as a whole and using Lemma 3.1, we write

$\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\left|m_{i}(\mathbf{x}\,|\,\mathbf{x}_{0})-\sum_{j=1}^{d_{q}}d_{ij}\cdot\sum_{k=1}^{q+1}\gamma_{k}\cdot\sigma\left(\beta_{k}(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}+u_{\sigma}\right)\right|$
(A.11) $\displaystyle\leq$
$\displaystyle\sum_{j=1}^{d_{q}}|d_{ij}|\cdot\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\left|[(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}]^{q}-\sum_{k=1}^{q+1}\gamma_{k}\cdot\sigma\left(\beta_{k}(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}+u_{\sigma}\right)\right|$
$\displaystyle=$ $\displaystyle O\left(h^{q+1}\right),$

where the last line follows from Lemma 3.1. Also, $\boldsymbol{\beta}$,
$\boldsymbol{\gamma}$ and $u_{\sigma}$ are known as discussed in Remark 3.1.
Thus, we can further write

$\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\left|p(\mathbf{x}\,|\,\mathbf{x}_{0},\boldsymbol{\lambda})-\sum_{j=1}^{d_{q}}\boldsymbol{\lambda}^{\top}\mathbf{d}_{j}\cdot\sum_{k=1}^{q+1}\gamma_{k}\cdot\sigma\left(\beta_{k}(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}+u_{\sigma}\right)\right|$
(A.12) $\displaystyle=$ $\displaystyle\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\left|\boldsymbol{\lambda}^{\top}\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})-\sum_{j=1}^{d_{q}}\boldsymbol{\lambda}^{\top}\mathbf{d}_{j}\cdot\sum_{k=1}^{q+1}\gamma_{k}\cdot\sigma\left(\beta_{k}(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}+u_{\sigma}\right)\right|$
$\displaystyle\leq$
$\displaystyle\|\boldsymbol{\lambda}\|\cdot\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0},h}}\left\|\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})-\sum_{j=1}^{d_{q}}\mathbf{d}_{j}\cdot\sum_{k=1}^{q+1}\gamma_{k}\cdot\sigma\left(\beta_{k}(1,\mathbf{x}^{\top}-\mathbf{x}_{0}^{\top})\boldsymbol{\alpha}_{j}+u_{\sigma}\right)\right\|$
$\displaystyle=$ $\displaystyle O\left(h^{q+1}\right),$

where the last line follows from the facts that $d_{q}$ is fixed and
$\|\boldsymbol{\lambda}\|=O(1)$.

Finally, let
$\widetilde{\boldsymbol{\lambda}}=(\widetilde{\lambda}_{1},\ldots,\widetilde{\lambda}_{d_{q}})^{\top}$
with $\widetilde{\lambda}_{j}=\boldsymbol{\lambda}^{\top}\mathbf{d}_{j}$.
Further, in view of the definitions of
$\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})$ and
$(\boldsymbol{\pi}_{1},\ldots,\boldsymbol{\pi}_{d_{q}(q+1)})$ in the body of
this lemma and (A.10), the proof is then completed. $\blacksquare$

Proof of Lemma 3.3:

Before starting the proof, we introduce a few notations to facilitate the
development. First, recall that in the body of this theorem, we have defined

$\displaystyle\widetilde{s}(\mathbf{x}\,|\,\widetilde{\boldsymbol{\Lambda}})=\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
s(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}),$

where
$s(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})=(\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}\otimes\boldsymbol{\gamma})^{\top}\boldsymbol{\sigma}(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})$
by the definition of Lemma 3.2. Second, note that the leading terms of the
$q^{th}$ order Taylor expansion of $g(\mathbf{x})$ at each
$\widetilde{\mathbf{x}}_{\mathbf{i}}\in(-a,a)^{d}$ can be written as follows:

$\displaystyle\sum_{0\leq|\mathbf{J}|\leq
q}\frac{1}{\mathbf{J}!}\cdot\frac{\partial^{|\mathbf{J}|}g(\widetilde{\mathbf{x}}_{\mathbf{i}})}{\partial\mathbf{x}^{\mathbf{J}}}(\mathbf{x}-\widetilde{\mathbf{x}}_{\mathbf{i}})^{\mathbf{J}}=\boldsymbol{\lambda}_{\mathbf{i}}^{\top}\mathbf{m}(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}}):=p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}}),$

where the definition of $\boldsymbol{\lambda}_{\mathbf{i}}$ should be obvious
in view of the definition of
$\mathbf{m}(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})$ according to
(3.2).

We are now ready to start the proof, and write

$\displaystyle
g(\mathbf{x})-\widetilde{s}(\mathbf{x}\,|\,\widetilde{\boldsymbol{\Lambda}})$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})-\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
s(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})$
$\displaystyle+g(\mathbf{x})-\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}}).$

Note that by Lemma 3.2 we choose
$\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$ which fulfils the
relationship:
$\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}=\mathbf{D}\boldsymbol{\lambda}_{\mathbf{i}}$.
It is worth mentioning that although
$(\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})$ vary
with respect to $\mathbf{i}$, the rotation matrix $\mathbf{D}$ in facts is
solely determined by $\mathbf{W}$ of Lemma 3.2. Therefore, without loss of
generality, we can fix $\mathbf{W}$ over $\mathbf{i}$, as it is user chosen.
Then $\mathbf{D}$ remains the same in view of the proof of Lemma 3.2.

Next, we write

$\displaystyle\left\|\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})-g(\mathbf{x})\right\|_{\infty}$
(A.13) $\displaystyle=$
$\displaystyle\left\|\sum_{\mathbf{i}\in[M]^{d}}[p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})-g(\mathbf{x})]\cdot
I_{\mathbf{i},h}(\mathbf{x})\right\|_{\infty}$ $\displaystyle\leq$
$\displaystyle\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0\mathbf{i}},h}}|p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})-g(\mathbf{x})|=O(h^{p}),$

where the inequality follows from the definition of
$I_{\mathbf{i},h}(\mathbf{x})$, and the last step follows from Lemma A.1.
Also, we can obtain that

$\displaystyle\left\|\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
s(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})-\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot
p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})\right\|_{\infty}$
(A.14) $\displaystyle\leq$
$\displaystyle\sum_{\mathbf{i}\in[M]^{d}}I_{\mathbf{i},h}(\mathbf{x})\cdot\sup_{\mathbf{x}\in
C_{\mathbf{x}_{0\mathbf{i}},h}}|s(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})-p(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}},\boldsymbol{\lambda}_{\mathbf{i}})|=O(h^{q+1}),$

where the inequality follows from the definition of
$I_{\mathbf{i},h}(\mathbf{x})$, and the last step follows from Lemma 3.2 by
letting
$\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}=\mathbf{D}\boldsymbol{\lambda}_{\mathbf{i}}$.

Therefore, based on (A.13) and (A.14), we obtain

$\displaystyle\|g(\mathbf{x})-\widetilde{s}(\mathbf{x}\,|\,\widetilde{\boldsymbol{\Lambda}})\|_{\infty}=O(h^{p}),$

where $p=q+s$. The proof is now completed. $\blacksquare$

A.5.1 Model (1.1)

Proof of Lemma A.3:

First, we expand the expression of
$\widehat{\boldsymbol{\theta}}_{\mathbf{i}}$ as follows:

$\displaystyle\widehat{\boldsymbol{\theta}}_{\mathbf{i}}-\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$
$\displaystyle=$
$\displaystyle\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}y_{t}-\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$
$\displaystyle=$
$\displaystyle\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})-\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}}$
$\displaystyle+\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})]$
$\displaystyle+\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\varepsilon_{t}$
$\displaystyle=$
$\displaystyle\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})]$
$\displaystyle+\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\varepsilon_{t},$

where the third equality follows from the fact that
$I_{\mathbf{i},h}(\mathbf{x}_{t})I_{\mathbf{j},h}(\mathbf{x}_{t})=0$ for
$\mathbf{i}\neq\mathbf{j}$, and the definition of
$\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})$. Below,
we consider the terms on the right hand side one by one.

Recall the argument of Remark 3.5, so we immediately obtain that

$\displaystyle\|\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\mathbf{D}^{\top}\mathbf{H}-I_{\mathbf{i},h}(\mathbf{x}_{t})\mathbf{m}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\cdot\mathbf{H}\|=O(h),$
(A.15)

where the last line follows from the definition of $\mathbf{H}$ and
$\max_{j}|\mathbf{n}_{j}|=q$. For notational simplicity, we further denote

$\displaystyle
I_{\mathbf{i},h}(\mathbf{x}_{t})\mathbf{H}\cdot\mathbf{m}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}}):=\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})=(\widetilde{m}_{1}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}}),\ldots,\widetilde{m}_{d_{q}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}}))^{\top}.$
(A.16)

First, we consider
$\frac{1}{Th^{d}}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}$,
and write

$\displaystyle\frac{1}{Th^{d}}\sum_{t=1}^{T}\mathbf{H}\mathbf{D}E\left[\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right]\mathbf{D}^{\top}\mathbf{H}$
(A.17) $\displaystyle=$
$\displaystyle\frac{1}{h^{d}}\mathbf{H}\mathbf{D}E[\widetilde{\mathbf{x}}_{\mathbf{i},1}\widetilde{\mathbf{x}}_{\mathbf{i},1}^{\top}]\mathbf{D}^{\top}\mathbf{H}$
$\displaystyle=$
$\displaystyle\frac{1}{h^{d}}E[I_{\mathbf{i},h}(\mathbf{x}_{1})\mathbf{H}\cdot\mathbf{m}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\mathbf{m}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\cdot\mathbf{H}]\cdot(1+o(1))$
$\displaystyle=$ $\displaystyle\frac{1}{h^{d}}\int_{\mathbf{x}\in
C_{\mathbf{x}_{0\mathbf{i}},h}}\mathbf{H}\cdot\mathbf{m}(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\mathbf{m}(\mathbf{x}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\cdot\mathbf{H}\cdot
f_{\mathbf{x}}(\mathbf{x})\mathrm{d}\mathbf{x}\cdot(1+o(1))$ $\displaystyle=$
$\displaystyle
f_{\mathbf{x}}(\widetilde{\mathbf{x}}_{\mathbf{i}})\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}\cdot(1+o(1))>0$

where the first equality follows from Assumption 2.1, the second equality
follows from (A.15), and the last step follows from Assumption 2.2 and the
definition of $\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})$ given in (3.2).

Next, we consider the second moment. Define
$\mathbf{X}_{\mathbf{i}}=E[\mathbf{H}\mathbf{D}\widetilde{\mathbf{x}}_{t,\mathbf{i}}\widetilde{\mathbf{x}}_{t,\mathbf{i}}^{\top}\mathbf{D}^{\top}\mathbf{H}]$
for notational simplicity, and by the development similar to those given for
(A.19) below it is easy to see that

$\displaystyle
E\left\|\frac{1}{Th^{d}}\sum_{t=1}^{T}\mathbf{H}\mathbf{D}\widetilde{\mathbf{x}}_{t,\mathbf{i}}\widetilde{\mathbf{x}}_{t,\mathbf{i}}^{\top}\mathbf{D}^{\top}\mathbf{H}-\mathbf{X}_{\mathbf{i}}\right\|^{2}=o(1),$
(A.18)

so we omit the details for now.

Also, we note

$\displaystyle\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})]$
$\displaystyle=$
$\displaystyle\mathbf{D}^{\top}\mathbf{H}\left(\mathbf{H}\mathbf{D}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\mathbf{D}^{\top}\mathbf{H}\right)^{-1}\mathbf{H}\mathbf{D}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})]$
$\displaystyle=$
$\displaystyle\mathbf{D}^{\top}\mathbf{H}(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\Delta\mathbf{G},$

where

$\displaystyle\widetilde{\mathbf{X}}_{\mathbf{i}}=(\widetilde{\mathbf{x}}_{\mathbf{i},1},\ldots,\widetilde{\mathbf{x}}_{\mathbf{i},T})^{\top},$
$\displaystyle\Delta\mathbf{G}=(I_{\mathbf{i},h}(\mathbf{x}_{1})[g(\mathbf{x}_{1})-\widetilde{s}(\mathbf{x}_{1}\,|\,\widetilde{\boldsymbol{\Lambda}})],\ldots,I_{\mathbf{i},h}(\mathbf{x}_{T})[g(\mathbf{x}_{T})-\widetilde{s}(\mathbf{x}_{T}\,|\,\widetilde{\boldsymbol{\Lambda}})])^{\top}.$

By Lemma 3.3, it is easy to see that

$\displaystyle\frac{1}{Th^{d}}E\|\Delta\mathbf{G}\|^{2}=\frac{1}{h^{d}}E\left[I_{\mathbf{i},h}(\mathbf{x}_{1})[g(\mathbf{x}_{1})-\widetilde{s}(\mathbf{x}_{1}\,|\,\widetilde{\boldsymbol{\Lambda}})]^{2}\right]=O(h^{2p}).$

Then we can write

$\displaystyle\|(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\Delta\mathbf{G}\|^{2}$
$\displaystyle=$
$\displaystyle\Delta\mathbf{G}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H}(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\Delta\mathbf{G}$
$\displaystyle\leq$
$\displaystyle\frac{1}{Th^{d}}\lambda_{\text{max}}\left\\{\left(\frac{1}{Th^{d}}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H}\right)^{-1}\right\\}\Delta\mathbf{G}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H}(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\Delta\mathbf{G}$
$\displaystyle\leq$ $\displaystyle
O_{P}(1)\lambda_{\text{max}}\\{\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H}(\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\widetilde{\mathbf{X}}_{\mathbf{i}}\mathbf{D}^{\top}\mathbf{H})^{-1}\mathbf{H}\mathbf{D}\widetilde{\mathbf{X}}_{\mathbf{i}}^{\top}\\}\cdot\|\Delta\mathbf{G}\|^{2}/(Th^{d})$
$\displaystyle=$ $\displaystyle O_{P}(h^{2p}),$

where the first inequality follows from the exercise 5 on page 267 of Magnus
and Neudecker (2007), and the second inequality follows from (A.17) and
(A.18).

Based on the above development, we can conclude that

$\displaystyle\left\|\mathbf{H}^{-1}\mathbf{D}^{\top,-1}\left(\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}\widetilde{\mathbf{x}}_{\mathbf{i},t}^{\top}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbf{x}}_{\mathbf{i},t}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\,\widetilde{\boldsymbol{\Lambda}})]\right\|=O_{P}(h^{p}).$

Finally, in order to establish the asymptotic distribution, we just need to
focus on
$\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}$
in view of (A.15). Write

$\displaystyle
E\left[\left(\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}\right)\left(\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}\right)^{\top}\right]$
(A.19) $\displaystyle=$
$\displaystyle\frac{1}{Th^{d}}\sum_{t=1}^{T}\sum_{s=1}^{T}E[\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{s}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{t}\varepsilon_{s}]$
$\displaystyle=$
$\displaystyle\frac{1}{Th^{d}}\sum_{t=1}^{T}\sigma_{\varepsilon}^{2}E[\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}]$
$\displaystyle+\frac{1}{h^{d}}\sum_{t=1}^{T-1}(1-t/T)E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]$
$\displaystyle+\frac{1}{h^{d}}\sum_{t=1}^{T-1}(1-t/T)E[\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}],$

where the last two terms are the same up to a transpose operation.

The term
$\frac{1}{h^{d}}\sum_{t=1}^{T-1}(1-t/T)E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]$
on the right hand side can be bounded as follows.

$\displaystyle\|E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]\|$
$\displaystyle\leq$ $\displaystyle
O(1)\alpha(t)^{\nu/(2+\nu)}\left\\{E\|\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{1}\|^{2+\nu}\right\\}^{2/(2+\nu)}$
$\displaystyle\leq$ $\displaystyle
O(1)\alpha(t)^{\nu/(2+\nu)}\left\\{E\|\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{1}\|^{2+\nu}]\frac{1}{h^{d}}\right\\}^{2/(2+\nu)}\cdot(h^{d})^{2/(2+\nu)}$
$\displaystyle=$ $\displaystyle
O((h^{d})^{2/(2+\nu)})\alpha(t)^{\nu/(2+\nu)}.$

Thus, we have

$\displaystyle\left\|\frac{1}{h^{d}}\sum_{t=1}^{T-1}(1-t/T)E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]\right\|$
(A.20) $\displaystyle\leq$ $\displaystyle
O(1)\frac{1}{h^{d}}\sum_{t=1}^{d_{T}}\left\|E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]\right\|$
$\displaystyle+O(1)\frac{1}{h^{d}}\sum_{t=d_{T}+1}^{T}\left\|E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]\right\|$
$\displaystyle\leq$ $\displaystyle
O(1)h^{d}\sum_{t=1}^{d_{T}}E\left\|\frac{1}{h^{d}}I_{\mathbf{i},h}(\mathbf{x}_{1})\frac{1}{h^{d}}I_{\mathbf{i},h}(\mathbf{x}_{1+t})\varepsilon_{1}\varepsilon_{1+t}\right\|$
$\displaystyle+O(1)\frac{1}{h^{d}}\sum_{t=d_{T}+1}^{T}\left\|E[\widetilde{\mathbf{m}}(\mathbf{x}_{1}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{1+t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}\varepsilon_{1}\varepsilon_{1+t}]\right\|$
$\displaystyle\leq$ $\displaystyle
O(1)h^{d}d_{T}+O(1)\frac{(h^{d})^{2/(2+\nu)}}{h^{d}}\sum_{t=d_{T}+1}^{T}\alpha^{\nu/(2+\nu)}(t)$
$\displaystyle=$ $\displaystyle O(1)h^{d}d_{T}+O(1)\frac{1}{h^{\frac{\nu
d}{2+\nu}}}\sum_{t=d_{T}+1}^{T}\alpha^{\nu/(2+\nu)}(t)=o(1),$

where $\nu>0$ is defined in Assumption 2.1, and the last step follows from
Assumption 2.1 that we can choose $d_{T}$ to ensure

$\displaystyle d_{T}h^{d}\to
0\quad\text{and}\quad\sum_{t=d_{T}+1}^{T}\alpha^{\nu/(2+\nu)}(t)=o\left(h^{\frac{\nu
d}{2+\nu}}\right),$ (A.21)

which can be achieved by choosing $d_{T}=\lfloor T^{c_{1}}\rfloor$,
$h=T^{-c_{2}}$ and $\alpha(t)=t^{-c_{3}}$ for some suitable $c_{1}>0$,
$c_{2}>0$ and $c_{3}>0$, for example.

Thus, we can conclude that

$\displaystyle
E\left[\left(\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}\right)\left(\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}\right)^{\top}\right]$
(A.22) $\displaystyle=$
$\displaystyle\frac{1}{Th^{d}}\sum_{t=1}^{T}\sigma_{\varepsilon}^{2}E[\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})^{\top}]+o(1)$
$\displaystyle\to$
$\displaystyle\sigma_{\varepsilon}^{2}f_{\mathbf{x}}(\widetilde{\mathbf{x}}_{\mathbf{i}})\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}.$

Below, we further use small-block and large-block to prove the normality. To
employ the small-block and large-block arguments, we partition the set
$\\{1,\ldots,T\\}$ into $2k_{T}+1$ subsets with large blocks of size $l_{T}$
and small blocks of size $s_{T}$ and the last remaining set of size
$T-k_{T}(l_{T}+s_{T})$, where $l_{T}$ and $s_{T}$ are selected such that

$\displaystyle s_{T}\to\infty,\quad\frac{s_{T}}{l_{T}}\to
0,\quad\frac{l_{T}^{1+\nu}}{(Th^{d})^{\frac{\nu}{2}}}\to
0,\quad\text{and}\quad
k_{T}\equiv\left\lfloor\frac{T}{l_{T}+s_{T}}\right\rfloor,$

and $\nu$ is defined in Assumption 2.1.

For $j=1,\ldots,k_{T}$, define

$\displaystyle\boldsymbol{\xi}_{j}=\sum_{t=(j-1)(l_{T}+s_{T})+1}^{jl_{T}+(j-1)s_{T}}\frac{1}{\sqrt{h^{d}}}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t},\quad\boldsymbol{\eta}_{j}=\sum_{t=jl_{T}+(j-1)s_{T}+1}^{j(l_{T}+s_{T})}\frac{1}{\sqrt{h^{d}}}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t},$
$\displaystyle\boldsymbol{\zeta}=\sum_{t=k_{T}(l_{T}+s_{T})+1}^{T}\frac{1}{\sqrt{h^{d}}}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}.$

Note that $\alpha(T)=o(1/T)$ and $k_{T}s_{T}/T\to 0$. By direct calculation,
we immediately obtain that

$\displaystyle\frac{1}{T}E\left\|\sum_{j=1}^{k_{T}}\boldsymbol{\eta}_{j}\right\|^{2}\to
0\quad\text{and}\quad\frac{1}{T}E\left\|\boldsymbol{\zeta}\right\|^{2}\to 0.$

Therefore,

$\displaystyle\frac{1}{\sqrt{Th^{d}}}\sum_{t=1}^{T}\widetilde{\mathbf{m}}(\mathbf{x}_{t}\,|\,\widetilde{\mathbf{x}}_{\mathbf{i}})\varepsilon_{t}=\frac{1}{\sqrt{T}}\sum_{j=1}^{k_{T}}\boldsymbol{\xi}_{j}+o_{P}(1).$

By Proposition 2.6 of Fan and Yao (2003), we have as $T\to 0$

$\displaystyle\left|E\left[\exp\left(\frac{iw}{\sqrt{T}}\sum_{j=1}^{k_{T}}\boldsymbol{\xi}_{j}\right)\right]-\prod_{j=1}^{k_{T}}E\left[\exp\left(\frac{iw\boldsymbol{\xi}_{j}}{\sqrt{T}}\right)\right]\right|$
$\displaystyle\leq$ $\displaystyle 16(k_{T}-1)\alpha(s_{T})\to 0,$

where $i$ is the imaginary unit. In connection with (A.19)-(A.22), the Feller
condition is fulfilled as follows:

$\displaystyle\frac{1}{T}\sum_{j=1}^{k_{T}}E[\boldsymbol{\xi}_{j}\boldsymbol{\xi}_{j}^{\top}]\to\sigma_{\varepsilon}^{2}f_{\mathbf{x}}(\widetilde{\mathbf{x}}_{\mathbf{i}})\int_{[-1,1]^{d}}\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})\mathbf{m}(\mathbf{x}\,|\,\mathbf{0})^{\top}\mathrm{d}\mathbf{x}.$

Also, we note that

$\displaystyle E[\|\boldsymbol{\xi}_{1}\|^{2}\cdot
I(\|\boldsymbol{\xi}_{1}\|\geq\epsilon\sqrt{T})]$ $\displaystyle\leq$
$\displaystyle\\{E\|\boldsymbol{\xi}_{1}\|^{2\cdot\frac{2+\nu}{2}}\\}^{\frac{2}{2+\nu}}\left\\{E[I(\|\boldsymbol{\xi}_{1}\|\geq\epsilon\sqrt{T})]\right\\}^{\frac{\nu}{2+\nu}}$
$\displaystyle\leq$
$\displaystyle\\{E\|\boldsymbol{\xi}_{1}\|^{2+\nu}\\}^{\frac{2}{2+\nu}}\left\\{\frac{E\|\boldsymbol{\xi}_{1}\|^{2+\nu}}{\epsilon^{2+\nu}T^{\frac{2+\nu}{2}}}\right\\}^{\frac{\nu}{2+\nu}}$
$\displaystyle=$