High Dimensional Time Series Regression Models: Applications to Statistical Learning Methods

High Dimensional Time Series Regression Models: Applications to Statistical Learning Methods

Introduction

These lecture notes series present a unified framework for recent developments of deep learning methodologies in time series econometrics and statistics. Specifically, statistical machine learning methods and theory are now widespread in social sciences, economics\footnote{The idea of representing the dynamic interactions of economic agents similar to the structure of brain neural networks dates back to ashby1957introduction. Furthermore the application of neural networks for modelling purposes dates back to the fundamental works of hebb1949organization and rosenblatt1958perceptron.}, finance and mathematical modelling. Neural networks can be considered as a type of linear sieve estimation where the basis functions themselves are flexibly learned from the data by optimizing over many combinations of simple functions (see, farrell2021deep). Deep reinforcement learning has gain attention in recent years especially when modelling dynamic interactions of economic agents. In particular deep neural nets contain many hidden layers of neurons between the input and output layers, which has been found to exhibit superior performance across a variety of contexts. Our goal is to consolidate the theoretical and practical underpinnings of machine learning methods when modelling time series data for estimation, inference and forecasting purposes.

In order to express precise conditions for ergodicity, we turn to the Lyapounov exponent, a concept which is well-known by those studying stability of dynamical systems. Therefore, in the context of nonlinear time series, nonstability means explosive behaviour. Thus, the Lyapounov exponent, as we define it for the state space model $\left\{ X_t \right\}$ of a time series, is given by

which measures the drift of the process when the sample size is large. Define with

The sets $\ell_1, \ell_2$ and $\ell_{ \infty }$ are all vector spaces and it holds that $\left\lVert ( x_n ) \right\rVert_1 = \sum_{ n \geq 1 } | x_n |$, $\left\lVert ( x_n ) \right\rVert_2 = \sum_{ n \geq 1 } \left( \sum_{ n \geq 1 } x_n^2 \right)^{ 1 / 2 }$ define norms on $\ell_1$ and $\ell_2$ spaces respectively. For any $x \in \mathbb{R}^n$, $\left\lVert x \right\rVert = \sqrt{ \sum_{i=1}^n x_i^2 }$, $\left\lVert x \right\rVert_{ \ell_1 } = \sum_{ i=1 }^n | x_i |$ and $\left\lVert x \right\rVert_{ \infty } = \underset{ 1 \leq i \leq n }{ \mathsf{max} } | x_i |$, denote $\ell_2$, $\ell_1$ and $\ell_{ \infty }-$norms respectively. We also denote with $\left\lVert x \right\rVert_{ \ell_0 } = \sum_{ i=1 }^n \mathbf{1} \left\{ x_i \neq 0 \right\}$ which is simply the number of non-zero entries of $x$. For any $n \times n$ matrix $M$, $\left\lVert M \right\rVert_{ \ell_{ \infty } } = \mathsf{max}_{ 1 \leq i \leq n } \sum_{ j = 1 }^n | M_{i,j} |$ denotes the induced $\ell_{ \infty }$ matrix norm. Detailed statistical theory for high dimensional statistics applications are presented in wainwright2019high (see also vershynin2018high). Let $\left\{ X_{t,n} \right\}_{ t = 1}^n$, $n \in \mathbb{N}$ be a random array and $\left\{ \mathcal{F}_t \right\}_{t = - \infty}^{ \infty }$ a filtration such that $X_{t,n}$ is $\mathcal{F}_t-$measurable for all $t$ and $n$. Based on the aforementioned notation we discuss important theoretical results from the statistics and econometrics literature.

Limit theory for High Dimensional Dependent Data

Limit theory of Covariance Matrices for Linear Processes

Covariance and Precision matrix estimation for high dimensional time series

Following the framework of Chen2013covariance, suppose we have $n$ temporally observed $p-$dimensional vectors $\left( \mathbf{z}_i \right)_{ i = 1}^n$ having mean zero and covariance matrix $\Sigma_i = \mathbb{E} \left( \mathbf{z}_i, \mathbf{z}_i \right)$ whose dimension is $p \times p$. Our goal is to estimate the covariance matrices $\Sigma_i$ and their inverses $\Omega_i = \Sigma_i^{-1}$ based on the data matrix $Z_{ p \times n} = \left( \mathbf{z}_1,..., \mathbf{z}_n \right)$. In the classical situation where $p$ is fixed, $n \to \infty$ and $\mathbf{z}_i$ are mean zero independent and identically distributed $\textit{i.i.d}$ random vectors, it is well known that the sample covariance matrix

is a consistent and well behaved estimator of $\Sigma$, and $\hat{\Omega}_n = \hat{\Sigma}_n^{-1}$ is a natural and good estimator of $\Omega$.

However, when the dimensionality $p$ grows with $n$, random matrix theory asserts that $\hat{\Sigma}_n$ is no longer a consistent estimate of $\Sigma$ in the sense that its eigenvalues do not converge to those of $\Sigma$, as the Marcenko-Pastur law states. Moreover, it is clear that $\hat{\Omega}_n$ is not defined when $\hat{\Sigma}_n$ is not invertible in the high-dimensional case with $p >> n$.

Furthermore, if $\hat{\lambda}_j \leq 0$, since $\mu_i \geq 0$, we have that $\left| \hat{\lambda}_j \right| \leq \left| \hat{\lambda}_j - \mu_i \right|$. Then,

Notice that the eigenvalues of $\tilde{S}_v$ are bounded below by $v$, and thus it is positive definite such that

The same positive-definization procedure also applies to the spectral norm and its rate can be similarly preserved.

\paragraph{Sparsity Conditions}

We begin by presenting the commonly used sparsity condition defined in terms of the strong $\ell^q-$ball such that

In particular, when $r = 0$, the sparsity condition implies that $\underset{ 1 \leq k \leq p }{ \text{max} } \sum_{ j = 1 }^p \mathbf{1} \left\{ \sigma_{jk} \neq 0 \right\} \leq \tilde{M}$.

Take the AR(1) process, $\mathbf{z}_i = A \mathbf{z}_{i-1} + \mathbf{e}_i$, where $A$ is the real matrix with spectral norm $\rho(A) < 1$, it is of the form as in (ref) with $A_j = A^j$. and the functional dependence measure $\theta_{ i, 2q, j} = \mathcal{O} \left( \rho(A)^i \right)$.

Furthermore, the structure of the proposed risk matrix has different properties in contrast to the covariance matrix. For instance, the estimation procedure is not invariant to variable permutations. In contrast the covariance matrix itself has the property of being invariant to variable permutations. Therefore, this property allows the implementation of high dimensional techniques in inference problems. For instance, in the case of the covariance matrix ledoit2004well proposed a way to compute an optimal linear combination of the sample covariance with the identity matrix, which also results in shrinkage of eigenvalues. Furthermore, shrinkage estimators are invariant to variable permutations, but they do not affect the eigenvectors of the covariance, only the eigenvalues, and it has been shown that the sample eigenvectors are also not consistent when $p$ is large. Therefore, developing also estimation methods which are invariant to variable permutations is important.

Main Results on Probability Bounds

Following the framework of Chen2013covariance we focus in how to establish the convergence theory for covariance matrix estimates, we shall use the functional dependence measure of Wu (2005). Recall that $Z_{ji} = g_j ( \mathcal{F}_i ), 1 \leq j \leq p$, where $g_j(.)$ is the $j$th coordinate projection of the $\mathbb{R}^p-$valued measurable function $\mathbf{g}$. For $w > 0$, the functional dependence measure of $Z_{ji}$ is defined by

where $Z_{ji}^{\prime} = g_j ( \mathcal{F}_i^{\prime} ), \mathcal{F}_i^{\prime} = \left( ..., \mathbf{e}_{-1}, \mathbf{e}_{0}^{\prime}, \mathbf{e}_{1},..., \mathbf{e}_{i} \right)$ and $\mathbf{e}_{0}^{\prime}$ and $\mathbf{e}_{0}^{\prime}$, $\mathbf{e}_{\ell}$, $\ell \in \mathbb{Z}$, are i.i.d. In other words, $Z_{ji}^{\prime}$ is a coupled version of $Z_{ji}$ with $\mathbf{e}_{0}$ in the latter replaced by an i.i.d copy of $\mathbf{e}_{0}^{\prime}$.

\paragraph{Proof of Theorem 2.1} (Chen2013covariance) We first assume that $\alpha > 1/2 - 1/q$. Notice that

Since the functional dependence measure for the product process $\left( Z_{ji} Z_{ki} \right)_i$ has a bounded inequality, under the decay condition $\Theta_{ m, 2q } \leq C m^{- \alpha }, \alpha > 1/2 - 1/q$. we have by Theorem 2(ii) in Wu (2013) that

holds for all $v > 0$. Using integration by parts we obtain

Moreover, some useful limit results are also presented in the paper of xiao2013asymptotic. In particular, the sample covariance between columns $x_i$ and $x_j$ is defined as

In high-dimensional covariance inference, a fundamental problem is to establish an asymptotic distributional theory for the maximum deviation

\paragraph{Proofs of (iii) and (iv)}

The first step is the truncation step.

\paragraph{A normal comparison principle}

Suppose that for each $n \geq 1$, $\left( X_{n,i} \right)_{ i \in \mathcal{I}_n }$ is a Gaussian random vector whose entries have mean zero and variance one, where $\mathcal{I}_n$ is an index set with cardinality $\left| \mathcal{I}_n \right| = s_n$. Let $\mathcal{\Sigma}_n \left( r_{n,i,j} \right)_{ i,j \in \mathcal{I}_n }$ be the covariance matrix of $\left( X_{n,i} \right)_{ i \in \mathcal{I}_n }$. Assume that $s_n \to \infty$ as $n \to \infty$. We impose either of the following two conditions.

Limiting Laws of Coherence of Random Matrices with Applications

Furthermore, following the framework of cai2011limiting, we present some useful lemmas below.

Therefore, we have that

Suppose we consider an $\mathbb{R}^d-$valued time series process $\left\{ \underline{X}_t, t \in \mathbb{Z} \right\}$ with $\underline{X}_t = \left( X_{1,t},...., X_{d,t} \right)^{\top}$, and we have data $\underline{X}_1,..., \underline{X}_n$ at hand to use for estimation and inference purposes.

Large deviations for quadratic forms

Large deviations for quadratic forms of stationary processes have been extensively studied in the literature. In particular, we need an upper bound for the tail probability under less restrictive conditions. In this section, we follow the framework proposed by xiao2012covariance. More specifically, we prove a result on probabilities of large deviations of quadratic forms of stationary processes, which take the form

The coefficients $a_{s,t} = a_{T,s,t}$ may depend on $T$. Throughout this section we assume that sup$_{s,t} \left| a_{s,t} \right| \leq 1$, and $a_{s,t} = 0$ when $| s - t | > B_T$, where $B_T \to \infty$, and $B_T = \mathcal{O} ( T^{\gamma} )$ for some $0 < \gamma < 1$.

Asymmetry Helps

In this section, we present the main results obtained by chen2021asymmetry. We being by considering the low-rank matrix completion of chen2021asymmetry.

\paragraph{Low-rank matrix completion} Suppose that $\mathbf{M}$ is generated using random partial entries of $\mathbf{M}^{\star}$ as

where $p$ denotes the fraction of the entries of $\mathbf{M}^{\star}$ being revealed. We have that $\mathbf{H} = \mathbf{M} - \mathbf{M}^{\star}$ is zero-mean and obeys $\left| H_{ij} \right| \leq \frac{ \mu }{ np } := B$ and Var$\left( H_{ij} \right) \leq \frac{ \mu^2 }{ p n^2 }$. For instance, in the case where $p \approx \frac{ \mu^2 \text{log} (n) }{ n }$, then invoking Corollaries 1-3 of chen2021asymmetry yields

which gives that

with high probability, where $\mathbf{a} \in \mathbf{R}^n$ is any fixed unit vector. Empirically, eigen-decomposition outperforms SVD in estimating both the leading eigenvalue and eigenvector of $\mathbf{M}^{\star}$

\paragraph{Why asymmetry helps?} According tochen2021asymmetry, if we consider their Theorem 3, focusing on the case with $\lambda^{\star} = 1$ for simplicity, then the key ingredient is the Neumann trick stated in Theorem 2. Specifically, in the rank-1 case we can expand

We obtain that

where the last inequality holds since (i) $\left| \mathbf{u}^{ \star \top} \mathbf{u} \right| \leq 1$, and (ii) $\lambda$ is real-valued and obeys $\lambda \approx 1$ if $ \left\lVert \mathbf{H} \right\rVert << 1$. As a result, the perturbation can be well controlled as long as $\left| \mathbf{a}^{\top} \mathbf{H}^s \mathbf{u}^{\star} \right|$ is small for every $s \leq 1$.

\paragraph{Pertrubation Analysis for the rank-r case}

\paragraph{Eigenvalue perturbation for the rank-r case} The eigenvalue perturbation analysis can be extended to accommodate the case where $\mathbf{M}^{\star}$ is symmetric and rank$-r$. As before, we assume that the $r$ nonzero eigenvalues of $\mathbf{M}^{\star}$ obey that $| \lambda_1 | \geq ... \geq | \lambda_r |$

for some sufficiently small constants $c_1 > 0$. Then for any fixed unit vector $\mathbf{a} \in \mathbb{R}^n$ and any $1 \leq \ell \leq r$, with probability at least $1 - \mathcal{O} \left( n^{-10} \right)$ one has that

The particular condition allows to control the perturbation of the linear form of eigenvectors. Thus, the perturbation upper bound grows as either the rank $r$ or the condition number $\kappa$ increases.

In comparison, the Bauer-Fike theorem (Lemma 2) together with Lemma 1 gives a perturbation bound

For the low-rank case where $r << \sqrt{n}$, the eigenvalue pertrubation bound derived in Corollary 5 can be much sharper than the Bauer-Fike theorem.

Consequently, by taking $\mathbf{a} = \mathbf{e}_i$ for $1 \leq i \leq n$ in Corollary 6, we arrive at the following statement regarding the alternative definition of the incoherence of the eigenvalue matrix $\mathbf{U}$.

\paragraph{Proof of Theorem (ref)}

Without loss of generality, we shall assume that $\lambda_{\text{max}}^{\star} = \lambda_1^{\star} = 1$ throughout the proof. To begin with, Lemma 2 implies that for all $1 \leq \ell \leq r$,

as long as $\left\lVert \mathbf{H} \right\rVert < 1 / ( 2 \kappa )$. In view of the Neumann trick (Theorem 2), we can derive that

where the third line follows since $\sum_{ j = 1}^r \left| \mathbf{u}_j^{ \star } \mathbf{u}_{\ell} \right|^2 \leq \left\lVert \mathbf{u}_{\ell} \right\rVert_2^2 = 1$, and the last inequality makes use of (49). Apply Corollary 4 in chen2021asymmetry to obtain a bound as below

with the provision that $\left| \lambda_{\ell} \right| > 1 / ( 2 \kappa )$ and $\text{max} \left\{ B \text{log}(n), \sqrt{n \sigma^2 \text{log} (n) } \right\} \leq c_1 / \kappa$ for some sufficiently small constant $c_1 > 0$. The condition $| \lambda_{\ell} | > 1 / (2 \kappa )$ follows immediately by combining Lemma 2, Lemma 1 and the condition (34).

\paragraph{Appendix}

We justify the existence and uniqueness of an eigenvalue in $\mathcal{B} \left( 1, \left\lVert \mathbf{H} \right\rVert \right)$. Denote with $\Lambda \left( \textbf{M} \right) = \left\{ \lambda_1,..., \lambda_n \right\}$ and define a set of auxiliary matrices

As we can see the set of eigenvalues of $\mathbf{M}_t$ depends continuously on $t$, we can write as below

with each $\lambda_j (t)$, $1 \leq j \leq n$ being a continuous function in $t$. Meanwhile as $\left\lVert \mathbf{H} \right\rVert < 1 / 2$ and $0 \leq t \leq 1$, the two disks $\mathcal{B} \left( 1, t \left\lVert \mathbf{H} \right\rVert \right)$ and $\mathcal{B} \left( 0, t \left\lVert \mathbf{H} \right\rVert \right)$ are always disjoint sets. Therefore, the continuity of the spectrum with respect to $t$ requires that $\lambda_j(t)$ to always stay within the same disk where $\lambda_j(0) \in \left\{ 0, 1 \right\}$ lies, namely, such that

Given that $\mathbf{M}^{\star}$ has $n - 1$ eigenvalues equal to 0 and one eigenvalue equal to 1, we establish the lemma for the rank-1 case.

\paragraph{ Rank-r } Follow the above argument for the case where we have a rank$-1$ matrix, we can immediately show that: if $\left\lVert \mathbf{H} \right\rVert < \lambda_r^{ \star } / 2$, then (i) there are exactly $n - r$ eigenvalues lying within $\mathcal{B} \left( 0, \left\lVert \mathbf{H} \right\rVert \right)$, (ii) all other eigenvalues lie within $\cup_{ 1 \leq j \leq r} \mathcal{B} \left( \lambda_j^{\star} , \left\lVert \mathbf{H} \right\rVert \right)$, which are exactly top$-r$ leading eigenvalues of $\mathbf{M}$.

\paragraph{Proof of Theorem 2 }[chen2021asymmetry] Following the definition of eigenvectors, we have that

When $\left\lVert \mathbf{H} \right\rVert_2 < \left| \lambda_{\ell} \right|$, one can invert $\left( \mathbf{I} - \frac{1}{ \lambda_{\ell} } \mathbf{H} \right)$ and therefore we obtain that

where the last term follows by rearranging terms. Finally, replacing $\left( \mathbf{I} - \frac{1}{ \lambda_{\ell}} \mathbf{H} \right)^{-1}$ with the Neumann series we obtain that $\sum_{ s = 0}^{\infty} \frac{1}{ \lambda^s } \mathbf{H}^s$, we establish the theorem.

\paragraph{Proof of Lemma 3}[chen2021asymmetry] We start with the rank-1 case. Towards this, we resort to the Neumann trick in Theorem 2, which in the rank-1 case gives the following expression

From Lemma 2, we know that $\lambda$ is real-valued and that $\lambda > 1 - \left\lVert \mathbf{H} \right\rVert \geq 3/4 > \left\lVert \mathbf{H} \right\rVert$ under our assumption. This together yields that

where the last inequality holds since $\lambda \geq 3 / 4$ and $\lambda - \left\lVert\mathbf{H}\right\rVert \geq 1 - 2 \left\lVert\mathbf{H}\right\rVert \geq 1 / 2 $. Next, by decomposing $\mathbf{u}$ into two orthogonal components

we obtain that

the particular inequality holds since $\left( \mathbf{u}^{\star \top} \mathbf{u} \right) \mathbf{u}^{\star}$ is orthogonal projection of $\mathbf{u}$ onto the subspace spanned by $\mathbf{u}^{\star}$, and hence

Moreover, since $\mathbf{u}$ is real-valued we obtain the following bound:

Next, we consider the rank$-r$ case, since for any $1 \leq r \leq r$ we have that

where the inequality arises since $\sum_{j=1}^r \left( \mathbf{u}_j^{\star \top} \mathbf{u}_{\ell} \right) \mathbf{u}_j^{\star}$ is the Euclidean projection of $\mathbf{u}_{\ell}$ onto the space of $\left\{ \mathbf{u}_{1}^{\star},..., \mathbf{u}_{r}^{\star} \right\}$. Furthermore, we observe that

This taken collectively with Theorem 2 leads to

The last two lines follow since, when $\left\lVert \mathbf{H} \right\rVert < 1 / (4 \kappa )$. We can show that $| \lambda_{\ell} | > | \lambda_{ \text{min} }^{ \star } | - \left\lVert \mathbf{H} \right\rVert \geq 3 / (4 \kappa )$ and (ii) $\left| \lambda_{\ell} \right| - \left\lVert \mathbf{H } \right\rVert \geq \left| \lambda_{ \text{min} }^{ \star } \right| - 2 \left\lVert \textbf{H} \right\rVert \geq 1 / (2 \kappa)$.

\paragraph{Proof of Lemma 5 on Edge and Vertices}[chen2021asymmetry]

To establish this lemma, we exploit entrywise independence of $\mathbf{H}$ and develop a combinatorial trick. To begin with, we expand the quantity of interest as below

In this section, we use that $ \mathcal{I} := \left\{ i_t^{(b)} | 0 < t < s, 1 < b < k \right\} \in [ n ]^{ (s+1) k } $, to denote such a collection of $( s + 1 ) k$ indices. Therefore, one can write that

Further details of these proofs can be found in the Appendix of chen2021asymmetry.

Limit theory for Functional Data

Simultaneous Diagonalization

In this section we briefly study the framework proposed by yuan2010reproducing. In particular, before studying the asymptotic properties of the regularized estimators, we first investigate the relationship between the eigen structures of the covariance operator for $X(.)$ and the reproducing kernel of the functional space $\mathcal{H}$. As observed in earlier studies, eigenvector structures play prominent roles in determining the nature of the estimation problem in functional linear regression.

Recall that $K$ is the reproducing kernel of $\mathcal{H}_1$. Because $K$ is continuous and square integrable, it follows from Mercer's theorem that $K$ admits the following spectral decomposition

Here $\rho_1 \geq \rho_2 \geq ... $ are the eigenvalues of $K$, and $\left\{ \psi_1, \psi_2,... \right\}$ are the corresponding eigenfunctions,

Moreover,

where $\delta_{ij}$ is the Kronecker's delta.

\paragraph{Convergence rates}

We now turn to the asymptotic properties of the smoothness regularized estimators. To fix ideas, in what follows, we shall focus on the squared error loss. Recall in this case we have that

Furthermore, we have that

Observe that

Write

For instance, we have the following decomposition

we refer to the terms above on the right hand side as the stochastic error and deterministic error respectively.

\paragraph{Deterministic Error}

Write

Then, Theorem 3 above implies that

Therefore,

Common Functional Principal Components

According to benko2009common, the elements of the risk matrix $\Gamma$ are functionals of the model estimates which can be estimated with asymptotically negligible bias and a parametric rate of convergence $T_i^{- 1 / 2}$. For example when the data are generated from a balanced, equidistance design, then it can be easily seen that for $i \neq j$ this rate of convergence is achieved by the following estimator. Notice that a bias corrected estimator may yield negative eigenvalues. However, in practise these values are small and can be interpreted close to zero. For instance, when evaluating eigenvalues and eigenfunctions of the empirical covariance operator of nonparametrically estimated curves $\hat{X}_i$, then for fixed $r \leq r_0$ the above rate of convergence for the estimated eigenfunctions may well be achieved for a suitable choice of smoothing parameters. For example, when using standard methods it does not seem to be possible to obtain a corresponding rate of convergence, since any smoothing bias $\left| \mathbb{E} \left[ \hat{X}_i \right] - X_i(t) \right|$ will invariably affect the quality of the corresponding estimate of $\hat{\lambda}_r$.

\paragraph{Main Results}

In the following, we have that $\left\lVert v \right\rVert = \left( \int_0^1 v(t)^2 dt \right)^{1 / 2}$ will denote the $L^2-$norm for any square integrable function $v$. Consider that $\hat{\Sigma}_m$ is the $m \times m$ matrix with elements given by

Let $\lambda_1 \left( \hat{\Sigma}_m \right) \geq \lambda_2 \left( \hat{\Sigma}_m \right) \geq ... \geq \lambda_m \left( \hat{\Sigma}_m \right)$ and $\hat{\zeta}_{1,m},...,\hat{\zeta}_{m,n}$ denote eigenvalues and corresponding eigenvectors of $\hat{\Sigma}_m$. Therefore, some straightforward algebra shows that

We will use $\Sigma_m$ to represent the $m \times m$ diagonal matrix with diagonal entries given by $\lambda_1 \leq ... \leq \lambda_m$. Then, the corresponding eigevectors are given by the $m-$dimensional unit vectors denoted by $e_{1,m},..., e_{m,m}$. Then, by Lemma A of Kneip and Utikal we obtain that the differences between the eigenvalues and eigenvectors of $\Sigma_m$ and $\hat{\Sigma}_m$ can be bounded by

Moreover, we have that

where we denote with $S_{r,m} = \sum_{ s \neq r } \frac{1}{ \lambda_s - \lambda_r } e_{s,m} e_{s,m}^{\top}$.

Assumption 1 implies that $\mathbb{E} \left( \hat{ \beta }_r \right) = 0$ and Var$\left( \hat{ \beta }_r \right) = \frac{ \lambda_r }{ r }$ and with $\delta_{ii} = 1$ as well as $\delta_{ij} = 0$ for $i \neq j$, we obtain that

for all $m$. Furthermore, since

Therefore after applying the central limit theorem we obtain that

Obviously the event $\hat{\lambda}_{r-1} > \hat{\lambda}_{r} > \hat{\lambda}_{r+1}$ occurs with probability 1.

Limit theorems for Dependent Sequences

In this Section, we consider some important results presented in the literature for dependent sequences. We consider a "high-dimensional" probability space $\left( \Omega, \mathcal{F}, \mathbb{P} \right)$. For functions of dependent random variables $( X_i )_{ i \in \mathbb{N} }$ the main challenge is often to quantify and bound the dependence among random variables $X_i$, in terms of various types of mixing coefficients (see, kontorovich2008concentration). To establish concentration bounds a sufficiently rapid decay of the mixing coefficients is necessary.

Method of bounded martingale differences

The particular methodology presented in the paper of kontorovich2008concentration provides a suitable mechanism of establishing concentration inequalities (see also doukhan1995invariance).

Concentration Inequalities for Dependent Random Variables

We consider the concept of pointwise inequalities, that is, inequalities that hold uniformly for any parameter $\theta \in \Theta$. Define the function (see, van2002hoeffding)

For a real-valued random variable $\xi$, we define with

Moreover, we write that $\xi \in \mathcal{L}^q$ for some $q > 0$ if it holds that

Central limit theorems for high dimensional dependent data

Following the framework proposed by chang2021central , recall that we define with $S_{n,x} = n^{-1/2} \sum_{t=1}^n X_t$. Let $\mathcal{G} \sim \mathcal{N} (0, \Xi)$ where $\Xi := \mathsf{Cov} \left( n^{-1/2} \sum_{t=1}^n X_t \right)$. Without loss of generality we assume that $\mathcal{G}$ is independent of $\mathcal{X} = \left\{ X_1,..., X_n \right\}$. We write with $X_t = \left( X_{t,1},..., X_{t,p} \right)^{\prime}$. Then, the long-run variance of the $j-$th coordinate marginal sequence $\left\{ X_{t,j} \right\}_{t=1}^n$ is defined as

Therefore, in order to determine the convergence rate of $\rho_n$ for the $\alpha-$mixing sequence $\left\{ X_t \right\}$, we impose additional regularity conditions. The above condition assumes that the partial sum $\frac{1}{\sqrt{n} } \sum_{t=1}^n X_{t,j}$ is non-degenerated which is necessary to bound the probability of a Gaussian vector taking values in a small region. When $\left\{ X_{t,j} \right\}_{ t \geq 1 }$ is stationary, then it holds that

where $\Gamma_k(k) = \mathsf{Cov} \left( X_{1,j}, X_{k+1,j} \right)$ is the autocovariance of $\left\{ X_{t,j} \right\}_{ t \geq 1}$ at lag $k$.

Due to the fact that $x \mapsto \mathsf{exp} \left[ \left( x \big/ {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert X \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{ e^{\gamma} } \right)^{\gamma} \right]$ is non-decreasing then

It holds that,

where $K_{\gamma} := \frac{ ( \mathsf{exp} a_{\gamma}^{\gamma} - 1) }{ a_{\gamma} }$ and $a_{\gamma}$ is defined as below

Moreover, it holds that

We define with $\phi(d)$ the following quantity

Sub-Weibull random vectors under $\beta-$mixing

Most analysis on lasso assume that data that have sub-Gaussian or subexponential tails. These assumptions ensure that the moment generating function exists, at least for some values of the free parameter. Nonexistence of the moment generating function is often taken as a definition of having a heavy tail.

A natural realization which allows for heavier tails is as follows. Fix some $\gamma > 0$, and require

Notice that the condition above requires that the tail is no heavier than that of a Weibull random variable with parameter $\gamma$.

Approximation Theorems for Strongly Mixing Random Variables

Discussion

The various mixing properties of time series (either contemporaneously or temporally dependent) are important for understanding how these particular features affect the asymptotic theory for model estimators, test statistics and the relevant properties of model selection methods based on dependent data which exhibit these features. Some potential data specific features such as the presence of serial dependence among errors and regressors, heteroscedasticity and fat tails can lead to misspecified time series models which especially in a high dimensional setting require to establish validity via error bound determination (see, adamek2023lasso and wong2020lasso among others).

In particular, shrinkage methodologies assume a certain structure on the unknown parameter vector of interest. Generally, the sparsity condition implies that a small but unknown subset of the high dimensional vector of covariates is considered to have "significantly different than zero" coefficients, while the remaining subset of covariates have negligible, or even exactly zero, coefficients. In other words, for Lasso shrinkage methodologies to render meaningful inference the penalty function exploits the underline sparsity condition. Specifically, adamek2023lasso obtain novel theoretical results for both point estimation and inference via the desparsified lasso. Furthermore, they consider a general time series framework where the regressors and error terms are allowed to be non-Gaussian, serially correlated and heteroscedastic and the number of variables can grow faster than the time dimension.

Time Series Regression Models with Many Covariates

A recent growing literature develops econometric frameworks for estimation and inference in regression models with many covariates. Relevant studies include among others cattaneo2018inference, karmakar2022long and wei2023inference (see, also farrell2015robust although not in a time series setting).

A forecasting Application

Let $\hat{y}_t$ be the forecast of $y_t$ based upon information up to $t-1$. Notice that when the interested of the researcher is the one-period ahead forecast, $\left( y_t - \hat{y}_t \right)^2$ is the cost to be minimized. However, there are two situations where the accumulated cost function, denoted with $\sum_{ j = 1 }^t \left( y_t - \hat{y}_t \right)^2$ is more appropriate. For example, in the sequential forecast case, the forecaster are updated sequentially over many periods and therefore the accumulated cost function is the target to be minimized.

Consider for example the one-period expected loss function $\mathbb{E}\left( y_{T+1} - \hat{y}_{T+1} \right)^2$. For instance, for the AR(1) model, under the assumption that $\mathbb{E} \left( \epsilon_t^2 | \mathcal{F}_{t-1} \right) = \sigma^2$ almost surely for all $t$, then it can be shown that under appropriate assumptions that $\frac{1}{T} \sum_{t=1}^T \left( y_{T+1} - \hat{y}_{T+1} \right)^2$ almost surely.

Let $\hat{\beta}_t$ be the Least square estimate of $\beta$, then we have that

where $Y_t = \left\{ y_1,...,y_t \right\}^{\prime}$, then $\hat{y}_t = \hat{\beta}_{t-1}^{\prime} Y_{t-1}$ is the least square prediction of $y_t$ at time $t-1$.

Best Linear Predictor

The main tool of prediction, or forecasting, in time series relies on the concept of the best linear predictor which is often employed in order to obtain the limiting distribution of an arbitraty size vector of autocovariance estimators (see, Dhrymes2013mathematics).

The first order conditions are given by

Thus, by rearranging we get that

or $\alpha = C_n(i-j) c_{nh}$.

Forecast Combination

Assume that we are interested to construct a forecast combination criterion based a specific econometric model of interest. Then, a simple approach would be to consider the construction of a criterion which uses OOS forecast combinations. Let $\mathbf{f}_t$ be a sequence of out-of-sample forecasts for $y_{t+1}$ derived using a set $m$ different econometric models. The combination forecast is then defined as $f( \mathbf{w} ) = \mathbf{w}^{'} \mathbf{f}$.

Therefore, the OOS empirical MSFE is given by

Furthermore, the Granger-Ramanathan forecast combination method implies to select $\mathbf{w}$ to minimize the OOS MSFE. However, minimization over $\mathbf{w}$ is equivalent to the least-squares regression over $y_t$ on the forecasts $y_{t+1} = \mathbf{w}^{'} \mathbf{f} + \epsilon_{t+1}$. Then, the unrestricted least-squares gives the following vector of weights

Notice that the above unrestricted least-squares approach can produce weights which are far outside $[0,1]$, thus an alternative representation can be constructed by subtracting $y_{t+1}$ from each side, such that $0 = \mathbf{w}^{'} \mathbf{f} - y_{t+1} + \epsilon_{t+1}$, and by defining $\mathbf{e}_{t+1} = \left( \mathbf{f}_t - y_{t+1} \right)$ to be the negative forecast errors. Then, it holds that $ 0 = \mathbf{w}^{'} \mathbf{e}_{t+1} + \epsilon_{t+1}$, which is the regression of 0 on the forecast errors. Therefore, the constrained GR weights solve the following problem.

where $\mathbf{A} = \displaystyle \sum_t \mathbf{e}_{t+1} \mathbf{e}^{'}_{t+1}$, the $ M \times M $ matrix of forecast error empirical variances/covariances.

Relevant studies include hansen2007least, hansen2008least and cheng2015forecasting while a recent approach related to forecasting can be found in boot2019forecasting and swanson2020predicting. For nonstationary time series models relevant studies include hansen2010averaging, coroneo2020comparing and kejriwal2021generalized among others.

Statistical Theory for Lasso Regression Models

Statistical Principles

Recall that the Euclidean norm is given by: $\left\lVertb\right\rVert_2 = \sqrt{ b^{\prime} b }$. In other words, the inverse of $\left( \boldsymbol{X}^{\prime} \boldsymbol{X} \right)$ is regularized by $n \lambda \boldsymbol{I}_k$ (called Tikhonov regularization for ill-posed problems).

Special cases of Ridge Regression include:

Therefore, it holds that

Based on the above results we can introduce the model selection idea with an $\ell_0-$penalty. Thus, in the literature various studies mention that one of the main advantages of the Lasso shrinkage is that it "bets on sparsity". Therefore, using an $\ell_0$ penalty term we have that

which implies that

where $\lambda \geq 0$ is the tuning (penalty) parameter which needs to be chosen. Similarly, using an $\ell_1$ penalty

where $\lambda \geq 0$ is the tuning (penalty) parameter which needs to be chosen.

The lasso shrinkage approach was proposed in the seminal study of tibshirani1996regression. Various studies examine its statistical properties such as fan2001variable and zhang2010nearly among others, while many variants of the lasso shrinkage have been proposed in the literature such as zou2005regularization, zou2006adaptive, huang2008adaptive, park2008bayesian and james2009dasso. A relevant empirical study is present by katsouris2021forecast. Recently, the literature focuses on applications of the lasso shrinkage approach with group structure (see, also huang2012selective and bing2022inference).

In particular, klau2018priority propose the priority lasso, which considers the use of different penalty across a baseline group of covariates (see, breheny2015group for case with grouped predictors) and the remaining set of variables while campbell2017within propose a variable selection method using the exclusive Lasso. Now, the Lasso estimates of the slope coefficients in a linear regression model solve the $\ell_1-$penalized least regression problem:

or, equivalently we have that

where $\boldsymbol{\beta} = \left( \beta_1,..., \beta_p \right)^{\prime}$, and $s$ and $\lambda$ are tuning parameters. The lasso uses a constraint in the form of $\ell_1-$norm: $\sum_{j=1}^p \left| \beta_j \right| \leq s$. Therefore, by using the $\ell_1-$penalty, the Lasso achieves variable selection and shrinkage simulatneously, while $\lambda$ controls the amount of regularization.

Denote with $\boldsymbol{\beta}^0$ the true vector of parameters. There are three types of errors of interest in LASSO regression, that is, (i) the prediction error: $ \left\lVert \boldsymbol{X} \left( \widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}^0 \right) \right\rVert_2^2$, (ii) the parameter estimation error: $ \left\lVert \left( \widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}^0 \right) \right\rVert_2^2$, and (iii) the model selection error: $\mathbb{P} \left( \mathsf{supp} \left( \widehat{\boldsymbol{\beta}} \right) = \left( \boldsymbol{\beta}^0 \right) \right)$. In other words, the prediction problem implies that given a random sample $\left( y_i, x_i \right), i = 1,...,n$ and the covariates $x_{n+1}$ of an additional observation, we aim to predict the unobserved outcome $y_{n+1}$. Choice of $\lambda$ is achieved using Cross-Validation method:

where $\widehat{\beta}^{(-i)}(\lambda)$ is the estimator $\widehat{\beta}^{\ell_1 } (\lambda)$ obtained after dropping observation $i$ from the sample. Notice that cross-validation is usually used to choose the tuning parameter $\lambda$ for in-sample model estimation and variable selection. However, there are theoretical results that show that LASSO performs well for prediction even for large number of covariates $p$, as long as the true parameter $\beta$ is sparse. Under appropriate regularity conditions with high probability (for large $n$) we have that

The cross-validation is usually used to choose the tuning parameter $\lambda$ for in-sample model estimation and variable selection. Furthermore, a key property for model selection implies that under appropriate regularity conditions with high probability (large $n$) one can show the following relation

Under our assumptions, $\frac{ \boldsymbol{U}^{\prime} \boldsymbol{X}_j }{T}$, $j = 1,...,p$, is Gaussian with mean zero and variance bounded by $\frac{ \mathcal{C} \sigma^2 }{n}$. Then, it holds that

The aim is to make (bound) the above probability as small as possible. Moreover, to obtain convergence rates and bounds we can use the following definitions $z = \sigma \sqrt{ \frac{2 \mathsf{log} \left( \frac{ep}{\delta} \right) }{n} }$, where $\delta > 0$ and $e := \mathsf{exp}(1)$. Set also $\mathcal{C} = 1$. Therefore, we obtain that

Therefore, we get that

Another important aspect related to the statistical properties of the lasso shrinkage approach is the sign consistency. Specifically, Strong Sign Consistency: implies that one can use a pre-selected $\lambda$ to achieve consistent model selection via the LASSO. On the other hand, a General Sign Consistency: means that for a random realization there exists a correct amount of regularization that selects the true model. Further discussion and related asymptotic results on the consistency and regularization property\footnote{Professor Marcelo C. Medeiros gave a seminar with title: "Bridging Factor and Sparse Models" (see, fan2021bridging) at the Department of Economics, University of Southampton on the 6th of October 2021.} under different distributional assumptions can be found in medeiros2016.

Dantzig Selector

A large stream of literature has focused on the $\ell_1$ penalized LASSO estimator of parameters in high-dimensional linear regression when the number of variables can be much larger than the sample size. We consider the linear regression with many covariates such as

where $X$ is the $n \times M$ deterministic design matrix, with $M$ possibly much larger than $n$, and $\boldsymbol{ \epsilon }$ is a vector of $\textit{i.i.d}$ standard normal random variables. In particular, we are interested in the case of sparsity parameters, which implies that the high-dimensional vector $\beta$ has coefficients that are mostly $0$. Denote with $\widehat{S} \left( \beta \right)$, the residual sum of squares such that

for all $\beta \in \mathbb{R}^{M}$. Define the Lasso solution such that $\widehat{\beta}_L = \left( \widehat{\beta}_{1,L},...., \widehat{\beta}_{M,L} \right)$ by the following expression

where $r > 0$ is some tuning constant. Then, the corresponding Lasso estimator

A necessary and sufficient condition of the minimizer is that 0 belongs to the subdifferential of the convex function $\beta \mapsto n^{-1} \left| y - X \beta \right|_2^2 + 2 r \left| D^{1 / 2} \beta \right|_1$.

This implies that the Lasso selector $\widehat{ \beta }_{L}$ satisfies the constraint

Various studies in the literature present relevant applications and statistical theory of the Dantzig selector such as osborne2000lasso, candes2007dantzig, bickel2009simultaneous, koltchinskii2009dantzig and james2009dasso among others. We say that $\beta \in \mathbb{R}^M$ satisfies the Dantzig constraint if $\beta$ belongs to the set

The Dantzig estimator is defined by the following expression

where $\widehat{\beta}_D = \left( \widehat{\beta}_{1,D},..., \widehat{\beta}_{M,D} \right)$ is the Dantzig selector. By the definition of Dantzig selector, we have that $\left| \widehat{\beta}_{D} \right|_1 \leq \left| \widehat{\beta}_{L} \right|_1$. Notice that the Dantzig selector is computationally feasible, since it reduces to a linear programming problem. Finally, for any $n \geq 1$, $M \geq 2$, we consider the Gram matrix as below

and let $\phi_{\text{max}}$ denote the maximal eigenvalue of $\psi_n$.

Oracle inequalities for prediction loss

In this section, we prove sparsity oracle inequalities for the prediction loss of the Lasso and Dantzig estimators. These inequalities allow us to bound the difference between the prediction errors of the estimators and the best sparse approximation of the regression function (e.g., by an oracle that knows the truth but is constrained by sparsity). We demonstrate that the distance between the prediction losses of the Dantzig and Lasso estimators is of the same order as the distances between them and their oracle approximations. Recall that an estimator has the oracle property if it is variable selection consistent and the limiting distribution of its subvector corresponding to the non-zero coefficients is the same as if their set were known prior to estimation (see, giurcanu2016thresholding).

Parameter Estimation and Model Selection Consistency Properties

Consider the linear regression model as below (see, amann2018uniform)

where $y \in \mathbb{R}^n$ is the response vector, $X \in \mathbb{R}^{ n \times p}$ the nonstochastic regressor matrix which is assumed to have full column rank, $\beta \in \mathbb{R}^p$ the unknown parameter vector, and $\epsilon \in \mathbb{R}^n$ the unobserved stochastic error term consisting of i.i.d distributed components with mean zero and finite second moments, defined on some probability space $\left( \Omega, \mathcal{F}, \mathbb{P} \right)$. Model consistency results can be found in shibata1986consistency and potscher1991effects. Moreover, to define the adaptive Lasso estimator, we consider the following expression

where $\hat{\beta}$ denotes the OLS estimator. In particular, we assume that the event $\left\{ \hat{\beta}_j = 0 \right\}$ to have zero probability, for all $j = 1,...,p$ and thus we do not consider this event occurring in the subsequent analysis.

Consistency Properties

Consider any sequence $\left( \beta_n \right)_{ n \in \mathbb{N} } \subset \mathbb{R}^p$ converging to $\phi$ and let

Then, by the Portmanteau Theorem, we have that

Then, for any $\omega \in \Omega$, we have that

In other words, while the limit of $\sqrt{ \frac{ n }{ \lambda^* } } \left( \hat{\beta} - \beta_n \right)$ will in general be random, the set $\mathcal{M}$ is not random. In particular, Proposition 8 shows that, for any $\omega$, the union of limits over all possible sequences of unknown parameters is always given by the same compact set $\mathcal{M}$.

Specifically, this observation is central for the construction of confidence regions in the following section. It also shows that while in general, a stochastic component will survive in the limit, it is always restricted to have bounded support that depends on the regressor matrix and the tuning parameter through the matrix $C$ and the quantities $\psi$ and $\lambda^0$. Interestingly, $\mathcal{M}$ only depends on $\psi$ for the components where $\psi_j = \infty$, in which case the set $\mathcal{M}$ loses a dimension. In particular, this can be seen as a result of the $j-$th component being penalized much less than the maximal one, so that the scaling factor used in Theorem 7 is not large enough for this component to survive in the limit.

Consider the following expression (see, amann2018uniform)

where $\hat{\beta}_{AL} \neq 0$. Notice that the left-hand side is bounded by $L$, whereas the right-hand side converges to $\frac{ c }{ | \beta_j | }$ in probability. We therefore get $\mathbb{P}_{ \beta } \left( \hat{\beta}_{AL,j} = 0 \right) \to 1$, for all $\beta_j \in \mathbb{R}$, satisfying the condition $\left| \beta_j \right| < \frac{ c }{ L }$. The set $\mathcal{M} = \mathcal{M}_1$ acts as a benchmark for confidence sets in the sense that if we take a "slightly larger" set, multiplied with the appropriate factor and centered at the adaptive Lasso estimator, we get a confidence region with minimal asymptotic coverage probability equal to 1. But if we base the region on a "slightly larger" set than $\mathcal{M}$, we end up with a confidence set of asymptotic minimal coverage 0.

A Lasso-based Time Series Regression Model

We denote with $\theta = \left( \rho^{*}, \beta ^{* \prime} \right)^{\prime}$ the set of model parameters. The set of active variables is denoted with

Moreover, denote with $\mathcal{S} = \text{diag} \left( T, \sqrt{T}, ..., \sqrt{T} \right)$ denotes a $\left( p + 1 \times p + 1 \right)$. Denote with $\hat{ \theta } = \left( \hat{\rho}, \hat{\beta}^{\prime} \right)^{\prime}$ denote the minimizer of the objective function. Furthermore, denoting $\mathcal{M}_0$ the true model and $\hat{\mathcal{M} }$ the estimated model, then we can say that the shrinkage methodology is consistent if for all $\left( \rho^{*}, \hat{\beta}^{*} \right)$ it holds that $\mathbb{P} \left( \hat{\mathcal{M}} = \mathcal{M}_0 \right) \to 1$. Thus, the shrinkage selection methodology is said to be conservative if for all $\left( \rho^{*}, \hat{\beta}^{*} \right)$ it holds that $\mathbb{P} \left( \hat{ \mathcal{M} }_0 \subset \hat{ \mathcal{M} } \right) \to 0$, which implies that the probability of excluding relevant variables tend to zero.

Oracle Property

We discuss the oracle property of the adaptive Lasso for stationary and nonstationary autoregressions.

Lasso shrinkage with long memory regression errors

In many problems of practical interest in which the lasso shrinkage is applied, such as for example when modelling and forecasting Realized Volatility measures, it is reasonable to consider the existence of long memory errors. We investigate the asymptotic behaviour of Lasso in regression models with long memory errors. Consider $X_i = \left( x_{i1},...., x_{ip} \right)^{\prime}$ for $i = 1,...,n$ to be the vector of design matrices and $Y_i$ to denote its response variable. Thus we have the model (see, kaul2014lasso)

The errors $\epsilon_i$ are assumed to be long memory moving average with i.i.d innovations, that is,

where $a_k = c_0 k^{-1 + d}, \forall k \geq 1$, $0 \leq d \leq \frac{1}{2}$ and some constant $c_0 > 0$, and $a_k = 0$ for $k \leq 0$. Moreover, we have that $\zeta_j, j \in \mathbb{Z} := \left\{ 0, \pm 1, \pm 2,... \right\}$, are i.i.d RV's with mean zero and variance $\sigma_{ \zeta }^2$. Without loss of generality we assume that $\sigma_{ \zeta }^2 = 1$. We denote with $X = \left( x_{ij} \right)_{ n \times p}$ as the design matrix, and $\epsilon := \left( \epsilon_1,..., \epsilon_n \right)^{\prime}$. Moreover, $\left\{ \epsilon_i, i \in \mathbb{Z} \right\}$ is a stationary process with autocovariance function

Moreover, the Lasso estimate of $\beta$ is defined as follows

where $Y = \left( Y_1,..., Y_n \right)^{\prime}$ and $\left\lVert \beta \right\rVert_1 := \sum_{j=1}^p \left| \beta_j \right|$ denotes $\ell_1$ norm of $\beta = \left( \beta_1,..., \beta_p \right)^{\prime}$.

Notice that the literature in the area of regularized estimation with dependence considerations is scarce. In this paper, we we investigate the asymptotic behaviour of Lasso under strong dependence structure and less restrictive model assumptions. In particular, we assign a long memory structure on the model errors $\epsilon$, that is, $\sum_{ k = 1}^{ \infty } \left| \gamma_{\epsilon} (k) \right| = \infty$. We provide restrictions on the rate of increase of the design variables as well as the rate of increase of the dimension of $p$ in order to obtain the corresponding finite sample error bounds. Furthermore, we allow the design variables to grow with the restriction $\sum_{1 \leq i \leq n} x_{ij}^2 = \mathcal{O}(n)$, and hence the results obtained can also easily be extended to the case of Gaussian random designs.

Notice that the framework proposed by babii2022machine considers a machine learning application for high-frequency time series panel data. The particular time series filters are constructed based on mixed-frequency data and therefore assumptions on the dependence structure of the errors such as persistence, mixing and long-memory are indeed plausible in such economic and finance studies.

Results with finite sample

In this Section we prove a finite sample oracle inequality for the Lasso solution when the design matrix is non-random. We define with

where we have that

Moreover, we denote with

Therefore, we shall prove that with an appropriate choice of $\lambda_n$, the Lasso solution obeys the following oracle inequality in the long memory case, for any $n \geq 1$,

where $\lambda_n = \left( \mathcal{O}(1) \right) \text{log} (p) / n^{1/2 - d}$, under some conditions on the design matrix. Moreover, $s_0$ denotes the cardinality of the set of non-zero components of $\beta$ and $\phi_0$ is a constant depending on the design matrix $X$. Therefore, in order to prove the result we need to obtain a probability bound for the set as below

for a proper choice of $\lambda_{0n}$. Thus, once this probability bound is obtained, the oracle inequality follows by deterministic arguments.

Notice that the notion of weak dependence for stationary time series is measured in terms of covariance functions. The following lemma given by gupta2012note is useful.

For the proof see gupta2012note. The argument is to let $\omega \in A$, then it follows that

Residual empirical process based on the ALasso

Moreover, the framework proposed by chatterjee2015residual provides regularity conditions and assumptions for deriving the asymptotic behaviour the adaptive Lasso estimator. Specifically, the Alasso estimator of $\beta$ is defined as the minimizer of the weighted $\ell_1-$penalized least squares criterion function

where, $\lambda_n > 0$ is a regularization parameter, $\gamma > 0$ and$\tilde{\beta}_{j,n}$ is the $j$th component of $\widetilde{ \boldsymbol{\beta} }_n$, a consistent preliminary estimator of $\boldsymbol{\beta}$.

Main results

\paragraph{Asymptotic uniform linearity in high dimensions}

For $t \in [0,1]$, we define

We set $J(t) = f \left( F^{-1} (t) \right)$. The first result of this section proves the AUL property of the ALASSO based residual empirical distribution function for p.

Notice that Theorem (ref) shows that the AUL property holds for the empirical distribution function of the residuals based on the ALASSO fit even for $p >>n$, provided the regularity conditions hold and the regression model has enough sparsity.

\paragraph{Functional oracle property of the ALASSO}

The AUL property has various important applications in the context of statistical inference in high dimensional regression. For instance, it allows to establish the asymptotic distribution of the residual edf which, in turn, can be used to carry the goodness of fit tests on $F$ and set confidence bands for $F$. Under standard weak convergence theory we have that

where $\widehat{F}_n^{OR}$ is the edf of the residuals based in OLS under the oracle property.

Moreover, $\mathcal{Z}_{\infty}$ is a zero mean Gaussian process with covariance function as below

We also define the ALASSO based residual empirical process as below

Sparsity in high dimensional estimation problems

The Lasso estimator appears to be suitable for high dimensional inference problems such as the case when $p$ is much larger than the time series observations $n$. Some useful results are as below:

\paragraph{Regression in the iid case} Under the usual assumption that the $\epsilon_i$ are iid and subGaussian,

for some known $\sigma^2$, then we have that

Lasso shrinkage for Predictive Regression

In this model, the unknown coefficients $\beta_n^{*}$ can be obtained from the data by running OLS

The asymptotic behaviour of the OLS estimator has been extensively examined in the time series econometrics literature (see, koo2020high, lin2020robust, yousuf2021boosting, lee2022lasso). A relevant framework for lasso time series regression is proposed by adamek2023lasso.

In particular, the following functional central limit theorem holds

To represent the asymptotic distribution of the OLS estimator, define $u_i^{+} = u_i - \Sigma_{eu}^{\prime} \Sigma_{ee} e_i^{\prime}$. By definition, Cov$\left( e_{ij}, u_i^{+} \right) = 0$ for all $j$ so that

which is the sum of a mixed normal random vector and a non-standard random vector. The OLS limit distribution is then given by

where $\Omega: = \displaystyle \int_0^1 B_e(r) B_e (r)^{\prime} dr$. This implies that when we inflate $\widehat{ \beta }_j^{OLS}$ by the factor $n^{ \delta_j }$ so that its magnitude is comparable to the constant $\beta_j^{0*}$, we attain consistency in that

Moreover, if $\beta_j^{0*} \neq 0$, when $\delta_j$ is close to 1, the signal of $x_j$ is weak, and therefore the rate of convergence is slow. However, notice that some true coefficients $\beta_j^{*}$ could be exactly zero, where the associated predictors would be redundant (inactive) in the regression. Let $M^{*} = \left\{ j : \beta_j^{0*} \neq 0 \right\}$ be the index set of regressors relevant to the regression, $p^{*} = \left| M^{*} \right|$, and $M^{* c} = \left\{ 1,..., p \right\} M^{*}$ be the set of redundant regressors. For simplicity, we refer to $M^{*}$ as the active set, meaning that it plays an active role in the regression, and we call $M^{* c}$ the inactive set. Define $\widehat{\beta}^{ora} = \left( \widehat{\beta}^{ \text{ora} \prime}_{ M^{*} }, \widehat{\beta}^{ \text{ora} \prime}_{ M^{*}} \right)$, where

and the complement of the oracle estimator is just the remaining coefficients which are identical equal, that is, $\widehat{\beta}^{ \text{ora} }_{ M^{*} c } = 0$. Notice that the oracle estimator is the infeasible oracle information of $M^{*}$. The asymptotic distribution of the estimator is given by

where $\Omega_{ M^{*} }$ is the $p^{*} \times p^{*}$ submatrix $\left( \Omega_{j,j^{\prime} } \right)_{ j,j^{\prime} \in M^{*} }$ and $\zeta_{ M^{*} }$ is the $p^{*} \times 1$ subvector $\left( \zeta_j \right)_{ j \in M^{*} }$.

Lasso estimation with a structural break

Relevant studies with frameworks for lasso estimation robust to the presence of structural breaks in the covariates include chan2014group, cheng2016shrinkage, qian2016shrinkage and smith2019variable. For any $n-$dimensional vector $W = \left( W_1,..., W_n \right)^{\prime}$, define the empirical norm as $\left\lVert W \right\rVert_n := \left( \frac{1}{n} \sum_{i=1}^n W_i^2 \right)^{1/2}$.

Then, we define the prediction risk as below

Define with $\alpha_0 = \left( \beta_0^{\prime}, \delta_0^{\prime} \right)^{\prime}$. Then, we can rewrite the model as below $Y_i = \boldsymbol{X}_i ( \uptau_0 )^{\prime} \alpha_0 + U_i$, for $i = 1,...,n$. Let $\boldsymbol{y} \equiv \left( Y_1,..., Y_n \right)^{\prime}$. Then, for any fixed $\uptau \in \mathcal{T}$ we consider the residual sum of squares given by

Classical Shrinkage Type Estimation Approach

In this section, we present the Score function approach proposed in the framework of sen1987preliminary, although the particular methodology doesn't correspond to a high-dimensional environment, it has many possible applications such as a "Shrinkage type Structural Break Estimation in High-Dimensional Predictive Regressions" or "Testing for Predictability in High Dimensional Environments via a Shrinkage type Structural Break Estimator".

Consider the score function given by $\psi = \left\{ \psi(x), x \in \mathbb{R} \right\}$ required for the definition of M-estimators. Moreover, suppose that $\psi (x) = \big\{ \psi_1 (x) + \psi_2(x) \big\}$, where $\psi_1(.)$ and $\psi_2(.)$ are both non-decreasing and skew-symmetric such that, $\psi_j (x) + \psi_j (-x) = 0, \forall \ x \in \mathbb{R}, j = 1,2$. Denote with $\boldsymbol{A}_T = \left( \boldsymbol{a}_1,.., \boldsymbol{a}_T \right)$ and for every $\boldsymbol{b} \in \mathbb{R}^p$ we define the following $p-$dimensional functional with $T$ elements in each dimension

Then, under the assumption of stationarity then, the RME for $\hat{\boldsymbol{\beta}}_{1T}$ of the model parameter $\boldsymbol{\beta}_{1}$ is a solution to the following equation

We introduce a suitable $M-$test statistic for testing the null hypothesis $\mathbb{H}_0: \boldsymbol{\beta}_2 = \boldsymbol{0}$,

where $\hat{\boldsymbol{\beta}}_{1T}$ is the restrained M-estimator (RME) of $\boldsymbol{\beta}_1$. Notice that under stationarity, it can be proved that the restrained M-estimator (RME) generally performs better than the unrestrained M-estimator (URE) especially when $\boldsymbol{\beta}_2$ is close to $\boldsymbol{0}$ (see, sen1987preliminary). Furthermore, denote with

Therefore, an appropriate aligned $M-$test is given by

Under the null hypothesis, $\mathbb{H}_0$, $\mathcal{L}_T$ has asymptotically the chi-square distribution with $p_2$ degrees of freedom. Therefore, based on the significance level $\alpha$ such that $0 < \alpha < 1$, the preliminary test under the null, $\mathbb{H}_0$, is rejected based on the critical values.

Considering for example one is interested to the conditional mean specification of a predictive regression model, focusing on the OLS against the IVX estimation approaches based on the full set of regressors in the model. The, the corresponding expressions are:

and

In other words, the authors of the particular framework focus on preliminary test M-estimation (PTME) formulation, thus the $\hat{\boldsymbol{\beta}}_{1n}^{PT}$ is chosen as the RME or UME, according as this preliminary test leads to the acceptance or rejection of the null hypothesis $\mathbb{H}_0$. Moreover, the shrinkage $M-$estimator (SME), based on the usual James-Stein rule incorporates the same test statistic in a smoother manner. When $\boldsymbol{\beta}_2$ is very close to $\boldsymbol{0}$, generally both the PTME and SME perform better than the UME, but the RME may still be better than either of them. On the other hand, for $\boldsymbol{\beta}_2$ away from $\boldsymbol{0}$, the RME may perform rather poorly, while both the PTME and SME are robust. This relative picture on the performance characteristics of all four versions of M-estimators can be best studies within an asymptotic framework. One can consider the notion of asymptotic distributional risk (ADR) as well as the asymptotic risk efficiency (ARE) for the various versions of the M-estimators (see, sen1987preliminary).

Limiting Distributional Risk

In the multivariate location model sen1985some have pointed out that shrinkage estimation works out well only in a shrinking neighbourhood of the pivot. Since for $\boldsymbol{\beta}_2$ the pivot is taken as $\boldsymbol{0}$, we consider a shrinking neighbourhood of $\boldsymbol{0}$ and we consider a sequence $\left\{ \mathcal{K}_T \right\}$ of alternatives,

which implies that the null hypothesis $\mathbb{H}_0$ reduces to $\mathbb{H}_0: \boldsymbol{\xi} = \boldsymbol{0}$. Using a suitable estimator $\boldsymbol{\beta}^{\star}_{1T}$ of $\boldsymbol{\beta}_{1}$, we denote by $\boldsymbol{G}^{\star} ( \boldsymbol{\mathsf{u}} ) = \underset{ T \to \infty }{ \mathsf{lim} } \ \mathbb{P} \bigg( T^{- 1/ 2} \big( \boldsymbol{\beta}^{\star}_{1T} - \boldsymbol{\beta}_{1} \big) \leq \boldsymbol{\mathsf{u}} \big| \mathcal{K}_T \bigg), \ \ \ \boldsymbol{\mathsf{u}} \in \mathbb{R}^{p_1}$, where we assume that $\boldsymbol{\mathsf{u}}$ is nondegenerate. Then, with a quadratic error loss such that $T \big( \boldsymbol{\beta}^{\star}_{1T} - \boldsymbol{\beta}_{1} \big)^{\prime} \boldsymbol{W} \big( \boldsymbol{\beta}^{\star}_{1T} - \boldsymbol{\beta}_{1} \big)$, for a suitable probability distribution function $\boldsymbol{W}$, the ADR of $\boldsymbol{\beta}^{\star}_{1T}$ is defined as below

where $\boldsymbol{V}^{\star}$ is the dispersion matrix for the asymptotic distribution of $\boldsymbol{G}^{\star} ( . )$.

Asymptotic representation of estimators

We follow the framework presented in mukherjee1998preliminary.

\paragraph{Notation} The $i-$th row and the $j-$th column of any matrix $\boldsymbol{S}$ are denoted by $s_{i \bullet}$ and $s_{ \bullet j}$, respectively. Then if $\boldsymbol{S}$ has $p$ columns, $\boldsymbol{D}_s$ denotes the maximal diagonal $\big( \left\lVerts_{ \bullet 1} \right\rVert, ..., \left\lVerts_{ \bullet p} \right\rVert \big)$, where $\left\lVerts_{ \bullet j} \right\rVert$ is the usual Euclidean norm of the $j-$th column of $\boldsymbol{S}$. Moreover, the usual partitions, $\boldsymbol{S}_{ii.j}$ stands for the matrix $\boldsymbol{S}_{ii} - \boldsymbol{S}_{ij} \boldsymbol{S}_{jj} \boldsymbol{S}_{ji}$. Then, we define the following matrices which appear as the scaling factors in the asymptotic distributions of the different estimators. We denote with $\boldsymbol{A}_x = n^{-d} \boldsymbol{D}_x$, $\boldsymbol{B}_x = n^{d} \boldsymbol{D}_x$, $\boldsymbol{A}_{x1} = n^{-d} \boldsymbol{D}_{x1}$. Moreover, we have that $\underset{ n \to \infty }{ \mathsf{lim} } n^{-1} \boldsymbol{X}^{\prime} \boldsymbol{X} \ \ \text{exists and equals} \ \boldsymbol{C}$, with $\boldsymbol{C}$ being positive definite. Consequently,

Denote with $\boldsymbol{\Sigma}_n$ to denote the dispersion matrix of $\boldsymbol{\epsilon}$. Then, the dispersion of $\boldsymbol{B}_x^{-1} \boldsymbol{X}^{\prime} \boldsymbol{\Sigma}_n \boldsymbol{X} \boldsymbol{B}_x^{-1}$ and it follows that the limit of $\boldsymbol{E}_n$ exists and equals $\boldsymbol{E}$, say. As in sen1987preliminary, the shrinkage least-squares estimator (SLSE) $\hat{\beta}_s$ can be motivated as a smoothed version of a preliminary test estimator. Thus, for each $c > 0$, this is defined as below

For large values of $\mathcal{L}_T$, $\hat{\beta}_s$ yields to $\tilde{\beta}_1$, whereas it performs quite differently from $\hat{\beta}_{pt}$ form small values of $\mathcal{L}_T$. Define the dispersion matrix

Then, it can be proved that the second term converges to $\boldsymbol{R}_{22.1} \boldsymbol{D}_{c2} \boldsymbol{\eta}$. Therefore, we define with

implying that $\boldsymbol{H}_n$ converging to $\boldsymbol{H}$ and $\boldsymbol{\xi} := \boldsymbol{D}_{c2} \boldsymbol{\eta} \in \mathbb{R}^{p_2}$. Then, $\mathcal{L}_n$ is defined as

which converges to a noncentral chi-square random variable with $p_2$ degrees of freedom and noncentrality parameter $\boldsymbol{\delta} := \boldsymbol{\xi}^{\prime} \boldsymbol{L}_{22}^{-1} \boldsymbol{\xi}$. Therefore, under the null hypothesis $\boldsymbol{\eta} = \boldsymbol{0}$ or equivalently, $\boldsymbol{\xi} = \boldsymbol{0}$, the test statistic $\mathcal{L}_n$ converges to a chi-square random variable with $p_2$ degrees of freedom. Define the projection matrix $\boldsymbol{P}_{x1} = \boldsymbol{X}_1 \left( \boldsymbol{X}_1^{\prime} \boldsymbol{X}_1 \right) \boldsymbol{X}_1^{\prime}$. Let $\boldsymbol{D}_n := \boldsymbol{R}_n^{-1}$ and $\boldsymbol{D} := \boldsymbol{R}^{-1}$. Then, the covariance matrices of the random vectors $\boldsymbol{D}_n \boldsymbol{U}_n$ and $\boldsymbol{D} \boldsymbol{U}$ are denoted by $\boldsymbol{L}_n$ and $\boldsymbol{L}$ respectively.

Lasso Inference for High Dimensional Time Series under NED

Following the framework proposed by adamek2023lasso, consider the linear time series regression

where $\boldsymbol{x}_t = \left( x_{1,t},..., x_{N,t} \right)^{\prime}$ is an $(N \times 1)$ vector of explanatory variables and we assume that we have a high-dimensional time series model where $N$ can be larger than $T$.

The framework proposed by adamek2023lasso is the first to utilize the NED framework for establishing uniformly valid high-dimensional inference. For instance, wong2020lasso consider time series models with $\beta-$mixing errors, which has the advantage of allowing for general forms of dynamic misspecification resulting in serially correlated error terms, but, rules out several relevant data generating processes, and is in addition typically difficult to verify. On the other hand, masini2022regularized use an m.d.s assumption on the innovations in combination with sub-Weibull tails and a mixing assumption on the conditional covariance matrix. However, the particular m.d.s assumption does not allow for dynamic misspecification of the full model. Importantly, the NED assumption on $u_t$ does allow for misspecified models as well, in which case we view $\boldsymbol{\beta}^0$ as the coefficients of the pseudo-true model when restricting the class of models to those linear in $\boldsymbol{x}_t$.

Therefore, this econometric intuition allows to view the high-dimensional time series model of adamek2023lasso as simply the linear projection of $y_t$ on all the variables in $\boldsymbol{x}_t$, with $\boldsymbol{\beta}^0$ in that case representing the best linear projection coefficients, which implies that $\mathbb{E} [ u_t ] = 0$ and $\mathbb{E} [ u_t x_{j,t} ] = 0$. On the other hand, $u_t$ is not likely to be an $\textit{m.d.s}$. Thus, allowing for misspecified dynamics is crucial for developing the theory for the nodewise regressions underlying the desparsified lasso. The NED-order $m$ and sequence of size $-d$ play a key role in the development of the limit results in which these parameters enter as the corresponding asymptotic rates of related quantities. Obviously, there is a trade-off between the thickness of the tails allowed for and the moment of dependence - which is measures via the mixing rate.

Notice that the above examples, provide cases in which although the classical martingale difference sequence framework fails to apply, the more general NED framework extends the applicability of certain time series processes.

Next, we focus on the related assumptions/conditions on the sparsity properties of the parameter vector. Specifically, the key condition to apply the lasso successfully is that the parameter vector $\boldsymbol{\beta}_0$ is (at least approximately) sparse. This is given by the following assumption.

Inference on low-dimensional parameters

In this section we establish the uniform asymptotic normality of the desparsified lasso focusing on low-dimensional parameters of interest (see, adamek2023lasso). Specifically, we consider testing $P$ joint hypotheses of the form $\boldsymbol{R}_N \boldsymbol{\beta}^0 = \boldsymbol{q}$, via a Wald statistic, where $\boldsymbol{R}_N$ is an appropriate $P \times N$ matrix whose non-zero columns are indexed by the set

with cardinality $h:= |H|$. Notice that we can allow for $h$ to increase in $N$ (and therefore $T$). Focusing on inference on finite set of parameters such that we can apply a standard central limit theorem. Therefore, given our time series setting, the long-run covariance matrix is

where $\boldsymbol{w}_t = \big( v_{1,t} u_t,..., v_{N,t} u_t \big)^{\prime}$, enters the asymptotic distribution. Then, $\boldsymbol{\Omega}_{N,T}$ can equivalently be written

In order to estimate the asymptotic variance $\boldsymbol{\Psi}$, we suggest to estimate $\boldsymbol{\Omega}_{N,T}$ with the long-run variance kernel estimator such that

with $\hat{w}_{j,t} = \hat{v}_{j,t} \hat{u}_t$, the kernel $K(.)$ can be taken as the Bartlett Kernel $K ( \frac{\ell}{ Q_T} ) = \left( 1 - \frac{\ell}{Q_T} \right)$.

Error Bound and the Consistency for the Lasso

In this section, we consider a new error bound for the lasso in a high-dimensional time series model. Thus, the Lasso estimator of the parameter vector $\boldsymbol{\beta}^0$ in the model (see, adamek2023lasso)

where $\boldsymbol{y} = \left( y_1,..., y_T \right)^{\prime}$ is the $T \times 1$ response vector, $\boldsymbol{X} = \big( \boldsymbol{x}_1,...., \boldsymbol{x}_T \big)^{\prime}$ is the $( T \times N )$ design matrix and $\lambda > 0$ a tuning parameter. In other words, the above optimization problem adds a penalty term to the least squares objective to penalize parameters that are different from zero. Thus, when deriving this error bound, we typically require that $\lambda$ is chosen sufficiently large to exceed the empirical process such that $\underset{ j }{ \mathsf{max} } \ \left| \frac{1}{T} \sum_{t=1}^T x_{j,t} u_t \right|$, with high probability. Moreover, we define the following set such that

and establish the conditions under which $\mathbb{P} \big( \mathcal{E}_T \left( \frac{T \lambda }{4} \right) \big) \to 1$.

Uniformly Valid Inference via the Disparsified Lasso

Following the framework of adamek2023lasso, we consider an application of the disparsified Lasso.

High Dimensional Feature Selection Methods

High dimensional statistical problems address the challenge of conducting robust estimation and inference when the number of unknown parameters, $p$, is much larger than the number of observations, $n$. The literature on high-dimensional feature selection methods include the conditional sure independence screening\footnote{Note that the covariate screening approach is widely used in the applied statistics literature, see aschard2017covariate.} approach as discussed by fan2008sure, wang2009forward, ke2014covariance and barut2016conditional as well as he2013quantile and kong2019screening in the case of a conditional quantile functional form is employed. The covariate screening approach as a tool for variable selection works well in ultra high-dimensional linear regression models with stationary covariates, without requiring distributional restrictions. Statistical estimation methodologies in high-dimensional settings include belloni2011, belloni2014inference, wang2012quantile, chernozhukov2015post, fan2017estimation as well as sun2015sprem and wang2016high for the cases of projection screening.

A fundamental research question remains in most cases when considering modelling a high-dimensional vector of regressors is determining the statistical significance of either individual covariates or the model selection mechanism. lockhart2014significance propose a test for detemining statistical significance with the Lasso, while chudik2018one propose a methodology for selecting covariates in a high-dimensional environment using the multiple testing approach. Furthermore, battey2018distributed considers hypothesis testing and parameter estimation in the context of the divide-and-conquer algorithm. In particular, in a unified likelihood-based framework, we propose new test statistics and point estimators obtained by aggregating various statistics from $k$ subsamples of size $n / k$, where $n$ is the sample size.

Specifically, in low dimensional and sparse high dimensional settings, battey2018distributed examine how large $k$ can be, as $n$ grows large, such as the loss of efficiency due to the divide-and-conquer algorithm is negligible. In other words, the resulting estimators have the same inferential efficiencies and estimation rates as an oracle with access to the full sample. While hypothesis testing in a low dimensional context is straightforward, in the sparse high dimensional setting, nuisance parameters introduce a nonnegligible bias, causing classical low dimensional theory to break down. Therefore, in their high dimensional Wald construction, the phenomenon is remedied through debiasing of the estimator, which gives rise to a test statistic with tractable limiting distribution. Thus, they find that the theoretical upper bound on the number of subsamples guaranteeing the same inferential or estimation efficiency as the whole-sample procedure is $k = o( ( s \mathsf{log} d )^{-1} \sqrt{n} )$ in the linear model, where $s$ is the sparsity of the parameter vector.

On the other hand, in time series regression models incorporating high dimensionality features requires to develop further tools such as Gaussian approximations (e.g., see chernozhukov2014gaussian and zhang2017gaussian) and consistent model selection techniques robust to the presence of heteroscedasticity (see, halunga2017heteroskedasticity) and nonstationarity. In particular in the time series econometrics literature many open problems remain such as robust model selection methodologies for ultra-high dimensional environments with heterogenous data. We return back to Vector Autoregression Processes, although not necessarily a high-dimensional VAR$(p)$ process. In practise, we consider the case where the dimension of the time-series vector is less than the sample size $n$. Formally, a $p-$dimensional vector-valued stationary time series $X_t = ( X_{1t},..., X_{pt} )$, with $t=1,...,n$, can be modelled using a VAR representation of lag $d$ with serially uncorrelated Gaussian errors, which takes the following form

where $A_1,..., A_d$ are (p x p) matrices and et is a p-dimensional vector of possibly correlated innovation shocks. Therefore, the main objective in VAR models is to estimate the transition matrices $A_1,..., A_d$, together with the order of the model $d$, based on time series realizations $( X_0, X_1,..., X_n )$. Then, the structure of the transition matrices provides insights into the complex temporal relationships amongst the p time series and the particular representation provides a way to apply forecasting techniques (see, michailidis2013autoregressive, basu2015regularized and basu2021graphical).

Ultra-high dimensionality under dependence

According to yousuf2018variable, it is possible to achieve the sure screening property in the ultrahigh dimensional setting with dependent errors and covariates. However, in order to do that we need to make stronger assumptions on the moments of both the error and covariate processes. Specifically, if the error and covariate processes are assumed to follow a stronger moment condition, such as, $\Delta_{0,q} (\varepsilon) < \infty$ and $\Phi_{0,q} ( \boldsymbol{x} ) < \infty$ for arbitrary $q > 0$, we can then achieve a much larger range of $p_n$ which will cover the ultra-high dimensional case. Specifically, we need a condition that implies that the tails of the covariate and error processes are exponentially tight. A wide range of processes satisfy the above condition.

Suppose that $\varepsilon_i$ is a linear process such that

with $\sum_{\ell = 0}^{\infty} | f_{\ell} | < \infty$ then it holds that

Thus, if we assume that $e_0$ is sub-Gaussian, then $\tilde{\alpha}_{\epsilon} = \frac{1}{2}$ since $\left\lVert e_0 \right\rVert_{ q } = \mathcal{O}_p ( \sqrt{q} )$ then, when $e_i$ is sub-exponential then it holds that $\tilde{\alpha}_{\epsilon} = 1$.

More generally, for $e_i := \sum_{j=0}^{\infty} f_j e_{i-j}^p$ and thus if $e_i$ is sub-exponential, we have that $\tilde{\alpha}_{\epsilon} = p$, where $p$ is the number of components of the high-dimensional vector. The high dimensionality as well as the non-parametric aspect can slow down the rate of convergence and thus especially the variance of the estimator can have slower convergence to the corresponding variance of the process.

Thus, for any fixed $q$, we are not placing additional assumptions on the temporal decay rate of the covariate error processes besides requiring $\Delta_{0,q} (\varepsilon), \Phi_{0,q} (\varepsilon) < \infty$. On the other hand, the ultrahigh dimensional setting requires geometrically decaying strong mixing coefficients, in addition to requiring sub-exponential tails for the response. Consider for example the case that $\varepsilon_i = \sum_{j=0}^{\infty} f_j e_{i-j}$, geometrically decaying strong mixing coefficients would require the coefficients, $f_j$, to decay geometrically.

SIS with dependent observations

Sure Independence Screening, is a method of variable screening based on ranking the magnitudes of the $p_n$ marginal regression estimates. Under appropriate conditions, this simple procedure is shown to possesses the sure screening property. Define with

Therefore, $\hat{\rho}_j$ is the OLS estimate of the linear projection of $Y_t$ onto $X_{tj}$. Let

and let $| \mathcal{M}_{*} | = s_n << n$ be the size of the true sparse model. Moreover, we sort the elements of $\hat{\boldsymbol{\rho}}$ by their magnitudes. Thus, for any given $\gamma_n$, define a sub-model such that

and let $\left| \hat{\mathcal{M}}_{\gamma_n} \right| = d_n$ be the size of the selected model. Furthermore, notice that the screening property states that for an appropriate choice of $\gamma_n$, we have that

Experimental design examples

\paragraph{Experimental Design I: Uncorrelated Features}

Consider the model for the covariate process we have

where $A_1 = \mathsf{diag} \left\{ \gamma \right\}$ and $\gamma \sim \mathsf{Unif} [ 0.4, 0.6 ]$. Moreover, we set with $\boldsymbol{\eta}_t \sim \mathcal{N} (\boldsymbol{0},\boldsymbol{\Sigma}_{\eta})$. In particular, when we are dealing with uncorrelated predictors we set $\Sigma_{\eta} = I_{ p_n }$ which represents the no correlation case.

\paragraph{Experimental Design II: Correlated Features}

One can compare the performance of SIS and GLSS in the case of correlated predictors. The covariate process is generated with $A_1 = \mathsf{diag} \left\{ 0.4^{ | i - j | + 1 } \right\}_{ i,j \leq p_n }$ such that $\boldsymbol{\eta}_t \sim \mathcal{N} ( 0, \Sigma_{\eta} )$ and $\boldsymbol{\eta}_t \sim t_5 ( 0, \boldsymbol{V} )$ , with

Uniform Asymptotic Inference and Model Selection

This is another important topic which we cannot cover here extensively (see, tibshirani2018uniform).

Uniform-in-Submodel Bounds

Following the framework proposed by kuchibhotla2021uniform, suppose that $\mathcal{M} = \left\{ M_1,...,M_L \right\}$ denotes a collection of submodels, where $M_j$ represents a subset of covariates for $1 \leq j \leq L$. Also, denote with $\widehat{\beta}_{M_j}$ to represent the least-squares estimator for the linear regression of the response on the covariates of the model $M_j$. Then, by simultaneous consistency, we mean the existence of target vectors $\left\{ \beta_{M_j} : 1 \leq j \leq L \right\}$ such that the following result holds

More specifically, if $\widehat{M}$ is a selected model, then one can perform inference of $\beta_{ \widehat{M} }$ by estimating the distribution of $\widehat{ \beta }_{ \widehat{M} }$. According to kuchibhotla2021uniform, in practise, even though various model-selection criteria like $C_p$ (AIC), (BIC) and lasso have been recommended for covariate selection in linear regression, developing more general asymptotic results to the asymptotic uniform linear representation in the special case of the least-squares linear regression estimator are still an ongoing research field. The framework of variable selection and linear regression is often considered in the context of high-dimensional linear regressions. Although, the procedure uses only a reduced set of variables in the final regression, it uses all variables in the proceding covariate selection stage of estimation.

Suppose that $\widehat{M} \in \mathcal{M}$ is the final selected submodel, where $\mathcal{M}$ is some finite and countable collection of models, with $\widehat{ \beta }_{ \widehat{M} }$ the corresponding least-squares estimator. The estimator $\widehat{ \beta }_{ \widehat{M} }$ is known as the post-regularization estimator in the high-dimensional statistics literature when $\widehat{M}$ is obtained from a regularized least-squares procedure. A relevant question regarding the asymptotic behaviour of the particular estimator is "What does $\widehat{ \beta }_{ \widehat{M} }$ estimate consistently?". These considerations are relevant to the literature of variable and model selection, which includes: methodologies for determining the statistical significance of predictors in settings with many covariates (see, chudik2018one) as well as the uniform inference post-model selection which includes among others the studies of belloni2016post and farrell2015robust).

Here we focus on the framework proposed by kuchibhotla2021uniform and present key results. Suppose that for any $M \subset \left\{ 1,2,...,p \right\}$ the ordinary least-squares (OLS) empirical risk (or objective) function is

Then, by expanding the square function, it is clear that

Furthermore, define the following quantities

Then, the least-squares linear regression estimator $\hat{\beta}_{n,M}$ is defined as below

The notation arg min$_{ \theta } f( \theta )$ denotes the minimized of $f( \theta )$. Based on the quadratic expansion of the empirical objective function $\hat{R}_n \left( \theta ; M \right)$, the estimator $\widehat{\beta}_{n,M}$ is given by the closed form expression below

assuming nonsingularity of $\hat{ \Sigma }_n (M)$. Therefore, it is clear that $\widehat{\beta}_{n,M}$ is a smooth (nonlinear) function of two averages $\hat{ \Sigma }_n (M)$ and $\hat{ \Gamma}_n(M)$. Then, assuming that the random vectors $( X_i, Y_i )$ are i.i.d with finite fourth moments, it follows that $\hat{ \Sigma }_n (M)$ and $\hat{ \Gamma}_n(M)$ converge in probability to their expectations. Thus, we can define the expected value in terms of matrices and vectors as below

Following the notation above $\widehat{\beta}_{n,M} = \left[ \hat{ \Sigma }_n (M) \right]^{-1} \hat{ \Gamma}_n (M)$, and if $\big( \hat{ \Sigma }_n - \hat{ \Sigma }_n , \hat{ \Gamma}_n - \Gamma_n \big) \overset{ p }{ \to } 0$ as $n \to \infty$, using a Slutsky-type argument, it follows that $ \left( \widehat{\beta}_{n,M} - \beta_{n,M} \right) \overset{ p }{ \to } 0$, as $n \to \infty$, where $\beta_{n,M}$ is:

Thus, the convergence statement above only concerns a single submodel $M$ and is not uniform over $M$. By uniform-in-submodel $\left\lVert.\right\rVert_2-$norm consistency of $\hat{ \beta }_{n,M}$ to $\beta_{n,M}$, for $M \in \mathcal{M} (k)$, we mean that

Thus, converges of $\widehat{\beta}_{ n,M } $ to $ \beta_{ n, M }$ only requires convergence of $\widehat{\Sigma}_n(M)$ to $\Sigma_n(M)$ and $\widehat{\Gamma}_n(M)$ to $\Gamma_n(M)$.

Non-nested Regressions and Variable Selection Specification Testing

Following the framework proposed by hagemann2012simple, suppose that we observe covariates

that give rise to $M \geq 2$ different possible linear regression models for $y := ( y_1, ..., y_n )^{\top} \in \mathbb{R}^n$, such that

where $X_m := \big( x_{1,m},..., x_{n,m} \big)^{\top} \in \mathbb{R}^{ n \times d_m }$ is the design matrix of model $m$.

Then, the matrices $X_1,...., X_M$ are assumed to be non-nested, which implies that for any two matrices with index $m \neq \ell \in \mathcal{M}$, no matrix can be obtained by another by a linear transformation. Furthermore, the particular assumption does not rule out the possibility that some of the columns of $X_m$ and $X_{\ell}$ are identical or that they may be nonlinear transformations of another. In particular, suppose that there is an unobserved design matrix $X_{ m^{*} } := \big( x_{ 1, m^{*} },..., x_{n, m^{*} } \big)^{*}$ associated with the correct model $m^{*}$. Suppose that $\mathcal{F} := \big\{ X_m : m \in \mathcal{M} \cup \left\{ m^{*} \right\} \big\}$ has the following properties:

Next, we impose a condition for the existence of a correct model specification.

The proposed $J$ test presumes that for some predetermined $m \in \mathcal{M}$, the null hypothesis is formulated such as $m = m^{*}$ against the alternative hypothesis $ m \neq m^{*}$ in the presence of non-nested alternatives $\ell \in \mathcal{M} \backslash \left\{ m \right\}$. To do this, we artificially nesting the models via an additional parameter vector $\alpha_{m} := \left( \alpha_{\ell, m} \right)_{\ell \in \mathcal{M} \backslash \left\{ m \right\} } \in \mathbb{R}^{M-1}$ such that

Therefore, since the vectors $\left( \alpha_{\ell, m}, \beta_{\ell} \right)_{ \ell \in \mathcal{M} \backslash \left\{ m \right\} }$ of the nesting model may not be identified, we can replace the $\beta_{\ell}$ by the OLS estimates such that $\hat{\beta}_{\ell} = \left( X_{\ell}^{\top} X_{\ell} \right)^{-1} X_{\ell}^{\top} y$. (see, hagemann2012simple).

After redefining the error term appropriately, this yields that

Therefore, a Wald test for testing the null hypothesis that $\alpha_m = 0$ is a $J$test for the validity of model $m$ in the presence of the alternatives $\mathcal{M} \backslash \left\{ m \right\}$. Thus, to construct the test statistic, let

are the usual projection matrices and $\hat{\Omega}_{n,m}$ is an estimate of the conditional expectation of $\mathbb{E} \left( u_m^{*} u_m^{*\top} | \mathcal{F} \right)$. The $J$ test statistic for model $m$ is then given by

In other words, the hypothesis that the model $m$ is the true model is rejected for large values of $J_{n,m}$. The following asymptotic results holds for obtaining the empirical size and power of the test statistic.

Divide and Conquer Variable Selection Algorithm

Further algorithmic procedures for model selection purposes in high dimensional linear regression models includes the framework of meinshausen2010stability propose a variable selection methodology where the algorithm repeatedly employs a subset selection and apply the backward algorithm. Then, the final stage provides all the significantly selected covariates with a pre-specified control rate. Therefore, the particular stability selection procedure is designed to address the issue of choosing the amount of regularization such that a certain familywise type I error rate in multiple testing can be conservatively controlled for finite sample size. The issue of correlated covariates can affect the behaviour of stability selection for high correlated designs. Specifically, the stability selection algorithm puts a large emphasis on avoiding false positive selections and, as a consequence, might miss important variables if they are highly correlated with irrelevant variables. Moreover, battey2018distributed (see, also fan2015power) consider the so-called "divide and conquer" variable selection algorithm as we briefly explain below.

On each subset $\mathcal{D}_j$, we compute the debiased estimator of $\boldsymbol{\beta}^{*}$ such that

where $d$ is used to indicate the debiased version of the estimator, such that $M^{(j)} = \big( \boldsymbol{m}_1^{(j)},..., \boldsymbol{m}_d^{(j)} \big)^{\top}$ and $\boldsymbol{m}_v$ is the solution of

The above procedure implies a divide-and-conquer Wald statistic of the following form

for $\beta_v^{*}$, where $\bar{\sigma}$ is an estimator for $\sigma$ based on the $k$ subsamples. Consider the desparsified estimator for sub-sample $\mathcal{D}_j$ is given by

where $\widehat{\Theta}^{(j)}$ is a regularized inverse of the Hessian matrix of second-order derivatives of $\ell_{ n_k }^{ (j) } ( \boldsymbol{\beta} )$ evaluated at $\boldsymbol{\beta}^{ \lambda } ( \mathcal{D}_j )$, denote by

The proposed approach for estimating $\widehat{J}^{ (j) }$, reduces to the empirical covariance of the design matrix in the case of the linear model. The, the aggregated debiased estimator over the $k$ subsamples is

Thus, to approximate the nodewise Lasso requires to approximately invert $\widehat{J}^{(j)}$ via $L_1 -$regularization. The basic idea is to find the regularized invert row via a penalized $L_1-$regression, which is the same as regressing the variable $X_{ \nu }$ on $\boldsymbol{X}_{ - \nu }$ but expressed in the sample covariance form. In particular for each row $\nu \in \left\{ 1,..., d \right\}$, consider the optimization as below:

where $\widehat{J}_{ \nu, - \nu }^{(j)}$ denotes the $\nu-$th row of $\widehat{J}^{(j)}$ without the $( \nu, \nu )-$th diagonal element and $\widehat{J}_{ -\nu, - \nu }^{(j)}$ is the principal submatrix without the $\nu-$th row and $\nu-$th column.

Introduce the following matrix

Denote with

Moreover, the estimation of the matrix $\widehat{ \Theta }^{(j)}$ is given by

Therefore, we consider establishing the limit distribution of the following term

for any $\nu \in \left\{ 1,..., d \right\}$ under the null hypothesis that $H_0: \beta_{\nu} = \beta_{ \nu }^H$.

\paragraph{Open Problems}

An open problem in econometrics and statistics when considering model comparison and selection methodologies remains the aspect of robust statistical inference for nonnested environments. More specifically, various studies in the literature have considered methodologies for establishing the consistency of the Bayes factor in these settings. A particular interesting solution is the approach of converting a nonnested problem to a nested statistical problem which can provide a pseudo-distance between the base model and the full model. Another approach is to have on overlapping window where the comparisons between the nested models is made (see, berger1996intrinsic, berger1999default). Moreover, the framework proposed by vuong1989likelihood relates the probabilistic model selection approach to the classical nested-hypothesis testing situation. Lastly, in the next section we discuss an application from the statistics literature on sample-splitting to distinguish it from subsampling. Although the particular implementation corresponds to i.i.d data, an extension to dependent data (e.g., time series data with a weak form of dependence) could be an interesting fruitful avenue for future research.

Sample-Splitting and Variable Importance Algorithm

In this section, we follow and discuss the framework proposed by rinaldo2019bootstrapping. Specifically, consider a distribution-free regression framework, where the pair $Z = (X,Y) \in \mathbb{R}^d \times \mathbb{R}$ of $d-$dimensional covariates and response variables has an unknown distribution belonging to a large non-parametric class of $\mathcal{Q}$. Then, using minimal assumptions (assumption-lean inference) on the regression function

where $\mu(x)$ describes the relationship between the vector of covariates and the expected value of the response variable. We observe the set of data observations $\mathcal{D}_n = \left( Z_1,...,Z_n \right)$ where $Z_i = ( X_i, Y_i ) \in \mathbb{R}^{d+1}$, where $i = 1,...,n$ and the class $\mathcal{Q} = \mathcal{Q}_n$, which may depend on the sample size. Then, we apply to the data a procedure $w_n$, which returns both a subset of the covariates and an estimator of the regression function over the selected covariates. Formally, we have that

where $\widetilde{S}$ is the selected model, is a random, nonempty subset of $\left\{ 1,...,d \right\}$ and $\widehat{\mu}_{\widetilde{S}}$ is an estimator of the regression function $x \in \mathbb{R}^d \mapsto \mathbb{E} \left[ Y | X_{\widetilde{S}} = x_{\widetilde{S}} \right]$ restricted to $\widehat{S}$, where $(X,Y) \sim P$ independent of $\mathcal{D}_n$ and, for a vector $x = \left( x(1),..., x(d)\right) \in \mathbb{R}^d$, we set $x_{\widehat{S}} = \left( x(j), j \in \widehat{S} \right)$. The only assumption that is imposed on $w_n$ is that the maximum size of the selected model is controlled by the experimenter, that is, $1 \leq | \widehat{S} | \leq k$, for a predefined positive integer $k \leq d$, where $k$ and $d$ can both increase with the sample size. Therefore, rinaldo2019bootstrapping propose a methodology for defining $\widehat{S}$ such that it contains any optimal model. In particular, their framework allows for arbitrary procedures including sparse variable selection, and stepwise-forward regression. Thus, the goal of the framework proposed by rinaldo2019bootstrapping, is to provide statistical guarantees for various measures of variable importance applied to the covariates in $\widehat{S}$, uniformly over the choice of $w_n$ and over all the distributions $P \in \mathcal{Q}_n$. Then confidence sets are considered for four random parameters taking values in $\mathbb{R}^{\widehat{S}}$, each providing a different assessment of the level of statistical significance of the variables in $\widehat{S}$ from a purely predictive standpoint. All the random parameters under consideration are function of the data generating distribution $P$, of the sample $\mathcal{D}_n$ and its size $n$, as well as of the mechanism for model selection and the estimation procedure of $w_n$.

The projection parameter $\beta_{\widehat{S}}$

Consider a linear estimator of the form $x \mapsto \widehat{\mu}_{\widehat{S}} (x) = \widehat{\beta}^{\top}_{\widehat{S}} x_{\widehat{S}}$, where $\widehat{\beta}_{\widehat{S}}$ is any estimator of the linear regression coefficients for the selected predictors based on OLS. Then, the linear projection parameter $\widehat{\beta}_{\widehat{S}}$ is defined to be the vector of coefficients of the best linear predictor of $Y$ using $X_{\widehat{S}}$, given by

where $\mathbb{E}_{X,Y}$ denote the expectation with respect to the joint distribution of the random variables $\left( X,Y \right)$.

In other words, based on the framework proposed by rinaldo2019bootstrapping, these projection parameters correspond to the fact that $X^{\top} \beta_{\widehat{S}}$ is the $L_2$ projection of $Y$ into the linear space of all random variables that can be obtained as linear functions of $X_{\widehat{S}}$. The projection parameter is well-defined even though the true regression function $\mu$ is not linear. Indeed, it is immediate that

where

The LOCO parameters $\gamma_{\widehat{S}}$ and $\phi_{\widehat{S}}$

A commonly used measure of the importance of the selected covariates is $\beta_{\widehat{S}}$, but there are of course other ways to quantify variable significance. We consider two parameters of variable importance, which we refer to as Leave Out COvariate Inference-or LOCO-parameters (see, rinaldo2019bootstrapping).

The first LOCO parameters is $\gamma_{\widehat{S}} = \left( \gamma_{\widehat{S}} (j) : j \in \widehat{S} \right)$, where

In the last expression, $\widehat{\beta}_{\widehat{S}}$ is any estimator of the projection parameter $\beta_{\widehat{S}}$ and $\widehat{S}(j)$ and $\widehat{\beta}_{\widehat{S}(j)}$ are obtained by rerunning the model selection and estimation procedure after removing the $j-$th predictor. In other words, for each $j \in \widehat{S}, \widehat{S}(j)$ is a subset of size at most $k$ of $\left\{ 1,...,d \right\} \ \left\{j\right\}$. Notice that the selected model can be different when the $j-$th covariate is held out from the data, so that the intersection between $\widehat{S}(j)$ and $\widehat{S}$ can be smaller than $k-1$. The interpretation of $\gamma_{\widehat{S}}(j)$ is simple: it is the increase in prediction error by not including the $j-$th predictor in the model. For instance, it is easy to extend the definition of this parameter by leaving out several variables from $\widehat{S}$ at once without additional conceptual difficulties.

Moreover, the parameter $\gamma_{\widehat{S}}$ has advantages over the projection parameter $\beta_{\widehat{S}}$ and this is because it refers directly to prediction error and as also demonstrate the accuracy of the Normal approximation and the bootstrap is much higher. The second type of LOCO parameters that we consider are the median LOCO parameters, $\phi_{\widehat{S}} = \left( \phi_{\widehat{S}}(j), j \in \widehat{S} \right)$ with

As with $\gamma_{\widehat{S}}$, we may leave out multiple predictors at the same time.

The prediction parameter $\rho_{\widehat{S}}$

Another interesting measure of variable importance is an omnibus parameter that measures how well the selected model will predict future observations. To this end, we define the future prediction error as

where $\widehat{\beta}_{\widehat{S}}$ is computed based on $\mathcal{D}_n$. In particular, the main idea of the approach proposed by rinaldo2019bootstrapping relies on sample splitting: assuming for notational convenience that the sample size is $2n$, we randomly split the data $\mathcal{D}_{2n}$ into two halves, $\mathcal{D}_{1,n}$ and $\mathcal{D}_{2,n}$. Next, we run the model selection and estimation procedure $w_n$ on $\mathcal{D}_{1,n}$, obtaining both $\widehat{S}$ and $\widehat{\mu}_{\widehat{S}}$. We then use the second half of the sample $\mathcal{D}_{2,n}$ to construct an estimator $\widehat{\theta}_{\widehat{S}}$ and a confidence set $\widehat{C}_{\widehat{S}}$ for $\theta_{\widehat{S}}$ satisfying the following properties:

Confidence sets for the projection parameters: The bootstrap

Following rinaldo2019bootstrapping, the confidence set based on the Normal approximation require the evaluation of both the matrix $\widehat{\Gamma}_{\widehat{S}}$ and the quantile $\widehat{t}_{\alpha}$ which may be computationally inconvenient. The authors show that the paired bootstrap can be deployed to construct analogous confidence sets, centered at $\widehat{\beta}_{\widehat{S}}$, without knowledge of $\widehat{\Gamma}_{\widehat{S}}$. Thus, the bootstrap distribution corresponds to the empirical probability measure associated to the subsample $\mathcal{D}_{2,n}$ and conditionally on $\mathcal{D}_{1,n}$ and the outcome of the sample splitting procedure. Denote $\widehat{\beta}_{\widehat{S}}^{*}$ the estimator of the projection parameters $\beta_{\widehat{S}}$ arising from i.i.d sample of size $n$ drawn from the bootstrap distribution. For a given $\alpha \in (0,1)$, let $\widehat{t}^{*}_{\alpha}$ be the smallest positive number

Next, let $\left( \widetilde{t}_j^{*}, j \in \widehat{S} \right)$ be such that

By the union bound, each $\widetilde{t}^{*}_j$ can be chosen to be the largest positive number such that

Consider the following two bootstrap confidence sets:

It is immediate that $\widehat{C}^{*}_{ \widehat{S} }$ and $\widetilde{C}^{*}_{ \widehat{S}}$ are just the bootstrap equivalent of the confidence sets above. Therefore, the case of sparse fitting where $k = \mathcal{O}(1)$ implies that the size of the selected model is not allowed to increase with $n$. The standard central limit theorem shows that

Note that $\Gamma$ can be consistently estimated by the sandwich estimator $\widehat{\Gamma} = \widehat{\Sigma}^{-1} A \widehat{\Sigma}^{-1}$, where $A = n^{-1} \mathbb{X}^{\top} R \mathbb{X}$, $\mathbb{X}_{ij} = X_i(j)$, $R$ is the $k \times k$ diagonal matrix with $R_{ii} = \left( Y_i - X_i^{\top} \widehat{\beta} \right)^2$. By Slutsky's theorem, valid asymptotic confidence sets can be based on the Normal distribution with $\widehat{\Gamma}$ in place of $\Gamma$.

However, as rinaldo2019bootstrapping points out there is a clear prediction/accuracy trade-off when employing the sample splitting approach for variable selection. In other words, the selected model may be less accurate because only part of the data are used to select the model. Thus, although splitting creates gains in accuracy and robustness for inference it with some loss of prediction accuracy.

\paragraph{Open Problems} An interesting application would be to consider the feasibility of using sample splitting as an estimation and inference approach in time series regression models. For instance as in the case, when quantile regression models are employed for modelling risk measures such as the VaR and the CoVaR. In that case, sample splitting can be seen as a methodology for out-of-sample forecasting of the unknown risk quantities. Thus, developing robust methodologies for accommodating not only the uncertainty induced by the selection of the optimal size of the out-of-sample period but also to be able to capture correctly the persistence properties of regressors included in the model. Another example, consider the block bootstrap which is commonly which as a resampling method to preserve the dependence structure in predictive regression models. Thus, in the case of time-varying predictive regression models, essentially one models the effect of time-varying persistence using a rolling window. Combining these sets of information in a meaningful way is crucial in understanding how the persistence properties of predictors affect the asymptotic efficiency of statistics such as the bootstrap confidence set or other methods of variable selection in high-dimensional settings with dependent data.

Multiple Testing Procedure and Variable Selection

In this section, we consider an alternative methodology proposed in the statistical literature especially, since the pioneered work of Abraham Wald on Sequential tests of statistical hypotheses. Towards the direction of the multiple testing approach applied to high dimensional regression models for variable selection and statistical inference purposes, a relevant framework is proposed by chudik2018one. Therefore, in this section we discuss their proposed statistical methodology, so-called One Covariate at a Time Multiple Testing (OCMT) (see, also zhang1993model). Specifically, the particular procedure is computationally simple and fast even for extremely large data sets (see, also romano2005exact).

Thus, when a Lasso penalization is allowed then the estimator of $\boldsymbol{\beta}$ is obtained via

A certain degree of sparsity is required in order to apply standard penalized linear models in high-dimensional settings. However, the methodology does not require that the regressor vector $\boldsymbol{x}_{nt}$ to have a sparse covariance matrix, and therefore it is still applicable even if the covariance matrix of the noise variables, is not sparse. The particular algorithm can be thought as a model selection device since the main idea of this procedure is to test the statistical significance of the net contribution of all $n$ available potential covariates in explaining $y_t$ individually, while accounting for the multiple testing nature of the problem under consideration. However, the OCMT procedure is not sequential and selects in a single step all covariates whose $t-$ratios exceed a given threshold. Moreover, using the OCMT procedure post-selection is only applied when there are still covariates whose net contribution to $y_t$ is zero, despite the fact that they belong to the true model for $y_t$. From the statistics perspective sequential model selection is discussed by fithian2015selective. These approaches construct regression models by selecting variables from active sets, based on a sequence of sets.

Rejection Principle of Familywise Error Control

From the statistical theory perspective the $\textit{partitioning principle}$ is a powerful tool in multiple decision theory. The following theorem explains the concept of familywise error rejection, (FWER).

Closed testing and partitioning are recognized as fundamental principles of familywise error control, (FWE). In practice, various multiple testing procedures that control the FWE are often sequential, in the sence that rejection of some of the hypotheses may make rejection of the remaining hypothesis easier. Sequential rejective procedures in the literature include that of romano2005exact (see, also politis1999subsampling). Moreover, the paper of goeman2010sequential presents a unified approach to the class of sequentially rejective multiple testing procedures, emphasizing the sequential aspect. Thus, the authors consider a general sequentially rejective procedure as a sequence of single-step methods, determined by a rule for setting the rejection regions for each null hypothesis based on the current collection of unrejected regions for each null hypothesis based on the current collection of unrejected null hypotheses. Thus, the sequential rejection principle, implies that a single-step familywise error controlling procedure is turned into a sequential one which is a general principle of familywise error control.

Subgroup Selection Methodology

Lastly, we briefly discuss the subgroup selection methodology proposed by reeve2021optimal\footnote{Dr. Henry Reeve gave a seminar with title: "Subgroup Selection in nonparametric regimes", at the S3RI Departmental Seminar Series at the University of Southampton on the 10th of November 2022.}. In subgroup selection, the main objective of the statistician is to leverage the underline structure dynamics in data to identify a subset $\hat{\mathcal{A}}$ of the population to treat. Suppose that we have a distribution $P$ on a covariate-response pairs $(X,Y)$ in $\mathbb{R}^d \times \mathbb{R}$. Let $\mu := \mu_P$ denote the marginal distribution of the covariate $X \in \mathbb{R}^d$. Let $\eta : \mathbb{R}^d \to \mathbb{R}$ be the regression function defined by $\eta(x) := \mathbb{E} \big[ Y | X = x \big]$ for some $x \in \mathbb{R}^d$. Then, one would like to select a subgroup $A \subset \mathbb{R}^d$ such that $\eta$ is above a user-specified threshold $\tau \in \mathbb{R}$ on $A$. Hence, we are interested in subsets of the $\tau-$super level set such that the following set holds

In other words, the statistical mechanism of reeve2021optimal implies that the user chooses a data-dependent subgroup $\hat{A} \equiv \hat{A} ( \mathcal{A} )$, which is a random subset of $\mathbb{R}^d$, taking values in $\mathcal{A}$. Specifically, the practitioner has access to a sample $\mathcal{D} := \big\{ (X_1,Y_1),..., (X_n,Y_n) \big\} \overset{ \textit{i.i.d} }{ \sim } P$. Therefore, from the objective function perspective the aim is to select a data-dependent subgroup $\hat{A} ( \mathcal{D} ) \subset \mathcal{X}_{\tau} (\eta)$, with high-probability where $\mathcal{X}_{\tau} (\eta) := \big\{ x \in \mathbb{R}^d: \eta(x) \geq x \big\}$.

Therefore, our objective is to choose $\hat{A} \equiv \hat{A} ( \mathcal{D} )$ which minimizes regret $R_{\tau} ( \hat{A} )$, subject to

Further details on the implementation of this methodology as well as related statistical theory and numerical illustrations are presented in reeve2021optimal. More recently the isotonic subgroup selection approach is proposed by muller2023isotonic.

Thus, the statistical problem corresponds to minimizing the regret that controls the type I error. Although one in practice needs to check the exact set that corresponds to the minimized regret. Therefore, the following expression holds:

Then, the Type I error guarantee implies that the statistical subgroup selection procedure is constructed by sequential partial ordering of multiple covariates. Then, the statistician chooses only the one null hypothesis with the highest partial ordering. In other words, using a graph with acyclic structure allows to expand the multiple testing procedures (the null hypothesis of all ancestors is true). Thus, if the null at that node is false, then the null hypothesis associated with all ancestors is also false.

Model Selection in Cointegrating Regressions

According to mendes2011model, usually under the presence of cointegration in the VAR system, which implies that there are $(k-m)$ linear combinations of $X_t$ which are $I(0)$, a bias analysis is necessary to obtain a robust econometric estimation methodology. To do this, we consider the $m$ common stochastic trends and at most $( k - m )$ cointegrating relations amongst the $I(1)$ components are represented by the same generating process as $\Delta X_{t-1} = ( A - I_k ) X_{t-1} + \varepsilon_t$, where $\epsilon_t \sim \mathcal{N} ( \boldsymbol{0}, \boldsymbol{\Omega} )$. In particular, these $m$ components are represented by $m$ unit roots, and the $(k-m)$ stable roots of the stochastic difference equation represent the $I(0)$ components of $X_t$ as well as cointegrating relations between the $I(1)$ variates of $X_t$ (see, abadir1999influence).

Model Selection Consistency and Oracle Property

We focus on showing that under certain conditions on $n, p$ and $\lambda'$s the Adaptive Lasso selects the correct subset of variables - sign consistency and it has the oracle property, meaning that our estimate has the same asymptotic distribution of the OLS as if we knew a priori the selected model and at optimal rate.

Hence,

Further aspects of consideration in the more recent literature include the development of robust frameworks for time series regressions within a high-dimensional environment for the purpose of estimation, inference and forecasting (see, gupta2019robust, baillie2022robust). Moreover, the development of a framework that allows for the use of an ultra-high dimensional environment in nonstationary time series models can be useful when considering statistical properties such as model selection consistency and hypothesis testing accuracy. An econometric framework for high-dimensional quantile predictive regressions is proposed by fan2023predictive. Various studies consider methods for variable screening in high dimensional linear models, however less attention is paid on how these methodologies perform in the case of nonstationary time series models, especially when for the conditional quantile functional form. Relevant research aspects include the development of a statistical mechanism for dimension reduction or screening methodology in the presence of a high-dimensional vector of possibly nonstationary predictors, which can improve forecasting performance (see, pitarakis2023novel and gonzalo2023out). Moreover, choosing variables from irrelevant cointegrating relations is another important issue (see, khalaf2020simulation and richard2023model).

Model Selection and Rank of a Matrix

In many econometric and statistic applications knowing the rank of a matrix is crucial to ensure robust inference and testing. For example, from the financial economics literature one can test the implications of Arbitrage Pricing Theory by testing the corresponding rank restrictions. Moreover, the identification of parameters in econometric models depends on conditions and restrictions regarding the rank of related moment matrices such as the Jacobian matrix.

Specifically, various tests for the rank of a matrix can be found in the literature. In particular, these testing methodologies are constructed under the assumption that the unrestricted matrix estimator has a kronecker covariance matrix, however this approach can be sensitive to the ordering of the variables (see, kleibergen2006generalized). Furthermore, cragg1997inferring propose suitable testing procedures for determining the rank of a matrix. In particular, the authors examine the use of model selection criterion and sequential hypothesis testing methods to estimate the rank consistently. On the other hand, kleibergen2006generalized propose a novel rank statistic which uses a $\sqrt{n}-$consistent estimator of the unrestricted matrix which does not have to have a kronecker covariance matrix. To do this, the authors decompose the estimator of the unrestricted matrix using the spectral vector decomposition. Testing for the rank of the matrix can be also used in cointegration testing for non-stationary time series models. Specifically, in this case the limiting distribution is found to be functional of Brownian motions and is equal to the asymptotic distribution of the johansen1991estimation trace test.

Statistical Learning Methods in Time Series Analysis

Motivation

In this section we discuss some key applications of statistical learning methodologies in time series analysis. In particular the use of Neural Networks is now widely spread in the finance, economics and actuarial statistics fields. Some applications worth mentioning include:

Furthermore, on the aspect of modelling heterogeneity in mortality, more specifically the statistician can employ a functional form represented by $f$ ("biometric" function) to represent the age pattern of mortality. Then, the following representation follows:

where $\big\{ f^{(1)}, f^{(2)},..., f^{(N)} \big\}$ represent risk factors (such as health status, occupation etc.) in order to capture the non-linear effects of mortality over time. In particular, these risk factors are considered as contributing factors that worsens the mortality level or probability of an individual that can be incorporated into a rating system. Therefore, splitting individuals with similar risk into groups allows to gain more information about uncertainty in heterogeneity\footnote{Professor Ermanno Pitacco gave a seminar with title: "Heterogeneity in mortality: A survey with an actuarial focus" at the S3RI Departmental Seminar Series at the University of Southampton on 24 of September 2018.} (which is a well-known principle in the statistical literature of credibility theory and risk premium models). Moreover, according to pitacco2019heterogeneity a rating system can be constructed based on the form $q_{t + x}^{spec} = q_{t+h} \cdot \left( 1 + \sum_{j=1}^r \rho^{(j)} \right)$, in which case peak mortality is achieved after entering disable stage, such that $q_{t + x}^{spec} = q_{t+x} \cdot \Delta ( x + h; \alpha, \gamma )$. Then, the unobservable risk factors can be captured via the form $\mu_{x}^{spec} = \Phi \left( \mu_{x+t}; z_{x,t} \right)$; known as random heterogeneity.

The aforementioned examples provide some brief illustrations of the various applications from the econometrics and statistics literature that can motivate the investigation of relevant research questions both from the theoretical as well as the empirical perspective when we consider statistical learning methods for time series analysis purposes. Further resources related to deep learning theory and applications can be found in goodfellow2016deep (see, also lecun2015deep).

Non-Asymptotic Probability Theory

One of the main goals of learning theory is the development of stability bounds of algorithmic procedures. Following the framework of mohri2010stability we employ the following definitions.

Specifically, the use of stability allow us to derive generalization bounds and an exponential concentration bounds of the following form

where the probability is over a sample of size $m$ and where $\frac{\tau}{m}$ is the Lipschitz parameter of $\Phi$, with $\tau$ a function of $m$. In the first section below, we consider the problem of shallow-vs-deep expressiveness from the perspective of approximation theory and general spaces of functions having derivatives up to certain order (Sobolev-type spaces). Specifically, in this framework the problem of expressiveness is very well studied in the case of shallow networks with a single hidden layer, where it is known, in particular, that to approximate a $C^n-$function on a $d-$dimensional set with infinitesimal error $\epsilon$ one needs a network of size about $\epsilon^{ - d / n}$, assuming a smooth activation function.

Shallow Neural Network Estimate learned by Gradient Descent

Consider functions of the form

where $w_0,...,w_d \in \mathbb{R}$ the weights of the neuron and $\sigma: \mathbb{R} \to \mathbb{R}$ the activation function. Shallow Neural Networks have only one hidden layer where a simple linear combination of neurons is used to define a function $f: \mathbb{R}^d \to \mathbb{R}$ by

Denote the weight between neuron $j$ in layer $(s-1)$ and neuron $i$ in layer $s$ by $w_{i,j}^{(s)}$. This leads to the following recursive definition of a neural network with $\mathcal{L}$ layers and $k_s$ neurons in layer $s \in \left\{ 1,..., \mathcal{L} \right\}$ such that the following representation applies:

for some $w_{1,0}^{ ( \mathcal{L} ) },..., w_{1, k_L }^{ ( \mathcal{L} ) }$ and for $f_i^{ ( \mathcal{L} ) }$'s recursively defined by

Main results in the literature show that NNs can achieve dimension reducion provided the regression function is a composition of sums of functions, where the input dimension of each of the functions is at most $d^{*} < d$. Denote by $f_{\mathsf{net},w}$ the neural network with weight vector $\boldsymbol{w} = \left( w_{j,k}^{(s)} \right)_{ s= 0,.., \mathcal{L}, j=1,..., k_{s+1}, k = 0,..., k_s }$ and set the following function

for some randomly chosen initial vector $\boldsymbol{v}$.

Then the optimization problem can be written as below:

Learning of linear penalized least squares estimates by gradient descent

Let $( x_1,y_1 ),..., ( x_n,y_n ) \in \mathbb{R}^d \times \mathbb{R}$, let $K \in \mathbb{N}$ and let $B_1,..., B_K : \mathbb{R}^d \to \mathbb{R}$, with $c_1 > 0$. We consider the problem to minimize as below

where $\boldsymbol{a} = \left( a_1,..., a_K \right)^{\top}$ \ \ \ and \ \ \ $\left\lVert \boldsymbol{a} \right\rVert^2 = \sum_{j=1}^K a_j^2$, by gradient descent.

To obtain the solution of the optimization problem, we choose $\boldsymbol{a}^{(0)} \in \mathbb{R}^K$ and set with

for some chosen $\lambda_n > 0$.

Consider that a Lipschitz continuity condition holds

Then, we have that

Consider the Cauchy-Schawrz inequality and since the activation function $\sigma$ is Lipschitz continuous

Deep Neural Network Estimate learned by Gradient Descent

Based on the framework proposed by shen2021deep, we present the following results which are useful to investigate the properties of Deep Neural Networks learned by gradient descent with respect to their dimensionality and complexity. In particular, the concept of pseudo-dimension is considered as a measure of complexity (see, mohri2008rademacher). Thus, the framework of shen2021deep considers the implementation of deep neural network for a high-dimensional quantile regression. An application to forecasting in time series is presented by chronopoulos2023forecasting. A key ingredient to derive excess risk bounds for the deep quantile regression framework is to consider the properties of the functional classes within which identification and estimation holds.

Notice that the denominator can be improved to $n$ if the response $Y$ is assumed to be sub-exponentially distributed, that is, there exists a constant $\sigma_Y > 0$ such that $\mathbb{E} \big[ \mathsf{exp} \big( \sigma_Y | Y | \big) \big]$.

Moreover, further interesting aspects here include the technical results with respect to the DNN architecture (e.g., composition), the relation of continuous functionals with the underline quantile regression model and the number of variables as well as linearity of the model and how it relates with the dimensionality of the problem. In particular, one can use numerical experiments to evaluate the performance of the standard kernel-based method for quantile regressions against the DQR estimation approach.

Furthermore, for the excess risk of the DQR estimator and the error bounds for the models, based on appropriately specified network parameters (depth, width and size of the network), we have the following upper bound for the excess risk,

Then, the convergence rate which is part of the error bound, $n^{ - \left( 1 - \frac{1}{p} \right) \frac{2 \alpha^{*} }{ 2 \alpha^{*} + t^{*} } }$, is determined by the number of moments $p$ of the response $Y$, the smoothness index of the composite function $\alpha^{*}$ and the intrinsic dimension of the model $t^{*}$. Other relevant aspects include whether imposing further distributional assumptions (such as heavy-tailed errors versus Gaussian errors) can affect the limit theory of estimators as well as its min-max optimality properties. Next, we employ the risk function (i.e., the expectation of the loss function) such that $R_{\tau} (f) := \mathbb{E} \big[ \rho_{\tau} ( Y - f(x)) \big]$ as well as the empirical risk, which is the minimizer of the empirical risk function (i.e., it allows to restrict the functional class). For the proof of the above lemma (Lemma 1 in shen2021deep) the following two steps are necessary and sufficient:

Deep Neural Network Architecture Approximations

Related studies on the theoretical aspects of network achitecture include farrell2021deep, zeng2021deep and dung2021deep. The approximation by deep ReLU neural networks of functions having a mixed smoothness is related to the high-dimensional sparse-grid approach which was introduced by Zenger for numerical solving partial differential equations (dung2021deep). Next, we consider the mathematical analysis for deriving error bounds of DNNs (see also yarotsky2017error and yarotsky2018optimal).

Error Bounds of DNNs

Consider the following $\mathsf{g}_s$ function with $2^{s-1}$ uniformly distributed "triangles". The key observation here is that the function $f(x) = x^2$ can be approximated by linear combinations of the functions $\mathsf{g}_s$. Specifically, $f_m$ is a piece-wise linear interpolation of $f$ with $2^m + 1$ uniformly distributed knots (breakpoints) $\frac{k}{2^m}$, $k \in \left\{ 0,..., 2^m \right\}$, such that

Thus, the function $f_m$ approximates $f$ with the error $\epsilon_m = 2^{-2(m+1)}$. Applying the linear interpolation method from $f_{m-1}$ to $f_m$ amounts to adjusting it by a function proportional to a sawtooth function given by

Hence, it holds that

Topological Space for DNN Architecture

Then, we focus based on Assumptions 1 and 2 above, on deriving high-probability bounds.

Statistical Inference with Stochastic Gradient Descent

Literature Review

Notice that both for the estimation of Shallow as well as Deep Neural Networks, useful optimization algorithms under examination are the Gradient Descent and Stochastic Gradient Descent. Related studies that discuss the implementation and applications of these algorithms are presented by chen2020statistical. Furthermore, toulis2017asymptotic study the asymptotic theory analysis and finite-sample properties of these estimators (see also tran2015stochastic). Moreover, a selective overview of deep learning estimation methodologies are discussed in the study of fan2021selective and braun2021smoking.

Following chen2020statistical, assume that we want to estimate the true parameter $\theta_0 \in \mathbb{R}^p$ of a distribution $f$ from $\textit{i.i.d}$ data points $( X_i, Y_i )$ such that conditional on covariate $X_i \in \mathbb{R}^p$ outcome $Y_i \in \mathbb{R}^d$ is distributed according to $f( Y_i; X_i, \theta_0 )$. Such statistical problems reduce to optimization. More specifically, suppose that our aim is to estimate the true parameter $\theta_0 \in \mathbb{R}^p$ of a distribution $f$ from i.i.d data points $( X_i, Y_i )$ such that conditional on covariate $X_i \in \mathbb{R}^p$ outcome $Y_i \in \mathbb{R}^d$ is distributed according to $f( Y_i; X_i, \theta_0 )$. Such statistical problems reduce to optimization.

where $\gamma_n$ is the learning rate sequence.

Consider the equivalent moment condition below

where the expectation corresponds to the true conditional distribution of outcome $Y$ given covariate $X$.

Therefore, it holds that

Recall the Fisher Information matrix is denoted $\hat{\mathcal{I}}_n ( \theta ) = - \nabla^2 \ell \left( X_n^{\prime} \theta; Y_n \right)$. Using a Taylor approximation of the gradient $\nabla \mathsf{log} \left[ f \left( Y_n ; X_n, \textcolor{red}{\theta_n^{im}} \right) \right]$ yields

where $\Delta \textcolor{red}{\theta_n^{im}} = \textcolor{red}{\theta_n^{im}} - \theta_0$ and $\Delta \theta_n^{sgd} = \theta_n^{sgd} - \theta_0$.

Preliminary Theory

Throughout, we use $\left\lVertx\right\rVert_p$ to denote the $\ell_p-$norm of $x$, $\left\lVertX\right\rVert$ the matrix operator norm of $X$ and $\left\lVert X \right\rVert_{\infty} = \mathsf{max}_{i,j} | X_{i,j} |$ the elementwise $\ell_{\infty}-$norm of $X$. For any sequences $\left\{ a_n \right\}$ and $\left\{ b_n \right\}$ of positive numbers, we write with $a_n \geq b_n$ and $a_n \leq b_n$.

Consider the SGD method where the iteration is given by (see, chen2021statistical)

where $\xi_n := \nabla F( x_{n-1} ) - \nabla f ( x_{n-1}, \zeta_n )$. In particular, the above formulation decomposes the descent into two parts: $\nabla F( x_{n-1} )$ which represents the direction of population gradient which is the major driving force behind the convergence of SGD and $\xi_n$ is a martingale difference sequence under the above assumption. Furthermore, it holds that

Notice that $\mathbb{E}_n (.)$ denotes the conditional expectation $\mathbb{E}_n (. | \mathcal{F}_n )$, where $\mathcal{F}_n$ is the $\sigma-$algebra generated by $\left\{ \zeta_1,..., \zeta_n \right\}$. Let $\Delta_n := x_n - x^*$ be the error of the $n-$th iterate. An equivalent expression gives us

Given the SGD recursion and under suitable assumptions it can be shown that when the step size sequence $n_i = \eta i^{ -\alpha }$, for $i \in \left\{ 1,..., n \right\}$ with $\alpha \in (1/2, 1)$, we have that

where $\bar{\Delta}_n = \frac{1}{n} \sum_{i=1}^n \Delta_i = \bar{x}_n - x^{*}$. \color{black}

Assumptions and error bounds

Further assume that $\nabla^2 F( x^{*} )$ exists, and $\nabla F(x)$ is Lipschitz continuous with a constant $L_F$, that is for any $x_1$ and $x_2$

Estimators for asymptotic covariance

Consider two consistent estimators, the plug-in estimator and the batch-means estimator (see, chen2021statistical). Consider the sample estimate as below

Consider the following nodewise Lasso approach

where $D_{.,j}$ is the $j-$th column of the design matrix $D$ and $D_{.,-j}$ is the design submatrix without the $j-$th column (see also guo2022doubly). Further, one can estimate $\Omega_{j,j}$ by

Averaged SGD Algorithms

Consider the statistical problem of searching for the minimum point $\theta_0$ of a smooth function $f( \theta )$ where $\theta \in \Theta \subset \mathbb{R}^d$. The stochastic GD method provides a direct way to solve the minimization problem. The algorithm is given as below. Let $\theta_0 \in \mathbb{R}^d$ be the initial value and for $n \geq 1$, we update $\theta_n$ using

where $\theta_n = \frac{1}{n} \sum_{i=0}^{ n-1 } \theta_i$.

Sieve Estimation Applications

Consider the knowledge transmission through neural networks as a recursive process that focuses on the diffusion of processes. Similar to neurobiological processes in the brain neural networks can activate certain neuros in the brain for transmission of understanding.

Stochastic approximations recursively approximate the zeros of an unknown function $\Psi ( \theta )$, say $\theta^{*}$,

where $a_t$ is a "learning rate" tending to zero and $\psi ( Z_t, \theta )$ is a measurement of $\Psi ( \theta )$ at time $t$, influenced by random variables $Z_t$. Furthermore, when $\Psi ( \theta ) \equiv \mathbb{E} \big[ \psi ( Z_t, \theta ) \big]$ this method yields a recursive implementation of the method of m-estimation of Huber (see, kuan1994artificial). The particular methodology can be used to estimate recursively the parameters of nonlinear regression models, such as those arising in neural network applications. Specifically, offline nonparametric estimation methods can usually be applied to the case of mixing processes using results for the methods of sieves. On the other hand, online nonparametric estimation methods require convergence to a global optimum of the underlying least squares problem, not just the local optimum.

Asymptotic Properties of Sieve Estimation

Another application in which constraint optimization is used is when estimating and simulating neural networks with sieves. In particular, based on a construction of the neural network sieve estimators, in each sieve space $\mathcal{F}_{r_n}$, there is a constraint on the $\ell_1$ norm for which $\sum_{i=0}^{ r_n } | \alpha_i | \leq V_n$. Therefore, finding the nearly optimal function in $\mathcal{F}_{r_n}$ for $\mathbb{Q}_n (f)$ is in fact a constrained optimization problem. A classical way to conduct this optimization is through introducing a Lagrange multiplier for each constraint. On the other hand, due to the difficulty in finding an explicit connection between the Lagrange multiplier and the upper bound in the inequality constraint we use instead the subgradient method (see, shen2023asymptotic). The main idea is to update the parameter $\left\{ \alpha_0,...., \alpha_{r_n} \right\}$ through

This allows us to investigate the asymptotic properties, including consistency, rate of convergence and asymptotic normality for neural network sieve estimators with one hidden layer.

Using the Generalized Dominated Convergence Theorem, we have that

Sieve Estimation for Panel Data

A growing interest in the estimation of panel data models with cross-section dependence but most of the literature focuses on the linear specification of the regression relationship. In particular, let $y_{it}$, $t = 1,...,n$ and $t = 1,...,T$ be the $i-$th cross section unit at time $t$. We suppose that $y_{it}$ is generated according to the following semiparametric panel data generating process (see, su2012sieve)

where $x_{it} \in \mathcal{X}_i \in \mathbb{R}^d$ is a vector of observed individual-specific regressors on the $i-$th cross-section unit at time $t$. In this section, we focus on the sieve estimation of semiparametric panel data models with multi-factor error structure. We develop the asymptotic theory under fairly general conditions when both the cross-section and time-dimensions are large. For instance, if only homogeneous regression relationships are of interest, the time dimension need not pass to infinity. Moreover one can consider testing the constancy of the nonparametric relationship over individuals in the presence of multi-factor error structure. In addition the individual specific regressors have the following structure

In practice, one may also be interested in estimating a restricted submodel of the following form

Sieve M inference on irregular parameters

We follow the framework proposed by chen2014sieve who consider plug-in sieve M estimators. More precisely, we assume that the data $\left\{ Z_i \right\}_{i=1}^n$ is a random sample from the distribution of $Z$ defined on an underlying complete probability space. Let $\mathcal{L} (., .) : \mathcal{Z} \times \mathcal{A} \mapsto \mathbb{R}$ be a measurable function and $\mathbb{E} \left[ \mathcal{L} \left( Z, \alpha \right) \right]$ be a population criterion. For simplicity we assume that there is a unique $\alpha_0 \in \left( \mathcal{A}, d_{A} \right)$ such that $\mathbb{E} \left[ \mathcal{L} \left( Z, \alpha_0 \right) \right] > \mathbb{E} \left[ \mathcal{L} \left( Z, \alpha \right) \right]$ for all $\alpha \in \left( \mathcal{A}, d_{A} \right)$ with $d_A \left( \alpha, \alpha_0 \right) > 0$. Different models in economics correspond to different choices of the criterion function $\mathbb{E} \left[ \mathcal{L} \left( Z, \alpha \right) \right]$ and the parameter space $\left( \mathcal{A}, d_{A} \right)$. A model does not need to be correctly specified and $\alpha_0$ could be pseudo true parameter. In this paper, we are interested in the estimation of and inference of a functional $f \left( \alpha_0 \right)$ via the method of sieves.

Let $\mathcal{A}_n$ be a sieve space for the whole parameter space $\mathcal{A}$. Then an approximate sieve M estimator $\widehat{a}_n \in \mathcal{A}_n$ of $\alpha_0$ solves

In the example above, we have that

with $\alpha(Y) = X_{0,i}^{\prime} \theta + \sum_{j=1}^q h_{j,0}\left( X_{j,i} \right)$. For any $\alpha \in \mathcal{A}$, we define a strong metric $\left\lVert.\right\rVert_s$ as below

By the definitions of the metrics, we have that for any $\alpha \in \mathcal{A}$, (see, chen2014sieve)

Let $v^{*}_n = \left( v^{*}_{\theta, n}, v^{*}_{h, n} \right)$ be the Riesz representer of the functional $\frac{ \partial f \left( \alpha_0 \right) }{ \partial \alpha } [v]$ on $\mathcal{V}_n$. Denote with

Then, the variance of the plug-in sieve $M$ estimator $f \left( \widehat{\alpha}_n \right)$ of $f \left( \alpha_0 \right)$ is expressed as below

An application to specification testing in econometrics using invariance principles on Sobolev spaces is proposed by kuersteiner2019invariance. Furthermore, an application of general sieves for inference purposes in nonparametric time series regression models identified with conditional moment restrictions is studied by chen2022inference. In particular, nonlinear sieve learning employs a more general functional class of sieves that can approximate nonlinear functions of high dimensional variables (such as in the case of an infinite-dimensional parameter space).

Further Statistical Algorithms in Economic Applications

This section is motivated from the Seminar of Dennis Kristensen, Professor of Economics at the Department of Economics of University College London (UCL), who presented the paper with title: "Iterative estimation of structural models with an application to perturbed utility models", in May 2023 at the Department of Economics, University of Exeter Business School. Specifically, the particular steam of literature is based on the framework proposed by hotz1994simulation and is relevant to estimators for dynamic models of discrete choice, that is, a general class of iterative estimators.

Take for example, the classical likelihood estimation, which is a statistical methodology that iteratively solves the statistical problem until convergence to an estimate that is close to the true unknown parameter of interest. In the same spirit, without sacrificing efficient statistical estimation we assume that such a sequence of estimators converge to the full solution of the optimization problem. These optimization steps are summarized as below:

A different stream of literature considers doubly robust estimation methods which have been developed in various semiparametric problems, including partially linear models, instrumental variable analysis, and dimension reduction among others. Therefore, as a somewhat under-appreciated limit result, we point out that the familiar least-squares estimator for each individual coefficient in linear regression is doubly robust in the context of a partially linear model. This result is also closely related to debiased Lasso estimation in high-dimensional linear regression (e.g., see guo2022doubly).