EconBase
← Back to paper

Estimation and Inference for a Class of Generalized Hierarchical Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

68,492 characters · 11 sections · 52 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\newtheorem{corollary}{Corollary} \newtheorem{definition}{Definition} \newtheorem{lemma}{Lemma} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newtheorem{theorem}{Theorem} \newtheorem{example}{Example} \newtheorem{assumption}{Assumption}

\numberwithin{corollary}{section} \numberwithin{definition}{section} \numberwithin{equation}{section} \numberwithin{lemma}{section} \numberwithin{proposition}{section} \numberwithin{remark}{section} \numberwithin{theorem}{section}

\allowdisplaybreaks[4]

titlepage\begin{center} { \bf Estimation and Inference for \\a Class of Generalized Hierarchical Models} { {\sc Chaohua Dong$^{\ast}$, Jiti Gao$^{\dag}$, Bin Peng$^{\dag}$ and Yayi Yan$^\sharp$} $^{\ast}$Zhongnan University of Economics and Law\\ $^{\dag}$Monash University\\ $^\sharp$Shanghai University of Finance and Economics} \today \begin{abstract} In this paper, we consider estimation and inference for the unknown parameters and function involved in a class of generalized hierarchical models. Such models are of great interest in the literature of neural networks (such as BK2019). We propose a rectified linear unit (ReLU) based deep neural network (DNN) approach, and contribute to the design of DNN by i) providing more transparency for practical implementation, ii) defining different types of sparsity, iii) showing the differentiability, iv) pointing out the set of effective parameters, and v) offering a new variant of rectified linear activation function (ReLU), etc. Asymptotic properties are established accordingly, and a feasible procedure for the purpose of inference is also proposed. We conduct extensive numerical studies to examine the finite-sample performance of the estimation methods, and we also evaluate the empirical relevance and applicability of the proposed models and estimation methods to real data. {\em Keywords}: Estimation Theory; Deep Neural Network; Hierarchical Model; ReLU {\em JEL classification}: C14, C45, G12 \end{abstract} \end{center}

Introduction

In recent decades, there has been a notable emphasis on deep neural networks (DNNs) in the literature of econometrics, (e.g., cw1999,AC2003,SL2004,FLM2021). Initially applied in machine learning, DNNs have since expanded into various fields, such as economics, finance, social sciences, among others. Related to the applications of DNN, LBH2015 offer a comprehensive overview of practical topics, while Athey2019 discusses its capacity in social science. Additionally, BMR2021 and FMZ2021 offer reviews on recent methodological advancements.

As the most important part of DNN, a variety of activation functions have been proposed and investigated theoretically and numerically (DUBEY202292). The rectified linear activation function (ReLU) sees its popularity due to its simplicity and partial linearity:

eqnarray*[eqnarray* omitted — 77 chars of source]

ReLU is a piecewise linear function that will output the input directly if it is positive, otherwise, it will output zero. Compared to Sigmoid functions, ReLU has a low computational cost, which makes it efficient for large-scale neural networks practically. SH2020 and FG2022 for example establish some fundamental results with respect to using ReLU.

However, there are still properties related to ReLU that remain unknown. To be more specific, we now define a simple DNN using ReLU for activation function, and then briefly review the relevant literature.

definition[Simple DNN] For $\forall \mathbf{u},\mathbf{v} \in \mathbb{R}^n$, define the shifted activation function $\pmb{\sigma}_{\mathbf{v}}: \mathbb{R}^n \to \mathbb{R}^n$ as \begin{eqnarray*} \pmb{\sigma}_{\mathbf{v}} (\mathbf{u})=(\sigma (u_1 - v_1),\ldots, \sigma (u_n- v_n))^\top , \end{eqnarray*} where $u_j$ and $v_j$ stand for the $j^{th}$ elements of $\mathbf{u}$ and $\mathbf{v}$ respectively. A simple DNN with $m$ hidden layers that realizes the mapping $\mathbb{U}\, (\subseteq \mathbb{R}^{c_1})\mapsto \mathbb{R}$ is defined as follows: \begin{figure}[H] \begin{tikzpicture} \tikzstyle{annot} = [text width=10em, text centered] \node[rectangle, rounded corners, minimum size = 6mm, fill=orange!30] (Input-1) at (-0.15,-2) {$\mathbf{u} \in \mathbb{U}$}; \node[rectangle, rounded corners, minimum size = 6mm, fill=teal!10] (Hidden-1) at (2,-2) {$\mathbf{w}_1\, \pmb{\sigma}_{\mathbf{v}_{1}}(\mathbf{u})$}; \node[rectangle, rounded corners, minimum size = 6mm, fill=teal!10] (Hidden-3) at (4,-2) {$\cdots$}; \node[rectangle, rounded corners, minimum size = 6mm, fill=teal!10] (Hidden-4) at (6,-2) {$\mathbf{w}_m\, \pmb{\sigma}_{\mathbf{v}_{m}}(\cdot)$}; \node[rectangle, rounded corners, minimum size = 6mm, fill=purple!30] (Output-1) at (8.8,-2) {$\text{A scalar output}$}; \draw[->, shorten >=1pt] (Input-1) -- (Hidden-1); \draw[->, shorten >=1pt] (Hidden-1) -- (Hidden-3); \draw[->, shorten >=1pt] (Hidden-3) -- (Hidden-4); \draw[->, shorten >=1pt] (Hidden-4) -- (Output-1); \node[annot,above of=Input-1, node distance=1.1cm] (hl) {Input layer}; \node[annot,above of=Hidden-3, node distance=1.1cm] {Hidden layers}; \node[annot,above of=Output-1, node distance=1.1cm] {Output layer}; \end{tikzpicture} \end{figure} Mathematically, it is written as \begin{eqnarray*} \operatorname*{\mathcal{N}}(\mathbf{u}\, |\, \mathbf{W}_{m})&:=&\mathbf{w}_m \, \pmb{\sigma}_{\mathbf{v}_m} \, \cdots \,\mathbf{w}_{1}\, \pmb{\sigma}_{\mathbf{v}_{1}}(\mathbf{u}), \end{eqnarray*} where $\mathbf{W}_{m}:= \{\mathbf{v}_1,\ldots, \mathbf{v}_m; \mathbf{w}_1,\ldots, \mathbf{w}_m \}$, and the weighting matrices and shift vectors have the following dimensions: \begin{eqnarray*} &&\mathbf{w}_j is \left\{ \begin{array}{ll} c_1\times c_1 & for j=1,\\ c_j\times c_{j-1} & for 2\le j \le m-1 , \\ 1\times c_{m-1} & for j=m, \end{array} \right. \quadand\quad \mathbf{v}_j is \left\{ \begin{array}{ll} c_1\times 1 & \text{for }j = 1,\\ c_{j-1}\times 1 & \text{for }j\ge 2. \end{array}\right. \end{eqnarray*}

The current literature agrees that ReLU is designed to offer sparsity, which leads to computational efficiency (e.g., hoefler2021sparsity,glorot11a). However, there have only been few efforts to explain how sparsity should be defined and why it occurs. To the best of our knowledge, there are already two types of sparsity involved: (1) non-active neurons and (2) parameters that are not effective. Additionally, the literature implicitly agrees that $\mathbf{W}_{m}$ can be estimated through a minimization process. It is noteworthy that ReLU is piecewise linear, and it is not yet clear how to handle the accumulated (non-)differentiability through layers in both theory and practice. While the concern raised here actually exists in some well known software packages, to the best of our knowledge, no satisfactory treatment has been offered. We summarize a very brief survey about existing software packages in Appendix (ref).

Moving on to our discussion about modelling, when it comes to practical analysis using DNN based models and methods, the existing literature of model building primarily focuses on fully nonparametric models (e.g., cw1999 and many extensions since then), with only a few mentions of semiparametric settings (e.g., KK2017, BK2019, and references therein). It is not clear how to estimate and recover the parameters of interest for the generalized hierarchical models mentioned in Definition 1 of BK2019. There is a lack of investigations in this area of research, and it is also related to the (non-)differentiability of ReLU. As far as we know, these questions have not been thoroughly investigated.

Meanwhile, the current literature heavily focuses on independent and identically distributed (i.i.d.) data (e.g., BK2019, FLM2021), while largely neglecting the implications of asymptotic properties when dealing with dependent data. This is especially significant for applications in the fields of finance and economics, such as those studied by Gu2020, CLMZ2022, and FALLAHGOUL2024105574, where accounting for dependence can pose challenges in constructing inference. The literature on this topic dates at least back to NW1987, with a comprehensive review provided by Shao2015. In this paper, we aim to address this gap by training DNN with time series data, establishing asymptotic properties, and offering valid inference.

In what follows, in order to address the aforementioned concerns collectively, we follow the Definition 1 of BK2019 and consider a class of generalized hierarchical models of the form: for $1\leq t\leq T$,

eqnarray[eqnarray omitted — 218 chars of source]

where $f_\star(\mathbf{x})$ is an unknown $r$-dimensional function of $\mathbf{x}=\left(x_1, \cdots, x_r\right)^{\top}$, $\pmb{\theta}_{\star j}$ is a vector of unknown parameters of interest, each $\mathbf{z}_{jt}$ is a vector of $d_j\times 1$, for $d_j\ge 2$, observed time series, $\mathbf{z}_{jt}$ for $1\leq j\leq r$ may contain some common components, but $x_{jt}=\mathbf{z}_{jt}^\top\, \pmb{\theta}_{\star j}$ are all different for $1\leq j\leq r$, and $\varepsilon_t$ is an idiosyncratic error term.

When no misunderstanding arises, we write (ref) as

eqnarray[eqnarray omitted — 93 chars of source]

for notational simplicity, where $\mathbf{z}_t=\operatorname*{\normalfont\textrm{diag}}\{\mathbf{z}_{1t}^\top,\ldots,\mathbf{z}_{rt}^\top \}$ represents a block-wise diagonal matrix dataset of $r\times d$ dimensionality, and $\pmb{\theta}_\star=(\pmb{\theta}_{\star 1}^\top, \ldots, \pmb{\theta}_{\star r}^\top )^\top$ is a $d\times 1$ vector, in which $d=\sum_{j=1}^r d_j$. As discussed below Assumption (ref), meanwhile, a vector of lagged dependent variables of the form $\mathbf{Y}_{t-1} = \left(y_{t-1}, \cdots, y_{t-p}\right)^{\top}$ can be part of $\mathbf{z}_{t}$ when strict stationarity on $\mathbf{z}_{t}$ is imposed.

Throughout the rest of this paper, we suppose that $d_j$'s and $r$ are finite, although both $r$ and $d_j$ may be very large. We also assign the script $\star$ to the true parameters and the true function. For the purpose of identifiability and estimability of $\pmb{\theta}_\star$, we choose $\|\pmb{\theta}_{\star j} \|=1$ for all $j$'s, and assume that the first elements of $\pmb{\theta}_{\star j}$'s is positive.

Note that model (ref) is in the same spirit as Definition 1 of BK2019 and Eq (1.1) of Xia2002. Note also that Xia2002 haven't imposed the same set of restrictions as ours due to a different purpose, and BK2019 haven't focused on such identifiability and estimability issues that we are interested in this paper.

To the best of our knowledge, there have been no attempts to infer $(\theta_{\star 1}, \cdots, \theta_{\star r})$ in the relevant DNN literature, although KK2017 and BK2019 draw to hierarchical models. The main goals of this paper are to estimate both $f_\star(\cdot)$ and $\pmb{\theta}_{\star j}$'s jointly by proposing a ReLU based DNN approach. Moreover, the proposed DNN based estimation method offers a unified way to deal with the case where the dimensionalities, $d_j$ for $1\leq j\leq r$, $r$ and $d$, can all be large, although being fixed (as reported in Table 2 of Section 5 below for the case of $r=8$ and $d=16$). By contrast, the existing nonparametric methods suffer from the so-called “curse of dimensionality" issue when $r\geq 3$ and the sample size $T$ is not large enough.

In summary, this paper adds the following main contributions to the relevant literature:

enumerate• We enhance the design of DNN by i) providing more transparency for practical implementation, ii) defining different types of sparsity, iii) showing the differentiability, iv) pointing out the set of effective parameters, and v) offering a new variant of ReLU, etc. • We investigate a class of generalized hierarchical models, and propose a ReLU based DNN approach for the unified estimation of both the parameters of interest and the unknown function. • We allow our models and methods to be applicable to dependent time series data and then establish a valid implementable procedure for the purpose of inference. • We conduct extensive numerical results to validate the theoretical findings before we also demonstrate the empirical relevance and applicability of the proposed model and estimation method to real data.

The remainder of this paper is structured as follows. Section (ref) presents the design of DNN, and establishes some basic results which can be applied to a wide class of nonparametric models. Section (ref) considers the estimation of model (ref), and derives the necessary asymptotic theory accordingly. In Section (ref), we point out a few possible extensions. Section (ref) provides extensive numerical studies to examine the theoretical findings by both simulated and real datasets. We conclude in Section (ref) with a few remarks. Appendix (ref) presents some useful plots, notations, and preliminary lemmas. Due to the space limitation, we provide extra simulation results and plots in Appendix (ref) and Appendix (ref) receptively, and present all the proofs of the main results in Appendix (ref) of the online supplementary file.

Before proceeding further, we introduce a few symbols which will be repeatedly used throughout this paper.

Symbols & basic operations --- For $\forall w\in \mathbb{R}$, we let $\lfloor w\rfloor$ and $\lceil w\rceil$ be the largest and smallest integers satisfying $\lfloor w\rfloor \le w$ and $\lceil w\rceil\ge w$ respectively. For $\forall n\in \mathbb{N}$, we let $\mathbf{I}_n$, $\mathbf{1}_n$, and $[n]$ be a $n\times n$ identify matrix, a $n\times 1$ vector of ones, and a set $\{1,2,\ldots,n\}$ respectively. For $\pmb{\alpha} \in \mathbb{N}_0^r$ with $\mathbb{N}_0 =0 \cup \mathbb{N}$ and $\mathbf{x}\in \mathbb{R}^r$, we let

eqnarray*[eqnarray* omitted — 438 chars of source]

For a matrix $\mathbf{A}$, we let $\|\mathbf{A} \|$ and $\| \mathbf{A}\|_2$ define its Frobenius norm and Spectral norm respectively. Throughout, we write

eqnarray*[eqnarray* omitted — 183 chars of source]

Function operations & monomials --- Let $f(\mathbf{x})$ be a sufficiently smooth function defined on $\mathbb{U}\subseteq \mathbb{R}^r$, and define

eqnarray*[eqnarray* omitted — 477 chars of source]

We define a space of monomials:

eqnarray*[eqnarray* omitted — 131 chars of source]

of which the dimension is $\text{dim}\mathscr{P}_n =\binom{r+n}{r} :=r_n$ by direct calculation. Denote the basis of $\mathscr{P}_n$ by $\{\psi_1(\mathbf{x} ),\ldots, \psi_{r_n}(\mathbf{x})\}$, and let

eqnarray*[eqnarray* omitted — 107 chars of source]

For $\forall\mathbf{x}_0\in \mathbb{U}$, define the re-centred basis by $\{ \psi_1(\mathbf{x}\, |\, \mathbf{x}_0),\ldots, \psi_{r_n}(\mathbf{x}\, |\, \mathbf{x}_0) \}$, and let

eqnarray*[eqnarray* omitted — 163 chars of source]

Having these notation and symbols in hand, we are now ready to start our investigation.

DNN via ReLU

In this section, we present the design of DNN, and establish some basic results, which can be applied to a wider class of non- and semi-parametric models. Specifically, Section (ref) provides some basic definitions, while Section (ref) presents the detailed design.

Basic Definitions

Recall that we have defined a simple DNN (i.e., $\operatorname*{\mathcal{N}}(\cdot\, |\, \mathbf{W}_m)$) in Definition (ref). Building on it, we further define a pair-wise hierarchical DNN (referred to as HDNN hereafter), which plays an important role in what follows.

definition[HDNN] Define a mapping $\pmb{\ell}_{x,y}:(x,y)\mapsto \mathbb{U}$, where $x$ and $y$ are scalars. For $\forall\mathbf{u} =(u_1,\ldots, u_{2^q})^\top$ with $q\in \mathbb{N}$, the HDNN (written as $\operatorname*{\mathcal{N}}_{\pmb{\ell}} (\mathbf{u}\, |\, \mathbf{W}_m)$) is implemented as follows: \begin{enumerate} • Step 1 -- Divide $\mathbf{u}$ into pairs, and calculate $\operatorname*{\mathcal{N}} (\pmb{\ell}_{u_1,u_2} \, |\, \mathbf{W}_m),\ldots,\operatorname*{\mathcal{N}} ( \pmb{\ell}_{u_{2^q-1},u_{2^q}} \, |\, \mathbf{W}_m)$; • Step $n\ (2\le n \le q)$ -- Apply $\operatorname*{\mathcal{N}} (\pmb{\ell}_{x,y} \, |\, \mathbf{W}_m) $ to each pair of the outcomes from Step $n-1$. \end{enumerate}

To better see Definition (ref), we plot Figure (ref) for the purpose of visualization.

In Figure (ref) below, it is obvious that after $q$ steps, there is only one scalar left, which is the output of HDNN. Thus, the total number of hidden layers is $mq$. Although $\operatorname*{\mathcal{N}}_{\pmb{\ell}} (\mathbf{u}\, |\, \mathbf{W}_m)$ may appear complex, its parameters are entirely determined by $\mathbf{W}_m$ and $\pmb{\ell}_{x,y}$. Consequently, the number of effective parameters is significantly less than it appears. Finally, the ordering of the elements in $\mathbf{u}$ does not matter.

figure[figure omitted — 3,315 chars of source]

We then regulate the unknown function that is to be estimated.

definition[Smoothness] Let $p=\vartheta+s$ for some $\vartheta\in \mathbb{N}$ and $0<s\le 1$, and $\mathbb{U}\subseteq \mathbb{R}^r$. A function $f\,:\, \mathbb{U}\mapsto \mathbb{R}$ is called $(p, \mathscr{C})$-smooth, if for $\forall\boldsymbol{\alpha} \in \mathbb{N}_0^r$ with $\|\boldsymbol{\alpha}\|_1=\vartheta$ the partial derivative $f^{(\boldsymbol{\alpha})} (\mathbf{x})$ exists and satisfies that \begin{eqnarray*} \|f^{(\boldsymbol{\alpha})} (\mathbf{x})-f^{(\boldsymbol{\alpha})} (\mathbf{z}) \|_{\infty}^{\mathbb{U}} \le \mathscr{C} \| \mathbf{x} - \mathbf{z}\|^s , \end{eqnarray*} where $\mathscr{C}$ is a constant.

Definition (ref) is adopted from BK2019, and has different names in the literature (e.g., H{\"o}lder smoothness in SH2020).

The Design

We present the design of neural network in this subsection, and establish some fundamental results regarding function approximation. First, we present the following lemma building on Definition (ref).

lemmaFor $\forall(x,y)\in [0,1]^2$, construct a DNN with $m+3$ hidden layers: \begin{eqnarray*} \operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )\quadwith\quad \mathbf{W}^{\star}_{m+3} = \{\mathbf{v}_1^\star,\ldots, \mathbf{v}_{m+3}^\star; \mathbf{w}_1^\star,\ldots, \mathbf{w}_{m+3}^\star\}, \end{eqnarray*} where \begin{eqnarray*} &&\pmb{\ell}_{x,y}^m=\frac{1}{2}\left( \frac{x-y+1}{2}, \ x-y+1, \ x+y+2^{-m} , \ \frac{x+y}{2}, \ x+y,\ \frac{1}{2} \right)^\top ,\nonumber \\ && \mathbf{v}_k^\star =\left\{\begin{array}{ll} \mathbf{1}_2\otimes (0,2^{1-2k},0)^\top & for k\in [m+1]\\ -1 & for k\in\{ m+2,m+3\} \end{array}\right. ,\nonumber\\ && \mathbf{w}_k^\star =\left\{\begin{array}{ll} \mathbf{I}_2\otimes \widetilde{\mathbf{w}} & for k\in [m]\\ (-1,1,-1,1-1,1) &for k=m+1 \\ -1&for k=m+2\\ 1 &\text{for }k=m+3 \end{array}\right.,\nonumber \\ && \widetilde{\mathbf{w}} = \begin{pmatrix} \frac{1}{2} & -\frac{1}{2} & 0 \\ 1 & -1 & 0 \\ 1 & -1 & 1 \end{pmatrix} . \end{eqnarray*} Here, $\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )$ is piecewise linear in $x$ and $y$, and $\frac{\partial }{\partial x} [\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )]$ is defined accordingly. Let $\mathbb{C}_m=[0,1-2^{-m}]\times [0,1]$. The following results hold: \begin{enumerate} • $0\le \operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} ) \le 1$ on $[0,1]^2$, and $\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )=1$ at $(x,y)=(1,1)$, • $0\le \operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} ) -xy \le 2^{-m}$ on $\mathbb{C}_m$, • $\left\|\frac{\partial }{\partial x} [\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )] -y \right\|_{\infty}^{\mathbb{C}_m} \le 2^{-m-1}$. \end{enumerate}

In Lemma (ref), $x$ and $y$ are symmetric in the sense that one can interchange $x$ and $y$ without violating the above results. In addition, Lemma (ref) actually offers a generic result. For example, one can replace $y$ with a generic function, say $f(y)$, and the result still holds with obvious modification. Thus, $\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )$ offers a way to approximate different monomials practically, of which the space is the key to carry on nonparametric regression. To see this numerically, we plot Figure (ref) in Section (ref). Finally, it is worth mentioning that provided $\mathbb{C}_m$, $\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3} )$ always approximates $xy$ from positive side.

Based on Lemma (ref), we are then able to further estimate different monomials. Specifically, we provide the following lemma.

lemmaUsing $\operatorname*{\mathcal{N}} (\pmb{\ell}^m_{x,y} \, |\, \mathbf{W}^{\star}_{m+3})$ of Lemma (ref), define $\mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3})$ according to Definition (ref), where $r\ge 2$, $\mathbf{x}\in [0,h]^r$, $h\le 1 -\lceil\log_2 r\rceil \cdot 2^{-m}$, and $\pmb{\alpha} \in \mathbb{N}_0^r$. Then the following results hold: \begin{enumerate} • $0\le \mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3})\le 1$ uniformly on $\mathbf{x}\in [0,h]^r$, • $0\le \mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3}) -\mathbf{x}^{\pmb{\alpha}} \le 3^{\lceil \log_2 r\rceil-1} 2^{-m} $ uniformly on $\mathbf{x}\in [0,h]^r$, • $\| \frac{\partial}{\partial x_i}[\mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3}) -\mathbf{x}^{\pmb{\alpha}}]\|_\infty^{[0,h]^r} \le 3^{\lceil \log_2 r\rceil-1} \cdot \|\pmb{\alpha}\|_1\cdot 2^{-m} $ for $\forall i\in [r]$. \end{enumerate}

Lemma (ref) offers a specific range (i.e., $[0,h]^r$) in which DNN can approximate monomials reasonably well. It is interesting to note that the range includes non-negative quantities only, and $\mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3}) $ converges to $\mathbf{x}^{\pmb{\alpha}} $ from the positive side, which are consistent with Lemma (ref). The construction of $\mathcal{N}_{\pmb{\ell}^m} (\pmb{\ell}_{\mathbf{x}\,|\, {\pmb{\alpha}}} \, |\, \mathbf{W}^{\star}_{m+3})$ and Lemma (ref).3 together offer a theoretical justification for Adam and Keras algorithms in which differentiability is required. In Section (ref), we plot Figure (ref) and Figure (ref) for the purpose of demonstration.

We are now ready to consider the estimation of $f_\star (\mathbf{x})$. To proceed, we impose the following assumption to facilitate the development.

assumptionLet $f_\star(\mathbf{x})$ of (ref) be $(p, \mathscr{C})$-smooth on $\mathbf{x} \in [-a,a]^r$, where $a> 0$ and $r\geq 1$ are fixed.

As explained under Definition (ref), Assumption (ref) is commonly adopted in the literature. To proceed, we recall the notation and symbols defined at the end of Section (ref), and present the following lemma that approximates $f_\star(\mathbf{x})$.

lemmaUnder Assumption 1, for $\forall \mathbf{x}_0\in [-a,a)^r$, there exits $\pmb{\beta}_\star$ such that \begin{eqnarray*} \| f_\star (\mathbf{x}) - \mathbf{N} (\mathbf{x} \, | \, \mathbf{x}_0)^\top \pmb{\beta}_{\star} \|_{\infty}^{C_{\mathbf{x}_0,h}}=O(h^p +2^{-m}), \end{eqnarray*} where $ C_{\mathbf{x}_0,h} =\{\mathbf{x} \, |\, \mathbf{x} -\mathbf{x}_0\in [0,h]^r\}$ with $h\to 0$, $\pmb{\alpha}_{j}$'s are the power vectors of $\mathscr{P}_{r_\vartheta}$, and \begin{eqnarray*} \mathbf{N} (\mathbf{x} \, | \, \mathbf{x}_0 ) &=& \begin{pmatrix} \mathcal{N}_{\pmb{\ell}^m}(\pmb{\ell}_{\mathbf{x}-\mathbf{x}_0\, |\, \pmb{\alpha}_1 }\, |\, \mathbf{W}^\star_{m+3}) \\ \vdots \\ \mathcal{N}_{\pmb{\ell}^m}(\pmb{\ell}_{\mathbf{x}-\mathbf{x}_0\, |\, \pmb{\alpha}_{r_\vartheta}}\, |\, \mathbf{W}^\star_{m+3}) \end{pmatrix} . \end{eqnarray*}

It should be understood that in Lemma (ref), $\pmb{\beta}_\star$ varies with respect to $\mathbf{x}_0$. Although $h\to 0$ is not required in Lemma (ref), it is essential to include the condition here from the perspective of function approximation. Building on Lemma (ref), for $\forall M\in \mathbb{N}$, we can subdivide $[-a, a]^r$ into $M^r$ cubes of side length $h = \frac{2a}{M}$. For comprehensibility, we label these cubes by $C_{\mathbf{x}_{\mathbf{i}}}$ with $\mathbf{i}\in [M]^r$, where $\mathbf{x}_{\mathbf{i}}$ represents the point at the bottom left corner of each cube. Mathematically, $C_{\mathbf{x}_{\mathbf{i}}}$ is expressed as

eqnarray[eqnarray omitted — 130 chars of source]

Using the partition, we further present the first theorem of this paper.

theoremLet Assumption (ref) hold and $h\to 0$. There exists $\mathbf{B}_\star=\{\pmb{\beta}_{\star\mathbf{i}}\, |\, \mathbf{i}\in [M]^r\}$ such that \begin{eqnarray*} \| f_\star (\mathbf{x})-\mathscr{N} (\mathbf{x} \, | \, \mathbf{B}_\star)\|_{\infty}^{[-a,a]^r}=O(h^{p} +2^{-m}), \end{eqnarray*} where $\mathscr{N} (\mathbf{x} \, | \, \mathbf{B}_\star) =\sum_{\mathbf{i}} I(\mathbf{x}\in C_{\mathbf{x}_{\mathbf{i}}})\cdot \mathbf{N} (\mathbf{x} \, | \, \mathbf{x}_{\mathbf{i}})^\top \pmb{\beta}_{\star \mathbf{i}}$.

Theorem (ref) shows that we can recover $f_\star(\mathbf{x})$ on $[-a,a]^r$ via HDNN. Having Theorem 1 in hand, we are ready to work on the estimation of (ref) using data. In both theory and practice, the partial derivatives of an unknown function are also of great interest, as they allow us to further calculate marginal effects of some important variables. Along this line, we present some useful discussion and Corollary (ref) in Section (ref) later.

Estimation and an Asymptotic Theory

In this section, we consider the estimation of (ref). To facilitate the development, we define a few more symbols. Let

eqnarray[eqnarray omitted — 226 chars of source]

where $\theta_{j,1}$ stands for the $1^{st}$ element of $\pmb{\theta}_j$. As in (ref), we write $(\mathbf{z}_{1t}^\top \,\pmb{\theta}_1 ,\ldots, \mathbf{z}_{rt}^\top \, \pmb{\theta}_r)^\top = \mathbf{z}_t\, \pmb{\theta}$ when no misunderstanding arises.

Still, we partition $[-a,a]^r$ into $M^r$ cubes with side length $h =\frac{2a}{M}$, and work with $\{\mathbf{x}_{\mathbf{i}} \, | \, \mathbf{i}\in [M]^r\}$ as in (ref). To accommodate the structure of (ref), for $\forall\pmb{\theta}$ we map $\{\mathbf{x}_{\mathbf{i}} \, | \, \mathbf{i}\in [M]^r\}$ to $\{ \mathbf{z}_{\mathbf{i}}\, | \, \mathbf{z}_{\mathbf{i}}= \operatorname*{\normalfont\textrm{diag}}\{\mathbf{x}_{\mathbf{i}}\}\pmb{\theta}_{\operatorname*{\normalfont\textrm{diag}}}, \mathbf{i}\in [M]^r\}$, where $\pmb{\theta}_{\operatorname*{\normalfont\textrm{diag}}}=\operatorname*{\normalfont\textrm{diag}}\{\pmb{\theta}_1^\top, \ldots, \pmb{\theta}_r^\top\}$, and the dimension of $\mathbf{z}_{\mathbf{i}}$ is obviously consistent with $\mathbf{z}_t$. We then group $\mathbf{z}_t$'s using the following sets:

eqnarray[eqnarray omitted — 426 chars of source]

where the second equality follows from $\pmb{\theta}_{\operatorname*{\normalfont\textrm{diag}}}\,\pmb{\theta}_{\operatorname*{\normalfont\textrm{diag}}}^\top =\mathbf{I}_r$ and $\pmb{\theta}_{\operatorname*{\normalfont\textrm{diag}}}^\top \mathbf{1}_r =\pmb{\theta}$ by (ref). Here, $C_{\mathbf{i}}$'s are equivalent to $C_{\mathbf{x}_{\mathbf{i}}}$'s of (ref).

With these symbols in hand, we conduct the following minimization procedure:

eqnarray[eqnarray omitted — 137 chars of source]

where $Q_T (\pmb{\theta}, \mathbf{B}) = \frac{1}{T}\sum_{t=1}^T [y_t - \mathscr{N} (\mathbf{z}_t\, \pmb{\theta}\, |\, \mathbf{B}) ]^2$, $\pmb{\theta}\in \pmb{\Theta}$, and $\mathbf{B} = \{\pmb{\beta}_{\mathbf{i}} \, |\, \| \pmb{\beta}_{\mathbf{i}}\|<\infty \}$. By Theorem (ref), for $\forall \mathbf{x}_{0}\in [-a,a]^r$, $f_\star(\mathbf{x}_0)$ is naturally estimated by

eqnarray[eqnarray omitted — 115 chars of source]

To facilitate the theoretical development, we impose the following conditions.

assumption\begin{enumerate} • $\{(\mathbf{z}_t, \varepsilon_t )\, |\, t\in [T] \}$ are strictly stationary and $\alpha$-mixing with mixing coefficient \begin{eqnarray*} \alpha(k) = \sup_{A\in \mathcal{F}_{-\infty}^0, B\in \mathcal{F}_k^\infty} |\Pr(A)\Pr(B) -\Pr(A\cap B)| \ \ for $k\geq 1$ \end{eqnarray*} satisfying $\sum_{k= 1}^{\infty} \alpha(k)^{\nu/(2+\nu)} <\infty$ for some $\nu>0$, where $\mathcal{F}_{-\infty}^0$ and $\mathcal{F}_{k}^\infty$ are the $\sigma$-algebras generated by $\{(\mathbf{z}_{s}, \varepsilon_s )\, |\, s \leq 0\}$ and $\{(\mathbf{z}_{s}, \varepsilon_s )\, |\, s \geq t\}$, respectively. In addition, suppose that almost surely $E[\varepsilon_1\, |\,\mathbf{z}_1]=0$, $E[\varepsilon_1^2\, |\, \mathbf{z}_1]=\sigma_\varepsilon^2$, and $E[|\varepsilon_1|^{2+\nu} \, |\, \mathbf{z}_1]<\infty$. • $E[f_\star(\mathbf{z}_1\,\pmb{\theta}_\star)-f_\star(\mathbf{z}_1\,\pmb{\theta})]^2$ is uniquely minimized on $\pmb{\Theta}$, and $\inf_{ \pmb{\Theta}\times [-a, a]^r}\phi_{\pmb{\theta}} (\mathbf{w})\ge c_0>0$, where $\phi_{\pmb{\theta}} (\mathbf{w})$ defines the density function of $\mathbf{z}_1\,\pmb{\theta}$, and is Lipschitz continuous on $[-a, a]^r$. \end{enumerate}

Assumption (ref).1 is standard (see, for example, Chapter 2 of FanYao; Chapter 2 and Appendix A of gao2007), and it requires strict stationarity. As alluded below model ((ref)), under the strict stationarity, model ((ref)) covers the dynamic case where a vector of lagged dependent variables of the form $\mathbf{Y}_{t-1} = \left(y_{t-1}, \cdots, y_{t-p}\right)^{\top}$ can be part of $\mathbf{z}_{t}$. In such dynamic cases, the relevant literature (see, Howell1990; FanYao; gao2007; ttg2010, for example) shows that $y_t$ of model ((ref)) can be strictly stationary when $f_{\star}(\cdot)$ satisfies $\lim_{\|\mathbf{x}\|\rightarrow \infty} \frac{\left|f_{\star}(\mathbf{x})\right|}{\|\mathbf{x}\|}<1$.

In Assumption (ref).2, the condition about $E[f_\star(\mathbf{z}_1\,\pmb{\theta}_\star)-f_\star(\mathbf{z}_1\,\pmb{\theta})]^2$ is necessary even in the case $f_\star(\cdot)$ is fully known. As $f_\star(\cdot)$ also needs to be estimated, we impose one more condition on $\phi_{\pmb{\theta}}$, which can be easily justified. For example, if $\mathbf{z}_1\,\pmb{\theta}$ follows a multivariate normal/$t$ distribution, the condition automatically holds. Notably, although we only infer $f_{\star}(\cdot)$ on $[-a,a]^r$, it does not mean that $\mathbf{z}_1\,\pmb{\theta}$ has to belong to a compact set. In Section (ref), we discuss how to relax the condition on $a$ being finite.

Using Assumption (ref), we present the consistency in the following lemma.

lemmaUnder Assumptions (ref)-(ref), as $(h, Th^r)\to (0,\infty)$, \begin{enumerate} • $ \frac{1}{M^r}\sum_{\mathbf{i}\in [M]^r}\|\mathbf{H}(\widehat{\pmb{\beta}}_{\mathbf{i}}-\pmb{\beta}_{\star\mathbf{i}} ) \|^2=o_P(1)$, • $\|\widehat{\pmb{\theta}}-\pmb{\theta}_\star\|=o_P(1)$, \end{enumerate} where $\mathbf{H}=\operatorname*{\normalfont\textrm{diag}}\{ h^{\|\pmb{\alpha}_1\|_1},\ldots, h^{\|\pmb{\alpha}_{r_\vartheta}\|_1} \}$ with $\pmb{\alpha}_j$'s being defined in Lemma (ref).

After we have established Lemma (ref), we recall the notation introduced in Section (ref) before we establish the following asymptotic distribution in the second theorem of this paper.

theoremSuppose that $\pmb{\Sigma}_{11} +\pmb{\Sigma}_{12}+\pmb{\Sigma}_{12}^\top$ and $\pmb{\Sigma}_{11}$ are positive definite, where \begin{eqnarray*} \pmb{\Sigma}_{11} &=& \sigma_\varepsilon^2 E [I_{a,1}\, \mathbf{f}^{(1)}_\star (\mathbf{z}_1\, \pmb{\theta}_\star) \, \widetilde{\mathbf{z}}_1 \widetilde{\mathbf{z}}_1^\top \mathbf{f}^{(1)}_\star(\mathbf{z}_1\, \pmb{\theta}_\star) ],\nonumber \\ \pmb{\Sigma}_{12}&=&\lim_{T} \sum_{t=1}^{T-1}(1-t/T ) E[\varepsilon_1 \varepsilon_{1+t} \,I_{a,1}I_{a,1+t}\, \mathbf{f}^{(1)}_\star (\mathbf{z}_1\, \pmb{\theta}_\star)\, \widetilde{\mathbf{z}}_1 \widetilde{\mathbf{z}}_{1+t}^\top \,\mathbf{f}^{(1)}_\star(\mathbf{z}_{1+t}\, \pmb{\theta}_\star) ], \end{eqnarray*} in which $\widetilde{\mathbf{z}}_t = \mathbf{z}_t^{\top} \mathbf{1}_r$. Under Assumptions (ref) and (ref), as $(h, Th^r)\to (0,\infty)$, for $\forall\mathbf{x}_0\in [-a, a]^r$ \begin{equation} \pmb{\Sigma}_{\mathbf{x}_{0}}^{-1/2}\mathbf{D}_T\left[\begin{pmatrix} \widehat{\pmb{\theta}}-\pmb{\theta}_\star\\ \widehat{f}(\mathbf{x}_0) -f_\star(\mathbf{x}_0) \end{pmatrix}+\mathbf{c}_{bias} +O_P( 2^{-m}) \right]\to_D N(\mathbf{0}_{d+1}, \mathbf{I}_{d+1}), \end{equation} where $\|\mathbf{c}_{\rm bias}\|=O_P(h^p)$ with its detailed form available from (ref) of the supplementary appendix, $\mathbf{D}_T = \operatorname*{\normalfont\textrm{diag}}\{ \sqrt{T}\mathbf{I}_{d}, \sqrt{Th^r}\}$, and \begin{eqnarray*} \pmb{\Sigma}_{\mathbf{x}_{0}} &=&\operatorname*{\normalfontdiag}\{\mathbf{I}_d , \pmb{\psi}_{r_\vartheta}(\mathbf{x}_0)^\top\mathbf{H}^{-1}\}\cdot\pmb{\Sigma}\cdot\operatorname*{\normalfontdiag}\{\mathbf{I}_d , \mathbf{H}^{-1}\pmb{\psi}_{r_\vartheta}(\mathbf{x}_0)\} ,\nonumber \\ \pmb{\Sigma}&=&\operatorname*{\normalfontdiag}\{ \pmb{\Sigma}_{11}^{-1}( \pmb{\Sigma}_{11} +\pmb{\Sigma}_{12}+\pmb{\Sigma}_{12}^\top)\pmb{\Sigma}_{11}^{-1}, \sigma_\varepsilon^2 \mathbf{I}_{r_{\vartheta}}\}. \end{eqnarray*}

There are two bias terms involved in Theorem (ref). The term $2^{-m}$ arises due to the use of ReLU, and it can be negligible as long as $m$ is sufficiently large. Also, $m$ is proportional to the number of layers, so it explains why increasing the layers of DNN can improve estimation accuracy exponentially. The term $\mathbf{c}_{\rm bias}$ comes from the nonparametric nature of DNN. It is pointed out that the bias terms remain with the estimate of $\pmb{\theta}_\star$ as well, because the approximation error term involved in the above Theorem 1 is not necessarily negligible asymptotically, unless a type of under-smoothing condition: $\sqrt{T} \left(h^p + 2^{-m}\right)=o(1)$, is imposed.

In addition to the involvement of two bias terms in equation ((ref)), a long-run covariance matrix is involved in $\pmb{\Sigma}_{12}$. These terms make the asymptotic distribution in ((ref)) infeasible in practice. Therefore, we propose a bootstrap procedure for the purpose of inference for $\pmb{\theta}_\star$.

enumerate• For each bootstrap replication, we draw $\ell$-dependent time series $\{ \eta_t\, |\, t\in [T]\}$. where $E[ \eta_t]=0$, $E[ \eta_t^2]=1$, $E[ \eta_t^4]<\infty$, $E[\eta_t\eta_s]=a\left(\frac{t-s}{\ell}\right)$, $(\frac{1}{\ell},\frac{\ell}{\sqrt{T}})\to (0,0)$, $a(\cdot)$ is a symmetric kernel defined on $[-1,1]$ satisfying that $a(0)=1$ and $K_a(x)=\int_{\mathbb{R}} a(u)e^{-iux}du\ge 0$ for $x\in \mathbb{R}$. • Construct $y_t^* = \mathscr{N} (\mathbf{z}_t\, \widehat{\pmb{\theta}} \, |\, \widehat{\mathbf{B}}) + \widehat{\varepsilon}_t \cdot \eta_{t}$, where $\widehat{\varepsilon}_t = y_t - \mathscr{N} (\mathbf{z}_t\, \widehat{\pmb{\theta}} \, |\, \widehat{\mathbf{B}})$. Conduct estimation using $\{y_t^*,\mathbf{z}_t \}$ as in (ref) to obtain $\widehat{\pmb{\theta}}^*$. • Repeat Steps 1 and 2 $R$ times, where $R$ is sufficiently large.

The condition $K_a(x)=\int_{\mathbb{R}} a(u)e^{-iux}du\ge 0$ for $x\in \mathbb{R}$ essentially regulates the kernel function, which together with other restrictions imposed on $a(\cdot)$ are satisfied by a few commonly used kernels, such as the Bartlett and Parzen kernels. More choices of the kernel function can be found in Andrews1991 and shao2010 for example. To get $\eta_{t}$'s, one may use $N(\mathbf{0}, \pmb{\Sigma}_{\eta})$ with $\pmb{\Sigma}_{\eta} =\{a\left(\frac{t-s}{\ell} \right) \}_{T\times T}$ for ease of implementation.

corollaryLet the conditions of Theorem (ref) hold. If, in addition, $\sqrt{T}(h^{p}+2^{-m})\to 0$ as $(T, m)\rightarrow (\infty, \infty)$, then we have for $\forall \mathbf{x}_0\in [-a, a]^r$, \begin{eqnarray*} \sup_{\mathbf{w}}\left| \normalfont Pr\left\{ \sqrt{T}(\widehat{\pmb{\theta}}-\pmb{\theta}_\star)\le \mathbf{w} \right\} - \normalfont Pr^*\left\{ \sqrt{T}(\widehat{\pmb{\theta}}^*-\widehat{\pmb{\theta}})\le \mathbf{w}\right\} \right| =o_P(1), \end{eqnarray*} where $\text{\normalfont Pr}^*$ is the probability measure induced by the bootstrap procedure.

Because the biases come from diverse and intricate resources, we require a under-smoothing condition (i.e., $\sqrt{T}(h^{p}+2^{-m})\to 0$) in Corollary (ref). Up to this point, we have completed our theoretical development for model (ref).

Extensions

Before we conduct simulation results and discuss real data analysis, we discuss about how to smooth ReLU and also how to estimate marginal effects in this section.

Smoothed ReLU --- In order to improve the finite sample performance of ReLU, many variants have been proposed, such as Leaky ReLU, parametric ReLU, Gaussian-error linear unit (GELU), etc. The rise of these variants is most likely due to the (non)--differentiability. See DUBEY202292 for a comprehensive review. Here, although it is not the main focus of the paper, we also offer a new variant, i.e., smoothed ReLU.

We first introduce an assumption.

assumptionSuppose that $\phi(\cdot)$ is nonnegative, $\int \phi(u)dv=1$, and $ \int |u| \phi(u)du<\infty$.

It is readily seen that conventional density functions satisfy this assumption. Then for positive integer $s$, we define

equation[equation omitted — 108 chars of source]

The function sequence $\{\sigma_s(u)\}$ is called a regular sequence of $\sigma(u)$. For $\sigma_s(u)$, the following lemma holds.

lemmaLet Assumption (ref) hold. \begin{enumerate} • $\sigma_s(u), \ s=1,2,\cdots,$ are smooth, and $$\|\sigma_s(u)-\sigma(u)\|_{\infty}^\mathbb{R}\le c \cdot s^{-1}$$ for some absolute constant $c$. • Furthermore, let $\phi(\cdot)$ be defined on $[-1,1]$, and be symmetric. Then we obtain that \begin{eqnarray*} \left\{\begin{array}{ll} 0\le \sigma_s(u)-\sigma(u) \le O(1)s^{-1}, & for |u|\le s^{-1},\\ \\ \sigma_s(u)-\sigma(u) =0, & for |u|\ge s^{-1}. \end{array} \right. \end{eqnarray*} \end{enumerate}

It is worth emphasizing that $\sigma_s(u)$'s of Lemma (ref) are completely independent of the results in Sections (ref) and (ref). Therefore, using triangle inequality and letting $s$ be sufficiently large. We can easily extend the results of Sections (ref) and (ref) via this lemma. Here, we basically trade smoothness with an additional approximation error due to the use of $\sigma_s(u)$'s. In Section (ref), we further draw Figure (ref) to demonstrate this lemma.

Marginal Effects --- Lemma (ref) shows the feasibility of derivatives, which may enable us to compute the marginal effects of different key variables. For example, we can have the following corollary immediately.

corollaryLet $h\to 0$ and $p>r$. Under Assumption 1, for $\forall \mathbf{x}_0\in [-a,a)^r$, there exits $\pmb{\beta}_\star$ such that for $\forall \pmb{\delta}=(\delta_{1},\ldots, \delta_{r})^\top$ with $\delta_{j}\in\{ 0,1\}$ and $j\in [r]$, \begin{eqnarray*} \| f_\star^{(\pmb{\delta})}(\mathbf{x}) - \mathbf{N}_{\pmb{\delta}} (\mathbf{x} \, | \, \mathbf{x}_0)^\top \pmb{\beta}_{\star}\|_{\infty}^{C_{\mathbf{x}_0,h}}=O(h^{p-\|\pmb{\delta}\|_1} +2^{-m}), \end{eqnarray*} where \begin{eqnarray*} &&\mathbf{N}_{\pmb{\delta}} (\mathbf{x} \, | \, \mathbf{x}_0 ) =\begin{pmatrix} \pmb{\alpha}_1^{\pmb{\delta}}\, \mathcal{N}_{\pmb{\ell}^m}(\pmb{\ell}_{\mathbf{x}-\mathbf{x}_0,\, \overline{\pmb{\alpha}}_1 }\, |\, \mathbf{W}^\star_{m+3}) \\ \vdots \\ \pmb{\alpha}_{r_\vartheta}^{\pmb{\delta}}\, \mathcal{N}_{\pmb{\ell}^m}(\pmb{\ell}_{\mathbf{x}-\mathbf{x}_0, \, \overline{\pmb{\alpha}}_{r_\vartheta}}\, |\, \mathbf{W}^\star_{m+3}) \end{pmatrix} \quad and\quad\overline{\pmb{\alpha}}_{j} = \pmb{\sigma}_{\pmb{\delta}} (\pmb{\alpha}_j). \end{eqnarray*}

In a similar fashion to the establishment of Theorem (ref), we can estimate the marginal effects consistently.

Finite-Sample Evaluation

In this section, we conduct extensive simulation studies to examine our theoretical findings in Section (ref) before we demonstrate the applicability and the finite-sample properties and advantages of the estimation method using a real dataset in Section (ref).

Simulation Results for Section (ref)

In this subsection, we conduct simulations to examine the theoretical findings for model (ref). Specifically, the data generating process is as follows. Let

eqnarray*[eqnarray* omitted — 190 chars of source]

where $\varepsilon_t = \rho_\varepsilon \varepsilon_{t-1}+v_t$ and $f_\star(\mathbf{x}) =\frac{2}{r}\left(\sum_{j\text{ is odd}} (5x_j+\sin(2x_j)) +\sum_{j\text{ is even}}\exp(2.5 x_j)\right)$ with $\mathbf{x} = \left(x_1, \cdots, x_r\right)^{\top}$. We consider the case of $r=2$ before we also have a look at the cases of $r=4,8$ in the second part of this simulation study. It is well understood that when $r\ge 4$, a typical nonparametric method such as kernel method will break down. When $r=8$, understandably, our NN method takes much longer time to compute than that used for the case of $r=2$.

That said, in what follows, we let

eqnarray*[eqnarray* omitted — 266 chars of source]

where $z_{it,j}$ stands for the $j^{th}$ element of $\mathbf{z}_{it}$. It is easy to know that as we set $a=0.9$, some observations will apparently be left out. Therefore, the design matches the framework of Sections (ref) and (ref).

For each generated dataset, we conduct estimation as in (ref). After $J$ replications, we report the following measures:

eqnarray*[eqnarray* omitted — 439 chars of source]

where $\widehat{\pmb{\theta}}_j$ and $\widehat{f}_j(\cdot)$ stand for the estimates of $\pmb{\theta}_\star$ and $f_\star(\cdot)$ in the $j^{th}$ replication, and $\widehat{\theta}_{j,s} $ and $ \theta_{\star,s}$ are the $s^{th}$ elements of $\widehat{\pmb{\theta}}_j$ and $\pmb{\theta}_\star$. In addition, $\text{CI}_{j,s}$ is the 95% confidence interval for the $s^{th}$ element in the $j^{th}$ replication based on the bootstrap draws, and $\{\mathbf{x}_l^{\ast}\, |\, l\in 0\cup [L]\}$ are some pre-determined points: $\mathbf{x}_l^{\ast}= (-a + l\cdot 2a/L) \mathbf{1}_r.$ To implement the bootstrap procedure, we choose the Bartlett kernel (i.e., $y=(1-|x|)I(|x|\le 1)$), and follow the suggestion from such as PSU2011 to let $\ell =\lfloor 1.75T^{1/3} \rfloor$ for simplicity. We limit the number of bootstrap replications to 100 due to time and computation constraints. Also we let $T\in \{500, 1000, 2000\}$, $h\asymp T^{-1/4}$, and $m\in\{3,4,5\}$, where $m$ controls the number of layers of in each hierarchy. Finally, to exam our argument under Lemma (ref) we also repeat the above procedure using $\sigma_{s}(\cdot)$ by letting $s=32$, so $\sigma_{s}(\cdot)$ is sufficiently close to $\sigma(\cdot)$, and is smooth.

table[table omitted — 1,067 chars of source]

As shown in Table (ref), it is clear that both $\text{RMSE}_{\pmb{\theta}}$ and $\text{RMSE}_{f}$ converge to 0, and $\text{RMSE}_{\pmb{\theta}}$ converges at a much faster rate. It is not surprising, as $\widehat{\pmb{\theta}}$ enjoys a parametric rate. Additionally, we see that as the sample size goes up, $\text{CR}_{\pmb{\theta}}$ converges to 95% which is the nominal rate. Overall, $m=5$ offers better finite sample performance, which is also consistent with the finding of Section (ref). Due to the complexity of DNN and the slow rate of $\text{RMSE}_{f}$, we need reasonably large sample size, which is also required by some of the existing simulation designs involving i.i.d. data (e.g., DFLSP2021). Last but not least, we see $\sigma_s(\cdot)$ and $\sigma(\cdot)$ offer equivalent finite sample performance. Although we see differences in the third digits, but in view of the number of replications of the Monte Carlo design, we should not emphasize these differences too much. Again, to understand the differences between $\sigma_s(\cdot)$ and $\sigma(\cdot)$, one should conduct a systematic comparison as in HG2023.

As has been mentioned above, it is well understood that when $r\ge 4$, a typical nonparametric method such as kernel method will break down. Of course, nothing comes for free. When $r=4,8$, our NN method takes long time to compute due to large dimension involved in the calculation. Take $r=8$ for example, it takes a few hours to get one estimation done for $T\in \{1000, 2000\}$. As a result, it is impossible for us to calculate coverage rates using bootstrap method plus a large number of simulation replications. In practice, as one does not need to repeatedly implement the bootstrap draws, calculating confidence intervals is still feasible. That said, in what follows, we focus on bias and standard deviation, which are the focus of majority of DNN based studies anyway. Therefore, the following numerical results do not lose any generality. Additionally, we examine the case with $r=16$ in the online supplement Appendix (ref).

We now introduce some new measures as follows: {

eqnarray[eqnarray omitted — 521 chars of source]

where $\overline{\pmb{\theta}}= \frac{1}{J}\sum_{j=1}^J\widehat{\pmb{\theta}}_j$ and $ \overline{f}_j(\mathbf{x}_l)= \frac{1}{J}\sum_{j=1}^J\widehat{f}_j(\mathbf{x}_l)$}. Here $\mathbf{x}_l$'s are selected in the same way as in the main text. Note that ideally, we would like to select say 10 points on each dimension of $\mathbf{x}$ when evaluating the performance of $\widehat{f}$. However, it is not feasible practically, because $10^r$ with $r=8$ requires us to evaluate an extreme large number of points just for one simulation replication, which will definitely create lots of computation burden.

We summarize the results in Table (ref). As we can see, the method presented in the paper still provides reasonable finite sample performance when $r=4, 8$. A few facts emerge: 1) It seems that $\text{Bias}_{\pmb{\theta}}$ has less bias via $\sigma(\cdot)$, and $\text{Std}_{\pmb{\theta}}$ is smaller via $\sigma_{32}(\cdot)$. Thus, introducing $\sigma_s$ seems to be balancing bias and standard deviation if we just focus on the estimate of $\pmb{\theta}_\star$. 2) The same pattern can be found for the estimate of $f_\star$. In this case, the biases associated with $\sigma_{32}(\cdot)$ seem to be much larger, but they do not increase much when $r$ changes from 4 to 8. And 3) $m$ does not seem to play any serious role over all.

{

table[table omitted — 2,711 chars of source]

}

A Real Data Analysis

In this subsection, we examine the bond return predictability. Understanding the predictability of bond returns has always been a central topic in finance (e.g., LN2009, AEMS2020,Betal2023). The literature seems to agree that excess bond returns are forecastable by financial indicators such as yield spreads (AEMS2020), while also acknowledges the nonlinearity from modelling perspective (Betal2023). In addition, to avoid any variable selection problem, one often conducts principal component analysis (PCA) to collect some key factors from a large group of macro variables to be the key predictors of the forecasting model (LN2009). The above features naturally fit our example of Section (ref).

We consider two benchmark models in this section:

eqnarray[eqnarray omitted — 142 chars of source]

where $y_{t+1}$ denotes the 1-month log excess holding period return on a 2-month zero-coupon treasury bond. Here (ref) is a simple constant mean model, while (ref) is a typical one step ahead linear model including some key predictors. Specifically, $\mathbf{x}_t$ includes a set of variables adopted in Betal2023 and covers the period from December 1961 to December 2018: 1) yields spreads computed as the difference between the yield on a bond with 2 months to maturity and the implied yield on a one-month treasury bill obtained from CRSP; 2) forward spreads computed as the difference between the 2-month forward rate and the one-month yield; 3) the first principal component (PC) of yields obtained using 12, 24, 36, 48, and 60 month yields; 4) a linear combination of forward rates obtained from projecting 12, 24, 36, 48, and 60 month forward rates onto the mean excess bond return across the maturity spectrum; and 5) a linear combination of macroeconomic factors obtained using the FRED-MD database and estimated analogously to CP. As the first two predictors are observed directly, they form $\mathbf{z}_{1t}$. The last three predictors are obtained from estimation, so they form $\mathbf{z}_{2t}$.

Having the dataset ready, we first partition the data into a training set $S_{tr}$ (data from 1961 to 2011) and a test set $S_{te}$ (data from 2012 and 2018) (referred to as Partition I). For each model, we run regression using training set only to obtain the parameter estimates, and calculate $\widehat{y}_{t}$ for the test set. For the DNN approach, we also let $m=3,4,5$ to examine the impact of different layers practically.

To evaluate the performance of different methods, we calculate

eqnarray[eqnarray omitted — 228 chars of source]

where $\sharp S_{te}$ stands for the cardinality of $S_{te}$. The first measure is a typical root mean squared errors based on the test set, and the second measure examines the sign prediction and reports the percentage of correct sign prediction. There is a vast literature regarding sign prediction (e.g., CD2006 and many extensions since then), which argues that return is not predictable in general, but one may better forecast the sign of return. Although predicting the sign is not the focus of the paper, it is interesting to see how DNN performs practically along this line of research.

table[table omitted — 512 chars of source]

Alternatively, we consider a rolling window idea to partition the data (referred to as Partition II). The window size is always 50 years. For each window, we run regression using each model, and then use the estimated parameters to forecast the following year (i.e., the test set always includes data of one year). Again, we report RMSE and CS based on our forecasts.

figure[figure omitted — 145 chars of source]

We summarize the results in Table (ref). As shown in the table, DNN approach in general outperforms Model (ref) and Model (ref) in terms of both RMSE and CS. Here, smaller RMSE implies less forecast errors, while higher CS means better sign prediction. To further see the differences when using different layers (i.e., $m=3,4,5$) for DNN approach, we plot the point forecasts and their confidence intervals using Partition I for example in Figure (ref). In the figure, the actual $y_t$'s are marked using “x", the estimated values are represented using blue solid lines, and the red dotted-hash lines stand for the 95% confidence intervals. When $m=5$, DNN approach seems to offer a better coverage for the period from 2012 to 2015, while for the period from 2016 to 2018, it is quite hard to cover all the points regardless the value of $m$. This finding is consistent with Betal2023, in which the authors argue that the predictability varies with respect to time and economic conditions. Certainly, one may conduct a more comprehensive investigation to be added to the relevant literature. We do not pursue these results further in order not to deviate from our main goal, but would like to leave this for future research.

Conclusion with Discussion

Before concluding, we provide a few useful remarks. Note that we require $f_\star(\mathbf{x})$ to be defined on a compact set, but do not impose restriction on the range of $\{\mathbf{z}_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{z}_t$ follows a sub-Gaussian distribution. In this case, after some algebra we can relax the condition on $a$ to $\sqrt{\log (Th^r)} \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 LTG2016 for example, so we do not further elaborate it here.

The total number of layers that we require for Theorem (ref) is

eqnarray*[eqnarray* omitted — 141 chars of source]

The width of most layers in the hierarchy is 6. The sparsity occurs naturally in view of Definition (ref) and Lemma (ref).

Connection with Some Existing Studies --- Recently, DFLSP2021 and FLM2021 apply DNN to study treatment effects using micro datasets, and KN2020 apply DNN to study climate data. Our proposed model and method can be extendable to related topics, and the proposed estimation method offers a simple implementational procedure that complements existing studies. Specifically, we have a clear understanding of the nonlinear minimization process.

To conclude, we consider the estimation of both the parameters of interest and the unknown function involved in a class of hierarchical models in the literature of DNN. We contribute to the design of DNN by i) providing more transparency for practical implementation, ii) defining different types of sparsity, iii) showing the differentiability, iv) pointing out the set of effective parameters, and v) offering a new variant of rectified linear activation function (ReLU), etc. The asymptotic properties of the proposed estimates are derived accordingly, and they can be applied to a wider class of non- and semi-parametric models. Finally, we conduct extensive numerical studies to examine the theoretical findings.

Acknowledgment

The first author was supported by the National Natural Science Foundation of China (Grant 72073143) and the Fundamental Research Funds for the Central Universities at Zhongnan University of Economics and Law (2722022EG001). The second author was supported by Australian Research Council Discovery Project Program under Grant Number DP200102769. The third author was supported by Australian Research Council Discovery Project Program under Grant Number DP210100476. The fourth author was supported by the National Natural Science Foundation of China (Grant Number: 72303142) and Fundamental Research Funds for the Central Universities (Grant Numbers: 2022110877 & 2023110099).