# Debiased and threshold ridge regression for linear model with
heteroskedastic and dependent error

Yunyi Zhang and Dimitris N. Politis

###### Abstract

Focusing on a high dimensional linear model $y=X\beta+\epsilon$ with
dependent, non-stationary, and heteroskedastic errors, this paper applies the
debiased and threshold ridge regression method that gives a consistent
estimator for linear combinations of $\beta$; and derives a Gaussian
approximation theorem for the estimator. Besides, it proposes a dependent wild
bootstrap algorithm to construct the estimator’s confidence intervals and
perform hypothesis testing. Numerical experiments on the proposed estimator
and the bootstrap algorithm show that they have favorable finite sample
performance.

Research on a high dimensional linear model with dependent(non-stationary)
errors is sparse, and our work should bring some new insights to this field.

##  1 Introduction

Linear regression with independent and identically distributed(i.i.d.) errors
is a fundamental topic in statistical inference. The classical setting assumes
the dimension of parameters in a linear model is constant. Under this setting,
research has been proposed on estimation, e.g., Zou and Hastie [1] and Zou
[2]; confidence intervals construction/hypothesis testing, e.g., Chatterjee
and Lahiri [3, 4]; and prediction, e.g., Stine [5] and Zhang and Politis [6].
We also refer Seber and Lee [7] for an extensive introduction.

In reality, however, errors in a linear model can be dependent or have
different distributions. For example,-as Vogelsang [8] and Petersen [9]
suggested-, heteroscedasticity, autocorrelation and spatial correlation
present in panel data. In this case confidence intervals developed for i.i.d.
errors may fail to capture the correct probability. New tools are developed to
adapt to the non-i.i.d. errors. Andrews [10] and Kim and Sun [11] considered
estimating the ordinary least square estimator’s covariance matrix; Kelejian
and Prucha [12] and Vogelsang [8] proposed test statistics; Sun and Wang [13]
and Conley et al. [14] worked on inference and hypothesis testing, etc.
Despite great success, their works focus on the classical situation, i.e., the
dimension of parameters is considered to be fixed and does not change as the
sample size increases.

In the modern era, observations may have a comparable or even larger dimension
than the number of samples. So the dimension of parameters cannot be
considered as fixed. In order to perform consistent estimation, statisticians
need to assume the underlying parameters are sparse(i.e., the parameters
contain many zeros), and proceed with statistical inference based on this
assumption. Lasso is a suitable algorithm for this setting since it conducts
an implicit model selection, i.e., zeroing out parameters that are not
significant, see Tibshirani [15]. More recent work includes Zhao and Yu [16],
Meinshausen and Bühlmann [17] and Meinshausen and Yu [18] for model selection;
Zhang and Zhang [19], Zhang and Cheng [20] and Chatterjee and Lahiri [4, 3]
for statistical inference and hypothesis testing; Greenshtein and Ritov [21]
for prediction and Zou [2] for algorithm improvement. We refer Bühlmann and
van de Geer [22] for a comprehensive overview of the Lasso method on high
dimensional data set.

Lasso is not the only choice for fitting a high-dimensional linear model. Fan
and Li [23] introduced a new penalty function, called SCAD, that is
continuously differentiable and maintains the sparsity of the underlying
model. Lee et al. [24], Liu and Yu [25] and Tibshirani et al. [26] considered
Post-selection inference, i.e., performing model selection with Lasso, then
fitting ordinary least square regression on the selected parameters. Shao and
Deng [27] applied thresholding on ridge regression to recover the sparsity of
a linear model; Zhang and Politis [28] made a further improvement through
debiasing and thresholding. After debiasing and thresholding, the ridge
regression estimator had a comparable performance to the complex models like
threshold Lasso or post-selection inference. Moreover, the estimator had a
closed-form formula, making it easy to derive theoretical guarantees.

This paper focuses on fitting a high dimensional linear model, constructing
confidence intervals, and performing hypothesis testing, with the presence of
dependent and heteroskedastic errors. Since non-i.i.d. errors may appear in
fixed dimensional linear models, they may also appear in high dimensional
linear models. Unfortunately, performing statistical inference for a high
dimensional linear model with dependent and heteroskedastic errors is
challenging. Wu and Wu [29] proposed an oracle inequality; Han and Tsay [30]
proved the consistency of the Lasso estimator. However, their works relied on
a stationary assumption, e.g., definition 1.3.3 in Brockwell and Davis [31].

Suppose a high dimensional sparse linear model $y=X\beta+\epsilon$ with
dependent, non-stationary, heteroskedastic errors
$\epsilon=(\epsilon_{1},...,\epsilon_{n})^{T}$. Here $\epsilon_{i},i=1,...,n$
are not necessary to be linear processes. This paper focuses on estimating
linear combinations of parameters $\gamma=M\beta$ with $M$ a given matrix. We
adopt the debiased and threshold ridge regression estimator proposed in Zhang
and Politis [28]. After selecting a suitable ridge parameter and a threshold,
this method has a comparable performance to Lasso and is easily analyzed.

Our work includes proving the model selection consistency as well as the
consistency of the estimator; deriving a Gaussian approximation theorem for
the estimator and constructing a confidence interval for $\gamma$. We are also
interested in testing the statistical hypothesis

$\text{null: }M\beta=z\ \text{versus the alternative: }M\beta\neq z$ (1)

with $z$ a given vector. To achieve this goal, we adopt a dependent wild
bootstrap(Shao [32]) and provide its theoretical guarantee. Since there is
little research on statistical inference for a high dimensional linear model
with non-i.i.d. errors, this paper should shed some light on this field.

The remainder of this paper is organized as follows: section 2 introduces
frequently used notations and assumptions. Section 3 presents the consistency
results and the Gaussian approximation theorem for the proposed estimator.
Section 4 constructs a confidence interval for $\gamma=M\beta$, and tests the
null hypothesis $\gamma=z$ versus the alternative hypothesis $\gamma\neq z$
via dependent wild bootstrap. Section 5 provides numerical experiments, and
section 6 makes conclusions. Technical proofs are deferred to the Appendix.

##  2 Preliminarily

Our work considers the fixed design linear model

$y=X\beta+\epsilon$ (2)

where the unknown parameter vector
$\beta=(\beta_{1},...,\beta_{p})^{T}\in\mathbf{R}^{p}$, and the $n\times
p(p<n)$ fixed(nonrandom) design matrix $X=(x_{ij})_{i=1,...,n,j=1,...,p}$ is
assumed to have rank $p$. Denote $y=(y_{1},...,y_{n})^{T}$ and the errors
$\epsilon=(\epsilon_{1},...,\epsilon_{n})^{T}$. From the thin singular value
decomposition(theorem 7.3.2 in [33]), $X=P\Lambda Q^{T}$ with
$P=(p_{ij})_{i=1,2,...,n,j=1,2,...,p},Q=(q_{ij})_{i,j=1,2,...,p}$ respectively
being the $n\times p,p\times p$ orthonormal matrix, i.e.,
$P^{T}P=Q^{T}Q=QQ^{T}=I_{p}$. Here $I_{p}$ denotes the $p\times p$ identity
matrix. And $\Lambda=diag(\lambda_{1},...,\lambda_{p})$,
$\lambda_{i},i=1,2,...,p$ are singular values of $X$. Define the set
$\mathcal{N}_{b}=\left\\{i=1,2,...,p:|\beta_{i}|>b\right\\},\forall b>0$. We
are interested in constructing the confidence region for $\gamma=M\beta$, and
testing the statistical hypothesis (1). Here
$M=(m_{ij})_{i=1,2,...,p_{1},j=1,2,...,p}$ is a given $p_{1}\times p$ linear
combination matrix, and $z=(z_{1},...,z_{p_{1}})^{T}$ is the given expected
value of $M\beta$.

After choosing a suitable threshold $b$, define
$c_{ij},i=1,2,...,p_{1},j=1,2,...,p$ as
$c_{ij}=\sum_{k\in\mathcal{N}_{b}}m_{ik}q_{kj}$, and
$\mathcal{M}=\left\\{i=1,2,...,p_{1}:\sum_{j=1}^{p}c^{2}_{ij}>0\right\\}$.
Define

$\tau_{i}=\sqrt{\frac{1}{n}+\sum_{j=1}^{p}c^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}},i=1,2,...,p_{1}$.
Define $\Sigma=(\sigma_{ij})_{i,j=1,2,...,n}$ with
$\sigma_{ij}=\mathbf{E}\epsilon_{i}\epsilon_{j}$ as the covariance matrix of
errors.

To quantify the dependency among random variables, Wu [34](also see Wu and Wu
[35]) introduced the concept ‘physical dependence’. However, this concept was
designed for stationary random variables. We extend this concept to adapt to
non-stationary errors $\epsilon$: suppose $e_{i},i=...,-1,0,1,...$ are
independent(not necessarily identically distributed) random variables, and
$\epsilon_{i}=g_{i}(...,e_{i-1},e_{i})$, here $g_{i}$ are measurable
functions, $i=1,2,...,n$. Define $\mathcal{F}_{i}$ as the $\sigma$-field
generated by $...,e_{i-1},e_{i}$, so $\epsilon_{i}$ is $\mathcal{F}_{i}$
measurable. Suppose independent random variables
$e_{i}^{\dagger},i=...,-1,0,1,...$ are independent with
$e_{j},j=...,-1,0,1,...$; and $e_{i}^{\dagger}$ has the same distribution as
$e_{i}$, $\forall i$. Define
$\epsilon_{i,j}=g_{i}(...,e_{i-j-2},e_{i-j-1},e_{i-j}^{\dagger},e_{i-j+1},...,e_{i})$;
and the filter $\mathcal{F}_{i,j}$ as the $\sigma$-field generated by
$e_{i-j},e_{i-j+1},...,e_{i}$. Here $i\in\mathbf{Z},j\geq 0$. For a constant
$m\geq 1$, define the norm $\|\cdot\|_{m}=(\mathbf{E}|\cdot|^{m})^{1/m}$. Then
define $\delta_{i,j,m}=\|\epsilon_{i}-\epsilon_{i,j}\|_{m}$.

This paper applies the standard order notations $O(.),o(.),O_{p}(.)$ and
$o_{p}(.)$: for two numerical sequence $a_{n},b_{n}$, we say $a_{n}=O(b_{n})$
if $\exists$ a constant $C>0$ such that $|a_{n}|\leq C|b_{n}|$ for all $n$;
and $a_{n}=o(b_{n})$ if $\lim_{n\to\infty}a_{n}/b_{n}=0$. For two random
variable sequences $X_{n},Y_{n}$, we say $X_{n}=O_{p}(Y_{n})$ if for any
$0<\varepsilon<1$, $\exists$ a constant $C_{\varepsilon}$ such that
$\sup_{n}Prob(|X_{n}|\geq C_{\varepsilon}|Y_{n}|)\leq\varepsilon$; and
$X_{n}=o_{p}(Y_{n})$ if $X_{n}/Y_{n}\to_{p}0$. See definition 1.9 and chapter
1.5.1 of Shao [36]. All order notations and convergences in this paper are
understood to hold as the sample size $n\to\infty$.

For a set $A$, we use $|A|$ to indicate the number of elements in $A$. For a
vector $a=(a_{1},...,a_{n})^{T}$, define its $q$ norm as
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|a\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{q}=(\sum_{i=1}^{n}|a_{i}|^{q})^{1/q}$.
For a matrix $T$, define(and mainly use) the operator norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|T\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}=\max_{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|a\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{2}=1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Ta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$.
The notation $\to$ and $\to_{p}$ respectively indicates convergence in
$\mathbf{R}$, and convergence in probability. $\exists$ and $\forall$
respectively means ‘there exists’ and ‘for all’. Define
$Prob^{*}(\cdot)=Prob(\cdot|y)$ and $\mathbf{E}^{*}\cdot=\mathbf{E}(\cdot|y)$,
i.e., the probability and the expectation in the bootstrap world.

This paper uses the following assumptions:

Assumptions

1\. The fixed design matrix $X$ has rank $p\leq n$. There exists constants
$c_{\lambda},\ C_{\lambda}>0,1/2\geq\eta>0$ such that

$C_{\lambda}n^{1/2}\geq\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{p}\geq
c_{\lambda}n^{\eta},\forall n$ (3)

Besides, $\max_{i=1,...,n,j=1,...,p}|x_{ij}|=O(1)$. $p_{1}=O(1)$. Here $p_{1}$
is the number of simultaneous linear combinations in (1).

2\. $\mathbf{E}\epsilon_{i}=0,i=1,2,...,n$. There exists constants
$m>\frac{3}{\eta}$ and $\alpha_{\epsilon}>1$ such that

$\sup_{k=0,1,...}(k+1)^{\alpha_{\epsilon}}\sum_{j=k}^{\infty}\max_{i=1,2,...,n}\delta_{i,j,m}=O(1)\
\text{and }\max_{i=1,2,...,n}\|\epsilon_{i}\|_{m}=O(1)$ (4)

Besides, $\exists$ a constant $c_{\Sigma}>0$ such that $\Sigma$’s minimum
eigenvalue is greater than $c_{\Sigma}$ for $\forall n$; here
$\Sigma=\mathbf{E}\epsilon\epsilon^{T}$.

3\.
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}=O(n^{\alpha_{\beta}})$
with a constant $\alpha_{\beta}$ such that $0\leq\alpha_{\beta}<3\eta$, and
$|\mathcal{N}_{b_{n}}|=o(n^{\eta-1/m})$.

4\. $\rho_{n}=O(n^{2\eta-\delta})$ with a constant $\delta$ such that
$\frac{\eta+\alpha_{\beta}}{2}<\delta<2\eta$. Besides, there exists constants
$C_{b}>0$ and $0<\nu_{b}<\eta-\frac{3}{m}$ such that $b_{n}=C_{b}n^{-\nu_{b}}$
for $\forall n$. Assume $\exists$ a constant $0<c_{b}<1$ such that
$\max_{i\not\in\mathcal{N}_{b_{n}}}|\beta_{i}|\leq c_{b}\times b_{n}$ and
$\min_{i\in\mathcal{N}_{b_{n}}}|\beta_{i}|\geq b_{n}/c_{b}$.

5\. $\mathcal{M}$ is not empty. And $\exists$ constants $0<c_{\mathcal{M}}\leq
C_{\mathcal{M}}<\infty$ such that
$c_{\mathcal{M}}\leq\sum_{k\in\mathcal{N}_{b_{n}}}m^{2}_{ik}(=\sum_{j=1}^{p}c^{2}_{ij})\leq
C_{\mathcal{M}}$ for $\forall i\in\mathcal{M}$ and $\forall n$. Also assume

$\max_{i=1,...,p_{1}}|\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}|=o(1/\sqrt{n})\text{
and
}\max_{i=1,...,n}|\sum_{j\not\in\mathcal{N}_{b_{n}}}x_{ij}\beta_{j}|=o(n^{\frac{1/m-\eta}{2}})$
(5)

6\. The matrix $(c_{ij})_{i\in\mathcal{M},j=1,2,...,p}$ has rank
$|\mathcal{M}|$ and

$c^{*}=\max_{i\in\mathcal{M},l=1,2,...,n}\frac{1}{\tau_{i}}|\sum_{j=1}^{p}c_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)|=o(n^{-1/4}\log^{-z}(n))$
(6)

Here $z=\max(\frac{9}{2},\frac{3\alpha_{\epsilon}}{2\alpha_{\epsilon}-2})$.

7\. Suppose a function $K:\mathbf{R}\to[0,\infty)$ is symmetric, continuously
differentiable, $K(0)=1$, $\int_{\mathbf{R}}K(x)dx<\infty$, and $K$ is
decreasing on $[0,\infty)$. Define the Fourier transformation of $K$ as
$\mathcal{F}K(x)=\int_{\mathbf{R}}K(t)\exp(-2\pi\rm{i}tx)dt$. Assume
$\mathcal{F}K(x)\geq 0,\forall x\in\mathbf{R}$ and
$\int_{\mathbf{R}}\mathcal{F}K(x)dx<\infty$. Define a bandwidth parameter
$k_{n}$ satisfying $\lim_{n\to\infty}k_{n}=\infty$,
$k_{n}=o(n^{2\eta-2\nu_{b}-4/m})$, and
$k_{n}=O(n^{\frac{\eta}{2}-\frac{3}{2m}})$.

In assumption 6, $(c_{ij})_{i\in\mathcal{M},j=1,2,...,p}$ denotes the sub-
matrix formed by selecting $\mathcal{M}$ rows from the matrix
$(c_{ij})_{i=1,...,p_{1},j=1,...,p}$.

Our work mainly focuses on fixed design, i.e., no randomness involves in the
design matrix $X$. In the case of random design, results in this paper hold
true after conditioning on the design matrix $X$, i.e., replacing
$Prob(\cdot)$ by $Prob(\cdot|X)$; $\mathbf{E}\cdot$ by $\mathbf{E}\cdot|X$;
$Prob^{*}(\cdot)$ by $Prob(\cdot|X,y)$ and $\mathbf{E}^{*}\cdot$ by
$\mathbf{E}\cdot|X,y$.

Example 1 introduces a situation in which assumption 1 is satisfied.

###### Example 1

Suppose $n>p$ and $p/n\to c\in(0,1)$. Choose
$X=(x_{ij})_{i=1,...,n,j=1,...,p}$ such that each $x_{ij}$ is a realization of
i.i.d. random variables with mean $0$, variance $1$, and finite fourth order
moment. According to Bai and Yin [37], the smallest eigenvalue of
$\frac{1}{n}X^{T}X$ converges to $(1-\sqrt{c})^{2}$ almost surely as
$n\to\infty$. So the smallest singular value of $X$, being the square root of
$X^{T}X$’s smallest eigenvalue, is greater than $\frac{1-\sqrt{c}}{2}\sqrt{n}$
for sufficiently large $n$ almost surely. On the other hand, the largest
eigenvalue of $\frac{1}{n}X^{T}X$ converges to $(1+\sqrt{c})^{2}$ as
$n\to\infty$. So the largest singular value of $X$ has order $O(\sqrt{n})$
almost surely.

###### Remark 1

1\. Under stationary assumptions, e.g., Wu [34]

$\delta_{i,j,m}=\|\epsilon_{0}-\epsilon_{0,j}\|_{m}=\delta_{0,j,m}\Rightarrow\sum_{j=k}^{\infty}\max_{i=1,...,n}\delta_{i,j,m}=\sum_{j=k}^{\infty}\delta_{0,j,m}$
(7)

and (4) coincides with (2.8) in Wu and Wu [35], i.e., the dependence adjusted
norm condition. Therefore, assumption 2 can be recognized as the dependence
adjusted norm condition for non-stationary random variables.

2\. According to Shao [32] and Fourier inversion theorem(theorem 8.26 in
[38]), $\forall x=(x_{1},...,x_{n})^{T}\in\mathbf{R}^{n}$,

$\displaystyle\sum_{s=1}^{n}\sum_{j=1}^{n}x_{s}x_{j}K\left(\frac{s-j}{k_{n}}\right)=\int_{\mathbf{R}}\sum_{s=1}^{n}\sum_{j=1}^{n}x_{s}x_{j}\mathcal{F}K(z)\exp\left(2\pi\rm{i}z\frac{s-j}{k_{n}}\right)dz$
(8)
$\displaystyle=\int_{\mathbf{R}}FK(z)|\sum_{s=1}^{n}x_{s}\exp(\frac{2\pi\rm{i}zs}{k_{n}})|^{2}dz\geq
0$

Therefore, the matrix
$\left\\{K\left(\frac{s-j}{k_{n}}\right)\right\\}_{s,j=1,2,...,n}$ is
symmetric positive semi-definite. One example of $K$ satisfying assumption 7
is
$K(x)=\exp(-x^{2}/2)\Rightarrow\mathcal{F}K(x)=\sqrt{2\pi}\exp(-2\pi^{2}x^{2})$.

3\. Since

$\sum_{l=1}^{n}\left(\sum_{j=1}^{p}c_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\right)^{2}=\sum_{j=1}^{p}c^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}<\tau_{i}$
(9)

Assumption 6 requires no single element in the array
$\left\\{\sum_{j=1}^{p}c_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\right\\}_{l=1,...,n}$
dominates others. We add $1/n$ in $\tau_{i}$ to prevent the normalizing
parameters from being $0$.

4\. Assumption 2 implies a polynomial decay of $\epsilon$’s covariance, i.e.,
$\max_{1\leq i<j\leq
n}(1+j-i)^{\alpha_{\epsilon}}\times|\mathbf{E}\epsilon_{i}\epsilon_{j}|=O(1)$.
We will prove this in appendix A. Panel data may require different types of
dependency, e.g., clustered dependency or user-defined spatial
dependency(Vogelsang [8]). Our work needs to be adjusted if those dependencies
show up.

##  3 Consistency and Gaussian approximation theorem

This paper applies the ridge regression estimator introduced in Zhang and
Politis [28]. Use notations in section 2. Suppose the classical ridge
regression estimator $\widetilde{\beta}^{\star}$ and the de-biased estimator
$\widetilde{\beta}$ as

$\displaystyle\widetilde{\beta}^{\star}=(X^{T}X+\rho_{n}I_{p})^{-1}X^{T}y$
(10)
$\displaystyle\widetilde{\beta}=\widetilde{\beta}^{\star}+\rho_{n}Q(\Lambda^{2}+\rho_{n}I_{p})^{-1}Q^{T}\widetilde{\beta}^{\star}$

Define
$\widehat{\mathcal{N}}_{b_{n}}=\left\\{i=1,2,...,p:|\widetilde{\beta}_{i}|>b_{n}\right\\}$,
here $b_{n}>0$ is a given number(related to the sample size $n$). Define
$\widehat{\beta}=(\widehat{\beta}_{1},...,\widehat{\beta}_{p})^{T}$ such that
$\widehat{\beta}_{i}=\widetilde{\beta}_{i}\times\mathbf{1}_{i\in\widehat{\mathcal{N}}_{b_{n}}},i=1,2,...,p$;
and
$\widehat{\gamma}=(\widehat{\gamma}_{1},...,\widehat{\gamma}_{p_{1}})^{T}=M\widehat{\beta}$.
$\mathbf{1}_{i\in\widehat{\mathcal{N}}_{b_{n}}}=1$ if
$i\in\widehat{\mathcal{N}}_{b_{n}}$ and $0$ otherwise. This paper will use
$\widehat{\gamma}$ to estimate
$\gamma=(\gamma_{1},...\gamma_{p_{1}})^{T}=M\beta$.

###### Remark 2

Remark 2 in Zhang and Politis [28] explained why $\widetilde{\beta}$ helps
decrease the bias as well as the estimation error. Generally speaking, the
bias of a classical ridge regression estimator $\widetilde{\beta}^{\star}$ is
$-\rho_{n}Q(\Lambda^{2}+\rho_{n}I_{p})^{-1}Q^{T}\beta$. We can estimate the
bias by
$-\rho_{n}Q(\Lambda^{2}+\rho_{n}I_{p})^{-1}Q^{T}\widetilde{\beta}^{\star}$,
then subtract the estimated bias from $\widetilde{\beta}^{\star}$, yielding
the new estimator $\widetilde{\beta}$. Compared to
$\widetilde{\beta}^{\star}$, $\widetilde{\beta}$ increases the variance but
decreases the bias. So with a suitable choice of $\rho_{n}$, the estimation
error will decrease.

The consistency of $\widetilde{\beta}$ is not a free lunch. We need
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$
not to be very large(which is achieved if $\beta$ is sparse) and
$\rho_{n}/\lambda_{p}^{2}\to 0$ to maintain consistency.

The consistency of $\widetilde{\beta}$ consists of two aspects: the model
selection consistency, i.e.,
$\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$(defined in section 2) with
probability tending to $1$. And the consistency of $\widehat{\gamma}$, i.e.,
$\max_{i=1,...,p_{1}}|\widehat{\gamma}_{i}-\gamma_{i}|=o_{p}(1)$. Theorem 12
will prove both consistencies. Another important result in theorem 12 is that
the estimated errors $\widehat{\epsilon}=y-X\widehat{\beta}$ are close to the
underlying errors $\epsilon$. This result is the theoretical foundation for
bootstrap algorithm 19.

###### Theorem 1

Suppose assumption 1 to 5. Then

1.

$Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)=o(1)\
\text{and
}\max_{i=1,2,...,p_{1}}|\widehat{\gamma}_{i}-\gamma_{i}|=O_{p}(n^{-\eta})$
(11)

2\. Define
$\widehat{\epsilon}=(\widehat{\epsilon}_{1},...,\widehat{\epsilon}_{n})^{T}$
such that
$\widehat{\epsilon}_{i}=y_{i}-\sum_{j=1}^{p}x_{ij}\widehat{\beta}_{j}$, then

$\max_{i=1,2,...,n}|\widehat{\epsilon}_{i}-\epsilon_{i}|=o_{p}(n^{1/(2m)-\eta/2})$
(12)

(12) is not self-evident for a high dimensional linear model. We refer Mammen
[39] for a detailed explanation.

If the dimension of a parameter vector changes according to the sample size,
then the estimator $\widehat{\beta}$ does not have an asymptotic distribution,
and the classical central limit theorem fails. Nevertheless, statisticians
still can use normal random variables to approximate $\widehat{\beta}$’s
distribution. In reality, joint normal random variables can be generated
easily by a computer. So statisticians can perform Monte-Carlo simulations to
derive the asymptotic probability of an event or make a consistent confidence
interval. This idea is common for high dimensional statistics, see e.g.,
Chernozhukov et al. [40] and Zhang and Wu [41].

Following this idea, theorem 2 presents a Gaussian approximation theorem of
$\widehat{\gamma}$. Define $\xi=(\xi_{1},...,\xi_{n})^{T}$ as joint normal
random variables with $\mathbf{E}\xi=0,\mathbf{E}\xi\xi^{T}=\Sigma$(the
covariance matrix of errors $\epsilon$), and

$H(x)=Prob\left(\max_{i\in\mathcal{M}}\frac{1}{\tau_{i}}|\sum_{j=1}^{p}\sum_{l=1}^{n}c_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)p_{lj}\xi_{l}|\leq
x\right),\ x\in\mathbf{R}$ (13)

See section 2 for the meaning of notations. $\tau_{i}$ and $\mathcal{M}$ are
defined by choosing $b=b_{n}$. According to assumption 5, $\mathcal{M}$ is not
empty. So $H$ is well-defined.

###### Theorem 2

Suppose assumption 1 to 6. Then

$\sup_{x\in\mathbf{R}}|Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\sum_{j=1}^{p}m_{ij}\beta_{j}|}{\widehat{\tau}_{i}}\leq
x\right)-H(x)|=o(1)$ (14)

Here
$\widehat{\tau}_{i}=\sqrt{\frac{1}{n}+\sum_{j=1}^{p}\widehat{c}^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}}$
and $\widehat{c}_{ij}=\sum_{k\in\widehat{\mathcal{N}}_{b_{n}}}m_{ik}q_{kj}$;
$i=1,2,...,p_{1},j=1,2,...,p$.

Define $c_{1-\alpha}$ as the $1-\alpha$ quantile of $H(x)$, i.e.,
$c_{1-\alpha}=\inf\left\\{x\in\mathbf{R}:H(x)\geq 1-\alpha\right\\}$. Theorem
2 implies that the set

$\left\\{\gamma\in\mathbf{R}^{p_{1}}:\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\gamma_{i}|}{\widehat{\tau}_{i}}\leq
c_{1-\alpha}\right\\}$ (15)

is an asymptotic consistent confidence region for $\gamma$.

##  4 Dependent wild bootstrap

In order to calculate (15), we need a method that derives(or approximates)
$c_{1-\alpha}$. Since $\mathcal{M},\mathcal{N}_{b_{n}}$, and $\Sigma$ are
unknown a priori and $H$ does not have a closed-form formula, we develop a
Monte-Carlo algorithm(algorithm 19) that finds $c_{1-\alpha}$. The advantage
of algorithm 19 is that it avoids explicitly estimating $\Sigma$, which is
hard for non-stationary random variables $\epsilon$. Our method applies the
idea of the dependent wild bootstrap. It was first introduced by Shao [32],
and has been used in linear regression (Conley et al. [42]), time series
analysis (Fragkeskou and Paparoditis [43], Shao [32] and Doukhan et al. [44]),
high dimensional data analysis (Kurisu et al. [45]), etc.

Our target is to construct the confidence interval for $\gamma=M\beta$. And
Perform the statistical hypothesis testing for the null hypothesis $M\beta=z$
with $z=(z_{1},...,z_{p_{1}})^{T}$ a known vector versus the alternative
hypothesis $M\beta\neq z$. Notably, this test problem receives lots of
attentions in econometrics, e.g., Gonçalves and Vogelsang [46], Sun and Kim
[47], and Conley et al. [14].

######  Algorithm 1 (Bootstrap inference and hypothesis testing)

Input: Design matrix $X$ and dependent variables $y$ satisfying
$y=X\beta+\epsilon$; the linear combination matrix $M$. Ridge parameter
$\rho_{n}$, threshold $b_{n}$, bandwidth $k_{n}$, a kernel function
$K:\mathbf{R}\to[0,\infty)$. Nominal coverage probability $1-\alpha$, number
of bootstrap replicates $B$.

Additional input for testing: $z=(z_{1},...,z_{p_{1}})^{T}$.

1\. Derive the estimator $\widehat{\beta}$ defined in section 3 and
$\widehat{\tau}_{i},i=1,...,n$ defined in theorem 2. Then calculate
$\widehat{\gamma}=(\widehat{\gamma}_{1},...,\widehat{\gamma}_{p_{1}})^{T}=M\widehat{\beta},\
\widehat{\epsilon}=(\widehat{\epsilon}_{1},...,\widehat{\epsilon}_{n})^{T}=y-X\widehat{\beta}$.

2\. Generate joint normal random variables
$\varepsilon^{*}_{1},...,\varepsilon^{*}_{n}$ with mean $0$ and covariance
matrix $\left\\{K\left(\frac{i-j}{k_{n}}\right)\right\\}_{i,j=1,...,n}$. Then
define $\epsilon^{*}=(\epsilon^{*}_{1},...,\epsilon^{*}_{n})^{T}$ such that
$\epsilon^{*}_{i}=\widehat{\epsilon}_{i}\times\varepsilon^{*}_{i},i=1,2,...,n$.

3\. Calculate

$\displaystyle y^{*}=X\widehat{\beta}+\epsilon^{*}$ (16)
$\displaystyle\widetilde{\beta}^{\star*}=(X^{T}X+\rho_{n}I_{p})^{-1}X^{T}y^{*}$
$\displaystyle\widetilde{\beta}^{*}=(\widetilde{\beta}^{*}_{1},...,\widetilde{\beta}^{*}_{p})^{T}=\widetilde{\beta}^{\star*}+\rho_{n}Q(\Lambda^{2}+\rho_{n}I_{p})^{-1}Q^{T}\widetilde{\beta}^{\star*}$
$\displaystyle\widehat{\mathcal{N}}_{b_{n}}^{*}=\left\\{i=1,2,...,p:|\widetilde{\beta}_{i}^{*}|>b_{n}\right\\}$
$\displaystyle\text{and
}\widehat{\beta}^{*}=(\widehat{\beta}_{1}^{*},...,\widehat{\beta}^{*}_{p})\
\text{such that
}\widehat{\beta}^{*}_{i}=\widetilde{\beta}^{*}_{i}\times\mathbf{1}_{i\in\widehat{\mathcal{N}}_{b_{n}}^{*}}$

4\. Derive
$\widehat{\gamma}^{*}=(\widehat{\gamma}^{*}_{1},...,\widehat{\gamma}^{*}_{p_{1}})^{T}=M\widehat{\beta}^{*}$,
$\widehat{c}_{ij}^{*}=\sum_{k\in\widehat{\mathcal{N}}_{b_{n}}^{*}}m_{ik}q_{kj}$
and

$\widehat{\tau}_{i}^{*}=\sqrt{\frac{1}{n}+\sum_{j=1}^{p}\widehat{c}^{*2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}}$.
Then calculate

$\delta^{*}_{b}=\max_{i=1,2,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}$
(17)

5\. Repeat step 2 to 4 for $b=1,...,B$, and calculate the $1-\alpha$ sample
quantile $\widehat{c}^{*}_{1-\alpha}$ of $\delta^{*}_{b},b=1,2,...,B$.

6.a(For constructing the confidence region) The $1-\alpha$ confidence region
for the parameter of interest $\gamma=M\beta$ is given by the set

$\left\\{\gamma=(\gamma_{1},...,\gamma_{p_{1}})^{T}:\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\gamma_{i}|}{\widehat{\tau}_{i}}\leq\widehat{c}^{*}_{1-\alpha}\right\\}$
(18)

6.b(For hypothesis testing) Reject the null hypothesis if

$\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-z_{i}|}{\widehat{\tau}_{i}}>\widehat{c}^{*}_{1-\alpha}$
(19)

For a numerical sequence $\delta_{1}\leq\delta_{2}\leq...\leq\delta_{n}$, the
$1-\alpha$ sample quantile $C_{1-\alpha}$ is defined as

$C_{1-\alpha}=\delta_{i*}\text{ such that
}i*=\min\left\\{i=1,...,n:\frac{1}{n}\sum_{j=1}^{n}\mathbf{1}_{\delta_{j}\leq\delta_{i}}\geq
1-\alpha\right\\}$ (20)

Define the conditional quantile of $\delta^{*}_{b},b=1,...,B$ as

$c^{*}_{1-\alpha}=\inf\left\\{x\in\mathbf{R}:Prob^{*}\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}\leq
x\right)\geq 1-\alpha\right\\},\ 0<\alpha<1$ (21)

According to theorem 1.2.1 in Politis et al. [48], the sample quantile
$\widehat{c}^{*}_{1-\alpha}$ converges to $c^{*}_{1-\alpha}$ as $B\to\infty$.
So consistency of algorithm 19 is assured by showing

$Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\gamma_{i}|}{\widehat{\tau}_{i}}\leq
c^{*}_{1-\alpha}\right)\to 1-\alpha$. We derive this result in theorem 3.

###### Theorem 3

Suppose assumption 1 to 7. Then

$\displaystyle\sup_{x\in\mathbf{R}}|Prob^{*}\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}\leq
x\right)-H(x)|=o_{p}(1)$ (22) $\displaystyle\text{and
}Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\gamma_{i}|}{\widehat{\tau}_{i}}\leq
c^{*}_{1-\alpha}\right)\to 1-\alpha$

Here $H$ is defined in (13).

##  5 Numerical simulations

###  5.1 Selecting hyper-parameters

Statisticians need to fine-tune the ridge regression parameter $\rho_{n}$, the
threshold $b_{n}$, and the bandwidth $k_{n}$. $\rho_{n}$ and $b_{n}$ can be
chosen by (ten-fold) cross-validation, i.e., separate the design matrix $X$
and the dependent variables $y$ into disjoint training set
$(X_{train},y_{train})$ and validation set $(X_{valid},y_{valid})$; for each
choice of $\rho_{n}$ and $b_{n}$, use $(X_{train},y_{train})$ to fit
$\widehat{\beta}$; and calculate
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|y_{valid}-X_{valid}\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$.
The optimal $\rho_{n}$ and $b_{n}$ should minimize
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|y_{valid}-X_{valid}\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$.
This paper implements the ten-fold cross validation. See Arlot and Celisse
[49] for a further introduction on the cross validation methods.

Fine-tuning $k_{n}$ is more challenging. Politis and White [50] introduced an
automatic bandwidth selection algorithm; Shao [32] applied this algorithm for
the dependent wild bootstrap. Andrews [10] and Kim and Sun [11] considered
selecting bandwidth in the HAC estimation setting. Following Shao, this paper
applies Politis and White’s algorithm [50] on fitted residuals
$\widehat{\epsilon}=(\widehat{\epsilon}_{1},...,\widehat{\epsilon}_{n})^{T}=y-X\widehat{\beta}$
to select $k_{n}$. However, $\epsilon$ is not assumed to be stationary. So
this algorithm may result in a suboptimal bandwidth.

###  5.2 Generating data

The numerical experiment applies the linear model $y=X\beta+\epsilon$. $X$ is
generated by i.i.d. normal random variables with mean $0$ and variance $1$,
and is fixed in each experiment. $\beta=(\beta_{1},...,\beta_{p})^{T}$ is
generated by the following strategy

$\beta_{i}=0.1\times(i+6)\ \text{for }i=1,2,...,10\ \text{and }0\
\text{otherwise}$ (23)

Define $e_{i},i\in\mathbf{Z}$ as i.i.d. normal random variables with mean $0$
and variance $1$. Choose $a=(a_{1},a_{2},...,a_{n})^{T}$ such that
$a_{i}=e^{2}_{i}e^{2}_{i-1}-1$. Then define
$\epsilon=(\epsilon_{1},...,\epsilon_{n})^{T}=(1/4)\times Ha$(the factor $1/4$
avoids $\epsilon$’s variances being too large). Here
$H=(h_{ij})_{i,j=1,...,n}$, $h_{ij}=0$ for $j>i$, $h_{ii}=1$, $h_{ij}$ is
generated by uniform distribution in $[0.6,0.9]$ for $i-10\leq j<i$ and
$h_{ij}=s_{ij}/(i-j)^{z}$ for $j<i-10$. $s_{ij}$ is generated by uniform
distribution on $[-1,1]$. $H$ is fixed in each experiment. For
$\mathbf{E}a_{i}=\mathbf{E}e^{2}_{i}\times\mathbf{E}e^{2}_{i-1}-1=0$, we have
$\mathbf{E}\epsilon=0$. Moreover, $a$ is not white noise since
$\mathbf{E}a_{i}a_{i-1}=2$. $\epsilon$ is not linear because of $a$, and is
not stationary because of $H$. Figure 1 plots an observation of the errors
$\epsilon$, and the first $10\times 10$ elements of $\epsilon$’s covariance
matrix. In figure 1, the errors demonstrate strong dependency, i.e.,
$\epsilon_{i+1}$ is likely to be large if $\epsilon_{i}$ is large. In figure
1, $\epsilon$’s covariance matrix is not a Toeplitz matrix, so the
distribution of $\epsilon$ is not stationary.

Figure 1: Figure 1 plots an observation of the errors $\epsilon$, and figure 1
plots the heatmap for the first $10\times 10$ elements of $\epsilon$’s
covariance matrix. Values in each grid represent the corresponding covariance.
The covariance matrix is calculated by simulating 40000 samples of the random
vector $\epsilon=(\epsilon_{1},...,\epsilon_{n})^{T}$.

The linear combination matrix $M$ is generated by i.i.d. normal random
variables with mean 0.5 and variance $0.25$, and is fixed in each experiment.
The hyper-parameters $\rho_{n},b_{n},k_{n}$ are tuned by methods described in
section 5.1. The sample size $n$, dimension $p$ and $z$ vary in each
experiment. We store the information about experiments in table 1.

Table 1: Experiment parameters. ‘No.’ abbreviates ‘the experiment number’, $n$ is the sample size. $p$ is the dimension of parameters. $p_{1}$ is the number of linear combinations. $z$ is defined in section 5.2. $\rho_{n},b_{n},k_{n}$ respectively is the selected ridge parameter, the selected threshold, and the selected bandwidth defined in section 2. $\lambda_{p}$ is the smallest singular value of the design matrix $X$. We use R-Package ‘np’[51] to facilitate choosing $k_{n}$. No. | $n$ | $p$ | $p_{1}$ | $z$ | $\rho_{n}$ | $b_{n}$ | $k_{n}$ | $\lambda_{p}$
---|---|---|---|---|---|---|---|---
1 | 500 | 250 | 20 | 2.5 | 23.87 | 0.56 | 17.67 | 6.95
2 | 500 | 400 | 20 | 2.5 | 34.97 | 0.66 | 23.42 | 2.47
3 | 500 | 250 | 60 | 2.5 | 69.22 | 0.47 | 16.75 | 6.66
4 | 800 | 640 | 60 | 2.5 | 31.30 | 0.61 | 17.99 | 3.32
5 | 3000 | 1500 | 20 | 2.5 | 117.10 | 0.31 | 46.80 | 13.46
6 | 500 | 250 | 20 | 2.2 | 79.07 | 0.48 | 15.76 | 6.53

###  5.3 Performance of linear regression algorithms

This section compares the performance of the debiased and threshold ridge
regression(thsDeb) to the classical linear regression algorithms including
Lasso, the Ridge regression(Ridge), the threshold Lasso(thsLas), the threshold
ridge regression(thsRid) and the elastic net. The evaluation indices consist
of the estimation
loss${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M\widehat{\beta}-M\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$
and the prediction loss
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|X\widehat{\beta}-X\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$.
The latter one is frequently used in evaluating linear regression algorithms,
see e.g. Dalalyan et al. [52].

We are also interested in the model selection performance(i.e., whether an
algorithm can recover the sparsity of the underlying linear model or not) of
the linear regression algorithms. Following the idea of Fithian et al. [53],
we evaluate this through the frequency of model misspecification(i.e.,
$\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}$),
$P(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}})$; the average size of
model misspecification
$|\widehat{\mathcal{N}}_{b_{n}}\Delta\mathcal{N}_{b_{n}}|$ ($\Delta$ denotes
the symmetric difference, i.e. $A\Delta B=(A-B)\cup(B-A)$); and the average
false discovery rate
$|\widehat{\mathcal{N}}_{b_{n}}-\mathcal{N}_{b_{n}}|/\max(1,|\widehat{\mathcal{N}}_{b_{n}}|)$.
Since Lasso, ridge regression and elastic net do not have thresholds, for
these algorithms we consider $i\in\widehat{\mathcal{N}}_{b_{n}}$ if
$|\widehat{\beta}_{i}|>0.01$.

According to figure 2, figure 3 and table 2, the debiased and threshold ridge
regression, the threshold Lasso and the threshold ridge regression have small
estimation loss and prediction loss compared to the other methods. Besides,
these three methods are able to recover the underlying sparsity of the linear
model correctly. Notably, even if Lasso is famous for generating sparse linear
regression estimators, further thresholding still significantly improves its
performance.

The necessity of debiasing can be illustrated in figure 2 and figure 3.
Although the threshold ridge regression has a small estimation loss when
$\rho_{n}$ is close to its optimal value, the estimation loss surges when
$\rho_{n}$ deviates from its optimal value. However, after debiasing, even if
a cross validation selects a sub-optimal $\rho_{n}$, the estimator’s
performance does not notably deteriorate. In reality, cross validation
algorithms cannot assure selecting the optimal ridge parameter $\rho_{n}$. So
an algorithm that is not sensitive to the slight fluctuations in hyper-
parameters is preferable.

Figure 2 also plots the relation between thresholds and the estimation losses.
In the experiments, the cross validation always selects sub-optimal
thresholds. However, the sub-optimal thresholds form wide-intervals.
Therefore, slight fluctuations in the threshold $b_{n}$ should not deteriorate
the performance of the regression algorithm as well.

(a) Experiment 1

(b) Experiment 2

(c) Experiment 3

(d) Experiment 4

Figure 2: Estimation losses of linear regression algorithms for case 1 to 4.
‘thsDeb’ represents the debiased and threshold ridge regression method;
‘thsLas’ represents the threshold Lasso method and ‘thsRid’ represents the
threshold ridge regression method. Dots reveal the optimal ridge parameters
for each method selected by 10-fold cross validation. y-axis records the error
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M\widehat{\beta}-M\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$,
here $\widehat{\beta}$ is the estimator of parameters $\beta$ by different
methods. We choose $l_{1}$ ratio $0.5$ in the ElasticNet. The small graphs
below each of the graphs plot the error of the debiased and threshold ridge
regression method with respect to different thresholds. $\rho_{n}$ and $b_{n}$
are chosen by 10-fold cross validation.

(a) Experiment 5

(b) Experiment 6

Figure 3: Estimation losses of linear regression algorithms for case 5 to 6. Meaning of notations coincide with figure 2. Table 2: Model selection performance of various linear regression estimators. Hyper-parameters are chosen by ten-fold cross validation. The overscore represents calculating the sample mean among 1000 simulations. ‘FDR’ abbreviates the false discovery rate. We omit experiment 3 here for it coincides with experiment 1. No. | Algorithm | $P\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)$ | $\overline{\widehat{\mathcal{N}}_{b_{n}}\Delta\mathcal{N}_{b_{n}}}$ | FDR | $\overline{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|X\widehat{\beta}-X\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}}$
---|---|---|---|---|---
1 | thsDeb | 0.211 | 0.334 | 0.010 | 11.335
| Lasso | 0.943 | 8.663 | 0.386 | 14.992
| thsLas | 0.190 | 0.261 | 0.005 | 10.212
| Ridge | 1.0 | 112.18 | 0.918 | 31.216
| thsRid | 0.323 | 0.411 | 0.006 | 12.369
| Elastic net | 1.0 | 27.71 | 0.706 | 18.992
2 | thsDeb | 0.790 | 1.233 | 0.014 | 21.575
| Lasso | 1.0 | 86.28 | 0.892 | 25.873
| thsLas | 0.142 | 0.323 | 0.015 | 10.797
| Ridge | 1.0 | 181.10 | 0.948 | 35.431
| thsRid | 0.801 | 1.214 | 0.019 | 27.727
| Elastic net | 1.0 | 122.05 | 0.923 | 30.638
4 | thsDeb | 0.377 | 0.568 | 0.009 | 17.365
| Lasso | 1.0 | 125.95 | 0.924 | 31.442
| thsLas | 0.079 | 0.259 | 0.011 | 10.667
| Ridge | 1.0 | 284.67 | 0.966 | 46.353
| thsRid | 0.409 | 0.585 | 0.006 | 18.101
| Elastic net | 1.0 | 184.87 | 0.948 | 38.673
5 | thsDeb | 0.031 | 0.045 | 0.004 | 9.781
| Lasso | 1.0 | 32.81 | 0.728 | 18.394
| thsLas | 0.093 | 0.248 | 0.017 | 9.703
| Ridge | 1.0 | 419.87 | 0.977 | 60.853
| thsRid | 0.035 | 0.055 | 0.004 | 9.933
| Elastic net | 1.0 | 79.88 | 0.879 | 27.446
6 | thsDeb | 0.138 | 0.236 | 0.012 | 10.856
| Lasso | 0.967 | 9.792 | 0.422 | 14.445
| thsLas | 0.079 | 0.157 | 0.010 | 9.229
| Ridge | 1.0 | 105.775 | 0.913 | 30.164
| thsRid | 0.173 | 0.261 | 0.010 | 10.874
| Elastic net | 1.0 | 28.313 | 0.711 | 18.617

###  5.4 Performance of bootstrap algorithm 19

This section evaluates the performance of algorithm 19 as well as the
bootstrap algorithms designed for independent errors, -including Efron’s
bootstrap(Efron [54] and Chatterjee and Lahiri [3]) and the wild
bootstrap(Mammen [55])-, for the linear model with dependent errors. The
results are demonstrated in table 3.

In table 3, the bootstrap algorithms have similar performance and are able to
capture the correct $1-\alpha$ confidence intervals. Meanwhile, the dependent
wild bootstrap tends to generate a wider confidence interval than Efron’s
bootstrap and the wild bootstrap. Notably, the number of simultaneous linear
combinations $p_{1}$ affects the performance of bootstrap algorithms. If
$p_{1}$ is too large, then the bootstrap algorithms may generate an
unnecessarily wide confidence interval. Numerical results show that the
Efron’s bootstrap and the wild bootstrap are robust to the dependent and
heteroskedastic errors. However, the consistency of those algorithms cannot be
assured.

Table 3: Performance of various bootstrap algorithms. The number of bootstrap replicates is 1000 and the nominal coverage probability is $90\%$. We use $P$ to abbreviate the real coverage probability and $C$ to abbreviate the average $90\%$ quantile of $\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}$. $P$ and $C$ are derived through 1000 simulations. ‘Efron’ means ‘Efron’s bootstrap’ and ‘Wild’ means the wild bootstrap.  No | Algorithm 19 $P$ | Algorithm 19 $C$ | Efron $P$ | Efron $C$ | Wild $P$ | Wild $C$
---|---|---|---|---|---|---
1 | $90.7\%$ | 9.32 | $88.0\%$ | 8.83 | $88.4\%$ | 8.86
2 | $89.1\%$ | 17.14 | $87.7\%$ | 16.71 | $87.2\%$ | 16.76
3 | $95.1\%$ | 11.77 | $94.4\%$ | 11.26 | $94.0\%$ | 11.31
4 | $89.8\%$ | 14.38 | $89.8\%$ | 13.89 | $88.8\%$ | 13.86
5 | $90.2\%$ | 4.72 | $89.1\%$ | 4.62 | $89.7\%$ | 4.64
6 | $92.8\%$ | 11.38 | $92.2\%$ | 10.70 | $92.1\%$ | 10.75

##  6 Conclusion

This paper applies the debiased and threshold ridge regression method to a
linear model $y=X\beta+\epsilon$. Then it derives a dependent wild bootstrap
algorithm that constructs confidence intervals for linear combinations of
parameters $\gamma=M\beta$; or tests the statistical hypothesis

$M\beta=z\textit{ versus }M\beta\neq z$ (24)

Here $M$ is a known matrix, and $z$ is the given expected value of $M\beta$.

Our work suits a high dimensional setting, i.e., the dimension of parameters
$\beta$ may grow as the sample size increases. Besides, the adopted regression
method and the bootstrap procedure are consistent for dependent, non-
stationary, heteroskedastic errors $\epsilon$. $\epsilon$ is not necessary to
be a linear process as well.

Numerical simulations indicate that the ridge regression estimator has a
comparable model selection and estimation performance to complex methods like
the threshold Lasso. Furthermore, it is robust to the non-optimal choice of
the ridge parameter $\rho_{n}$ as well as the threshold $b_{n}$. Therefore, it
can be considered as a practical method to handle real-life problems.

Our work generalizes results in Conley et al. [14] to a more realistic
setting. We have not yet seen many works on a high-dimensional linear model
with non-stationary errors. So this work should bring new insights to this
field.

## References

  * [1] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2):301–320, 2005.
  * [2] Hui Zou.  The adaptive lasso and its oracle properties.  Journal of the American Statistical Association, 101(476):1418–1429, 2006.
  * [3] A. Chatterjee and S. N. Lahiri.  Bootstrapping lasso estimators.  Journal of the American Statistical Association, 106(494):608–625, 2011.
  * [4] A. Chatterjee and S. N. Lahiri.  Asymptotic properties of the residual bootstrap for lasso estimators.  Proceedings of the American Mathematical Society, 138(12):4497–4509, 2010.
  * [5] Robert A. Stine.  Bootstrap prediction intervals for regression.  Journal of the American Statistical Association, 80(392):1026–1031, 1985.
  * [6] Yunyi Zhang and Dimitris N. Politis.  Bootstrap prediction intervals with asymptotic conditional validity and unconditional guarantees.  (2005.09145), 2021.
  * [7] George A. F. Seber and Alan J. Lee.  Linear Regression Analysis.  John Wiley $\&$ Sons, 2 edition, 2003.
  * [8] Timothy J. Vogelsang.  Heteroskedasticity, autocorrelation, and spatial correlation robust inference in linear panel models with fixed-effects.  Journal of Econometrics, 166(2):303–319, 2012.
  * [9] Mitchell A. Petersen.  Estimating Standard Errors in Finance Panel Data Sets: Comparing Approaches.  The Review of Financial Studies, 22(1):435–480, 06 2008.
  * [10] Donald W. K. Andrews.  Heteroskedasticity and autocorrelation consistent covariance matrix estimation.  Econometrica, 59(3):817–858, 1991.
  * [11] Min Seong Kim and Yixiao Sun.  Spatial heteroskedasticity and autocorrelation consistent estimation of covariance matrix.  Journal of Econometrics, 160(2):349–371, 2011.
  * [12] Harry H. Kelejian and Ingmar R. Prucha.  Hac estimation in a spatial framework.  Journal of Econometrics, 140(1):131–154, 2007.
  * [13] Yixiao Sun and Xuexin Wang.  An asymptotically f-distributed chow test in the presence of heteroscedasticity and autocorrelation.  Econometric Reviews, 0(0):1–35, 2021.
  * [14] Timothy Conley, Sílvia Gonçalves, Min Seong Kim, and Benoit Perron.  Bootstrap inference under cross sectional dependence.  unpublished, 2019.
  * [15] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288, 1996.
  * [16] Peng Zhao and Bin Yu.  On model selection consistency of lasso.  J. Mach. Learn. Res., 7:2541– 2563, 2006.
  * [17] Nicolai Meinshausen and Peter Bühlmann.  High-dimensional graphs and variable selection with the Lasso.  The Annals of Statistics, 34(3):1436 – 1462, 2006.
  * [18] Nicolai Meinshausen and Bin Yu.  Lasso-type recovery of sparse representations for high-dimensional data.  The Annals of Statistics, 37(1):246 – 270, 2009.
  * [19] Cun-Hui Zhang and Stephanie S. Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 76(1):217–242, 2014.
  * [20] Xianyang Zhang and Guang Cheng.  Simultaneous inference for high-dimensional linear models.  Journal of the American Statistical Association, 112(518):757–768, 2017.
  * [21] Eitan Greenshtein and Ya’Acov Ritov.  Persistence in high-dimensional linear predictor selection and the virtue of overparametrization.  Bernoulli, 10(6):971 – 988, 2004.
  * [22] Peter Bühlmann and Sara van de Geer.  Statistics for High-Dimensional Data.  Springer-Verlag Berlin Heidelberg, Berlin, Heidelberg, 1 edition, 2011\.
  * [23] Jianqing Fan and Runze Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360, 2001.
  * [24] Jason D. Lee, Dennis L. Sun, Yuekai Sun, and Jonathan E. Taylor.  Exact post-selection inference, with application to the lasso.  The Annals of Statistics, 44(3):907 – 927, 2016.
  * [25] Hanzhong Liu and Bin Yu.  Asymptotic properties of lasso+mls and lasso+ridge in sparse high-dimensional linear regression.  Electron. J. Statist., 7:3124–3169, 2013.
  * [26] Ryan J. Tibshirani, Alessandro Rinaldo, Rob Tibshirani, and Larry Wasserman.  Uniform asymptotic inference and the bootstrap after model selection.  The Annals of Statistics, 46(3):1255 – 1287, 2018.
  * [27] Jun Shao and Xinwei Deng.  Estimation in high-dimensional linear models with deterministic design matrices.  Ann. Statist., 40(2):812–831, 04 2012.
  * [28] Yunyi Zhang and D. N. Politis.  Ridge regression revisited: Debiasing, thresholding and bootstrap.  (arXiv:2009.08071), 2020.
  * [29] Wei-Biao Wu and Ying Nian Wu.  Performance bounds for parameter estimates of high-dimensional linear models with correlated errors.  Electronic Journal of Statistics, 10(1):352 – 379, 2016.
  * [30] Yuefeng Han and Ruey S. Tsay.  High-dimensional linear regression for dependent data with application to nowcasting.  Statistica Sinica, 30:1797–1827, 2020.
  * [31] Peter J Brockwell and Richard A Davis.  Time Series: Theory and Methods.  Springer-Verlag, Berlin, Heidelberg, 1986.
  * [32] Xiaofeng Shao.  The dependent wild bootstrap.  Journal of the American Statistical Association, 105(489):218–235, 2010.
  * [33] Roger A. Horn and Charles R. Johnson.  Matrix Analysis.  Cambridge University Press, 2012.
  * [34] Wei Biao Wu.  Nonlinear system theory: Another look at dependence.  Proceedings of the National Academy of Sciences, 102(40):14150–14154, 2005.
  * [35] Wei-Biao Wu and Ying Nian Wu.  Performance bounds for parameter estimates of high-dimensional linear models with correlated errors.  Electron. J. Statist., 10(1):352–379, 2016.
  * [36] Jun Shao.  Mathematical Statistics.  Springer-Verlag New York, 2003.
  * [37] Z. D. Bai and Y. Q. Yin.  Limit of the Smallest Eigenvalue of a Large Dimensional Sample Covariance Matrix.  The Annals of Probability, 21(3):1275 – 1294, 1993.
  * [38] Gerald B. Folland.  Real Analysis.  John Wiley $\&$ Sons, 2 edition, 2007.
  * [39] Enno Mammen.  Empirical process of residuals for high-dimensional linear models.  The Annals of Statistics, 24(1):307 – 335, 1996.
  * [40] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  Ann. Statist., 41(6):2786–2819, 12 2013.
  * [41] Danna Zhang and Wei Biao Wu.  Gaussian approximation for high dimensional time series.  The Annals of Statistics, 45(5):1895 – 1919, 2017.
  * [42] Timothy Conley, Sílvia Gonçalves, Min Seong Kim, and Benoit Perron.  Bootstrap inference under cross sectional dependence.  unpublished, 2019.
  * [43] Maria Fragkeskou and Efstathios Paparoditis.  Extending the range of validity of the autoregressive (sieve) bootstrap.  Journal of Time Series Analysis, 39(3):356–379, 2018.
  * [44] Paul Doukhan, Gabriel Lang, Anne Leucht, and Michael H. Neumann.  Dependent wild bootstrap for the empirical process.  Journal of Time Series Analysis, 36(3):290–314, 2015.
  * [45] Daisuke Kurisu, Kengo Kato, and Xiaofeng Shao.  Gaussian approximation and spatially dependent wild bootstrap for high-dimensional spatial data.  (arXiv:2103.10720), 2021.
  * [46] Sílvia Gonçalves and Timothy J. Vogelsang.  Block bootstrap hac robust tests: the sophistication of the naive bootstrap.  Econometric Theory, 27(4):745 – 791, 2011.
  * [47] Yixiao Sun and Min Seong Kim.  Asymptotic f-test in a gmm framework with cross-sectional dependence.  The Review of Economics and Statistics, 97(1):210–223, 2015.
  * [48] Dimitris N. Politis, Joseph P. Romano, and Michael Wolf.  Subsampling.  Springer-Verlag New York, 1999.
  * [49] Sylvain Arlot and Alain Celisse.  A survey of cross-validation procedures for model selection.  Statistics Surveys, 4:40 – 79, 2010.
  * [50] Dimitris N. Politis and Halbert White.  Automatic block-length selection for the dependent bootstrap.  Econometric Reviews, 23(1):53–70, 2004.
  * [51] Tristen Hayfield and Jeffrey S. Racine.  Nonparametric econometrics: The np package.  Journal of Statistical Software, 27(5), 2008.
  * [52] Arnak S. Dalalyan, Mohamed Hebiri, and Johannes Lederer.  On the prediction performance of the Lasso.  Bernoulli, 23(1):552 – 581, 2017.
  * [53] William Fithian, Dennis Sun, and Jonathan Taylor.  Optimal inference after model selection.  (arXiv:1410.2597), 2017.
  * [54] B. Efron.  Bootstrap Methods: Another Look at the Jackknife.  The Annals of Statistics, 7(1):1 – 26, 1979.
  * [55] Enno Mammen.  Bootstrap and Wild Bootstrap for High Dimensional Linear Models.  The Annals of Statistics, 21(1):255 – 285, 1993.
  * [56] Bernt øksendal.  Stochastic Differential Equations.  Springer-Verlag Berlin Heidelberg, 2003.
  * [57] D. L. Burkholder, B. J. Davis, and R. F. Gundy.  Integral inequalities for convex functions of operators on martingales.  In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory, pages 223–240, Berkeley, Calif., 1972. University of California Press.
  * [58] P. Whittle.  Bounds for the moments of linear and quadratic forms in independent variables.  Theory of Probability & Its Applications, 5(3):302–305, 1960.
  * [59] Mengyu Xu, Danna Zhang, and Wei Biao Wu.  Pearson’s chi-squared statistics: approximation theory and beyond.  Biometrika, 106(3):716–723, 04 2019.
  * [60] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Comparison and anti-concentration bounds for maxima of gaussian random vectors.  Probability Theory and Related Fields, 162:47–70, 2015.

##  Appendix A Some useful results

Notations in the appendix should coincide with section 2 unless we specify
their meanings. This section starts with some useful corollaries of assumption
2. From corollary C.9 in [56],
$\epsilon_{i}=\lim_{j\to\infty}\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j}$
almost surely; and
$\lim_{j\to\infty}\|\epsilon_{i}-\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j}\|_{m}=0,\forall
i=1,2,...,n$. For

$\|\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j}-\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j-1}\|_{m}=\|\mathbf{E}(\epsilon_{i}-\epsilon_{i,j})|\mathcal{F}_{i,j}\|_{m}\leq\|\epsilon_{i}-\epsilon_{i,j}\|_{m}\leq\max_{i=1,2,...,n}\delta_{i,j,m}$
(25)

Therefore, assumption 2 implies $\sup_{i=1,2,...,n;k\geq
0}(k+1)^{\alpha_{\epsilon}}\sum_{j=k}^{\infty}\|\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j}-\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,j-1}\|_{m}=O(1)$.
For any $a_{i}\in\mathbf{R},i=1,2,...,n$ and $s\in\mathbf{Z}^{+}$, define
$M_{i,s}=\sum_{j=n+1-i}^{n}a_{j}(\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s-1})$,
then $M_{i,s}$ is measurable in the filter $\mathcal{F}_{n,s+i-1}$. Besides,
$M_{i+1,s}-M_{i,s}=a_{n-i}\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s}-a_{n-i}\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s-1}$.
Apply $\pi-\lambda$ theorem to the

$\lambda-$system
$\left\\{A\in\mathcal{F}_{n,s+i-1}\Big{|}\mathbf{E}(\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s})\mathbf{1}_{A}=\mathbf{E}(\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s-1})\mathbf{1}_{A}\right\\}$
and the $\pi-$system $\\{A_{n}\times A_{n-1}\times...\times A_{n-s-i+1}\\}$.
Here $A_{i}$ is generated by $e_{i}$. Then
$\mathbf{E}(\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s}-\mathbf{E}\epsilon_{n-i}|\mathcal{F}_{n-i,s-1})\Big{|}\mathcal{F}_{n,s+i-1}=0$
almost surely. In other words, $\\{M_{i,s}\\}_{i=1,2,...,n}$ is a martingale.
From Burkholder’s inequality(theorem 1.1 in [57]),

$\|M_{n,s}\|_{m}\leq
C\sqrt{\|\sum_{j=1}^{n}a^{2}_{j}(\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s-1})^{2}\|_{m/2}}\leq
C\sqrt{\sum_{i=1}^{n}a^{2}_{i}}\times\max_{i=1,2,...,n}\delta_{i,s,m}$ (26)

Here $C$ is independent of $s$. In particular, from theorem 2 in [58] and
assumption 2,

$\displaystyle\|\sum_{i=1}^{n}a_{i}\epsilon_{i}\|_{m}\leq\|\sum_{i=1}^{n}a_{i}\mathbf{E}\epsilon_{i}|\mathcal{F}_{i,0}\|_{m}+\sum_{s=1}^{\infty}\|M_{n,s}\|_{m}$
(27) $\displaystyle\leq
C\sqrt{\sum_{i=1}^{n}a^{2}_{i}\|\epsilon_{i}\|^{2}_{m}}+C\sqrt{\sum_{i=1}^{n}a^{2}_{i}}\times\sum_{s=1}^{\infty}\max_{i=1,2,...,n}\delta_{i,s,m}=O\left(\sqrt{\sum_{i=1}^{n}a^{2}_{i}}\right)$

Here $C$ is a constant. Since $m>2$, (27) implies $\Sigma$’s largest
eigenvalue has order $O(1)$. For $\forall 1\leq i<j\leq n$, assumption 2 and
(25) imply

$\displaystyle|\mathbf{E}\epsilon_{i}\epsilon_{j}|=|\mathbf{E}\epsilon_{i}(\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,j-i-1})|\leq\|\epsilon_{i}\|_{2}\times\|\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,j-i-1}\|_{2}$
(28)
$\displaystyle\leq\|\epsilon_{i}\|_{m}\times\sum_{s=j-i-1}^{\infty}\|\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s+1}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}\|_{m}\Rightarrow\max_{1\leq
i<j\leq
n}(1+j-i)^{\alpha_{\epsilon}}\times|\mathbf{E}\epsilon_{i}\epsilon_{j}|=O(1)$

In other words, assumption 2 implies a polynomial decay on the error’s
covariance.

For $\forall\tau,\psi>0,z\in\mathbf{R}$, define
$F_{\tau}(x_{1},...,x_{n})=\frac{1}{\tau}\log(\sum_{i=1}^{n}\exp(\tau
x_{i}))$;
$G_{\tau}(x_{1},...,x_{n})=\frac{1}{\tau}\log(\sum_{i=1}^{n}\exp(\tau
x_{i})+\sum_{i=1}^{n}\exp(-\tau
x_{i}))=F_{\tau}(x_{1},...,x_{n},-x_{1},...,-x_{n})$;

$g_{0}(x)=(1-\min(1,\max(x,0))^{4})^{4}$; $g_{\psi,z}(x)=g_{0}(\psi(x-z))$;
and
$h_{\tau,\psi,z}(x_{1},...,x_{n})=g_{\psi,z}\left(G_{\tau}(x_{1},...,x_{n})\right)$.
From lemma A.2 and (8) in [40] and (S1) to (S5) in [59],
$g_{*}=\sup_{x\in\mathbf{R}}(|g^{{}^{\prime}}_{0}(x)|+|g^{{}^{\prime\prime}}_{0}(x)|+|g^{{}^{\prime\prime\prime}}_{0}(x)|)<\infty$;
$\mathbf{1}_{x\leq z}\leq g_{\psi,z}(x)\leq\mathbf{1}_{x\leq z+1/\psi}$; and
$\sup_{x,z\in\mathbf{R}}|g^{{}^{\prime}}_{\psi,z}(x)|\leq g_{*}\psi$,
$\sup_{x,z\in\mathbf{R}}|g^{{}^{\prime\prime}}_{\psi,z}(x)|\leq
g_{*}\psi^{2}$,
$\sup_{x,z\in\mathbf{R}}|g^{{}^{\prime\prime\prime}}_{\psi,z}(x)|\leq
g_{*}\psi^{3}$. Besides, $\frac{\partial F_{\tau}}{\partial x_{i}}\geq 0$;
$\sum_{i=1}^{n}\frac{\partial F_{\tau}}{\partial x_{i}}=1$;
$\sum_{i=1}^{n}\sum_{j=1}^{n}|\frac{\partial^{2}F_{\tau}}{\partial
x_{i}\partial x_{j}}|\leq 2\tau$;
$\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}|\frac{\partial^{3}F_{\tau}}{\partial
x_{i}\partial x_{j}\partial x_{k}}|\leq 6\tau^{2}$; and
$F_{\tau}(x_{1},...,x_{n})-\frac{\log(n)}{\tau}\leq\max_{i=1,2,...,n}x_{i}\leq
F_{\tau}(x_{1},...,x_{n})$. Therefore,

$G_{\tau}(x_{1},...,x_{n})-\frac{\log(2n)}{\tau}\leq\max_{i=1,2,...,n}|x_{i}|\leq
G_{\tau}(x_{1},...,x_{n})$ (29)

Since $\frac{\partial G_{\tau}}{\partial x_{i}}=\frac{\partial
F_{\tau}}{\partial x_{i}}-\frac{\partial F_{\tau}}{\partial x_{i+n}}$, we get
$\sum_{i=1}^{n}|\frac{\partial G_{\tau}}{\partial x_{i}}|\leq 1$. For
$\frac{\partial^{2}G_{\tau}}{\partial x_{i}\partial
x_{j}}=\frac{\partial^{2}F_{\tau}}{\partial x_{i}\partial
x_{j}}-\frac{\partial^{2}F_{\tau}}{\partial x_{i}\partial
x_{j+n}}-\frac{\partial^{2}F_{\tau}}{\partial x_{i+n}\partial
x_{j}}+\frac{\partial^{2}F_{\tau}}{\partial x_{i+n}\partial x_{j+n}}$, we have
$\sum_{i=1}^{n}\sum_{j=1}^{n}|\frac{\partial^{2}G_{\tau}}{\partial
x_{i}\partial x_{j}}|\leq 2\tau$. For

$\displaystyle\frac{\partial^{3}G_{\tau}}{\partial x_{i}\partial x_{j}\partial
x_{k}}=\left(\frac{\partial^{3}F_{\tau}}{\partial x_{i}\partial x_{j}\partial
x_{k}}+\frac{\partial^{3}F_{\tau}}{\partial x_{i}\partial x_{j+n}\partial
x_{k+n}}+\frac{\partial^{3}F_{\tau}}{\partial x_{i+n}\partial x_{j}\partial
x_{k+n}}+\frac{\partial^{3}F_{\tau}}{\partial x_{i+n}\partial x_{j+n}\partial
x_{k}}\right)$ (30) $\displaystyle-\left(\frac{\partial^{3}F_{\tau}}{\partial
x_{i}\partial x_{j}\partial x_{n+k}}+\frac{\partial^{3}F_{\tau}}{\partial
x_{i}\partial x_{j+n}\partial x_{k}}+\frac{\partial^{3}F_{\tau}}{\partial
x_{i+n}\partial x_{j}\partial x_{k}}+\frac{\partial^{3}F_{\tau}}{\partial
x_{i+n}\partial x_{j+n}\partial x_{k+n}}\right)$

we have
$\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}|\frac{\partial^{3}G_{\tau}}{\partial
x_{i}\partial x_{j}\partial x_{k}}|\leq 6\tau^{2}$. Since $\frac{\partial
h_{\tau,\psi,z}}{\partial
x_{i}}=g^{{}^{\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\frac{\partial
G_{\tau}}{\partial x_{i}}$, we get $\sum_{i=1}^{n}|\frac{\partial
h_{\tau,\psi,z}}{\partial x_{i}}|\leq g_{*}\psi$.

$\displaystyle\frac{\partial^{2}h_{\tau,\psi,z}}{\partial x_{i}\partial
x_{j}}=g^{{}^{\prime\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\frac{\partial
G_{\tau}}{\partial x_{i}}\frac{\partial G_{\tau}}{\partial
x_{j}}+g^{{}^{\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\frac{\partial^{2}G_{\tau}}{\partial
x_{i}\partial x_{j}}$ (31)
$\displaystyle\Rightarrow\sum_{i=1}^{n}\sum_{j=1}^{n}|\frac{\partial
h_{\tau,\psi,z}}{\partial x_{i}\partial x_{j}}|\leq
g_{*}\psi^{2}+2g_{*}\psi\tau$
$\displaystyle\frac{\partial^{3}h_{\tau,\psi,z}}{\partial x_{i}\partial
x_{j}\partial
x_{k}}=g^{{}^{\prime\prime\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\frac{\partial
G_{\tau}}{\partial x_{i}}\frac{\partial G_{\tau}}{\partial
x_{j}}\frac{\partial G_{\tau}}{\partial
x_{k}}+g^{{}^{\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\frac{\partial^{3}G_{\tau}}{\partial
x_{i}\partial x_{j}\partial x_{k}}$
$\displaystyle+g^{{}^{\prime\prime}}_{\psi,z}(G_{\tau}(x_{1},...,x_{n}))\left(\frac{\partial^{2}G_{\tau}}{\partial
x_{i}\partial x_{k}}\frac{\partial G_{\tau}}{\partial
x_{j}}+\frac{\partial^{2}G_{\tau}}{\partial x_{j}\partial x_{k}}\frac{\partial
G_{\tau}}{\partial x_{i}}+\frac{\partial^{2}G_{\tau}}{\partial x_{i}\partial
x_{j}}\frac{\partial G_{\tau}}{\partial x_{k}}\right)$
$\displaystyle\Rightarrow\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}|\frac{\partial^{3}h_{\tau,\psi,z}}{\partial
x_{i}\partial x_{j}\partial x_{k}}|\leq
g_{*}\psi^{3}+6g_{*}\psi\tau^{2}+6g_{*}\tau\psi^{2}$

##  Appendix B Proof of theorem 12 and theorem 2

We first prove necessary lemmas.

###### Lemma 1

Suppose random variables $\epsilon_{i},i=1,2,...,n$ satisfy assumption 2, and
$\\{\gamma_{ij}\\}_{i=1,2,...,k,j=1,2,...,n}\in\mathbf{R}^{k\times n}$ satisfy

$\max_{i=1,2,...,k}\sum_{j=1}^{n}\gamma^{2}_{ij}\leq D^{2}$. Then

$\max_{i=1,2,...,k}|\sum_{j=1}^{n}\gamma_{ij}\epsilon_{j}|=O_{p}(k^{1/m}D)$
(32)

###### Proof.

Form (27), $\forall\xi>0$,

$\displaystyle
Prob\left(\max_{i=1,2,...,k}|\sum_{j=1}^{n}\gamma_{ij}\epsilon_{j}|>\xi\right)\leq\frac{1}{\xi^{m}}\sum_{i=1}^{k}\|\sum_{j=1}^{n}\gamma_{ij}\epsilon_{j}\|_{m}^{m}=O\left(\frac{k}{\xi^{m}}\times\left(\max_{i=1,2,...,k}\sum_{j=1}^{n}\gamma^{2}_{ij}\right)^{m/2}\right)$
(33)

(32) is proved by choosing $\xi=Ck^{1/m}D$ with a sufficiently large constant
$C$. ∎

Lemma 32 lays the theoretical foundation for the threshold ridge regression
estimator $\widehat{\beta}$’s (model selection) consistency. Lemma 35 is a
corollary of Chernozhukov et al. [60]. It introduces some properties of a
Gaussian random vector.

###### Lemma 2

1\. Suppose $\epsilon=(\epsilon_{1},...,\epsilon_{n})^{T}$ are joint normal
random variables, $\mathbf{E}\epsilon=0$,
$\mathbf{E}\epsilon\epsilon^{T}=\Sigma=(\sigma_{ij})_{i,j=1,...,n}$ is non-
singular; and $\exists 0<c_{0}\leq C_{0}<\infty$ such that
$c_{0}\leq\sigma_{ii}\leq C_{0},i=1,2,...,n$. Then

$\sup_{x\in\mathbf{R}}|Prob(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x+\delta)-Prob(\max_{i=1,2,...,n}|\epsilon_{i}|\leq x)|\leq
C\delta(1+\sqrt{\log(n)}+\sqrt{|\log(\delta)|})$ (34)

Here $C$ only depends on $c_{0}$ and $C_{0}$.

2\. In addition suppose $\xi=(\xi_{1},...,\xi_{n})^{T}$ are joint normal
random variables with $\mathbf{E}\xi=0$,
$\mathbf{E}\xi\xi^{T}=\Sigma^{*}=(\sigma^{*}_{ij})_{i,j=1,2,...,n}$, and
$\Delta=\max_{i,j=1,2,...,n}|\sigma_{ij}-\sigma^{*}_{ij}|<1$. Then

$\displaystyle\sup_{x\in\mathbf{R}}|Prob(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x)-Prob(\max_{i=1,2,...,n}|\xi_{i}|\leq x)|$ (35) $\displaystyle\leq
C^{*}\left(\Delta^{1/3}(1+\log^{3}(n))+\frac{\Delta^{1/6}}{1+\log^{1/4}(n)}\right)$

for $n=1,2,...$. Here $C^{*}$ only depends on $c_{0}$ and $C_{0}$.

We emphasize that $\Sigma^{*}$ can be singular. From theorem 4.1.5 in [33]
$\Sigma^{*}=Q\Lambda
Q^{T},\Lambda=diag(\lambda_{1},...,\lambda_{r},0,0,...,0),0\leq r<n$, and
$Q=(q_{ij})_{i,j=1,...,n}$ satisfies $QQ^{T}=Q^{T}Q=I_{n}$. We define
$\tau=Q^{T}\xi=(\tau_{1},...,\tau_{r},0,...,0)^{T}$ almost surely. For any
continuously differentiable function $f$ such that
$\mathbf{E}\sum_{i=1}^{n}|\frac{\partial f}{\partial x_{i}}|<\infty$, lemma 2
in [40] implies

$\mathbf{E}\xi_{i}f(\xi)=\sum_{j=1}^{r}q_{ij}\mathbf{E}\tau_{j}f(Q\tau)=\sum_{j=1}^{r}\sum_{k=1}^{r}\sum_{l=1}^{n}q_{ij}q_{lk}\mathbf{E}(\tau_{j}\tau_{k})\times\mathbf{E}\frac{\partial
f}{\partial x_{l}}(Q\tau)=\sum_{l=1}^{n}\mathbf{E}\frac{\partial f}{\partial
x_{l}}(\xi)\mathbf{E}\xi_{i}\xi_{l}$ (36)

This observation assures that (35) works for degenerated $\xi$.

###### Proof.

Since
$|\epsilon_{i}|=\max(\epsilon_{i},-\epsilon_{i})\Rightarrow\max_{i=1,2,...,n}|\epsilon_{i}|=\max(\max_{i=1,...,n}\epsilon_{i},\max_{i=1,2,...,n}-\epsilon_{i})$,
and $-\epsilon$ has the same distribution as $\epsilon$,

$\displaystyle\sup_{x\in\mathbf{R}}\left(Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x+\delta\right)-Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x\right)\right)\leq\sup_{x\in\mathbf{R}}Prob\left(x<\max_{i=1,...,n}\epsilon_{i}\leq
x+\delta\right)$ (37)
$\displaystyle+\sup_{x\in\mathbf{R}}Prob\left(x<\max_{i=1,2,...,n}-\epsilon_{i}\leq
x+\delta\right)\leq
2\sup_{x\in\mathbf{R}}Prob\left(|\max_{i=1,...,n}\epsilon_{i}-x|\leq\delta\right)$

From theorem 3 and (18), (19) in Chernozhukov et al.[60](also see lemma 3 in
Zhang and Politis [28]), we define
$\underline{\sigma}=\min_{i=1,2,...,n}\sigma_{ii}$ and
$\overline{\sigma}=\max_{i=1,2,...,n}\sigma_{ii}$,

$\displaystyle\sup_{x\in\mathbf{R}}Prob\left(|\max_{i=1,2,...,n}\epsilon_{i}-x|\leq\delta\right)\leq\frac{\sqrt{2}\delta}{\underline{\sigma}}\left(\sqrt{\log(n)}+\sqrt{\max(1,\log(\underline{\sigma})-\log(\delta))}\right)$
(38)
$\displaystyle+\frac{4\sqrt{2}\delta}{\underline{\sigma}}\times\left(\frac{\overline{\sigma}}{\underline{\sigma}}\sqrt{\log(n)}+2+\frac{\overline{\sigma}}{\underline{\sigma}}\sqrt{\max(0,\log(\underline{\sigma})-\log(\delta))}\right)$
$\displaystyle\leq\frac{\sqrt{2}\delta}{c_{0}}\left(\sqrt{\log(n)}+\sqrt{1+|\log(c_{0})|+|\log(C_{0})|}+\sqrt{|\log(\delta)|}\right)$
$\displaystyle+\frac{4\sqrt{2}\delta
C_{0}}{c_{0}^{2}}\left(\sqrt{\log(n)}+2+\sqrt{|\log(c_{0})|+|\log(C_{0})|}+\sqrt{|\log(\delta)|}\right)$
$\displaystyle\leq\left(\frac{\sqrt{2\times(1+|\log(c_{0})|+|\log(C_{0})|)}}{c_{0}}+\frac{4\sqrt{2}C_{0}}{c_{0}^{2}}(2+\sqrt{|\log(c_{0})|+|\log(C_{0})|})\right)$
$\displaystyle\times\delta\left(\sqrt{\log(n)}+1+\sqrt{|\log(\delta)|}\right)$

and we prove (34).

Without loss of generality, assume $\epsilon$ is independent of $\xi$. Similar
to Chernozhukov et al.[60], for any $0\leq t\leq 1$, we define random
variables $Z_{i}(t)=\sqrt{t}\epsilon_{i}+\sqrt{1-t}\xi_{i}$. According to
(31), (36), and theorem 2.27 in Folland [38],

$\displaystyle\mathbf{E}h_{\tau,\psi,x}(\epsilon_{1},...,\epsilon_{n})-\mathbf{E}h_{\tau,\psi,x}(\xi_{1},...,\xi_{n})$
(39)
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}\int_{0}^{1}t^{-1/2}\mathbf{E}(\frac{\partial
h_{\tau,\psi,x}(Z_{1}(t),...Z_{n}(t))}{\partial
x_{i}}\epsilon_{i})dt-\frac{1}{2}\sum_{i=1}^{n}\int_{0}^{1}(1-t)^{-1/2}\mathbf{E}(\frac{\partial
h_{\tau,\psi,x}(Z_{1}(t),...,Z_{n}(t))}{\partial x_{i}}\xi_{i})dt$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}\sum_{k=1}^{n}(\sigma_{ik}-\sigma^{*}_{ik})\int_{0}^{1}\mathbf{E}\frac{\partial^{2}h_{\tau,\psi,x}(Z_{1}(t),...,Z_{n}(t))}{\partial
x_{i}\partial x_{k}}dt$
$\displaystyle\Rightarrow|\mathbf{E}h_{\tau,\psi,x}(\epsilon_{1},...,\epsilon_{n})-\mathbf{E}h_{\tau,\psi,x}(\xi_{1},...,\xi_{n})|\leq\Delta\times(g_{*}\psi^{2}+g_{*}\psi\tau)$

For any $x\in\mathbf{R}$ and given $\tau,\psi>0$, define
$t=\frac{1}{\psi}+\frac{\log(2n)}{\tau}$, then

$\displaystyle Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x\right)-Prob\left(\max_{i=1,2,...,n}|\xi_{i}|\leq x\right)$ (40)
$\displaystyle\leq Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x-t\right)+Ct(\sqrt{\log(n)}+\sqrt{|\log(t)|}+1)-Prob\left(\max_{i=1,2,...,n}|\xi_{i}|\leq
x\right)$
$\displaystyle\leq\mathbf{E}h_{\tau,\psi,x-\frac{1}{\psi}}(\epsilon_{1},...,\epsilon_{n})-\mathbf{E}h_{\tau,\psi,x-\frac{1}{\psi}}(\xi_{1},...,\xi_{n})+Ct(\sqrt{\log(n)}+\sqrt{|\log(t)|}+1)$
$\displaystyle Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x\right)-Prob\left(\max_{i=1,2,...,n}|\xi_{i}|\leq x\right)$
$\displaystyle\geq Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x+t\right)-Ct(\sqrt{\log(n)}+\sqrt{|\log(t)|}+1)-Prob\left(\max_{i=1,2,...,n}|\xi_{i}|\leq
x\right)$
$\displaystyle\geq\mathbf{E}h_{\tau,\psi,x+\frac{\log(2n)}{\tau}}(\epsilon_{1},...,\epsilon_{n})-\mathbf{E}h_{\tau,\psi,x+\frac{\log(2n)}{\tau}}(\xi_{1},...,\xi_{n})-Ct(\sqrt{\log(n)}+\sqrt{|\log(t)|}+1)$
$\displaystyle\Rightarrow\sup_{x\in\mathbf{R}}|Prob\left(\max_{i=1,2,...,n}|\epsilon_{i}|\leq
x\right)-Prob\left(\max_{i=1,2,...,n}|\xi_{i}|\leq x\right)|$
$\displaystyle\leq\Delta\times(g_{*}\psi^{2}+g_{*}\psi\tau)+Ct(\sqrt{\log(n)}+\sqrt{|\log(t)|}+1)$

We choose $\tau=\psi=\left(1+\log^{3/2}(n)\right)/\Delta^{1/3}$, then
$\exists$ a constant $C_{1}>0$ such that
$\frac{1}{C_{1}}\frac{\Delta^{1/3}}{1+\log^{1/2}(n)}\leq
t=\Delta^{1/3}\left(\frac{1+\log(2)}{1+\log^{3/2}(n)}\right)+\frac{\Delta^{1/3}\times\log(n)}{1+\log^{3/2}(n)}\leq\frac{C_{1}\Delta^{1/3}}{1+\log^{1/2}(n)}$
for $n=1,2,...$; and we prove (35). ∎

Lemma 3 approximates the distribution of the linear combinations of errors
$\epsilon$ by a multivariate normal random vector. This lemma can be used to
derive the Gaussian approximation theorem for the estimator $\widehat{\beta}$.

######  Lemma 3 (Gaussian approximation theorem)

Suppose $\epsilon_{i},i=1,2,...,n$ satisfy assumption 2, and
$(\gamma_{ij})_{i=1,...,p_{1},j=1,...,n}\in\mathbf{R}^{p_{1}\times n}$ has
rank $p_{1}$, $p_{1}=O(1),\ p_{1}\leq n$. Besides, suppose $\exists$ constants
$0<c_{\gamma}\leq C_{\gamma}<\infty$ such that
$c_{\gamma}\leq\sum_{j=1}^{n}\gamma^{2}_{ij}\leq C_{\gamma}$ for $\forall
i=1,2,...,p_{1}$ and $n=1,2,...$. Define
$\gamma^{*}=\max_{i=1,...,p_{1},j=1,...,n}|\gamma_{ij}|$ and assume
$\gamma^{*}=o(n^{-1/4}\log^{-z}(n))$, here
$z=\max(\frac{9}{2},\frac{3\alpha_{\epsilon}}{2\alpha_{\epsilon}-2})$. Then

$\sup_{x\in\mathbf{R}}|Prob\left(\max_{i=1,2,...,p_{1}}|\sum_{j=1}^{n}\gamma_{ij}\epsilon_{j}|\leq
x\right)-Prob\left(\max_{i=1,2,...,p_{1}}|\sum_{j=1}^{n}\gamma_{ij}\xi_{j}|\leq
x\right)|=o(1)$ (41)

Here $\xi=(\xi_{1},...,\xi_{n})^{T}$ has multivariate normal distribution,
$\mathbf{E}\xi=0$ and $\mathbf{E}\xi\xi^{T}=\Sigma$, the covariance matrix of
$\epsilon$.

###### Proof.

For any given $\psi,\tau>0$, define
$t=\frac{1}{\psi}+\frac{\log(2p_{1})}{\tau}$. Assumption 2, (27), and (40)
imply

$\displaystyle\sup_{x\in\mathbf{R}}|Prob\left(\max_{i=1,2,...,p_{1}}|\sum_{j=1}^{n}\gamma_{ij}\epsilon_{j}|\leq
x\right)-Prob\left(\max_{i=1,2,...,p_{1}}|\sum_{j=1}^{n}\gamma_{ij}\xi_{j}|\leq
x\right)|$ (42)
$\displaystyle\leq\sup_{x\in\mathbf{R}}|\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\epsilon_{j},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\epsilon_{j})-\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\xi_{j},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\xi_{j})|$
$\displaystyle+O\left(t(\sqrt{\log(p_{1})}+\sqrt{|\log(t)|}+1)\right)$

For any integer $s>0$, (26) implies

$\displaystyle\|\sum_{j=1}^{n}\gamma_{ij}(\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s})\|_{m}\leq\sum_{k=s}^{\infty}\|\sum_{j=1}^{n}\gamma_{ij}(\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,k+1}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,k})\|_{m}$
(43) $\displaystyle\leq
C\sqrt{\sum_{j=1}^{n}\gamma^{2}_{ij}}\sum_{k=s}^{\infty}\max_{j=1,2,...,n}\delta_{j,k,m}$
$\displaystyle\Rightarrow\max_{i=1,2,...,p_{1}}\|\sum_{j=1}^{n}\gamma_{ij}(\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s})\|_{m}=O\left(\frac{1}{(s+1)^{\alpha_{\epsilon}}}\right)$

Here $C$ is a constant. Therefore,

$\displaystyle\sup_{x\in\mathbf{R}}|\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\epsilon_{j},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\epsilon_{j})-\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s})|$
(44) $\displaystyle\leq
g_{*}\psi\mathbf{E}\max_{i=1,2,...,p_{1}}|\sum_{j=1}^{n}\gamma_{ij}(\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s})|$
$\displaystyle\leq
g_{*}\psi\sum_{i=1}^{p_{1}}\|\sum_{j=1}^{n}\gamma_{ij}(\epsilon_{j}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s})\|_{m}=O\left(\frac{\psi}{(s+1)^{\alpha_{\epsilon}}}\right)$

For any integer $k>s$, define the big block
$S_{il}=\sum_{j=(l-1)\times(k+s)+1}^{(l-1)\times(k+s)+k}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}$
and the small block
$s_{il}=\sum_{j=(l-1)\times(k+s)+k+1}^{l\times(k+s)}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}$
for $l=1,2,...,r=\lfloor\frac{n}{k+s}\rfloor$. Here $\lfloor x\rfloor$ denotes
the largest integer that is smaller than or equal to $x$. If $(k+s)r+k\geq n$,
we define
$S_{i(r+1)}=\sum_{j=(k+s)r+1}^{n}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}$
and $s_{i(r+1)}=0$. If $(k+s)r+k<n$, we define
$S_{i(r+1)}=\sum_{j=(k+s)r+1}^{j=(k+s)r+k}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}$
and
$s_{i(r+1)}=\sum_{j=(k+s)r+k+1}^{n}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}$.
$(S_{1l},...,S_{p_{1}l})^{T},\ l=1,2,...,r+1$ are independent with each other;
$(s_{1l},...,s_{p_{1}l})^{T},l=1,2,...,r+1$ are independent with each other;
and
$\sum_{j=1}^{n}\gamma_{ij}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}=\sum_{j=1}^{r+1}(S_{ij}+s_{ij})$.
For
$\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s}=\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,0}+\sum_{l=1}^{s}\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,l}-\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,l-1}$,
(26) and (27) imply

$\displaystyle|\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{r+1}S_{1j}+s_{1j},...,\sum_{j=1}^{r+1}S_{p_{1}j}+s_{p1j})-\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{r+1}S_{1j},...,\sum_{j=1}^{r+1}S_{p_{1}j})|$
(45) $\displaystyle\leq
g_{*}\psi\sum_{i=1}^{p_{1}}\|\sum_{j=1}^{r+1}s_{ij}\|_{m}=O\left(\psi\max_{i=1,...,p_{1}}\sqrt{\sum_{j=1}^{r+1}\sum_{l=(j-1)\times(k+s)+k+1}^{j\times(k+s)}\gamma^{2}_{il}}\right)$
$\displaystyle=O\left(\psi\times\gamma^{*}\sqrt{\frac{ns}{k}}\right)$

Here we define $\gamma_{ij}=0$ for $i=1,2,...,p_{1},j>n$.

Define $\xi_{j,s},j=1,2,...,n$ as the joint normal random variables with mean
$0$ and the same covariance matrix as
$\mathbf{E}\epsilon_{j}|\mathcal{F}_{j,s},\ j=1,2,...,n$; and are independent
of $e_{i},i=...,-1,0,1,...$. Define
$S^{*}_{il}=\sum_{j=(l-1)(k+s)+1}^{(l-1)(k+s)+k}\gamma_{ij}\xi_{j,s}$ for
$l=1,...,r$ and
$S^{*}_{ir+1}=\sum_{j=r(k+s)+1}^{\min(n,r(k+s)+k)}\gamma_{ij}\xi_{j,s}$.
Define $H_{il}=\sum_{j=1}^{l-1}S_{ij}+\sum_{j=l+1}^{r+1}S^{*}_{ij}$, then
$H_{il}+S_{il}=H_{il+1}+S_{il+1}^{*}$. $S^{*}_{il},i=1,...,p_{1},l=1,...,r+1$
has joint normal distribution(we treat a constant as a degenerated normal
random variable) and $\mathbf{E}S^{*}_{il}S^{*}_{jk}=\mathbf{E}S_{il}S_{jk}=0$
for $l\neq k$, so $(S^{*}_{1l},...,S^{*}_{p_{1}l})^{T}$ are independent with
each other. From Taylor’s theorem,

$\displaystyle|\mathbf{E}\left(h_{\tau,\psi,x}(H_{1j}+S_{1j},...,H_{p_{1}j}+S_{p_{1}j})-h_{\tau,\psi,x}(H_{1j}+S^{*}_{1j},...,H_{p_{1}j}+S^{*}_{p_{1}j})\right)|H_{1j},...,H_{p_{1}j}|$
(46) $\displaystyle\leq|\sum_{l=1}^{p_{1}}\frac{\partial
h_{\tau,\psi,x}(H_{1j},...,H_{p_{1}j})}{\partial
x_{l}}\mathbf{E}(S_{lj}-S^{*}_{lj})|$
$\displaystyle+|\frac{1}{2}\sum_{l_{1}=1}^{p_{1}}\sum_{l_{2}=1}^{p_{1}}\frac{\partial^{2}h_{\tau,\psi,x}(H_{1j},...,H_{p_{1}j})}{\partial
x_{l_{1}}\partial
x_{l_{2}}}\mathbf{E}(S_{l_{1}j}S_{l_{2}j}-S^{*}_{l_{1}j}S^{*}_{l_{2}j})|$
$\displaystyle+g_{*}(\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\mathbf{E}(\max_{i=1,2,...,p_{1}}|S_{ij}|^{3}+\max_{i=1,2,...,p_{1}}|S_{ij}^{*}|^{3})$
$\displaystyle\leq
g_{*}(\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\sum_{i=1}^{p_{1}}\left(\|S_{ij}\|_{m}^{3}+\|S^{*}_{ij}\|_{m}^{3}\right)$
$\displaystyle\leq
C\left((\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\times\max_{i=1,2,...,p_{1}}\left(\sum_{l=(j-1)(k+s)+1}^{(j-1)(k+s)+k}\gamma^{2}_{il}\right)^{3/2}\right)$

Here $C$ is a constant. Since
$\max_{i=1,...,p_{1}}\sum_{j=1}^{n}\gamma^{2}_{ij}\leq C_{\gamma}$,

$\displaystyle\sup_{x\in\mathbf{R}}|\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{r+1}S_{1j},...,\sum_{j=1}^{r+1}S_{p_{1}j})-\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{r+1}S_{1j}^{*},...,\sum_{j=1}^{r+1}S_{p_{1}j}^{*})|$
(47)
$\displaystyle\leq\sum_{j=1}^{r+1}\sup_{x\in\mathbf{R}}|\mathbf{E}h_{\tau,\psi,x}(H_{1j}+S_{1j},...,H_{p_{1}j}+S_{p_{1}j})-\mathbf{E}h_{\tau,\psi,x}(H_{1j}+S^{*}_{1j},...,H_{p_{1}j}+S^{*}_{p_{1}j})|$
$\displaystyle\leq
C(\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\sum_{j=1}^{r+1}\sum_{i=1}^{p_{1}}\left(\sum_{l=(j-1)(k+s)+1}^{(j-1)(k+s)+k}\gamma^{2}_{il}\right)^{3/2}$
$\displaystyle=O\left((\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\times\max_{i=1,...,p_{1},j=1,...,r+1}\sqrt{\sum_{l=(j-1)(k+s)+1}^{(j-1)(k+s)+k}\gamma^{2}_{il}}\right)$

which has order
$O\left((\psi^{3}+\tau\psi^{2}+\psi\tau^{2})\times\gamma^{*}\sqrt{k}\right)$.
We define $s^{*}_{il}=\sum_{j=(l-1)(k+s)+k+1}^{l(k+s)}\gamma_{ij}\xi_{j,s}$
for $l=1,2,...,r$, and
$s^{*}_{ir+1}=\sum_{j=r(k+s)+k+1}^{n}\gamma_{ij}\xi_{j,s}$. For $s^{*}_{il}$
has normal distribution,

$\displaystyle|\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{r+1}S^{*}_{1j},...,\sum_{j=1}^{r+1}S^{*}_{p_{1}j})-\mathbf{E}h_{\tau,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\xi_{j,s},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\xi_{j,s})|$
(48) $\displaystyle\leq
g_{*}\psi\sum_{i=1}^{p_{1}}\|\sum_{j=1}^{r+1}s_{ij}^{*}\|_{m}=O\left(\psi\max_{i=1,2,...,p_{1}}\|\sum_{j=1}^{r+1}s_{ij}\|_{2}\right)$

According to (45), this has order
$O\left(\psi\times\gamma^{*}\sqrt{\frac{ns}{k}}\right)$. From (39), (28), and
0.9.7 in [33],

$\displaystyle|\mathbf{E}h_{\beta,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\xi_{j,s},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\xi_{j,s})-\mathbf{E}h_{\beta,\psi,x}(\sum_{j=1}^{n}\gamma_{1j}\xi_{j},...,\sum_{j=1}^{n}\gamma_{p_{1}j}\xi_{j})|$
(49)
$\displaystyle\leq(g_{*}\psi^{2}+g_{*}\psi\tau)\times\max_{i,j=1,2,...,p_{1}}|\sum_{k_{1}=1}^{n}\sum_{k_{2}=1}^{n}\gamma_{ik_{1}}\gamma_{jk_{2}}(\mathbf{E}(\mathbf{E}\epsilon_{k_{1}}|\mathcal{F}_{k_{1},s}\times\mathbf{E}\epsilon_{k_{2}}|\mathcal{F}_{k_{2},s})-\sigma_{k_{1}k_{2}})|$
$\displaystyle\leq(g_{*}\psi^{2}+g_{*}\psi\tau)\max_{i,j=1,...,p_{1}}\sum_{|k_{1}-k_{2}|\leq
s}|\gamma_{ik_{1}}\gamma_{jk_{2}}|(\|\mathbf{E}\epsilon_{k_{1}}|\mathcal{F}_{k_{1},s}\|_{2}\|\mathbf{E}\epsilon_{k_{2}}|\mathcal{F}_{k_{2},s}-\epsilon_{k_{2}}\|_{2}+\|\mathbf{E}\epsilon_{k_{1}}|\mathcal{F}_{k_{1},s}-\epsilon_{k_{1}}\|_{2}\|\epsilon_{k_{2}}\|_{2})$
$\displaystyle+C(g_{*}\psi^{2}+g_{*}\psi\tau)\max_{i,j=1,...,p_{1}}\sum_{|k_{1}-k_{2}|>s}\frac{|\gamma_{ik_{1}}\gamma_{jk_{2}}|}{(1+|k_{1}-k_{2}|)^{\alpha_{\epsilon}}}$
$\displaystyle=O\left((\psi^{2}+\psi\tau)\times\left(\max_{i=1,2,...,p_{1}}\sum_{k=1}^{n}\gamma_{ik}^{2}\right)\times\left(\frac{1}{(1+s)^{\alpha_{\epsilon}-1}}\right)\right)$

Here $\sigma_{ij}=\mathbf{E}\epsilon_{i}\epsilon_{j}$.

Define $V$ such that $\frac{1}{V}=\gamma^{*}\times n^{1/4}\log^{z}(n)\to 0$.
Choose $k=\lfloor\sqrt{n}\rfloor$,
$\psi=\tau=V^{(\alpha_{\epsilon}-1)/(3\alpha_{\epsilon}+3)}$, and $s=\lfloor
V^{2/(\alpha_{\epsilon}+1)}\log^{3/(\alpha_{\epsilon}-1)}(n)\rfloor$. Then we
prove (41). ∎

###### proof of theorem 12.

$\displaystyle\max_{i=1,2,...,p}\sum_{l=1}^{n}\left(\sum_{j=1}^{p}q_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\right)^{2}$
(50)
$\displaystyle=\max_{i=1,2,...,p}\sum_{j=1}^{p}q^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}\leq\frac{4}{\lambda^{2}_{p}}$

Define $s=(s_{1},...,s_{p})^{T}=Q^{T}\beta$, from lemma 32

$\displaystyle\widetilde{\beta}_{i}-\beta_{i}=-\rho^{2}_{n}\sum_{j=1}^{p}\frac{q_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}+\sum_{j=1}^{p}\sum_{l=1}^{n}q_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\epsilon_{l}$
(51)
$\displaystyle\Rightarrow\max_{i=1,...,p}|\widetilde{\beta}_{i}-\beta_{i}|$
$\displaystyle\leq\max_{i=1,...,p}\sqrt{\sum_{j=1}^{p}q^{2}_{ij}}\times\frac{\rho^{2}_{n}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}{\lambda^{4}_{p}}+\max_{i=1,2,...,p}|\sum_{j=1}^{p}\sum_{l=1}^{n}q_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\epsilon_{l}|$
$\displaystyle=O_{p}\left(n^{1/m-\eta}\right)$

Therefore,

$\displaystyle
Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)\leq
Prob\left(\min_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|\leq
b_{n}\right)+Prob\left(\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|>b_{n}\right)$
(52) $\displaystyle\leq
Prob\left(\min_{i\in\mathcal{N}_{b_{n}}}|\beta_{i}|-\max_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}-\beta_{i}|\leq
b_{n}\right)$
$\displaystyle+Prob\left(\max_{i\not\in\mathcal{N}_{b_{n}}}|\beta_{i}|+\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}-\beta_{i}|>b_{n}\right)$

From assumption 4,
$\min_{i\in\mathcal{N}_{b_{n}}}|\beta_{i}|-b_{n}\geq(1/c_{b}-1)b_{n}$ and
$b_{n}-\max_{i\not\in\mathcal{N}_{b_{n}}}|\beta_{i}|\geq(1-c_{b})b_{n}$. So
$Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)=o(1)$
is proved by (51).

If $\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$, we define
$\\{c_{ij}\\}_{i=1,...,p_{1},j=1,...,p}$ as in section 2 with $b=b_{n}$ and
$\gamma=(\gamma_{1},...,\gamma_{p_{1}})^{T}=M\beta$. Then

$\displaystyle\widehat{\gamma}_{i}-\gamma_{i}=\sum_{j\in\mathcal{N}_{b_{n}}}m_{ij}(\widetilde{\beta}_{j}-\beta_{j})-\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}$
(53)
$\displaystyle=-\rho^{2}_{n}\sum_{j=1}^{p}\frac{c_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}+\sum_{j=1}^{p}\sum_{l=1}^{n}c_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\epsilon_{l}-\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}$

For
$\max_{i=1,...,p_{1}}|\rho^{2}_{n}\sum_{j=1}^{p}\frac{c_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}|\leq\frac{\rho^{2}_{n}\sqrt{C_{\mathcal{M}}}}{\lambda^{4}_{p}}\times{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}=O(n^{\alpha_{\beta}-2\delta})$,
and

$\displaystyle\max_{i=1,2,...,p_{1}}\sum_{l=1}^{n}\left(\sum_{j=1}^{p}c_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\right)^{2}$
(54)
$\displaystyle=\max_{i=1,2,...,p_{1}}\sum_{j=1}^{p}c^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)^{2}\leq\frac{4C_{\mathcal{M}}}{\lambda^{2}_{p}}$

lemma 32 implies
$\max_{i=1,2,...,p_{1}}|\widehat{\gamma}_{i}-\gamma_{i}|=O_{p}(n^{-\eta})$.

If $\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$, we define
$u_{ij}=\sum_{k\in\mathcal{N}_{b_{n}}}x_{ik}q_{kj},i=1,...,n,j=1,...,p$,

$\displaystyle\widehat{\epsilon}_{i}-\epsilon_{i}=-\sum_{j\in\mathcal{N}_{b_{n}}}x_{ij}(\widetilde{\beta}_{j}-\beta_{j})+\sum_{j\not\in\mathcal{N}_{b_{n}}}x_{ij}\beta_{j}$
(55)
$\displaystyle=\rho_{n}^{2}\sum_{j=1}^{p}\frac{u_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}-\sum_{j=1}^{p}\sum_{l=1}^{n}u_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\epsilon_{l}+\sum_{j\not\in\mathcal{N}_{b_{n}}}x_{ij}\beta_{j}$

For
$\max_{i=1,...,n}\rho_{n}^{2}|\sum_{j=1}^{p}\frac{u_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}|\leq\frac{\rho^{2}_{n}{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\beta\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{2}}{\lambda^{4}_{p}}\times\max_{i=1,..,n}\sqrt{\sum_{j=1}^{p}u^{2}_{ij}}=O(n^{\alpha_{\beta}-2\delta}\sqrt{|\mathcal{N}_{b_{n}}|})$,
and

$\displaystyle\max_{i=1,2,...,n}\sum_{l=1}^{n}\left(\sum_{j=1}^{p}u_{ij}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)p_{lj}\right)^{2}$
(56)
$\displaystyle\leq\frac{4}{\lambda^{2}_{p}}\max_{i=1,2,...,n}\sum_{j\in\mathcal{N}_{b_{n}}}x^{2}_{ij}=O\left(|\mathcal{N}_{b_{n}}|\times
n^{-2\eta}\right)$

from lemma 32,

$\max_{i=1,...,n}|\widehat{\epsilon}_{i}-\epsilon_{i}|=o_{p}(n^{(1/m-\eta)/2})+O_{p}(n^{\alpha_{\beta}-2\delta}\sqrt{|\mathcal{N}_{b_{n}}|})+O_{p}(n^{1/m-\eta}\times\sqrt{|\mathcal{N}_{b_{n}}|})$
(57)

and we prove (12). ∎

###### proof of theorem 2.

Define
$v_{il}=\frac{1}{\tau_{i}}\sum_{j=1}^{p}c_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)p_{lj}$,
then $v_{il}=0$ if $i\not\in\mathcal{M}$. If
$\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$, then
$\widehat{\tau}_{i}=\tau_{i}$. From (53) and Cauchy-Schwarz inequality,

$\displaystyle\max_{i=1,...,p_{1}}|\frac{1}{\widehat{\tau}_{i}}\left(\widehat{\gamma}_{i}-\sum_{j=1}^{p}m_{ij}\beta_{j}\right)-\sum_{l=1}^{n}v_{il}\epsilon_{l}|$
(58)
$\displaystyle\leq\max_{i\in\mathcal{M}}\frac{\rho^{2}_{n}}{\tau_{i}}|\sum_{j=1}^{p}\frac{c_{ij}s_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}|+\max_{i=1,...,p_{1}}\frac{1}{\tau_{i}}|\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}|$
$\displaystyle\leq\max_{i\in\mathcal{M}}\frac{\rho^{2}_{n}}{\tau_{i}}\sqrt{\sum_{j=1}^{p}\frac{c^{2}_{ij}\lambda^{2}_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}}\sqrt{\sum_{j=1}^{p}\frac{s^{2}_{j}}{\lambda^{2}_{j}(\lambda^{2}_{j}+\rho_{n})^{2}}}+\sqrt{n}\max_{i=1,...,p_{1}}|\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}|$
$\displaystyle\leq\frac{\rho^{2}_{n}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}{\lambda^{3}_{p}}+\sqrt{n}\max_{i=1,...,p_{1}}|\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}|$

which has order $O(n^{\alpha_{\beta}+\eta-2\delta})+o(1)$. Besides, if
$i\in\mathcal{M}$, then for sufficiently large $n$

$\displaystyle\sum_{l=1}^{n}v^{2}_{il}=\frac{1}{\tau^{2}_{i}}\sum_{j=1}^{p}c^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)^{2}\leq
1$ (59)
$\displaystyle\sum_{l=1}^{n}v^{2}_{il}=\frac{1}{1+\frac{1}{n\sum_{j=1}^{p}c^{2}_{ij}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)^{2}}}\geq\frac{1}{1+\frac{1}{n\sum_{j=1}^{p}c^{2}_{ij}/(4\lambda^{2}_{j})}}\geq\frac{1}{1+\frac{4C^{2}_{\lambda}}{c_{\mathcal{M}}}}>0$

Assumption 6 implies that the matrix
$\\{v_{il}\\}_{i\in\mathcal{M},l=1,...,n}$ has rank $|\mathcal{M}|$. Define
$t=\frac{\rho^{2}_{n}{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\beta\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{2}}{\lambda^{3}_{p}}+\max_{i=1,...,p_{1}}\sqrt{n}|\sum_{j\not\in\mathcal{N}_{b_{n}}}m_{ij}\beta_{j}|$,
lemma 35 and lemma 3 imply

$\displaystyle
Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\sum_{j=1}^{p}m_{ij}\beta_{j}|}{\widehat{\tau}_{i}}\leq
x\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x\right)$ (60) $\displaystyle\leq
Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)+|Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\epsilon_{l}|\leq
x+t\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x+t\right)|$ $\displaystyle+Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})$
$\displaystyle
Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\sum_{j=1}^{p}m_{ij}\beta_{j}|}{\widehat{\tau}_{i}}\leq
x\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x\right)$ $\displaystyle\geq-
Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)-|Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\epsilon_{l}|\leq
x-t\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x-t\right)|$ $\displaystyle-Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})$
$\displaystyle\Rightarrow\sup_{x\in\mathbf{R}}|Prob\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}_{i}-\sum_{j=1}^{p}m_{ij}\beta_{j}|}{\widehat{\tau}_{i}}\leq
x\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x\right)|$ $\displaystyle\leq
Prob\left(\widehat{\mathcal{N}}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)+\sup_{x\in\mathbf{R}}|Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\epsilon_{l}|\leq
x\right)-Prob\left(\max_{i\in\mathcal{M}}|\sum_{l=1}^{n}v_{il}\xi_{l}|\leq
x\right)|$ $\displaystyle+Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})$

and we prove (14). Here $C$ is a constant. ∎

##  Appendix C Proof of theorem 3

Errors in the linear model (2) can be dependent, non-stationary and
heteroskedastic. So it is hopeless to estimate the errors’ variance and
covariance. However, it is still possible to make a consistent estimation of
the estimator $\widehat{\gamma}=M\widehat{\beta}$’s variances. Lemma 62
justifies this idea.

###### Lemma 4

Suppose random variables $\epsilon_{i},i=1,2,...,n$ satisfy assumption 2;
$(\gamma_{ij})_{i=1,...,p_{1},j=1,...,n}\in\mathbf{R}^{p_{1}\times n}$ satisfy
conditions in lemma 3; the kernel function $K$ satisfies assumption 7; and
$k_{n}>0$ is a chosen bandwidth. Then

$\displaystyle\max_{i,j=1,...,p_{1}}|\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}\gamma_{is_{1}}\gamma_{js_{2}}K\left(\frac{s_{1}-s_{2}}{k_{n}}\right)\epsilon_{s_{1}}\epsilon_{s_{2}}-\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}\gamma_{is_{1}}\gamma_{js_{2}}\sigma_{s_{1}s_{2}}|$
(61) $\displaystyle=o_{p}(k_{n}\times n^{-1/4}\log^{-z}(n))+O_{p}(v_{n})$

Here $\sigma_{s_{1}s_{2}},s_{1},s_{2}=1,...,n$ are defined in section 2 and
$z=\max(\frac{9}{2},\frac{3\alpha_{\epsilon}}{2\alpha_{\epsilon}-2})$,

$v_{n}=\begin{cases}k_{n}^{1-\alpha_{\epsilon}}\ \text{if
}1<\alpha_{\epsilon}<2\\\ \log(k_{n})/k_{n}\ \text{if }\alpha_{\epsilon}=2\\\
1/k_{n}\ \text{if }\alpha_{\epsilon}>2\end{cases}$ (62)

###### Proof.

From (28) and section 0.9.7 in [33], there exists a constant $C>0$ such that

$\displaystyle|\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}\left(1-K\left(\frac{s_{1}-s_{2}}{k_{n}}\right)\right)\gamma_{is_{1}}\gamma_{js_{2}}\sigma_{s_{1}s_{2}}|$
(63) $\displaystyle\leq
C\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}|\gamma_{is_{1}}\gamma_{js_{2}}|\times\frac{\left(1-K\left(\frac{s_{1}-s_{2}}{k_{n}}\right)\right)}{(1+|s_{1}-s_{2}|)^{\alpha_{\epsilon}}}$
$\displaystyle\leq
2C\sqrt{\sum_{s=1}^{n}\gamma^{2}_{is}}\times\sqrt{\sum_{s=1}^{n}\gamma^{2}_{js}}\times\sum_{s=0}^{\infty}\frac{1-K(s/k_{n})}{1+s^{\alpha_{\epsilon}}}$

$K$ is continuously differentiable, so

$\displaystyle\sum_{s=0}^{\infty}\frac{1-K(s/k_{n})}{1+s^{\alpha_{\epsilon}}}\leq\frac{\max_{x\in[0,1]}|K^{{}^{\prime}}(x)|}{k_{n}}\sum_{s=0}^{k_{n}}\frac{s}{1+s^{\alpha_{\epsilon}}}+\sum_{s=k_{n}+1}^{\infty}\frac{1}{1+s^{\alpha_{\epsilon}}}$
(64)
$\displaystyle=O\left(\frac{1}{k_{n}}\left(1+\int_{[1,k_{n}]}x^{1-\alpha_{\epsilon}}dx\right)+\int_{[k_{n},\infty)}x^{-\alpha_{\epsilon}}dx\right)$

we have
$\max_{i,j=1,...,p_{1}}|\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}\left(1-K\left(\frac{s_{1}-s_{2}}{k_{n}}\right)\right)\gamma_{is_{1}}\gamma_{js_{2}}\sigma_{s_{1}s_{2}}|=O(v_{n})$.

On the other hand, define
$\zeta_{i,k}=\epsilon_{i}\epsilon_{i+k}-\sigma_{ii+k}=h_{i,i+k}(...,e_{i+k-1},e_{i+k})$(i.e.,
$\zeta_{i,k}$ is $\mathcal{F}_{i+k}$ measurable), and define
$\zeta_{i,k,t}=h_{i,i+k}(...,e_{i+k-t-1},e_{i+k-t}^{\dagger},e_{i+k-t+1},...,e_{i+k})$.
Here $k\geq 0,i\geq 1,i+k\leq n,t\geq 0$. $e_{i},e^{\dagger}_{i},\sigma_{ij}$
are defined in section 2. Then $\mathbf{E}\zeta_{i,k}=0$ and

$\displaystyle\psi_{i,k,t,m/2}=\|\zeta_{i,k}-\zeta_{i,k,t}\|_{m/2}=\|\epsilon_{i}\epsilon_{i+k}-\epsilon_{i,t-k}\epsilon_{i+k,t}\|_{m/2}$
(65)
$\displaystyle\leq\|\epsilon_{i}-\epsilon_{i,t-k}\|_{m}\|\epsilon_{i+k}\|_{m}+\|\epsilon_{i,t-k}\|_{m}\|\epsilon_{i+k}-\epsilon_{i+k,t}\|_{m}$
$\displaystyle\Rightarrow\max_{i=1,...,n-k}\psi_{i,k,t,m/2}\leq
C\max_{i=1,...,n}\delta_{i,t-k,m}+C\max_{i=1,...,n}\delta_{i,t,m}\ \text{if
$t-k\geq 0$}$ $\displaystyle\text{ and }C\max_{i=1,...,n}\delta_{i,t,m}\
\text{if $t-k<0$}$

Here $\epsilon_{i,t}=\epsilon_{i}$ if $t<0$ and $C$ is a constant. For fixed
$i,j$, define

$N_{s,k,t}=\sum_{l=n-k+1-s}^{n-k}\gamma_{il}\gamma_{jl+k}\left(\mathbf{E}\zeta_{l,k}|\mathcal{F}_{l+k,t}-\mathbf{E}\zeta_{l,k}|\mathcal{F}_{l+k,t-1}\right)$,
then $N_{s,k,t}$ is $\mathcal{F}_{n,s+t-1}$ measurable. Apply $\pi-\lambda$
theorem to the $\lambda-$ system

$\\{A\in\mathcal{F}_{n,s+t-1}\Big{|}\mathbf{E}\left(\mathbf{E}\zeta_{n-k-s,k}|\mathcal{F}_{n-s,t}\right)\times\mathbf{1}_{A}=\mathbf{E}\left(\mathbf{E}\zeta_{n-k-s,k}|\mathcal{F}_{n-s,t-1}\right)\times\mathbf{1}_{A}\\}$
(66)

and the $\pi-$system $\\{A_{n}\times A_{n-1}\times...\times A_{n-s-t+1}\\}$,
$A_{i}$ is generated by $e_{i}$. Then we know that $N_{s,k,t},s=1,2,...,n-k$
form a martingale for any given $k,t$. From (25), (26) and (27)

$\displaystyle\|\sum_{l=1}^{n-k}\gamma_{il}\gamma_{jl+k}\zeta_{l,k}\|_{m/2}\leq\|\sum_{l=1}^{n-k}\gamma_{il}\gamma_{jl+k}\mathbf{E}\zeta_{l,k}|\mathcal{F}_{l+k,0}\|_{m/2}+\sum_{t=1}^{\infty}\|N_{n-k,k,t}\|_{m/2}$
(67) $\displaystyle\leq
C\sqrt{\sum_{l=1}^{n-k}\gamma^{2}_{il}\gamma_{jl+k}^{2}\|\zeta_{l,k}\|_{m/2}^{2}}+C\sqrt{\sum_{l=1}^{n-k}\gamma^{2}_{il}\gamma_{jl+k}^{2}}\times\sum_{t=1}^{\infty}\max_{l=1,...,n-k}\psi_{l,k,t,m/2}$

Here $C$ is a constant independent with $i,j$. Since
$\|\zeta_{l,k}\|_{m/2}\leq\|\epsilon_{i}\|_{m}\times\|\epsilon_{i+k}\|_{m}+|\sigma_{ii+k}|$
and

$\sum_{t=1}^{\infty}\max_{l=1,...,n-k}\psi_{l,k,t,m/2}\leq
2C\sum_{t=0}^{\infty}\max_{i=1,...,n}\delta_{i,t,m}\ \text{with $C$ a
constant}$ (68)

we have

$\max_{i,j=1,...,p_{1}}\frac{\|\sum_{l=1}^{n-k}\gamma_{il}\gamma_{jl+k}\zeta_{l,k}\|_{m/2}}{\sqrt{\sum_{l=1}^{n-k}\gamma^{2}_{il}\gamma^{2}_{jl+k}}}=O(1)$
(69)

Therefore

$\displaystyle\|\sum_{s_{1}=1}^{n}\sum_{s_{2}=1}^{n}\gamma_{is_{1}}\gamma_{js_{2}}K\left(\frac{s_{1}-s_{2}}{k_{n}}\right)\times(\epsilon_{s_{1}}\epsilon_{s_{2}}-\sigma_{s_{1}s_{2}})\|_{m/2}$
(70) $\displaystyle\leq
2\sum_{l=0}^{n-1}K\left(\frac{l}{k_{n}}\right)\|\sum_{s=1}^{n-l}\gamma_{is}\gamma_{js+l}\zeta_{s,l}\|_{m/2}$
$\displaystyle\leq
C\sum_{l=0}^{\infty}K\left(\frac{l}{k_{n}}\right)\times\max_{i=1,...,p_{1},j=1,...,n}|\gamma_{ij}|\
\text{with $C$ a constant}$

For $K$ is decreasing on $[0,\infty)$,

$\sum_{l=1}^{\infty}K(l/k_{n})\leq\sum_{l=1}^{\infty}\int_{[l-1,l]}K(x/k_{n})dx=k_{n}\int_{[0,\infty)}K(x)dx=O(k_{n})$
(71)

From (64) and (70), we prove (61). ∎

###### proof of theorem 3.

Recall $\mathbf{E}^{*}\cdot=\mathbf{E}\cdot|y$,
$Prob^{*}(\cdot)=Prob(\cdot|y)$, and define $\varepsilon^{*}_{i},i=1,...,n$ as
in algorithm 19. $\forall a=(a_{1},...,a_{n})^{T}\in\mathbf{R}^{n}$,
$\sum_{i=1}^{n}a_{i}\varepsilon^{*}_{i}$ has normal distribution. Section
0.9.7 in [33] and (71) implies

$\mathbf{E}^{*}|\sum_{i=1}^{n}a_{i}\varepsilon^{*}_{i}|^{m}=C\left(\sum_{i=1}^{n}\sum_{j=1}^{n}a_{i}a_{j}K\left(\frac{i-j}{k_{n}}\right)\right)^{m/2}=O\left(k_{n}^{m/2}\times{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|a\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}^{m}\right)$
(72)

Here $C=\mathbf{E}|Y|^{m}$, $Y$ has normal distribution with mean $0$ and
variance $1$. Therefore, $\forall\xi>0$,

$\displaystyle
Prob^{*}\left(\max_{i=1,...,p}|\sum_{j=1}^{p}\sum_{l=1}^{n}q_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\epsilon^{*}_{l}|>\xi\right)$
(73)
$\displaystyle\leq\frac{1}{\xi^{m}}\sum_{i=1}^{p}\mathbf{E}^{*}|\sum_{j=1}^{p}\sum_{l=1}^{n}q_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\epsilon^{*}_{l}|^{m}$
$\displaystyle=\frac{Ck_{n}^{m/2}}{\xi^{m}}\sum_{i=1}^{p}\left(\sum_{l=1}^{n}\widehat{\epsilon}_{l}^{2}\left(\sum_{j=1}^{p}q_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda_{j}^{2}+\rho_{n})^{2}}\right)\right)^{2}\right)^{m/2}$
$\displaystyle\leq\frac{Ck_{n}^{m/2}}{\xi^{m}}\times\max_{l=1,...,n}|\widehat{\epsilon}_{l}|^{m}\times\frac{2^{m}p}{\lambda^{m}_{p}}$

Here $C$ is a constant. From theorem 12,
$\max_{i=1,...,n}|\widehat{\epsilon}_{i}|^{m}\leq
2^{m}\max_{i=1,...,n}|\epsilon_{i}|^{m}+2^{m}\max_{i=1,...,n}|\widehat{\epsilon}_{i}-\epsilon_{i}|^{m}=O_{p}(n)$.
If $\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$, (51) implies
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq
2\sqrt{\sum_{i\in\mathcal{N}_{b_{n}}}(\widetilde{\beta}_{i}-\beta_{i})^{2}}+2{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}=O_{p}(n^{\alpha_{\beta}})$.
Define
$\widehat{s}=(\widehat{s}_{1},...,\widehat{s}_{p})^{T}=Q^{T}\widehat{\beta}$.
If $\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$, then for any given
$0<a<1$, there exists a constant $C>0$ such that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq
Cn^{\alpha_{\beta}}$ and $\max_{l=1,...,n}|\widehat{\epsilon}_{i}|^{m}\leq Cn$
with probability at least $1-a$. From (73), $\forall c>0$,

$\displaystyle\widetilde{\beta}^{*}_{i}-\widehat{\beta}_{i}=-\rho^{2}_{n}\sum_{j=1}^{p}\frac{q_{ij}\widehat{s}_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}+\sum_{j=1}^{p}\sum_{l=1}^{n}q_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\epsilon^{*}_{l}$
(74) $\displaystyle\Rightarrow
Prob^{*}\left(\max_{i=1,...,p}|\widetilde{\beta}^{*}_{i}-\widehat{\beta}_{i}|>2cn^{-\nu_{b}}\right)$
$\displaystyle\leq
Prob^{*}\left(\frac{n^{\nu_{b}}\rho^{2}_{n}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}{\lambda^{4}_{p}}>c\right)+\frac{C^{{}^{\prime}}}{c^{m}}\times\left(\sqrt{k_{n}}\times
n^{2/m-\eta+\nu_{b}}\right)^{m}$

Here $C^{{}^{\prime}}$ depends on $C$. If
$\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$,

$\displaystyle
Prob^{*}\left(\widehat{\mathcal{N}}_{b_{n}}^{*}\neq\mathcal{N}_{b_{n}}\right)\leq
Prob^{*}\left(\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}^{*}_{i}|>b_{n}\right)+Prob^{*}\left(\min_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}^{*}_{i}|\leq
b_{n}\right)$ (75) $\displaystyle\leq
Prob^{*}\left(\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}^{*}_{i}-\widehat{\beta}_{i}|>b_{n}-\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|\right)+Prob^{*}\left(\max_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}^{*}_{i}-\widehat{\beta}_{i}|\geq\min_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|-b_{n}\right)$

From assumption 4 and (51), with probability tending to $1$,
$\max_{i\not\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|\leq\frac{1+c_{b}}{2}b_{n}$
and
$\min_{i\in\mathcal{N}_{b_{n}}}|\widetilde{\beta}_{i}|\geq\frac{2b_{n}}{1+c_{b}}$.
(74) implies
$Prob^{*}\left(\widehat{\mathcal{N}}_{b_{n}}^{*}\neq\mathcal{N}_{b_{n}}\right)=o_{p}(1)$.

If
$\widehat{\mathcal{N}}_{b_{n}}^{*}=\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$,
define $c_{ij},i=1,...,p_{1},j=1,...,p$ as in section 2,

$\displaystyle\widehat{\gamma}^{*}_{i}-\widehat{\gamma_{i}}=\sum_{j\in\mathcal{N}_{b_{n}}}m_{ij}(\widetilde{\beta}^{*}_{j}-\widehat{\beta}_{j})$
(76)
$\displaystyle=-\rho_{n}^{2}\sum_{j=1}^{p}\frac{c_{ij}\widehat{s}_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}+\sum_{j=1}^{p}\sum_{l=1}^{n}c_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda^{2}_{j}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right)\epsilon^{*}_{l}$

Besides, $\widehat{\tau}^{*}_{i}=\tau_{i}$, and
$\widehat{c}^{*}_{ij}=c_{ij},i=1,...,p_{1},j=1,...,p$. Define
$v_{il}=\frac{1}{\tau_{i}}\sum_{j=1}^{p}c_{ij}p_{lj}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\rho_{n}}+\frac{\rho_{n}\lambda_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}\right),i=1,...,p_{1},l=1,...,n$;
and $z=(z_{1},...,z_{p_{1}})^{T}$ such that
$z_{i}=\sum_{l=1}^{n}v_{il}\epsilon^{*}_{l}$. We have
$\mathbf{E}^{*}z_{i}z_{j}=\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\widehat{\epsilon}_{l_{1}}\widehat{\epsilon}_{l_{2}}K\left(\frac{l_{1}-l_{2}}{k_{n}}\right)$.
Form section 0.9.7 in [33], lemma 62, and theorem 12, $\forall
i,j\in\mathcal{M}$,

$\displaystyle|\mathbf{E}^{*}z_{i}z_{j}-\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\sigma_{l_{1}l_{2}}|\leq|\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}K\left(\frac{l_{1}-l_{2}}{k_{n}}\right)\widehat{\epsilon}_{l_{1}}(\widehat{\epsilon}_{l_{2}}-\epsilon_{l_{2}})|$
(77)
$\displaystyle+|\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}K\left(\frac{l_{1}-l_{2}}{k_{n}}\right)(\widehat{\epsilon}_{l_{1}}-\epsilon_{l_{1}})\epsilon_{l_{2}}|$
$\displaystyle+|\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\epsilon_{l_{1}}\epsilon_{l_{2}}K\left(\frac{l_{1}-l_{2}}{k_{n}}\right)-\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\sigma_{l_{1}l_{2}}|$
$\displaystyle\leq
2\sum_{l=0}^{\infty}K\left(\frac{l}{k_{n}}\right)\times\left(\sqrt{\sum_{l=1}^{n}v^{2}_{il}\widehat{\epsilon}_{l}^{2}}\sqrt{\sum_{l=1}^{n}v^{2}_{jl}(\widehat{\epsilon}_{l}-\epsilon_{l})^{2}}+\sqrt{\sum_{l=1}^{n}v^{2}_{il}(\widehat{\epsilon}_{l}-\epsilon_{l})^{2}}\sqrt{\sum_{l=1}^{n}v^{2}_{jl}\epsilon^{2}_{l}}\right)$
$\displaystyle+|\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\epsilon_{l_{1}}\epsilon_{l_{2}}K\left(\frac{l_{1}-l_{2}}{k_{n}}\right)-\sum_{l_{1}=1}^{n}\sum_{l_{2}=1}^{n}v_{il_{1}}v_{jl_{2}}\sigma_{l_{1}l_{2}}|$
$\displaystyle=o_{p}(k_{n}\times
n^{\frac{3}{2m}-\frac{\eta}{2}})+o_{p}(k_{n}\times
n^{-1/4}\log^{-z}(n))+O_{p}(v_{n})=o_{p}(1)$

$v_{n}$ is defined in lemma 62 and $z$ is defined in assumption 6. From lemma
35, assumption 2, assumption 6, (27), and (59),
$\sup_{x\in\mathbf{R}}|Prob^{*}\left(\max_{i\in\mathcal{M}}|z_{i}|\leq
x\right)-H(x)|=o_{p}(1)$. Since

$\displaystyle\max_{i=1,...,p_{1}}\frac{\rho_{n}^{2}}{\tau_{i}}|\sum_{j=1}^{p}\frac{c_{ij}\widehat{s}_{j}}{(\lambda^{2}_{j}+\rho_{n})^{2}}|\leq\max_{i\in\mathcal{M}}\frac{\rho_{n}^{2}}{\tau_{i}}\sqrt{\sum_{j=1}^{p}\frac{c^{2}_{ij}\lambda_{j}^{2}}{(\lambda^{2}_{j}+\rho_{n})^{2}}}\sqrt{\sum_{j=1}^{p}\frac{\widehat{s}^{2}_{j}}{\lambda^{2}_{j}(\lambda^{2}_{j}+\rho_{n})^{2}}}$
(78)
$\displaystyle\leq\frac{\rho^{2}_{n}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\beta}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}{\lambda^{3}_{p}}=o_{p}(1)$

Define
$t=\frac{\rho^{2}_{n}{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\widehat{\beta}\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{2}}{\lambda^{3}_{p}}$
and assume $\widehat{\mathcal{N}}_{b_{n}}=\mathcal{N}_{b_{n}}$,

$\displaystyle
Prob^{*}\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}\leq
x\right)-H(x)\leq
Prob^{*}\left(\widehat{\mathcal{N}}^{*}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)$
(79) $\displaystyle+\left(Prob^{*}\left(\max_{i\in\mathcal{M}}|z_{i}|\leq
x+t\right)-H(x+t)\right)+Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})$
$\displaystyle\text{and
}Prob^{*}\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}\leq
x\right)-H(x)\geq-
Prob^{*}\left(\widehat{\mathcal{N}}^{*}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)$
$\displaystyle+\left(Prob^{*}\left(\max_{i\in\mathcal{M}}|z_{i}|\leq
x-t\right)-H(x-t)\right)-Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})$
$\displaystyle\Rightarrow\sup_{x\in\mathbf{R}}|Prob^{*}\left(\max_{i=1,...,p_{1}}\frac{|\widehat{\gamma}^{*}_{i}-\widehat{\gamma}_{i}|}{\widehat{\tau}^{*}_{i}}\leq
x\right)-H(x)|\leq
Prob^{*}\left(\widehat{\mathcal{N}}^{*}_{b_{n}}\neq\mathcal{N}_{b_{n}}\right)$
$\displaystyle+\sup_{x\in\mathbf{R}}|Prob^{*}\left(\max_{i\in\mathcal{M}}|z_{i}|\leq
x\right)-H(x)|+Ct(1+\sqrt{\log(p_{1})}+\sqrt{|\log(t)|})=o_{p}(1)$

For the non-degenerated joint normal distribution is absolutely continuous
with respect to Lebesgure measure, $\forall 0<\alpha<1$, $\exists
c_{1-\alpha}\in\mathbf{R}$ such that $H(c_{1-\alpha})=1-\alpha$. $\forall
0<\tau<\min(\alpha/2,(1-\alpha)/2)$, assign
$x=c_{1-\alpha+\tau},c_{1-\alpha-\tau}$ in (79), then $c_{1-\alpha-\tau}\leq
c^{*}_{1-\alpha}\leq c_{1-\alpha+\tau}$ with probability tending to $1$.
Therefore, we have
$Prob\left(\max_{i=1,...,p_{1}}\frac{1}{\widehat{\tau}_{i}}|\widehat{\gamma}_{i}-\gamma_{i}|\leq
c^{*}_{1-\alpha}\right)\leq
Prob\left(c^{*}_{1-\alpha}>c_{1-\alpha+\tau}\right)+Prob\left(\max_{i=1,...,p_{1}}\frac{1}{\widehat{\tau}_{i}}|\widehat{\gamma}_{i}-\gamma_{i}|\leq
c_{1-\alpha+\tau}\right)$; and
$Prob\left(\max_{i=1,...,p_{1}}\frac{1}{\widehat{\tau}_{i}}|\widehat{\gamma}_{i}-\gamma_{i}|\leq
c^{*}_{1-\alpha}\right)\geq-
Prob(c^{*}_{1-\alpha}<c_{1-\alpha-\tau})+Prob(\max_{i=1,...,p_{1}}\frac{1}{\widehat{\tau}_{i}}|\widehat{\gamma}_{i}-\gamma_{i}|\leq
c_{1-\alpha-\tau})$. By setting $\tau\to 0$, theorem 2 implies (22).

∎