Model selection for high-dimensional linear regression with dependent observations

TL;DR

This paper studies the effectiveness of the orthogonal greedy algorithm (OGA) combined with a high-dimensional AIC for model selection in high-dimensional linear regression with dependent data, achieving optimal prediction rates.

Contribution

It introduces a novel combination of OGA and HDAIC that attains optimal convergence rates without prior sparsity knowledge.

Findings

OGA with HDAIC achieves optimal convergence rates.

The method prevents overfitting in high-dimensional settings.

Applicable to models with dependent observations.

Abstract

We investigate the prediction capability of the orthogonal greedy algorithm (OGA) in high-dimensional regression models with dependent observations. The rates of convergence of the prediction error of OGA are obtained under a variety of sparsity conditions. To prevent OGA from overfitting, we introduce a high-dimensional Akaike's information criterion (HDAIC) to determine the number of OGA iterations. A key contribution of this work is to show that OGA, used in conjunction with HDAIC, can achieve the optimal convergence rate without knowledge of how sparse the underlying high-dimensional model is.

Equations

$$y_t=\sum _{j=1}^p{\beta }_j{x}_{tj}+{\epsilon }_t,t=1,\dots ,n,$$

$$C_1{j}^{-\gamma }\le |{\beta }_j|\le C_2{j}^{-\gamma },0<C_1\le C_2<\infty ,\gamma >1,$$

$$C_3\exp (-\beta j)\le |{\beta }_j|\le C_4\exp (-\beta j),0<C_3\le C_4<\infty ,\beta >0,$$

$$L{j}^{-\gamma }\le |{\beta }_{(j)}^{\ast }|\le U{j}^{-\gamma },$$

$$L_1\exp (-\beta j)\le |{\beta }_{(j)}^{\ast }|\le U_1\exp (-\beta j),$$

$$\sum _{j=1}^p|{\beta }_j^{\ast }{|}^{1/\gamma }<M_4,\text{mbox}forsome\gamma \ge 1,0<M_4<\infty ,$$

$$\begin{array}{r}{\hat{J}}_m={\hat{J}}_{m-1}\bigcup \{{\hat{j}}_m\},\end{array}$$

$${\hat{j}}_m=\arg \underset{1\le j\le p,j\notin {\hat{J}}_{m-1}}{\max }|{\hat{\mathbf{\mu }}}_{{\hat{J}}_{m-1},j}|,$$

$$\begin{array}{r}P(\underset{1\le j\le p}{\max }|n^{-1}\sum _{i=1}^n{z}_{ij}{\epsilon }_i|\ge c_1^{\ast }(\log p{)}^{1/2}/n^{1/2})=o(1).\end{array}$$

$$\begin{array}{r}P(\underset{1\le k,l\le p}{\max }|n^{-1}\sum _{i=1}^n{z}_{ik}{z}_{il}-{\rho }_{kl}|\ge c_2^{\ast }(\log p{)}^{1/2}/n^{1/2})=o(1),\end{array}$$

$$x_{tl}=\sum _{j=0}^{\infty }{\mathbf{w}}_j^{\top }(l){\mathbf{\delta }}_{t-j},{\epsilon }_t=\sum _{j=0}^{\infty }{\mathbf{w}}_j^{\top }(0){\mathbf{\delta }}_{t-j},$$

$$\mathrm{E}({\mathbf{\delta }}_t)=0,\underset{1\le i\le q}{\max }\mathrm{E}({\delta }_{ti}^2)<\bar{c}<\infty ,{\lambda }_{\min }(\mathrm{E}({\mathbf{\delta }}_t{\mathbf{\delta }}_t^{\top }))>{\underline{c}}_0>0,$$

$$\underset{0\le l\le p}{\max }\sum _{j=0}^{\infty }\Vert {\mathbf{w}}_j(l){\Vert }_1<M_2<\infty ,\underset{0\le l\le p}{\min }\Vert {\mathbf{w}}_0(l)\Vert >{\underline{c}}_1>0.$$

$$\text{mbox}(A1)and(A2)holdtrueunder\text{eqref}aug2--\text{eqref}aug3and\text{eqref}aug5,$$

$$p\to \infty \text{mbox}asn\to \infty ,\frac{\log p}{n}=o(1).$$

$$y_t=\sum _{j=1}^{q_0}{a}_j{y}_{t-j}+\sum _{l=1}^{p_1}\sum _{j=1}^{r_l}{\beta }_j^{(l)}{x}_{t-j+1}^{(l)}+{\epsilon }_t,$$

$$\Vert {\delta }_{tk}{\Vert }_{{\psi }_2}\le L,k=1,\dots ,q,$$

$$\sum _{j\in J}|{\beta }_j^{\ast }|\le C_{\gamma }(\sum _{j\in J}{\beta }_j^{{\ast }^2}{)}^{(\gamma -1)/(2\gamma -1)}.$$

$$\begin{array}{rl}& \sum _{j\in J}|{\beta }_j^{\ast }|\le (\sum _{j\in J}|{\beta }_j^{\ast }{|}^{1/\gamma }{)}^{\gamma /(2\gamma -1)}(\sum _{j\in J}{\beta }_j^{{\ast }^2}{)}^{(\gamma -1)/(2\gamma -1)}\\ & \le M_4^{\gamma /(2\gamma -1)}(\sum _{j\in J}{\beta }_j^{{\ast }^2}{)}^{(\gamma -1)/(2\gamma -1)},\end{array}$$

$$\sum _{j\in J}|{\beta }_j^{\ast }|\le M_1\underset{j\in J}{\max }|{\beta }_j^{\ast }|.$$

$$\underset{1\le ♯(J)\le \bar{D}(n/\log p{)}^{1/2},i\notin J}{\max }\Vert {\mathbf{\Gamma }}^{-1}(J){\mathbf{g}}_i(J){\Vert }_1<M,$$

$$\Vert {\mathbf{\beta }}^{\ast }-{\mathbf{\beta }}^{\ast }(J){\Vert }_1\le C\sum _{j\notin J}|{\beta }_j^{\ast }|,$$

$$\Vert {\mathbf{\beta }}^{\ast }-{\mathbf{\beta }}^{\ast }(J){\Vert }_1\le C{C}_{\gamma }(\sum _{j\notin J}{\beta }_j^{{\ast }^2}{)}^{(\gamma -1)/(2\gamma -1)},$$

$$\Vert {\mathbf{\beta }}^{\ast }-{\mathbf{\beta }}^{\ast }(J){\Vert }_1\le \Vert {\mathbf{\beta }}^{\ast }(J)-{\mathbf{\beta }}_J^{\ast }{\Vert }_1+\sum _{j\notin J}|{\beta }_j^{\ast }|,$$

$$\begin{array}{rl}& \mathbf{\Gamma }(J{)}^{-1}\mathbf{\Gamma }(J)({\mathbf{\beta }}^{\ast }(J)-{\mathbf{\beta }}_J^{\ast })=\mathbf{\Gamma }(J{)}^{-1}(\sum _{l\notin J}{\rho }_{j_{1}l}{\beta }_l^{\ast },\dots ,\sum _{l\notin J}{\rho }_{j_{q}l}{\beta }_l^{\ast }{)}^{\top }\\ & =\sum _{l\notin J}{\beta }_l^{\ast }\mathbf{\Gamma }(J{)}^{-1}{\mathbf{g}}_l(J).\end{array}$$

$$k_0=♯(N_n)<M_7,$$

$$\begin{array}{rl}& {\mathrm{E}}_n\{y(\mathbf{x})-{\hat{y}}_m(\mathbf{x}){\}}^2={\mathrm{E}}_n\{y(\mathbf{x})-y_{{\hat{J}}_m}(\mathbf{x}){\}}^2+{\mathrm{E}}_n\{y_{{\hat{J}}_m}(\mathbf{x})-{\hat{y}}_m(\mathbf{x}){\}}^2,\end{array}$$

$$\displaystyle\lambda_{\min}(\mbox{\boldmath$\Gamma$})\geq\lambda_{1}>0,$$

$$\log p=o(n),K_n=\bar{\delta }{\left(\frac{n}{\log p}\right)}^{1/2}$$

$$\bar{\tau }=\sup \mathbf{\tau }\equiv \sup \{\tau :\tau >0,\underset{n\to \infty }{\mathrm{limsup}}\frac{\tau c_2^{\ast }}{\underset{♯(J)\le \tau (n/\log p{)}^{1/2}}{\min }{\lambda }_{\min }(\mathbf{\Gamma }(J))}\le 1\}.$$

Full text

63M30,

62F07,

62F12,

Best $m$-term approximations,

high-dimensional Akaike’s information criterion,

orthogonal greedy algorithm,

sparsity conditions,

time series,

arXiv:0000.0000 \startlocaldefs

\endlocaldefs

t1Research was supported in part by the Science Vanguard Research Program of the Ministry of Science and Technology, Taiwan

1 Introduction

Model selection for high-dimensional regression models has been one of the most vibrant topics in statistics over the past decade. It also has broad applications in a variety of important fields such as bioinformatics, quantitative finance, image processing, and advanced manufacturing; see Negahban et al. (2012) and Ing et al. (2017) for further discussion. A typical high-dimensional regression model takes the following form:

[TABLE]

where $n$ is the sample size, $x_{t1},\cdots,x_{tp}$ are predictor variables, $\varepsilon_{t}$ are mean-zero random disturbance terms, and $p=p_{n}$ is allowed to be much larger than $n$. There are computational and statistical difficulties in estimating the regression function by standard regression methods owing to $p\gg n$. However, by assuming sparsity conditions on $\beta_{j}$, eigenvalue conditions on the covariance (correlation) matrix of the predictor variables, and distributional conditions on $\varepsilon_{t}$ or $x_{tj}$, it has been shown that consistent estimation of the regression function or optimal prediction is still possible either through penalized least squares methods (see Zhao and Yu (2006), Candes and Tao (2007), Bickel et al. (2009), and Zhang (2010)) or through greedy forward selection algorithms (see Bühlmann (2006), Chen and Chen (2008), Wang (2009), Fan and Lv (2008), and Ing and Lai (2011)).

The vast majority of studies on model (1.1), however, have focused on situations where $\mathbf{x}_{t}=(x_{t1},\ldots,x_{tp})^{\top}$ are nonrandom and $\varepsilon_{t}$ are independently and identically distributed (i.i.d.) or $(\mathbf{x}_{t},\varepsilon_{t})$ are i.i.d., which regrettably preclude most serially correlated data. In fact, (1.1) can encompass a broad array of time series models if these restrictions are relaxed. For example, it becomes the well-known autoregressive (AR) model when $x_{tj}=y_{t-j}$. Since the predictor variables in AR models have a natural ordering, a commonly used sparsity condition is

[TABLE]

in which $|\beta_{j}|$ decay polynomially, or

[TABLE]

in which $|\beta_{j}|$ decay exponentially (see Shibata (1980) and Ing (2007)). Moreover, the model selection problem in the AR case is simplified to an order selection one, which has been well explored in the literature (see Shibata (1980)). When $x_{tj},j=1,\ldots,p$, do not have a natural ordering, e.g., the autoregressive exogenous (ARX) model, (1.2) and (1.3) can be generalized as

[TABLE]

and

[TABLE]

respectively, where $0<L\leq U<\infty$, $0<L_{1}\leq U_{1}<\infty$, and $|\beta^{*}_{(1)}|\geq|\beta^{*}_{(2)}|\geq\cdots\geq|\beta^{*}_{(p)}|$ is a rearrangement of $\{|\beta^{*}_{j}|\}$ in decreasing order with $\beta^{*}_{j}=\sigma_{j}\beta_{j}$ and $\sigma_{j}^{2}=\mathrm{E}(x_{tj}^{2})$. However, unlike the order selection problem, the model selection problem in (1.1) with dependent observations and with coefficients satisfying (1.4) or (1.5) seems to be seldom investigated. The problem becomes more challenging when $\beta_{j}$ may obey either one of (1.4), (1.5), or $k_{0}\ll n$, but it is unclear which of the three is true. Here, $k_{0}$ denotes the number of nonzero coefficients in model (1.1), and $k_{0}\ll n$ is referred to as the strong sparsity condition.

In this paper, we assume that the $(\mathbf{x}_{t},\varepsilon_{t})$ in model (1.1) is a time series obeying concentration inequalities (2.2) and (2.3). We also assume that the $\beta_{j}$ in model (1.1) follow one of the following sparsity conditions: (i) (A3), (ii) (A4), or (iii) $k_{0}\ll n$, where (A3) and (A4) are defined in Section 2.1. Note that (A3) includes (1.4) and

[TABLE]

as special cases, whereas (A4) contains (1.5). We use the orthogonal greedy algorithm (OGA) (Temlyakov, 2000) to sequentially include candidate variables and introduce a high-dimensional Akaike’s information criterion (HDAIC) to determine the number of OGA iterations. This model selection procedure is denoted by OGA+HDAIC. A key contribution of this paper is to show that OGA+HDAIC achieves the optimal convergence rate without knowing which sparsity condition among (i), (ii), and (iii) would follow, thereby alleviating the dilemma mentioned in the previous paragraph.

Following this introduction, the rest of the paper is organized as follows. In Section 2.1, we introduce OGA and the assumptions required for our asymptotic analysis of the algorithm. Section 2.2 derives an error bound for OGA, which is the sum of an approximation error and a term accounting for the sampling variability. Since the approximation error decreases as the number $m$ of iterations increases and the sampling variability increases with $m$, the optimal $m$ can be determined by equating the two terms in the error bound for OGA. This approach, however, is infeasible because not only does the solution involve the unknown parameters in (A3) or (A4), but it is unknown which kind of sparsity among (i), (ii), and (iii) holds true. To overcome this difficulty, Theorem 3.1 in Section 3.1 proposes using HDIC to determine the number of iterations, and shows that OGA+HDAIC is rate optimal regardless of which sparsity condition is true. In Section 3.2, we offer a comprehensive comparison of our results with those in Negahban et al. (2012) and Ing and Lai (2011), in which the statistical properties of Lasso (Tibshirani, 1996) and OGA, respectively, are explored under model (1.1) with independent observations. In this connection, Section 3.2 also discusses the papers by Basu and Michailidis (2015) and Wu and Wu (2016), which investigate the performance of Lasso under sparse high-dimensional time series models. The proof of Theorem 3.1 is given in Section 3.3. We conclude in Section 4. An appendix consisting of some technical results is given at the end of the paper. A simulation study to illustrate the performance of OGA+HDAIC, along with further technical details, is deferred to the supplementary material.

2 Asymptotic Theory of OGA in Weakly Sparse Models

This section aims at establishing the convergence rate of OGA under sparse high-dimensional regression models with dependent observations. The definition of OGA and the assumptions required for our analysis of OGA are given in Section 2.1. The main result of this section is stated and proved in Section 2.2.

2.1 Models and Assumptions

We assume that $\{(\mathbf{x}_{t},\varepsilon_{t})\}$ in model (1.1) is a zero-mean stationary time series satisfying $\mathrm{E}(\mathbf{x}_{t}\varepsilon_{t})=\mathbf{0}$. The OGA is a recursive procedure that selects variables from the set of predictor variables in (1.1) one at a time. Define $\mathbf{X}_{i}=(x_{1i},\ldots,x_{ni})^{\top}$, $\mathbf{Z}_{i}=(z_{1i},\ldots,z_{ni})^{\top}=\mathbf{X}_{i}/\sigma_{i}$, and $\mathbf{Y}=(y_{1},\ldots,y_{n})^{\top}$. The algorithm is initialized by setting $\hat{J}_{0}=\emptyset$, where $\hat{J}_{m}$ denotes the index set of the variables chosen by OGA at the $m$-th iteration. For $m\geq 1$, $\hat{J}_{m}$ is recursively updated by

[TABLE]

where

[TABLE]

with $\hat{\bm{\mu}}_{J,i}=\mathbf{Z}^{\top}_{i}(\mathbf{I}-\mathbf{H}_{J})\mathbf{Y}/(n^{1/2}\|\mathbf{Z}_{i}\|)$, $\|\mathbf{a}\|$ denoting the $L_{2}$-norm of vector $\mathbf{a}$, and $\mathbf{H}_{J},J\subseteq{\cal P}\equiv\{1,\ldots,p\}$, being the orthogonal projection matrix onto the linear span of $\{\mathbf{Z}_{i},i\in J\}$ ($\mathbf{H}_{\emptyset}=\mathbf{0}$).

To investigate the performance of OGA, we make the following distributional assumptions:

(A1)

There exists $c_{1}^{*}>0$ such that

[TABLE]

(A2)

There exists $c_{2}^{*}>0$ such that

[TABLE]

where $\rho_{kl}=\mathrm{E}(z_{1k}z_{1l})$.

The following examples help illustrate (A1) and (A2). Let $\lambda_{\min}(A)$ ($\lambda_{\max}(A)$) denote the minimum (maximum) eigenvalue of matrix $A$ and $\|\mathbf{a}\|_{1}$ the $L_{1}$-norm of vector $\mathbf{a}$.

Example 1 (Gaussian linear processes.) Let

[TABLE]

where $\bm{\delta}_{t}=(\delta_{t1},\ldots\delta_{tq})^{\top}$ are i.i.d. Gaussian random vectors satisfying

[TABLE]

and $\bm{w}_{j}(l)$ obey

[TABLE]

Then, by making use of the Hanson-Wright inequality (see Theorem 1.1 of Rudelson and Vershynin (2013)), it is shown in Section S1 of the supplementary material that

[TABLE]

where (2.8) is given by

[TABLE]

As an application, we consider a high-dimensional ARX model,

[TABLE]

in which $p=q_{0}+\sum_{l=1}^{p_{1}}r_{l}$ satisfies (2.8), $1-\sum_{j=1}^{q_{0}}a_{j}z^{j}\neq 0$ for all $|z|\leq 1+\iota$ with $\iota$ being some positive constant, $\sum_{j=1}^{q_{0}}|a_{j}|+\sum_{l=1}^{p_{1}}\sum_{j=1}^{r_{l}}|\beta_{j}^{(l)}|<M_{5}<\infty$, $x_{t}^{(l)}=\epsilon_{t}^{(l)}+\sum_{j=1}^{\infty}b_{j}^{(l)}\epsilon_{t-j}^{(l)}$, with $\sum_{j=1}^{\infty}|b_{j}^{(l)}|<M_{6}<\infty$ for all $1\leq l\leq p$, and $\bm{\delta}_{t}=(\epsilon_{t}^{(1)},\ldots,\epsilon_{t}^{(p)},\varepsilon_{t})^{\top}$ are i.i.d. $(p+1)$-dimensional Gaussian random vectors obeying (2.5) with $q=p+1$. It is not difficult to see that (2.4) and (2.6) are fulfilled by the regressor variables and the error term in (2.9). Hence (A1) and (A2) are applicable to model (2.9).

Example 2 (Linear processes with sub-Gaussian innovations.) Suppose that (2.4)–(2.6) and (2.8) are satisfied except that the Gaussianity of $\bm{\delta}_{t}$ is replaced by

[TABLE]

where $\lVert\cdot\rVert_{\psi_{2}}$ denotes the $\psi_{2}$ Orlicz norm and $L$ is some positive number. We note that (2.10) is fulfilled by sub-Gaussian random variables. Assume $q=p^{s}$ for some $0\leq s<\infty$. Then by making use of the concentration inequality given in Theorem 1.4 of Adamczak and Wolff (2015), it can be shown that (A1) and (A2) hold for some large $c_{1}^{*}$ and $c_{2}^{*}$. For more details, see Huang and Ing (2019). In addition, the regressor variables and the error term in (2.9) still obey (A1) and (A2), provided assumption (2.10) is used in place of the Gaussian assumption in Example 1.

We also need a sparsity condition on regression coefficients:

(A3)

There is $0<\bar{M}_{0}<\infty$ such that $\sum_{j=1}^{p}\beta_{j}^{*^{2}}\leq\bar{M}_{0}$. In addition, there exist $\gamma\geq 1$ and $0<C_{\gamma}<\infty$ such that for any $J\subseteq{\cal P}$,

[TABLE]

When $\gamma=1$, (2.11) and (1.6) are equivalent. However, (2.11) is weaker than (1.6) for $\gamma>1$. To see this, note that if (1.6) is true for some $\gamma>1$, then by Hölder’s inequality,

[TABLE]

implying that (2.11) holds for $C_{\gamma}=M_{4}^{\gamma/(2\gamma-1)}$. In view of the connection between (2.11) and (1.6), the parameter $\gamma$ in (2.11) can be understood as an index to describe the degree of sparseness in the underlying high-dimensional models. The larger the $\gamma$ is, the sparser the model is. Although assumptions similar to (1.6) are quite popular for high-dimensional regression analysis (see, e.g., Wang et al. (2014)), there is a subtle difference between (2.11) and (1.6). To see this, assume that (1.4) holds for some $\gamma>1$. Then, (2.11) holds for the same $\gamma$ (see Lemma A1.2 in the Appendix), whereas (1.6) is violated due to $L(1+\log p)\leq\sum_{j=1}^{p}|\beta^{*}_{j}|^{1/\gamma}\leq U(1+\log p)$. It is worth mentioning that (1.4) not only plays an important role in time series modeling, it also allows us to demonstrate that the approximation error of the population counterpart of OGA (which is defined at the beginning of Section A1 and is referred to as the population OGA) is almost as small as that of the best $m$-term approximation (see (3.24) and Lemma A1.3 in the Appendix). In the sequel, we refer to (2.11) as the ‘polynomial decay’ case, owing to its connection with (1.4). To broaden OGA’s applications, we also consider a coefficient condition sparser than (2.11):

(A4)

There exists $0<M_{0}<\infty$ such that $\max_{1\leq j\leq p}|\beta_{j}^{*}|\leq M_{0}$. Moreover, there exists $M_{1}>1$ such that for any $J\subseteq{\cal P}$,

[TABLE]

Assumption (A4) is referred to as the ‘exponential decay’ case because (1.5) is included by (2.12).

The following assumption on the covariance structure of $\mathbf{z}_{t}=(z_{1},\ldots,z_{p})^{\top}$ is frequently used throughout the paper. Define $\bm{\Gamma}(J)=\mathrm{E}\{\mathbf{z}_{t}(J)\mathbf{z}_{t}^{\top}(J)\}$ and $\bm{g}_{i}(J)=\mathrm{E}(z_{ti}\mathbf{z}_{t}(J))$, where $J\subseteq{\cal P}$ and $\mathbf{z}_{t}(J)=(z_{ti},i\in J)^{\top}$.

(A5)

For some positive numbers $\bar{D}$ and $M$,

[TABLE]

where $\sharp(J)$ denotes the cardinality of $J$ .

Since $\bm{\Gamma}^{-1}(J)\bm{g}_{i}(J)=\arg\min_{\bm{c}\in R^{\sharp(J)}}\mathrm{E}(z_{ti}-\bm{c}^{\top}\mathbf{z}_{t}(J))^{2}$, (2.13) essentially says that the regression coefficients for $z_{ti}$ on $\mathbf{z}_{t}(J)$ with all $i\notin J$ and $\sharp(J)\leq\bar{D}(n/\log p)$ are $L_{1}$ bounded. This condition holds even when $z_{t1},\cdots,z_{tp}$ are highly correlated; see Section S3 of the supplementary document. Let $\bm{g}_{y}(J)=\mathrm{E}(y_{t}\mathbf{z}_{t}(J))$ and \bm{\beta}^{*}(J)=\bm{\Gamma}^{-1}(J)\bm{g}_{y}(J)$$=\arg\min_{\bm{c}\in R^{\sharp(J)}}\mathrm{E}(y_{t}-\bm{c}^{\top}\mathbf{z}_{t}(J))^{2}, which is the regression coefficients for $y_{t}$ on $\mathbf{z}_{t}(J)$. By making use of (2.13), we will show later that for any $J\subseteq{\cal P}$ with $\sharp(J)\leq\bar{D}(n/\log p)^{1/2}$, there exists $0<C<\infty$ such that

[TABLE]

where $\bm{\beta}^{*}=(\beta_{1}^{*},\ldots,\beta^{*}_{p})^{\top}$ and $\bm{\beta}^{*}(J)$ here is regarded as a $p$-dimensional vector with undefined entries set to 0. Inequality (2.14) is referred to as the uniform Baxter’s inequality. (For more details on Baxter’s inequality in autoregressive modeling, see Baxter (1962), Berk (1974), and Pourahmadi (1989).) This inequality can be used together with (2.11) to yield, for all $\sharp(J)\leq\bar{D}(n/\log p)^{1/2}$,

[TABLE]

which is one of the key ingredients in our asymptotic analysis of OGA+HDAIC.

To derive (2.14) from (2.13), we may assume without loss of generality that $J=\{j_{1},\ldots,j_{q}\}$ for some $1\leq q\leq\bar{D}(n/\log p)^{1/2}$, where $j_{i},i=1,\ldots q$, are distinct elements in ${\cal P}$. Note first that

[TABLE]

where $\bm{\beta}^{*}_{J}=(\beta^{*}_{j_{1}},\ldots,\beta^{*}_{j_{q}})^{\top}$. Denote $\bm{g}_{y}(J)$ by $(\gamma_{y,j_{1}},\ldots,\gamma_{y,j_{q}})^{\top}$. Then, it follows that $\gamma_{y,j_{i}}=\sum_{l=1}^{p}\rho_{j_{i}l}\beta_{l}^{*}$, $i=1,\ldots,q$, and hence

[TABLE]

Taking the $L_{1}$-norm on both sides, (2.14) (with $C=M+1$) follows from (2.16) and (2.13).

Before closing this section, we remark that (2.12) can be viewed as a limiting case of (2.11). To see this, note that (2.12) implies that for any $J\subseteq{\cal P}$, $\sum_{j\in J}|\beta_{j}^{*}|\leq M_{1}(\sum_{j\in J}\beta_{j}^{*^{2}})^{\lim_{\gamma\to\infty}(\gamma-1)/(2\gamma-1)}$. In addition, the strong sparsity condition,

[TABLE]

where $N_{n}=\{j:\beta_{j}^{*}\neq 0,1\leq j\leq p\}$ and $M_{7}$ is some positive integer, is also a limiting case of (2.11) because (2.17) yields that for any $J\subseteq{\cal P}$, $\sum_{j\in J}|\beta_{j}^{*}|\leq M_{7}^{1/2}(\sum_{j\in J}\beta_{j}^{*^{2}})^{\lim_{\gamma\to\infty}(\gamma-1)/(2\gamma-1)}$.

2.2 Rates of Convergence of the OGA

Let $\mathbf{x}=(x_{1},\ldots x_{p})^{\top}$ be independent of and have the same covariance structure as $\{\mathbf{x}_{t}\}$ and $y(\mathbf{x})=\sum_{j=1}^{p}\beta_{j}x_{j}$. Then, $y(\mathbf{x})$ can be predicted by $\hat{y}_{m}(\mathbf{x})=\mathbf{x}^{\top}(\hat{J}_{m})\hat{\bm{\beta}}(\hat{J}_{m})$, where $\mathbf{x}(J)=(x_{i},i\in J)$ and $\hat{\bm{\beta}}(J)=(\sum_{t=1}^{n}\mathbf{x}_{t}(J)\mathbf{x}_{t}^{\top}(J))^{-1}\sum_{t=1}^{n}\mathbf{x}_{t}(J)y_{t}$. Note also that $\hat{y}_{m}(\mathbf{x})=\mathbf{z}^{\top}(\hat{J}_{m})\hat{\bm{\beta}}^{*}(\hat{J}_{m})$, where $\mathbf{z}(J)=(z_{i},i\in J)$ with $z_{i}=x_{i}/\sigma_{i}$, and $\hat{\bm{\beta}}^{*}(\hat{J}_{m})=(\sum_{t=1}^{n}\mathbf{z}_{t}(J)\mathbf{z}_{t}^{\top}(J))^{-1}\sum_{t=1}^{n}\mathbf{z}_{t}(J)y_{t}$. One of the most natural performance measures for $\hat{y}_{m}(\mathbf{x})$ is the conditional mean squared prediction error (CMSPE),

[TABLE]

where $\mathrm{E}_{n}(\cdot)=\mathrm{E}(\cdot|y_{1},\mathbf{x}_{1},\cdots,y_{n},\mathbf{x}_{n})$, and y_{J}(\mathbf{x})$$=\bm{\beta}^{*^{\top}}(J)\mathbf{z}(J). A convergence rate of the left-hand side of (LABEL:aug14) is established in the next theorem.

Theorem 2.1**.**

Suppose that (1.1), (A1)–(A3), (A5),

[TABLE]

and

[TABLE]

hold, where \mbox{\boldmath\Gamma}=\mathrm{E}(\mathbf{z}\mathbf{z}^{\top}) and $0<\bar{\delta}<\min\{\bar{\tau},\bar{D}\}$, with $\bar{D}$ defined in assumption (A5) and

[TABLE]

Then,

[TABLE]

Moreover, if (A4) holds instead of (A3), then

[TABLE]

where $G_{3}$ is some positive constant given in (A1.2) in the Appendix.

Proof. We first prove (2.22). Recall $\hat{J}_{k}=\{\hat{j}_{1},\cdots,\hat{j}_{k}\}$ and define

[TABLE]

where $(\hat{y}_{1;J},\ldots,\hat{y}_{n;J})^{\top}=\mathbf{H}_{J}\mathbf{Y}$. Moreover, let

[TABLE]

and

[TABLE]

where $\mu_{J,i}=\mathrm{E}[(y(\mathbf{x})-y_{J}(\mathbf{x}))z_{i}]$, $s>0$ is some large constant, and $\tilde{\xi}=2/(1-\xi)$ with $0<\xi<1$ being arbitrarily given.

By an argument similar to that of (3.10) in Ing and Lai (2011), it follows that for all $1\leq q\leq m$,

[TABLE]

This and (A1.1) in the Appendix, which gives an error bound for the population OGA under (A3), lead to

[TABLE]

Moreover, (A3) and (2.19) imply that for any $0\leq i\leq m-1$,

[TABLE]

and hence

[TABLE]

By (2.25),

[TABLE]

Since $A_{n}(m)$ decreases as $m$ increases, (2.24) and (2.26) yield that for all $1\leq m\leq K_{n}$ and some $C_{2}>0$,

[TABLE]

We show in Section S1 of the supplementary material that

[TABLE]

where $\hat{\bm{\Gamma}}(J)=n^{-1}\sum_{t=1}^{n}\mathbf{z}_{t}(J)\mathbf{z}^{\top}_{t}(J)$ and

[TABLE]

noting that the positiveness of the denominator is ensured by (2.20) and (2.21). With the help of (2.28), (A1), (A2), and (A5), it is shown in the same section that there exists a sufficiently large $s$ such that

[TABLE]

which, together with (2.27), yields

[TABLE]

Moreover, we have

[TABLE]

which is also proved in Section S1 of the supplementary material. In view of (2.31), (2.32), and the fact that $(\log p/n)^{1-(2\gamma)^{-1}}\leq m^{-2\gamma+1}$ if $m\leq(n/\log p)^{(2\gamma)^{-1}}$ and $(\log p/n)^{1-(2\gamma)^{-1}}\leq n^{-1}m\log p$ if $m\geq(n/\log p)^{(2\gamma)^{-1}}$, the desired conclusion (2.22) follows. Equation (2.23) follows from (A1.2) in the Appendix (which gives an error bound for the population OGA under (A4)) and an argument similar to that used to prove (2.22). We skip the details in order to save space. $\Box$

Remark 1. It is easy to see that the $\bm{\tau}$ defined in (2.21) is nonempty. In particular, $x\in\bm{\tau}$ for any $x\in(0,\lambda_{1}/c_{2}^{*}]$. It is also not difficult to see that $\bar{\tau}<a/c_{2}^{*}$ for any $a>1$.

In view of (2.22), to strike a suitable balance between squared bias and variance, one should choose $m\approx(n/\log p)^{1/2\gamma}$ in the polynomial decay case, which yields a rate of convergence, $(n^{-1}\log p)^{1-(2\gamma)^{-1}}$. Similarly, (2.23) suggests that the best convergence rate one can expect in the exponential decay case is $n^{-1}\log n\log p$, which is ensured by selecting $m\approx\log n/G_{3}$. The optimality of the rates, $(n^{-1}\log p)^{1-(2\gamma)^{-1}}$ and $n^{-1}\log n\log p$, will be discussed further in Section 3.2. In most practical situations, however, not only do we not know what $\gamma$ or $G_{3}$ is, we do not even know which of (A3) and (A4) is true. To attain the aforementioned optimal convergence rates without knowing the degree of sparseness, a data-driven method to determine the number of OGA iterations is called for. In the next section, we show that HDAIC (see (3.1)) can fulfill this need.

Finally, we note that if (2.2) and (2.3) are weakened to

[TABLE]

and

[TABLE]

respectively, where $0\leq\bar{c}_{1},\bar{c}_{2}<\infty$ are some constants, and (2.20) is strengthened to

[TABLE]

where $\bar{c}=\max\{\bar{c}_{1},\bar{c}_{2}\}$ and $\underline{\delta}$ is some positive constant, then (2.22) and (2.23) become

[TABLE]

and

[TABLE]

respectively. While (2.33) and (2.34) are satisfied by a broader class of time series models (see Wu and Wu (2016) for a detailed discussion), to determine the optimal $m$ in (2.36) or (2.37), the HDAIC must also be corrected according to the value of $\bar{c}$. This kind of correction, however, is hardly implemented in practice because $\bar{c}$ is in general unknown.

3 Analysis of OGA+HDAIC

In Section 3.1, the rate of convergence of OGA+HDAIC is established under various sparsity conditions; see Theorem 3.1. Comparisons of Theorem 3.1 and related existing results are given in Section 3.2. The proof of Theorem 3.1 is provided in Section 3.3.

3.1 Error bounds for OGA+HDAIC

Define

[TABLE]

where $\hat{\sigma}_{J}^{2}=n^{-1}\mathbf{Y}^{\top}(\mathbf{I}-\mathbf{H}_{J})\mathbf{Y}$ and $s_{a}$ is some positive constant, and define

[TABLE]

noting that $\hat{J}_{k}$ is defined in (2.1).

Theorem 3.1**.**

Suppose that (1.1), (A1), (A2), (A5), (2.19), (2.20), and

[TABLE]

hold. Then, for

[TABLE]

where $\bar{B}$ is defined in (2.29), we have

(i)

[TABLE]

provided (A3) is true,

(ii)

[TABLE]

provided (A4) is true and $\log p=o(n/(\log n)^{2})$,

(iii)

[TABLE]

provided $\mathrm{E}(y_{t}^{2})$ is bounded above by a finite constant and

[TABLE]

Remark 2. The sparsity condition (3.7) implies k_{0}=o\big{(}(n/\log p)^{1/3}\big{)}, allowing $k_{0}$ to grow to $\infty$ slowly with $n$. Moreover, (3.6) also holds when (2.19) is weakened to

[TABLE]

for some $\eta>0$; see Section S2 in the supplementary document. However, since it is unknown which kind of sparsity condition is true among those described in (i), (ii), and (iii) of Theorem 3.1, and since (2.19) appears to be indispensable for the proofs of (3.4) and (3.5), the latter assumption is still adopted in our unified theory.

Remark 3. We briefly discuss extensions of Theorems 2.1 and 3.1 to the following multivariate time series models,

[TABLE]

where $\bm{y}_{t}$, $\bm{e}_{t}$, and $\bm{b}_{j}$ are $d$-dimensional vectors, $d$ is allowed to grow to infinity with $n$, and $\{(\bm{e}^{\top}_{t},\mathbf{x}^{\top}_{t})^{\top}\}$ is a zero-mean stationary time series satisfying $\mathrm{E}(\mathbf{x}_{t}\bm{e}^{\top}_{t})=\mathbf{0}$. Define $\hat{\bm{\psi}}_{J,i}=\|{\cal Y}^{\top}(\mathbf{I}-\mathbf{H}_{J})\mathbf{Z}^{\top}_{i}\|/(n^{1/2}\|\mathbf{Z}_{i}\|)$, where ${\cal Y}=(\bm{y}_{1},\cdots,\bm{y}_{n})^{\top}$, and $\mathbf{H}_{J}$ and $\mathbf{Z}_{i}$ are defined as in Section 2.1. A multivariate version of OGA, MOGA, is initialized by $\hat{L}_{0}=\emptyset$. For $m\geq 1$, $\hat{L}_{m}$ is recursively updated by

[TABLE]

where $\hat{l}_{m}=\arg\max_{1\leq l\leq p,l\notin\hat{L}_{m-1}}\hat{\bm{\psi}}_{\hat{L}_{m-1},l}$. Consider multivariate extensions of the sparsity conditions (A3) and (A4):

(A3${}^{{}^{\prime}}$)

There is $0<\bar{M}_{0}<\infty$ such that

[TABLE]

Moreover, there exist $\gamma\geq 1$ and $0<C_{\gamma}<\infty$ such that for any $J\subseteq{\cal P}$,

[TABLE]

where $\bm{b}_{j}^{*}=\sigma_{j}\bm{b}_{j}$.

(A4${}^{{}^{\prime}}$)

There is $0<M_{0}<\infty$ such that

[TABLE]

Moreover, there exists $M_{1}>1$ such that for any $J\subseteq{\cal P}$,

[TABLE]

Moreover, a natural generalization of (A1) under model (3.9) is

(A1${}^{{}^{\prime}}$)

There exists $c_{1}^{*}>0$ such that

[TABLE]

where $(\varepsilon_{t1},\ldots,\varepsilon_{td})^{\top}=\bm{e}_{t}$.

Let $\mathbf{x}=(x_{1},\ldots,x_{p})^{\top}$ be defined as in Section 2.2 and $\bm{y}(\mathbf{x})=\sum_{j=1}^{p}\bm{b}_{j}x_{j}$. Then, $\bm{y}(\mathbf{x})$ can be predicted by $\hat{\bm{y}}_{m}(\mathbf{x})=\hat{\bm{B}}(\hat{L}_{m})^{\top}\mathbf{x}(\hat{L}_{m})$, where $\hat{\bm{B}}(J)=(\sum_{t=1}^{n}\mathbf{x}_{t}(J)\mathbf{x}_{t}^{\top}(J))^{-1}\sum_{t=1}^{n}\mathbf{x}_{t}(J)\bm{y}_{t}^{\top}$. Suppose that

[TABLE]

for some $\zeta>0$. Then, under (3.10) and the assumptions of Theorem 2.1, with (A1), (A3), and (A4) replaced by (A1${}^{{}^{\prime}}$), (A3${}^{{}^{\prime}}$), and (A4${}^{{}^{\prime}}$), it can be shown that

[TABLE]

and for some $G_{4}>0$,

[TABLE]

To choose a suitable number of MOGA iterations, one may consider a multivariate extension of HDAIC (MHDAIC),

[TABLE]

where \hat{\bm{\Sigma}}_{J}=(nd)^{-1}tr\big{(}{\cal Y}^{\top}(\mathbf{I}-\mathbf{H}_{J}){\cal Y}\big{)} and $\iota_{a}$ is some positive constant, and define

[TABLE]

We conjecture that $d^{-1}\mathrm{E}_{n}\|\bm{y}(\mathbf{x})-\hat{\bm{y}}_{\hat{m}_{n}}(\mathbf{x})\|^{2}$ is of order $O_{p}((\log pd/n)^{1-1/(2\gamma)})$, $O_{p}(\log n\log pd/n)$, or $O_{p}(k_{0}\log pd/n)$ under (A3${}^{{}^{\prime}}$), (A4${}^{{}^{\prime}}$), or a strong sparsity condition resembling (3.7), respectively. However, the rigorous proof of this result and those of (3.11) and (3.12) are out of the scope of this paper, and are left for future work.

3.2 Some comparisons with existing results

It would be interesting to compare (3.4) with Corollary 3 of Negahban et al. (2012), which provides an error bound for Lasso in the following high-dimensional regression model,

[TABLE]

where $\{\epsilon_{t}\}$ is a sequence of i.i.d. $N(0,\sigma^{2})$ random variables and $\{x_{tj}\}$ are non-random constants satisfying $n^{-1}\sum_{t=1}^{n}x_{tj}^{2}\leq 1,1\leq j\leq p$, and the restricted eigenvalue condition defined in (31) of their paper. When

[TABLE]

for some $\gamma\geq 1$, it is shown in the corollary that

[TABLE]

where $\hat{\bm{\beta}}_{\lambda_{n}}$ is the Lasso estimate of $\bm{\beta}^{*}$ with $\lambda_{n}=4\sigma(\log p/n)^{1/2}$. On the other hand, (3.4) implies that under model (1.1),

[TABLE]

In addition to allowing for serially correlated data, (3.16) may lead to a faster convergence rate than (3.15). In particular, the bound on the right-hand side of (3.15) is larger than that on the right-hand side of (3.16) by a factor of $\log p$ as $p\to\infty$ when (1.4), with $\gamma>1$, and (3.14) follows.

Assuming that the $\{\mathbf{x}_{t}\}$ and $\{\epsilon_{t}\}$ in (3.13) are generated according to independent, centered, Gaussian stationary time series, Proposition 3.3 of Basu and Michailidis (2015) establishes for Lasso the following bounds:

[TABLE]

and

[TABLE]

where $p\to\infty$, $k_{0}=O(n/\log p)$, and $\lambda_{n}\geq c^{*}(\log p/n)^{1/2}$ for some $c^{*}>0$. By (3.6) and an argument used in Section S2 of the supplementary document, it can be shown that under model (1.1),

[TABLE]

and

[TABLE]

Although (3.17)–(3.20) suggest that Lasso and OGA+HDAIC share the same error rate in the case of $k_{0}\ll n$, they are obtained under somewhat different assumptions. Note first that unlike (3.17) and (3.18), (3.19) and (3.20) do not require that $\{\mathbf{x}_{t}\}$ and $\{\varepsilon_{t}\}$ are independent, and hence are applicable to ARX models. Moreover, (3.17) and (3.18) are established under

[TABLE]

and

[TABLE]

where $0<\underline{s}\leq\bar{S}<\infty$ and $f_{\mathbf{x}}(\theta)=[1/(2\pi)]\sum_{l=-\infty}^{\infty}\bm{\Gamma}_{\mathbf{x}}(l)\exp(-il\theta)$ with $\bm{\Gamma}_{\mathbf{x}}(l)=\mathrm{E}(\mathbf{x}_{t}\mathbf{x}_{t+l}^{\top})$. Assumption (3.22) is comparable to (2.19) (which assumes that $\lambda_{\min}(\bm{\Gamma})$ is bounded away from zero and is needed for proving (3.19) and (3.20)), but is more stringent than the latter because

[TABLE]

Maximum eigenvalue assumptions like (3.21) are not required for (3.19) and (3.20). This type of assumption can be easily violated when the components of $\mathbf{x}_{t}$ are highly correlated, as illustrated by an ARX example in Section S3 of the supplementary document, in which $\lambda_{\max}(\bm{\Gamma})\to\infty$ as $p\to\infty$ and hence $\mathrm{ess\,sup}_{\theta\in[-\pi,\pi]}\lambda_{\max}(f_{\mathbf{x}}(\theta))\geq[1/(2\pi)]\lambda_{\max}(\bm{\Gamma})$. On the other hand, while (3.19) and (3.20) are obtained under the beta-min condition given in (3.7), (3.17) and (3.18) do not assume any beta-min condition. Wu and Wu (2016) also investigate the performance of Lasso under (3.13) with $k_{0}\ll n$ and $\{x_{ti}\}$ being nonrandom and obeying the restricted eigenvalue condition defined in (4.2) of their paper. They allow $\{\epsilon_{t}\}$ to be a stationary process following some general moment and dependence conditions. The error rates that they derive for Lasso, however, are usually larger than those in (3.17)–(3.20).

In fact, it can be argued that all error bounds obtained in Theorem 3.1 are rate optimal. To see this, let $\hat{J}(m),1\leq m\leq K_{n}$, be a sequence of nested models chosen from $p$ candidate variables in a data-driven fashion, where $\sharp(\hat{J}(m))=m$. The CMSPE of model $\hat{J}(m)$ is $\mathrm{E}_{n}(y(\mathbf{x})-\hat{y}_{\hat{J}(m)}(\mathbf{x}))^{2}=\mathrm{E}_{n}(y(\mathbf{x})-y_{\hat{J}(m)}(\mathbf{x}))^{2}+\mathrm{E}_{n}(\hat{y}_{\hat{J}(m)}(\mathbf{x})-y_{\hat{J}(m)}(\mathbf{x}))^{2}$, where $\hat{y}_{J}(\mathbf{x})=\mathbf{x}^{\top}(J)\hat{\bm{\beta}}(J)$. It is not difficult to show that the squared bias terms obey

[TABLE]

where $y_{J_{m}^{*}}(\mathbf{x})$, satisfying $\sharp(J_{m}^{*})=m$ and

[TABLE]

is called the best $m$-term approximation of $y(\mathbf{x})$. In addition, an argument similar to that used to prove (2.32) implies that the variance terms satisfy

[TABLE]

In view of (3.23) and (3.25), the best possible rate that can be achieved by a forward inclusion method accompanied by a stopping criterion is the same as that of

[TABLE]

According to Lemma A1.3, (A1.11), and $\mathrm{E}(y(\mathbf{x})-y_{J_{m}^{*}}(\mathbf{x}))^{2}=0$ if $m\geq k_{0}$, the convergence rate of $\bar{L}_{n}(m_{n}^{*})$ under (1.4), (1.5), or (2.17) is $(\log p/n)^{1-1/2\gamma}$, $\log n\log p/n,$ or $k_{0}\log p/n$, which coincides with that of (3.4), (3.5), or (3.6), respectively. We therefore conclude that the bounds obtained in Theorem 3.1 are rate optimal. In this connection, we also note that when (1.1) is a stationary AR($p$) model with $p\gg n$, the set of candidate models are usually given by AR(1)$,\ldots,$ AR($K_{n}$), with $K_{n}$ approaching $\infty$ at a rate slower than $n$. Unlike $\hat{J}(m),1\leq m\leq K_{n}$, the candidate set in this case is not determined by any data-driven methods, and hence the corresponding variance terms can get rid of the variance inflation factor $\log p$ (see (3.25)), which is introduced by data-dependent selection of the candidate set from all $p$ variables. As a result, the optimal rate that can be attained by an order selection criterion is equivalent to that of

[TABLE]

see Shibata (1980) for more details. Under (1.2), (1.3), or (2.17) with $N_{n}=\{1,\ldots,k_{0}\}$, the convergence rate of (3.27) is $(1/n)^{1-1/2\gamma}$, $\log n/n$, or $k_{0}/n$, which differs by a factor of $(\log p)^{1-1/2\gamma}$ from that of $\bar{L}_{n}(m_{n}^{*})$ under (1.4), (1.5), or (2.17), respectively.

We would also like to point out the differences between the current paper and the paper by Ing and Lai (2011), which investigates the performance of OGA under (1.1) with $(\mathbf{x}_{t},\varepsilon_{t})$ being i.i.d. and obeying sub-Gaussian or subexponential distributions. Note first that Theorem 1 of Ing and Lai (2011) can be understood as a special case of Theorem 2.1 when $\gamma=1$ and observations are independent over time. However, since the former theorem only focuses on the case of $\gamma=1$, its proof does not involve the approximation errors of the population OGA under general sparsity conditions such as those given in Lemma A1.1 in the Appendix. Moreover, when $\gamma=1$ is known, the optimal rate, $(\log p/n)^{1/2}$, can be achieved by choosing $m=(n/\log p)^{1/2}$, without recourse to any data-driven method to help determine the number of iterations. Alternatively, Theorem 3.1 encompasses a much wider class of sparsity conditions, and demonstrates that HDAIC can automatically choose a suitable $m$, leading to the optimal balance between the squared bias term and the variance term, without knowing the degree of sparseness. Indeed, Theorem 4 of Ing and Lai (2011) has suggested using a high-dimensional information criterion (whose penalty is heavier than that of HDAIC) to decide the number of OGA iterations when the regression coefficients satisfy the strong sparsity condition, (2.17), and a beta-min condition. Theorem 5 of Ing and Lai (2011) further introduces a backward elimination method based on the aforementioned information criterion to remove possible redundant variables surviving the first two (variable) screening stages, and shows that the resultant set of variables is equivalent to $N_{n}$ with probability tending to 1. Although the approaches adopted in both papers can be considered similar to a certain extent, their goals are entirely different. In particular, whereas Ing and Lai (2011) aim to establish selection consistency under the strong sparsity condition, this paper focuses on prediction efficiency under much more general sparsity conditions, which include the strong sparsity one as a special case. From a technical point of view, the main differences between the two papers are: (i) serial correlation is not allowed in Ing and Lai (2011); and (ii) the squared bias term in Theorems 4 (or Theorem 5) of Ing and Lai (2011) completely vanishes along the OGA path in the sense that

[TABLE]

which is ensured by the sure screening property of OGA under the strong sparsity condition (see Theorem 3 of Ing and Lai (2011)), but the squared bias term in Theorem 3.1 decays at a variety of unknown rates and can never be zero along the OGA path, making it much harder to pursue the bias-variance tradeoff along this data-driven path.

We close this section by mentioning that while condition (3.3) on $s_{a}$ involves unknown parameters, we have introduced a data-driven method for determining $s_{a}$ in Section S3 of the supplementary document, which is of practical relevance.

3.3 Proof of Theorem 3.1

We only prove (3.4). The proof of (3.5) is similar to that of (3.4), and hence is omitted. The proof of (3.6) is slightly different, and is deferred to the supplementary material because of space constraints. In the rest of the proof, a weaker restriction on the penalty term,

[TABLE]

is used instead of (3.3), although the latter one is required in the proof of (3.6).

By making use of (2.11), (2.14) (which is ensured by (2.13)), and (2.19), we show in Section A2 in the Appendix that for any $1\leq m\leq K_{n}$,

[TABLE]

where $C_{M,\gamma,\lambda_{1}}=(M+1)^{2}C^{2}_{\gamma}\lambda_{1}^{-(2\gamma-2)/(2\gamma-1)}$, with $M$ defined in (2.13), $R_{1,p}=\max_{1\leq i,l\leq p}|n^{-1}\sum_{t=1}^{n}z_{ti}z_{tl}-\rho_{il}|$, $\varepsilon_{t}(J)=y_{t}-\varepsilon_{t}-\bm{\beta}^{*^{\top}}(J)\mathbf{z}_{t}(J)$, and $\varepsilon(J)=y(\mathbf{x})-y_{J}(\mathbf{x})=y(\mathbf{x})-\bm{\beta}^{*^{\top}}(J)\mathbf{z}(J)$. In addition, it is shown in Section S2 of the supplementary material that

[TABLE]

[TABLE]

and

[TABLE]

where $R_{2,p}=\max_{1\leq i\leq p}|n^{-1}\sum_{t=1}^{n}z_{ti}\varepsilon_{t}|$ and $\|\bm{\nu}\|^{2}_{\bm{A}}=\bm{\nu}^{\top}\bm{A}\bm{\nu}$ for vector $\bm{\nu}$ and non-negative definite matrix $\bm{A}$.

Let $m^{*}_{n}=\min\{(n/\log p)^{1/2\gamma},K_{n}\}$ and

[TABLE]

in which $G\gg C_{2}$ and $C2$ is defined in (2.27). Using (3.29)–(3.32), we next show that

[TABLE]

Since (2.27) implies $\mathrm{E}_{n}(\varepsilon^{2}(\hat{J}_{m_{n}^{*}}))\leq C_{2}m_{n}^{*^{-2\gamma+1}}\leq Gm^{*^{-2\gamma+1}}_{n}$ on $A_{n}(K_{n})$, it follows that $m^{*}_{n}\geq\tilde{k}_{n}$ on $A_{n}(K_{n})$. By (2.30), one obtains

[TABLE]

where

[TABLE]

By (2.2), (2.3), and (2.28),

[TABLE]

where

[TABLE]

Moreover, (3.29)–(3.32), (2.28), and (2.27) imply that for $1\leq k\leq\tilde{k}_{n}-1$ and all large $n$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

By (3.37)–(3.42), it follows that for large enough $G$ in (3.33), there exists $0<\iota<1/2$ such that for all large $n$,

[TABLE]

In addition, (2.2), (2.3), (A3), and $\log p/n\leq m^{*^{-2\gamma+1}}_{n}$ ensure that there exists $\bar{M}_{2}>0$ such that

[TABLE]

By (2.30), (3.36), (3.44), and selecting $G$ in (3.43) larger than $2\bar{M}_{2}$, we obtain the desired conclusion (3.34).

Using (3.29)–(3.32) again, it is shown in Section S2 of the supplementary material that

[TABLE]

where $V$ is a sufficiently large constant to be specified in the proof of (3.45). With the help of (3.34) and (3.45), the desired conclusion follows if one can show that for $\gamma>1$,

[TABLE]

and for $\gamma=1$,

[TABLE]

To show (3.46), note first that

[TABLE]

where $\mathbf{L}(J)=n^{-1}\sum_{t=1}^{n}\mathbf{z}_{t}(J)(\varepsilon_{t}+\varepsilon_{t}(J))$. By (A3), (2.2), (2.3), (2.13), and straightforward algebraic manipulations, it holds that

[TABLE]

Moreover, we have

[TABLE]

and

[TABLE]

where the equality is ensured by (2.3) and (2.20). Consequently, (3.46) follows from (3.48)–(3.51), (2.30), and (2.28). The proof of (3.47) is similar to that of (3.46). The details are omitted. $\Box$

4 Conclusions

This paper has addressed the important problem of selecting high-dimensional linear regression models with dependent observations when knowledge is lacking about the degree of sparseness of the true model. When the true model is known to be an AR model or a regression model whose predictor variables have been ranked a priori based on their importance, this type of problem has been tackled in the past; see, e.g., Ing (2007), Yang (2007), Zhang and Yang (2015), and Ding et al. (2018). These authors have proposed various ways to combine the strengths of AIC and BIC and shown that their methods achieve the optimal rate without knowing whether (1.2), (1.3), or (2.17), with $N_{n}=\{1,\ldots,k_{0}\}$, is true. Their approaches, however, are not applicable to situations where the predictor variables have no natural ordering or their importance ranks are unknown. To alleviate this difficulty, we first use OGA to rank predictor variables, and then choose along the OGA path the model that has the smallest HDAIC value. Our approach is not only computationally feasible, but also rate optimal without the need for knowing how sparse the underlying time series model is.

Compared to a similar attempt made in Negahban et al. (2012), in which Lasso is used instead of OGA+HDAIC, the novelty of this paper is threefold: first, the validity of OGA+HDAIC is established not only for independent data, but also for time series data; second, the advantage of OGA+HDAIC is obtained in the important special case (1.5), which is seldom discussed in the high-dimensional literature; third, in another important special case (1.4), it is shown that OGA+HDAIC can have a faster convergence rate than Lasso. Finally, we note that OGA is exclusive for linear models. The counterpart of OGA in nonlinear models is the Chebyshev greedy algorithm (CGA) (Temlyakov, 2015). Investigating the performance of CGA+HDAIC in high-dimensional nonlinear time series models would be an interesting topic for future research.

Appendix

A1 Rates of Convergence of the Population OGA

In this section, we consider the population counterpart of OGA, whose convergence rate plays a crucial role in the analysis of the first term on the right-hand side of (LABEL:aug14). Let $0<\xi\leq 1$ be given. The algorithm initializes $J_{\xi,0}=\emptyset$. For $m\geq 1$, $J_{\xi,m}$ is recursively updated by

[TABLE]

where $j_{\xi,m}$ is any element $l$ in ${\cal P}$ satisfying

[TABLE]

with $u_{0}=y(\mathbf{x})$ and $u_{m}=y(\mathbf{x})-y_{J_{\xi,m}}(\mathbf{x})$ if $m\geq 1$. Because the algorithm is implemented based on the ‘population’ correlations of $\mathbf{x}$, it is referred to as the population OGA when $\xi=1$, and the population weak OGA when $0<\xi<1$. The following lemma provides a rate of convergence of the $\mathrm{E}(u^{2}_{m})$ under (A3) or (A4).

Lemma A1.1**.**

Assume (2.11) and (2.19). Then, there exists $G_{1}>0$ such that

[TABLE]

Moreover, if (2.12) holds instead of (2.11), then there exist $G_{2},G_{3}>0$ such that

[TABLE]

Proof. Straightforward calculations yield

[TABLE]

recalling $\mu_{J,i}=\mathrm{E}[(y(\mathbf{x})-y_{J}(\mathbf{x}))z_{i}]$. In addition, (2.19) implies

[TABLE]

By (A1.3), (A1.4) and (2.11), it follows that

[TABLE]

and hence

[TABLE]

In view of (A1.6), one has

[TABLE]

The desired conclusion (A1.1) follows from (A1.7) and Lemma 1 of Gao et al. (2013).

To show (A1.2), note first that (2.12), (A1.3) and (A1.4) yield

[TABLE]

This and an argument similar to that used in (A1.7) imply

[TABLE]

Since $M_{1}>1$, $0<\lambda_{1}\leq 1$, and $0<\xi\leq 1$, (A1.8) leads directly to (A1.2). $\Box$

Lemma A1.2 shows that (1.4) is a special case of (2.11). Using Lemmas A1.1 and A1.2, Lemma A1.3 demonstrates that the rate $m^{-2\gamma+1}$ obtained in (A1.1) cannot be improved under (1.4). More specifically, recall the best $m$-term approximation, $y_{J^{*}_{m}}(\mathbf{x})$, of $y(\mathbf{x})$ (see (3.24)). Lemma A1.3 asserts that when (1.4) and (2.19) hold true, the approximation errors of $y_{J_{\xi,m}}(\mathbf{x})$ and $y_{J^{*}_{m}}(\mathbf{x})$ only differ by a positive constant.

Lemma A1.2**.**

Suppose that (1.4) is true for some $\gamma>1$. Then (2.11) holds for the same $\gamma$.

Proof. See Section S1 of the supplementary document.

Lemma A1.3**.**

Suppose that (1.4) holds for some $\gamma>1$ and (2.19) is true. Then, for all $1\leq m\leq(1-\epsilon)p$, where $\epsilon$ is an arbitrarily small positive constant, there exist $D_{1},D_{2}$ and $D_{3}$ such that

[TABLE]

and

[TABLE]

Proof. By Lemmas A1.1 and A1.2, it follows that for all $1\leq m\leq(1-\epsilon)p$,

[TABLE]

where $J^{o}_{m}$ is the index set corresponding to $\{\beta^{2}_{(1)},\cdots,\beta^{2}_{(m)}\}$ and $\underline{d}>0$ depends only on $\gamma$ and $\epsilon$. These inequalities lead immediately to (A1.9) and (A1.10). $\Box$

Remark A.1. Theorem 2.1 of Temlyakov (1998) shows that a near best $m$-term approximation can be realized by a greedy-type algorithm under a basis $L_{p}$-equivalent to the Haar basis. Since the Haar basis yields an identity correlation matrix, our correlation assumption, (2.19), appears to be substantially weaker. The performance of the $m$-term approximation of OGA has been investigated by Tropp (2004) under a noise-free underdetermined system and a condition on the cumulative coherence function, which requires that the atoms in the dictionary are ‘nearly’ uncorrelated. His approximation error for OGA is larger than that of the best $m$-term approximation by a factor of $(1+6m)^{1/2}$. Suppose that (1.5) holds. Then,

[TABLE]

which, together with (1.5) and Minkowski’s inequality, yields

[TABLE]

where $C_{1,\beta}\leq C_{2,\beta}$ are some positive constants depending on $\beta$. On the other hand, the argument used to prove (A1.2) leads to

[TABLE]

where $f_{{\rm oga}}=\xi^{2}\lambda_{1}(L1/U_{1})^{2}(1-{\rm exp}(-\beta))^{2}<2\beta$. Equations (A1.11) and (A1.12) suggest that the population OGA and the best $m$-term approximation in general do not share the same convergence rate in the exponential decay case. To be as efficient as the best $m$-term approximation, the population OGA needs to run for another $m(2\beta/f_{{\rm oga}}-1)$ iterations, which is still of order $m$.

A2 Proof of (3.29)

Recall that (2.13) implies (2.14) with $C=M+1$. This, (2.11), and (2.19) yield

[TABLE]

Thus, (3.29) follows.

SUPPLEMENTARY MATERIAL

Supplement to “Model selection for high-dimensional linear regression with dependent observations”. The supplementary material contains the proofs of (2.7), (2.28), (2.30), (2.32), (3.6), (3.30)–(3.32), (3.45), and Lemma A1.2, and a simulation study to demonstrate the performance of OGA+HDAIC under a high-dimensional ARX model whose $\bm{\Gamma}$ obeys $\lambda_{\max}(\bm{\Gamma})\to\infty$ and (2.19).

Acknowledgments

I thank Hai-Tang Chiou, Hsueh-Han Huang, and Tze Leung Lai for helpful suggestions on earlier versions of this paper. I would also like to thank an Associate Editor and two anonymous referees for insightful and constructive comments.