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.
80,941 characters · 24 sections · 52 citation commands
Conditional nonparametric variable screening by neural factor regression
\def\spacingset#1{ {#1}} \spacingset{1.0} \def\r#1{\textcolor{red}{\bf #1}} \def\b#1{\textcolor{blue}{\bf #1}}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Neural networks, factor model, non-parametric regression, non-parametric tests, functional of derivatives, high-dimensionality.
\spacingset{1.2}
Variable screening is a powerful tool to expeditiously identify the set of predictors that potentially affect the regression outcome FanLv2008. It can reduce a very large number of predictors to a smaller and more manageable set. Then, on this reduced set, one can apply some more refined but computationally demanding variable selection methods such as Lasso, SCAD, Danzig selector Tibshirani1996,FanLi2001,CandesTao2007,FanLiZhangZou2020 and their non-parametric counterpart FAST-NN FanGu2023factor. Conditional marginal screening in parametric regression barut2016conditional further augments the screening by conditioning on a known set of useful predictors to reduce the impact of the correlations among the predictors, thus making the important {predictors more visible and reducing the false positive and false negative rates in the vanilla screening method}. Despite the aforementioned advances, conducting variable screening or selection in non-parametric regression with high-dimensional inputs remains a challenge due to curse of dimensionality, and deep learning offers a promising solution thanks to its ability to adapt to unknown low-dimensional structures in multivariate non-parametric regression.
Deep learning has achieved tremendous empirical successes in numerous applications LeCunBengioHinton2015deep, GoodfellowBengioCourville2016deep, for instance, in high-dimensional problems such as image recognition SimonyanZisserman2015, KrizhevskySutskeverHinton2017imagenet, deep reinforcement learning mnih2015human, and large language models kasneci2023chatgpt,thirunavukarasu2023large. There is now also a growing literature justifying theoretically the benefit of depth in deep neural networks Telgarsky2016benefits, Yarotsky2017, ElbrachterPerekrestenkoGrohsBolcskei2021 and their power in alleviating the curse of dimensionality in non-parametric regression via algorithmic learning of unknown low-dimensional structures within complex functions schmidt2020nonparametric, KohlerLanger2021. Now being a component of the standard toolbox for statisticians, deep neural networks may nevertheless not be efficient if the dimension of the predictors is very high due to the fundamental limit of multivariate non-parametric regression. Recently FanGu2023factor proposed a factor-augmented sparse throughput regression model that simultaneously leveraged the aforementioned adaptivity of deep neural networks and a factor model on the predictors to facilitate variable selection. To complement the variable selection effort in FanGu2023factor, DinhHo2020, and HoRichardsonTran2023, in this paper, we investigate the issue of conditional variable screening with deep neural networks when facing potentially very high-dimensional inputs.
We will assess conditional contribution of a candidate variable for screening by examining its partial derivative in the multivariate non-parametric regression function with a given set of variables and construct our screening test statistics on the moment generating function (MGF) of the smoothed partial derivative of the regression function estimator. We chose a deep neural network as our regression estimator due to its aforementioned algorithmic adaptation to the low-dimensional structure. Thus, quantifying and improving the performance of derivatives of deep neural networks are integral to our study, and also form an interesting topic of its own right given the importance of derivative estimation in non-parametric regression across diverse research domains and practical applications GijbelsGoderniaux2005, RondonottiMarronPark2007, horel2020significance. Traditional non-parametric derivative estimation in general follows one of the following three methods: empirical derivative-type estimation MullerStadtmullerSchmitt1987, DeBrabanterDeBrabanterGijbelsDeMoor2013, LiuDeBrabanter2020, kernel/local polynomial-type estimation GasserMuller1984,fan1996local, and series/spline-type estimation Stone1985,ZhouWolfe2000. However, the first method heavily relies on the existence of an order among the predictor samples and hence naturally applies when the predictor dimension is just one, and the last two methods, just like their original regression estimation counterparts, suffer from the curse of dimensionality problem {when facing predictors of a moderate dimension}.
Despite their deteriorated performances posed by high dimensionality, the traditional derivative estimation methods are relatively amenable to theoretical analysis due to their closed-form solutions. Excluding the empirical derivative method, which avoids explicitly fitting a regression function, the closed-form solutions for the last two methods can be attributed to the {close} connection between estimating the original regression functions and their accompanying derivatives. For instance, in the kernel and the spline methods, a {closed-form} derivative estimator can simply be obtained as the derivative of the original regression function estimator (see, for instance, the discussion between Eqs. (2) and (3) in ZhouWolfe2000). Thus, unsurprisingly, in these methods the quality of the derivative estimation closely follows the quality of the original regression function estimation. However, for deep neural networks, estimating the derivative of a regression function can be quite different from the task on the regression function itself due to its smoothness.
Take for example any candidate regression function estimator $m$ within the canonical neural network class, precisely defined in (ref) later, built from the popular ReLU activation function. The first order partial derivatives of $m$ are necessarily piecewise constant due to the piecewise linear nature of the ReLU function, and the second order partial derivatives of $m$ are necessarily zero almost everywhere, irrespective of the underlying truth that the function $m$ may attempt to recover. Numerically, this is easily demonstrable through standard software packages such as PyTorch. In addition, the first order partial derivatives of $m$ could exhibit convergence behaviors qualitatively different from the convergence behavior of $m$ itself as we will explain in Section (ref). Such phenomena are in sharp contrast to the traditional non-parametric estimators, and can in part be attributed to the features of deep neural networks: they are highly non-linear and lack easily interpretable closed-form expressions. {Moreover, existing asymptotic results on functionals acting on deep neural networks almost invariably assume some continuity of the functionals with respect to their inputs. However, when such functionals in effect act on the derivatives of neural networks, the assumed continuities {could} break down due to the aforementioned different convergence behaviors of the said derivatives.} We will explain and address this discrepancy as we progress through the paper.
Our main contribution in this paper is a non-parametric conditional variable screening test using deep neural networks when the ambient dimension of the inputs is potentially very high. In addition to complementing the variable selection methods, our contribution is also an advance over the aforementioned paper by barut2016conditional in that we work with non-parametric regression and our conditioning variables are not known in advance but instead are extracted from the inputs with the help of the factor model as in fan2022learning and FanGu2023factor. We note that horel2020significance has also conducted screening test based on partial derivatives (in low dimensions). However their study focused on single-layer neural networks, and hence does not benefit from the power nor reveal the intricacies of deep neural networks in derivative estimation.
In addition to the general procedures outlined above, we also make the following contributions which could be of independent interest:
Through simulation studies, we illustrate the favorable size and power performance of our test statistics, and the benefit of the smoothing operation in generating accurate derivative estimators. We also highlight the potential of our test statistics as a viable tool for model specification in nonlinear factor models via two empirical applications.
Let $\mathbb{N}$ denote the set of positive integers. For a vector $\mathbf{v} = (v_1, \ldots, v_d)^\top \in \mathbb{R}^d$, we let $\|\mathbf{v}\|_p = (\sum_{k=1}^d |v_i|^p)^{1/p}$ be the $\ell_p$ norm of $\mathbf{v}$. For a matrix $\mathbf{A} = (a_{j,k})_{1\le j\le l, 1\le k \le m}\in\mathbb{R}^{l\times m}$, we define the operator norm $\|\mathbf{A}\|_{\textup{op}} = \max_{\mathbf{v}\in\mathbb{R}^m: \|\mathbf{v}\|_2 = 1} \|\mathbf{A}\mathbf{v}\|_2$. We let $C$ denote an (absolute) constant that may change for each occurrence, and let “$\lesssim$” denote an inequality that holds up to such a multiplicative factor $C$; moreover, let $c$, $C$ and $M$ with super/subscripts denote constants with particular (though often non-specified) values. Limits are taken as $n\to\infty$ unless otherwise stated. For positive number sequences $(a_n: n\ge 1)$ and $(b_n: n\ge 1)$, we denote $a_n\lesssim b_n$ if there exists a positive constant $C$ such that $a_n/b_n\le C$ (for all $n$), and denote $a_n={\scriptstyle{\mathcal{O}}}(b_n)$ (resp. $a_n\sim b_n$) if $a_n/b_n\to 0$ (resp. $a_n\lesssim b_n$ and $b_n\lesssim a_n$). We use “$\rightarrow_{d}$” to denote convergence in distribution. Finally, let $\|\cdot\|_{L_2}=\|\cdot\|_{L_2(\mathbb{P})}$ denote the $L_2$ norm of the argument function with respect to the measure $\mathbb{P}$ to be formally introduced in Section (ref), so $\|f\|_{L_2}= \{ \int f^2 {\rm d}\mathbb{P} \}^{1/2}$. Sections in the supplement are labelled alphabetically.
We organize the remainder of our paper as follows. Section (ref) specifies our regression model and screening test, and develop our derivative estimators and test statistics. Section (ref) depicts the accompanying theoretical properties of the derivative estimators and test statistics from Section (ref). Section (ref) presents a simulation study. Section (ref) applies our test to an empirical example in asset pricing and another in macroeconomics time series. Section (ref) concludes and suggests several extensions. Additional results for the empirical application, proofs and supporting details are deferred to the supplementary materials.
Our screening test aims to tackle the high-dimensional regime where the ambient dimension $d$ of our observed predictors ${\bm{X}}\in\mathbb{R}^d$ can increase with the sample size $n$. Even with the remarkable capacity of deep neural networks to represent complex functions, our task is still infeasible if $d$ grows too fast. To address this issue, we assume that ${\bm{X}}$ admits the following factor model FanGu2023factor:
where ${\bm{F}} \in \mathbb{R}^r$ is a vector of latent factors, $\mathbf{B}\in\mathbb{R}^{d\times r}$ (usually $d\gg r$ for the high-dimensional ${\bm{X}}$) is the factor loading matrix, and ${\bm{u}}$ is the vector of the idiosyncratic noises. We let $(Y,{\bm{F}},{\bm{X}})$ have joint distribution $\mathbb{P}$ which also determines the distribution of ${\bm{u}}$ by (ref).
The conditional marginal screening is to see whether a component $X_j$ of ${\bm{X}}$ has additional contributions to the response variable $Y$ given ${\bm{F}}$. Then, the problem involves the working regression function $m_0(\cdot):\mathbb{R}^{r+1}\rightarrow\mathbb{R}$ defined through
The additional contribution of variable $X_j$ is measured through the partial derivative
We impose a blanket notational convention that a function with subscript $j$ denotes the partial derivative with respect to the last argument, which will almost always be $x_j$, while keeping the other argument (in this case, ${\bm{f}}$) fixed. Then, to examine whether $X_j$ has additional contribution, we propose the null and alternative hypotheses
We will assume that $(Y,{\bm{F}},{\bm{X}})$ is a $\mathbb{R}\times[-b,b]^{r+d}$ valued random vector. The function $m_0$ is akin to $\mathbb{E}[Y|X_j]$ in condition C in FanFengSong2011 on unconditional non-parametric screening, but here we exercise finer control over the potential role of $X_j$ by introducing ${\bm{f}}$ in $m_0=m_0({\bm{f}},x_j)$. Note that {working with} $m_0$ involves estimating the latent factors ${\bm{F}}$, which will be discussed in Section (ref). By Remark (ref) in the supplement, the screening hypotheses for $X_j$ in (ref) and (ref) are also equivalent to a significance test for the $j$-th idiosyncratic term $u_j$. We deliberately avoid working with the full conditional expectation of $Y$, namely $m_0^*({\bm{F}},{\bm{X}}) \equiv \mathbb{E}[Y|{\bm{F}},{\bm{X}}]$, because this is likely infeasible if the dimension of ${\bm{X}}$ is too high {without} additional structure. In contrast, the total number of variables $r+1$ in our working regression function $m_0$ and its derivative $m_{0,j}$ is much less than $d$, which allows us to circumvent the problem of testing a potentially very large number of coordinates of ${\bm{X}}$ simultaneously FanGu2023factor. Nevertheless, our model retains its validity even when the regression outcome $Y$ involves more coordinates of ${\bm{X}}$, as we will explain in Remark (ref).
Although we shall focus on the high-dimensional ${\bm{X}}$ case with growing $d$, our screening test also easily accommodates the low-dimensional ${\bm{X}}$ case, as we will comment in Section (ref). At the other end of the spectrum, our findings can be generalized to testing ${\bm{X}}$ over a fixed or expanding set of coordinates that form a subset of $\{1,\cdots,d\}$, as we will briefly outline in Section (ref).
We focus on multi-layer feed-forward neural networks that are fully connected between adjacent layers AnthonyBartlett1999. In this paper we exclusively consider the Rectified Linear Unit (ReLU) activation function defined as $\sigma=\sigma_{\text{ReLU}}(x)=\max \{x, 0\}$ due to its widespread popularity GoodfellowBengioCourville2016deep.
For now, we consider neural network functions with a generic input dimension ${\overline{d}}$. We finalize the structure of our neural networks by specifying the tuple $(L, \mathbf{k})$ where $L\in\mathbb{N}$ represents the number of hidden layers and the width vector $\mathbf{k}=({\overline{d}}, k_1, \ldots, k_L, k_{L+1}=1) \in \mathbb{N}^{L+2}$ specifies the number of nodes (i.e., neurons) in each hidden layer and the input/output dimensions. More precisely, such a neural network function $f:\mathbb{R}^{\overline{d}}\rightarrow\mathbb{R}$ with architecture $(L, \mathbf{k})$ is given by
here $L_\ell(\bm{z})=\mathbf{w}_\ell \bm{z} +\bm{b}_\ell$ is an affine map with weight matrix {$\mathbf{w}_{\ell}\in\mathbb{R}^{k_\ell \times k_{\ell-1}}$ if $\ell\geq 2$, $\mathbf{w}_{1}\in\mathbb{R}^{k_1 \times \overline{d}}$} and bias vector $\bm{b}_\ell\in\mathbb{R}^{k_\ell}$, and the function ${\overline{\sigma}}_\ell:\mathbb{R}^{k_\ell}\rightarrow\mathbb{R}^{k_\ell}$ applies the ReLU activation function $\sigma$ entry-wise. Additionally, we will truncate the (univarite) output of $f$ at a pre-specified constant level $M>0$ with the truncation operator $T_M(x) \equiv \textup{sgn}(x) \min\{|x|, M\}$. We denote such a collection of truncated $f$ by $\mathcal{F}_n({\overline{d}})$ where for brevity of notation we only retain the dependence on the input dimension ${\overline{d}}$:
where both $L$ and $\mathbf{k}$ can scale with the sample size $n$.
Next, we address the latency issue of the unobserved factor ${\bm{F}}$. For a given diversified weight $\mathbf{w} \in\mathbb{R}^d$, fan2022learning note from (ref) that within the projection $\mathbf{w}^T {\bm{X}}$, the projected idiosyncratic terms $\mathbf{w}^T {\bm{u}}$ will be negligible in high-dimension due to the law of the averages, so $\mathbf{w}^T {\bm{X}}$ yields an approximate linear combination of ${\bm{F}}$. To yield $r$ factors, we need at least $r$ projections. Since $r$ is {unknown}, one specifies an upper bound ${\overline{r}}$ {of} $r$.
Let $\mathbf{W}\in\mathbb{R}^{d\times {\overline{r}}}$ be a pre-trained diversified projection matrix as termed by Definition 3 in FanGu2023factor, and let $\mathbf{H}\equiv d^{-1}\mathbf{W}^\top \mathbf{B}\in\mathbb{R}^{{\overline{r}}\times r}$. Then, by (ref), we have $$ \widetilde{\bm{F}} \equiv d^{-1}\mathbf{W}^\top {\bm{X}} = \mathbf{H} {\bm{F}} + d^{-1}\mathbf{W}^\top {\bm{u}} \approx \mathbf{H} {\bm{F}} \in\mathbb{R}^{{\overline{r}}}, $$ by the law of the averages. Hence, $\mathbf{H}^\dagger \widetilde{\bm{F}} \approx {\bm{F}}$ under some appropriate conditions, where $\mathbf{H}^\dagger$ is the pseudo-inverse of $\mathbf{H}$. We call $\widetilde{\bm{F}}$ as the diversified factor, which is observable and a proxy of the latent factor ${\bm{F}}$. In practice, acquiring $\mathbf{W}$ in advance is necessary, either through domain expertise or data-driven methods. For instance, Proposition 1 in FanGu2023factor proposes estimating $\mathbf{W}$ via pretraining, where a tiny portion of size approximately $\log(n)$ of the samples is reserved to {extract the top ${\overline{r}}$ principal components in order to} construct a $\mathbf{W}$ that will satisfy Assumption (ref) with high probability. Henceforth we shall assume that $\mathbf{W}$ is pre-determined, exogenous and fixed.
With the proxies of the latent factors, for each given $j$, the coordinate of the candidate variable for screening, we can compute $\{Y_i, \widetilde{\bm{F}}_i, X_{i,j}\}_{i=1}^n$ based on the observed data, where $\widetilde{\bm{F}}_i= d^{-1}\mathbf{W}^\top {\bm{X}}_i$. Then we fit the neural network regression
where $\ell(y,\widetilde{\bm{f}},x_j;g)=\{y-g(\widetilde{\bm{f}},x_j)\}^2$ is the square loss and $\mathbb{P}_n$ is the empirical distribution. Intuitively, $\widehat g_n$ is a neural network estimation of
Indeed, Lemma (ref) shows that $g_0$ approximates $m_0$ well.
To test our screening hypotheses (ref) and (ref), we start from $\widehat g_{n,j}$, which by the convention in Section (ref) is the partial derivative of $\widehat g_n$ in (ref) with respect to $x_j$: $\widehat g_{n,j}(\widetilde{\bm{f}},x_j)=\partial\widehat g_n(\widetilde{\bm{f}},x_j) / \partial x_j$. Now, $\widehat g_{n,j}$ can be regarded as an estimator of the derivative $m_{0,j}$. We then consider the following initial MGF/exponentially tilted test statistic:
Partly owing to the fact that a MGF uniquely characterizes the distribution of the underlying random variable (specifically when considering the MGF over an interval including zero), MGF-based tests are popular in the literature EppsSingletonPulley1982,BaringhausEbnerHenze2017. For instance, they have been extensively employed in testing (multivariate) normality EbnerHenze2020. Under the null hypothesis $H_0: m_{0,j}=0$, we expect $\widehat g_{n,j}$ to be close to zero and so $\eta_t(\widehat g_n)$ is also centered around zero.
To enhance the quality of both the derivative estimator $\widehat g_n$ and the test statistic $\eta_t(\widehat g_n)$, we will further conduct the following refinements sequentially:
We focus solely on the smoothing step (ref) in this section, and will introduce step (ref) in Section (ref), and both (ref) and (ref) in Section (ref).
We first describe the smoothing operation applied to a generic function $g:\mathbb{R}^{\bar{r}+1}\rightarrow\mathbb{R}$ (or analogously $m:\mathbb{R}^{r+1}\rightarrow\mathbb{R}$). Let $K$ be a univariate continuously differentiable kernel supported on $[-1,1]$ with derivative $\dot K$. Then, for the generic function $g$, denote its smoothed version in the variable $x_j$ by $g^\textup{s}=g^\textup{s}(\widetilde{\bm{f}}, x_j) = \int g(\widetilde {\bm{f}}, z) K_h(x_j - z)dz$, where $K_h(\cdot)=K(\cdot/h) / h$ with a bandwidth parameter $h$. Now, the smoothed function $g^\textup{s}$ becomes differentiable in $x_j$ everywhere with the partial derivative given by
Accordingly, we let our initial smoothed derivative estimator based on $\widehat g_n$ be $\widehat g_{n,j}^\textup{s}$. In practice, we can select the bandwidth $h$ through cross-validation; see our Remark (ref).
To test the screening hypotheses (ref) and (ref), we proceed to estimate the score function corresponding to the statistic $\eta_t^\textup{s}(\widehat g_n)\equiv\mathbb{P}_n\{ \exp(t\,\widehat g_{n,j}^\textup{s}) - 1 \}$ now refined over (ref). This score function $\alpha_{t,n}^*:\mathbb{R}^{{\overline{r}}+1}\rightarrow\mathbb{R}$ satisfies, under the null and for any $L_2(\mathbb{P})$-integrable $\alpha$, the equality $\int \alpha\alpha_{t,n}^* {\rm d}\mathbb{P} = \int_{\Omega_h} \alpha_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P}$, where $\alpha_j^\textup{s}$ is the smoothed derivative of $\alpha$ and $\Omega_h=\{\omega\in\Omega: X_j(\omega)\in\mathcal{B}_h\}$ {is the interior sample space corresponding to} the interior set $\mathcal{B}_h=[-b+h,b-h]$. (Such an $\alpha_{t,n}^*$ exists because it is the Riesz representer of the directional derivative functional of our test statistic, as we will formally explain in Section (ref).) Therefore, the mean-square error (MSE) can be expressed as
Discarding the last term $\mathbb{E}[\alpha_{t,n}^{*2}(\widetilde{\mathbf{F}},X_j)]$, which does not depend on $\alpha$, the MSE can be estimated by the empirical loss
where $\mathcal{I}_h=\{i\in\{1,\dots,n\}: X_{i,j}\in\mathcal{B}_h\}$ is the interior index set. Hence, we estimate $\alpha_{t,n}^*$ by
which is a neural network estimator of the score function.
We tailor the loss function $\widehat{R}_n^{\textup{null}}$ and the estimator $\widehat\alpha_n$ to the null hypothesis (ref). One notable computational advantage is that $\widehat{R}_n^{\textup{null}}$ and hence $\widehat\alpha_n$ are independent of both the response and the value of $t$: A single optimization suffices to yield $\widehat\alpha_n$ for all $t$. However, the potential trade-off is a bias of $\widehat\alpha_n$ with respect to $\alpha_{t,n}^*$ induced under the alternative hypothesis, which will be characterized in Theorem (ref). Numerically, computing the right-hand side in (ref) is straightforward, as detailed in Lemma (ref) in the supplement.
Our refined statistic $\eta_t^\textup{s}(\widehat g_n)$ needs to be further centered to ensure asymptotic normality. To achieve this, ideally we adjust $\widehat g_n$ in the direction of the score function $\alpha_{t,n}^*$ to attain a smaller loss than in (ref). With the estimate $\widehat\alpha_n$ of $\alpha_{t,n}^*$ given by (ref), one naturally uses
to arrive at a debiased statistic $\eta_t^\textup{s}(\widecheck g_n)$, where
is chosen to minimize $\mathbb{P}_n \ell(\cdot;\widehat g_n + \delta_\textup{t} \widehat\alpha_n)$ in $\delta_\textup{t}$, as derived in Section (ref). We can interpret this update step as a form of the targeted machine learning VanderLaanRose2018 specialized to our deep neural network context. In Section (ref) we will show that $\widecheck g_n$ indeed satisfies an approximate minimization condition: for an infinitesimal stepsize $r_{\textup{inf},n}$,
for a small tolerance $b_n$. This will in turn center our test statistics.
To alleviate the technical possibility of the unboundnessed of our derivative estimator, we let $\Psi: \mathbb{R} \rightarrow \mathbb{R}$ be a bounded truncation function that also satisfies Assumption (ref); we provide an example of such a truncation function below the assumption. We then define the final (fixed-$t$) test statistic to be $\widecheck\eta_t^\textup{s}(\widecheck g_n)$ where the operator $\widecheck\eta_t^\textup{s}(\cdot)$ acts on functions $g:\mathbb{R}^{{\overline{r}}+1}\rightarrow\mathbb{R}$ as
Its standard error is estimated as $t \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}= t [ \sum_{i=1}^n \{ \widehat\epsilon_i \widehat\alpha_n(\widetilde{\bm{F}}_i,X_{i,j}) \}^2 / n ]^{1/2}$, where $\widehat\epsilon_i=Y_i-\widecheck g_n(\widetilde{\bm{F}}_i,X_i)$ is the $i$-th residual. This leads to a level-$\alpha$ test of the null hypothesis (ref) with the critical region
where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard normal distribution; the choice of this critical value is supported by the first part of (ref) in Theorem (ref).
To increase the power of our test, we aggregate the ensemble of test statistics $\widecheck\eta_t^\textup{s}(\widecheck g_n)$ over a range of $t$. We consider two types of aggregation: a sup-statistic and an integrated weighted square-statistic (or simply a square-statistic). In both cases, let $\delta > 0$ be a small but fixed constant, let $T > \delta$ be a potentially large but also fixed constant, and define $\mathcal{T}_\delta= [-T, -\delta] \cup [\delta, T]$. The two aforementioned aggregations of test statistics and their associated critical regions are as follows.
The critical values are supported by the last two parts of (ref) in Theorem (ref).
The major implementation steps of our partial derivative estimation and conditional screening test procedures in Section (ref) are summarized in Algorithm (ref). { \spacingset{1.5}
}
We first collect the assumptions most relevant to the construction and rates of our estimators $\widehat g_n$ and $\widehat g_{n,j}^\textup{s}$; we defer the more technical assumptions to Section (ref). In particular, Assumption (ref) imposes mild conditions on the data generating processes, and Assumption (ref) specifies the function class for $m_0$ and its associated neural networks.
In Theorem (ref) we establish the convergence of $\widehat g_n-m_0$ for regression function estimation and of $\widehat g_{n,j}^\textup{s} - m_{0,j}$ for derivative estimation. The following rates will appear:
where $c_{\textup{p}}$ is the constant appearing in (ref), $d_u$ appears in Assumption (ref), and $\delta_{\textup{f}} \to 0$ as $d\to \infty$. Moreover, because the target of the smoothed derivative estimator $\widehat g_{n,j}^\textup{s}$ is the smoothed derivative $m_{0,j}^\textup{s}$, not $m_{0,j}$, naturally a bias
occurs. However, this bias vanishes under the null in (ref) because there $m_{0,j}=0$ everywhere.
The proof of Theorem (ref) is deferred to Section (ref). With the rate $r_{m,n}$ in (ref) for the regression estimator, the rate (ref) for the smoothed derivative estimator follows from Lemma (ref) that applies to the smoothed derivatives of all estimators. Next, the rate $r_{m,n}$ consists of two parts: $\nu_n$ represents the combined effect of stochastic error and neural network approximation bias, while $\delta_{\textup{f}}$ represents the error induced by relying on the diversified factors in place of the latent factors. Theorem (ref) also covers as a special example the low-dimensional ${\bm{X}}$ case in Section (ref) by simply setting $\delta_{\textup{f}}=0$.
In this section, to complete our discussion in Section (ref), we formally introduce the score function $\alpha_{t,n}^*$ of our test statistics and the bias of $\widehat\alpha_n$ relative to $\alpha_{t,n}^*$. Define the population version of the operator $\widecheck\eta_t^\textup{s}(g)$ in (ref) as
and its associated directional derivative (in the direction $v=v(\widetilde{\bm{f}},x_j)$)
By (ref) in the proof of Lemma (ref), $\frac{\partial \widetilde\eta_t^\textup{s}(g_0)}{\partial g}[v] / t$ is a bounded linear functional of $v$. Hence, there exists a Riesz representer that becomes our $\alpha_{t,n}^*$ and that satisfies $\frac{\partial \widetilde\eta_t\left(g_0\right)}{\partial g}[v] / t = {\int v \alpha_{t,n}^* {\rm d}\mathbb{P}}$ for all $v$ Conway1990. In particular, by (ref), the effect of $\alpha_{t,n}^*$ is captured explicitly as
Because the loss $\widehat{R}_n^{\textup{null}}$ relies on the last step of (ref), $\widehat\alpha_n$ consistently estimates $\alpha_{t,n}^*$ under the null, as confirmed by Theorem (ref). However $\widehat\alpha_n$ may not be consistent for $\alpha_{t,n}^*$ under the alternative. Instead, it is not hard to show that it is consistent for a population limit $\alpha_n^{\textup{null}}$ which exists and satisfies $\int v \alpha_n^{\textup{null}} {\rm d}\mathbb{P} = \int_{\Omega_h} v_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P}$ for all $v$. See the remark below Eq. (ref) for details.
To obtain an intuitive idea of what the Riesz representer $\alpha_{t,n}^*$ could look like, assume that $\widetilde{\bm{F}}, X_j$ admit a joint density {$p(\widetilde{\bm{f}},x_j)$} that is differentiable in $x_j$ with the derivative being {$p_j(\widetilde{\bm{f}},x_j)$}. Then, by Lemma (ref), under the null hypothesis and in the limit $h\rightarrow 0$, $\alpha_{t,n}^*(\widetilde{\bm{f}},x_j)= - p_j(\widetilde{\bm{f}},x_j) / p(\widetilde{\bm{f}},x_j)$; in this case, if further $\widetilde{\bm{F}}, X_j$ are jointly Gaussian {on their support}, then $\alpha_{t,n}^*(\widetilde{\bm{f}},x_j)$ is simply a linear function in $\widetilde{\bm{f}}$ and $x_j$.
To ensure the rate of $\widehat\alpha_n$, our Assumption (ref) mirrors our earlier Assumptions (ref) and (ref) for estimating the regression function $m_0$ with deep neural networks and in particular imposes a hierarchical composition model on $\alpha_n^{\textup{null}}$ that takes the diversified predictors $(\widetilde{\bm{F}},X_j)$ as arguments. Then, Assumption (ref) places a mild condition on the smoothing bandwidth. As in Section (ref), we defer the more technical assumptions to Section (ref).
To characterize the bias of $\alpha_n^{\textup{null}}$ with respect to $\alpha_{t,n}^*$ under the alternative, define the random variable that represents the signal of the alternative as
Note that $Z_{t,j,h}$ always equals zero under the null hypothesis.
The proof of Theorem (ref) is deferred to Section (ref). Compared with the rate (ref) in Theorem (ref) for the regression estimator $\widehat g_n$, the rates for the score function estimator $\widehat\alpha_n$ mainly differ in two aspects: first, a factor $h^{-1}$ precedes $\nu_n$ which is the consequence of the smoothing operation in the loss function $\widehat{R}_n^{\textup{null}}$ in (ref); second, under the alternative hypothesis a bias $h^{-1} \| Z_{t,j,h} \|_{L_2}$ in (ref) is induced relative to $\alpha_{t,n}^*$.
Recall from Section (ref) that the adjusted estimator $\widecheck g_n$ in (ref) was introduced to center our test statistics through a suitable minimization condition (ref). In this section Proposition (ref) first shows that a small $\widehat{\delta}_\textup{t}$ is sufficient to arrive at $\widecheck g_n$ that satisfies (ref) with a small $b_n$.
The proof of Proposition (ref) is deferred to Section (ref). In the proposition, under the null the tolerance $b_n=r_{m,n} r_{\alpha,n}^{\textup{null}}$, which is faster than $n^{-1/2}$ under appropriate conditions and will ensure the asymptotic normality of our test statistics. Under the alternative, $b_n$ is not necessarily faster than $n^{-1/2}$, but is still fast enough to ensure {consistency under the local alternatives}.
We will start working mostly with the adjusted estimator $\widecheck g_n$ from now on. The next theorem is the counterpart of Theorem (ref) and shows that $\widecheck g_n$ and its smoothed derivative estimator $\widecheck g_{n,j}^\textup{s}$ maintain convergence rates similar to the unadjusted ones.
The proof of Theorem (ref) is deferred to Section (ref).
This section gives the result on the asymptotic null distributions of the test statistics. The results for consistency against the local alternatives are given in Section (ref).
The proof of Theorem (ref) is deferred to Section (ref). The extra condition (ref) in Theorem (ref) is natural due to the presence of the $\sqrt{n}$ scaling factor in the our tests and is mild. Note that if the bandwidth $h$ is held as a constant, then condition (ref) simply reduces to $\delta_{\textup{f}}={\scriptstyle{\mathcal{O}}}(n^{-1/2})$, and $\nu_n={\scriptstyle{\mathcal{O}}}(n^{-1/4})$ which is in turn implied by $\kappa>1/2$ for $\kappa$ in Assumption (ref) on the hierarchical composition model.
In this section, we conduct simulations to assess the size and power of our conditional screening tests across different scenarios. Recall that we summarize the implementations of these tests in Algorithm (ref). We consistently employ the Quartic/biweight kernel.
We choose $t=1$ for the fixed-$t$ test statistic, and simply set $\mathcal{T}_\delta = \{-1.25,-0.5,0.5,1.25\}$ in the sup statistic and the square statistic in (ref) and (ref) respectively. Then, we implement the sup statistic in (ref) as $\mbox{sup}_{t \in T_{\delta}}| {\widehat\eta_t^\textup{s}(\widecheck g_n)} |/t$, and the square statistic as $\widehat\chi^2_{\delta} = n \sum_{t\in \mathcal{T}_\delta} \left\{ {\widehat\eta_{t}^\textup{s}(\widecheck g_n)} \right\}^2 / t^2$ (we take $w(t)=1$ in (ref)).
We set $d=200$ and $d=400$ when $n=256$ and $n=512$ respectively, and the full regression model as $Y=m_0^*({\bm{F}},X_1,X_3)+\epsilon^*$ (see $m_0^*$ defined below (ref)). To illustrate the performance of our conditional screening test, we have designed $m_0^*$ so that a significant portion of its variation is accounted for by the factors. We further consider a nonlinear and a linear model of $m_0^*$. Specifically, under the null hypothesis, we set:
For our screening test, we always select $X_3$ as the variable of interest. Accordingly, under the alternative hypothesis we add signals in $X_3$ to $m_0^{*(null)}$ above:
The signal under the linear model is very weak and allows us to discern different test settings. While the signal under the nonlinear model may not appear weak at first, here the average derivative $\mathbb{E} m_{0,j}=0$, so detecting departure from the null critically depends on the higher-order, nonlinear effect of our MGF test statistics. Under the nonlinear model $r=4$ while under the linear model $r=5$. In both models we set ${\overline{r}}=r$.
Under the nonlinear model, $u_1$ is incorporated for a richer structure that specifically leads to the working regression function $m_0({\bm{F}},X_3)$ being different from the full regression function $m_0^*({\bm{F}},X_1,X_3)$ by design; see Remark (ref). We incorporate $u_1 = X_1 - \mathbf{B}_{1\cdot}{\bm{F}}$ instead of $X_1$ directly because $u_1$ will be enforced to be independent of $X_3$ under screening. Moreover, under the null hypothesis, the selection of $X_3$ is purely for clarity: under the nonlinear model, we could conduct our test on any $X_j$ for any $j\in\{2,\dots,d\}$, and under the linear model, for any $j\in\{1,\dots,d\}$.
For the high-dimensional predictor ${\bm{X}}$ in the factor model (ref), the factor loading matrix $\mathbf{B}$ is generated with i.i.d.\,Unif$[-\sqrt{3},\sqrt{3}]$ entries; the factor ${\bm{F}}$ and the idiosyncratic terms ${\bm{u}}$ both have i.i.d. $\mathcal{N}(0,0.6)$ entries. We further draw the noise as $\epsilon^*\sim\mathcal{N}(0,0.3)$. The quantities $\mathbf{B}$, ${\bm{F}}$, ${\bm{u}}$ and $\epsilon^*$ are all drawn independently. We pre-train $\mathbf{W}$ with samples of size $100$.
We set the hyper-parameters for the neural network fitting as follows: for the regression estimator $\widehat g_n$ in (ref), we employ neural networks with $L=5$ hidden layers, a common width of $k_{\textup{0}}=16$ per layer, and the ReLU activation function. Training proceeds over $800$ epochs, utilizing a batch size of $256$ and a constant learning rate of $0.005$. For the score function estimator $\widehat\alpha_n$ in (ref), we maintain the same neural network fitting parameters, except that we set $L=2$, batch size to $64$ and the number of epochs to $400$. We employ early stopping as the only regularization technique, and terminate training if there's no improvement after 20 epochs ($\text{patience}=20$) on a validation set.
We conduct our tests at significance levels of either $5\%$ or $10\%$ under both the nonlinear and the linear models, resulting in a total of four combinations summarized in Tables (ref) to (ref). In each table we present in alternating rows the sizes and powers (the latters in parentheses) of the tests in (ref), (ref) and (ref), and of the same tests but without centering (that is, tests employing the non-adjusted estimator $\widehat g_n$), under different sample sizes. The columns, arranged from left to right, correspond to the four different bandwidths. Each entry is calculated based on $500$ Monte Carlo repetitions.
We make a few observations from the results in the tables. First, centering generally improves the size of our tests without affecting their power. Moreover, under the null, sizes arrive at the nominal level if we further increase the bandwidth; this is especially evident for the sup test, which naturally tends to reject more often (under both the null and alternative). This shows the improvement of the power of the centered, smoothed test compared to the initial test (ref) if both tests are calibrated to have the same size. Of course, if we increase $h$ more than $1.5$, we expect the biases will arrive {under the alternative}. In the nonlinear case (Tables (ref) and (ref)), the ensemble tests by aggregating various $t$ also outperform the fixed-$t$ test in power.
In our first empirical analysis, we examine a comprehensive dataset containing both returns and specific characteristics of firms. Twelve monthly returns in $2021$ for the 100 largest financial institutions are drawn from the CRSP database as our response variable, along with a set of $d=48$ firm-specific characteristics as predictors; the resulting sample size is {$n=1200$}. The dataset's foundation is credited to chen2021open.
We select the smoothing bandwidth from the sequence $\{0.4, 0.6, \dots, 1.8, 2.0\}$ using the bandwidth selection algorithm suggested in Remark (ref). We set ${\overline{r}}=5$ by the conventional practice in finance, the hidden size $k_{\textup{0}}=24$, and the learning rate $\tau=10^{-3}$; we keep the other neural network fitting parameters identical to Section (ref).
In Figure (ref) we plot for each of the $44$ continuous predictors the intervals whose right ends are the test statistics normalized by their corresponding critical values in (ref), (ref) and (ref). Specifically, we plot the intervals $$ [ 0, n^{1/2} |\widecheck\eta_t^\textup{s}(\widecheck g_n)| / ( t z_{1-\alpha/2} \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)} ) ], \quad [ 0, \widehat Z/(z_{1-\alpha/2}\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)})], \quad [ 0, \sqrt{\widehat\chi^2}/(\sqrt{\chi^2_{1,1-\alpha}}\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)})],$$ for the fixed-$t$, sup, and square tests respectively. We also plot the same intervals but for the non-centered versions of these tests. In all cases, an interval covering one indicates statistical significance. The results suggest that the presence of some idiosyncratic terms related to tail risk, such as ReturnSkewCAPM and DownsideBeta, exert additional influence on returns and complement the explanatory power of the five factors.
In this section, to substantiate Section (ref), we illustrate our conditional screening test with another empirical application, this time on the macroeconomics dataset FRED-MD introduced in mccracken2016fred and also later studied in FanGu2023factor. The dataset collects $d=127$ monthly U.S.\,macroeconomic variables, such as unemployment rate and real personal income, starting from 1959/01. It is shown in mccracken2016fred that these variables can be explained well by several latent factors.
Our analysis setup is in general similar to Section C in FanGu2023factor. Our target variables are UEMP15T26, TB3SMFFM or TB6SMFFM, and we aim to identify which variables contribute to predicting the target variables beyond the latent factors. The variable UEMP15T26 represents the civilians unemployed for $15-26$ weeks. The variable TB3SMFFM (TB6SMFFM) measures the 3-month (6-month) treasury bill rate minus the effective federal funds rate. For each target response variable $y_{t+1}$ in \{UEMP15T26, TB3SMFFM, TB6SMFFM\}, we regress $y_{t+1}$ on $\bm{x}_t\in\mathbb{R}^{127}$ where $\bm{x}_t$ is the vector of all variables at the previous month. We choose the $n=330$ valid sample pairs $\{\bm{x}_t, y_{t+1}\}$ between January 1980 and July 2022. We employ the same neural network fitting parameters as Section (ref) except that we revert back to a hidden size $k_0=16$ as in our simulation studies.
For each variable, we plot the same six intervals derived from our tests as in Section (ref) and we recall that an interval covering one indicates statistical significance. Our analysis, based on the variables categorized by mccracken2016fred, underscores the predominant impact of a limited set of variables on the three response variables. As there are 127 variables, it takes four sub-figures to display a figure that is similar Figure (ref) for each given response variable. Figures (ref) {to} (ref) in the appendix display the intervals, with the response variable being UEMP15T26 (civilian unemployment). Notably, most variables fail to reach the $5\%$ significance level, except for some variables such as RETAILx, NDMANEMP, and COMPAPFFx. Additional significant variables for the response variables TB3SMFFM and TB6SMFFM are illustrated from Figures (ref) to (ref). As for TB3SMFFM, there are many significant variables. For example, the variables such as IPDMAT in the output and income variable group and DTCTHFNM in the prices group contribute additionally to the factors. TB6SMFFM demonstrates slightly higher susceptibility to variables compared to TB3SMFFM, which, in turn, is notably influenced by variables related to the stock market, employment, interest rates, and money and credit. Our overall findings validate our regression model based on a latent factors plus sparse idiosyncratic terms structure.
We have introduced a conditional variable screening test for non-parametric regression using deep neural networks; the inputs to the networks are obtained with the help of a factor model that further enables us to handle very high-dimensional predictors. In our test statistics, we employ high-quality estimators of the partial derivatives of the non-parametric regression function, which could be of independent interest. To demonstrate the versatility of our test, we apply it to assess the adequacy of non-parametric factor regression. An intriguing avenue for further exploration involves extending this framework to dependent data, other statistical machine learning losses, and simultaneous testing for multiple variables. Relevant examples for the latter direction include the $\ell_{\infty}$ statistics proposed by chen2022inference and the $\ell_2-\ell_{\infty}$ type statistics discussed by li2022ell.
This paper focuses on the popular feed-forward networks with the ReLU activation function. The derivative irregularity is not limited to the ReLU activation function, and hence, addressing this issue can benefit various neural network classes. Moreover, in our approach, we apply the smoothing procedure after we have obtained a neural network estimator. This is different from other derivative regularizations, for instance, Sobolev training Czarneck2017sobolev, that employ roughness penalties during optimization, which could form a potential future topic. Last but not least, we could also potentially conduct different tests for testing higher- or mixed-derivatives. For instance, testing the monotonicity in a variable in the regression is equivalent to testing the sign of the associated partial derivative, a task already studied in the context of an one-dimensional predictor by, for instance, BowmanJonesGijbels1998testing, GhosalArusharka2000testing and HallHuang2001nonparametric.
{ \spacingset{1.0} }