$\epsilon:=\left(\epsilon_{1},...,\epsilon_{n}\right)^{\prime}$. Moreover,
$\left\\{\epsilon_{i},i\in\mathbb{Z}\right\\}$ is a stationary process with
autocovariance function

$\displaystyle\gamma_{\epsilon}(k)=\sum_{j=1}^{\infty}a_{j}a_{j+k}=k^{-1+2d}B(d,1-2d)$
(3.56)

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

$\displaystyle\hat{\beta}^{n}(\lambda)=\underset{\beta}{\text{arg
min}}\left\\{\frac{1}{n}\left\lVert
Y-X^{\prime}\beta\right\rVert_{2}^{2}+\lambda_{n}\left\lVert\beta\right\rVert_{1}\right\\},\
\ \lambda>0,$ (3.57)

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 Babii et al.,  (2022) 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.

####  3.5.1. 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

$\displaystyle
W_{nj}=n^{-(1/2+d)}\sum_{i=1}^{n}x_{ij}\epsilon_{i}=n^{-(1/2+d)}\sum_{i=1}^{n}\sum_{s=-\infty}^{i}x_{ij}a_{i-s}\zeta_{s}=\sum_{s=-\infty}^{n}c_{ns,j}\zeta_{s},$
(3.58)

where we have that

$\displaystyle c_{ns,j}$
$\displaystyle:=n^{-(1/2+d)}\sum_{i=1}^{n}x_{ij}a_{i-s},\ \ \
s\in\mathbb{Z},j=1,...,p,$ (3.59) $\displaystyle c_{n,j}$
$\displaystyle:=\underset{-\infty<s\leq n}{\text{sup}}\left|c_{ns,j}\right|,\
\ \ c_{n}=\underset{1\leq j\leq p}{\text{max}}c_{n,j}.$ (3.60)

Moreover, we denote with

$\displaystyle\sigma_{n,j}^{2}:=\text{Var}\left(W_{nj}\right),\ \
\sigma_{n}^{2}=\underset{1\leq j\leq p}{\text{max}}\sigma_{n,j}^{2}.$ (3.61)

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$,

$\displaystyle\frac{1}{n}\left\lVert
X\left(\hat{\beta}-\beta\right)\right\rVert_{2}^{2}+\lambda_{n}\left\lVert\hat{\beta}-\beta\right\rVert_{1}\leq\frac{4\lambda_{n}^{2}s_{0}}{\phi_{0}^{2}}$
(3.62)

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

$\displaystyle\Lambda=\left\\{\underset{1\leq j\leq
p}{\text{max}}\frac{2}{n}\left|\sum_{i=1}^{n}x_{ij}\epsilon_{i}\right|\leq\lambda_{0n}\right\\},$
(3.63)

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

######  Remark 11.

The presence of long memory errors can affect the convergence rates of
estimators as well as the variable selection procedure. Moreover, the presence
of serial correlation can be captures by modeling the error term such that
$\varepsilon_{t}=\rho\varepsilon_{t-1}+u_{t}$. In particular, there are cases
in which serial correlation can manifest as structural breaks in high
dimensional models (see, Kapetanios and Zikes,  (2018)). On the other hand,
the long memory property is a common feature of mean-revering processes, which
implies the significance of sample autocorrelations at large lags. In other
words, long memory processes implies that when predicting future values this
will depend persistently from the past observations. For instance, $d=0.5$ is
refereed to the fractional or long-memory parameter of the stochastic process.

######  Theorem 9.

For the long memory regression model suppose that the design variables satisfy
the model assumptions. Further, suppose that the tuning parameter
$\lambda_{n}$ is such that $\lambda_{n}\to\lambda_{0}\geq 0$, then we have
that

$\displaystyle\hat{\beta}_{n}$
$\displaystyle\overset{p}{\to}\mathsf{argmin}\big{(}Z(\phi)\big{)}$ (3.64)
$\displaystyle Z(\phi)$
$\displaystyle=\left(\phi-\beta\right)^{\prime}C\left(\phi-\beta\right)+\lambda_{0}\sum_{j=1}^{p}\left|\phi_{j}\right|,\
\ \phi\in\mathbb{R}^{p}$ (3.65)

Thus, if $\lambda_{n}=o(1)$ then $\mathsf{arg\
min}_{\phi}\left(Z(\phi)\right)=\beta$ and $\hat{\beta}_{n}(\lambda_{n})$ is
consistent for the unknown parameter $\beta$.

###### Proof.

To prove this theorem we consider the following objective function

$\displaystyle
Z_{n}(\phi)=\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-X_{i}^{\prime}\phi\right)^{2}+\lambda_{n}\sum_{j=1}^{n}\left|\phi_{j}\right|,$
(3.66)

then $Z_{n}(\phi)$ is convex. Therefore, we need to show the pointwise
convergence (in probability) of $Z_{n}(\phi)$ to $Z(\phi)+k^{2}$ for some
constant. Clearly, it holds that

$\displaystyle\lambda_{n}\sum_{j=1}^{n}\left|\phi_{j}\right|\to\lambda_{0}\sum_{j=1}^{n}\left|\phi_{j}\right|.$
(3.67)

Consider expanding the expression,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-X_{i}^{\prime}\phi\right)^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left[\epsilon_{i}-X_{i}^{\prime}\left(\phi-\beta\right)\right]^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}^{2}+\frac{1}{n}\sum_{i=1}^{n}\left(\phi-\beta\right)^{\prime}X_{i}X_{i}^{\prime}\left(\phi-\beta\right)-\frac{2}{n}\left(\phi-\beta\right)^{\prime}\sum_{i=1}^{n}X_{i}\epsilon_{i},$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}^{2}+\frac{1}{n}\sum_{i=1}^{n}\left(\phi-\beta\right)^{\prime}X_{i}X_{i}^{\prime}\left(\phi-\beta\right)-\frac{2}{n}\left(\phi-\beta\right)^{\prime}X^{\prime}\epsilon,$

Notice that the first term in the above expression converges to $k^{2}$ by the
ergodic theorem, since the error sequence $\left\\{\epsilon_{i}\right\\}$
forms a stationary ergodic sequence. The second term converges to
$\left(\phi-\beta\right)^{\prime}C\left(\phi-\beta\right)$ and the last term
converges to zero in probability, since we have that
$\frac{1}{n^{0.5+d}}X^{\prime}\epsilon$ converges in distribution. ∎

######  Theorem 10.

For the long memory regression model assume that the design variables satisfy
the model assumptions. Suppose that $n^{1/2-d}\lambda_{n}\to\lambda_{0}\geq 0$
as $n\to\infty$, then

$\displaystyle n^{1/2-d}\left(\hat{\beta}^{n}-\beta\right)$
$\displaystyle\to_{D}\underset{u}{\mathsf{arg\ min}}\ V(u),$ (3.68)
$\displaystyle V(u)$
$\displaystyle=-2u^{\prime}\mathcal{W}+u^{\prime}Cu+\lambda_{0}\sum_{j=1}^{p}\left[u_{j}\text{sign}\left(\beta_{j}\right)\mathbf{I}_{\left[\beta_{j}\neq
0\right]}+\left|u_{j}\right|\mathbf{I}_{\left[\beta_{j}\neq 0\right]}\right]$
(3.69)

such that $\mathcal{W}$ is an $\mathcal{N}\left(0,\Sigma\right)$ random
variable.

###### Proof.

We define the following expression

$\displaystyle
V_{n}(u)=n^{1-2d}\left\\{\sum_{i=1}^{n}\frac{1}{n}\left[\left(\epsilon_{i}-\frac{X_{i}^{\prime}u}{n^{\frac{1}{2}-d}}\right)^{2}-\epsilon_{i}^{2}\right]+\lambda_{n}\sum_{j=1}^{p}\left[\left|\beta_{j}+\frac{u_{j}}{n^{\frac{1}{2}-d}}\right|\right]-\left|\beta_{j}\right|\right\\}.$
(3.70)

Furthermore, we denote the first term of the above expression by (I) and the
second term by (II). Then, we obtain that

$\displaystyle(I)=n^{1-2d}\left\\{\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}^{2}-2\frac{\sum_{i=1}^{n}X_{i}^{\prime}u\epsilon_{i}}{nn^{1-2d}}+\frac{u^{\prime}\sum_{i=1}^{n}X_{i}X_{i}^{\prime}u}{nn^{\frac{1}{2}-d}}-\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}^{2}\right\\}$
$\displaystyle=\left\\{\frac{u^{\prime}\sum_{i=1}^{n}X_{i}X_{i}^{\prime}u}{n}-2\frac{\sum_{i=1}^{n}X_{i}^{\prime}u\epsilon_{i}}{n^{\frac{1}{2}+d}}\right\\}$
$\displaystyle\to u^{\prime}Cu-2u^{\prime}\mathcal{W},\ \ \text{as}\
n\to\infty,$

where $\mathcal{W}\sim\mathcal{N}\left(0,\Sigma\right)$. Moreover, we have
that

$\displaystyle(II)=n^{\frac{1}{2}-d}\lambda_{n}\sum_{j=1}^{p}\left[\left|n^{\frac{1}{2}-d}\beta_{j}+u_{j}\right|-n^{\frac{1}{2}-d}\left|\beta_{j}\right|\right]\to\lambda_{0}\sum_{j=1}^{p}\left[u_{j}\text{sign}\left(\beta_{j}\right)\mathbf{I}_{\left[\beta_{j}\neq
0\right]}+\left|u_{j}\right|\mathbf{I}_{\left[\beta_{j}\neq 0\right]}\right].$

∎

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

######  Lemma 9.

For each fixed $n$, let

$\displaystyle A$
$\displaystyle:=\left\\{\left|\sum_{i=-\infty}^{n}Y_{ni}\right|>r\right\\},\ \
\ B_{m}=\left\\{\left|\sum_{i=-m}^{n}Y_{ni}\right|>r-\delta\right\\},\ \
r>0,\delta>0,m=1,2,...$ $\displaystyle B$
$\displaystyle=\underset{m\to\infty}{\text{lim inf}}\ B_{m}.$

If $\left|\displaystyle\sum_{i=-\infty}^{n}Y_{ni}\right|<\infty$, almost
surely, then, for each fixed $n$, $A\subset B$.

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

$\displaystyle\left|\displaystyle\sum_{i=-\infty}^{n}Y_{ni}\right|>r,\ \ \
\left|\displaystyle\sum_{i=-\infty}^{n}Y_{ni}\right|<\infty.$

######  Remark 12.

Another relevant aspect for the lasso shrinkage approach regardless of the
properties of the error terms is the cross-validation consistency (see,
Chetverikov et al.,  (2021) and Yang,  (2007)). Moreover, although we discuss
various cases of the Lasso shrinkage algorithm, those correspond to a
different penalization property as well as possibly different-type of
thresholding, especially when one considers the pathwise selection approach
against the cross-validation approach.

###  3.6. Residual empirical process based on the ALasso

Moreover, the framework proposed by Chatterjee et al.,  (2015) 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

$\displaystyle\widehat{\boldsymbol{\beta}}_{n}=\underset{\mathbf{u}\in\mathbb{R}^{p}}{\text{arg
min}}\
\sum_{i=1}^{n}\left(y_{i}-\mathbf{x}_{i}^{\prime}\mathbf{u}\right)^{2}+\lambda_{n}\sum_{j=1}^{p}\frac{|u_{j}|}{\tilde{\beta}_{j,n}}^{\gamma},$
(3.71)

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}$.

####  3.6.1. Main results

##### Asymptotic uniform linearity in high dimensions

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

$\displaystyle\widehat{Z}_{n}(t)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\mathbf{1}\left(F(e_{i})\leq
t\right)-t\right]$ (3.72) $\displaystyle Z_{n}(t)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\mathbf{1}\left(F(\epsilon_{i})\leq
t\right)-t\right]$ (3.73)

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.

######  Theorem 11.

$\displaystyle\underset{t\in[0,1]}{\text{sup}}\
\left|\widehat{Z}_{n}(t)-\left[Z_{n}(t)-\bar{\mathbf{x}}_{n}^{\prime}\sqrt{n}\left(\widehat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}\right)J(t)\right]\right|=o_{p}(1).$
(3.74)

Notice that Theorem 11 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.

##### 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

$\displaystyle\mathcal{Z}_{n}^{OR}=n^{1/2}\left(\widehat{F}_{n}^{OR}-F\right)\to_{w}\mathcal{Z}_{\infty}\
\ \text{on}\ \ L^{\infty}\left(\left[-\infty,+\infty\right]\right)$ (3.75)

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

$\displaystyle\tau(x,y)=\text{Cov}\left(\mathcal{Z}_{\infty}(x),\mathcal{Z}_{\infty}(y)\right),\
\ \ x,y\in\left(-\infty,+\infty\right)$ (3.76)

We also define the ALASSO based residual empirical process as below

$\displaystyle\mathcal{Z}_{n}(x)=\sqrt{x}\left[\widehat{F}_{n}(x)-F(x)\right],\
\ x\in\left(-\infty,+\infty\right).$ (3.77)

###### Proof.

We denote with

$\displaystyle\Gamma_{n}\left(t,\mathbf{u}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\mathbf{1}\left(F\left(\epsilon_{i}-n^{1/2}\mathbf{x}_{i}^{\prime}\mathbf{u}\right)\leq
t\right)\right],\ \ t\in[0,1],\ \text{and}\ \mathbf{u}\in\mathbb{R}^{p}.$
(3.78)

Next, we define the random vector
$\mathbf{W}_{n}=n^{-1/2}\sum_{i=1}^{n}\mathbf{x}_{i}^{\prime}\epsilon_{i}$ and
write
$\mathbf{W}_{n}=\left(\mathbf{W}_{n}^{(1)\prime},\mathbf{W}_{n}^{(2)\prime}\right)$.

∎

####  3.6.2. 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:

##### Regression in the iid case

Under the usual assumption that the $\epsilon_{i}$ are iid and subGaussian,

$\displaystyle\forall\
s,\mathbb{E}\left[\text{exp}\left(s\epsilon_{i}^{2}\right)\right]\leq\text{exp}\left(\frac{s^{2}\sigma^{2}}{2}\right)$
(3.79)

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

$\displaystyle\mathbb{P}\left(\left|\frac{1}{n}\sum_{i=1}^{n}W_{i}^{(j)}\right|\geq\frac{t}{\sqrt{n}}\right)\leq\psi(t)=\text{exp}\left(-\frac{t^{2}}{2\sigma^{2}}\right).$
(3.80)

###  3.7. Lasso shrinkage for Predictive Regression

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

$\displaystyle\widehat{\beta}^{OLS}=\underset{\beta}{\text{arg
min}}\left\lVert y-X\beta\right\rVert^{2},$ (3.81)

The asymptotic behaviour of the OLS estimator has been extensively examined in
the time series econometrics literature (see, Koo et al.,  (2020), Lin and Tu,
(2020), Yousuf and Ng,  (2021), Lee et al.,  (2022)). A relevant framework for
lasso time series regression is proposed by Adamek et al.,  (2023).

In particular, the following functional central limit theorem holds

$\displaystyle\frac{1}{\sqrt{n}}\sum_{j=1}^{\lfloor
nr\rfloor}\begin{pmatrix}e_{j.}^{\prime}\\\
u_{j}\end{pmatrix}\Rightarrow\begin{pmatrix}B_{e}(r)\\\
B_{u}(r)\end{pmatrix}\equiv\text{BM}\left(\Sigma\right).$ (3.82)

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

$\displaystyle\frac{X^{\prime}u}{n}\Rightarrow\zeta:=\int_{0}^{1}B_{e}(r)dB_{u^{+}}(r)+\int_{0}^{1}B_{e}(r)\Sigma^{\prime}_{eu}\Sigma_{ee}^{-1}dB_{e}(r)^{\prime}$
(3.83)

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

$\displaystyle
n\left(\hat{\beta}^{OLS}-\beta_{n}^{*}\right)=\left(\frac{X^{\prime}X}{n^{2}}\right)^{-1}\frac{X^{\prime}u}{n}\Rightarrow\Omega^{-1}\zeta,$
(3.84)

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

$\displaystyle
n^{\delta_{j}}\left(\widehat{\beta}^{OLS}-\beta_{jn}^{*}\right)=n^{\delta_{j}}\widehat{\beta}^{OLS}-\beta_{j}^{0*}=\mathcal{O}_{p}\left(n^{\delta_{j}-1}\right)=o_{p}(1)\
\ \text{for all}\ j\leq p.$ (3.85)

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

$\displaystyle\widehat{\beta}^{\text{ora}}_{M^{*}}=\underset{\beta}{\text{arg
min}}\left\lVert y-\sum_{j\in M^{*}}x_{j}\beta_{j}\right\rVert^{2}$ (3.86)

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

$\displaystyle
n\left(\widehat{\beta}^{\text{ora}}-\beta_{n}^{*}\right)_{M^{*}}\Rightarrow\Omega_{M^{*}}\zeta_{M^{*}},$
(3.87)

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^{*}}$.

###  3.8. Lasso estimation with a structural break

Relevant studies with frameworks for lasso estimation robust to the presence
of structural breaks in the covariates include Chan et al.,  (2014), Cheng et
al.,  (2016), Qian and Su,  (2016) and Smith et al.,  (2019). 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}$.

$\displaystyle f_{(\alpha,\uptau)}(x,q)$
$\displaystyle:=x^{\prime}\beta+x^{\prime}\delta\mathbf{1}\big{\\{}q<\uptau\big{\\}}$
(3.88) $\displaystyle f_{0}(x,q)$
$\displaystyle:=x^{\prime}\beta_{0}+x^{\prime}\delta_{0}\mathbf{1}\big{\\{}q<\uptau_{0}\big{\\}}$
(3.89) $\displaystyle\hat{f}_{(\alpha,\uptau)}(x,q)$
$\displaystyle:=x^{\prime}\hat{\beta}+x^{\prime}\hat{\delta}\mathbf{1}\big{\\{}q<\hat{\uptau}\big{\\}}$
(3.90)

Then, we define the prediction risk as below

$\displaystyle\left\lVert\hat{f}-f_{0}\right\rVert_{n}:=\left[\frac{1}{n}\sum_{i=1}^{n}\bigg{(}\hat{f}(X_{i},Q_{i})-f_{0}(X_{i},Q_{i})\bigg{)}^{2}\right]^{1/2}.$
(3.91)

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

$\displaystyle S_{n}(\alpha,\uptau)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\bigg{(}Y_{i}-X_{i}^{\prime}\beta-
X_{i}^{\prime}\mathbf{1}\left\\{Q_{i}<\uptau\right\\}\bigg{)}^{2}\equiv\left\lVert\boldsymbol{y}-\boldsymbol{X}(\uptau)\alpha\right\rVert_{n}^{2}.$
(3.92)

######  Remark 13.

The examples discussed in the previous sections are concerned with the use of
the shrinkage approach under the assumption of sparsity in high dimensional
time series regression model for estimation purposes (see, also Basu and
Michailidis,  (2015) and García-Donato and Paulo,  (2022)). Consider that we
have $p-$variate time series then we can quantify the stability of the process
based on the largest eigenvalue of the transition matrix. Since the largest
eigenvalue is bounded, then a condition to roll out the cross-sectional and
serial correlation on the data is ensured. The literature on high-dimensional
time series modelling allows for dominant cross-sectional autocorrelations to
be accounted for, using common factors. In particular, Cho et al.,  (2023)
propose the factor-adjusted VAR model, in which case the time series is
decomposed into two latent components. A relevant framework for robust
estimation of high-dimensional VAR Models is given by Wang and Tsay,  (2023).
The frameworks mentioned in Section 3.8. are suitable for statistical
estimation in high dimensional regression models under the presence of
structural breaks555Dr. Haeran Cho (University of Bristol) gave a seminar
titled: ”High-dimensional time series segmentation via factor-adjusted VAR
modelling”, at the S3RI Departmental Seminar Series, University of Southampton
on the 24th of March 2022.
. One may be interested to develop an econometric framework for structural
break detection within these data-dependent model specifications such as as
the piecewise stationary factor-adjusted VAR models (see, Breitung and
Eickmeier,  (2011)). An alternative algorithmic procedure666Professor Degui Li
(University of York) gave a seminar titled: ”Detection of Multiple Structural
Breaks in Large Covariance Matrices”, at the S3RI Departmental Seminar Series,
University of Southampton on the 2nd of December 2021. for structural change
detection in high dimensional time series models is proposed by Li et al.,
(2021).

###  3.9. Classical Shrinkage Type Estimation Approach

In this section, we present the Score function approach proposed in the
framework of Sen and Saleh,  (1987), 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

$\displaystyle\boldsymbol{M}_{T}\left(\boldsymbol{b}\right)=\big{(}\boldsymbol{M}_{T1}\left(\boldsymbol{b}\right),...,\boldsymbol{M}_{Tp}\left(\boldsymbol{b}\right)\big{)}^{\prime}=\sum_{t=1}^{T}\boldsymbol{a}_{t}\psi\big{(}\boldsymbol{y}_{t}-\boldsymbol{a}_{t}^{\prime}\boldsymbol{b}\big{)},\
\ \boldsymbol{b}\in\mathbb{R}^{p}.$ (3.93)

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

$\displaystyle\boldsymbol{M}_{T(1)}\big{(}\boldsymbol{b}_{1},\boldsymbol{0}\big{)}=\boldsymbol{0}.$
(3.94)

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

$\displaystyle\hat{\boldsymbol{M}}_{T(2)}=\boldsymbol{M}_{T(2)}\big{(}\hat{\boldsymbol{\beta}}_{1},\boldsymbol{0}\big{)},$
(3.95)

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, Sen and Saleh,  (1987)). Furthermore, denote
with

$\displaystyle\mathcal{S}_{T}^{2}$
$\displaystyle=T^{-1}\sum_{t=1}^{T}\psi^{2}\bigg{(}\boldsymbol{y}_{t}-\hat{\boldsymbol{\beta}}_{1T}^{\prime}\boldsymbol{x}_{(j)t-1}\bigg{)},\
\ \ \
\boldsymbol{x}_{t-1}=\left(\boldsymbol{x}_{(1)t-1},\boldsymbol{x}_{(2)t-1}\right),$
(3.96) $\displaystyle\boldsymbol{C}_{ii.j}$
$\displaystyle=\boldsymbol{C}_{ii}-\boldsymbol{C}_{ij}\boldsymbol{C}_{jj}^{-1}\boldsymbol{C}_{ji},\
\ \text{for}\ \ i\neq j=\left\\{1,2\right\\}.$ (3.97)

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

$\displaystyle\mathcal{L}_{T}=S_{T}^{-2}\bigg{\\{}\hat{\boldsymbol{M}}_{T(2)}\boldsymbol{C}_{22.1}\hat{\boldsymbol{M}}_{T(2)}\bigg{\\}}.$
(3.98)

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:

$\displaystyle\mathcal{S}_{T}^{2}=T^{-1}\sum_{t=1}^{T}\psi^{2}\bigg{(}\boldsymbol{y}_{t}-\hat{\boldsymbol{\beta}}_{1T,ols}^{\prime}\boldsymbol{x}_{(j)t-1}\bigg{)},\
\ \ \
\boldsymbol{x}_{t-1}=\left(\boldsymbol{x}_{(1)t-1},\boldsymbol{x}_{(2)t-1}\right),\
j\in\left\\{1,2\right\\},$ (3.99)

and

$\displaystyle\widetilde{\mathcal{S}}_{T}^{2}=T^{-1}\sum_{t=1}^{T}\psi^{2}\bigg{(}\boldsymbol{y}_{t}-\hat{\boldsymbol{\beta}}_{1T,ivx}^{\prime}\boldsymbol{x}_{(j)t-1}\bigg{)},\
\ \ \
\boldsymbol{x}_{t-1}=\left(\boldsymbol{x}_{(1)t-1},\boldsymbol{x}_{(2)t-1}\right),\
j\in\left\\{1,2\right\\}.$ (3.100)

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, Sen and Saleh,
(1987)).

####  3.9.1. Limiting Distributional Risk

In the multivariate location model Sen and Saleh,  (1985) 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,

$\displaystyle\mathcal{K}_{T}:\boldsymbol{\beta}_{2}=\boldsymbol{\beta}_{2(T)}=T^{-1/2}\boldsymbol{\xi},\
\ \
\boldsymbol{\xi}=\big{(}\xi_{p_{1}+1},...,\xi_{p}\big{)}^{\prime}\in\mathbb{R}^{p_{2}},$
(3.101)

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

$\displaystyle\mathcal{R}\big{(}\boldsymbol{\beta}^{\star}_{1T};\boldsymbol{W}\big{)}=\mathsf{trace}\left\\{\boldsymbol{W}\int_{\mathbb{R}^{p}}...\int_{\mathbb{R}^{p}}\boldsymbol{\mathsf{u}}\boldsymbol{\mathsf{u}}^{\prime}d\boldsymbol{G}^{\star}(\boldsymbol{\mathsf{u}})\right\\}\equiv\mathsf{trace}\left\\{\boldsymbol{W}\boldsymbol{V}^{\star}\right\\},$
(3.102)

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

####  3.9.2. Asymptotic representation of estimators

We follow the framework presented in Mukherjee,  (1998).

##### 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\lVert s_{\bullet 1}\right\rVert,...,\left\lVert
s_{\bullet p}\right\rVert\big{)}$, where $\left\lVert s_{\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,

$\displaystyle\boldsymbol{R}_{n}:=\boldsymbol{D}_{x}^{-1}\boldsymbol{X}^{\prime}\boldsymbol{X}\boldsymbol{D}_{x}^{-1}\overset{p}{\to}\boldsymbol{R}:=\boldsymbol{D}_{c}^{-1}\boldsymbol{C}\boldsymbol{D}_{c}^{-1}$
(3.103)

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 Sen and Saleh,  (1987), 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

$\displaystyle\hat{\beta}_{s}=\hat{\beta}_{1}+\left(1-\frac{\kappa}{\mathcal{L}_{T}}\right)\left(\tilde{\beta}_{1}-\hat{\beta}_{1}\right),c>0.$
(3.104)

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

$\displaystyle\boldsymbol{H}:=\big{[}\boldsymbol{0}\ \vdots\
\boldsymbol{R}_{22.1}\big{]}\boldsymbol{R}^{-1}\boldsymbol{E}\boldsymbol{R}^{-1}\big{[}\boldsymbol{0}\
\vdots\ \boldsymbol{R}_{22.1}\big{]}$ (3.105)

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

$\displaystyle\boldsymbol{H}_{n}:=\big{[}\boldsymbol{0}\ \vdots\
\boldsymbol{R}_{n22.1}\big{]}\boldsymbol{R}_{n}^{-1}\boldsymbol{E}_{n}\boldsymbol{R}_{n}^{-1}\big{[}\boldsymbol{0}\
\vdots\ \boldsymbol{R}_{n22.1}\big{]}$ (3.106)

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

$\displaystyle\mathcal{L}_{n}:=\left\lVert\boldsymbol{H}_{n}^{-1/2}\boldsymbol{B}_{x2}^{-1}\boldsymbol{X}_{2}^{\prime}\big{(}\boldsymbol{Y}-\boldsymbol{X}_{1}\hat{\boldsymbol{\beta}}_{1}\big{)}\right\rVert^{2}$
(3.107)

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.

###  3.10. Lasso Inference for High Dimensional Time Series under NED

Following the framework proposed by Adamek et al.,  (2023), consider the
linear time series regression

$\displaystyle y_{t}=\boldsymbol{x}_{t}^{\prime}\boldsymbol{\beta}^{0}+u_{t},\
\ \ \ t=1,...,T,$ (3.108)

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$.

######  Assumption 4.

Let
$\boldsymbol{z}_{t}=\left(\boldsymbol{x}_{t}^{\prime},u_{t}\right)^{\prime}$
and let there exists some constants $\bar{m}>m>2$ and
$d\geq\mathsf{max}\left\\{1,\frac{\frac{\bar{m}}{m-1}}{\bar{m}-2}\right\\}$,

$\displaystyle\mathbb{E}\big{[}\boldsymbol{z}_{t}\big{]}=\boldsymbol{0},\ \ \
\mathbb{E}\big{[}\boldsymbol{x}_{t}u_{t}\big{]}=\boldsymbol{0}\ \ \
\text{and}\ \ \ \underset{1\leq j\leq N+1,1\leq t\leq T}{\mathsf{max}}\
\mathbb{E}|z_{j,t}|^{2\bar{m}}\leq C.$ (3.109)

######  Assumption 5.

Let $\boldsymbol{s}_{T,t}$ denote a $k(T)-$dimensional triangular array that
is $\alpha-$mixing of size

$\displaystyle\frac{d}{\left(\frac{1}{m}-\frac{1}{\bar{m}}\right)},\ \ \
\text{with}\ \ \sigma-\text{field}\
\mathcal{F}_{t}^{s}:=\sigma\left\\{\boldsymbol{s}_{T,t},\boldsymbol{s}_{T,t-1},...\right\\}$
(3.110)

such that $\boldsymbol{z}_{t}$ is $\mathcal{F}_{t}^{s}-$measurable. The
process $\left\\{z_{j,t}\right\\}$ is $L_{2m}-$near-epoch dependent (NED) of
size $-d$ on $\boldsymbol{s}_{T,t}$ with positive bounded NED constants,
uniformly over $j\in\left\\{1,...,N+1\right\\}$.

######  Remark 14.

The first assumption above ensures that the error terms are contemporaneously
uncorrelated with each of the regressors, and that the process has finite and
constant unconditional moments. The second assumption above implies that
$\boldsymbol{s}_{T,t}$ is an underlying shock process driving the regressors
and errors in $\boldsymbol{z}_{t}$, where we assume $\boldsymbol{z}_{t}$ to
depend almost nearly entirely on the "near epoch" of $\boldsymbol{s}_{T,t}$.
In other words, the near epoch dependence of the vector $\boldsymbol{z}_{t}$
can be interpreted as $\boldsymbol{z}_{t}$ being "approximately" mixing, in
the sense that it can be well-approximated by a mixing process. The NED
framework allows for general forms of dependence and that are often
encountered in econometrics applications including, but not limited to, strong
mixing processes, linear processes including ARMA models, various types of
stochastic volatility and GARCH models as well as nonlinear processes (see,
Davidson,  (2002)).

The framework proposed by Adamek et al.,  (2023) is the first to utilize the
NED framework for establishing uniformly valid high-dimensional inference. For
instance, Wong et al.,  (2020) 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, Masini et al.,  (2022) 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 Adamek et al.,  (2023) 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 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.

######  Example 13 (ARCL Model with GARCH errors, see Adamek et al.,  (2023)).

Consider the autoregressive distributed lag (ARDL) model with GARCH errors
such that

$\displaystyle y_{t}$
$\displaystyle=\sum_{i=1}^{p}\rho_{i}y_{t-i}+\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}\boldsymbol{w}_{t-j}+u_{t}=\boldsymbol{x}_{t}^{\prime}\boldsymbol{\beta}^{0}+u_{t}$
(3.111) $\displaystyle u_{t}$ $\displaystyle=\sqrt{h}_{t}\varepsilon_{t},\ \ \
\ \varepsilon_{t}\sim\mathcal{N}(0,1),$ (3.112) $\displaystyle h_{t}$
$\displaystyle=\beta_{0}+\beta_{1}h_{t-1}+\beta_{2}u_{t-1}^{2},$ (3.113)

where the roots of the lag polynomial $\rho(z)=1-\sum{j=1}^{p}\rho_{j}z^{j}$
are outside the unit circle. Then, $u_{t}$ is a strictly stationary
geometrically $\beta-$mixing process. Moreover, we assume that the vector of
exogenous variables $\boldsymbol{w}_{t}$ is stationary and geometrically
$\beta-$mixing as well with $2\bar{m}$ finite moments.

Given the invertability of the lag polynomial, we may then write

$\displaystyle y_{t}=\rho^{-1}(L)v_{t},\ \ \ \
v_{t}=\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}\boldsymbol{w}_{t-j}+u_{t},$
(3.114)

and the inverse lag polynomial $\rho^{-1}(z)$ has geometrically decaying
coefficients. Then, it follows directly that $y_{t}$ is NED on $v_{t}$, where
$v_{t}$ is strong mixing of size $-\infty$ as its components are geometrically
$\beta-$mixing, and the sum inherits the mixing properties.

######  Example 14 (Equation-by-Equation VAR, see Adamek et al.,  (2023)).

Consider the VAR model below

$\displaystyle\boldsymbol{y}_{t}=\sum_{j=1}^{p}\boldsymbol{\Phi}_{j}\boldsymbol{y}_{t-j}+\boldsymbol{u}_{t},$
(3.115)

where $\boldsymbol{y}_{t}$ is a $K\times 1$ vector of dependent variables and
$\mathbb{E}|u_{t}|^{2\bar{m}}\leq C$ and the $K\times K$ matrices
$\boldsymbol{\Phi}_{i}$ satisfy appropriate stationarity and summability
conditions. The equivalent equation-by-equation representation

$\displaystyle
y_{k,t}=\sum_{j=1}^{p}\big{[}\Phi_{k,1,j},...,\Phi_{k,K,j}\big{]}\boldsymbol{y}_{t-j}+u_{k,t}=\big{[}\boldsymbol{y}_{t-1}^{\prime},...,\boldsymbol{y}_{t-p}^{\prime}\big{]}\boldsymbol{\beta}_{k}+u_{k,t},\
\ \ k\in\left\\{1,...,K\right\\}.$ (3.116)

Therefore, assuming a well-specified model with
$\mathbb{E}\big{[}\boldsymbol{u}_{t}|\boldsymbol{y}_{t-1},...,\boldsymbol{y}_{t-p}\big{]}=\boldsymbol{0}$,
then the NED conditions are then satisfied trivially.

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.

######  Example 15 (Misspecified AR Model, see Adamek et al.,  (2023)).

Consider an autoregressive AR model of order 2 as

$\displaystyle y_{t}=\rho_{1}y_{t-1}+\rho_{2}y_{t-2}+v_{t},\ \ \
v_{t}\sim\mathcal{N}(0,1).$ (3.117)

where the roots of the characteristic polynomial $(1-\rho_{1}L-\rho_{2}L^{2})$
are outside the unit circle. Moreover, define the misspecified model
$y_{t}=\tilde{\rho}y_{t-1}+u_{t}$, where

$\displaystyle\tilde{\rho}=\underset{\rho}{\mathsf{argmin}}\
\mathbb{E}\left[\left(y_{t}-\rho
y_{t-1}\right)^{2}\right]=\frac{\mathbb{E}\big{[}y_{t}y_{t-1}\big{]}}{\mathbb{E}\left[y_{t-1}^{2}\right]}=\frac{\rho_{1}}{1-\rho_{2}}$
(3.118)

and $u_{t}$ is autocorrelated.

Therefore, an m.d.s assumption would be inappropriate in this case, since

$\displaystyle\mathbb{E}\big{[}u_{t}|\sigma\left\\{y_{t-1},y_{t-2},....\right\\}\big{]}=\mathbb{E}\big{[}y_{t}=\tilde{\rho}y_{t-1}+u_{t}|\sigma\left\\{y_{t-1},y_{t-2},....\right\\}\big{]}=-\frac{\rho_{1}\rho_{2}}{1-\rho_{2}}y_{t-1}+\rho_{2}y_{t-2}\neq
0.$ (3.119)

On the other hand, it can be shown that $\big{(}y_{t-1},u_{t}\big{)}^{\prime}$
satisfies the NED-condition by considering the moving average representation
of $y_{t}$ and thus by extension, of $u_{t}=y_{t}-\tilde{\rho}y_{t-1}$. In
addition, since the coefficients are geometrically decaying, $u_{t}$ is thus
clearly NED on $v_{t}$ and therefore the NED-condition is satisfied.

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.

######  Assumption 6.

For some $0\leq r\leq 1$ and sparsity level $s_{r}$, define the
$N-$dimensional sparse compact parameter space such that

$\displaystyle\mathcal{B}_{N}(r,s_{r}):=\big{\\{}\boldsymbol{\beta}\in\mathbb{R}^{N}:\left\lVert\boldsymbol{\beta}\right\rVert_{r}^{r}\leq
s_{r},\left\lVert\boldsymbol{\beta}\right\rVert_{\infty}\leq C,\exists\
C<\infty\big{\\}},$ (3.120)

and assume that $\boldsymbol{\beta}^{0}\in\mathcal{B}_{N}(r,s_{r})$

######  Example 16 (Infinite Order AR,see Adamek et al.,  (2023)).

Consider an infinite order autoregressive model as below

$\displaystyle y_{t}=\sum_{j=1}^{\infty}\rho_{j}y_{t-j}+\varepsilon_{t},$
(3.121)

where $\varepsilon_{t}$ is a stationary m.d.s with sufficient moments
existing, and the lag polynomial
$\left(1-\sum_{j=1}^{\infty}\rho_{j}L^{j}\right)$ is invertible and satisfies
the summability condition $\sum_{j=1}^{\infty}j^{\alpha}|\rho_{j}|<\infty$ for
some $\alpha\geq 0$.

####  3.10.1. 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,
Adamek et al.,  (2023)). 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

$\displaystyle\mathcal{H}:=\left\\{j:\sum_{p=1}^{P}|r_{N,p,j}|>0\right\\}$
(3.122)

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

$\displaystyle\boldsymbol{\Omega}_{N,T}=\mathbb{E}\left[\frac{1}{T}\left(\sum_{t=1}^{T}\boldsymbol{w}_{t}\right)\left(\sum_{t=1}^{T}\boldsymbol{w}_{t}^{\prime}\right)\right],$
(3.123)

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

$\displaystyle\boldsymbol{\Omega}_{N,T}$
$\displaystyle\equiv\boldsymbol{\Xi}(0)+\sum_{\ell=1}^{T-1}\big{[}\boldsymbol{\Xi}(\ell)+\boldsymbol{\Xi}^{\prime}(\ell)\big{]},$
(3.124) $\displaystyle\boldsymbol{\Xi}(\ell)$
$\displaystyle\equiv\frac{1}{T}\sum_{t=\ell+1}^{T}\mathbb{E}\left[\boldsymbol{w}_{t}\boldsymbol{w}_{t-\ell}^{\prime}\right].$
(3.125)

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

$\displaystyle\boldsymbol{\Omega}$
$\displaystyle\equiv\widehat{\boldsymbol{\Xi}}(0)+\sum_{\ell=1}^{Q_{T}-1}K\left(\frac{\ell}{Q_{T}}\right)\big{[}\boldsymbol{\Xi}(\ell)+\boldsymbol{\Xi}^{\prime}(\ell)\big{]},$
(3.126) $\displaystyle\boldsymbol{\Xi}(\ell)$
$\displaystyle=\frac{1}{T-\ell}\sum_{t=\ell+1}^{T}\hat{\boldsymbol{w}}_{t}\hat{\boldsymbol{w}}_{t}^{\prime}$
(3.127)

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)$.

####  3.10.2. 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, Adamek et al.,  (2023))

$\displaystyle\hat{\boldsymbol{\beta}}:=\underset{\beta\in\mathbb{R}^{N}}{\mathsf{argmin}}\
\left\\{\frac{1}{T}\left\lVert\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right\rVert_{2}^{2}+2\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{1}\right\\}.$
(3.128)

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

$\displaystyle\mathcal{E}_{T}(z):=\left\\{\underset{j\leq N,\ell\leq
T}{\mathsf{max}}\ \left|\sum_{t=1}^{\ell}u_{t}x_{j,t}\leq z\right|\right\\}$
(3.129)

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

####  3.10.3. Uniformly Valid Inference via the Disparsified Lasso

Following the framework of Adamek et al.,  (2023), we consider an application
of the disparsified Lasso.

######  Example 17 (Sparse Factor Models).

Consider the factor model

$\displaystyle y_{t}$
$\displaystyle=\boldsymbol{\beta}^{0\prime}\boldsymbol{x}_{t}+u_{t},\ \ \
u_{t}\sim\mathcal{N}(0,1)$ (3.130) $\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{\Lambda}\boldsymbol{f}_{t}+\boldsymbol{v}_{t},\ \ \
\boldsymbol{v}_{t}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{v}),\ \
\ \boldsymbol{f}_{t}\sim(\boldsymbol{0},\boldsymbol{\Sigma}_{f}).$ (3.131)

Moreover, we assume that
$\boldsymbol{\Sigma}=\boldsymbol{\Lambda}\boldsymbol{\Sigma}_{f}\boldsymbol{\Lambda}^{\prime}+\boldsymbol{\Sigma}_{v}$.

######  Example 18 (Sparse VAR(1)).

Consider a stationary VAR(1) model for
$\boldsymbol{z}_{t}=\left(y_{t},\boldsymbol{x}_{t}^{\prime}\right)^{\prime}$

$\displaystyle\boldsymbol{z}_{t}=\boldsymbol{\Phi}\boldsymbol{z}_{t-1}+\boldsymbol{u}_{t},\
\ \
\mathbb{E}\left[\boldsymbol{u}_{t}\boldsymbol{u}_{t}^{\prime}\right]:=\boldsymbol{\Omega},\
\
\mathbb{E}\left[\boldsymbol{u}_{t}\boldsymbol{u}_{t}^{\prime}\right]=\boldsymbol{0},\
\ \forall\ \ell\neq 0.$ (3.132)

Then, the regression
$y_{t}=\boldsymbol{\phi}_{1}\boldsymbol{z}_{t-1}+u_{1,t}$, where
$\boldsymbol{\phi}_{j}$ is the $j-$th row of $\boldsymbol{\Phi}$, i.e., the
first line of the VAR.

######  Example 19.

Consider the population nodewise regressions defined by the linear projections
below

$\displaystyle
x_{j,t}=\boldsymbol{x}^{\prime}_{-j,t}\boldsymbol{\gamma}_{j}^{0}+v_{j,t},\ \
\ \
\boldsymbol{\gamma}_{j}^{0}:=\underset{\boldsymbol{\gamma}}{\mathsf{argmin}}\
\mathbb{E}\left[\frac{1}{T}\sum_{t=1}^{T}\left(x_{j,t}-\boldsymbol{x}_{-j,t}^{\prime}\boldsymbol{\gamma}\right)^{2}\right],\
\ \ \ j=1,....,N.$ (3.133)

with
$\tau_{j}^{2}:=\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[v_{j,t}^{2}\right]$.
Specifically, if we let $\boldsymbol{\Phi}\phi\boldsymbol{I}$, with
$|\phi|<1$, and let $\boldsymbol{\Omega}$ have a Toeplitz structure, which
implies that the elements of the matrix are defined such that
$\omega_{i,j}=\rho^{|i-j|},|\rho|<1$. Therefore, the nodewise regression
parameter $\boldsymbol{\gamma}_{j}^{0}$ is only weakly sparse, in that it
contains no zeroes, but its entries follow a geometrically decaying
pattern777Notice that in high dimensional settings relevant conditions that
capture features such as sub-Gaussianity, heavy tailed observations with slowy
decaying temporal dependence are essential (see, Baek et al.,  (2021) and
Guillaumin et al.,  (2022). Dr. Adam Sykulski gave a seminar with title: ”The
Debiased Spatial Whittle Likelihood”, at the S3RI Departmental Seminar Series
at the University of Southampton on the 5th of May 2022., meaning that
$\underset{j}{\mathsf{max}}\left\lVert\boldsymbol{\gamma}_{j}^{0}\right\rVert_{r}^{r}\leq
C$.

##  4\. 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 screening888Note that the covariate screening approach is widely
used in the applied statistics literature, see Aschard et al.,  (2017).
approach as discussed by Fan and Lv,  (2008), Wang,  (2009), Ke et al.,
(2014) and Barut et al.,  (2016) as well as He et al.,  (2013) and Kong et
al.,  (2019) 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 Belloni and
Chernozhukov,  (2011), Belloni et al.,  (2014), Wang et al.,  (2012),
Chernozhukov et al.,  (2015), Fan et al.,  (2017) as well as Sun et al.,
(2015) and Wang and Leng,  (2016) 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. Lockhart et al.,  (2014) propose a test for detemining
statistical significance with the Lasso, while Chudik et al.,  (2018) propose
a methodology for selecting covariates in a high-dimensional environment using
the multiple testing approach. Furthermore, Battey et al.,  (2018) 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, Battey
et al.,  (2018) 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.

######  Lemma 10 (Battey et al.,  (2018)).

Assume that
$\Sigma=\mathbb{E}\big{(}\boldsymbol{X}_{i}\boldsymbol{X}_{i}^{\top}\big{)}$
satisfies
$C_{\mathsf{min}}<\lambda_{\mathsf{min}}(\Sigma)\leq\lambda_{\mathsf{min}}(\Sigma)\leq
C_{\mathsf{max}}$ as well as
$\left\lVert\Sigma^{-1/2}\boldsymbol{X}_{1}\right\rVert_{\psi_{2}}=\kappa$,
then it holds that

$\displaystyle\mathbb{P}\left(\underset{j=1,...,k}{\mathsf{max}}\left\lVert
M^{(j)}\widehat{\Sigma}^{(j)}-I\right\rVert_{\mathsf{max}}\leq\alpha\sqrt{\frac{\mathsf{log}d}{n}}\right)\geq
1-2kd^{-c_{2}},\ \ \
c_{2}=\frac{\alpha^{2}C_{\mathsf{min}}}{24e^{2}\kappa^{4}C_{\mathsf{max}}}-2.$
(4.1)

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 Chernozhukov et al.,  (2014) and Zhang and Wu,
(2017)) and consistent model selection techniques robust to the presence of
heteroscedasticity (see, Halunga et al.,  (2017)) 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

$\displaystyle X_{t}=A_{1}X_{t-1}+...+A_{d}X_{t-d}+e_{t},\ \ \
e_{t}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{e})$ (4.2)

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, Michailidis and d’Alché Buc,  (2013), Basu and Michailidis,
(2015) and Basu and Rao,  (2021)).

###  4.1. Ultra-high dimensionality under dependence

According to Yousuf,  (2018), 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

$\displaystyle\varepsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j},\ \ \
e_{i}\sim\textit{i.i.d},$

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

$\displaystyle\Delta_{0,q}(\varepsilon)(\epsilon_{\ell})=\left\lVert
e_{0}-e_{0}^{*}\right\rVert_{q}\sum_{\ell=0}^{\infty}|f_{\ell}|$ (4.3)

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.

######  Theorem 12 (Theorem 2 in Yousuf,  (2018)).

Define with $\alpha=\frac{2}{1+4\tilde{\alpha}_{x}}$, then it holds that for
any $c_{2}>0$

$\displaystyle\mathbb{P}\left(\underset{j\leq
p_{n}}{\mathsf{max}}\left|\hat{\rho}_{j}-\rho_{j}\right|>c_{2}n^{-k}\right)\leq\mathcal{O}_{p}\left(s_{n}p_{n}\mathsf{exp}\left\\{-\frac{n^{1/2-\kappa}}{\nu_{x}^{2}s_{n}}\right\\}^{\tilde{\alpha}}\right)+\mathcal{O}_{p}\left(p_{n}\mathsf{exp}\left\\{-\frac{n^{1/2-\kappa}}{\nu_{x}\nu_{\epsilon}}\right\\}^{\tilde{\alpha}^{\prime}}\right)$
(4.4)

####  4.1.1. 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

$\displaystyle\hat{\boldsymbol{\rho}}=\left(\hat{\rho}_{1},...,\hat{\rho}_{p_{n}}\right),\
\ \text{where}\ \
\hat{\rho}_{j}=\left(\sum_{t=1}^{n}X_{tj}^{2}\right)^{-1}\left(\sum_{t=1}^{n}X_{tj}Y_{t}\right).$
(4.5)

Therefore, $\hat{\rho}_{j}$ is the OLS estimate of the linear projection of
$Y_{t}$ onto $X_{tj}$. Let

$\displaystyle\mathcal{M}_{*}=\big{\\{}1\leq i\leq p_{n}:\beta_{i}\neq
0\big{\\}}$ (4.6)

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

$\displaystyle\hat{\mathcal{M}}_{\gamma_{n}}=\big{\\{}1\leq i\leq
p_{n}:\left|\hat{\rho}_{i}\right|\geq\gamma_{n}\big{\\}}$ (4.7)

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

$\displaystyle\mathbb{P}\left(\mathcal{M}_{*}\subset\hat{\mathcal{M}}_{\gamma_{n}}\right)\to
1.$

####  4.1.2. Experimental design examples

##### Experimental Design I: Uncorrelated Features

Consider the model for the covariate process we have

$\displaystyle\boldsymbol{x}_{t}=\boldsymbol{A}_{1}\boldsymbol{x}_{t-1}+\boldsymbol{\eta}_{t}$
(4.8)

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.

##### 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

$\displaystyle\Sigma_{\eta}=\left\\{0.3^{|i-j|}\right\\}$ (4.9)

###### Proof.

Recall that it holds that
$\sum_{k=1}^{p_{n}}\mathbf{1}\left\\{\left|\beta_{k}\right|>0\right\\}=s_{n}$.
Hence, we obtain the following expression

$\displaystyle\mathbb{P}\left(\left|S_{1}-\mathbb{E}(S_{1})\right|>\frac{c_{2}n^{-\kappa}}{2}\right)\leq\
\sum_{k\in
M_{*}}\mathbb{P}\left(\left|\frac{X_{tj}\left(X_{tk}\beta_{k}\right)}{n}-\beta_{k}\mathbb{E}\left[X_{tj}X_{tk}\right]\right|>\frac{c_{2}n^{-\kappa}}{2s_{n}}\right)$
(4.10)

Moreover, it holds that $\left\lVert
X_{ij}\right\rVert_{r}\leq\Delta_{0,r}(\boldsymbol{X}_{j})\leq\Phi_{0,r}(\boldsymbol{x})$.
Using this we compute the cumulative functional dependence measure of
$X_{tk}X_{tj}$ such that

$\displaystyle\sum_{t=m}^{+\infty}\left\lVert
X_{tj}X_{tk}-X_{tj}^{*}X_{tk}^{*}\right\rVert_{r/2}$
$\displaystyle\leq\sum_{t=m}^{\infty}\big{(}\left\lVert
X_{tj}\right\rVert_{r}\left\lVert X_{tk}-X_{tk}^{*}\right\rVert_{r}\left\lVert
X_{tj}-X_{tj}^{*}\right\rVert_{r}\big{)}$ $\displaystyle\leq
2\Phi_{0,r}(x)\Phi_{m,r}(x)=\mathcal{O}_{p}\left(m^{-\alpha_{x}}\right).$

Therefore, we obtain that

$\displaystyle\underset{m}{\mathsf{sup}}\
(m+1)^{\alpha_{x}}\sum_{t=m}^{\infty}\left\lVert
X_{tj}\right\rVert_{r}\left\lVert X_{tk}-X_{tk}^{*}\right\rVert_{r/2}\leq
2K_{x,r}^{2}.$ (4.11)

Using the above results we get the following probability bound

$\displaystyle\mathbb{P}\left(\left|S_{1}-\mathbb{E}\left[S_{1}\right]\right|>\frac{c_{2}n^{-\kappa}}{2}\right)\leq
Cs_{n}\times\left\\{\frac{n^{\omega}K_{x,r}^{r}}{(n/s_{n})^{r/2-rk/2}}+\mathsf{exp}\left(-\frac{n^{1-2\kappa}}{s_{n}^{2}K_{x,r}^{4}}\right)\right\\}$
(4.12)

Hence, by our choice of $\gamma_{n}$, we obtain that
$\mathbb{P}\left(\mathcal{M}_{*}\subset\hat{\mathcal{M}}_{\gamma_{n}}\right)>\mathbb{P}\left(\mathcal{A}_{n}\right)$.
In other words, this allow us to obtain a probability bound with respect to
the spectral radius of the corresponding covariance matrix

$\displaystyle\sum_{k=1}^{p_{n}}\rho_{k}^{2}=\mathcal{O}_{p}\left(\lambda_{\mathsf{max}}(\Sigma)\right).$
(4.13)

Then, we extend the above result to the corresponding set defined by
$\mathcal{B}_{n}$ such that

$\displaystyle\mathcal{B}_{n}:=\left\\{\underset{k\leq
p_{n}}{\mathsf{max}}\left|\hat{\rho}_{k}-\rho_{k}\right|\leq
c_{4}n^{-\kappa}\right\\}$ (4.14)

Therefore, we obtain a measurability condition for the number of
$\left\\{k:\left|\hat{\rho}_{k}\right|>2c_{4}n^{-4}\right\\}$ cannot exceed
the number of $\left\\{k:\left|\rho_{k}\right|>c_{4}n^{-k}\right\\}$ which is
bounded by $\mathcal{O}_{p}\left(n^{2k}\lambda_{\mathsf{max}}(\Sigma)\right)$.
Thus, by choosing $c_{4}=c_{3}/2$ we can obtain the correct probability bound
such that

$\displaystyle\mathbb{P}\left(\left|\hat{\mathcal{M}}_{\gamma_{n}}\right|<\mathcal{O}_{p}\left(n^{2k}\lambda_{\mathsf{max}}(\Sigma)\right)\right)>\mathbb{P}\left(\mathcal{B}_{n}\right)$
(4.15)

We are interested to bound the following probability event

$\displaystyle\mathcal{A}_{1}:=\mathbb{P}\left\\{\left|\hat{\gamma}_{i,k}-\gamma_{i,k}\right|>\frac{cn^{-k}}{\ell_{n}}\right\\}$
(4.16)

where a probability bound for the event $\mathcal{A}_{1}$, can be obtained as
below:

$\displaystyle\mathcal{A}_{1}$
$\displaystyle\leq\mathbb{P}\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}\varepsilon_{t+|i|,k}-\mathbb{E}\left[\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}\varepsilon_{t+|i|,k}\right]\right|\right)$
$\displaystyle+\mathbb{P}\left(\left|\mathbb{E}\left[\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}\varepsilon_{t+|i|,k}\right]-\gamma_{i,k}\right|>cn^{-k}/4\ell_{n}\right)$
$\displaystyle+\mathbb{P}\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t+|i|,k}X_{t+|i|,k}\left(\frac{\sum_{j=1}^{n}X_{jk}\varepsilon_{j,k}/n}{\sum_{j=1}^{n}X_{jk}^{2}/n}\right)\right|>cn^{-k}/4\ell_{n}\right)$
$\displaystyle+\mathbb{P}\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}X_{tk}X_{t+|i|,k}\left(\frac{\sum_{j=1}^{n}X_{jk}\varepsilon_{j,k}/n}{\sum_{j=1}^{n}X_{jk}^{2}/n}\right)^{2}\right|>cn^{-k}/4\ell_{n}\right)$

Furthermore, we can determine the bias of the term as below

$\displaystyle\left|\mathbb{E}\left[\sum_{t=1}^{n-|i|}\frac{\varepsilon_{t,k}\varepsilon_{t+|i|,k}}{n}\right]-\gamma_{i,k}\right|\leq\frac{i\gamma_{i,k}}{n}.$
(4.17)

Moreover, it holds that

$\displaystyle\mathbb{P}$
$\displaystyle\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}X_{t+|i|,k}>M\right|\right)\leq\mathbb{P}\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}X_{t+|i|,k}-\mathbb{E}\left[\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}X_{t+|i|,k}\right]\right|>M-\left|\mathbb{E}\left[\frac{1}{n}\sum_{t=1}^{n-|i|}\varepsilon_{t,k}X_{t+|i|,k}\right]\right|\right)$
$\displaystyle M>\underset{k\leq p_{n}}{\mathsf{max}}\
\underset{i\leq\ell_{n}}{\mathsf{max}}\
2\left|\mathbb{E}\left(\varepsilon_{t,k}X_{t+|i|,k}\right)\right|+\varepsilon,\
\ \text{for some}\ \varepsilon>0.$ (4.18)

∎

###  4.2. Uniform Asymptotic Inference and Model Selection

This is another important topic which we cannot cover here extensively (see,
Tibshirani et al.,  (2018)).

####  4.2.1. Uniform-in-Submodel Bounds

Following the framework proposed by Kuchibhotla et al.,  (2021), 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

$\displaystyle\underset{M\in\mathcal{M}}{\text{sup}}\left\lVert\widehat{\beta}_{M}-\widehat{\beta}\right\rVert=o_{p}(1),\
\ \text{as}\ \ n\to\infty,\ \ \ \ \text{for some}\ \left\lVert.\right\rVert.$
(4.19)

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 Kuchibhotla et al.,  (2021), 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, Chudik et al.,  (2018)) as
well as the uniform inference post-model selection which includes among others
the studies of Belloni et al.,  (2016) and Farrell,  (2015)).

Here we focus on the framework proposed by Kuchibhotla et al.,  (2021) 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

$\displaystyle\hat{R}_{n}\left(\theta;M\right):=\frac{1}{n}\sum_{i=1}^{n}\left\\{Y_{i}-X_{i}^{\top}\left(M\right)\theta\right\\}^{2},\
\ \text{for}\ \ \theta\in\mathbb{R}^{|M|}.$ (4.20)

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

$\displaystyle\hat{R}_{n}\left(\theta;M\right)=\frac{1}{n}\sum_{i=1}^{n}Y_{i}^{2}-\frac{2}{n}\sum_{i=1}^{n}Y_{i}X_{i}^{\top}(M)\theta+\theta^{\top}\left(\frac{1}{n}\sum_{i=1}^{n}X_{i}(M)X_{i}^{\top}(M)\right)\theta.$
(4.21)

Furthermore, define the following quantities

$\displaystyle\hat{\Sigma}_{n}:=\frac{1}{n}\sum_{i=1}^{n}X_{i}X_{i}^{\top}\in\mathbb{R}^{p\times
p},\ \ \ \text{and}\ \ \
\hat{\Gamma}_{n}:=\frac{1}{n}\sum_{i=1}^{n}X_{i}Y_{i}\in\mathbb{R}^{p\times
p}.$ (4.22)

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

$\displaystyle\hat{\beta}_{n,M}=\underset{\theta\in\mathbb{R}^{|M|}}{\text{arg
min}}\
\hat{R}_{n}\left(\theta;M\right)=\underset{\theta\in\mathbb{R}^{|M|}}{\text{arg
min}}\left\\{\theta^{\top}\hat{\Sigma}_{n}(M)\theta-2\theta^{\top}\hat{\Gamma}_{n}(M)\right\\}.$
(4.23)

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

$\displaystyle\widehat{\beta}_{n,M}=\left[\hat{\Sigma}_{n}(M)\right]^{-1}\hat{\Gamma}_{n}(M),$
(4.24)

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

$\displaystyle\Sigma_{n}:=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[X_{i}X_{i}^{\top}\right]\in\mathbb{R}^{p\times
p},\ \ \ \text{and}\ \ \
\Gamma_{n}:=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[X_{i}Y_{i}\right]\in\mathbb{R}^{p\times
p}.$ (4.25)

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:

$\displaystyle\beta_{n,M}:=\left[\Sigma_{n}(M)\right]^{-1}\Gamma_{n}(M)\equiv\underset{\theta\in\mathbb{R}^{|M|}}{\text{arg
min}}\left\\{\theta^{\top}\Sigma_{n}(M)\theta-2\theta^{\top}\Gamma_{n}(M)\right\\}.$
(4.26)

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

$\displaystyle\underset{M\in\mathcal{M}(k)}{\text{sup}}\left\lVert\widehat{\beta}_{n,M}-\beta_{n,M}\right\rVert=o_{p}(1)\
\ \text{as}\ \ n\to\infty.$ (4.27)

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)$.

######  Remark 15.

A relevant stream of literature to uniform inference for submodel selection,
is the testing methodologies which consider the nested property between a
comparison of model as a mechanism for specification testing. In particular,
Hagemann,  (2012) considers a non-nested specification test (see, also
MacKinnon,  (1983)) while Rinaldo et al.,  (2019) develops a framework for
model selection based on sample-splitting for high dimensional assumption lean
inference purposes.

###  4.3. Non-nested Regressions and Variable Selection Specification Testing

Following the framework proposed by Hagemann,  (2012), suppose that we observe
covariates

$\displaystyle\left\\{\big{(}x_{i,1}^{\top},...,x_{i,M}^{\top}\big{)}^{\top}\in\mathbb{R}^{d_{1}+...+d_{M}}:i=1,...,n\right\\}$

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

$\displaystyle y=X_{m}\beta_{m}+u_{m},\ \
m\in\mathcal{M}:=\left\\{1,...,M\right\\},$ (4.28)

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:

######  Assumption 7.

Consider that
$\left\\{y_{i},(x_{i,m})^{\top}_{m\in\mathcal{M}\cup\left\\{m^{*}\right\\}}:i\geq
1\right\\}$ is a sequence of i.i.d random vectors. We have that
$\mathbb{E}|x_{i,m}|^{4}<\infty$ for all
$m\in\mathcal{M}\cup\left\\{m^{*}\right\\}$, where the number of elements of
$\mathcal{M}$ (number of covariates) does not depend on $n$ (sample size). For
all $m\in\mathcal{M}$, the matrices
$\mathbb{E}\left(x_{i,m}x_{i,m}^{\top}\right)$ are positive definite.

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

######  Assumption 8.

Model $m^{*}$ satisfies
$\mathbb{E}\big{(}y|\mathcal{F}\big{)}=X_{m^{*}}\beta_{m^{*}}$. Define with
$u_{i,m^{*}}:=y_{i}-x_{i,m}^{\top}\beta_{m^{*}}$ for all $i\geq 1$ and
$\mathbb{E}u^{4}_{i,m^{*}}<\infty$.

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

$\displaystyle y$
$\displaystyle=X_{m}b_{m}+\sum_{\ell\in\mathcal{M}\backslash\left\\{m\right\\}}\alpha_{\ell,m}X_{\ell}\beta_{\ell}+u,$
(4.29) $\displaystyle b_{m}$
$\displaystyle:=\left(1-\sum_{\ell\in\mathcal{M}\backslash\left\\{m\right\\}}\alpha_{\ell,m}\right)\beta_{m}.$
(4.30)

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, Hagemann,  (2012)).

######  Remark 16.

In particular, the minimum J-test proposed by Hagemann,  (2012) does not
require the correct model to be among the considered specifications and avoids
ambiguous test outcomes. More specifically, the MJ test determines with
asymptotically correct size if the correct model is among the specifications
under consideration. Furthermore, if the correct model is present, it is
chosen with probability approaching one as the sample size becomes large.

After redefining the error term appropriately, this yields that

$\displaystyle
y=X_{m}b_{m}+\sum_{\ell\in\mathcal{M}\backslash\left\\{m\right\\}}\alpha_{\ell,m}X_{\ell}\hat{\beta}_{\ell}+u,$
(4.31)

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

$\displaystyle\lambda_{n,m}$
$\displaystyle:=n^{-1/2}\left(y^{\top}P_{\ell}M_{m}y\right)_{\ell\in\mathcal{M}\backslash\left\\{m\right\\}}$
(4.32) $\displaystyle\hat{\Sigma}_{n,m}$
$\displaystyle:=n^{-1}\left(y^{\top}P_{\ell}M_{m}\hat{\Omega}M_{m}P_{\ell^{\prime}}y\right)_{\ell,\ell^{\prime}\in\mathcal{M}\in\mathcal{M}\backslash\left\\{m\right\\}}$
(4.33) $\displaystyle P_{m}$
$\displaystyle:=X_{m}\left(X_{m}^{\top}X_{m}\right)^{-1}X_{m}^{\top}\ \ \
\text{and}\ \ \ M_{m}:=I_{n}-\left(P_{m}\right)_{{\color[rgb]{1,0,0}n\times
n}}$ (4.34)

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

$\displaystyle
J_{n,m}:=\lambda_{n,m}^{\top}\hat{\Sigma}_{n,m}^{-1}\lambda_{n,m}$ (4.35)

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.

######  Lemma 11 (Hagemann,  (2012)).

Suppose that Assumptions are satisfied. Then, it holds that

  * (i).

If $m^{*}\in\mathcal{M}$, then $J_{n,m^{*}}\Rightarrow\chi^{2}_{M-1}$.

  * (ii).

For every $m\in\mathcal{M}\backslash\left\\{m\right\\}$ and every
$B\in\mathbb{R}$, we have that
$\mathsf{lim}_{n\to\infty}\mathbb{P}\big{(}J_{n,m}>B\big{)}=1$.

######  Lemma 12 (Hagemann,  (2012)).

Suppose that the conditions of the Theorem hold. If $m^{*}\in\mathcal{M}$,
then

$\displaystyle\underset{n\to\infty}{\mathsf{plim}}\
\underset{x\in\mathbb{R}}{\mathsf{sup}}\
\left|\mathbb{P}^{*}\big{(}MJ_{n}^{*}\leq
x\big{)}-\mathbb{P}^{*}\big{(}MJ_{n}^{*}\leq x\big{)}\right|=0.$ (4.36)

###### Proof.

Let $\hat{m}:=\mathsf{argmin}\mathcal{J}_{n}^{*}$. We first show that if
$m^{*}\in\mathcal{M}$, then $\hat{m}^{*}$ approximates $m^{*}$.

Without loss of generality, fix any $0<\epsilon<1$, then

$\displaystyle\mathbb{P}\big{(}\boldsymbol{1}\left\\{\exists\
m\in\mathcal{M}\backslash\left\\{\hat{m}\right\\}:J_{n,m}\leq
B\right\\}>\epsilon\big{)}\leq\mathbb{P}\big{(}\exists\
m\in\mathcal{M}\backslash\left\\{\hat{m}\right\\}:J_{n,m}\leq
B\big{)}+\mathbb{P}\big{(}\hat{m}\neq m^{*}\big{)},$ (4.37)

which converges to zero as $n\to\infty$.

In summary, it can be proved that

$\displaystyle\underset{n\to\infty}{\mathsf{plim}}\
\mathbb{P}^{*}\big{(}J^{*}_{n,\hat{m}^{*}}\leq x\big{)}\equiv K(x),$ (4.38)

where $K(x)$ is the distribution function of a $\chi^{2}_{M-1}$ variable. ∎

######  Remark 17.

Notice that although the methodology proposed by Hagemann,  (2012) considers
testing for correct model specification it differs from the literature of
conditional subvector testing. Specifically, the subvector testing method
considers the roots of the following characteristic polynomial

$\displaystyle\left|\hat{\kappa}I_{p}-n^{-1}\widehat{G}_{n}^{-1/2}\big{(}\bar{Y}_{0},W\big{)}^{\prime}Z\widehat{H}_{n}Z^{\prime}\big{(}\bar{Y}_{0},W\big{)}\widehat{G}_{n}^{-1/2}\right|=0.$
(4.39)

######  Remark 18.

A different stream of literature considers the classification of non-
overlapping and overlapping models and discuss how the relationship between
candidate models affects the asymptotic distributions of the test statistics.
The particular stream of literature was initiated by Vuong,  (1989). This
approach considers measuring the distance from the model under investigation
to the true distribution, which requires to solve the minimization problem
$\mathsf{inf}_{P\in\mathcal{P}}d(P,\mu)$.

###  4.4. Divide and Conquer Variable Selection Algorithm

Further algorithmic procedures for model selection purposes in high
dimensional linear regression models includes the framework of Meinshausen and
Bühlmann,  (2010) 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, Battey et al.,  (2018) (see, also Fan et al.,  (2015))
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

$\displaystyle\widehat{\boldsymbol{\beta}}^{d}=\widehat{\boldsymbol{\beta}}_{Lasso}^{d}(\mathcal{D}_{j})+\frac{1}{n_{k}}M^{(j)}\left(X^{(j)}\right)^{\top}\left(Y^{(j)}-X^{(j)}\widehat{\boldsymbol{\beta}}_{Lasso}^{d}(\mathcal{D}_{j})\right),$
(4.40)

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

$\displaystyle\boldsymbol{m}_{v}^{(j)}=\underset{\boldsymbol{m}}{\mathsf{argmin}}\
\boldsymbol{m}^{\top}\widehat{\boldsymbol{\Sigma}}^{(j)}\boldsymbol{m}\ \ \
\text{s.t}\ \ \
\left\lVert\widehat{\boldsymbol{\Sigma}}^{(j)}\boldsymbol{m}-\boldsymbol{e}_{v}\right\rVert_{\infty}\leq\vartheta_{1},\left\lVert\widehat{\boldsymbol{\Sigma}}^{(j)}\boldsymbol{m}\right\rVert\leq\vartheta_{2}.$
(4.41)

######  Theorem 13.

Suppose that $\mathbb{E}\big{[}\varepsilon_{1}^{4}\big{]}<\infty$ and choose
$\vartheta_{1},\vartheta_{2}$ and $k$ such that
$\theta_{1}\approx\sqrt{k\mathsf{log}d/n},\theta_{2}n^{-1/2}=o(1)$ and
$k=o\big{(}\sqrt{n}(s\mathsf{log}d)^{-1}\big{)}$. For any
$v\in\left\\{1,...,d\right\\}$,

$\displaystyle\frac{\sqrt{n}}{k}\sum_{j=1}^{k}\frac{\widehat{\beta}^{d}_{v}(\mathcal{D}_{j})-\beta_{v}^{*}}{\widehat{Q}_{v}^{(j)}}\Rightarrow\mathcal{N}(0,\sigma^{2}),\
\ \ \text{where}\ \ \
\widehat{Q}_{v}=\left(\boldsymbol{m}_{v}^{(j)\top}\widehat{\boldsymbol{\Sigma}}^{(j)}\boldsymbol{m}_{v}^{(j)}\right)^{1/2}$
(4.42)

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

$\displaystyle\bar{S}_{n}=\frac{\sqrt{n}}{k}\sum_{j=1}^{k}\frac{\widehat{\beta}^{d_{v}}(\mathcal{D}_{j})-\beta_{v}^{H}}{\bar{\sigma}\left(\boldsymbol{m}_{v}^{(j\top)}\widehat{\boldsymbol{\Sigma}}^{(j)}\boldsymbol{m}_{v}^{(j)}\right)^{1/2}}.$
(4.43)

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

$\displaystyle\widehat{\boldsymbol{\beta}}^{d}(\mathcal{D}_{j})=\widehat{\boldsymbol{\beta}}^{\lambda}(\mathcal{D}_{j})-\widehat{\Theta}^{(j)}\nabla\ell_{n_{k}}^{(j)}\left(\boldsymbol{\beta}^{\lambda}(\mathcal{D}_{j})\right),$
(4.44)

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

$\displaystyle\widehat{J}^{(j)}=\nabla^{2}\ell_{n_{k}}^{(j)}\left(\widehat{\boldsymbol{\beta}}^{\lambda}(\mathcal{D}_{j})\right)$
(4.45)

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

$\displaystyle\bar{\boldsymbol{\beta}}^{d}:=\frac{1}{k}\sum_{j=1}^{k}\widehat{\boldsymbol{\beta}}^{d}\left(\mathcal{D}_{j}\right).$
(4.46)

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:

$\displaystyle\widehat{\boldsymbol{\kappa}}_{\nu}\left(\mathcal{D}_{j}\right)=\underset{\boldsymbol{\kappa}\in\mathbb{R}^{d-1}}{\mathsf{argmin}}\left(\widehat{J}_{\nu\nu}^{(j)}-2\widehat{J}_{\nu,-\nu}^{(j)}\boldsymbol{\kappa}+\boldsymbol{\kappa}^{\top}\widehat{J}_{-\nu,-\nu}^{(j)}\boldsymbol{\kappa}+2\lambda_{\nu}\left\lVert\boldsymbol{\kappa}\right\rVert_{1}\right),$
(4.47)

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.

######  Remark 19.

Notice that Meinshausen and Bühlmann,  (2006) established a link between the
nodewise regression and the optimal linear prediction of excess asset returns
under the assumption that the returns are normally distributed. Relevant
studies include Callot et al.,  (2021) and Caner et al.,  (2023).

Introduce the following matrix

$\displaystyle\widehat{C}:=\begin{pmatrix}1&-\widehat{\kappa}_{1,2}(\mathcal{D}_{j})&\ldots&-\widehat{\kappa}_{1,d}(\mathcal{D}_{j})\\\
-\widehat{\kappa}_{2,1}(\mathcal{D}_{j})&1&\ldots&-\widehat{\kappa}_{1,d}(\mathcal{D}_{j})\\\
\vdots&\vdots&\ddots&\vdots\\\
-\widehat{\kappa}_{d,1}(\mathcal{D}_{j})&-\widehat{\kappa}_{d,2}(\mathcal{D}_{j})&\ldots&1\end{pmatrix}.$
(4.48)

Denote with

$\displaystyle\widehat{\Xi}^{(j)}=\mathsf{diag}\big{(}\mathcal{\tau}_{1}(\mathcal{D}_{j}),...,\mathcal{\tau}_{d}(\mathcal{D}_{j})\big{)},\
\ \ \text{where}\ \ \
\widehat{\mathcal{\tau}}_{\nu}\left(\mathcal{D}_{j}\right)^{2}\equiv\widehat{J}_{\nu\nu}^{(j)}-\widehat{J}_{\nu,-\nu}^{(j)}\widehat{\boldsymbol{\kappa}}_{\nu}(\mathcal{D}_{j}).$
(4.49)

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

$\displaystyle\widehat{\Theta}^{(j)}:=\left(\widehat{\Xi}^{(j)}\right)^{-2}\widehat{C}^{(j)},$
(4.50)

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

$\displaystyle\bar{S}_{n}:=\frac{\sqrt{n}}{k}\sum_{j=1}^{k}\frac{\beta_{\nu}^{d}(\mathcal{D}_{j})-\beta_{\nu}^{H}}{\sqrt{\Theta^{*}_{\nu\nu}}}$
(4.51)

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

######  Remark 20.

Therefore, the above statistical procedure provides the basis for the
statistical testing based on divide and conquer. On the other hand, the
particular methodology is only implemented for estimation and inference
purposes in high dimensional sparse models. In the sparse high dimensional
setting, nuisance parameters introduce a nonnegligible bias, causing classical
low dimensional theory to break down. On the other hand, in their high
dimensional Wald construction, the phenomenon is remedied through a debiasing
of the estimator, which gives rise to a test statistic with tractable limiting
distribution.

##### 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, Berger and Pericchi,
(1996), Berger and Mortera,  (1999)). Moreover, the framework proposed by
Vuong,  (1989) 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.

###  4.5. Sample-Splitting and Variable Importance Algorithm

In this section, we follow and discuss the framework proposed by Rinaldo et
al.,  (2019). 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

$\displaystyle x\in\mathbb{R}^{d}\mapsto\mu(x):=\mathbb{E}\big{[}Y|X=x\big{]}$

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

$\displaystyle\mathcal{D}_{n}\mapsto
w_{n}\left(\mathcal{D}_{n}\right)=\left(\widetilde{S},\widetilde{\mu}_{\widetilde{S}}\right),$
(4.52)

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, Rinaldo et al.,  (2019) 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 Rinaldo et al.,  (2019), 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}$.

####  4.5.1. 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

$\displaystyle\widehat{\beta}_{\widehat{S}}=\underset{\beta\in\mathbb{R}^{\widehat{S}}}{\mathsf{arg\
min}}\ \mathbb{E}_{X,Y}\left(Y-\beta^{\top}X_{\widehat{S}}\right)^{2},$ (4.53)

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 Rinaldo et al.,  (2019),
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

$\displaystyle\beta_{\widehat{S}}=\Sigma^{-1}_{\widehat{S}}Q_{\widehat{S}},\ \
\ \text{with}\ \ \
Q_{\widehat{S}}:=\left(Q_{\widehat{S}},j\in\widehat{S}\right)$ (4.54)

where

$\displaystyle
Q_{\widehat{S}}=\mathbb{E}_{X,Y}\left[YX_{\widehat{S}}(j)|\mathcal{D}_{n}\right]\
\ \ \text{and}\ \ \
\Sigma_{\widehat{S}}=\mathbb{E}\left[X_{\widehat{S}}X_{\widehat{S}}^{\top}|\mathcal{D}_{n}\right]$

####  4.5.2. 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, Rinaldo et
al.,  (2019)).

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

$\displaystyle\gamma_{\widehat{S}}(j)=\mathbb{E}\left[\big{|}Y-\widehat{\beta}^{\top}_{\widehat{S}(j)}X_{\widehat{S}(j)}\big{|}-\big{|}Y-\widehat{\beta}^{\top}_{\widehat{S}}X_{\widehat{S}}\big{|}\bigg{|}\mathcal{D}_{n}\right].$
(4.55)

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

$\displaystyle\phi_{\widehat{S}}(j)=\mathsf{median}\left[\big{|}Y-\widehat{\beta}^{\top}_{\widehat{S}(j)}X_{\widehat{S}(j)}\big{|}-\big{|}Y-\widehat{\beta}^{\top}_{\widehat{S}}X_{\widehat{S}}\big{|}\bigg{|}\mathcal{D}_{n}\right].$
(4.56)

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

####  4.5.3. 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

$\displaystyle\rho_{\widehat{S}}=\mathbb{E}\left[\big{|}Y-\widehat{\beta}^{\top}_{\widehat{S}}X_{\widehat{S}}\big{|}\bigg{|}\mathcal{D}_{n}\right],$
(4.57)

where $\widehat{\beta}_{\widehat{S}}$ is computed based on $\mathcal{D}_{n}$.
In particular, the main idea of the approach proposed by Rinaldo et al.,
(2019) 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:

  * •

Concentration:

$\displaystyle\underset{n\to\infty}{\mathsf{lim\
sup}}\underset{w_{n}\in\mathcal{W}_{n}}{\mathsf{sup}}\underset{P\in\mathcal{Q}_{n}}{\text{sup}}\mathbb{P}\left(\big{\|}\widehat{\theta}_{\widehat{S}}-\theta_{\widehat{S}}\big{\|}>r_{n}\right)\to
0,$ (4.58)

  * •

Coverage Validity:

$\displaystyle\underset{n\to\infty}{\mathsf{lim\
inf}}\underset{w_{n}\in\mathcal{W}_{n}}{\mathsf{inf}}\underset{P\in\mathcal{Q}_{n}}{\mathsf{inf}}\mathbb{P}\left(\theta_{\widehat{S}}\in\widehat{C}_{\widehat{S}}\right)\geq
1-\alpha,$ (4.59)

  * •

Accuracy:

$\displaystyle\underset{n\to\infty}{\mathsf{lim\
sup}}\underset{w_{n}\in\mathcal{W}_{n}}{\mathsf{sup}}\underset{P\in\mathcal{Q}_{n}}{\mathsf{sup}}\mathbb{P}\left(\nu\left(\widehat{C}_{\widehat{S}}\right)>\epsilon_{n}\right)\to
0,$ (4.60)

where $\alpha\in(0,1)$ is a prespecified level of significance,
$\mathcal{W}_{n}$ is the set of all the model selection and estimation
procedures on samples of size $n$, $r_{n}$ and $\epsilon_{n}$ both vanish as
$n\to\infty$ and $\nu$ is the volume (Lebsegue measure) of the set. The
probability statements above take into account both the randomness in the
sample $\mathcal{D}_{n}$ and the randomness associated to splitting it into
halves.

######  Remark 21.

The property that the coverage of $\widetilde{C}_{\widetilde{S}}$ is
guaranteed uniformly over the entire class $\mathcal{Q}_{n}$ is known as
(asymptotic) honesty. Moreover, the confidence sets are for the random
parameters (based on half the data) but the uniform coverage , accuracy and
concentration guarantee to hold with respect to the distribution of the entire
sample and the randomness associated to the splitting of the full sample.

######  Remark 22.

Since we are particularly interested to study suitable variable selection
methodologies for time series regression models, notice that usually the
subsampling approach in such modelling settings requires to consider the
dependence structure of the data. For example, the greater the number of
subsamples, the larger the bias due to the fact that consistent estimation of
long-memory parameters requires larger samples. Determining the length of
subsamples under the presence of long memory in time series data is another
aspect of concern.

####  4.5.4. Confidence sets for the projection parameters: The bootstrap

Following Rinaldo et al.,  (2019), 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

$\displaystyle\mathbb{P}\left(\sqrt{n}\left\|\widehat{\beta}_{\widehat{S}}^{*}-\widehat{\beta}_{\widehat{S}}\right\|\leq\widehat{t}^{*}_{\alpha}\big{|}\mathcal{D}_{2,n}\right)\geq
1-\alpha.$ (4.61)

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

$\displaystyle\mathbb{P}\left(\sqrt{n}\left|\widehat{\beta}_{\widehat{S}}^{*}(j)-\widehat{\beta}_{\widehat{S}}(j)\right|\leq\widetilde{t}^{*}_{j},\
\forall\ j\big{|}\mathcal{D}_{2,n}\right)\geq 1-\alpha.$ (4.62)

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

$\displaystyle\mathbb{P}\left(\sqrt{n}\left|\widehat{\beta}_{\widehat{S}}^{*}(j)-\widehat{\beta}_{\widehat{S}}(j)\right|>\widetilde{t}^{*}_{j},\
\forall\ j\big{|}\mathcal{D}_{2,n}\right)\leq\frac{\alpha}{k}.$ (4.63)

Consider the following two bootstrap confidence sets:

$\displaystyle\widehat{C}^{*}_{\widehat{S}}$
$\displaystyle=\left\\{\beta\in\mathbb{R}^{\widehat{S}}:\left|\beta-\widehat{\beta}_{\widehat{S}}\right|\leq\frac{\widehat{t}_{\alpha}^{*}}{\sqrt{n}}\right\\},$
(4.64) $\displaystyle\widetilde{C}^{*}_{\widehat{S}}$
$\displaystyle=\left\\{\beta\in\mathbb{R}^{\widehat{S}}:\left|\beta(j)-\widehat{\beta}_{\widehat{S}}(j)\right|\leq\frac{\widetilde{t}_{j}^{*}}{\sqrt{n}},\
\forall\ j\in\widehat{S}\right\\}.$ (4.65)

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

$\displaystyle\sqrt{n}\left(\widehat{\beta}-\beta\right)\to\mathcal{N}\left(0,\Gamma\right),\
\
\Gamma=\Sigma^{-1}\mathbb{E}\big{[}\left(Y-\beta^{\top}X\right)^{2}\big{]}\Sigma^{-1}.$
(4.66)

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 Rinaldo et al.,  (2019) 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.

##### 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.

###  4.6. 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 Chudik et al.,  (2018). Therefore, in this section we discuss
their proposed statistical methodology, so-called One Covariate at a Time
Multiple Testing (OCMT) (see, also Zhang,  (1993)). Specifically, the
particular procedure is computationally simple and fast even for extremely
large data sets (see, also Romano and Wolf,  (2005)).

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

$\displaystyle\hat{\boldsymbol{\beta}}=\underset{\boldsymbol{\beta}\in\mathbb{R}^{p}}{\mathsf{arg\
min}}\sum_{t=1}^{T}\left\\{\big{(}y_{t}-x_{nt}^{\prime}\boldsymbol{\beta}\big{)}^{2}+P_{\lambda}(\boldsymbol{\beta})\right\\}.$
(4.67)

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 Fithian et al.,  (2015). These approaches construct
regression models by selecting variables from active sets, based on a sequence
of sets.

####  4.6.1. Rejection Principle of Familywise Error Control

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

######  Theorem 14 (Finner and Strassburger,  (2002)).

Let $\alpha,\gamma\in(0,1)$. If the null hypothesis $H_{0,j}:\beta_{j}=0$ gets
rejected whenever $Q_{j}(\gamma)\leq\alpha$, then the FWER is asymptotically
controlled at level $\alpha$, that is,

$\displaystyle\underset{n\to\infty}{\mathsf{lim\ sup}}\
\mathbb{P}\bigg{(}\underset{j\in\mathbb{N}}{\mathsf{min}}\
Q_{j}(\gamma)\leq\alpha\bigg{)}\leq\alpha.$ (4.68)

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 Romano and Wolf,  (2005) (see, also Politis et al.,  (1999)).
Moreover, the paper of Goeman and Solari,  (2010) 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.

######  Theorem 15 (Sequential rejection principle, Goeman and Solari,
(2010)).

Suppose that for every $\mathcal{R}\subseteq\mathcal{S}\subset\mathcal{H}$,
almost surely, $\mathcal{N}(\mathcal{R})\subseteq(\mathcal{S})\cup\mathcal{S}$
and that for every $M\in\mathbb{M}$,

$\displaystyle\mathbb{P}\big{(}\mathcal{N}(\mathcal{F}(M))\subseteq\mathcal{F}(M)\big{)}\geq
1-\alpha.$ (4.69)

Then, for every $M\in\mathbb{M}$,
$\mathbb{P}\big{(}\mathcal{R}_{\infty}\subseteq\mathcal{F}(M)\big{)}\geq
1-\alpha$.

######  Remark 23.

The first condition - monotonicity condition, guarantees that no false
rejection in the critical case (during the single-step), which implies no
false rejection in situations with fewer rejections than in the critical case
so that type I error control in the critical case is sufficient for overall
FWE control of the sequential procedure. The second condition above - single-
step condition, guarantees FWE control when we have rejected all false null
hypotheses and none of the true ones. A general admissibility criterion is the
case of restricted combinations999A standard example concerns testing pairwise
equality of means in a one-way ANOVA model: if any single null hypothesis is
false, it is not possible that all other null hypotheses are simultaneously
true (see, also Vesely et al.,  (2021)). can be also constructed. Restricted
combinations occur if, for some $\mathcal{R}\subseteq\mathcal{H}$, there is no
model $M\in\mathbb{M}$ such that $\mathcal{R}=\mathcal{F}(M)$ (see, Goeman and
Solari,  (2010)).

  * •

Resampling-based multiple testing procedures use resampling techniques to let
the multiple testing procedure estimate or accommodate the actual dependence
structure between the test statistics. For example, resampling of a test
statistic under the complete null hypothesis, using permutations or the
bootstrap, can give consistent estimates of the desired quantiles. However,
regardless of the underlying assumptions, consistent estimation of the
quantiles of $\mathsf{max}_{H\in\mathcal{T}(M)}S_{H}$ only guarantees control
of the familywise error in an asymptotic sense. Resampling-based methods with
exact FWE control, imply that one can obtain control of the FWE by
generalizing the treatment of permutation testing to a multiple testing
procedure. In order to define a resampling-based sequentially rejective
multiple testing procedure with exact familywise error control, we choose a
set $\boldsymbol{\pi}=\left\\{\pi_{1},...,\pi_{r}\right\\}$ of $r$ functions
that we shall refer to as null-invariant transformations.

  * •

Graph-based procedures. Another application of the sequential rejection
principle, which is of interest in the statistics, econometrics and finance
literature, is the development of multiple testing101010Notice that for
multiple testing procedures researchers are also interested in reporting
multiplicity-adjusted p-values. Such multiplicity-adjusted p-values are
defined for each null hypothesis as the smallest $\alpha-$level that allows
rejection of that hypothesis (see also the method proposed by McKeague and
Qian,  (2015) and Huang et al.,  (2019)). procedures for graph-structured
hypotheses. Specific procedures for controlling the familywise error for
graph-structured hypotheses have been proposed by several authors.

######  Example 20 (Neighborhood Selection with the Lasso, Meinshausen and
Bühlmann,  (2006)).

Consider the $p-$dimensional multivariate normal distributed random variable
$X=\left(X_{1},...,X_{p}\right)\sim\mathcal{N}\left(\mu,\Sigma\right)$. This
includes the Gaussian linear models where $X_{1}$ is the response variable and
$\left\\{X_{k},2\leq k\leq p\right\\}$ are the predictor variables. The
conditional independence structure of the distribution can be represented by a
graphical model $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\left\\{1,...,p\right\\}$ is the set of nodes and $\mathcal{E}$
is the set of edges in $\mathcal{V}\times\mathcal{V}$. A pair $(a,b)$ is
contained in the edge set $\mathcal{E}$ if and only if $X_{a}$ is
conditionally dependent on $X_{b}$, given all remaining variables
$X_{\mathcal{V}\backslash\left\\{a,b\right\\}}$. Every pair of variables not
contained in the edge set is conditionally independent, given all remaining
variables, and corresponds to a zero entry in the inverse covariance matrix.
When predicting a variable $X_{a}$ with all remaining variables
$\left\\{X_{k};k\in\Gamma(n)\backslash\left\\{a\right\\}\right\\}$, the
vanishing Lasso coefficient estimates identify asymptotically the neighborhood
of a node $a$ in the graph. Let the $n\times p(n)-$dimensional matrix
$\mathbf{X}$ contain $n$ independent observations of $X$, so that the columns
$\mathbf{X}_{a}$ correspond for all $a\in\Gamma(n)$ to the vector of $n$
independent observations of $X_{a}$. Let $\langle.,.\rangle$ be the usual
inner product on $\mathbb{R}^{n}$ and $\left\lVert\ .\ \right\rVert_{2}$ the
corresponding norm. The Lasso estimate $\hat{\theta}^{a,\lambda}$ of
$\theta^{a}$ is given by

$\displaystyle\hat{\theta}^{a,\lambda}=\underset{\theta:\theta_{a}=0}{\mathsf{arg\
min}}\left(n^{-1}\left\lVert\mathbf{X}_{a}-\mathbf{X}\theta\right\rVert_{2}^{2}+\lambda\left\lVert\theta\right\rVert_{1}\right),\
\ \left\lVert\theta\right\rVert_{1}=\sum_{b\in\Gamma(n)}|\theta_{b}|,$ (4.70)

where $\ell_{1}-$norm of the coefficient vector. Furthermore, Salgueiro et
al.,  (2005) and Salgueiro et al.,  (2006) propose a formal statistical
framework for edge exclusion in graphical Gaussian models. Moreover,
statistical theory for constructing relevant matrix moments functions is
presented by Roverato and Whittaker,  (1998). Recently, Fan et al.,  (2020)
propose a conditional dependence measure in high-dimensional undirected
graphical models (UGM). In particular, the UGM approach examines the internal
conditional dependency structure of a multivariate random vector.

####  4.6.2. Subgroup Selection Methodology

Lastly, we briefly discuss the subgroup selection methodology proposed by
Reeve et al.,  (2021)111111Dr. 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

$\displaystyle\mathcal{X}_{\tau}(\eta):=\big{\\{}x\in\mathbb{R}^{d}:\eta(x)\geq\tau\big{\\}}.$
(4.71)

In other words, the statistical mechanism of Reeve et al.,  (2021) 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{\\}}$.

######  Proposition 3 (Type 1 error guarantee,Reeve et al.,  (2021)).

Let $\mathcal{P}$ be a family of distributions $P$ on pairs $(X,Y)$ and choose
a significance level $\alpha\in(0,1)$. We say that the data-dependent subgroup
$\hat{A}$ controls Type 1 error at the level $\alpha$ over the class
$\mathcal{P}$ if

$\displaystyle\underset{P\in\mathcal{P}}{\mathsf{inf}}\
\mathbb{P}_{P}\left(\hat{A}(\mathcal{D})\subseteq\mathcal{X}_{\tau}(\eta)\right)\geq
1-\alpha.$ (4.72)

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

$\displaystyle\underset{P\in\mathcal{P}}{\mathsf{inf}}\
\mathbb{P}_{P}\left(\hat{A}(\mathcal{D})\subseteq\mathcal{X}_{\tau}(\eta)\right)\geq
1-\alpha.$

######  Definition 7 (Holder Class, Reeve et al.,  (2021)).

Given that $(\beta,C_{S})\in(0,1]\times[1,\infty)$, we let
$\mathcal{P}_{Hol,\tau}(\beta,C_{S})$ denote the class of all distributions
$P$ on $\mathbb{R}^{d}\times[0,1]$ with marginal $\mu$ on $\mathbb{R}^{d}$
such that the regression function $\eta$ is $(\beta,C_{S})-$Holder on
$\mathcal{X}_{\tau}(\eta)\cap\mathsf{supp}(\mu)$ in the sence that

$\displaystyle\left|\eta(x^{\prime})-\eta(x)\right|\leq C_{X}.\left\lVert
x^{\prime}-x\right\rVert_{\infty}^{\beta},$ (4.73)

for all $x,x^{\prime}\mathcal{X}_{\tau}(\eta)\cap\mathsf{supp}(\mu)$.

Further details on the implementation of this methodology as well as related
statistical theory and numerical illustrations are presented in Reeve et al.,
(2021). More recently the isotonic subgroup selection approach is proposed by
Müller et al.,  (2023).

######  Remark 24.

Notice that according to Reeve et al.,  (2021), the subgroup selection
approach requires from the statistician to use both FWE control as well as the
selection step. In other words, the subgroup selection approach corresponds to
a data-dependent selection methodology, such that
$\hat{\mathcal{A}}=\hat{\mathcal{A}}(D)$. Therefore, a special class of
functions needs to be considered which allows to test for local null
hypothesis.

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:

$\displaystyle\underset{p\in P}{\mathsf{inf}}\
\mathbb{P}_{P}\left(\hat{\mathcal{A}}(D)\subset\mathcal{X}_{\uptau}(n)\right)\geq
1-\alpha.$ (4.74)

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.

###  4.7. Model Selection in Cointegrating Regressions

According to Mendes,  (2011), 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, Abadir et al.,  (1999)).

######  Example 21.

Consider the following cointegration regression model given by the expression

$\displaystyle
y_{t}=\alpha_{0}+\beta_{0}^{\prime}x_{t}+\gamma_{0}^{\prime}z_{t}+u_{t}$
(4.75)

where $\beta_{0}$ is an $n_{1}\times 1$ vector of parameters and $\gamma_{0}$
is an $n_{2}\times 1$ vector of parameters. Furthermore, the process
$\left\\{x_{t}\right\\}_{t=1}^{\infty}$ satisfies the following integrated
stochastic process

$\displaystyle x_{t}=x_{t-1}+v_{t},$ (4.76)

Therefore, the main aim shall be to control the number of $I(0)$ variables in
the model and we assume that $n\equiv n_{1}+n_{2}$ to possibly be grater than
$T$, but only a fraction of those coefficients are in fact nonzero.
Furthermore, we assume without loss of generality that each coefficient
vectors can be partitioned into zero and non-zero coefficients, such that
$\beta_{0}=\big{(}\beta_{0}(1)^{\prime},\beta_{0}(2)^{\prime}\big{)}^{\prime}$
and
$\gamma_{0}=\big{(}\gamma_{0}(1)^{\prime},\gamma_{0}(2)^{\prime}\big{)}^{\prime}$,
with all non-zero coefficients stacked first, where $\beta_{0}(1)$ is
$q_{1}\times 1$ and $\gamma_{0}(1)$ is $q_{2}\times 1$. In particular, we
assume that the number of non-zero coefficients, measured by $q_{1}$, is fixed
(does not depend on $T$), while the number of zero coefficients, measured by
$q_{2}$, may depend on $T$, also set $q=(q_{1}+q_{2})$.

Then, the Adaptive Lasso estimate is given by (see, also Medeiros and Mendes,
(2017))

$\displaystyle\left(\hat{\beta},\hat{\gamma}\right)=\underset{\beta,\gamma}{\mathsf{argmin}}\left\lVert
Y-X\beta-Z\gamma\right\rVert_{2}^{2}+\lambda_{1}\sum_{j=1}^{n}\lambda_{1j}|\beta_{j}|+\lambda_{2}\sum_{j=1}^{n}\lambda_{2j}|\gamma_{j}|,$
(4.77)

######  Lemma 13 (KKT Conditions).

The solutions
$\hat{\beta}=\big{(}\hat{\beta}(1)^{\prime},\hat{\beta}(2)^{\prime}\big{)}$
and
$\hat{\gamma}=\big{(}\hat{\gamma}(1)^{\prime},\hat{\gamma}(2)^{\prime}\big{)}$
to the minimization problem above exists if:

$\displaystyle\frac{\partial\left\lVert
Y-X\beta-Z\gamma\right\rVert_{2}^{2}}{\partial\beta_{j}}\bigg{|}_{\beta_{j}(1)=\hat{\beta}_{j}(1)}$
$\displaystyle=\mathsf{sgn}\left(\hat{\beta}_{j}(1)\right)\lambda_{1}\lambda_{1j}$
(4.78) $\displaystyle\frac{\partial\left\lVert
Y-X\beta-Z\gamma\right\rVert_{2}^{2}}{\partial\beta_{j}}\bigg{|}_{\gamma_{j}(1)=\hat{\gamma}_{j}(1)}$
$\displaystyle=\mathsf{sgn}\left(\hat{\gamma}_{j}(1)\right)\lambda_{2}\lambda_{2j}$
(4.79)

####  4.7.1. Model Selection Consistency and Oracle Property

We focus on showing that under certain conditions on $n,p$ and
$\lambda^{\prime}$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.

######  Lemma 14.

Let

$\displaystyle\Omega_{\infty}$
$\displaystyle=\begin{pmatrix}\Omega_{X,\infty}&\boldsymbol{0}\\\
\boldsymbol{0}^{\prime}&\Omega_{X,\infty}\end{pmatrix}$ (4.80)
$\displaystyle\Omega_{X,\infty}$
$\displaystyle=\int_{0}^{1}B_{X}(1)B_{X}(1)^{\prime}(r)dr\ \ \ \text{and}\ \ \
\Omega_{Z,\infty}=\Sigma_{Z(1)^{2}},$ (4.81)

where for any $0\leq r\leq 1$,
$B_{X(1)}(r)=\underset{T\to\infty}{\mathsf{lim}}\sum_{t=1}^{\lfloor
rT\rfloor}\boldsymbol{v}_{t}(1)$. Similarly, split the matrix $\Omega_{11}$
into the partition

$\displaystyle\Omega_{11}=\begin{pmatrix}\Omega_{X(1)^{2}}&\Omega_{Z(1)X(1)}\\\
\Omega^{\prime}_{Z(1)X(1)}&\Omega_{Z(1)^{2}}\end{pmatrix}=\begin{pmatrix}\frac{1}{T^{2}}X(1)^{\prime}X(1)&\frac{1}{T^{3/2}}Z(1)^{\prime}X(1)\\\
\frac{1}{T^{3/2}}X(1)^{\prime}Z(1)&\frac{1}{T}Z(1)^{\prime}Z(1)\end{pmatrix}.$
(4.82)

Hence,

$\displaystyle\mathbb{P}\big{(}\mathcal{A}^{c}_{T}(X)\big{)}$
$\displaystyle=\mathbb{P}\left(\left\\{\left[\big{|}T^{-1}\Omega_{X,\infty}^{-1}X(1)^{\prime}U\big{|}\right]_{j}>T|\beta_{0j}|\right\\},j=1,...,q_{1}\right)+o_{p}(1),$
$\displaystyle\leq\sum_{j=1}^{q_{1}}\mathbb{P}\left(\left[\big{|}T^{-1}\Omega_{X,\infty}^{-1}X(1)^{\prime}U\big{|}\right]_{j}>T|\beta_{0j}|\right)+o_{p}(1),$
$\displaystyle\leq\frac{q_{1}}{T^{2}\beta^{2}_{*}}\ \underset{1\leq j\leq
q_{1}}{\mathsf{max}}\
\mathbb{E}\left(\left[T^{-1}\big{|}\Omega_{X,\infty}^{-1}X(1)^{\prime}U\big{|}\right]^{2}_{j}\right)\to
0.$

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, Gupta and Seo,  (2019), Baillie et al.,  (2022)). 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 Fan et al.,  (2023). 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,
Pitarakis,  (2023) and Gonzalo and Pitarakis,  (2023)). Moreover, choosing
variables from irrelevant cointegrating relations is another important issue
(see, Khalaf and Richard,  (2020) and Richard,  (2023)).

####  4.7.2. 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, Kleibergen and Paap,  (2006)). Furthermore, Cragg and
Donald,  (1997) propose suitable testing procedures for determining the rank
of a matrix. In particular, the authors examine the use of model selection