Localized Neural Network Modelling of Time Series: A Case Study on US
  Monetary Policy
In this paper, we investigate a semiparametric regression model under the
context of treatment effects via a localized neural network (LNN) approach. Due
to a vast number of parameters involved, we reduce the number of effective
parameters by (i) exploring the use of identification restrictions; and (ii)
adopting a variable selection method based on the group-LASSO technique.
Subsequently, we derive the corresponding estimation theory and propose a
dependent wild bootstrap procedure to construct valid inferences accounting for
the dependence of data. Finally, we validate our theoretical findings through
extensive numerical studies. In an empirical study, we revisit the impacts of a
tightening monetary policy action on a variety of economic variables, including
short-/long-term interest rate, inflation, unemployment rate, industrial price
and equity return via the newly proposed framework using a monthly dataset of
the US.
Jiti Gao, Fei Liu, Bin Peng, Yanrong Yang
# A Localized Neural Network with Dependent Data: Estimation and Inference

###### Abstract

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
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:

\[y =g(\mathbf{x})+\varepsilon\quad\text{for continuous outcome}, \tag{1.1}\] \[y =\left\{\begin{array}{ll}1&g(\mathbf{x})-\varepsilon\geq 0 \\ 0&\text{otherwise}\end{array}\right.\quad\text{for binary outcome}, \tag{1.2}\]

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}\rightarrow\mathbb{R}\), where \(a\) is fixed1. The main goal of this paper is to infer \(g(\cdot)\) through a neural network architecture for both models. In addition, for Model (1.2) 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.2) 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 Krzyzak (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:

\[\{(y_{t},\mathbf{x}_{t})\,|\,t\in[T]\}, \tag{1.3}\]

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}\rightarrow[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. Provided a set of dependent time series data, how do we train a neural network, such as that in Figure 1?
2. Why are there only a small of number of neurons getting activated?
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. When the least absolute shrinkage and selection operator (i.e., LASSO) is employed, how do we define the set of true parameters?
5. Can any inference (such as a confidence interval) be established? and so forth.

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

Figure 1: A Neural Network with Sparsity
packages (e.g., Gunther 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 Gunther 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.2) 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.2) 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; \(\rightarrow_{P}\) and \(\rightarrow_{D}\) stand for convergence in probability and convergence in distribution, respectively. For \(g(\mathbf{x})\) defined in Model (1.1) and Model (1.2), 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}^{\boldsymbol{\alpha}_{1}}\cdots\partial x_{d}^{\boldsymbol{ \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}\),_

\[\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. _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. _Let Sigmoidal function_ \(\sigma(\cdot)\,:\,\mathbb{R}\to[0,1]\) _satisfy that_ 1. \(\sigma(\cdot)\) _is at least_ \(q+1\) _times continuously differentiable with bounded derivatives._ 2. _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.3).

**Assumption 2**.:
1. \(\{(\mathbf{x}_{t},\varepsilon_{t})\,|\,t\in[T]\}\) _are strictly stationary and_ \(\alpha\)_-mixing with mixing coefficient_ \[\alpha(t)=\sup_{A\in\mathcal{F}^{0}_{-\infty},B\in\mathcal{F}^{\infty}_{t}}|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}^{0}_{-\infty}\) _and_ \(\mathcal{F}^{\infty}_{t}\) _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. _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. _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\mathsf{c}\) _a.s._ _Model (_1.2_):_ \(\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

\[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.2), we may also introduce heterogeneity in a more general setup such as

\[\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.2).
### 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_

\[\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. 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. \(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. \(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. 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

\[\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\}, \tag{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

\[\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

\[\mathbf{m}(\mathbf{x}\,|\,\mathbf{x}_{0})=(m_{1}(\mathbf{x}\,|\,\mathbf{x}_{0 }),\ldots,m_{d_{q}}(\mathbf{x}\,|\,\mathbf{x}_{0}))^{\top}, \tag{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]\}, \tag{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\mathbf{c}\). Then a LNN of the type_

\[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_

\[\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_

\[(\boldsymbol{\pi}_{1},\ldots,\boldsymbol{\pi}_{d_{q}(q+1)})=\frac{\operatorname {diag}\{h,\mathbf{I}_{d}\}\mathbf{W}(\mathbf{I}_{d_{q}}\otimes\boldsymbol{ \beta}^{\top})}{d+1}+\begin{bmatrix}u_{\sigma}\\ \mathbf{0}_{d}\end{bmatrix}\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^{\text{th}}\) column of \(\mathbf{W}\)._

**Remark 3.2**.:
1. Note that we use vector form to rewrite each element of \(\boldsymbol{\sigma}(\mathbf{x}\,|\,\mathbf{x}_{0})\), i.e., \[\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. 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: \[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. 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. 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\to\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:

\[h=\frac{a}{M},\]

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

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

where \(x_{j}\) and \(x_{0i,j}\) are the \(j^{th}\) elements of \(\mathbf{x}\) and \(\mathbf{x}_{0i}\) 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:_

\[\|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_

\[\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}_{0i},\widetilde{\boldsymbol{\lambda}}_{\mathbf{i}})\quad\text{ with}\quad I_{\mathbf{i},h}(\mathbf{x})=I(\mathbf{x}\in C_{\mathbf{x}_{0i},h}).\]

**Remark 3.3**.:
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}_{0i}\).
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. 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.

### 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:

\[\mathcal{S}=\left\{\widetilde{s}(\mathbf{x}\,|\,\mathbf{\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\}, \tag{3.5}\]

in which \(\mathbf{\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(\mathbf{\Theta})=\sum_{t=1}^{T}[y_{t}-\widetilde{s}(\mathbf{x}_{t}\,|\, \mathbf{\Theta})]^{2}. \tag{3.6}\]

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

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

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

\[\widehat{g}(\mathbf{x}_{0})=\widetilde{s}(\mathbf{x}_{0}\,|\,\widehat{ \mathbf{\Theta}}). \tag{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 \(\mathbf{\widehat{\theta}_{\mathbf{i}}}\). To see this, we write

\[\frac{\partial Q(\mathbf{\Theta})}{\partial\mathbf{\theta}_{\mathbf{i}}} =-2\sum_{t=1}^{T}[y_{t}-\widetilde{s}(\mathbf{x}_{t}\,|\,\mathbf{ \Theta})]\cdot I_{\mathbf{i},h}(\mathbf{x}_{t})\cdot\frac{\partial s(\mathbf{x} _{t}\,|\,\mathbf{x}_{0\mathbf{i}},\mathbf{\theta}_{\mathbf{i}})}{\partial\mathbf{ \theta}_{\mathbf{i}}}\] \[=-2\sum_{t=1}^{T}[y_{t}-\widetilde{\mathbf{x}}_{\mathbf{i},t}^{ \top}\mathbf{\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\mathbf{\gamma}^{\top})\mathbf{\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

\[\mathbf{\widehat{\theta}_{\mathbf{i}}}=\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}. \tag{3.9}\]

**Remark 3.5**.: We note that by (A.10) and (A.11) of the online supplementary appendix:

\[\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\mathbf{\gamma}^{\top})\mathbf{\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

\[\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}\),_

\[\sqrt{Th^{d}}\widehat{\sigma}_{\mathbf{x}_{0}}^{-1}(\widehat{g}(\mathbf{x}_{0 })-g(\mathbf{x}_{0})+O_{P}(h^{p}))\rightarrow_{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}\mathbf{\Sigma}_{\mathbf{ x}_{0}}^{-1}\mathbf{H}\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})\), \(\mathbf{\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}=\mathrm{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. 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. 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. 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. By (3.8), we calculate \(\widehat{\varepsilon}_{t}=y_{t}-\widehat{g}(\mathbf{x}_{t})\) for \(t\in[T]\).
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. 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_

\[\sup_{w}|\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 \(\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\).

### On Model (1.2)

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

Direct calculation shows that

\[\Pr(y=1\,|\,\mathbf{x}) =\Phi_{\varepsilon}(g(\mathbf{x})),\] \[\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.3), the likelihood function is specified as follows:

\[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:

\[\log L(g) =\sum_{t=1}^{T}l_{t}(g(\mathbf{x}_{t}))\] \[=\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\}, \tag{3.10}\]

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:

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

which yields the following maximum likelihood estimator:

\[\widehat{\mathbf{\Theta}}=\operatorname*{argmax}_{\widetilde{s}\in\mathcal{S}} \log L(\widetilde{s}(\cdot\,|\,\mathbf{\Theta})). \tag{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{\mathbf{\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,_

\[\frac{1}{T}\sum_{t=1}^{T}[\Phi_{\varepsilon}(g(\mathbf{x}_{t}))-\Phi_{ \varepsilon}(\widetilde{s}(\mathbf{x}_{t}\,|\,\widehat{\mathbf{\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{\mathbf{\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

\[\frac{1}{T}\sum_{t=1}^{T}[g(\mathbf{x}_{t})-\widetilde{s}(\mathbf{x}_{t}\,|\, \widehat{\mathbf{\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.2).

**Theorem 3.3**.: _Suppose that Assumptions 1 and 2 hold. For \(\forall\mathbf{x}_{0}\in[-a,a]^{d}\),_

\[\sqrt{Th^{d}}\widetilde{\sigma}_{\mathbf{x}_{0}}^{-1}(\widetilde{g}(\mathbf{x }_{0})-g(\mathbf{x}_{0})+O_{P}(h^{p}))\rightarrow_{D}N(0,1),\]

_where \(\widetilde{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{\Sigma}_{\mathbf{x}_{0}}^{-1} \mathbf{H}\mathbf{m}(\mathbf{x}_{0}\,|\,\mathbf{x}_{0})\), and_

\[\widetilde{\mathbf{\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. For each bootstrap replication, we collect i.i.d. draws of \(\{\eta_{t}\,|\,t\in[T]\}\) from \(N(0,1)\), and calculate \[\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. 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.28) and (A.29) 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.2) 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_

\[\sup_{w}|\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 \(\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.

### 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

\[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.2) 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.

### Model (1.1)

Consider the following regression model:

\[y_{t}=g(\mathbf{x}_{t})+\varepsilon_{t}, \tag{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

\[(q,d,u_{\sigma}). \tag{4.2}\]

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

\[\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}, \tag{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 let2\(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.

Footnote 2: We 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.

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, \(\text{RMSE}_{g}\) converges to 0 and \(\text{CR}_{g}\) converges to 0.95. In particular, for \(\text{CR}_{g}\), 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, \(\text{RMSE}_{g}\) increases rather significantly, which is understood as the curse of dimensionality. The curse is also reflected in \(\text{CR}_{g}\), 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 \(\text{RMSE}_{g}\). However, the change in \(q\) does not alter \(\text{CR}_{g}\) 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\))
### Model (1.2)

Next, we consider the following data generating process:

\[y_{t}=\left\{\begin{array}{ll}1&g(\mathbf{x}_{t})-\varepsilon_{t}\geq 0\\ 0&\mbox{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

\[\operatorname*{argmin}_{\widetilde{\mathbf{x}}\in\mathcal{S}}[-\log L( \widetilde{\mathbf{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

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline  & & & \multicolumn{2}{c}{RMSE\({}_{g}\)} & \multicolumn{2}{c}{CR\({}_{g}\)} \\  & & \(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 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Simulation Results of Example 1.1
\(\mathbf{\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

\[\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 \(\text{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.2 (\(u_{\sigma}=-0.5,d=2\))
Figure 5: Simulation Results of Example 1.2 (\(u_{\sigma}=0.5,d=2\))
## 5 Empirical Study

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

### 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

\[\mathbf{x}_{t}=\mathbf{u}(\tau_{t})+\mathbf{s}_{t}+\mathbf{v}_{t}, \tag{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 Kaggle3. This dataset includes four time series (temperature, humidity, wind speed, atmospheric pressure)

\begin{table}
\begin{tabular}{c c r r r r} \hline \hline  & & \(T\) & RMSE\({}_{g}\) & RMSE\({}_{g}^{*}\) & CR\({}_{g}\) \\ \(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 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Simulation Results of Example 1.2
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}\)).

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

\[a=2,\quad q\in\{3,4,5,6\},\quad u_{\sigma}\in\{-0.5,0.5\} \tag{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

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

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

Figure 6: Climate Data of Delhi
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.
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://wildshares.com/about.php](https://wildshares.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

\[y_{t}=\left\{\begin{array}{ll}1&g(\mathbf{x}_{t})-\varepsilon_{t}\geq 0\\ 0&\mbox{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.

The entire time period we consider covers \(4^{th}\) January 2010 to \(31^{st}\) December 2019, in which we use the period from \(1^{st}\) July 2019 to \(31^{st}\) 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\)).

Figure 8: Plot of the Fitted Curve of the Test Period
After training our LNN, we calculate the mean absolute error (MAE) as follows:

\[\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.

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

\begin{table}
\begin{tabular}{l r r r} \hline \hline \(u_{\sigma}\setminus q\) & 3 & 4 & 5 \\ -0.5 & 0.5068 & **0.4384** & 0.5137 \\
0.5 & 0.5205 & **0.4110** & 0.4589 \\ \hline \hline \end{tabular}
\end{table}
Table 4: MAE for the Test Set

Figure 9: Clean Energy Index Data
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.

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.2) 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.

Figure 10: Plot of \(\widehat{g}(\mathbf{x})\) and associated Probabilities of the Test Period
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.
