# Lag selection and estimation of stable parameters for
multiple autoregressive processes through
convex programming

Somnath Chakrabortylabel=e1]somnath.chakraborty@ruhr-uni-
bochum.de0000-0002-5369-3053 [ Johannes Ledererlabel=e2]johannes.lederer@ruhr-
uni-bochum.de0000-0002-5369-3053 [ Rainer von
<EMAIL_ADDRESS>[ Fakultät für Mathematik, Ruhr-
Universität Bochum, 44801 Bochum, Deutschlandpresep=, ]e1,e2  Institut de
Statistique, Biostatistique et Sciences Actuarielles, LIDAM, UCLouvain, 1348
Louvain-la-Neuve, Belgiumpresep=, ]e3

###### Abstract

Motivated by a variety of applications, high-dimensional time series have
become an active topic of research. In particular, several methods and finite-
sample theories for individual stable autoregressive processes with known lag
have become available very recently. We, instead, consider multiple stable
autoregressive processes that share an unknown lag. We use information across
the different processes to simultaneously select the lag and estimate the
parameters. We prove that the estimated process is stable, and we establish
rates for the forecasting error that can outmatch the known rate in our
setting. Our insights on the lag selection and the stability are also of
interest for the case of individual autoregressive processes.

62M10,

62L99,

49M29,

regularised least square,

LASSO,

autoregressive process,

hierarchical-group norm,

dual norm,

stability,

sample complexity,

###### keywords:

[class=MSC]

###### keywords:

, and

##  1 Introduction

Today’s world of acquisition of complex data in areas such diverse as
macroeconomics and finance, everyday weather predictions, brain imaging, and
many more, has called for intelligent model approaches that avoid needing to
use potentially a (too) high number of model parameters per available sample
size. Moreover, often these data are of high dimensionality - as they arise
together in a panel or in the form of a multivariate vector. These stylized
facts render the purpose of predicting the evolution of these data into the
(near) future really challenging. To face this challenge choosing a data
generating model that assumes some common underlying structure relating the
different components of the observed multivariate data set will not only turn
out to be advantageous but reflects the observation that the different series
do not behave independently from each other - they might actually be driven by
latent (i.e. unobservable) mechanism (such as a leading economic indicator, or
a global climate trend, etc, often modelled by a latent factor model).
Moreover, we almost always observe serial correlation between present and past
observations, which traditionally has been modelled by assuming some sort of
weak dependence over time (translating into dynamic latent factor models,
e.g., Forni et al. (2000).

In this context, as factor modelling does not necessarily allow for component-
wise prediction, the approach of (parametric) vector autoregression (VAR) has
already for a long time become an overly prominent tool for modeling such
multivariate time series - with in particular the idea that the common serial
dependence is limited by the existence of a common maximal lag-order for all
components. However, as the number of component series is increased, VAR
models have the known tendency to become overparametrized. In the virtue of
having to do with a high-dimensional parameter estimation problem, more recent
possibilities to address this issue are regularized approaches, such as the
LASSO for estimating the parameters of these models (essentially by some kind
of regularised least-squares approach, see, for example, Nardi and Rinaldo
(2011)). This is in contrast to more traditional approaches (based mostly on
information criteria for lag-order selection such as AIC, BIC, etc.) which
address overparametrization by selecting a low lag order, based on the
assumption of short range dependence, assuming that a universal lag order
applies to all components. For a good forecast performance in a high-
dimensional context, these approaches turned out to fall behind the LASSO -
which, until recently, did however not incorporate the notion of lag order
selection. It has been only the recent work by Nicholson et al. (2020) that
proposed a class of hierarchical lag structures that embed the notion of lag
selection into a convex regularizer. The key modeling tool has been a group
LASSO with nested groups which guarantees that the sparsity pattern of lag
coefficients honors the VAR’s ordered structure. For more details on the
literature on dimension reduction methods which address the VAR’s
overparametrization problem we refer to Section 2 of the mentioned work by
Nicholson et al. (2020). A clear shortcoming, however, of this approach is the
necessity to model all components of the observed multivariate time series to
be of the same data length, a constraint in classical VAR-modelling that
cannot be circumvented.

Motivated by the approach of Nicholson et al. (2020), in this paper, we
propose a method to analyse multiple stable autoregressive processes of
(potentially) different lengths in the framework of regularized LASSO, where
the regularization is achieved via an overlapping group-norm that induces
sparsity at the group level. Moreover, we show that, even in absence of any
information on the maximum lag of the processes, the proposed framework
estimates the true lag and the coefficients of the AR model. Finally, we show
that the model fitted with the AR coefficients returned by this proposed
method is stable. As our results on statistical guarantees are essentially of
non-asymptotic nature - interesting even in the context of observing a single
time series - we first review the (sparse) literature on those non-asymptotic
results in a time series context, before we turn in more detail to the
similarities and differences between our and the approach of Nicholson et al.
(2020).

Most of the research in time series analysis — until recently — focused on
deriving asymptotic behaviour of the predictors. This severely restricts
applicability of these results, especially in the reign of low sample-to-
predictable ratio. Popular approaches to overcome this nuisance of dependency
have been using the assumption of stability (leading to stationarity) of the
data-generating process. For example, both Negahban and Wainwright (2011) and
Loh and Wainwright (2012) used stability in deriving the guarantees in small
sample regimes; however, these works established these results under the
condition that the coefficient matrices are severely norm-bounded (namely, the
sum of the operator-norms of the coefficient matrices is smaller than 1),
which is much stronger than stability of the process determined by those
coefficients. Recently, Basu and Michailidis (2015) made a big stride towards
understanding the effect of temporal and cross-sectional dependency in the
small sample regime. The underlying hypothesis in that work was that data be
amenable to modelling via stable vector autoregressive process; they tracked
the restrictions on the spectral domain — as enforced by stability of the
process — and derived non-asymptotic prediction guarantees for high
dimensional vector autoregressive process with Gaussian white noise
innovation. Several follow-up works (see e.g. Wong et al. (2020), Masini et-al
(2022), and the references there-in) then extended their results to the case
where the innovations are heavy-tailed, and moreover, they derived guarantees
assuming only the stationarity and finiteness of second moment of the
underlying process — conditions weaker than stability.

In the remainder of this Introduction we go now into more details about the
relation of our approach to the one of Nicholson et al. (2020). Essentially,
the latter contributed by finding out that algorithmically the method by
Mairal et al. (2011) (and its computationally faster amendment by Tseng
(2009)) can be used in such a setting to address, in the presence of prior
information on an upper bound $L$ of the true unknown lag order, estimation of
order $L$ vector-autoregressive models of dimension $M$, abbreviated
${\mbox{VAR}}_{M}(L)$ in the sequel. Our work can now be seen as fitting a
common (i.e. “diagonal vector”) autoregressive model to a multivariate time
series of dimension $M$, i.e. a panel of $M$ observed (univariate) time series
of in general not equal lengths $n_{m},1\leq m\leq M$. We address the
challenging question on how to choose, solely from the information available
in the model for the observed data, a common appropriate lag order $L$ that
allows us to phrase and solve our problem via a penalised LASSO approach. We
derive non-asymptotic bounds on the multivariate one-step ahead prediction
error and estimate the collection of the autoregressive coefficients
$\mathbf{\beta}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\beta^{m}_{1},\ldots,\beta^{m}_{L})_{1\leq
m\leq M}$ under the paradigm of sparseness. Assuming that there is a common
true unknown lag order $L_{0}$ that generated our $M$ time series, our
algorithm, akin Nicholson et al. (2020) (and Mairal et al. (2011)) for fitting
a common model is based on a modification of a hierarchical group LASSO
approach: We first determine an appropriate (minimal) upper bound $L\geq
L_{0}$ depending essentially on the sample size $n_{\min}$ of the shortest
observed time series component (and on $M$, of course) from a thorough
analysis of the theoretical complexity of our group-LASSO based approach. With
this appropriate $L$, necessary to embed our autoregression problem into the
framework of high-dimensional multivariate regression, we transfer existing
technology on LASSO estimation (with overlapping hierarchically constructed
groups) to our problem. We derive statistical guarantees for the estimators
$\mathbf{\hat{\beta}}$, solution of our aforementioned learning algorithm, and
for the estimator $\hat{L}_{0}$ (essentially taken from the support of
$\mathbf{\hat{\beta}}$). More specifically we deliver non-asymptotic bounds on
the multivariate one-step ahead prediction error, on the estimation error of
$\mathbf{\beta}$, on the false discoveries for the support of
$\mathbf{\beta}$, and quite innovatively on the stability of the fitted model.
For the latter, we show that the fitted autoregressive model of order
$\hat{L}_{0}$ with estimated coefficients $\mathbf{\hat{\beta}}$ fulfils the
conditions of the true model for stability (via a more explicit concept of
$\varepsilon$-stability that we introduce to asses the difficulty of the
statistical estimation problem).

In the following paragraph we are even more explicit about the exact nature of
our contributions motivated from the existing limitations of the current
approaches we found in the literature.

##### Current limitations and our contributions

Some of the questions that are not sufficiently addressed in recent existing
work on non-asymptotic time series analysis are as follows.

  * L1

Nicholson et al. (2020) illustrate that an approach based on LASSO with
overlapping groups (such as in Mairal et al. (2011)) can determine the
component-wise lag orders $L_{m}$ of stable, high-dimensional vector
autoregressive processes (including the ”cross-over” lag order $L_{ij}$ of the
dependence of the $i-$th component on the $j-$th component of a VAR model).
Their formulation based on (dual) convex programming is computationally
attractive and intriguing more generally, but it requires a uniform upper-
bound on the $L_{i,j}$ as a parameter both in their theoretical bounds and in
practice. The current literature either ignores this issue altogether or sets
those bounds based on model selection such as AIC or BIC or via Bayesian
shrinkage, whose suitability is unclear here.

  * C1

We mimic Nicholson et al. (2020)’s methodological approach, but we establish a
suitable upper bound for the lag orders. Our choice on this minimal upper
bound guarantees the following:

    1. 1.

the smallness of the empirical prediction risk (see Theorem 4.3 and the
discussion following the theorem),

    2. 2.

a resulting tuning parameter that is not too large (see Corollary 4.5), and

    3. 3.

the restricted eigenvalue property of the data matrix to hold (Theorem 4.6 and
the Corollary 4.7 immediately after that).

In this sense, our upper bound on the lag order is optimal for our setup.

  * L2

Real world applications often involve multiple univariate, decoupled time
series. This would translate to VAR modeling with a diagonal coefficient
matrix. However, purely transferring existing results from a
${\mbox{VAR}}_{M}(L)$ modeling approach can be cumbersome in practice for
situations in which the number of available samples for each individual
component time series is not the same: obviously needing to chop off the
samples in order to work with the minimal individual sample size could result
in possibly weaker theoretical guarantees and practical performances (see
below).

  * C2

Applying a hierarchical group norm enables us to derive prediction guarantees
that depend on the average number of samples per component (instead of the
minimum number of samples per model, as would be the case had we translated
naively to a $M$-dimensional VAR model). More specifically, availability of
perfect information on the true lag $L_{0}$ and
$n_{1}=L_{0}+T_{1},\cdots,n_{m}=L_{0}+T_{M}$ (respective) number of samples
for the $M$ individual components yields the following: the
${\mbox{VAR}}_{M}(L)$ translation would result in the following provable error
(see Nicholson et al. (2020)), stating that with high probability

$|\\!|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}|\\!|_{2}=O\left(\sqrt{\frac{\log(M^{2}L_{0})}{n_{\min}M}}\right)\,,$
(1)

whereas the error from the algorithm we describe is of the order (see equation
(33) below)

$|\\!|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}|\\!|_{2}=O\left(\sqrt{\frac{\log(ML)}{D}}\right)\,.$
(2)

Here $D=T_{1}+\cdots+T_{M}$ is the total number of ”postsamples”
($T_{m}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=n_{m}-L_{0}$).
Furthermore, our strategy results in weaker dependency of the error on the
behaviour of the reverse characteristic polynomial on the unit disk, unlike in
the relevant ${\mbox{VAR}}_{M}(L)$ model translation (compare with Basu and
Michailidis (2015)).

  * L3

Recall that a (univariate) autoregressive process
$X_{t}=a_{1}X_{t-1}+\cdots+a_{L}X_{t-L}+U_{t}$ is stable if the “reverse
characteristics polynomial” $1-a_{1}z-\cdots-a_{L}z^{L}$ has no complex roots
on the closed unit disk; it is known that stability implies stationarity (see
Section 2 below). But what about stability for fitted autoregressive models?
While this question has an affirmative answer in the special case of Yule-
Walker estimation of the coefficients of a univariate autoregressive process
(known however to be less efficient), the question of stability or
stationarity of the process reconstructed from the parameters estimated by
LASSO-based approaches does not seem to have been addressed in the literature.
However, starting with a stable process to have generated the input
observations of these devised algorithms, it is reasonable to expect that the
reconstructed process be stable (and thus, multi-step predictions be reliable
as well).

  * C3

We show, in Theorem 4.10, that the process reconstructed from the parameters
returned by our algorithm is stable when the samples available as input are
generated by stable processes. As Basu and Michailidis (2015) showed, a
measure of stability for an autoregressive process is, equivalently, a
boundedness criteria on the spectral density, and the boundedness in the
Fourier domain translates in a sense to ‘smoothness’ of the process in the
temporal domain. Thus, as, intuitively, the stable autoregressive processes
form a ‘smooth’ subclass, it is desirable that any algorithm for learning the
parameters of processes from this smooth subclass should return estimators
lying in this subclass; in this paper, this is ensured by Theorem 4.10.

It is important to mention here that the results in this paper demonstrate
that our overlapping group-lasso approach yields a stable process when the
underlying process is stable as well. The aim of the paper is not to choose
the most optimal tuning parameters or some absolute constants, but rather to
show existence of these proposing reasonable candidates for such
parameters/constants.

### Organization of the paper

The paper is organized as follows. In Section 2, we recall relevant
definitions from existing literature, and we set notations. In Section 3, (1)
we specify the model and formulate the learning problem as a regularized
group-LASSO with overlapping group norm, where the data matrix is a block-
diagonal matrix — each block of which consists of the data matrix that treats
a least-squares problem corresponding to the associated component time series;
(2) we present the learning algorithm based on the group-LASSO problem.
Section 4 contains the bulk of the technical contents, in particular, the
proof of the statements of the main results, already presented at the end of
Section 3. This Section 4 is divided into four subsections: Subsection 4.1
presents an oracle inequality bounding the one-step-ahead prediction error
(Theorem 4.2), as well as a high-probability bound on the effective noise of
the model (Theorem 4.3). Subsection 4.2 starts with restricted eigenvalue
bounds (Proposition 4.6) for the blocks of the data matrix, and goes on to
integrate the blockwise results to finally arrive at an estimate (Theorem 4.9)
of the error in estimating the AR-coefficients. Finally, combining the results
from these subsections, stability of the estimated AR model (Theorem 4.10) has
been established in Subsection 4.3. All proofs of auxiliary results are
deferred to a series of Appendices.

##  2 Preliminaries

###  2.1 General notations

  * •

$M_{d}(\mathbb{F})$ denotes the ring of $d\times d$ matrices with entries in
the field $\mathbb{F}\in\\{\mathbb{R},\mathbb{C}\\}$, and for $M\in
M_{d}(\mathbb{F})$, we write $M^{\top}$ for the transpose; if
$\mathbb{F}=\mathbb{C}$, then $M^{\star}$ denotes the conjugate transpose.

  * •

$\mathbb{D}$ is the complex closed unit disk $\mathbb{D}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{z\in\mathbb{C}:|z|\leq
1\\}$, and its boundary is $\partial\mathbb{D}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{z\in\mathbb{C}:|z|=1\\}$.

  * •

For any integer $d>0$, if $x\in\mathbb{R}^{d}$, then
$|\\!|x|\\!|_{2}=\sqrt{x^{\top}x}$, and $\mathbb{S}^{d-1}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{x\in\mathbb{R}^{d}:|\\!|x|\\!|_{2}=1\\}$;
in particular, under standard identification $\mathbb{C}=\mathbb{R}^{2}$, we
have $\mathbb{S}^{1}=\partial\mathbb{D}$.

  * •

For integer $n>0$, we will denote the set $\\{1,2,\cdots,n\\}$ by $[n]$.

  * •

For a vector
$\hat{\boldsymbol{\beta}}_{m}=(\hat{\boldsymbol{\beta}}_{m,1},\cdots,\hat{\boldsymbol{\beta}}_{m,L})$
and an integer $0<L_{0}\leq L$, we write

$\hat{\boldsymbol{\beta}}_{m}(L_{0})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\hat{\boldsymbol{\beta}}_{m,1},\cdots,\hat{\boldsymbol{\beta}}_{m,L_{0}})\,.$
(3)

###  2.2 Notations for autoregressive process

  1. Conventions:

We use $X,Y,Z,\dots$ to denote random variables, and
$\boldsymbol{X},\boldsymbol{Y},\boldsymbol{Z},\dots$ to denote random vectors.
On the other hand, we use $a,b,c,\dots$ to denote real (or complex) constants,
and $\boldsymbol{a},\boldsymbol{b},\boldsymbol{c},\dots$ for vector-valued
constants.

  2. Notations:

     * •

For $d$-dimensional autoregressive process

$\displaystyle X_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=A_{1}X_{t-1}+\cdots+A_{L}X_{t-L}+U_{t}\,,$
(4)

and $z\in\mathbb{C}$, we write

$\displaystyle{\bf\mathcal{A}}_{z}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=I-A_{1}z-\cdots-
A_{L}z^{L}\,.$ (5)

     * •

$L$ will denote an initially determined “ad-hoc” upper-bound on the true lag
of the process in (4), and $L_{0}$ the true lag.

######  Definition 2.1 (Weak stationarity).

A $d$-dimensional time series $\\{X_{t}\\}_{t\in\mathbb{Z}}$ is said to be
weakly stationary if the following holds:
$a)~{}\mathbb{E}[|\\!|X_{t}|\\!|^{2}_{2}]<\infty$ for all $t\in\mathbb{Z}$,
$b)~{}\mathbb{E}[X_{t}]=\mu$ for all $t\in\mathbb{Z}$, and $c)$
$\mathbb{E}[X_{t}X_{t-h}^{\top}]=\Gamma(h)$ for all $t,h\in\mathbb{Z}$.

######  Definition 2.2 (Strong stationarity).

A $d$-dimensional time series $\\{X_{t}\\}_{t\in\mathbb{Z}}$ is said to be
strongly stationary if for each integer $n>0$, and all integers
$t_{1},\cdots,t_{n},h$, the distributions of the vectors
$(X_{t_{1}},\cdots,X_{t_{n}})$ and $(X_{t_{1}+h},\cdots,X_{t_{n}+h})$ are
identical.

######  Definition 2.3 (Autoregressive time series).

A $d$-dimensional time series $\\{X_{t}\\}_{t\in\mathbb{Z}}$ is autoregressive
of lag at most $L>0$ if there are $d\times d$ matrices $A_{1},\cdots,A_{L}$
such that

$\displaystyle X_{t}=A_{1}X_{t-1}+\cdots+A_{L}X_{t-L}+U_{t}\,.$ (6)

holds for all $t\in\mathbb{Z}$, for some random white noise process
$\\{U_{t}\\}_{t\in\mathbb{Z}}$.

Associated to each $d$-dimensional lag-$L$ autoregressive process
$\\{X_{t}\\}_{t\in\mathbb{Z}}$ — as in equation (6) — is the associated
order-1 process ${\bf X}_{t}={\bf A}{\bf X}_{t-1}+{\bf U}_{t}$, where

$\displaystyle{\bf A}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\begin{pmatrix}A_{1\rightarrow
L},A_{L}\\\ {\bf I}_{dL-d}\,,{\bf 0}\end{pmatrix}\,,\hskip 28.45274pt{\bf
U}_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\begin{pmatrix}U_{t}\\\
{\bf 0}\end{pmatrix}\,,$ (7)

and $A_{1\rightarrow L}$ is the block matrix $(A_{1}\cdots A_{L-1})$.

######  Definition 2.4 (Stability).

A $d$-dimensional lag-$L$ autoregressive process
$\\{X_{t}\\}_{t\in\mathbb{Z}}$ — as in equation 6 — is said to be stable if
$\det({\bf I}-{\bf A}z)\neq 0$ for $|z|\leq 1$. Equivalently, the process is
stable if $\det({\mathcal{A}}_{z})\neq 0$ for $|z|\leq 1$.

######  Definition 2.5 (Reverse characteristic polynomial).

The polynomial $\det({\bf\mathcal{A}}_{z})$ is called the reverse
characteristic polynomial of the process in equation 6.

We note the equality $\det({\bf I}-{\bf A}z)=\det({\bf\mathcal{A}}_{z})$. In
particular, the process in equation (6) is stable if and only if every
eigenvalue of ${\bf A}$ is inside the open unit disk.

######  Definition 2.6 ($\epsilon$-stability).

A stable autoregressive process, as in equation (6), is said to be
$\epsilon$-stable for an $\epsilon\in(0,1)$ if the following holds:

$\displaystyle\epsilon\leq\min_{|z|=1}|\det({\bf\mathcal{A}_{z}})|\leq\max_{|z|=1}|\det({\bf\mathcal{A}_{z}})|\leq\epsilon^{-1}\,.$
(8)

######  Remark.

By maximum modulus principle, this is equivalent to saying that

$\epsilon\leq\min_{|z|\leq 1}|\det({\bf\mathcal{A}_{z}})|\leq\max_{|z|\leq
1}|\det({\bf\mathcal{A}_{z}})|\leq\epsilon^{-1}\,.$

The lower-bound here is a convenient quantification of the notion of
stability, which demands that $\min_{|z|\leq
1}|\det({\bf\mathcal{A}_{z}})|>0$. Note that we also require the upper bound
to derive our statistical guarantees.

A well-known fact about autoregressive processes is the following: see
Lütkepohl (2005, proposition 2.1) for details.

######  Lemma 2.7 (Stability implies weak stationarity).

A stable autoregressive process is weakly stationary.

##  3 Statistical Model and Estimator

This section introduces our statistical model and estimator, and presents an
algorithm to learn the parameters of the model from observed samples.

###  3.1 Statistical Model

We start with the model. Suppose that we observe time-samples generated by $M$
univariate autoregressive process, for which we know a uniform upper-bound $L$
of the true lag-order. Then, we can aggregate these $M$ univariate lag at most
$L$ autoregressive processes

$\displaystyle X^{1}_{t}~{}$
$\displaystyle=~{}\beta_{1}^{1}X_{t-1}^{1}+\cdots+\beta_{L}^{1}X_{t-L}^{1}+U^{1}_{t}\,;$
$\displaystyle~{}\,\vdots$ (9) $\displaystyle X^{M}_{t}~{}$
$\displaystyle=~{}\beta_{1}^{M}X_{t-1}^{M}+\cdots+\beta_{L}^{M}X_{t-L}^{M}+U^{M}_{t}\,.$

In this paper, we work under the simplified assumption that the true lag of
all the $M$ component processes is identical, namely, $L_{0}$, and that $L\geq
L_{0}$ is generic; neither $L_{0}$ nor $L$ is known _a priori_. Additionally,
we assume mean-zero, Gaussian white-noise innovations; that is, for each
$m\in[M]$ the set $\\{U^{m}_{t}\\}_{t\in\mathbb{Z}}$ consists of independent
mean-zero, univariate Gaussians with coordinate-wise standard deviation
$\sigma_{m}\in(0,\infty)$. Additionally, we assume that for each
$t\in\mathbb{Z}$, the noise variables $U^{1}_{t},\dots,U^{M}_{t}$ are
independent. We summarize the parameters of the model in a matrix
$\Theta\in\mathbb{R}^{M\times L}$ via $\Theta_{ml}\,\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\beta_{l}^{m}$
to refer to groups of parameters more easily later. Our goal is 1. to estimate
the parameters of all models simultaneously; and 2. to assess the lag $L_{0}$,
which is assumed to be the same over all $M$ processes.

We make the following assumption on the absolute value of the smallest
$\beta$-coefficient, which is widely known as _$\beta$ -min assumption_ in the
LASSO literature; see, for example, Bunea (2008). We note that this assumption
will only be needed to achieve the bound in Theorem 3.1 after
$\lambda$-thresholding; in particular, when no thresholding is employed, the
analysis in this paper does not require the assumption.

######  Assumption 1 (_$\beta$ -min assumption_).

There is an absolute constant $c_{\beta}>0$ such that the true autoregressive
coefficient vector $\boldsymbol{\beta}$ satisfies

$\boldsymbol{\beta}^{m}_{j}\neq
0~{}\Rightarrow~{}\boldsymbol{\beta}^{m}_{j}\geq c_{\beta}\,.$ (10)

Broadly speaking, this assumption ensures that the non-zero coefficients can
be detected in the first place. We now set out to define a regularizer. Let
$n_{j}$ denote the total number of samples from

$X_{t}^{j}=\beta^{j}_{1}X^{j}_{t-1}+\cdots+\beta^{j}_{L}X^{j}_{t-L}+U^{j}_{t}\,,$

and $T_{j}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=n_{j}-L$.
The main idea is as follows. Suppose that $L\geq L_{0}$ is some integer, and
for each $t\in\\{-L+1,\dots,1,\dots,T_{m}\\}$ and $m\in\\{1,\dots,M\\}$, we
have an observation $x_{t}^{m}$ of $X_{t}^{m}$. Denote
$\boldsymbol{\beta}_{m}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\beta^{m}_{1},\dots,\beta^{m}_{L})^{\top}$
for each $m\in[M]$, and let $\boldsymbol{\beta}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\boldsymbol{\beta}_{1}^{\top},\dots,\boldsymbol{\beta}_{m}^{\top})^{\top}$.
We define $\mathcal{G}_{1},\dots,\mathcal{G}_{L}\subset S_{ML}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{1,2,\dots,M\\}\times\\{1,2,\dots,L\\}$
by

$\mathcal{G}_{l}~{}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{1,2,\dots,M\\}\times\bigl{\\{}l,l+1\dots,L\bigr{\\}}$
(11)

for all $l\in\\{1,\dots,L\\}$. The groups are nested:
$\mathcal{G}_{1}\supset\dots\supset\mathcal{G}_{L}$. Let
$\boldsymbol{\beta}_{{\mathcal{G}}_{l}}\in\mathbb{R}^{M\times(L-l+1)}$ be the
submatrix of $\boldsymbol{\Theta}$, consisting of columns having index larger
or equal to $l$. We set the group norm to be

$\displaystyle|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}~{}$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sum_{l=1}^{L}\sqrt{M(L-l+1)}~{}|\\!|\boldsymbol{\beta}_{{\mathcal{G}}_{l}}|\\!|_{\mathbb{F}}\,,$
(12) $\displaystyle\mbox{where}\hskip
28.45274pt|\\!|\boldsymbol{\beta}_{\mathcal{G}_{l}}|\\!|_{\mathbb{F}}~{}$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sqrt{\sum_{m=1}^{M}\sum_{j=l}^{L}|\beta_{j}^{m}|^{2}}\,.$

is the Frobenius norm of $\boldsymbol{\beta}_{\mathcal{G}_{l}}$. We will
alternatively write ${\mathcal{N}}(\boldsymbol{\beta})$ for the group norm
$|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}$, for the sake of notational ease.

The overall post-sample size is denoted

$D\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=T_{1}+\cdots+T_{M}\,.$

In order to estimate the coefficient vector $\boldsymbol{\beta}$, we propose
solving the following constrained convex program

minimize
$\displaystyle\frac{1}{D}\sum_{m=1}^{M}\sum_{t=1}^{T_{m}}\bigl{(}x_{t}^{m}-\beta^{m}_{1}x_{t-1}^{m}-\cdots-\beta^{m}_{L}x_{t-L}^{m}\bigr{)}^{2}+\lambda|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}\,,$
(13)

with an appropriate tuning parameter $\lambda>0$.

The objective function can be put in a concise form. For this, we define the
vector $\boldsymbol{y}\in\mathbb{R}^{D}$, the matrix
$X\in\mathbb{R}^{D\times(ML)}$, and the parameter
$\boldsymbol{\beta}\in\mathbb{R}^{ML}$, as follows:

$\displaystyle\boldsymbol{y}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(x_{1}^{1},\dots,x_{T_{1}}^{1},x_{1}^{2},\dots,x_{T_{2}}^{2},\cdots,x_{1}^{M},\dots,x_{T_{M}}^{M})^{\top}\,;$
$\displaystyle X$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\begin{pmatrix}x_{1-1}^{1},\dots,x_{1-L}^{1}\\\
\vdots\\\ x_{T_{1}-1}^{1},\dots,x_{T_{1}-L}^{1}\\\
&x_{1-1}^{2},\dots,x_{1-L}^{2}\\\ &\vdots\\\
&x_{T_{2}-1}^{2},\dots,x_{T_{2}-L}^{2}\\\ &&\ddots\\\
&&&x_{1-1}^{M},\dots,x_{1-L}^{M}\\\ &&&\vdots\\\
&&&x_{T_{M}-1}^{M},\dots,x_{T_{M}-L}^{M}\\\ \end{pmatrix}\,;$ (14)
$\displaystyle\boldsymbol{\beta}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\beta_{1}^{1},\dots,\beta_{L}^{1},\beta^{2}_{1},\dots,\beta^{2}_{L},\cdots,\cdots,\beta_{1}^{M},\dots,\beta_{L}^{M})^{\top}\,.$

It is immediate that

$\sum_{m=1}^{M}\sum_{t=1}^{T_{m}}\bigl{(}x_{t}^{m}-\beta^{m}_{1}x_{t-1}^{m}-\cdots-\beta^{m}_{L}x_{t-L}^{m}\bigr{)}^{2}~{}=~{}|\\!|\boldsymbol{y}-X\boldsymbol{\beta}|\\!|_{2}^{2}\,.$

In conclusion, the above estimation program is equivalent to

$\displaystyle\widehat{\boldsymbol{\beta}}~{}\in~{}\operatornamewithlimits{argmin}_{\boldsymbol{\beta}\in\mathbb{R}^{ML}}\biggl{\\{}\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\beta}|\\!|_{2}^{2}+\lambda|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}\biggr{\\}}\,.$
(15)

Hence, the estimator can be cast as a modified group-lasso estimator, which
means that we can use established group-lasso algorithms that allow for
overlapping groups (Mairal et al., 2011). In essence, the above estimator
generalizes the _elementwise_ estimator $\mbox{HLag}^{\text{E}}$ in Nicholson
et al. (2020) to multiple time series. Note that, the ordinary LASSO
estimator—as well as any group-LASSO estimators with non-overlapping
groups—enforces sparsity by setting coefficients to zero without paying heed
to the fact that when only an upper-bound to the true-lag $L_{0}$ is an input
to the regression— _all_ coefficients indexed between $L_{0}+1$ and $L$ are
supposed to be zero before any coefficient with index smaller or equal to
$L_{0}$; however, the penalty obtained via the chained groups
$\mathcal{G}_{1}\supseteq\dots\supseteq\mathcal{G}_{L}$ precisely achieves
this feat.

######  Example (Regularizer).

We consider the case of two univariate lag (at most) three autoregressive
processes; that is, $M=2$ and $L=3$. The corresponding groups are the
following:

$\displaystyle{\mathcal{G}}_{1}$
$\displaystyle=\\{(1,1),(1,2),(1,3),(2,1),(2,2),(2,3)\\}\;;$
$\displaystyle{\mathcal{G}}_{2}$
$\displaystyle=\\{(1,2),(1,3),(2,2),(2,3)\\}\;;$
$\displaystyle{\mathcal{G}}_{3}$ $\displaystyle=\\{(1,3),(2,3)\\}\,.$

Thus, the group norm of $\boldsymbol{\beta}\in\mathbb{R}^{2\times 3}$ is

$\displaystyle|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}$
$\displaystyle=\sqrt{6}\sqrt{(\beta^{1}_{1})^{2}+(\beta^{2}_{1})^{2}+(\beta^{1}_{2})^{2}+(\beta^{2}_{2})^{2}+(\beta^{1}_{3})^{2}+(\beta^{2}_{3})^{2}}$
$\displaystyle\hskip
71.13188pt+~{}\sqrt{4}\sqrt{(\beta^{1}_{2})^{2}+(\beta^{2}_{2})^{2}+(\beta^{1}_{3})^{2}+(\beta^{2}_{3})^{2}}$
$\displaystyle\hskip
142.26378pt+~{}\sqrt{2}\sqrt{(\beta^{1}_{3})^{2}+(\beta^{2}_{3})^{2}}\,.$

Notice that, when the (regularized) LASSO sets a certain group (say the second
group above) to $\boldsymbol{0}$, it automatically sets all the following
groups to $\boldsymbol{0}$ as well. More specifically, when the (regularized)
LASSO sets $\boldsymbol{\beta}_{{\mathcal{G}}_{l}}\boldsymbol{0}$, then the
hierarchical structure
${\mathcal{G}}_{l}\supseteq{\mathcal{G}}_{l+1}\supseteq\cdots$ means that for
all $r>0$, each coordinate of $\boldsymbol{\beta}_{{\mathcal{G}}_{l+r}}$ comes
as a coordinate of $\boldsymbol{\beta}_{{\mathcal{G}}_{l}}$, thus ensuring
$\boldsymbol{\beta}_{{\mathcal{G}}_{l+r}}=\boldsymbol{0}$ for each $r\geq 0$.

In what follows, we use the following notations. For $m\in[M]$ and $l\in[L]$,
and $t\leq T_{m}$, we write

$\displaystyle\boldsymbol{U}^{(m)}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(U^{m}_{1},\cdots,U^{m}_{T_{m}})^{\top}\,;$
(16) $\displaystyle X^{(m,l)}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(X^{m}_{1-l},\dots,X^{m}_{T_{m}-l})^{\top}\,;$
$\displaystyle X_{t}^{(m)}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(X^{m}_{t-1},\dots,X^{m}_{t-L})\,;$
$\displaystyle X^{(m)}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(X^{(m,1)},\dots,X^{(m,L)})\,.$

Moreover, we will write $X^{(m)}_{<j}$ to denote any of the variables
$X^{m}_{j^{\prime}}$ for $j^{\prime}<j$.

###  3.2 Estimation Pipeline

Learning autoregressive coefficients of multiple time series of potentially
different lengths and identical true lag is more complex than just the usual
group-LASSO problem, where the groups form a partition of the index set. From
a methodological perspective, some immediate technical challenges are 1.
deciding on what $L$ should be taken in the formulation of the convex problem
(13), 2. how to disentangle the dual of the group norm (in order to apply
Hölder’s inequality to derive oracle prediction guarantees as in subsection
4.1 below); from a practical perspective, the challenge lies in incorporating
the varying number of samples into the convex problem.

We now give a high-level overview of our estimation pipeline (described below)
that takes as input the multiple time series, and forms the appropriate convex
problem in the form of (15), and solves this convex problem via stochastic
proximal gradient method as in Nicholson et al. (2020). In essence, the idea
is to start looking into the data set to first find the samples corresponding
to the component which has the minimum number of samples (breaking ties
arbitrarily). We then use this component to find the initial input lag, $L$,
as described in (38). In the next step, we solve the regularized least-squares
problem in (13) — where the penalty function is the group norm $\mathcal{N}$,
discussed further below (see (12)). The rate of convergence of the procedure
is quadratic in number of computational steps as discussed in Nicholson et al.
(2020).

Input : samples from the component stable processes, confidence parameters
$A\geq 1$ and $\delta>0$, stability parameter $\epsilon\in(0,1)$

Output : estimated lag $\hat{L}_{0}$, and estimated autoregression coefficient
vector $\tilde{\boldsymbol{\beta}}$

Let $n_{\min}$ be the minimum number of samples from the components. Solve for
$L$ (see equation (38)):
$n_{\min}=L+84Ae\zeta^{-2}L\log\left(\frac{ML}{\delta}\right)\,.$

Use the learning algorithm in Nicholson et al. (2020) as subroutine to solve
the convex problem in (15) — with $L$ as above; call the output
$\hat{\boldsymbol{\beta}}$.

If ${\boldsymbol{\beta}}_{m}={\boldsymbol{\beta}}_{m^{\prime}}$ for all
$m,m^{\prime}\in[M]$ (equivalently, the component time series are from
identical AR process), return $\hat{L}_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\max\\{j:|(\hat{\boldsymbol{\beta}}_{1})_{j}|>\lambda\\}$
and $\tilde{\boldsymbol{\beta}}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\hat{\boldsymbol{\beta}}^{\prime}_{0},\cdots,\hat{\boldsymbol{\beta}}^{\prime}_{0})$,
where (refer Equation (3) for notation)
$\hat{\boldsymbol{\beta}}_{0}^{\prime}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{M}\sum_{m=1}^{M}\hat{\boldsymbol{\beta}}_{m}(\hat{L}_{0})\,;$
else, return  $\displaystyle\hat{L}_{0}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\max\\{j:|\hat{\boldsymbol{\beta}}_{j}^{m}|>\lambda~{}\mbox{for
some}~{}m\in[M]\\}\,,$ $\displaystyle\text{and}\hskip
28.45274pt\tilde{\boldsymbol{\beta}}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\hat{\boldsymbol{\beta}}_{1}(\hat{L}_{0}),\cdots,\hat{\boldsymbol{\beta}}_{m}(\hat{L}_{0}))\,.$

Algorithm 1 AR Coefficient Estimation Pipeline

The following is the main theorem on the theoretical properties of the output
of the estimation pipeline. This theorem is essentially a summary of the
results contained in Section 4, and will be discussed in all detail in
subsections 4.2.2 and 4.3, as indicated below.

######  Theorem 3.1 (Main Theorem).

Let $\hat{\boldsymbol{\beta}}$ be the output of the group-LASSO in (15) above,
with

$\lambda={24(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\,,$

where $\zeta=6^{-3}\epsilon^{4}$, and $L\geq L_{0}$ satisfying

$n_{\min}=L+84Ae\zeta^{-2}L\log\left(\frac{ML}{\delta}\right)\,.$

Suppose that the $\beta$-min condition (10) holds with $c_{\beta}=\lambda$. If
the total number $D$ of post-samples satisfies

$D\geq
3^{9}\cdot(84Ae)C_{\sharp}^{5}(1+\epsilon^{-2}+\epsilon^{-4})^{2}\left(\frac{\sigma_{\max}}{\sigma_{\min}}\right)^{4}(\epsilon^{3}\zeta)^{-2}M^{a}L_{0}^{2}L^{3}\log\left(\frac{ML}{\delta}\right)\log(2L)\,,$

and if $T_{\min}\geq 84eA\zeta^{-2}L_{0}\log L$, then the following holds with
high probability:

  1. 1.

the estimation error is given by

$|\\!|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}|\\!|_{2}\leq\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}\,;$

  2. 2.

if $S_{\lambda}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{j:|\hat{\boldsymbol{\beta}}_{j}|>\lambda\\}$,
then the false discovery is bounded by the following inequality:

$|S_{\lambda}\setminus\operatorname{supp}({\boldsymbol{\beta}})|~{}\leq\frac{243(84Ae)^{\frac{1}{2}}LL_{0}^{\frac{3}{2}}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}\lambda}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}\,.$

  3. 3.

the AR-models — fitted with coefficients $\hat{\boldsymbol{\beta}}_{0}$
returned by the Algorithm 1 (”AR Coefficient Estimation Pipeline”). — are
stable, with high probability.

###### Proof of Theorem 3.1.

Subject to the stated $\beta$-min condition (10), this follows immediately
from Theorem (4.9) in combination with Theorem (4.10). ∎

In the pipeline above, it is necessary to consider two distinct stability
parameters $A\geq 1$ and $\delta$, as it is not possible to integrate them
into a single parameter — due mainly to the different number of samples from
the component processes in our set-up. Also, we separately mention the two
cases (of $\boldsymbol{\beta}_{m}$’s being identical (or not) for all $m$) in
order to specifically emphasize that in the first case, our algorithm requires
a smaller number of samples than in the later case (which allows savings of a
factor of $M$ in the sample complexity).

Decent algorithms like the proximal gradient method used in the subroutine
above might produce small non-zero valued parameters as numerical artifacts.
One could consider other types of algorithms instead, but the required number
of computational steps could be much higher. Moreover, those artifacts, and
statistical false positives more generally, can also be controlled by standard
$\lambda$-thresholding as mentioned in the algorithm.

##  4 Statistical Guarantees

This section contains the main theoretical results of this paper. We begin
with deriving bounds for the one-step ahead prediction error, first formulated
by an oracle inequality (see Theorem 4.2), which depends on the tuning
parameter $\lambda$ of our least-squares penalisation approach. Then we
control this tuning parameter by controlling the effective noise of our Lasso-
optimisation problem. Both things together will finally yield more explicit
rates of our one-step ahead prediction error. The second part of this section
treats control of the estimated autoregressive coefficients, including control
of false discovery for their support (Theorem 4.9). In the end we present our
result on stability of the estimated AR-model (Theorem 4.10).

To start with, we briefly recall the notion of the dual norm of our group-
LASSO norm $\mathcal{N}$ defined above.

Our underlying space is $\mathbb{R}^{ML}$. Observe that

$\boldsymbol{\beta}\mapsto{\mathcal{N}}(\boldsymbol{\beta})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sum\sqrt{|{\mathcal{G}_{l}}|}\cdot|\\!|\boldsymbol{\beta}_{{\mathcal{G}}_{l}}|\\!|_{2}$

is a norm for any
${\mathcal{G}}=\\{{\mathcal{G}}_{1},\cdots,{\mathcal{G}}_{l}\\}$ that covers
$[ML]$. However, this norm is singular at any point where all the coordinates
in a group vanish. Moreover, all the coordinates in all of the smaller-sized
groups will vanish, which is important in our analysis, since we will
basically care about sparse solutions. Thus, we can not appeal to differential
techniques to get bounds on the norm, but rather need to appeal to the dual
norm approach.

######  Definition 4.1 (Dual norm).

For a norm ${\mathcal{N}}$ on $\mathbb{R}^{ML}$, the dual norm
${\mathcal{N}}_{\star}(\boldsymbol{\alpha})$ of
$\boldsymbol{\alpha}\in\mathbb{R}^{ML}$ is the optimum solution of the
following convex program:

maximize $\displaystyle\langle\boldsymbol{\alpha},\boldsymbol{\beta}\rangle$
subject to $\displaystyle\mathcal{N}(\boldsymbol{\beta})\leq 1\,.$

The dual norm is used to encapsulate the effective noise of LASSO. More
explicitly, the dual norm shows up in the form
${\mathcal{N}}_{\star}(X^{\top}\boldsymbol{U})$, called the _effective noise_
of the LASSO in equation (15). Recall that the effective noise vector
$X^{\top}\boldsymbol{U}$ can be thought of as the “projection” of the noise
vector $\boldsymbol{U}$ on the column space of $X$, and thus,
${\mathcal{N}}_{\star}(X^{\top}\boldsymbol{U})$ is a measure of the “true”
noise present in the data.

In Appendix D, we obtain a generic bound on the dual norm, which will be used
to obtain the statistical guarantees of this paper.

###  4.1 Prediction Error

Here we now give the announced theoretical results on non-asymptotic bounds
for the one-step ahead prediction error, first formulated by an oracle
inequality (see Theorem 4.2), then in the following subsection including
control of the tuning parameter $\lambda$ by control of the effective noise of
our Lasso-optimisation problem. This enables us to formulate concrete rates of
the prediction error. Note that this approach delivers an explicit way of how
to select $L$ (via equation (24), and subsequently (38)) the input parameter
for Step 2 of Algorithm 1 (”AR Coefficient Estimation Pipeline”).

####  4.1.1 Oracle Prediction Error

The problem (15) has a solution because, given any specific realization of the
time series (thus, effectively, fixing $\boldsymbol{y}$ and $X$) and any
$\lambda>0$, the convex function

$\displaystyle
f_{\lambda}\left(\boldsymbol{\beta}\right)=\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\beta}|\\!|_{2}^{2}+\lambda|\\!|\boldsymbol{\beta}|\\!|_{\mathcal{G}}$
(17)

is continuous, and

$\lim_{|\\!|\boldsymbol{\beta}|\\!|_{2}\rightarrow\infty}f_{\lambda}(\boldsymbol{\beta})=\infty\
.$

This also shows that we can consider this as a convex program on a compact
domain, because we should (at least in theory) be able to restrict the domain
of minimization to be an $\ell^{2}$-ball of suitable radius, say
$R_{r,X,\boldsymbol{y}}\in(0,\infty)$. Let $\hat{\boldsymbol{\beta}}$ be a
solution of this convex program.

Now write the time series observations in the matrix form as $\mathcal{X}$. We
are interested in an estimate of the ’risk’, equivalently, the in-sample, one-
step-ahead mean squared forecast error
$\mathbb{E}[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}/D\mid\mathcal{X}]$.
In the derivations below, we follow a well-known approach for its control, as
appeared (for example) in (Lederer, 2021, Chapter 6). We defer the proof to
Appendix A.

######  Theorem 4.2 (Prediction Guarantee).

Suppose that
$\lambda\geq\frac{2}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})$, where
$\mathcal{N}(\boldsymbol{\alpha})=\sum_{l=1}^{L}\sqrt{|\mathcal{G}_{l}|}\cdot|\\!|\boldsymbol{\alpha}_{(\geq
l)}|\\!|_{2}$ as in Proposition D.3. Write
$\overline{\sigma}^{2}=D^{-1}(T_{1}\sigma_{1}^{2}+\cdots+T_{M}\sigma_{M}^{2})$;
then

$\displaystyle\frac{1}{D}\mathbb{E}\left[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]$
$\displaystyle\leq$
$\displaystyle\overline{\sigma}^{2}+\min_{\alpha\in\mathbb{R}^{ML}}\left(\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+2\lambda\mathcal{N}(\boldsymbol{\alpha})\right),$

In particular, the following inequality holds:

$\displaystyle\mathbb{E}\left[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]-\mathbb{E}\left[|\\!|\boldsymbol{y}-X{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]$
$\displaystyle\leq$
$\displaystyle\min_{\alpha\in\mathbb{R}^{ML}}\left(|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+2D\lambda\mathcal{N}(\boldsymbol{\alpha})\right)\,.$

If, moreover,
$\lambda\geq\frac{4}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})$, then

$\displaystyle\frac{1}{D}\mathbb{E}\left[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]$
$\displaystyle\leq\overline{\sigma}^{2}+\frac{\lambda}{2}\min\left\\{3\mathcal{N}(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}),3\mathcal{N}(\boldsymbol{\beta})-\mathcal{N}(\hat{\boldsymbol{\beta}})\right\\}.$
(18)

Consequently,

$\displaystyle\frac{1}{D}\mathbb{E}\left[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]\leq\overline{\sigma}^{2}+2\lambda\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}-\boldsymbol{\hat{\beta}}|\\!|_{{\mathcal{G}}_{\ell}}\
.$ (19)

This result is in the form of standard oracle inequalities in high-dimensional
statistics (Lederer, 2021, Chapter 6). It shows that the estimator minimizes
the one-step-ahead mean-squared forecast risk up to a complexity term that is
linear in the tuning parameter and the model complexity. Hence, the above
bound yields an upper bound for the rate convergence once we can control the
tuning parameter $\lambda$ via an upper bound on the effective noise.

Thanks to Theorem 4.2 and Proposition D.3, in order to find the smallest
tuning parameter $\lambda$ fulfilling $\lambda\geq
4\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})/D$ it suffices to derive a high-
probability upper-bound on
$D^{-1}L^{-\frac{1}{2}}|\\!|X^{\top}\boldsymbol{U}|\\!|_{\infty}$ (which is
precisely the bound on the dual norm of $X^{\top}U$).

####  4.1.2 Control of the Effective Noise for bounding the tuning parameter

We can prove the following high-probability bound on the tuning parameter
(equivalently, on the effective noise). Again, its proof appears in Appendix
B. Note that this is a major inequality, and while the arguments are well-
known, we have applied those arguments to the case where the sparsity
enforcing regularizer is induced by the overlapping group norm.

######  Theorem 4.3 (Bound on the Effective Noise).

Let $C_{\sharp}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=T_{\max}/T_{\min}$.
For any $\eta>0$ satisfying

$\eta\geq\frac{8C_{\sharp}\sigma^{2}_{\max}(1+\epsilon^{-2}+\epsilon^{-4})}{M}\,,$
(20)

and any $\delta>0$, if $D=T_{1}+\cdots+T_{m}$ satisfies

$D\geq\frac{8\sigma_{\max}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,,$
(21)

then the following inequality holds:

$\mathbb{P}\left[\frac{2}{D}{\mathcal{N}}_{\star}(X^{\top}\boldsymbol{U})\geq\frac{3\eta}{2\sqrt{L}}\right]~{}\leq~{}\delta\,.$
(22)

Here $c_{0}>0$ is the absolute constant from the Gaussian concentration
inequality in proposition E.1.

The interpretation of the above result is that for the LASSO oracle inequality
(Theorem 4.2) to hold with high probability, it is sufficient to choose
$\lambda$ to be just as large as $(3\eta)/(2\sqrt{L})$, but it is not
necessary to take it larger. Indeed, equation (22) shows that the probability
that $2/D\mathcal{N}_{\star}(X^{T}U)$ is larger than $3\eta/2\sqrt{L}$ is
small; thus, we may assume $2/D\mathcal{N}_{\star}(X^{T}U)<3\eta/2\sqrt{L}$,
to hold with probability $1-\delta$.

The above result bounds the tails of the effective noise. For any $A\geq 1$,
we will now set

$\eta=~{}C_{0}\sqrt{A}C_{\epsilon}C_{\sharp}^{\frac{3}{2}}\sigma_{\max}^{2}\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\,,$

for some absolute constant $C_{0}>0$ and parameter $C_{\epsilon}>0$ that
depends only on the stability parameter $\epsilon$. Henceforth, we take

$\eta\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=8(84Ae)^{\frac{1}{2}}\zeta^{-1}(1+\epsilon^{-2}+\epsilon^{-4})C_{\sharp}^{\frac{3}{2}}\sigma_{\max}^{2}\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\,,$
(23)

which is obtained by plugging-in values of $C_{0}$ and $C_{\epsilon}$ — as in
the proof of Lemma 4.4 below — into Equation (23); we do not attempt to
optimize these constants.

######  Lemma 4.4 (Data-dependent selection of $L$).

There is an absolute constant $C\in(0,\infty)$ and a parameter
$C_{\epsilon}>0$ that depends only on the stability parameter $\epsilon>0$
such that if $\eta$ is as in (23) and

$T_{\min}\leq 84Ae\zeta^{-2}L\log\left(\frac{ML}{\delta}\right)\,,$ (24)

then the inequality (20) holds.

This results shows that there is a suitable $\eta$ for our theories to hold.
The question of finding such an $\eta$ in practice will need to be discussed
in more applied future work.

###### Proof.

We set $C_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(12)^{3}\sqrt{84e}$,
and

$C_{\epsilon}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\epsilon^{-4}(1+\epsilon^{-2}+\epsilon^{-4})\,.$

Note that

$C_{\sharp}T_{\min}=T_{\max}\geq~{}\frac{D}{M}\,.$ (25)

Now, with $\eta$ as in (23), the inequality (20) is ensured by the following
sequence of inequalities:

$\displaystyle(84Ae)^{\frac{1}{2}}\zeta^{-1}\sqrt{{L\log\left(\frac{ML}{\delta}\right)}}$
$\displaystyle\geq\sqrt{T_{\min}}$ $\displaystyle\Rightarrow\hskip 14.22636pt$
$\displaystyle(84Ae)^{\frac{1}{2}}\zeta^{-1}C_{\sharp}^{\frac{1}{2}}\sqrt{{L\log\left(\frac{ML}{\delta}\right)}}$
$\displaystyle\geq\sqrt{\frac{D}{M}}$ $\displaystyle\Rightarrow\hskip
14.22636pt$ $\displaystyle
8(84Ae)^{\frac{1}{2}}\zeta^{-1}(1+\epsilon^{-2}+\epsilon^{-4})C_{\sharp}^{\frac{3}{2}}\sigma_{\max}^{2}\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}$
$\displaystyle\geq\frac{8C_{\sharp}\sigma^{2}_{\max}(1+\epsilon^{-2}+\epsilon^{-4})}{M}\,.$

∎

The inequality (24) provides for the theoretical support of our choice of $L$
in formulating the penalized least squares program in (15).

The explicit nature of the parameters in the proof above is crucial here, in
order to satisfy both (20) and (21) above, as well as remaining compatible
with the requirements involving the restricted eigenvalue property (as in
Corollary 4.7) below. Note that the current LASSO-based literature often
ignores combining the requirements coming from standard oracle inequality type
result (Theorem 4.2) and the restricted eigenvalue type results.

######  Remark (On choosing $L$ via equation (24)).

The use of such an upper-bound $L$ conforms with the (by now) well-known
restricted isometry property of sub-sampled Gaussian matrices (that is,
matrices whose entries are iid Gaussian) in the compressed-sensing literature;
for more details, see (for example) (Candes and Tao, 2006, section 1.E) and
the references therein.

For $D$ as in (21), this gets us the following bound:

$\displaystyle\mathbb{P}\left[{\mathcal{N}}_{\star}\left(\frac{2}{D}X^{\top}\boldsymbol{U}\right)>{12(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\right]$
$\displaystyle\leq\delta,$

In view of Theorem 4.3, we correspondingly set

$\lambda={24(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}$
(26)

to force (18) to be true.

######  Remark.

From the above, we note that $\lambda$ decreases when any of $M,D$ increases
(as the function $x^{-1}\log x\rightarrow 0$ when $x\rightarrow\infty$). This
is intuitive since larger $M$ requires that LASSO must not set too many
coefficients to 0, and larger the $D$ better the non-regularized estimator is
as approximation of the true coefficients.

Together with Theorem 4.2, the above yields the following bound:

######  Corollary 4.5 (Concrete rates of Prediction Error ).

Suppose that

$\eta=~{}8(84Ae)^{\frac{1}{2}}\zeta^{-1}(1+\epsilon^{-2}+\epsilon^{-4})C_{\sharp}^{\frac{3}{2}}\sigma_{\max}^{2}\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\,,$

as set in Equation (23) above, and

$D\geq\frac{8\sigma_{\max}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,.$

Then, the following inequality holds with probability at least $1-\delta$:

$\displaystyle\mathbb{E}\left[\frac{1}{D}|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid\mathcal{X}\right]-\overline{\sigma}^{2}$
$\displaystyle\leq~{}$
$\displaystyle\min_{\alpha\in\mathbb{R}^{ML}}\left(\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+{\frac{24}{\zeta}(84Ae)^{\frac{1}{2}}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\mathcal{N}(\boldsymbol{\alpha})\right)\,.$

###  4.2 Estimating error of parameters, false discoveries

We now deliver the treatment of assertions 1. and 2. on estimation error and
support (via control of false discoveries) of our Main Theorem 3.1. For this
we need to cope with the Restricted Eigenvalue Property as it typically arises
in LASSO analysis (see Basu and Michailidis (2015)).

####  4.2.1 Control of the restricted eigenvalue property

Our prediction guarantees stated so far did not impose restrictions on the
minimum number of samples from each component. For estimation and stability
guarantees, however, we need to demand a minimum number of those samples. In
particular, we will need the following proposition, on the restricted
eigenvalue property of the Gram matrix $X^{\top}X$. The proof presented in
Appendix C follows the same lines of arguments as in Basu and Michailidis
(2015). Because of the non-uniform nature of the block dimensions of the data
matrix in (3.1), as well as the individual weights assigned to the components
of the coefficient $\boldsymbol{\beta}$, we can only hope for a block-wise
result as stated below.

The following proposition delivers a lower bound on the error-expression
$|\\!|X_{m}\boldsymbol{v}_{m}|\\!|_{2}^{2}$ with ${\boldsymbol{v}_{m}}$ the
$m-$th component of $\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}$.

######  Proposition 4.6 (Bound on the Restricted Eigenvalue).

For any confidence parameter $\delta>0$ and all vectors
$\boldsymbol{v}=(\boldsymbol{v}_{1}^{\top},\dots,\boldsymbol{v}_{M}^{\top})^{\top}\in(\mathbb{R}^{L})^{M}$,
the following inequality holds with $\zeta\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{\epsilon^{4}}{216}$
and $s_{m}$ even positive integers:

$\displaystyle\mathbb{P}\left[\displaystyle\forall~{}{m\in[M]}\inf_{\boldsymbol{v}_{m}\in\mathbb{R}^{L}}\left(|\\!|X_{m}\boldsymbol{v}_{m}|\\!|_{2}^{2}-\frac{T_{m}\sigma_{m}^{2}\epsilon^{2}}{2}\left(|\\!|\boldsymbol{v}_{m}|\\!|_{2}^{2}-\frac{2}{s_{m}}|\\!|\boldsymbol{v}_{m}|\\!|_{1}^{2}\right)\right)\geq
0\right]$ (27) $\displaystyle\geq$ $\displaystyle
1-2\sum_{m=1}^{M}\exp\left(-\frac{T_{m}\min\\{\zeta,\zeta^{2}\\}}{2}+s_{m}\min\\{\log
L,\log(21eL/s_{m})\\}\right)\,.$

A proof of Proposition 4.6 appears in Appendix C.

Setting

$s_{m}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=2\Big{\lfloor}\frac{T_{m}\zeta^{2}}{8\log
L}\Big{\rfloor}\,,$ (28)

we get the following immediate corollary to Proposition 4.6:

######  Corollary 4.7 (Bound on the Restricted Eigenvalue for specific
$s_{m}$).

If

$T_{\min}\geq 84e\zeta^{-2}\log L\,,$ (29)

where $\zeta=6^{-3}\epsilon^{4}$, then the following holds:

$\displaystyle\mathbb{P}\left[\forall{m\in[M]}\inf_{\boldsymbol{v}_{m}\in\mathbb{R}^{L}}\left(|\\!|X_{m}\boldsymbol{v}_{m}|\\!|_{2}^{2}-\frac{T_{m}\sigma_{m}^{2}\epsilon^{2}}{2}|\\!|\boldsymbol{v}_{m}|\\!|_{2}^{2}\left(1-\frac{8L\log
L}{T_{m}\zeta^{2}}\right)\right)\geq 0\right]$ $\displaystyle\geq~{}$
$\displaystyle 1-2\sum_{m=1}^{M}e^{-\frac{T_{m}\zeta^{2}}{4}}\,.$ (30)

###### Proof.

We use $|\\!|\boldsymbol{v}_{m}|\\!|_{1}^{2}\leq
L|\\!|\boldsymbol{v}_{m}|\\!|_{2}^{2}$ in Proposition 4.6. ∎

####  4.2.2 Bounds on estimation error of the autoregressive coefficients, and
false discovery

This section contains the bound on the error of estimation of the
autoregressive coefficients, and the false discovery.

The following is a compact notation used in the proof of the theorem below
(also in the proof of Theorem 4.2).

######  Definition 4.8 (Notation).

${\mathcal{N}}_{\leq L_{0}}(\boldsymbol{v})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}\,.$
(31)

The theorem below bounds the $\ell^{2}$-error in estimating the autoregressive
coefficients $\boldsymbol{\beta}$, as well as the size of the set of false
positives, using the penalized LASSO formulation as in (15).

######  Theorem 4.9 (Bounds on the Estimation Error and False Discovery).

Let $\hat{\boldsymbol{\beta}}$ be the solution to the group-regularized LASSO
in (15) with $\lambda$ as in Equation (26) above — namely,

$\lambda={24(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}$

where $\zeta=6^{-3}\epsilon^{4}$; if

$\alpha\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\min_{m\in[M]}\frac{T_{m}\sigma_{m}^{2}}{D}\,,$
(32)

and

$D\geq\frac{8\sigma_{\max}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,,$

then, for any confidence parameter $\delta>0$, the inequality

$|\\!|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}|\\!|_{2}\leq\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}$
(33)

holds with probability at least

$(1-\delta)\left(1-2\sum_{m=1}^{M}L^{-\frac{T_{m}}{4\log L}}\right)\,,$ (34)

provided $T_{\min}\geq 84eA\zeta^{-2}L_{0}\log L$ with
$\zeta=6^{-3}\epsilon^{4}$.

Moreover, introducing the notation $S_{\lambda}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{j:|\hat{\boldsymbol{\beta}}_{j}|>\lambda\\}$,
then the false discovery is bounded by the following inequality holding with
probability as in (34) above:

$|S_{\lambda}\setminus\operatorname{supp}({\boldsymbol{\beta}})|~{}\leq\frac{243(84Ae)^{\frac{1}{2}}LL_{0}^{\frac{3}{2}}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}\lambda}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}\,.$
(35)

######  Remark.

In the above, the term inside the square root decreases asymptotically, and
the terms outside may change or stay fixed (depending on how the number of
samples for each component increases).

###### Proof of Theorem 4.9.

For ease of notation, write

$\boldsymbol{v}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\,.$

By Proposition 4.6 above (which can be applied since $T_{\min}\geq
84e\zeta^{-2}\log L$), one has

$\displaystyle\frac{1}{D}|\\!|X\boldsymbol{v}|\\!|_{2}^{2}$
$\displaystyle\geq$
$\displaystyle\frac{1}{D}\sum_{m=1}^{M}|\\!|X_{m}\boldsymbol{v}_{m}|\\!|_{2}^{2}$
$\displaystyle\geq$
$\displaystyle\frac{\epsilon^{2}}{2}\sum_{m=1}^{M}\frac{\sigma_{m}^{2}T_{m}}{D}|\\!|\boldsymbol{v}_{m}|\\!|_{2}^{2}\left(1-\frac{8L\log
L}{T_{m}\zeta^{2}}\right)$ $\displaystyle>$
$\displaystyle\frac{4\alpha\epsilon^{2}}{9}|\\!|\boldsymbol{v}|\\!|_{2}^{2}\,.$

To this, we now apply (see (45) in Appendix A) the inequality

$\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\leq
3\lambda{\mathcal{N}}_{\leq
L_{0}}(\boldsymbol{\beta}-\boldsymbol{\hat{\beta}})\,,$

to obtain

$\displaystyle|\\!|\boldsymbol{v}|\\!|_{2}^{2}$ $\displaystyle\leq$
$\displaystyle\frac{27\lambda}{8\alpha\epsilon^{2}}\mathcal{N}_{\leq
L_{0}}(\boldsymbol{v})$ $\displaystyle\leq$ $\displaystyle\frac{27\lambda
L_{0}\sqrt{ML}}{8\alpha\epsilon^{2}}|\\!|\boldsymbol{v}|\\!|_{2}$
$\displaystyle\Rightarrow\hskip 28.45274pt|\\!|\boldsymbol{v}|\\!|_{2}$
$\displaystyle\leq$ $\displaystyle\frac{27\lambda
L_{0}\sqrt{ML}}{8\alpha\epsilon^{2}}$

With the choice of $\lambda$ as in equation (26), namely,

$\lambda={24(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}\,,$

it now follows that

$\displaystyle|\\!|\boldsymbol{v}|\\!|_{2}$ $\displaystyle\leq$
$\displaystyle\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}\,,$

as claimed. To get the bound on the false discovery, write $S\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\operatorname{supp}({\boldsymbol{\beta}})$;
the bound on the false discovery can be obtained as follows:

$\displaystyle|S_{\lambda}\setminus S|$ $\displaystyle=\sum_{j\in[ML]\setminus
S}1_{|\hat{\boldsymbol{\beta}}_{j}|>\lambda}(\hat{\boldsymbol{\beta}})$
$\displaystyle=\sum_{j\in[ML]\setminus
S}1_{|\hat{\boldsymbol{v}}_{j}|>\lambda}(\hat{\boldsymbol{v}})$
$\displaystyle\leq\frac{1}{\lambda}\sum_{j\in[ML]\setminus
S}|\hat{\boldsymbol{v}}_{j}|$ by remark (Remark)
$\displaystyle\leq\frac{3}{\lambda}\sum_{j\in S}|\hat{\boldsymbol{v}}_{j}|$
$\displaystyle\leq\frac{3\sqrt{L_{0}}}{\lambda}|\\!|\hat{\boldsymbol{v}}|\\!|_{2}$

which yields the bound in (35) by inequality (33). ∎

######  Remark.

For the explicit choice of $\lambda$ as in equation (26), namely,

$\lambda={24(84Ae)^{\frac{1}{2}}\zeta^{-1}\sigma_{\max}^{2}}{C_{\sharp}}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})\sqrt{\frac{L\log\left(\frac{ML}{\delta}\right)}{DM}}$

with $\zeta=6^{-3}\epsilon^{4}$, the inequality (35) becomes

$|S_{\lambda}\setminus\operatorname{supp}({\boldsymbol{\beta}})|~{}\leq\frac{81(ML)^{\frac{1}{2}}L_{0}^{\frac{3}{2}}}{8\alpha\epsilon^{2}}\,.$
(36)

###  4.3 Stability of the estimated AR model

We finally prove that the coefficients estimated as per the overlapping-group-
LASSO in (15), with $\lambda$ as in (26), lie in the region for stability of
univariate lag-$L$ autoregressive processes, even when the number of post-
samples is ’not too large’. More specifically, we have the following theorem.

######  Theorem 4.10 (Stability Guarantee).

Let $A\geq 1$ be a confidence parameter. Let $\hat{\boldsymbol{\beta}}$ be the
output of the group-LASSO in (15) above, using $D$ post-samples, where

$D\geq
3^{9}\cdot(84Ae)C_{\sharp}^{5}(1+\epsilon^{-2}+\epsilon^{-4})^{2}\left(\frac{\sigma_{\max}}{\sigma_{\min}}\right)^{4}(\epsilon^{3}\zeta)^{-2}M^{a}L_{0}^{2}L^{3}\log\left(\frac{ML}{\delta}\right)\log(2L)\
.$ (37)

If $T_{\min}\geq 84eA\zeta^{-2}L_{0}\log L$, and $L$ satisfies

$n_{\min}=L+84Ae\zeta^{-2}L\log\left(\frac{ML}{\delta}\right)\,,$ (38)

with $\zeta=6^{-3}\epsilon^{4}$, then the AR-models — fitted with coefficients
$\hat{\boldsymbol{\beta}}_{0}$ returned by Algorithm 1 (”AR Coefficient
Estimation Pipeline”) — are stable, with probability as in (34) — in the
following scenarios:

  1. 1.

when all the time series are different realizations of a unique underlying
stochastic process, $a\geq 1$, and

$\hat{\boldsymbol{\beta}}_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{M}\sum_{m=1}^{M}\hat{\boldsymbol{\beta}}_{m}\,;$

  2. 2.

when all the time series are realizations of different underlying stochastic
processes (equivalently, all the $\boldsymbol{\beta}_{m}$’s are different),
$a\geq 2$, and $\hat{\boldsymbol{\beta}}_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\hat{\boldsymbol{\beta}}$.

The lower bound on $T_{\min}$ aligns with the upper bound in Equation (24) as
$L\geq L_{0}$. The stated value of $n_{\min}$ chooses the smallest value of
$L$ to satisfy both the upper and lower bound on $T_{\min}$.

###### Proof of Theorem 4.10.

  1. 1.

We start with the first case. The idea of the proof is as follows. All the
component $\hat{\boldsymbol{\beta}}_{m}$’s of $\hat{\boldsymbol{\beta}}$ are
approximations of the same underlying ${\boldsymbol{\beta}}_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}={\boldsymbol{\beta}}_{m}$
for all $m\in[M]$. Therefore, by the bound in Proposition 4.9 above and by
convexity of the square function, their mean must be $\ell^{2}$-close to
$\boldsymbol{\beta}_{0}$, with $D=O(ML_{0}^{2}L^{3}\log(ML)\log(2L))$ many
samples. Moreover, the $0.5L^{-1}\epsilon$ perturbation of the coefficients
preserves stability of the $\epsilon$-stable process.

More explicitly, we have

$\displaystyle\Bigg{|}\\!\Bigg{|}\left(\sum_{m=1}^{M}\frac{1}{M}\hat{\boldsymbol{\beta}}_{m}\right)-{\boldsymbol{\beta}}_{0}\Bigg{|}\\!\Bigg{|}_{2}^{2}$
$\displaystyle=\Bigg{|}\\!\Bigg{|}\frac{1}{M}\sum_{m=1}^{M}(\hat{\boldsymbol{\beta}}_{m}-{\boldsymbol{\beta}}_{m})\Bigg{|}\\!\Bigg{|}_{2}^{2}$
$\displaystyle\mbox{convexity}~{}\Rightarrow\hskip 28.45274pt$
$\displaystyle\leq\frac{1}{M}\sum_{m=1}^{M}|\\!|\hat{\boldsymbol{\beta}}_{m}-{\boldsymbol{\beta}}_{m}|\\!|_{2}^{2}$
$\displaystyle=\frac{1}{M}|\\!|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}|\\!|_{2}^{2},$

and by Theorem 4.9, this yields

$\displaystyle\Bigg{|}\\!\Bigg{|}\left(\sum_{m=1}^{M}\frac{1}{M}\hat{\boldsymbol{\beta}}_{m}\right)-{\boldsymbol{\beta}}\Bigg{|}\\!\Bigg{|}_{2}$
$\displaystyle\leq\frac{1}{\sqrt{M}}|\\!|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}|\\!|_{2}$
$\displaystyle\leq\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{MD}}$
(39)

with high probability. Note that

$\alpha=\min_{m\in[M]}\frac{T_{m}\sigma_{m}^{2}}{D}\geq\frac{T_{\min}\sigma_{\min}^{2}}{D}\geq\frac{T_{\min}\sigma_{\min}^{2}}{MT_{\max}}\geq\frac{\sigma_{\min}^{2}}{MC_{\sharp}}\,.$
(40)

When

$D\geq
3^{9}\cdot(84Ae)C_{\sharp}^{5}(1+\epsilon^{-2}+\epsilon^{-4})^{2}\left(\frac{\sigma_{\max}}{\sigma_{\min}}\right)^{4}(\epsilon^{3}\zeta)^{-2}ML_{0}^{2}L^{3}\log\left(\frac{ML}{\delta}\right)\log(2L)\,,$

this yields

$\displaystyle\Bigg{|}\\!\Bigg{|}\left(\sum_{m=1}^{M}\frac{1}{M}\hat{\boldsymbol{\beta}}_{m}\right)-{\boldsymbol{\beta}}\Bigg{|}\\!\Bigg{|}_{2}$
$\displaystyle\leq\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{MD}}$
$\displaystyle\leq\frac{\epsilon\sigma_{\min}^{2}}{\sqrt{3}C_{\sharp}M\alpha}\sqrt{\frac{1}{L\log(2L)}}$
$\displaystyle\eqref{alpha2}~{}\Rightarrow\hskip 28.45274pt$
$\displaystyle<\frac{\epsilon}{\sqrt{3}}\sqrt{\frac{1}{L\log(2L)}}\,.$

Since $\epsilon\in(0,1)$, the stability follows by the triangle inequality.
More explicitly, if $\hat{\boldsymbol{\beta}}_{0}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=M^{-1}(\hat{\boldsymbol{\beta}}_{1}+\cdots+\hat{\boldsymbol{\beta}}_{M})$,
then for every $z\in\mathbb{D}$, we have :

$\displaystyle|\boldsymbol{f}_{\hat{\boldsymbol{\beta}}_{0}}(z)|$
$\displaystyle=|1-\hat{\boldsymbol{\beta}}_{0}\cdot(z,z^{2},\cdots,z^{L})|$
$\displaystyle=|1-{\boldsymbol{\beta}}_{0}\cdot(z,z^{2},\cdots,z^{L})-(\hat{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}_{0})\cdot(z,z^{2},\cdots,z^{L})|$
$\displaystyle\geq|\boldsymbol{f}_{{\boldsymbol{\beta}_{0}}}(z)|-|\\!|\hat{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}_{0}|\\!|_{2}\cdot|\\!|(z,z^{2},\cdots,z^{L})|\\!|_{2}\
$ $\displaystyle>\epsilon-\frac{\epsilon}{\sqrt{3}}\sqrt{\frac{1}{\log(2L)}}$
$\displaystyle>0\,.$

Evidently, this shows that the autoregressive process with
${\hat{\boldsymbol{\beta}}_{0}}$ coefficients is stable.

  2. 2.

The arguments are similar in the second case. Since

$D\geq
3^{9}\cdot(84Ae)C_{\sharp}^{5}(1+\epsilon^{-2}+\epsilon^{-4})^{2}\left(\frac{\sigma_{\max}}{\sigma_{\min}}\right)^{4}(\epsilon^{3}\zeta)^{-2}M^{2}L_{0}^{2}L^{3}\log\left(\frac{ML}{\delta}\right)\log(2L)\,,$

by Theorem 4.9, for any $m\in[M]$ we have

$\displaystyle|\\!|\hat{\boldsymbol{\beta}}_{m}-{\boldsymbol{\beta}}_{m}|\\!|_{2}$
$\displaystyle\leq|\\!|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}|\\!|_{2}$
$\displaystyle\leq\frac{81(84Ae)^{\frac{1}{2}}LL_{0}\sigma_{\max}^{2}C_{\sharp}^{\frac{3}{2}}(1+\epsilon^{-2}+\epsilon^{-4})}{\zeta\alpha\epsilon^{2}}\sqrt{\frac{\log\left(\frac{ML}{\delta}\right)}{D}}$
$\displaystyle\eqref{alpha2}~{}\Rightarrow\hskip 28.45274pt$
$\displaystyle<\frac{\epsilon}{\sqrt{3}}\sqrt{\frac{1}{L\log(2L)}}\,.$

As in the first case, this implies stability for each of the components.

∎

##  5 Algorithmic Aspects

In this section we briefly introduce the main gazette of our algorithm —
namely, the proximal operator — to understand the procedure of solving LASSO
as was done in Nicholson et al. (2020). To sketch the outline of the standard
procedure for solving regularized LASSO penalized with an overlapping group-
norm, we start with the following definition (see Mairal et al. (2011), for
example).

######  Definition 5.1 (Proximal operator).

Given a norm $\mathcal{N}$ on $\mathbb{R}^{ML}$, and a tuning parameter
$\lambda$, the associated proximal operator
$\operatornamewithlimits{Prox}_{\mathcal{N},\lambda}$ is defined for every
$\boldsymbol{\alpha}\in\mathbb{R}^{ML}$ as the optimum value of the following
convex problem:

$\displaystyle\operatornamewithlimits{Prox}_{\mathcal{N},\lambda}(\boldsymbol{\alpha})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\operatornamewithlimits{argmin}_{\boldsymbol{\beta}}\left\\{\frac{1}{2}|\\!|\boldsymbol{\beta}-\boldsymbol{\alpha}|\\!|_{2}^{2}+\lambda\mathcal{N}(\boldsymbol{\beta})\right\\}\,.$

By Zhao et al (2009), the proximal operator
$\operatornamewithlimits{Prox}_{\mathcal{N},\lambda}$ — for $\mathcal{N}$ as
defined in (12) — is the composition

$\operatornamewithlimits{Prox}_{\mathcal{G}_{1},\lambda}\circ\cdots\circ\operatornamewithlimits{Prox}_{\mathcal{G}_{L},\lambda}$

of the proximal operators for the individual groups and can be computed
inductively — starting from
$\operatornamewithlimits{Prox}_{\mathcal{G}_{L},\lambda}$, which is the well-
known soft-thresholding — in $O(ML)$ computational steps. By Combettes and
Wajs (2005, Proposition 3.1(iii)b), the solutions to the regularized LASSO
problem in (15) are precisely the fixed points of the operator

$\boldsymbol{\beta}\longmapsto\operatornamewithlimits{Prox}_{\mathcal{N},\lambda}\left(\boldsymbol{\beta}+\frac{2\lambda}{D}\boldsymbol{X}^{\top}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\right)\,.$
(41)

Here we use an appropriate $\lambda$ (as given in (26) above). As Nicholson et
al. (2020) observed, the proximal operator can be evaluated via duality. The
proximal gradient method of Mairal et al. (2011) then finds the fixed point
(which exists and is unique by convexity of (15)) of this proximal operator.
For pseudo-code of this procedure, see Nicholson et al. (2020), where an
accelerated version of the proximal descent method was employed for achieving
quadratic convergence rate.

##  6 Conclusion

We have established a set-up (see Equation (15) and the discussion preceding
this equation) in which LASSO — regularized with a hierarchical group norm —
can be used to derive statistical guarantees in terms of the one-step ahead
prediction error (Theorem 4.2) in the realms of multiple $\epsilon$-stable
(Definition (2.6)) univariate autoregressive processes of different lengths
but identical true lag $L_{0}$. The results presented here assume no prior
knowledge of the true lag (or any upper-bound of the true lag); in fact, we
show that the sample size itself suggests a certain lag $\hat{L}$ to be used
for the group-LASSO, and given an appropriately large sample size, such that
$\hat{L}$ will be an upper-bound of the true lag $L_{0}$. Moreover, this
$\hat{L}$ will be of the order that is required for our theoretical guarantees
to hold. We proved that the group-LASSO formulated with a suitable tuning
parameter $\lambda$ estimates the AR coefficients with an arbitrarily high
degree of accuracy. We also showed the support of the estimated coefficient-
set approximately matches the support of the original parameters (Theorem
4.9). Finally, we proved that the fitted models with coefficients as estimated
by the group-LASSO are $\epsilon$-stable (Theorem 4.10), a property that is
known in the literature solely for Yule-Walker estimates of the parameters of
univariate autoregressive processes.

From a theoretical perspective, it will be interesting to investigate
adaptations of the group-LASSO method to the case of multiple decoupled AR
processes with multivariate components. We expect that this will require,
among others, (1) additional techniques to deal with the group-norm, and (2)
integrating the stability issues and the restricted eigenvalue issues; these
will be technically far more demanding in the multivariate components
settings. On the practical front, the most important question is to get a
better hold on the tuning parameter $\lambda$ (see Equation (26)), which
requires better constant/parameters than, for example, those appearing in the
proof of Lemma 4.4. A better control on these will enable a more realistic
estimate of $\hat{L}$ — to be used by the group-LASSO as an upper-bound on the
true lag $L_{0}$.

S.C. and J.L. acknowledge funding from the Deutsche Forschungsgemeinschaft
(DFG) under grant number 451920280. S.C. was partially supported by a
‘Research Assistantship Fellowship for Early Postdocs’ from Ruhr-Universität
Bochum Research School by means of the German Academic Exchange Service (DAAD)
STIBET funds.

## References

  * Basu and Michailidis (2015) Sumanta Basu and George Michailidis: Regularized estimation in sparse high-dimensional time series models. _Ann. Statist._ , 43(4): 1535–1567, 2015.
  * Bunea (2008) Florentina Bunea: Honest variable selection in linear and logistic regression models via $\ell_{1}$ and $\ell_{1}+\ell_{2}$ penalization. _Electron. J. Stat. 2_ , 2008.
  * Boyd and Vandenberghe (2004) Stephen Boyd, Lieven Vandenberghe: Convex Optimization. _Cambridge University Press_ , Cambridge, 2004.
  * Candes and Tao (2006) Emmanuel Candes, Terence Tao: Near-optimal signal recovery from random projections: universal encoding strategies? _IEEE Trans. Inform. Theory_ 52, no. 12, 2006
  * Combettes and Wajs (2005) Patrick Combettes, Valérie Wajs: Signal recovery by proximal forward-backward splitting. _Multiscale Model. Simul._ 4, no. 4, 1168–1200, 2005
  * Forni et al. (2000) Mario Forni, Marc Hallin, Marco Lippi, Lucrezia Reichlin: The Generalized Dynamic-Factor Model: Identification and Estimation. _The Review of Economics and Statistics_ , 82, no. 4, 540–554, 2000.
  * Lederer (2021) Johannes Lederer: _Fundamentals of High-dimensional Statistics—with exercises and R labs_. Springer Texts in Statistics. Springer, Cham, 2021.
  * Loh and Wainwright (2012) Po-Ling Loh and Martin J. Wainwright: High-dimensional regression with noisy and missing data: provable guarantees with nonconvexity. _Ann. Statist._ , 40(3):1637–1664, 2012.
  * Lütkepohl (2005) Helmut Lütkepohl: _New Introduction to Multiple Time Series Analysis_. Springer-Verlag, Berlin, 2005.
  * Mairal et al. (2011) Julien Mairal, Rodolphe Jenatton, Guillaume Obozinski, and Francis Bach: Convex and network flow optimization for structured sparsity. _J. Mach. Learn. Res._ , 12:2681–2720, 2011.
  * Masini et-al (2022) Ricardo Masini, Marcelo Medeiros, and Eduardo Mendes: Regularized estimation of high-dimensional vector autoregressions with weakly dependent innovations. _J. Time Ser. Anal._ , 43: 532–557, 2011.
  * Nardi and Rinaldo (2011) Yuval Nardi and Alessandro Rinald: Autoregressive process modeling via the Lasso procedure. _J. Multivariate Anal._ 102, no. 3, 528–549, 2011
  * Negahban and Wainwright (2011) Sahand Negahban and Martin J. Wainwright: Estimation of (near) low-rank matrices with noise and high-dimensional scaling. _Ann. Statist._ , 39(2):1069–1097, 2011.
  * Nicholson et al. (2020) William B. Nicholson, Ines Wilms, Jacob Bien, and David S. Matteson: High dimensional forecasting via interpretable vector autoregression. _J. Mach. Learn. Res._ , 21:Paper No. 166, 52p, 2020.
  * Tseng (2009) Paul Tseng: On accelerated proximal gradient methods for convex-concave optimization. 2008. available at _http://www. math. washington. edu/~ tseng/papers/apgm. pdf_.
  * Vershynin (2018) Roman Vershynin: _High-dimensional probability_. An Introduction with Applications in Data Science. With a foreword by Sara van de Geer. Cambridge Series in Statistical and Probabilistic Mathematics, 47. Cambridge University Press, Cambridge, 2018.
  * Wong et al. (2020) Kam Chung Wong, Zifan Li, and Ambuj Tewari: Lasso guarantees for $\beta$-mixing heavy-tailed time series. _Ann. Statist._ , 48(2), 2020.
  * Zhao et al (2009) Peng Zhao, Guilherme Rocha, Bin Yu: The composite absolute penalties family for grouped and hierarchical variable selection. _Ann. Statist._ 37, no. 6A, 3468–3497, 2009

##  Appendix A Proof of Theorem 4.2

###### Proof.

Denote the data by ${\mathcal{X}}$. Define the random vector
$\boldsymbol{U}\in\mathbb{R}^{D}$ as follows:

$\boldsymbol{U}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(U_{1}^{1},\dots,U_{T_{1}}^{1},U_{1}^{2},\dots,U_{T_{2}}^{2},\cdots,U_{1}^{M},\dots,U_{T_{M}}^{M})^{\top}\,.$

We note the following series of inequalities, all of which follow from
linearity of expectation: here, we write $\boldsymbol{\beta}$ for the true
coefficient vector of the model. One has

$\displaystyle\mathbb{E}\left[|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid
X\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[|\\!|\boldsymbol{y}-X{\boldsymbol{\beta}}|\\!|_{2}^{2}\mid
X\right]+\mathbb{E}\left[|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\mid
X\right]$ $\displaystyle\hskip
85.35826pt+2\mathbb{E}\left[\langle\boldsymbol{y}-X\boldsymbol{\beta},X\boldsymbol{\beta}-X\hat{\boldsymbol{\beta}}\rangle\mid
X\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[|\\!|\boldsymbol{U}|\\!|_{2}^{2}\mid
X\right]+\mathbb{E}\left[|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\mid
X\right]$ $\displaystyle\hskip
85.35826pt+2\mathbb{E}\left[\langle\boldsymbol{U},X\boldsymbol{\beta}-X\hat{\boldsymbol{\beta}}\rangle\mid
X\right]$ $\displaystyle=$ $\displaystyle
T_{1}\sigma_{1}^{2}+\cdots+T_{M}\sigma_{M}^{2}+\mathbb{E}\left[|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\mid
X\right]+2\langle\mathbb{E}\left[\boldsymbol{U}\mid
X\right],X\boldsymbol{\beta}-X\hat{\boldsymbol{\beta}}\rangle$
$\displaystyle=$ $\displaystyle
T_{1}\sigma_{1}^{2}+\cdots+T_{M}\sigma_{M}^{2}+|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\,,$

where the last equality follows from $\mathbb{E}\left[\boldsymbol{U}\mid
X\right]=\boldsymbol{U}$, and
$\boldsymbol{U}\perp(X\boldsymbol{\beta}-X\hat{\boldsymbol{\beta}})$.
Therefore, we need to derive an estimate of the error
$|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$. We notice
that

$\displaystyle|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{M}\sum_{t=1}^{T_{m}}\left(\sum_{l=1}^{L}\left(\hat{\beta}_{l}^{m}-\boldsymbol{\beta}_{l}^{m}\right)x_{t-l}^{m}\right)^{2}\,,$

where $\\{\boldsymbol{\beta}_{l}^{m}\\}$ denotes the true-parameters of the
models. Using the fact that $\hat{\boldsymbol{\beta}}$ is a minimizer of
$\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\alpha}|\\!|_{2}^{2}+\lambda\sum_{l=1}^{L}\mathcal{N}({\boldsymbol{\alpha}})$,
over $\boldsymbol{\alpha}\in\mathbb{R}^{ML}$, one derives

$\displaystyle\frac{1}{D}|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}+\lambda\mathcal{N}(\hat{\boldsymbol{\beta}})$
$\displaystyle\leq\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\alpha}|\\!|_{2}^{2}+\lambda\mathcal{N}({\boldsymbol{\alpha}})$
$\displaystyle\Rightarrow$
$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle\leq\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+\frac{2}{D}\langle\boldsymbol{U},X(\hat{\boldsymbol{\beta}}-\boldsymbol{\alpha})\rangle+\lambda\mathcal{N}(\boldsymbol{\alpha})-\lambda\mathcal{N}(\hat{\boldsymbol{\beta}})$
$\displaystyle\leq\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+\frac{2}{D}\langle
X^{\top}\boldsymbol{U},\hat{\boldsymbol{\beta}}-\boldsymbol{\alpha}\rangle+\lambda\mathcal{N}(\boldsymbol{\alpha})-\lambda\mathcal{N}(\hat{\boldsymbol{\beta}})$

for all $\boldsymbol{\alpha}\in\mathbb{R}^{ML}$. By the Cauchy-Schwarz
inequality, this implies

$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+\frac{2}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})\mathcal{N}(\hat{\boldsymbol{\beta}}-\boldsymbol{\alpha})+\lambda(\mathcal{N}({\boldsymbol{\alpha}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))$
$\displaystyle\leq$
$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\boldsymbol{\alpha})|\\!|_{2}^{2}+\lambda\mathcal{N}(\hat{\boldsymbol{\beta}}-{\boldsymbol{\alpha}})+\lambda(\mathcal{N}({\boldsymbol{\alpha}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))$

since $\lambda\geq\frac{2}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})$.
Finally, the triangle inequality
$\mathcal{N}(\hat{\boldsymbol{\beta}}-{\boldsymbol{\alpha}})\leq\mathcal{N}(\hat{\boldsymbol{\beta}})+\mathcal{N}({\boldsymbol{\alpha}})$
establishes the first part of the theorem.

In the following, we will choose a slightly different threshold for the tuning
parameter $\lambda$, as proposed in the statement of theorem. Let
$\boldsymbol{v}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\boldsymbol{\beta}+\hat{\boldsymbol{\beta}}$,
and notice that the true lag $L_{0}$ satisfies

$L_{0}=\max\\{\ell\in[L]:|\\!|\boldsymbol{\beta}|\\!|_{{\mathcal{G}}_{\ell}}\neq
0\\}\,.$

We have

$\displaystyle\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\boldsymbol{\beta}+\boldsymbol{v})$
$\displaystyle=\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}|\\!|_{{\mathcal{G}}_{\ell}}+\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\cancelto{0}{\boldsymbol{\beta}}|\\!|_{{\mathcal{G}}_{\ell}}$
$\displaystyle\hskip
28.45274pt-\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}+\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\cancelto{0}{\boldsymbol{\beta}}+\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}$
$\displaystyle=\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}+\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}$
$\displaystyle\leq\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}},$
(42)

where the last step is due to the triangle inequality:
$|\\!|\boldsymbol{\beta}+\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}\geq|\\!|\boldsymbol{\beta}|\\!|_{{\mathcal{G}}_{\ell}}-|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}$.
Suppose that
$\lambda\geq\frac{4}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})$. Using the
fact that $\hat{\boldsymbol{\beta}}$ is a minimizer of the objective
$\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\alpha}|\\!|_{2}^{2}+\lambda\mathcal{N}({\boldsymbol{\alpha}})$,
over $\boldsymbol{\alpha}\in\mathbb{R}^{ML}$, one derives

$\displaystyle\frac{1}{D}|\\!|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}|\\!|_{2}^{2}+\lambda\mathcal{N}(\hat{\boldsymbol{\beta}})$
$\displaystyle\leq\frac{1}{D}|\\!|\boldsymbol{y}-X\boldsymbol{\beta}|\\!|_{2}^{2}+\lambda\mathcal{N}({\boldsymbol{\beta}})$
$\displaystyle\Rightarrow$
$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle\leq\frac{2}{D}\langle
X^{\top}\boldsymbol{U},\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\rangle+\lambda\mathcal{N}(\boldsymbol{\beta})-\lambda\mathcal{N}(\hat{\boldsymbol{\beta}})\,.$

By the Cauchy-Schwarz inequality, this implies

$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle\leq\frac{2}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})\mathcal{N}(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})+\lambda(\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))$
$\displaystyle\leq\frac{\lambda}{2}\mathcal{N}(\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}})+\lambda(\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))\hskip
28.45274pt\because~{}\lambda\geq\frac{4}{D}\mathcal{N}_{\star}(X^{\top}\boldsymbol{U})$
$\displaystyle\leq\frac{\lambda}{2}\min\left\\{3\mathcal{N}(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}),3\mathcal{N}(\boldsymbol{\beta})-\mathcal{N}(\hat{\boldsymbol{\beta}})\right\\}\,.$
(43)

The inequalities

$0\leq\frac{2}{\lambda
D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}\leq\mathcal{N}(\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}})+2(\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))\,,$

and

$\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}})\leq\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}$

yield the following:

$\displaystyle 0$
$\displaystyle\leq\mathcal{N}(\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}})+2(\mathcal{N}({\boldsymbol{\beta}})-\mathcal{N}(\hat{\boldsymbol{\beta}}))$
$\displaystyle\leq\mathcal{N}(\hat{\boldsymbol{v}})+2\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-2\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}$
$\displaystyle=3\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}-\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{v}|\\!|_{{\mathcal{G}}_{\ell}}\,.$

That is, we have

$\displaystyle\sum_{\ell=L_{0}+1}^{L}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}-\boldsymbol{\hat{\beta}}|\\!|_{{\mathcal{G}}_{\ell}}$
$\displaystyle\leq
3\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}-\boldsymbol{\hat{\beta}}|\\!|_{{\mathcal{G}}_{\ell}}\,.$
(44)

By equation (A), we then have

$\displaystyle\frac{1}{D}|\\!|X(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})|\\!|_{2}^{2}$
$\displaystyle\leq
2\lambda\sum_{\ell=1}^{L_{0}}\sqrt{|{\mathcal{G}}_{\ell}|}\cdot|\\!|\boldsymbol{\beta}-\boldsymbol{\hat{\beta}}|\\!|_{{\mathcal{G}}_{\ell}}$
(45) $\displaystyle=2\lambda{\mathcal{N}}_{\leq
L_{0}}(\boldsymbol{\beta}-\boldsymbol{\hat{\beta}})\,.$

∎

######  Remark (A).

If, instead of $\mathcal{N}$, we take standard $\ell^{2}$ in $\mathbb{R}^{ML}$
(which correspond to the largest group in $\mathcal{N}$), the argument leading
to (44) yields

$\displaystyle\sum_{\ell=L_{0}+1}^{L}|\boldsymbol{\hat{\beta}}^{m}_{\ell}|$
$\displaystyle=\sum_{\ell=L_{0}+1}^{L}|\boldsymbol{\beta}^{m}_{\ell}-\boldsymbol{\hat{\beta}}^{m}_{\ell}|$
$\displaystyle\leq
3\sum_{\ell=1}^{L_{0}}|\boldsymbol{\beta}^{m}_{\ell}-\boldsymbol{\hat{\beta}}^{m}_{\ell}|\,.$

##  Appendix B Proof of Theorem 4.3

###### Proof.

Since Lemma D.3 implies

$\displaystyle{\mathcal{N}}_{\star}\left(\frac{2}{D}X^{\top}\boldsymbol{U}\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{L}}\left\|\frac{2}{D}X^{\top}\boldsymbol{U}\right\|_{\infty}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{L}}\max_{m\in[M]}\max_{l\in[L]}\left|\frac{2}{D}\sum_{t=1}^{T_{m}}U^{m}_{t}X^{m}_{t-l}\right|\,,$

it suffices, by an union bound argument, to find high-probability upper bound
of the inner-product (see (16) for notations)

$\displaystyle\tau_{m,l}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=2\langle\boldsymbol{U}^{m},X^{m,l}\rangle$
$\displaystyle=2\left(\sum_{t=1}^{T_{m}}U^{m}_{t}X^{m}_{t-l}\right)\,.$ (46)

We now present standard techniques, as in Basu and Michailidis (2015), for
example. One has

$\displaystyle\tau_{m,l}$ $\displaystyle=$
$\displaystyle\left\|\boldsymbol{U}^{(m)}+X^{(m,l)}\right\|^{2}_{2}-\left\|\boldsymbol{U}^{(m)}\right\|_{2}^{2}-\left\|X^{(m,l)}\right\|^{2}_{2}\,.$

Hence, for any $\eta>0$, an union bound argument, together with mutual
independence of the Gaussian random variable $U^{m}_{j}$ and the variables
$X^{m}_{<j}$ for each $j\in[T_{m}]$, implies

$\displaystyle\mathbb{P}\left[\left|\frac{\tau_{m,l}}{D}\right|>\frac{3}{2}\eta\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left[\left|(\boldsymbol{U}^{(m)})^{\top}(\boldsymbol{U}^{(m)})-\text{Var}(U^{m})\right|>{\frac{D}{2}\eta}\right]+\mathbb{P}\left[\left|(X^{(m,l)})^{\top}(X^{(m,l)})-\text{Var}(X^{m})\right|>{\frac{D}{2}\eta}\right]$
$\displaystyle\hskip
81.09052pt+~{}\mathbb{P}\left[\left|(\boldsymbol{U}^{(m)}+X^{(m,l)})^{\top}(\boldsymbol{U}^{(m)}+X^{(m,l)})-\text{Var}(X^{m}+U^{m})\right|>{\frac{D}{2}\eta}\right]\,.$

Here, $\text{Var}(Y)$ denotes the trace of the covariance matrix of the random
variable $Y$. We now estimate each summand separately, starting with

$p_{m,l,\eta}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{P}\left[\left|(\boldsymbol{U}^{(m)})^{\top}(\boldsymbol{U}^{(m)})-\text{Var}(U^{(m)})\right|>\frac{D}{2}\eta\right]\,.$

By the running assumptions, the vector $\boldsymbol{U}^{(m)}$ is
$T_{m}$-dimensional mean zero Gaussian having covariance matrix
$\sigma_{m}^{2}\mathbb{I}_{T_{m}}$; by the inequality in proposition E.1, one
has

$\displaystyle p_{m,l,\eta}$ $\displaystyle\leq
2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}}\min\\{1,\frac{D\eta}{8\sigma_{m}^{2}T_{m}}\\}}\,.$
(47)

Next, we consider

$q_{m,l,\eta}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{P}\left[\left|(\boldsymbol{X}^{(m,l)})^{\top}(\boldsymbol{X}^{(m,l)})-\text{Var}(X^{(m,l)})\right|>\frac{D}{2}\eta\right]\,.$

Recall that the random vector $X^{(m,l)}$ is mean-zero Gaussian with
covariance matrix $\Gamma^{(m)}$ mentioned in equation (59) below. One has

$\text{Var}(X^{(m,l)})=\text{tr}(\Gamma^{(m)})=T_{m}\mathbb{E}[(X_{0}^{m})^{2}]\,.$

From lemma E.2, one has
$|\\!|\Gamma^{(m)}|\\!|_{\text{op}}\leq\epsilon^{-2}\sigma_{m}^{2}$; thus,
proposition E.1 yields

$\displaystyle q_{m,l,\eta}$ $\displaystyle\leq
2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}\epsilon^{-2}}\min\\{1,\frac{D\eta}{8\sigma_{m}^{2}T_{m}\epsilon^{-2}}\\}}\,.$
(48)

Finally, we consider

$r_{m,l,\eta}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{P}\left[\left|(\boldsymbol{U}^{(m)}+X^{(m,l)})^{\top}(\boldsymbol{U}^{(m)}+X^{(m,l)})-\text{Var}(X^{m}+U^{m})\right|>\frac{D}{2}\eta\right]\,.$

Note that $\boldsymbol{U}^{(m)}+X^{(m,l)}$ is a mean zero Gaussian with
symmetric covariance matrix $\tilde{\Gamma}^{(m)}$, whose $(t,s)$-entry (for
$t\geq s$) is given by

$\displaystyle\tilde{\Gamma}^{(m)}_{t,s}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{E}[(X^{m}_{t-l}+U^{m}_{t})(X^{m}_{s-l}+U^{m}_{s})]$
$\displaystyle=\mathbb{E}[(X^{m}_{t-l}X^{m}_{s-l}]+\mathbb{E}[X^{m}_{s-l}U^{m}_{t}]+\mathbb{E}[X^{m}_{t-l}U^{m}_{s}]+\mathbb{E}[U^{m}_{t}U^{m}_{s}]$
$\displaystyle=\mathbb{E}[(X^{m}_{t-l}X^{m}_{s-l}]+\mathbb{E}[X^{m}_{t-l}U^{m}_{s}]+\sigma_{m}^{2}\boldsymbol{1}_{t=s}\,.$

Together with lemma E.3, this implies

$\displaystyle|\\!|\tilde{\Gamma}^{(m)}|\\!|_{\text{op}}$
$\displaystyle\leq\epsilon^{-2}\sigma_{m}^{2}+\epsilon^{-4}\sigma_{m}^{2}+\sigma_{m}^{2}$
$\displaystyle=\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})\,,$

and by proposition E.1, we fetch

$\displaystyle r_{m,l,\eta}$ $\displaystyle\leq
2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}\min\\{1,\frac{D\eta}{8\sigma_{m}^{2}T_{m}(1+\epsilon^{-2}+\epsilon^{-4})}\\}}\,.$
(49)

Now suppose that

$\eta\geq\frac{8C_{\sharp}\sigma^{2}_{\max}(1+\epsilon^{-2}+\epsilon^{-4})}{M}\,;$

then, the inequalities in (47), (48), and (49) are simplified as follows:

$\displaystyle p_{m,l,\eta}$ $\displaystyle\leq
2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}}}\,;$ (50) $\displaystyle q_{m,l,\eta}$
$\displaystyle\leq 2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}\epsilon^{-2}}}\;,$
$\displaystyle r_{m,l,\eta}$ $\displaystyle\leq
2e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}}\,.$

Of these, the right hand side is the largest in the bottom-most inequality
(50). Note that, the union bound implies

$\displaystyle\mathbb{P}\left[\sup_{m,l}\left|\frac{\tau_{m,l}}{D}\right|>\frac{3}{2}\eta\right]$
$\displaystyle\leq$
$\displaystyle\displaystyle\sum_{m,l}\mathbb{P}\left[\left|\frac{\tau_{m,l}}{D}\right|>\frac{3}{2}\eta\right]$
$\displaystyle\leq$
$\displaystyle\displaystyle\sum_{l\in[L]}\displaystyle\sum_{m\in[M]}\mathbb{P}\left[\left|\frac{\tau_{m,l}}{D}\right|>\frac{3}{2}\eta\right]$
$\displaystyle\eqref{eq:4.66}~{}\Rightarrow\hskip 28.45274pt$
$\displaystyle\leq$ $\displaystyle
6\displaystyle\sum_{l\in[L]}\displaystyle\sum_{m\in[M]}e^{-\frac{c_{0}D\eta}{8\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}}\,.$

Thus, for any $\delta>0$, in order to have the inequality

$\displaystyle\mathbb{P}\left[{\mathcal{N}}_{\star}\left(\frac{2}{D}X^{\top}\boldsymbol{U}\right)>\frac{3\eta}{2\sqrt{L}}\right]$
$\displaystyle\leq\delta\,,$

it suffices to have

$\exp\left(-\frac{c_{0}D\eta}{8\sigma_{m}(1+\epsilon^{-2}+\epsilon^{-4})}\right)\leq\frac{\delta}{ML}\,$
(51)

for each $m\in[M]$. Equivalently, it suffices to have the following for each
$m\in[M]$:

$\displaystyle{\frac{c_{0}D\eta}{8\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}}$
$\displaystyle\geq\log\left(\frac{ML}{\delta}\right)$
$\displaystyle\mbox{equivalently,}\hskip 28.45274ptD$
$\displaystyle\geq\frac{8\sigma_{m}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,.$

This is equivalent to

$\displaystyle D$
$\displaystyle\geq\frac{8\sigma_{\max}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,.$
(52)

It then follows from the argument above that

$\mathbb{P}\left[\frac{2}{D}{\mathcal{N}}_{\star}(X^{\top}\boldsymbol{U})\geq\frac{3\eta}{2\sqrt{L}}\right]~{}\leq~{}\delta$
(53)

holds, provided

$\displaystyle D$
$\displaystyle\geq\frac{8\sigma_{\max}^{2}(1+\epsilon^{-2}+\epsilon^{-4})}{c_{0}\eta}\log\left(\frac{ML}{\delta}\right)\,.$

∎

##  Appendix C Proof of Proposition 4.6

We first observe what happens in the $M=1$ case: this amounts to restricting
$X$ to the top left block

$X_{1}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\begin{pmatrix}x_{1-1}^{1}&\dots&x_{1-L}^{1}\\\
\cdot&\dots&\cdot\\\ \cdot&\dots&\cdot\\\
x_{T_{1}-1}^{1}&\dots&x_{T_{1}-L}^{1}\end{pmatrix}\,.$

######  Lemma C.1 (Blockwise concentration inequality).

For any integer $s_{1}>0$, write $\mathbb{K}(s_{1})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{B}_{0}(s_{1})\cap\mathbb{B}_{1}(1)$.
Then

$\displaystyle\mathbb{P}\left[\sup_{\boldsymbol{v}_{1}\in\mathbb{K}(s_{1})}\left|\boldsymbol{v}_{1}^{\top}(X_{1}^{\top}X_{1}-T_{1}\Gamma^{(1)})\boldsymbol{v}_{1}\right|\geq\frac{T_{1}^{2}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{2\mathfrak{m}(f)}\right]$
$\displaystyle\leq
2\exp\left(-\frac{T_{1}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}{2}+s_{1}\min\\{\log
L,~{}\log(21eL/s_{1})\\}\right)\,.$ (54)

###### Proof of Lemma C.1.

Let $\boldsymbol{v}_{1}\in\mathbb{R}^{L}$ be a fixed unit-normed vector. Then,
$\boldsymbol{u}_{1}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=X_{1}\boldsymbol{v}_{1}$
is a random mean-zero Gaussian vector in $\mathbb{R}^{T_{1}}$, with covariance
matrix $\Sigma_{\boldsymbol{v}_{1}}^{X_{1}}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{E}\left[X_{1}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}X_{1}^{\top}\right]$.
One has

$\displaystyle(\Sigma_{\boldsymbol{v}_{1}}^{X_{1}})_{r,s}$
$\displaystyle=\sum_{j,l=1}^{L}\mathbb{E}\left[X^{1}_{r-j}\boldsymbol{v}_{1,j}\boldsymbol{v}_{1,l}X^{1}_{s-l}\right]$
$\displaystyle=\sum_{j,l=1}^{L}\boldsymbol{v}_{1,j}\mathbb{E}\left[X^{1}_{r-j}X^{1}_{s-l}\right]\boldsymbol{v}_{1,l}$
$\displaystyle=\boldsymbol{v}_{1}^{\top}\Gamma_{r,s}\boldsymbol{v}_{1}\,,$

where $\Gamma_{r,s}$ is the covariance of the vectors $X^{(1)}_{r}$ and
$X^{(1)}_{s}$. Note that

$\displaystyle\text{tr}(\Sigma_{\boldsymbol{v}_{1}}^{X_{1}})$
$\displaystyle=\mathbb{E}\left[\text{tr}(X_{1}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}X_{1}^{\top})\right]$
$\displaystyle=T_{1}\mathbb{E}\left[(X_{1}\boldsymbol{v}_{1})_{11}^{2}\right]$
$\displaystyle=T_{1}\boldsymbol{v}_{1}^{\top}\Gamma^{(1)}\boldsymbol{v}_{1}\,,$

where $\Gamma^{(1)}$ is the autocovariance matrix of $\\{X^{1}_{t}\\}$.
Moreover, one has by Lemma E.2, the inequality

$\displaystyle|\\!|\Sigma_{\boldsymbol{v}_{1}}^{X_{1}}|\\!|_{\text{op}}\leq\frac{\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\,.$

By the inequality in Proposition E.1, one has

$\displaystyle\mathbb{P}\left[\left|\boldsymbol{u}_{1}^{T}\boldsymbol{u}_{1}-\text{tr}(\Sigma_{\boldsymbol{v}_{1}}^{X_{1}})\right|\geq
4T_{1}\zeta_{1}|\\!|\Sigma_{\boldsymbol{v}_{1}}^{X_{1}}|\\!|_{\text{op}}\right]$
$\displaystyle\leq 2e^{-\frac{T_{1}}{2}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}$

for any $\zeta_{1}>0$. In particular, this implies that for any fixed
$\boldsymbol{v}_{1}\in\mathbb{R}^{L}$, the inequality

$\displaystyle\mathbb{P}\left[\left|\boldsymbol{v}_{1}^{\top}(X_{1}^{\top}X_{1}-T_{1}\Gamma^{(1)})\boldsymbol{v}_{1}\right|\geq\frac{4T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\right]$
$\displaystyle\leq 2e^{-\frac{T_{1}}{2}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}$
(55)

holds for any $\zeta_{1}>0$. Applying Lemma E.4 with $G\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=X_{1}^{\top}X_{1}-T_{1}\Gamma^{(1)}$,
we conclude that

$\displaystyle\mathbb{P}\left[\sup_{\boldsymbol{v}_{1}\in\mathbb{K}(s_{1})}\left|\boldsymbol{v}_{1}^{\top}(X_{1}^{\top}X_{1}-T_{1}\Gamma^{(1)})\boldsymbol{v}_{1}\right|\geq\frac{4T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\right]$
$\displaystyle\leq$ $\displaystyle
2\exp\left(-\frac{T_{1}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}{2}+s_{1}\min\\{\log
L,~{}\log(21eL/s_{1})\\}\right)$ (56)

holds for any integer $s_{1}\geq 1$, and any $\zeta_{1}>0$. ∎

### Proof of Proposition 4.6

###### Proof.

When $s_{1}>0$ is even, it follows from Lemma E.5 that

$\displaystyle\mathbb{P}\left[\sup_{\boldsymbol{v}_{1}\in\mathbb{R}^{L}}\left|\boldsymbol{v}_{1}(X^{\top}_{1}X_{1}-T_{1}\Gamma^{(1)})\boldsymbol{v}_{1}\right|\leq\frac{108T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\left(|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}+\frac{2}{s_{1}}|\\!|\boldsymbol{v}_{1}|\\!|_{1}^{2}\right)\right]$
$\displaystyle\geq~{}$ $\displaystyle
1-2\exp\left(-\frac{T_{1}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}{2}+s_{1}\min\\{\log
L,~{}\log(21eL/s_{0})\\}\right)\,.$

An application of Lemma E.2 implies that for

$\delta\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=2e^{-\frac{T_{1}\min\\{\zeta_{1},\zeta_{1}^{2}\\}}{2}+s_{1}\min\\{\log
L,~{}\log(21eL/s_{1})\\}}\,,$

the conclusion of the following sequence of inequalities holds for all
$\boldsymbol{v}_{1}\in\mathbb{R}^{L}$ with probability at least $1-\delta$:

$\displaystyle|\\!|X_{1}\boldsymbol{v}_{1}|\\!|_{2}^{2}$ $\displaystyle\geq
T_{1}\lambda_{\min}(\Gamma^{(1)})|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}-\frac{108T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\left(|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}+\frac{2}{s_{1}}|\\!|\boldsymbol{v}_{1}|\\!|_{1}^{2}\right)$
$\displaystyle\geq\frac{T_{1}\Lambda_{\min}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{M}(f)}|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}-\frac{108T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f)}\left(|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}+\frac{2}{s_{1}}|\\!|\boldsymbol{v}_{1}|\\!|_{1}^{2}\right)$
$\displaystyle=T_{1}\left(\frac{\Lambda_{\min}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{M}(f_{1})}-\frac{108\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{\mathfrak{m}(f_{1})}\right)|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}-\frac{216T_{1}\zeta_{1}\Lambda_{\max}(\Sigma_{\boldsymbol{U}_{1}})}{s_{1}\mathfrak{m}(f_{1})}|\\!|\boldsymbol{v}_{1}|\\!|_{1}^{2}\,.$

We specialize to the case of univariate $\epsilon$-stable autoregressive time
series, with $\Sigma_{\boldsymbol{U}_{1}}=\sigma_{1}^{2}$, and (by
$\epsilon$-stability)

$\epsilon^{2}\leq\mathfrak{m}(f_{1})\leq\mathfrak{M}(f_{1})\leq\epsilon^{-2}.$

Letting $\zeta_{1}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=6^{-3}\epsilon^{4}$,
we obtain

$\displaystyle\mathbb{P}\left[\inf_{\boldsymbol{v}_{1}\in\mathbb{R}^{L}}\left(|\\!|X_{1}\boldsymbol{v}_{1}|\\!|_{2}^{2}-\frac{T_{1}}{2}\sigma_{1}^{2}\epsilon^{2}\left(|\\!|\boldsymbol{v}_{1}|\\!|_{2}^{2}-\frac{2}{s_{1}}|\\!|\boldsymbol{v}_{1}|\\!|_{1}^{2}\right)\right)\geq
0\right]$ (57) $\displaystyle\geq$ $\displaystyle
1-2e^{-\frac{T_{1}}{2}\zeta_{1}^{2}+s_{0}\min\\{\log
L,~{}\log(21eL/s_{0})\\}}\,.$

Finally, we consider all the blocks of $X^{\top}X$; the deviation inequality
above yields

$\displaystyle\mathbb{P}\left[\displaystyle\forall_{m\in[M]}\inf_{\boldsymbol{v}_{m}\in\mathbb{R}^{L}}\left(|\\!|X_{m}\boldsymbol{v}_{m}|\\!|_{2}^{2}-\frac{T_{m}\sigma^{2}\epsilon^{2}}{2}\left(|\\!|\boldsymbol{v}_{m}|\\!|_{2}^{2}-\frac{2}{s_{m}}|\\!|\boldsymbol{v}_{m}|\\!|_{1}^{2}\right)\right)\geq
0\right]$ (58) $\displaystyle\geq$ $\displaystyle
1-2\sum_{m=1}^{M}e^{-\frac{T_{m}\zeta^{2}}{2}+s_{m}\min\\{\log
L,~{}\log(21eL/s_{m})\\}}\,.$

This proves the proposition. ∎

##  Appendix D More on Dual Norms

Since ${\mathcal{N}}_{\star}(\boldsymbol{0})=0$, we will assume, in the
following computation of the dual norm, that
$\boldsymbol{\alpha}\neq\boldsymbol{0}$. This ensures
${\mathcal{N}}_{\star}(\boldsymbol{\alpha})>0$ too. Note that
$\mathcal{N}_{\star}(\boldsymbol{\alpha})=\mathcal{N}_{\star}(\boldsymbol{|\alpha|})$,
where $|\boldsymbol{\alpha}|=(|\alpha_{1}|,\cdots,|\alpha_{L}|)$; this is
because ${\mathcal{N}}$ is invariant under arbitrary sign changes of the
coordinates of its argument. Thus, for
$\mathcal{N}_{\star}(\boldsymbol{\alpha})=\boldsymbol{\alpha}\cdot\boldsymbol{\beta}$,
we may assume (without loss of generality) that
$|\boldsymbol{\alpha}|=\boldsymbol{\alpha}$ and
$|\boldsymbol{\beta}|=\boldsymbol{\beta}$.

The advantage of the hierarchical group norm
${\mathcal{N}}(\boldsymbol{\beta})$ — in the setting of sparse recovery via
regularized least square regression — is that, while it sets any group of
parameters to zero — because of the hierarchy in the group structure — all
variables in all groups that appear further down the order of hierarchy are
set to zero automatically. Thus, these norms are well-suited for determination
of the true lag order in the context of autoregressive process.

In order to use this group norm in our analysis, we need to collect some basic
facts on the norm. Note that the norm resembles $\ell^{1}$ at the group level,
while within each group it is the $\ell^{2}$ norm; from this, one might expect
that the dual of the norm should “resemble” the $\ell^{\infty}$-norm at the
group level, and within a group it should be $\ell^{2}$. We will prove that
this intuition goes quite well, in the sense that the actual dual norm can be
upper-bounded by this mixed $\ell^{\infty,2}$ norm.

######  Lemma D.1 (Dual norm is attained on the boundary).

$\displaystyle{\mathcal{N}}_{\star}(\boldsymbol{\alpha})$ $\displaystyle=$
$\displaystyle\max_{{\mathcal{N}}(\boldsymbol{\beta})=1}\langle\boldsymbol{\alpha},\boldsymbol{\beta}\rangle\,.$

###### Proof of Lemma D.1.

Note that — by continuity of
$\boldsymbol{\beta}\mapsto\langle\boldsymbol{\alpha},\boldsymbol{\beta}\rangle$,
compactness of the subset $C_{\boldsymbol{\beta}}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{\boldsymbol{\beta}:{\mathcal{N}}(\boldsymbol{\beta})\leq
1\\}$, and since ${\mathcal{N}}_{\star}(\boldsymbol{\alpha})\neq 0$ — there is
nonzero $\boldsymbol{\beta}_{0}\in C_{\boldsymbol{\beta}}$ such that
${\mathcal{N}}_{\star}(\boldsymbol{\alpha})=\langle\boldsymbol{\alpha},\boldsymbol{\beta}_{0}\rangle$.
Suppose, if possible, that $\boldsymbol{\beta}_{0}$ satisfies
${\mathcal{N}}(\boldsymbol{\beta}_{0})<1$. In $(0,+\infty)$, the function
$c\mapsto{\mathcal{N}}(c\boldsymbol{\beta}_{0})$ is strictly increasing
(continuous) function; thus, there is $c>1$ such that
$c{\mathcal{N}}(\boldsymbol{\beta}_{0})={\mathcal{N}}(\boldsymbol{c\beta}_{0})\leq
1$; one has

$\displaystyle\max\\{\langle\boldsymbol{\alpha},\boldsymbol{\beta}\rangle:{{\mathcal{N}}(\boldsymbol{\beta})\leq
1}\\}$
$\displaystyle\geq\langle\boldsymbol{\alpha},c\boldsymbol{\beta_{0}}\rangle$
$\displaystyle=c\cdot\langle\boldsymbol{\alpha},\boldsymbol{\beta}_{0}\rangle$
$\displaystyle>\langle\boldsymbol{\alpha},\boldsymbol{\beta}_{0}\rangle$
$\displaystyle={\mathcal{N}}_{\star}(\boldsymbol{\alpha})\,,$

which is a contradiction. ∎

In order to facilitate our computations of the dual norm in $\mathbb{R}^{ML}$,
we start with the special case of $M=1$; after Proposition D.2, we will extend
this to the general case. Now, we define

${\mathcal{N}}^{1}(\boldsymbol{\alpha})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sum_{l=1}^{L}\sqrt{(L-l+1)\sum_{j=l}^{L}\alpha_{j}^{2}}\,.$

Let us write $\boldsymbol{\alpha}_{\geq
j}=(\boldsymbol{0},\alpha_{j},\alpha_{j+1},\cdots,\alpha_{L})$ for $j\in[L]$.
From now on, we assume that

$\prod_{j\in[L]}\alpha_{j}\neq 0\,.$

We have the following proposition.

######  Proposition D.2 ($L^{\infty}$ norm bounds dual norm, $M=1$ case).

The following inequality holds for $\boldsymbol{\alpha}\in\mathbb{R}^{L}$:

$\displaystyle{\mathcal{N}}^{1}_{\star}(\boldsymbol{\alpha})$
$\displaystyle\leq L^{-\frac{1}{2}}|\\!|\boldsymbol{\alpha}|\\!|_{\infty}\,.$

###### Proof of Proposition D.2.

We apply the elementary inequalities

$\displaystyle(\sqrt{L-l+1}|\\!|\beta_{\geq
l}|\\!|_{2}-\sqrt{L-l+1}|\beta_{l}|)^{2}$
$\displaystyle\geq(L-l+1)|\\!|\beta_{\geq
l}|\\!|_{2}^{2}-(L-l+1)|\beta_{l}|^{2}$
$\displaystyle=(L-l+1)\sum_{j>l}\beta_{j}^{2}$
$\displaystyle=\left(\sqrt{(L-l+1)\sum_{j>l}\beta_{j}^{2}}\right)^{2}$
$\displaystyle\mbox{Cauchy-Schwarz}~{}\Rightarrow\hskip 28.45274pt$
$\displaystyle\geq\left(\sum_{j>l}\beta_{j}\right)^{2}\,,$

to derive

$\displaystyle{\mathcal{N}}^{1}(\boldsymbol{\beta})$
$\displaystyle=\sum_{l=1}^{L}\sqrt{L-l+1}|\\!|\beta_{\geq l}|\\!|_{2}$
$\displaystyle\geq\sum_{l=1}^{L}\left(\sqrt{L-l+1}|\beta_{l}|+\sum_{j>l}|\beta_{j}|\right)$
$\displaystyle=\sum_{l=1}^{L}\left(\sqrt{L-l+1}+l-1\right)|\beta_{l}|$
$\displaystyle\geq\sqrt{L}|\beta_{l}|\,.$

Since
${\mathcal{N}}^{1}_{\star}(\boldsymbol{\alpha})\leq|\\!|\boldsymbol{\alpha}|\\!|_{\infty}{\mathcal{N}}^{1}(\boldsymbol{1})$,
where $\boldsymbol{1}\in\mathbb{R}^{L}$ is the all-one vector, it suffices to
show that ${\mathcal{N}}^{1}_{\star}(\boldsymbol{1})\leq L^{-\frac{1}{2}}$.
Since

${\mathcal{N}}^{1}_{\star}(\boldsymbol{1})=\max_{\mathcal{N}^{1}(\boldsymbol{\beta})=1}\sum_{j\in[L]}\beta_{j}\,,$

and

$\displaystyle 1$ $\displaystyle=\mathcal{N}^{1}(\boldsymbol{\beta})$
$\displaystyle=\sum_{j\in[L]}\sqrt{L-l+1}|\\!|\beta_{\geq j}|\\!|_{2}$
$\displaystyle\geq\sqrt{L}\sum_{j\in[L]}|\beta_{j}|\,,$

the claim follows. ∎

With the above proposition dealing with the case $M=1$, we can formulate in
general the following proposition.

######  Proposition D.3 ($L^{\infty}$ norm bounds dual norm, general case).

Suppose $\boldsymbol{\alpha}\in\mathbb{R}^{ML}$, and

$\displaystyle{\mathcal{N}}(\boldsymbol{\alpha})$
$\displaystyle=\sum_{l=1}^{L}\sqrt{M(L-l+1)}|\\!|\boldsymbol{\alpha}_{(\geq
l)}|\\!|_{2}\,,$ $\displaystyle\mbox{where}\hskip
28.45274pt|\\!|\boldsymbol{\alpha}_{(\geq l)}|\\!|_{2}:$
$\displaystyle=\sqrt{\sum_{j=l}^{L}\sum_{m=1}^{M}\alpha^{2}_{(m-1)+j}}\,.$

Then ${\mathcal{N}}_{\star}(\boldsymbol{\alpha})\leq
L^{-\frac{1}{2}}|\\!|\boldsymbol{\alpha}|\\!|_{\infty}$.

###### Proof of Proposition D.3.

Again, it suffices to show that ${\mathcal{N}}(\boldsymbol{\alpha})\leq
L^{-\frac{1}{2}}$ when $\boldsymbol{\alpha}$ is the all-one vector.

Now, ${\mathcal{N}}_{\star}(\boldsymbol{\alpha})$ is the optimum value of the
problem

maximize $\displaystyle\sum_{j\in[ML]}\beta_{j}$ subject to
$\displaystyle{\mathcal{N}}(\boldsymbol{\beta})=1\,.$

Write $\boldsymbol{\beta}_{(l)}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(\beta_{l},\beta_{L+l},\cdots,\beta_{(M-1)L+l})^{\top}$
and $x_{l}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\sqrt{M}|\\!|\boldsymbol{\beta}_{(l)}|\\!|_{2}$
for $l\in[L]$; by Cauchy-Schwarz, we have

$\displaystyle\sum_{j\in[ML]}\beta_{j}$
$\displaystyle\leq\sum_{l=1}^{L}x_{l}\,.$

Moreover,
${\mathcal{N}}(\boldsymbol{\beta})={\mathcal{N}}^{1}(\boldsymbol{x})$, where
$\boldsymbol{x}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(x_{1},\cdots,x_{L})$;
thus, the lemma follows by D.2. ∎

##  Appendix E Some known results used in the proofs

The following tail bound on the norm of Gaussian random vectors is well-known;
see (Vershynin, 2018, Theorem 6.2.1) for details.

######  Proposition E.1 (Wright-Hansen inequality).

There is an absolute constant $c_{0}>0$ such that the following statement
holds. If $\boldsymbol{q}\sim N(\boldsymbol{0},\Sigma)$ in $\mathbb{R}^{d}$,
where $\Sigma$ is symmetric positive definite, then the following holds for
any $\tau>0$:

$\displaystyle\mathbb{P}\left[\frac{1}{d}\left|\boldsymbol{q}^{T}\boldsymbol{q}-\text{tr}(\Sigma)\right|\geq
4\tau|\\!|\Sigma|\\!|_{\text{op}}\right]$ $\displaystyle\leq
2e^{-c_{0}d\min\\{\tau^{2},\tau\\}}\,.$

We need to collect some information about the time series model in equation
(3.1). Let $T\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\min_{m\in[M]}T_{m}$,
and recall that we have written $D=T_{1}+\cdots+T_{M}$. Let $\Gamma^{(m,l)}$
be the covariance matrix of the random vector $X^{(m,l)}$; this is a symmetric
matrix, and for any $r,s\in\\{1-l,\dots,T_{m}-l\\}$ with $r\leq s$, the
stationarity of the time series yields

$\Gamma^{(m,l)}_{r,s}=\mathbb{E}[X^{(m,l)}_{r-l}X^{(m,l)}_{s-l}]=\mathbb{E}[X^{m}_{0}X^{m}_{s-r}]\,.$
(59)

Thus, we can write $\Gamma^{(m)}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Gamma^{(m,l)}$
for the “auto-covariance” matrix.

Now, let further, for any two integers $T>L>0$, ${}_{T,L}\Gamma$ denote the
$T_{1}L\times T_{1}L$ matrix, whose $(r,s)$-th block is the covariance of the
vectors $(X_{r-1},X_{r-2},\dots,X_{r-L})$ and
$(X_{s-1},X_{s-2},\dots,X_{s-L})$.

The following lemma gives eigenvalue bounds of these block covariance
matrices.

######  Lemma E.2 (Eigenvalue bound for Toeplitz matrices).

Let $f(z)$ denote the reverse characteristic polynomial of a stable univariate
AR-process $\\{X_{t}\\}$. Let $\mathfrak{m}(f)=\inf_{|z|\leq 1}|f(z)|^{2}$ and
$\mathfrak{M}(f)=\sup_{|z|\leq 1}|f(z)|^{2}$. The following inequality holds:

$\frac{\Lambda_{\min}(\Sigma_{\epsilon})}{\mathfrak{M}(f)}\leq\Lambda_{\min}(_{T,L}\Gamma)\leq\Lambda_{\max}(_{T,L}\Gamma)\leq\frac{\Lambda_{\max}(\Sigma_{\epsilon})}{\mathfrak{m}(f)}\,.$

###### Proof of Lemma E.2.

This is immediate from (Basu and Michailidis, 2015, Proposition 2.3) if we
consider the stable $\dim$-L VAR process $Y_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=(X_{t-1},\dots,X_{t-L})$.
∎

The following bound appeared as of (Basu and Michailidis, 2015, Proposition
2.4b).

######  Lemma E.3 (Basu-Michailidis).

Let $\Delta^{(m)}$ denote the $T_{m}$-dimensional matrix with
$\Delta_{t,s}^{(m)}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{E}[X^{m}_{t-l}U^{m}_{s}]$.
The following inequality holds:

$\Lambda_{\max}(\Delta^{(m)})\leq\frac{\Lambda_{\max}(\Sigma_{\epsilon_{m}})\mathfrak{M}(f_{m})}{\mathfrak{m}(f)}\,.$

The following appeared as (Basu and Michailidis, 2015, Lemma F.2). Recall that
$\mathbb{M}_{L}(\mathbb{R})$ denotes the algebra of all $L\times L$ real
matrices, and $\mathbb{S}^{L-1}$ the unit sphere in $\mathbb{R}^{L}$.

######  Lemma E.4 (Basu-Michailidis).

Suppose that $G\in\mathbb{M}_{L}(\mathbb{R})$ is a random symmetric matrix for
which the following inequality holds for every
$\boldsymbol{u}\in\mathbb{S}^{L-1}$, $T_{1}\in\mathbb{N}$, and $\eta>0$:

$\mathbb{P}\left[\left|\boldsymbol{u}^{T}G\boldsymbol{u}\right|\geq
C\eta\right]\leq 2e^{-cT_{1}\min\\{\eta,\eta^{2}\\}}\,,$

where $C,c>0$ are parameters independent of $\eta$ and $\boldsymbol{u}$. Then
the following inequality holds for any integer $s_{0}>0$:

$\displaystyle\mathbb{P}\left[\sup_{|\\!|\boldsymbol{u}|\\!|_{0}\leq
s_{0},~{}|\\!|\boldsymbol{u}|\\!|_{2}\leq
1}\left|\boldsymbol{u}^{T}G\boldsymbol{u}\right|\geq C\eta\right]$
$\displaystyle\leq 2e^{-cT_{1}\min\\{\eta,\eta^{2}\\}+s_{0}\min\\{\log
L,\log(21eL/s_{0})\\}}\,.$ (60)

The following result was obtained as (Loh and Wainwright, 2012, Lemma 12).
Recall that, for any integer $s_{0}>0$, we denote

$\mathbb{K}(s_{0})=\\{\boldsymbol{v}\in\mathbb{R}^{L}:|\\!|\boldsymbol{v}|\\!|_{2}\leq
1,~{}|\\!|\boldsymbol{v}|\\!|_{0}\leq s_{0}\\}\,.$

######  Lemma E.5 (Loh-Wainwright).

Suppose that $G\in\mathbb{M}_{L}(\mathbb{R})$ is a symmetric matrix for which
the following holds:

$\sup_{|\\!|\boldsymbol{u}|\\!|_{0}\leq
s_{0},~{}|\\!|\boldsymbol{u}|\\!|_{2}\leq
1}\left|\boldsymbol{u}^{T}G\boldsymbol{u}\right|\leq\delta\,.$

Then

$\sup_{\boldsymbol{u}\in\mathbb{R}^{L}}\left|\boldsymbol{u}^{T}G\boldsymbol{u}\right|\leq
27\delta\left(|\\!|\boldsymbol{u}|\\!|_{2}^{2}+\frac{2}{s_{0}}|\\!|\boldsymbol{u}|\\!|_{1}^{2}\right)\,.$
(61)