High-dimensional inference for single-index model with latent factors

The unknown smooth function $g(\cdot,\cdot)$ makes the estimation of the FASIM challenging. To address this concern, we reformulate the FASIM as a linear regression model with transfomed response. With this transformation, we estimate the parameters avoiding estimating the unknown function $g(\cdot,\cdot)$.

Without loss of generality, we assume that $\mathbb{E}(\bm{x})=\bm{0}$, and $\bm{\Sigma}=\mathrm{Cov}({\bm{x}})$ is a positive definite matrix. Let $\bm{\sigma}_{\bm{u}h}=\mathrm{Cov}\{\bm{u},h(Y)\}$ and $\bm{\sigma}_{\bm{f}h}=\mathrm{Cov}\{\bm{f},h(Y)\}$ for a given transformation function $h(\cdot)$ of the response. Define $\bm{\beta}_{h}=\mathbb{E}\left(\bm{u}\bm{u}^{\top}\right)^{-1}\bm{\sigma}_{{\bm{u}}h}$, $\bm{\gamma}_{h}=\mathbb{E}\left(\bm{f}\bm{f}^{\top}\right)^{-1}\bm{\sigma}_{\bm{f}h}$, $\bm{\eta}_{h}=(\bm{\beta}_{h}^{\top},\bm{\gamma}_{h}^{\top})^{\top}$, $\bm{\eta}=(\bm{\beta}^{\top},\bm{\gamma}^{\top})^{\top}$ and $\bm{v}=(\bm{u}^{\top},\bm{f}^{\top})^{\top}$. In the context of the SIM framework, the following linear expectation condition is commonly assumed.

Proposition 2.1 is from Theorem 2.1 of Li and Duan, (1989). The condition in this proposition is referred to as the linearity condition (LC) for predictors. It is satisfied when $\bm{v}$ follows an elliptical distribution, a commonly used assumption in the sufficient dimension reduction literature (Li,, 1991; Cook and Ni,, 2005). Hall and Li, (1993) showed that the LC holds approximately to many high-dimensional settings. Throughout the paper, we assume $\kappa_{h}\neq 0$, which is a relatively mild assumption. In fact, when $h(\cdot)$ is monotone and $g(\cdot,\cdot)$ is monotonic with respect to the first argument, this assumption is anatomically satisfied.

Following from the definition of $\bm{\beta}_{h}$ and $\bm{\gamma}_{h}$, our model can be rewritten as:

$\displaystyle h(Y)$

$\displaystyle=\bm{x}^{\top}\bm{\beta}_{h}+\bm{f}^{\top}\bm{\varphi}_{h}+{e}_{h},$

(2.1)

where $\bm{\varphi}_{h}=\bm{\gamma}_{h}-\bm{B}^{\top}\bm{\beta}_{h}$, and the error term $e_{h}$ satisfies $\mathbb{E}(e_{h})={{0}}$, $\mathbb{E}(e_{h}\bm{u})={\bm{0}_{p}}$ and $\mathbb{E}(e_{h}\bm{f})={\bm{0}_{K}}$. This implies that we now recast the FASIM as a FARM with transformed response. Actually for identification of the FASIM, it is usually assumed that $\|\bm{\eta}\|_{2}=1$ and its first element is positive. Thus the estimation of the direction of $\bm{\eta}$ is sufficient for the FASIM. The proportionality between $\bm{\eta}$ and $\bm{\eta}_{h}$, along with the previous mentioned transformed linear regression model (2.1), reduces the difficulty of analysis.

In practice, we need to choose an appropriate transformation function $h(\cdot)$. It’s worth mentioning that the procedure in Eftekhari et al., (2021) essentially works with $h(Y)=Y$. Motivated by Zhu and Zhu, (2009) and Rejchel and Bogdan, (2020), we set $h(Y)=F(Y)-1/2$, where $F(Y)$ is the distribution function of $Y$. This specific choice could make our procedures be robust against outliers and heavy tails. Actually with equation (2.1), given the widely imposed sub-Gaussian assumption on the predictors, the boundedness of $F(Y)-1/2$ would lead the transformed error term $e_{h}$ being sub-Gaussian, even if the original error term $\varepsilon$ comes from Cauchy distribution. The response-distribution transformation is preferred due to some additional reasons which will be discussed later.

Throughout the paper, we assume that the data $\{\bm{x}_{i},\bm{f}_{i},Y_{i},e_{hi}\}_{i=1}^{n}$ are independent and identically distributed (i.i.d.) copies of $\{\bm{x},\bm{f},Y,e_{h}\}$. Let $\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\top}$, $\bm{F}=(\bm{f}_{1},\ldots,\bm{f}_{n})^{\top}$, $\bm{Y}=(Y_{1},\ldots,Y_{n})^{\top}$ and $\bm{e}_{h}=({e}_{h1},\ldots,{e}_{hn})^{\top}$. We consider the high-dimensional scenario where the dimension $p$ of the observed covariate vector $\bm{x}$ can be much larger than the sample size $n$. For the FASIM, we first need to estimate the latent factor vector $\bm{f}$ since only the predictor vector ${\bm{x}}$ and the response $Y$ are observable. To address this issue, we impose an identifiability assumption similar to that in Fan et al., (2013). That is

$\displaystyle\mathrm{Cov}(\bm{f})=\bm{I}_{K},\ \ \text{and}\ \ \bm{B}^{\top}\bm{B}\ \ \text{is \ diagonal}.$

Consequently, the constrained least squares estimator of $(\bm{F},\bm{B})$ based on ${\bm{X}}$ is given as

	$\displaystyle(\widehat{\bm{F}},\widehat{\bm{B}})=$	$\displaystyle\arg\min\limits_{\bm{F}\in\mathbb{R}^{n\times K},~{}\bm{B}\in\mathbb{R}^{p\times K}}\|{\bm{X}}-\bm{F}\bm{B}^{\top}\|_{\mathbb{F}}^{2},$		(2.2)
	subjectto	$\displaystyle\ \frac{1}{n}{\bm{F}}^{\top}\bm{F}=\bm{I}_{K}\ \ \text{and}\ \ {\bm{B}}^{\top}\bm{B}\ \ \text{is}\ \ \text{diagonal}.$

Let $\bm{U}=(\bm{u}_{1},\ldots,\bm{u}_{n})^{\top}$, with $\bm{u}_{i}=\bm{x}_{i}-\bm{B}\bm{f}_{i},i=1,\ldots,n$. Denote $\bm{\Sigma}_{\bm{u}}=\mathrm{Cov}(\bm{u})$.

Elementary manipulation yields that the columns of $\widehat{\bm{F}}/\sqrt{n}$ are the eigenvectors corresponding to the largest $K$ eigenvalues of the matrix ${\bm{X}}{{\bm{X}}}^{\top}$ and $\widehat{\bm{B}}={\bm{X}}^{\top}{\widehat{\bm{F}}}({\widehat{\bm{F}}}^{\top}\widehat{\bm{F}})^{-1}=n^{-1}{\bm{X}}^{\top}{\widehat{\bm{F}}}$. Then the estimator of $\bm{U}$ is

$$\widehat{\bm{U}}={\bm{X}}-\widehat{\bm{F}}{\widehat{\bm{B}}}^{\top}=\left(\bm{I}_{n}-\frac{1}{n}\widehat{\bm{F}}{\widehat{\bm{F}}}^{\top}\right){\bm{X}},$$

see Fan et al., (2013). Since $K$ is related to the number of spiked eigenvalues of ${\bm{X}}{\bm{X}}^{\top}$, it is usually small. Therefore, we treat $K$ as a fixed constant as suggested by Fan et al., (2024). Additionally, let $\widehat{\bm{\Sigma}}$, $\widehat{\bm{\Lambda}}=\text{diag}\left(\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{K}\right)$ and $\widehat{\bm{\Gamma}}=\left(\widehat{\bm{\zeta}}_{1},\ldots,\widehat{\bm{\zeta}}_{K}\right)$ be suitable estimators of the covariance matrix $\bm{\Sigma}$ of $\bm{x}$, the matrix consisting of its leading $K$ eigenvalues $\bm{\Lambda}=\text{diag}(\lambda_{1},\ldots,\lambda_{K})$ and the matrix consisting of their corresponding $K$ orthonormalized eigenvectors $\bm{\Gamma}=(\bm{\zeta}_{1},\ldots,\bm{\zeta}_{K})$, respectively.

We proceed by presenting the regularity assumptions imposed in seminal works about factor analysis, such as Bai, (2003) and Fan et al., (2013, 2024).

We summarize the theoretical results related to consistent factor estimation in Lemma 1 in Supplementary Material, which directly follows from Proposition 2.1 in Fan et al., (2024) and Lemma 3.1 in Bayle and Fan, (2022).

In practice, the number of latent factors $K$ is often unknown, and determining $K$ in a data-driven way is a crucial challenge. Numerous methods have been introduced in the literature to estimate the value of $K$ (Bai and Ng,, 2002; Lam and Yao,, 2012; Ahn and Horenstein,, 2013; Fan et al.,, 2022). In this paper, we employ the ratio method for our numerical studies (Luo et al.,, 2009; Lam and Yao,, 2012; Ahn and Horenstein,, 2013). Let $\lambda_{k}({\bm{X}}{\bm{X}}^{\top})$ be the $k$-th largest eigenvalue of ${\bm{X}}{\bm{X}}^{\top}$. The number of factors can be consistently estimated by

$\displaystyle\widehat{K}=\arg\max\limits_{k\leq K_{\max}}\frac{\lambda_{k}({\bm{X}}{\bm{X}}^{\top})}{\lambda_{k+1}({\bm{X}}{\bm{X}}^{\top})},$

where $1\leq K_{\max}\leq n$ is a prescribed upper bound for $K$. In our subsequent theoretical analysis, we treat $K$ as known. All the theoretical results remain valid conditioning on that $\widehat{K}$ is a consistent estimator of $K$.

In this subsection, we develop a factor-adjusted score type test (FAST) and derive its Gaussian approximation result. For the adequacy testing problem, Fan et al., (2024) considered the maximum of the debiased lasso estimator of $\bm{\beta}$ under FARM. However, their procedure requires to estimate high-dimensional regression coefficients and also high-dimensional precision matrix, and thus is computationally expensive. While our proposed FAST does not need to estimate high-dimensional regression coefficients, nor high-dimensional precision matrix. Hence, Our procedure is straightforward to implement, saving both computation time and computational resources.

Under the null hypothesis in (3.1), we have

$$\mathbb{E}\left[\left\{F(Y)-{\frac{1}{2}}-\bm{f}^{\top}\bm{\gamma}_{h}\right\}\bm{u}\right]=\mathbb{E}[e_{h}\bm{u}]=\bm{0}.$$

While under alternative hypothesis $H_{1}$, we have

$$\mathbb{E}\left[\left\{F(Y)-{\frac{1}{2}}-\bm{f}^{\top}\bm{\gamma}_{h}\right\}\bm{u}\right]=\mathbb{E}\{\bm{u}(e_{h}+\bm{u}^{\top}\bm{\beta}_{h})\}=\bm{\Sigma}_{\bm{u}}\bm{\beta}_{h}\neq\bm{0}.$$

This observation motivates us to consider the following individual test statistic

$\displaystyle{T}_{nj}=\frac{1}{\sqrt{n}}\sum\limits_{i=1}^{n}\left[\left\{F_{n}(Y_{i})-\frac{1}{2}-\widehat{\bm{f}}_{i}^{\top}\widehat{\bm{\gamma}}_{h}\right\}\widehat{U}_{ij}\right].$

(3.2)

Here we utilize the empirical distribution $F_{n}(y)=n^{-1}\sum_{i=1}^{n}\mathbb{I}(Y_{i}\leq y)$ as an estimator of $F(Y)$. In the empirical distribution function, the term $\sum_{j=1}^{n}\mathbb{I}(Y_{j}\leq Y_{i})$ is the rank of $Y_{i}$. Since statistics with ranks such as Wilcoxon test and the Kruskall-Wallis ANOVA test, are well-known to be robust, this intuitively explains why our procedures with response-distribution transformation function would be robust with respect to outliers in response. We consider the least squares estimator $\widehat{\bm{\gamma}}_{h}=(\widehat{\bm{F}}^{\top}\widehat{\bm{F}})^{-1}\widehat{\bm{F}}^{\top}\left\{F_{n}(\bm{Y})-{1}/{2}\right\}$. Denote ${\bm{T}}_{n}=({T}_{n1},\ldots,{T}_{np})^{\top}$. To test the null hypothesis $H_{0}:\bm{\beta}_{h}=\bm{0}$, we consider $\ell_{\infty}$ norm of ${\bm{T}}_{n}$. That is,

$\displaystyle M_{n}=\|{\bm{T}}_{n}\|_{\infty}.$

(3.3)

It is clear that in the FAST statistic $M_{n}$ defined above, we only need to estimate a low-dimensional parameter $\bm{\gamma}_{h}$, and no high-dimensional parameters are required. In addition, since FAST does not rely on the estimation of the precision matrix, we can avoid assumptions about the $\ell_{\infty}$ norm of the precision matrix, as demonstrated in Fan et al., (2024). Further define $S_{nj}$ as follows:

$\displaystyle S_{nj}=\frac{1}{\sqrt{n}}\sum\limits_{i=1}^{n}\left[\left\{F(Y_{i})-\frac{1}{2}-\bm{f}_{i}^{\top}\bm{\gamma}_{h}\right\}U_{ij}+m_{j}(Y_{i})\right],$

(3.4)

where $m_{j}(y)=\mathbb{E}[(X_{1j}-\bm{f}_{1}^{\top}\bm{b}_{j})\{\mathbb{I}(Y\geq y)-F(Y)\}]$. Denote ${\bm{S}}_{n}=({S}_{n1},\ldots,{S}_{np})^{\top}$ and $\bm{\Omega}^{*}=\mathrm{Cov}({\bm{S}}_{n})$. Let $\sigma_{j}^{2}=\bm{\Omega}_{jj}^{*},\ j=1,\ldots,p$. It could be shown that ${T}_{nj}={S}_{nj}+o_{\mathbb{P}}(1)$.

Assumption 3 is mild and frequently employed in high-dimensional settings, which is the technical requirement to bound the difference between $M_{n}$ and $\|\bm{\mathrm{Z}}\|_{\infty}$, where $\bm{\mathrm{Z}}\sim\mathrm{N}_{p}(\bm{0},\bm{\Omega}^{*})$. Specifically, condition (i) imposes suitable restriction on the growth rate of $p$ and is commonly used in the high-dimensional inference literature, see Zhang and Cheng, (2017) and Dezeure et al., (2017). Condition (ii) imposes a boundedness restriction on the second moment of $S_{nj}$, which is also assumed in Ning and Liu, (2017).

Theorem 3.1 indicates that our test statistic can be approximated by the maximum of a high-dimensional Gaussian vector under mild conditions. Based on this, we can reject the null hypothesis $H_{0}$ at the significant level $\alpha$ if and only if $M_{n}>c_{1-\alpha}$, where $c_{1-\alpha}$ is the $(1-\alpha)$-th quantile of the distribution of $\|\bm{\mathrm{Z}}\|_{\infty}$.

As demonstrated by many authors, such as Cai et al., (2014) and Ma et al., (2021), the distribution of $\|\bm{\mathrm{Z}}\|_{\infty}$ can be asymptotically approximated by the Gumbel distribution. However, this generally requires restrictions on the covariance matrix of $\bm{\mathrm{Z}}$. For instance, it requires that the eigenvalues of $\bm{\Omega}^{*}$ are uniformly bounded. In addition, the critical value obtained from the Gumbel distribution may not work well in practice since this weak convergence is typically slow (Zhang and Cheng,, 2017). Instead of adopting Gumbel distribution, we consider using bootstrap to approximate the distribution of $M_{n}$.

Given that $\bm{\Omega}^{*}$ is unknown, we employ the Gaussian multiplier bootstrap to derive the critical value $c_{1-\alpha}$. The procedures and theoretical properties of the Gaussian multiplier bootstrap are outlined in this subsection.

• 1.

  Generate i.i.d. random variables $\mathrm{N}_{1},\ldots,\mathrm{N}_{n}\sim\mathrm{N}(0,1)$ independent of the observed dataset $\mathcal{D}=\{\bm{x}_{1},\ldots,\bm{x}_{n},Y_{1},\ldots,Y_{n}\}$, and compute

  $\displaystyle\widehat{\mathcal{G}}=\max\limits_{j\in[p]}\left|\frac{1}{\sqrt{n}}\sum\limits_{i=1}^{n}\left[\left\{F_{n}(Y_{i})-\frac{1}{2}-\widehat{\bm{f}}_{i}^{\top}\widehat{\bm{\gamma}}_{h}\right\}\widehat{U}_{ij}+\widehat{m}_{j}(Y_{i})\right]\mathrm{N}_{i}\right|.$

  (3.6)

  Here $\widehat{m}_{j}(y)={n}^{-1}\sum_{i=1}^{n}\left(X_{ij}-\widehat{\bm{f}}_{i}^{\top}\widehat{\bm{b}}_{j}\right)\{\mathbb{I}(Y_{i}\geq y)-F_{n}(Y_{i})\}$.

• 2.

  Repeat the first step independently for $B$ times to obtain $\widehat{\mathcal{G}}_{1},\ldots,\widehat{\mathcal{G}}_{B}$. Approximate the critical value $c_{1-\alpha}$ via the $(1-\alpha)$-th quantile of the empirical distribution of the bootstrap statistics:

  $\displaystyle\widehat{c}_{1-\alpha}=\inf\left\{t\geq 0:\frac{1}{B}\sum\limits_{b=1}^{B}\mathbb{I}(\widehat{\mathcal{G}}_{b}\leq t)\geq 1-\alpha\right\}.$

  (3.7)

• 3.

  We reject the null hypothesis $H_{0}$ if and only if

  $\displaystyle M_{n}\geq\widehat{c}_{1-\alpha}.$

  (3.8)

Theorem 3.2 demonstrates the validity of the proposed bootstrap procedure, grounded in the Gaussian approximation theory established in Theorem 3.1. Based on this result, we can define our decision rule as follows:

$\displaystyle\psi_{\infty,\alpha}=\mathbb{I}\left(M_{n}>\widehat{c}_{1-\alpha}\right).$

The null hypothesis $H_{0}$ is rejected if and only if $\psi_{\infty,\alpha}=1$.

We next consider the asymptotic power analysis of the $M_{n}$. To demonstrate the efficiency of the test statistic, we consider the following local alternative parameter family for $\bm{\beta}_{h}$ under $H_{1}$.

$\displaystyle\bm{\Delta}(c)=\left\{\bm{\beta}_{h}\in\mathbb{R}^{p}:\max\limits_{j\in[p]}\left|{\beta_{hj}}\right|\geq\sqrt{c\frac{\log p}{n}}\right\},$

(3.9)

where $c$ is a positive constant. The $s=\|\bm{\beta}_{h}\|_{0}$ quantifies the sparsity of the parameter $\bm{\beta}_{h}$. Denote $\mu_{j}=\mathbb{E}(U_{ij}^{2})$.

Assumption 4 outlines the diagonal structure of the covariance matrix $\bm{\Sigma}_{\bm{u}}$, a common assumption in factor analysis (Kim and Mueller,, 1978). We should also note that this diagonality of $\bm{\Sigma}_{\bm{u}}$ is only imposed here for simplify the power analysis.

Theorem 3.3 suggests that our test procedure maintains high power even when only a few components of $\bm{\beta}_{h}$ have magnitudes larger than $\sqrt{(2+\varrho_{0})(\log p)/{n}}$. Thus, our testing procedure is powerful against sparse alternatives. Notably, this separation rate represents the minimax optimal rate for local alternative hypotheses, as discussed in Verzelen, (2012), Cai et al., (2014), Zhang and Cheng, (2017), and Ma et al., (2021).

In the high-dimensional regime, we employ $\ell_{1}$ penalty (Tibshirani,, 1996) to estimate the unknown parameter vectors $\bm{\beta}_{h}$ and $\bm{\gamma}_{h}$ in model (2.1):

$\displaystyle(\widehat{\bm{\beta}}_{h},\widehat{\bm{\gamma}}_{h})=\arg\min\limits_{\bm{\beta}\in\mathbb{R}^{p},\bm{\gamma}\in\mathbb{R}^{K}}\left[\frac{1}{2n}\sum\limits_{i=1}^{n}\left\{F_{n}(Y_{i})-\frac{1}{2}-\widehat{\bm{u}}_{i}^{\top}\bm{\beta}-\widehat{\bm{f}}_{i}^{\top}\bm{\gamma}\right\}^{2}+\lambda\|\bm{\beta}\|_{1}\right],$

(4.1)

where $\lambda>0$ is a tuning parameter.

By utilizing pseudo response observations $F_{n}(Y_{i})$ instead of $Y_{i}$, our approach is inherently robust when confronting with outliers. We note that for the SIM without latent factors, the distribution function transformation has also been considered by many other authors, such as Zhu and Zhu, (2009), Wang et al., (2012), and Rejchel and Bogdan, (2020). However, in our model $\bm{f}$ and $\bm{u}$ are unobserved and must be estimated firstly. Thus we require additional efforts to derive the theoretical properties of $\widehat{\bm{\beta}}_{h}$.

Let $\widetilde{F}_{n}(\bm{Y})=(\bm{I}_{n}-\widehat{\bm{P}})\left\{F_{n}(\bm{Y})-1/2\right\}$ represent the residuals of the response vector $F_{n}(\bm{Y})-1/2$ after it has been projected onto the column space of $\widehat{\bm{F}}$, where $\widehat{\bm{P}}=n^{-1}\widehat{\bm{F}}\widehat{\bm{F}}^{\top}$ is the corresponding projection matrix. Since $\widehat{\bm{U}}=(\bm{I}_{n}-\widehat{\bm{P}})\bm{X}$, $\widehat{\bm{F}}^{\top}\widehat{\bm{U}}={\bm{0}_{K\times p}}$. Then direct calculations yield that the solution of (4.1) is equivalent to

	$\displaystyle\widehat{\bm{\beta}}_{h}$	$\displaystyle=\arg\min\limits_{\bm{\beta}\in\mathbb{R}^{p}}\left\{\frac{1}{2n}\left\|\widetilde{F}_{n}(\bm{Y})-\widehat{\bm{U}}\bm{\beta}\right\|_{2}^{2}+\lambda\|\bm{\beta}\|_{1}\right\},$		(4.2)
	$\displaystyle\widehat{\bm{\gamma}}_{h}$	$\displaystyle=(\widehat{\bm{F}}^{\top}\widehat{\bm{F}})^{-1}\widehat{\bm{F}}^{\top}\left\{F_{n}(\bm{Y})-\frac{1}{2}\right\}.$		(4.3)

For any vector $\widehat{\bm{\beta}}_{h}$ defined in (4.2), we have the following result.

Based on the empirical distribution function of the response, we establish the $\ell_{1}$ and $\ell_{2}$-error bounds for our introduced estimator $\widehat{\bm{\beta}}_{h}$ in Theorem 4.1. Compared with Zhu and Zhu, (2009), Wang et al., (2012), and Rejchel and Bogdan, (2020), in our model $\bm{f}$ and $\bm{u}$ are unobserved and must be estimated firstly. This adds new technical difficulty. Further compared with Fan et al., (2024), our procedures do not need any moment condition for the error term $\varepsilon$ in the model.

Due to the bias inherent in the penalized estimation, $\widehat{\bm{\beta}}_{h}$ is unsuitable for direct utilization in statistical inference. To overcome this obstacle, Zhang and Zhang, (2014) and Van de Geer et al., (2014) proposed the debiasing technique for Lasso estimator in linear regression. Eftekhari et al., (2021) extended the technique to SIM. Han et al., (2023) and Yang et al., (2024) further considered robust inference for high-dimensional SIM. However, their procedures cannot handle latent factors. To this end, we introduce debiased estimator for the FASIM. Motivated by Zhang and Zhang, (2014) and Van de Geer et al., (2014), we construct the following debiased estimator

$\displaystyle\widetilde{\bm{\beta}}_{h}=\widehat{\bm{\beta}}_{h}+\frac{1}{n}\widehat{\bm{\Theta}}\widehat{\bm{U}}^{\top}\left\{F_{n}(\bm{Y})-\frac{1}{2}-\widehat{\bm{U}}\widehat{\bm{\beta}}_{h}\right\},$

(4.4)

where $\widehat{\bm{\Theta}}$ is a sparse estimator of $\bm{\Sigma}_{u}^{-1}$, with $\bm{\Sigma}_{u}=\mathbb{E}\left({\bm{u}}_{i}{\bm{u}}_{i}^{\top}\right)$. Denote $\widehat{\bm{\Sigma}}_{u}=n^{-1}\sum_{i=1}^{n}\widehat{\bm{u}}_{i}\widehat{\bm{u}}_{i}^{\top}$. Then we construct ${\widehat{\bm{\Theta}}}$ by employing the approach in Cai et al., (2011). Concretely, ${\widehat{\bm{\Theta}}}$ is the solution to the following constrained optimization problem:

	$\displaystyle\widehat{\bm{\Theta}}$	$\displaystyle=\arg\min\limits_{\bm{\Theta}\in\mathbb{R}^{p\times p}}\|\bm{\Theta}\|_{\text{sum}},$
		$\displaystyle\text{s.t.}\ \|\bm{\Theta}\widehat{\bm{\Sigma}}_{u}-\bm{I}_{p}\|_{\max}\leq\delta_{n},$		(4.5)

where $\delta_{n}$ is a predetermined tuning parameter. In general, $\widehat{\bm{\Theta}}$ is not symmetric since there is no symmetry constraint in (4.2). Symmetry can be enforced through additional operations. Denote $\widehat{\bm{\Theta}}=\left(\widehat{\bm{\xi}}_{1},\ldots,\widehat{\bm{\xi}}_{p}\right)^{\top}=(\widehat{\xi}_{ij})_{1\leq i,j\leq p}$. Write $\widehat{\bm{\Theta}}^{\text{sym}}=\left(\widehat{\xi}_{ij}^{\text{sym}}\right)_{1\leq i,j\leq p}$, where $\widehat{\xi}_{ij}^{\text{sym}}$ is defined as:

$\displaystyle\widehat{\xi}_{ij}^{\text{sym}}=\widehat{\xi}_{ij}\mathbb{I}\left(|\widehat{\xi}_{ij}|\leq|\widehat{\xi}_{ji}|\right)+\widehat{\xi}_{ji}\mathbb{I}\left(|\widehat{\xi}_{ij}|>|\widehat{\xi}_{ji}|\right).$

Apparently, $\widehat{\bm{\Theta}}^{\text{sym}}$ is a symmetric matrix. For simplicity, we write $\widehat{\bm{\Theta}}$ as the symmetric estimator in the rest of the paper. Next, we consider the estimation error bound of $\widehat{\bm{\Theta}}$. To achieve this target, we need to introduce the following assumption on the inverse of the $\bm{\Sigma}_{u}$.

Assumption 5 necessitates that $\bm{\Sigma}_{u}^{-1}$ be sparse in terms of both its $\ell_{\infty}$-norm and matrix row space. Van de Geer et al., (2014), Cai et al., (2016) and Ning and Liu, (2017) also discussed similar assumptions on precision matrix estimation and more general inverse Hessian matrix estimation.

The estimation error bound of $\widehat{\bm{\Theta}}$ and the upper bound of $\left\|\widehat{\bm{\Theta}}-\bm{\Sigma}_{u}^{-1}\right\|_{\infty}$ are shown in Proposition 7 and Lemma 8 in the Supplementary Material, which are crucial for establishing theoretical results afterwards. Denote $\bm{m}=\{m_{1}(Y),\ldots,m_{p}(Y)\}^{\top}$, with $m_{j}(y)=\mathbb{E}[(X_{1j}-\bm{f}_{1}^{\top}\bm{b}_{j})\{\mathbb{I}(Y\geq y)-F(Y)\}]$. Let $\bm{m}_{i}$ be the i.i.d copies of $\bm{m}$. Namely, $\bm{m}_{i}=\{m_{1}(Y_{i}),\ldots,m_{p}(Y_{i})\}^{\top}$.

Theorem 4.2 indicates that the asymptotic representation of $\sqrt{n}(\widetilde{\bm{\beta}}_{h}-{\bm{\beta}}_{h})$ can be divided into two terms. The major term $\bm{Z}_{b}$ is associated with the error $\bm{e}_{h}$ and $\bm{m}$, while the remainder $\bm{E}$ vanishes as $n$ approaches to infinity.

Further for $j=1,\ldots,p$, the variance of ${Z}_{bj}$ is ${\sigma}_{zj}^{2}=\bm{e}_{j}^{\top}\bm{\Sigma}_{u}^{-1}\mbox{Var}\{\bm{u}e_{h}+\bm{m}\}\bm{\Sigma}_{u}^{-1}\bm{e}_{j}$, with $\bm{e}_{j}$ being an unit vector with only its $j$-th element being 1. Based on Theorem 4.2, the asymptotic variance of $\sqrt{n}(\widetilde{{\beta}}_{hj}-{{\beta}}_{hj})$ can be estimated by $\widehat{\sigma}_{zj}^{2}=\bm{e}_{j}^{\top}\widehat{\bm{\Theta}}\widetilde{\bm{\Sigma}}\widehat{\bm{\Theta}}\bm{e}_{j}$, with $\widetilde{\bm{\Sigma}}=(\widetilde{\sigma}_{lk})_{1\leq l,k\leq p}$, $\widetilde{\sigma}_{lk}={n}^{-1}\sum_{i=1}^{n}\left\{\widehat{U}_{il}\widetilde{e}_{hi}+\widehat{m}_{l}(Y_{i})\right\}\left\{\widehat{U}_{ik}\widetilde{e}_{hi}+\widehat{m}_{k}(Y_{i})\right\}$, $\widetilde{e}_{hi}=F_{n}(Y_{i})-{1}/{2}-\widehat{\bm{u}}_{i}^{\top}\widehat{\bm{\beta}}_{h}-\widehat{\bm{f}}_{i}^{\top}\widehat{\bm{\gamma}}_{h}$ and $\widehat{m}_{j}(y)={n}^{-1}\sum_{i=1}^{n}\left(X_{ij}-\widehat{\bm{f}}_{i}^{\top}\widehat{\bm{b}}_{j}\right)\{\mathbb{I}(Y_{i}\geq y)-F_{n}(Y_{i})\}$. Thus the $(1-\alpha)$ confidence interval for $\beta_{hj}$, $j\in[p]$ can be constructed as follows

$\displaystyle\mathcal{CI}_{\alpha}(\beta_{hj})=\left(\widetilde{\beta}_{hj}-\frac{\widehat{\sigma}_{zj}z_{1-\alpha/2}}{\sqrt{n}},\ \widetilde{\beta}_{hj}+\frac{\widehat{\sigma}_{zj}z_{1-\alpha/2}}{\sqrt{n}}\right),$

(4.6)

where $\widetilde{\beta}_{hj}$ is the $j$-th component of $\widetilde{\bm{\beta}}_{h}$ and $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-th quantile of standard normal distribution.

In this subsection, we set $n=200$, $K=2$, $p=200$ or $500$, and ${\bm{X}}=\bm{F}\bm{B}^{\top}+\bm{U}$. Here, the entries of $\bm{B}$ are generated from the uniform distribution $\mathrm{Unif}(-1,1)$ and every row of $\bm{U}$ follows from $\mathrm{N}(\bm{0},\bm{\Sigma}_{\bm{u}})$ with $\bm{\Sigma}_{\bm{u}}=(\sigma_{uij})_{1\leq i,j\leq p}$ and $\sigma_{uij}=0.5^{|i-j|}$. We consider two cases of $\bm{F}$ generation:
Case i. We generate $\bm{f}_{i}$ from standard normal distribution, that is,

$\displaystyle\bm{f}_{i}\sim\mathrm{N}(\bm{0},\bm{I}_{K}),\ i=1,\ldots,n.$

(5.3)

Case ii. We generate $\bm{f}_{i}$ from the AR(1) model:

$\displaystyle\bm{f}_{i}=\bm{\Phi}\bm{f}_{i-1}+\bm{\pi}_{i},\,i=2,\ldots,n,$

(5.4)

where $\bm{\Phi}=(\phi_{ij})\in\mathbb{R}^{K\times K}$ with ${\phi}_{ij}=0.4^{|i-j|+1},\ i,j\in[K]$. In addition, ${\bm{f}_{1}}$ and ${\bm{\pi}_{i}},\ i\geq 1$ are independently drawn from $\mathrm{N}(\bm{0},\bm{I}_{K})$. We generate ${\varepsilon}$ either from (a) $\mathrm{N}(0,0.25)$ or (b) Student’s t distribution with 3 degree of freedom, denoted as $\mathrm{t}_{3}$. We set $\gamma=(0.5,0.5)^{\top}$ and $\bm{\beta}=\omega\ast\left(\mathbf{1}_{3},\mathbf{0}_{p-3}\right)^{\top}$, $\omega\geq 0$. When $\omega=0$, it indicates that the null hypothesis holds, and the simulation results correspond to the empirical size. Otherwise, they correspond to the empirical power.