## 1

Hyeyoung Maeng, Idris Eckley and Paul Fearnhead
Lancaster University, United Kingdom

## Collective anomaly detection in High-dimensional VAR Models

Hyeyoung Maeng, Idris Eckley and Paul Fearnhead
Lancaster University, United Kingdom

###### Abstract

There is increasing interest in detecting collective anomalies: potentially
short periods of time where the features of data change before reverting back
to normal behaviour. We propose a new method for detecting a collective
anomaly in VAR models. Our focus is on situations where the change in the VAR
coefficient matrix at an anomaly is sparse, i.e. a small number of entries of
the VAR coefficient matrix change. To tackle this problem, we propose a test
statistic for a local segment that is built on the lasso estimator of the
change in model parameters. This enables us to detect a sparse change more
efficiently and our lasso-based approach becomes especially advantageous when
the anomalous interval is short. We show that the new procedure controls Type
1 error and has asymptotic power tending to one. The practicality of our
approach is demonstrated through simulations and two data examples, involving
New York taxi trip data and EEG data.

Keywords: Collective anomaly; high-dimensional time series; lasso; sparse
changes; epidemic change

##  2 Introduction

There is a growing need for modelling and analysis of high-dimensional time
series, as such series have become increasingly common in many application
areas. Applications include forecasting using a large panel of time series in
economics (De Mol et al.,, 2008; Bańbura et al.,, 2010), estimating causal
relationships among genes and constructing gene regulatory networks (Shojaie
and Michailidis,, 2010), identifying the monetary transmission mechanism from
macroeconomic time series (Bernanke et al.,, 2005), discovering causal
interactions in Neuroimaging (Smith,, 2012; Seth et al.,, 2015), analysing
housing markets for home-price estimation and forecasting (Rapach et al.,,
2007; Calomiris et al.,, 2008; Stock and Watson,, 2008) and analysing the
network structure of volatility interconnections in the Standard & Poor’s 100
data (Barigozzi and Hallin,, 2017).

The majority of existing methods are built on the assumption of stationary and
stable time series. If there is either a structural change or a period of
anomalous behaviour in a time series, detecting the location of the
change/anomaly is not only an important task in itself, but also useful for a
follow-up analysis after detection. It is indeed a problem of significant
interest in many applications. For example, Chen and Gupta, (1997) detects
multiple change-points in variance of weekly stock prices and Cribben and Yu,
(2017) study a network change-point detection problem for resting state
functional magnetic resonance imaging data. Ombao et al., (2005) propose a way
of segmenting multivariate nonstationary time series and analyse time-varying
electroencephalogram data that is recorded during an epileptic seizure. Other
examples include detecting changes that have occurred in a sparse subset of
time series (Cho and Fryzlewicz,, 2015; Cho,, 2016; Wang and Samworth,, 2018),
covariance change-point detection for multivariate or high-dimensional time
series (Inclan and Tiao,, 1994; Aue et al.,, 2009; Wang et al.,, 2017) and
detecting change-points under the factor model framework (Breitung and
Eickmeier,, 2011; Chen et al.,, 2014; Baltagi et al.,, 2017; Barigozzi et
al.,, 2018).

One of the most popular models for high-dimensional time series is the vector
autoregressive (VAR) model (Sims,, 1980; Lütkepohl,, 2005), due to its ability
to capture complex temporal and cross-sectional relationships. However, the
estimation of the coefficient matrix becomes challenging as the number of
parameters increases quadratically with the number of time series. To overcome
this, structured sparsity of the VAR coefficients is often assumed as this
assumption dramatically reduces the number of model parameters. For example,
Song and Bickel, (2011) use lasso type, that is $\ell_{1}$, penalties to
encourage sparsity in the estimates of the VAR coefficients. Davis et al.,
(2016) propose a two-stage approach to fit sparse VAR models and provide a
numerical evidence that a log-likelihood based loss function improves the
forecasting performance compared to a least squares based one as the former
includes information on the error covariance matrix. Basu and Michailidis,
(2015) investigate the theoretical properties of $\ell_{1}$-penalised
estimators for a Gaussian VAR model and show consistency results, while Lin
and Michailidis, (2017) generalise the results by considering a general norm
instead of being restricted to the $\ell_{1}$-norm for the penalty. Recently,
more complex structures have been studied in the literature: Basu et al.,
(2019) study the low-rank and structured sparse VAR model and Nicholson et
al., (2020) impose a hierarchical structure on VAR coefficient matrices
according to the lag order, thus addressing both the dimensionality and the
lag selection issues at the same time.

Despite the large body of literature on VAR models, detecting a structural
break has rarely been studied. Kirch et al., (2015) consider two scenarios,
detecting at-most-one-change and epidemic change in model parameters of
multivariate time series which is not restricted to VAR models. Safikhani and
Shojaie, (2020) consider the multiple change-point setting for the VAR
coefficient matrix under a high-dimensional regime and propose a three-stage
procedure that returns consistent estimators of both change-points and
parameters. Wang et al., (2019) also study the same setting (i.e. when the
model parameters have a form of piecewise constant over time) and use a
dynamic programming approach for localising change-points and improving the
corresponding error rates. Bai et al., (2020) study the multiple change-point
setting but assume the low-rank plus sparse structure on the VAR coefficient
matrices and consider the case where only the sparse structure changes over
time, while the low-rank parts remain constant. We will explain how our
proposal is different from these existing works later in this section.

In contrast to these earlier works, we focus on settings where we have plenty
of information about the current or normal behaviour of our time-series, and
wish to detect periods of different or abnormal behaviour. This can arise when
detecting collective anomalies (Fisch et al.,, 2018; Tveten et al.,, 2020) or
epidemic changepoints (Yao,, 1993) – where we have a, potentially short,
period of time where the behaviour of our model changes before it reverts back
to current behaviour. This also arises with sequential change detection (Lai,,
1995), when we observe data in real-time and wish to detect any change away
from the current behaviour as quickly as possible. For ease of presentation we
focus primarily on detecting collective anomalies/epidemic changes, and use
the terminology collective anomaly from now on. We show how our method can be
extended to the online framework in Section 6. The key feature of these
problems is that we have substantially more information about the current or
normal behaviour than about the anomaly. This suggests that we should
potentially use different procedures to estimate the parameters of the VAR
model for the normal behaviour than for the anomaly. We do this through making
an assumption that it is the change in VAR parameters that is sparse.

We focus on improving the detection power when the difference between the
coefficient matrices at anomaly point is sparse (i.e. a small number of
entries of the VAR coefficient matrix change). To tackle this problem, we
propose a test statistic for a local segment which is built on the lasso
estimator of the change in model parameters. This enables us to detect a
sparse change more efficiently, as the sparsity of change is considered in
establishing the test statistic. Moreover, our lasso-based approach become
more advantageous over, say, the standard likelihood-ratio test statistic for
shorter anomalous intervals: as for shorter intervals we have fewer
observations to estimate the new VAR coefficient matrix, and it becomes more
like a high-dimensional problem where the number of observations is similar to
or less than the number of parameters to estimate. In Section 5, our approach
is compared with a method that is built on estimating the change in VAR matrix
using ordinary least squares estimator, and the results show that our method
outperforms it when detecting sparse change. As we consider the setting where
a relatively longer region has a normal behaviour than the anomalous
behaviour, it is reasonable to assume that the underlying VAR coefficient
matrix is estimated well enough. Thus, we first develop our method when the
normal behaviour is assumed to be known and extend it to the case where an
appropriate estimator for the VAR coefficient is used instead. Our theory in
Section 4 shows the validity of this approach providing that the estimator for
the VAR coefficient is close enough to the true one. Although our main focus
is on single anomaly detection, we show that the new method can be extended
for detecting multiple anomalies in Section 3.1.

Among those relevant works already introduced earlier in this section, the
work of Safikhani and Shojaie, (2020) and Bai et al., (2020) are most closely
related to our work, in that they also control the change in VAR parameters
with a lasso penalty in their objective functions, however their approaches
are different from our method in several aspects. To obtain the initial
estimate of change-points before screening, Safikhani and Shojaie, (2020) use
a fused lasso penalty on a full model considering all time points being a
candidate for change-point. Thus their objective function controls the
sparsity of VAR parameters and the sparsity of its difference at the same
time. Bai et al., (2020) follows a similar procedure to Safikhani and Shojaie,
(2020) under the multiple change-point framework. They use a block fused lasso
penalty by assuming that the model parameters in a block is fixed, while our
objective function controls only the sparsity of change in building a test
statistic and search many segments to find an anomalous interval. Also,
Safikhani and Shojaie, (2020) and Bai et al., (2020) assume that the
$l_{2}$-norm of a change in VAR parameter is bounded away from zero, whereas
our assumption on the $l_{2}$-norm of a change is related to the sparsity of
change which is in line with the assumptions used in Wang et al., (2019).
Although those change-point detection methods are not exactly designed for the
anomaly setting we consider in this paper, we compare our performance with
theirs and present results in the supplementary material. Our method works
better especially when the underlying VAR coefficient matrix is dense but the
change is sparse, and surprisingly even in the case where the VAR coefficient
matrix has a low rank plus sparse structure and only a sparse component
changes. Full details can be found in the supplementary material.

The remainder of the article is organised as follows. Section 3 gives a full
description of our procedure and the relevant theoretical results are
presented in Section 4. The supporting simulation studies are described in
Section 5. Our methodology is illustrated through two datasets in Section 6
and we end with additional discussion in Section 7. The proofs of our main
theoretical results are in the supplementary material.

##  3 Methodology

###  3.1 Problem setting

We consider a zero-mean stationary ${p}$-dimensional multivariate time series
$\boldsymbol{x}_{t}=(x_{1t},\ldots,x_{pt})^{\prime}$ generated by a VAR(1)
model:

$\boldsymbol{x}_{t}=\boldsymbol{A}_{t}\boldsymbol{x}_{t-1}+{\boldsymbol{\varepsilon}}_{t},\quad\boldsymbol{\varepsilon}_{t}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}N(\boldsymbol{0},\Sigma_{\varepsilon}),\quad
t=1,\ldots,T,$ (1)

where $\\{\boldsymbol{A}_{t}\\}_{t=1}^{T}$ is a $p\times p$ matrix and
$\Sigma_{\varepsilon}$ is a positive definite matrix. We assume that the high-
dimensional VAR model shows an anomalous behaviour at
$t\in[\eta_{1},\eta_{2}]$ such that

$0=\eta_{0}<\eta_{1}<\eta_{2}<\eta_{3}=T,$ (2)

which gives the sets

$\boldsymbol{\mathfrak{x}}_{1}=\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{\eta_{1}-1}\\},\quad\boldsymbol{\mathfrak{x}}_{2}=\\{\boldsymbol{x}_{\eta_{1}},\ldots,\boldsymbol{x}_{\eta_{2}}\\},\quad\boldsymbol{\mathfrak{x}}_{3}=\\{\boldsymbol{x}_{\eta_{2}+1},\ldots,\boldsymbol{x}_{T}\\}$
(3)

and the sequence of $\\{\boldsymbol{A}_{t}\\}_{t=1}^{T}$ forms piecewise-
constant coefficient matrices as follows,

$\displaystyle\boldsymbol{A}^{(1)}=\boldsymbol{A}_{1}=\cdots=\boldsymbol{A}_{\eta_{1}-1},\quad\boldsymbol{A}^{(2)}=\boldsymbol{A}_{\eta_{1}}=\cdots=\boldsymbol{A}_{\eta_{2}},\quad\boldsymbol{A}^{(1)}=\boldsymbol{A}_{\eta_{2}+1}=\cdots\boldsymbol{A}_{T},$

where $\boldsymbol{A}^{(1)}\neq\boldsymbol{A}^{(2)}$ and
$\boldsymbol{A}^{(1)},\boldsymbol{A}^{(2)}\in\mathbb{R}^{p\times p}$.

The model in equation (1) can be represented as the following linear
regression,

$\begin{pmatrix}\boldsymbol{x}_{1}^{\prime}\\\ \boldsymbol{x}_{2}^{\prime}\\\
\vdots\\\ \boldsymbol{x}_{T}^{\prime}\end{pmatrix}_{T\times
p}=\begin{pmatrix}\boldsymbol{x}_{0}^{\prime}&0\\\ \vdots&\vdots\\\
\boldsymbol{x}_{\eta_{1}-2}^{\prime}&0\\\
\boldsymbol{x}_{\eta_{1}-1}^{\prime}&\boldsymbol{x}_{\eta_{1}-1}^{\prime}\\\
\vdots&\vdots\\\
\boldsymbol{x}_{\eta_{2}-1}^{\prime}&\boldsymbol{x}_{\eta_{2}-1}^{\prime}\\\
\boldsymbol{x}_{\eta_{2}}^{\prime}&0\\\ \vdots&\vdots\\\
\boldsymbol{x}_{T-1}^{\prime}&0\end{pmatrix}_{T\times
2p}\begin{pmatrix}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
{\boldsymbol{\theta}^{(2)}}^{\prime}\end{pmatrix}_{2p\times
p}+\begin{pmatrix}\varepsilon_{1}^{\prime}\\\ \varepsilon_{2}^{\prime}\\\
\vdots\\\ \varepsilon_{T}^{\prime}\end{pmatrix}_{T\times p},$ (4)

where

$\displaystyle\boldsymbol{\theta}^{(1)}=\boldsymbol{A}^{(1)},\;\boldsymbol{\theta}^{(2)}=\boldsymbol{A}^{(2)}-\boldsymbol{A}^{(1)}.$

The model, as written in equation (4), is a linear regression of the form
$\mathcal{Y}=\mathcal{X}\Theta+{E}.$ As such, it can be represented as

$\boldsymbol{Y}_{Tp\times 1}=\boldsymbol{X}_{Tp\times
2p^{2}}\boldsymbol{\Theta}_{2p^{2}\times 1}+\boldsymbol{E}_{Tp\times 1},$ (5)

where $\boldsymbol{X}=\mathit{I}_{p}\otimes\mathcal{X}$ and $\otimes$ is the
tensor product of two matrices.

Now our interest is in estimating the collective anomaly
$[\eta_{1},\eta_{2}]$. Our motivation is for scenarios where there is
substantial information about the normal or pre-change behaviour of the data.
Thus, for ease of presentation, we will first assume that
$\boldsymbol{\theta}^{(1)}$ in (4) is known. In practice we will use an
estimate of $\boldsymbol{\theta}^{(1)}$, and our theory shows that our
approach has good asymptotic properties if we plug-in a suitably accurate
estimate of $\boldsymbol{\theta}^{(1)}$ in the following procedure. We assume
that the change $\boldsymbol{\theta}^{(2)}$ is sparse in that it has small
number of nonzero entries which will be formulated in a later section.
Assuming the base coefficient matrix $\boldsymbol{A}^{(1)}$ is known, we can
rewrite the model as

$\begin{pmatrix}\boldsymbol{x}_{1}^{\prime}\\\ \boldsymbol{x}_{2}^{\prime}\\\
\vdots\\\ \boldsymbol{x}_{T}^{\prime}\end{pmatrix}_{T\times
p}-\begin{pmatrix}\boldsymbol{x}_{0}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
\vdots\\\
\boldsymbol{x}_{\eta_{1}-2}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
\boldsymbol{x}_{\eta_{1}-1}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
\vdots\\\
\boldsymbol{x}_{\eta_{2}-1}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
\boldsymbol{x}_{\eta_{2}}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
\vdots\\\
\boldsymbol{x}_{T-1}^{\prime}{\boldsymbol{\theta}^{(1)}}^{\prime}\end{pmatrix}_{T\times
p}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\ \boldsymbol{x}_{\eta_{1}-1}^{\prime}\\\
\vdots\\\ \boldsymbol{x}_{\eta_{2}-1}^{\prime}\\\ 0\\\ \vdots\\\
0\end{pmatrix}_{T\times
p}\begin{pmatrix}{\boldsymbol{\theta}^{(2)}}^{\prime}\end{pmatrix}_{p\times
p}+\begin{pmatrix}\varepsilon_{1}^{\prime}\\\ \varepsilon_{2}^{\prime}\\\
\vdots\\\ \varepsilon_{T}^{\prime}\end{pmatrix}_{T\times p},$ (6)

which can be represented as
$\mathcal{Y}-\mathcal{X}^{(1)}{\boldsymbol{\theta}^{(1)}}^{\prime}=\mathcal{X}^{(2)}{\boldsymbol{\theta}^{(2)}}^{\prime}+{E}.$
With slight abuse of notation as we are using different definitions of
$\boldsymbol{Y}$, $\boldsymbol{X}$ and $\boldsymbol{\Theta}$, this can be re-
written as

$\boldsymbol{Y}_{Tp\times 1}=\boldsymbol{X}_{Tp\times
p^{2}}\boldsymbol{\Theta}_{p^{2}\times 1}+\boldsymbol{E}_{Tp\times 1},$ (7)

where $\boldsymbol{X}=\mathit{I}_{p}\otimes\mathcal{X}^{(2)}$.

###  3.2 Lasso-based approach

To detect a collective anomaly we derive a test for whether data in an
interval of time is anomalous, and then apply this test to data from a set of
suitably chosen intervals, $\mathbb{J}_{T,p}(L)$. To help with the
presentation of theory in Section 4, we parameterise this set by the length,
$L$, of the smallest interval it contains. For any interval
$J\in\mathbb{J}_{T,p}(L)$, by extracting the corresponding rows from each
matrix in (6), the linear regression form can be rewritten as:
$\mathcal{Y}_{J}-\mathcal{X}^{(1)}_{J}{\boldsymbol{\theta}^{(1)}}^{\prime}=\mathcal{X}^{(2)}_{J}{\boldsymbol{\theta}^{(2)}}^{\prime}+{E}_{J},$
that can be vectorised in a form of

$\boldsymbol{Y}_{J}=\boldsymbol{X}_{J}\boldsymbol{\Theta}+\boldsymbol{E}_{J},$
(8)

as in (7).

One of the standard ways to detect change or epidemic changes in regression
models is to use a likelihood ratio test (Kim and Siegmund,, 1989; Siegmund
and Venkatraman,, 1995; Yau and Zhao,, 2016; Baranowski et al.,, 2019; Dette
and Gösmann,, 2020), and these methods can be applied in the VAR setting. To
detect a collective anomaly in a set of intervals, our procedure involves
calculating the likelihood ratio statistic for each interval
$J\in\mathbb{J}_{T,p}(L)$ as

$-2\bigg{\\{}\sum_{s\in
J}l_{s}\big{(}\boldsymbol{\Theta}=0,\Sigma_{\varepsilon}\big{)}-\sum_{s\in
J}l_{s}\big{(}\boldsymbol{\hat{\boldsymbol{\Theta}}},\Sigma_{\varepsilon}\big{)}\bigg{\\}},$
(9)

where $\boldsymbol{\hat{\boldsymbol{\Theta}}}$ is the maximum likelihood
estimator and the likelihood function has the form of

$\sum_{s\in
J}l_{s}\big{(}\boldsymbol{\boldsymbol{\Theta}},\Sigma_{\varepsilon}\big{)}=-\frac{1}{2}\bigg{\\{}|J|p\log(2\pi)+|J|\log|\Sigma_{\varepsilon}|+(\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}})^{\top}(\Sigma_{\varepsilon}^{-1}\otimes\mathit{I})(\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}})\bigg{\\}}.$

As we consider only ${\boldsymbol{\Theta}}$ varying, the first two terms are
constant and will cancel in the test statistic.

It is common to assume $\Sigma_{\varepsilon}$ is the identity matrix, in which
case the maximum likelihood estimator of $\boldsymbol{\Theta}$ is the ordinary
least squares (OLS) estimator. Alternatively we can estimate the variance from
the residuals obtained when estimating the parameters of the VAR model on
training data. For ease of presentation, we will assume $\Sigma_{\varepsilon}$
is an identity matrix from now on, but our theoretical results are still valid
if this assumption is not correct. Furthermore, the theory can be extended to
situations where we assume either $\Sigma_{\varepsilon}$ is any positive
identity matrix or an estimate of $\Sigma_{\varepsilon}$ is used. We now give
details of the likelihood ratio statistic and our suggested improvement based
on penalised estimation of the change in the VAR coefficients.

##### The OLS method

Before introducing the lasso-based approach, we consider the test statistic
based on the least squares estimator which we refer to as the OLS method. The
OLS estimator has been popularly used in the change point detection literature
e.g. in a linear model setup, CUSUM-type approaches built on the least squares
estimator are studied by Horváth et al., (2004), Aue et al., (2006), Chen and
Tian, (2010) and Fremdt, (2015). For any interval $J\in\mathbb{J}_{T,p}(L)$,
the test statistic of the OLS method takes the form,

$\displaystyle T(J)$
$\displaystyle=\|\boldsymbol{Y}_{J}\|_{2}^{2}-\min_{\boldsymbol{\Theta}}\big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}\big{\\}}$
$\displaystyle=\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2},$
(10)

that is the same as the likelihood ratio statistic in (9) when
$\Sigma_{\varepsilon}$ is the identity matrix. $T(J)$ has a $\chi^{2}_{p^{2}}$
distribution under the null, $\boldsymbol{\Theta}=\boldsymbol{0}$. The
classical least squares estimator $\hat{\boldsymbol{\Theta}}$ in (10) is not
able to be used when the dimension $p$ is greater than $T$. Note that
$\hat{\boldsymbol{\Theta}}$ also depends on $J$ but this is suppressed in the
notation for simplicity.

##### The Lasso method

To handle the case when $\boldsymbol{\Theta}$ is sparse more effectively, we
propose a test statistic based on a lasso estimator:

$T^{\text{lasso}}(J)=\|\boldsymbol{Y}_{J}\|_{2}^{2}-\min_{\boldsymbol{\Theta}}\big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}+\lambda\|\boldsymbol{\Theta}\|_{1}\big{\\}}.$
(11)

To detect a collective anomaly, we calculate this test statistic for a
collection of intervals, $\mathbb{J}_{T,p}(L)$. We detect an anomaly if the
maximum value of these test statistics is above a pre-determined threshold. If
we detect an anomaly, we estimate its location as the interval in
$\mathbb{J}_{T,p}(L)$ with the largest test-statistic value. The detailed
procedure is given in Algorithm 1.

INPUT: $\boldsymbol{X}$ matrix in (7), $L$, $\lambda^{\textsuperscript{thr}}$

  1. Step 1: Set a collection of intervals $\mathbb{J}_{T,p}(L)$ where $L$ is the minimum length of intervals.

  2. Step 2: For any interval $J\in\mathbb{J}_{T,p}(L)$, calculate $T^{\text{lasso}}(J)$ as in (11).

  3. Step 3: Using a pre-specified threshold $\lambda^{\textsuperscript{thr}}$, pick the candidate set

$\mathbb{I}^{*}=\Big{\\{}J\in\mathbb{J}_{T,p}(L):T^{\text{lasso}}(J)>\lambda^{\textsuperscript{thr}}\Big{\\}}.$

If $\mathbb{I}^{*}\neq\emptyset$, reject the null hypothesis (no anomaly
exists) and save the estimator of the anomaly interval,
$\hat{I}=\operatorname*{\arg\\!\max}_{J\in\mathbb{J}_{T,p}(L)}T^{\text{lasso}}(J).$
(12)

OUTPUT: $\hat{I}$.

Algorithm 1 Single anomaly detection

For setting the collection of intervals $\mathbb{J}_{T,p}(L)$ in Step 1, there
exist two general methods; randomly generated intervals (Fryzlewicz,, 2014;
Baranowski et al.,, 2019) and deterministic construction of intervals (Kovács
et al.,, 2020). In this paper, we use both construction methods and compare
their performance in Section 5.

###  3.3 Extension to detecting multiple anomalies

Following the ideas in Fryzlewicz, (2014) and Kovács et al., (2020), to deal
with multiple anomalies, we repeatedly update the candidate set by removing
the intervals that overlap with any detected anomalies. The detailed procedure
is given in Algorithm 2.

INPUT: $\boldsymbol{X}$ matrix in (7), $L$, $\lambda^{\textsuperscript{thr}}$

  1. Step 1: Set a collection of intervals $\mathbb{J}_{T,p}(L)$ where $L$ is the minimum length of intervals.

  2. Step 2: For any interval $J\in\mathbb{J}_{T,p}(L)$, calculate $T^{\text{lasso}}(J)$ as in (11).

  3. Step 3: Using a pre-specified threshold $\lambda^{\textsuperscript{thr}}$, pick the candidate set

$\mathbb{I}^{*}=\Big{\\{}J\in\mathbb{J}_{T,p}(L):T^{\text{lasso}}(J)>\lambda^{\textsuperscript{thr}}\Big{\\}}.$

If $\mathbb{I}^{*}\neq\emptyset$, reject the null hypothesis (no anomaly
exist). Set $\mathbb{I}^{(1)}=\mathbb{I}^{*}$, $j=1$ and proceed the following
steps.

while  _$\mathbb{I}^{(j)}\neq\emptyset$_ do

  1. Step 4: Save the estimator of the anomaly interval,

$\hat{I}_{j}=\operatorname*{\arg\\!\max}_{J\in\mathbb{I}^{(j)}}T^{\text{lasso}}(J),$

and update the candidate set as

$\mathbb{I}^{(j+1)}=\mathbb{I}^{(j)}\setminus\\{J:J\in\mathbb{I}^{(j)},J\cap\hat{I}_{j}\neq\emptyset\\}$

  2. Step 5: Set $j=j+1$.

end while

OUTPUT: $\hat{I}=\\{\hat{I}_{1},\hat{I}_{2},\cdots\\}$.

Algorithm 2 Multiple anomaly detection

###  3.4 Extension to VAR(q) model

The VAR process of order $1$ presented in Section 3.1 can simply be extended
to VAR(q) as follows,

$\boldsymbol{x}_{t}=\boldsymbol{A}_{t,1}\boldsymbol{x}_{t-1}+\cdots+\boldsymbol{A}_{t,q}\boldsymbol{x}_{t-q}+{\boldsymbol{\varepsilon}}_{t},\quad\boldsymbol{\varepsilon}_{t}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}N(\boldsymbol{0},\Sigma_{\varepsilon}),\quad
t=1,\ldots,T,$ (13)

where $\\{\boldsymbol{A}_{t,k}\\}_{t=1}^{T}$ is a $p\times p$ matrix for all
$k=1,\ldots,q$ and $\Sigma_{\varepsilon}$ is assumed to be a positive definite
matrix. With a slight abuse of notation, the piecewise-constant coefficient
matrices are as follows,

$\displaystyle\boldsymbol{A}^{(1)}=(\boldsymbol{A}_{t^{\prime},1},\ldots,\boldsymbol{A}_{t^{\prime},q})\in\mathbb{R}^{p\times
pq},\quad\text{ for any
}t^{\prime}=1,\ldots,{\eta_{1}-1},{\eta_{2}+1},\ldots,T$
$\displaystyle\boldsymbol{A}^{(2)}=(\boldsymbol{A}_{t^{\prime},1},\ldots,\boldsymbol{A}_{t^{\prime},q})\in\mathbb{R}^{p\times
pq},\quad\text{ for any }t^{\prime}={\eta_{1}},\ldots,{\eta_{2}},$

and the model (13) can be represented as

$\begin{pmatrix}\boldsymbol{x}_{q}^{\prime}\\\
\boldsymbol{x}_{q+1}^{\prime}\\\ \vdots\\\
\boldsymbol{x}_{T}^{\prime}\end{pmatrix}_{(T-q+1)\times
p}=\begin{pmatrix}\boldsymbol{x}_{q-1}^{\prime}&\cdots&\boldsymbol{x}_{0}^{\prime}&0&\cdots&0\\\
\vdots&&\vdots&\vdots&&\vdots\\\
\boldsymbol{x}_{\eta_{1}+q-3}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{1}-2}^{\prime}&0&\cdots&0\\\
\boldsymbol{x}_{\eta_{1}+q-2}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{1}-1}^{\prime}&\boldsymbol{x}_{\eta_{1}+q-2}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{1}-1}^{\prime}\\\
\vdots&&\vdots&\vdots&&\vdots\\\
\boldsymbol{x}_{\eta_{2}+q-2}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{2}-1}^{\prime}&\boldsymbol{x}_{\eta_{2}+q-2}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{2}-1}^{\prime}\\\
\boldsymbol{x}_{\eta_{2}+q-1}^{\prime}&\cdots&\boldsymbol{x}_{\eta_{2}}^{\prime}&0&\cdots&0\\\
\vdots&&\vdots&\vdots&&\vdots\\\
\boldsymbol{x}_{T-1}^{\prime}&\cdots&\boldsymbol{x}_{T-q}^{\prime}&0&\cdots&0\\\
\end{pmatrix}_{(T-q+1)\times
2pq}\begin{pmatrix}{\boldsymbol{\theta}^{(1)}}^{\prime}\\\
{\boldsymbol{\theta}^{(2)}}^{\prime}\end{pmatrix}_{2pq\times
p}+\begin{pmatrix}\varepsilon_{q}^{\prime}\\\ \varepsilon_{q+1}^{\prime}\\\
\vdots\\\ \varepsilon_{T}^{\prime}\end{pmatrix}_{(T-q+1)\times p},$ (14)

where
$\boldsymbol{\theta}^{(1)}=\boldsymbol{A}^{(1)},\boldsymbol{\theta}^{(2)}=\boldsymbol{A}^{(2)}-\boldsymbol{A}^{(1)}$.
With the larger dimension of the parameters, the same argument for the VAR
process of order $q$ can be achieved by following the logic from (5).

##  4 Theoretical results

In this section, we explore the asymptotic behaviour of the proposed method.
We show that our method controls the familywise error under the null (i.e.
when there exist no anomaly) with an appropriate threshold and give conditions
under which the asymptotic power of the method tends to 1. These results are
based upon the following assumptions.

###### Assumption 1

For each $j=1,2$, let $\Gamma_{j}(\ell)$ be the population version of the
lag-$\ell$ covariance matrix of $\boldsymbol{\mathfrak{x}}_{j}$ where
$\boldsymbol{\mathfrak{x}}_{j}$ is as in (3). For $\kappa\in[-\pi,\pi]$, there
exist the spectral density matrices,

$f_{j}(\kappa)=\frac{1}{2\pi}\sum_{l\in\mathbb{Z}}\Gamma_{j}(l)\exp^{-\sqrt{-1}\kappa
l}.$

In addition,

$\max_{j}\mathcal{M}(f_{j})=\max_{j}\Big{\\{}\text{ess}\sup_{\kappa\in[-\pi,\pi]}\Lambda_{\max}(f_{j}(\kappa))\Big{\\}}<+\infty,$

and

$\min_{j}\boldsymbol{\mathfrak{m}}(f_{j})=\min_{j}\Big{\\{}\text{ess}\inf_{\kappa\in[-\pi,\pi]}\Lambda_{\min}(f_{j}(\kappa))\Big{\\}}>0,$

where $\Lambda_{\max}(A)$ and $\Lambda_{\min}$ are the largest and the
smallest eigenvalues of the symmetric matrix $A$, respectively.

This first condition is needed to control the stability properties of the VAR
models. This is a spectral density condition that is not only valid for VAR
model but also holds for a large class of general linear process. Basu and
Michailidis, (2015) use the same assumption but for a stable VAR setting
without considering anomalies, while we extend it to the single collective
anomaly setting by assuming a spectral density function for each common and
anomalous segments separately.

In order to bound the power of our method we need conditions on the size and
length of any anomaly and the set of intervals we use – essentially we will
need at least one interval of sufficient length to be contained within the
anomaly. To this end we introduce the following:

###### Assumption 2

There exist at least one interval $J\in\mathbb{J}_{T,p}(L)$ such that
$J\subseteq[\eta_{1},\eta_{2}]$ and the choice of $L$ for a set of intervals
$\mathbb{J}_{T,p}(L)$ satisfies the following condition as
$T,p\rightarrow\infty$,

$\frac{log(T\lor p)}{L}\rightarrow 0$

where any interval $J\in\mathbb{J}_{T,p}(L)$ has length at least $L$.

###### Assumption 3

The sparsity of change is fixed; $\|\boldsymbol{\Theta}\|_{0}=d_{0}$.

###### Assumption 4

For any $\xi>0$, $L\cdot\|\boldsymbol{\Theta}\|_{2}^{2}>C_{2}\cdot
d_{0}^{2}\cdot\log^{1+\xi}{(T\lor p)}$, where $C_{2}>0$ is a constant.

Assumption 3 gives the condition on the number of nonzero entries of the
coefficient matrix, where the sparsity parameter $d_{0}$ affects the signal-
to-noise ratio condition in Assumption 4. Our Assumption 4 is similar to the
conditions required in other change-point problem in high-dimensional VAR
model. For example, Wang et al., (2019) study a multiple change point setting
and their signal-to-noise ratio assumption becomes equal to ours in the case
when single change-point is considered, while Safikhani and Shojaie, (2020)
assume $\|\boldsymbol{\Theta}\|_{2}$ is bounded away from zero.

Our final assumption is used to extend our results to the case where we
estimate $\boldsymbol{\theta}^{(1)}$.

###### Assumption 5

It holds for the estimator $\hat{\boldsymbol{\theta}}^{(1)}$ that
$\big{\|}\boldsymbol{\theta}^{(1)}-\hat{\boldsymbol{\theta}}^{(1)}\big{\|}_{\infty}<C\sqrt{\frac{\log({T\lor
p})}{L}}$ with probability approaching $1$ as $T\rightarrow\infty$ and
$p\rightarrow\infty$, where $C>0$ is a constant.

Assumption 5 states the condition on the estimation error bound in
$\ell_{\infty}$-norm. This is in line with the estimation error presented in
Proposition 4.1 of Basu and Michailidis, (2015) and Lemma 15 of Wang et al.,
(2019) in which the sparsity assumption is imposed on VAR coefficient
matrices. For instance, when $\boldsymbol{\theta}^{(1)}$ is assumed to be
sparse with the condition $\|\boldsymbol{\theta}^{(1)}\|_{0}=k$, then its
lasso estimator, $\hat{\boldsymbol{\theta}}^{(1)}$, satisfies
$\big{\|}\boldsymbol{\theta}^{(1)}-\hat{\boldsymbol{\theta}}^{(1)}\big{\|}_{2}\leq
c\sqrt{k}\sqrt{\frac{\log({T\lor p})}{T}}$, where
$\hat{\boldsymbol{\theta}}^{(1)}$ is obtained from a sample of size $T$. When
the sparsity $k$ is fixed, the estimation error bound in $\ell_{2}$-norm
implies Assumption 5.

We now present our main theoretical results. The following theorem gives
conditions on the lasso penalty to ensure the procedure asymptotically
controls the familywise error when there is no anomaly.

###### Theorem 1

Let Assumptions 1-2 hold. If there exist no anomaly, for a tuning parameter
$\lambda=C_{3}\sqrt{L(2\log{p}+\log{T})}$ with a constant $C_{3}$ large
enough, we have

$\displaystyle
P\bigg{(}\max_{J\in\mathbb{J}_{T,p}(L)}T^{\text{lasso}}(J)\leq\lambda^{\text{thr}}\bigg{)}$
$\displaystyle\geq
P\bigg{(}\max_{J\in\mathbb{J}_{T,p}(L)}T^{\text{lasso}}(J)=0\bigg{)}$
$\displaystyle\geq 1-C_{4}\exp(-C_{5}(2\log{p}+\log{T})),$

where $C_{4},C_{5}>0$, $\lambda^{\text{thr}}$ is strictly positive and
$\lambda$ is a tuning parameter controlling the penalty term in lasso
regression in (11).

In Theorem 1, it is clear that our result applies to any positive threshold
$\lambda^{\text{thr}}$. In the proof of Theorem 1, we show that the familywise
error is controlled under an appropriate tuning parameter $\lambda$ and the
proof can be found in the supplementary material. We now turn to the
asymptotics of the test statistic under the alternative.

###### Theorem 2

Let Assumptions 1-4 hold. If there exist an anomaly, with a tuning parameter
$\lambda=C_{2}\sqrt{L(2\log{p}+\log{T})}$ for a large enough $C_{2}$, as
$T\rightarrow\infty$,

$\displaystyle
P\bigg{(}\max_{J\in\mathbb{J}_{T,p}(L)}T^{\text{lasso}}(J)\leq\lambda^{\text{thr}}\bigg{)}\rightarrow
0$

and

$P(\hat{I}\cap[\eta_{1},\eta_{2}]\neq\emptyset)\rightarrow 1,$

where the threshold $\lambda^{\text{thr}}$ has the order of
$\sqrt{L\cdot\log(p\lor T)}$, the estimated anomaly $\hat{I}$ is as in (12)
and $\lambda$ is a tuning parameter controlling the penalty term in lasso
regression in (11).

Theorem 2 states that the test statistic corresponding to the intervals in the
candidate set is greater than the pre-specified threshold if the interval is
located within the true anomaly. In other words, it shows that the individual
test has asymptotic power one. The following theorem shows that our method has
larger power to detect a sparse collective anomaly.

###### Theorem 3

Assume that $\boldsymbol{x}_{t}$ follows (6) and let Assumptions 1-4 hold. Let
the null hypothesis hold, then for any
$\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\}$, the test
statistic of the OLS method in (10) follows a $\chi^{2}_{p^{2}}$ distribution.
Consequently, we have an asymptotic level $\alpha$ test if the null hypothesis
is rejected for $T(J)>\chi^{2}_{p^{2};(1-\alpha)}$, where
$\chi^{2}_{p^{2};(1-\alpha)}$ is the $(1-\alpha)$-quantile of chi-square
distribution with $p^{2}$ degrees of freedom.

Under the alternative, for any $J\in\mathbb{J}_{T,p}(L)$ such that
$J\subseteq[\eta_{1},\eta_{2}]$, the upper bound on the power of the OLS
method is given as

$\frac{E\big{(}\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}\big{)}}{W_{p}},$
(15)

where $W_{p}=O_{p}(p)$.

Note that $W_{p}$ in (15) is linked to the false positive rate as it is the
approximation of $\chi^{2}_{p^{2};(1-\alpha)}-p^{2}$. See the proof in the
supplementary material for further details.

Theorem 3 shows the asymptotic behaviours of the test statistic of the OLS
method under both the null and the alternative hypotheses. Furthermore,
Theorem 3 implies that the test statistic built on the lasso estimator can
detect weaker anomalies than using the OLS estimator when the change is
sparse. The intuition behind this is that the test statistic of the OLS method
in (10) can be written as

$\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}+\big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2}\big{\\}},$
(16)

and
$E(\|\boldsymbol{Y}\|_{2}^{2}-\|\boldsymbol{Y}-\boldsymbol{X}{\boldsymbol{\Theta}}\|_{2}^{2})$
needs to be at least as large as $O_{p}(p)$ to have high power. By comparison,
if we denote the lasso estimator of $\boldsymbol{\Theta}$ by
$\hat{\boldsymbol{\Theta}}$, then the test statistic of the lasso method in
(11) can be written as

$\displaystyle\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}-\lambda\|\boldsymbol{\Theta}\|_{1}+\big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}+\lambda\|\boldsymbol{\Theta}\|_{1}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2}-\lambda\|\hat{\boldsymbol{\Theta}}\|_{1}\big{\\}}.$
(17)

Noting that the term in $\\{\\}$s in (17) is positive, the lasso-based test
statistic requires that
$\|\boldsymbol{Y}\|_{2}^{2}-\|\boldsymbol{Y}-\boldsymbol{X}{\boldsymbol{\Theta}}\|_{2}^{2}$
should at least as large as $O_{p}(\lambda\|\boldsymbol{\Theta}\|_{1})$ and
$\lambda=C_{2}\sqrt{L(2\log{p}+\log{T})}$.

The following two corollaries state that the assertions in Theorems 1-2 remain
true if the $\boldsymbol{\theta}^{(1)}$ is replaced by an estimator
${\hat{\boldsymbol{\theta}}^{(1)}}$ that satisfies the condition in Assumption
5.

###### Corollary 1

Theorem 1 holds with a different constant if
${\hat{\boldsymbol{\theta}}^{(1)}}$ is used in calculating the test statistic
instead of the true parameter $\boldsymbol{\theta}^{(1)}$, where
${\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}$ is an estimator fulfilling
Assumption 5.

###### Corollary 2

Theorem 2 holds with a different constant if
${\hat{\boldsymbol{\theta}}^{(1)}}$ is used in calculating the test statistic
instead of the true parameter $\boldsymbol{\theta}^{(1)}$, where
${\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}$ is an estimator fulfilling
Assumption 5.

The proofs of Theorems 1-3 and Corollaries 1-2 can be found in the
supplementary material.

##  5 Simulation study

###  5.1 Preliminaries

We compare the performance of our lasso-based approach with the OLS method
described in Section 3.2. Whilst there are other methods for detecting changes
in a VAR model, such as those of Safikhani and Shojaie, (2020) and Bai et al.,
(2020), they are not designed for the collective anomaly setting that we
consider. For completeness, we compare their performances with ours, and the
details can be found in the supplementary material. Perhaps due to not being
designed for the collective anomaly setting, we find these alternative methods
perform substantially worse than ours, particularly when the underlying matrix
$A^{(1)}$ is dense but the change is sparse.

In practice, the underlying parameter $A^{(1)}$ is often unknown and needs to
be estimated. In this case, as the accuracy of our method depends on how
accurately we can estimate $A^{(1)}$, considering two extreme cases gives
upper and lower bounds on our method: $A^{(1)}$ is known and $A^{(1)}$ is
estimated from a relatively small amount of data. The threshold of each test
is selected by choosing the $99\%$ quantile of the test statistics obtained
through the 100 simulation runs performed under the null. For the error
variance, we set $\Sigma_{\varepsilon}$ to be the identity matrix. In the
following sections, we report the results when $\Sigma_{\varepsilon}$ is
known. The results for the case when $\Sigma_{\varepsilon}$ is estimated can
be found in the supplementary material.

We also look at how the choice of the set of intervals, $\mathbb{J}_{T,p}(L)$,
affects performance. We vary both the number of intervals which we denote by
$s$, and the way we choose the intervals, randomly or deterministically, with
a pre-determined minimum length of interval. For the deterministic
construction of intervals, we use the technique proposed in Definition 1 of
Kovács et al., (2020) with the decay parameter $1/a=1.1,1.2$. Regardless of
the way of choosing the intervals, we force the minimum length intervals to be
greater than $p$ in order to compare our approach with the OLS method. In the
following sections, we present the simulation results for two scenarios: (1)
$A^{(1)}$ is dense and (2) $A^{(1)}$ is sparse; where the number of non-zero
elements is large in $(1)$ and small in $(2)$. Note that when $A^{(1)}$ is
assumed to be unknown, it is estimated from the null region with ridge or
lasso penalty depending on the given sparsity of $A^{(1)}$.

###  5.2 Dense $A^{(1)}$

We first consider the case when all entries of $A^{(1)}$ are non-zero. The
coefficient matrix is randomly generated by using the algorithm proposed by
Ansley and Kohn, (1986) and implemented in R package gmvarkit which forces the
resulting VAR model to be stationary, where the range of the entries of
$A^{(1)}$ is obtained as $[-0.67,0.58]$. We set $T=500$ and $p=10$, where only
a few (five or ten) entries in the VAR coefficient matrix undergo change in
anomalous interval. We investigate both cases: (1) $A^{(1)}$ is assumed to be
known and (2) $A^{(1)}$ is estimated from the training data with a ridge
penalty. In the latter case, the training data contains the same amount of the
test data which we examine for detecting an anomaly. Our lasso-based method is
implemented by using the tuning parameter $\lambda$ presented in Theorems 1-2
with the constant $C=0.15$. In the following sections, we consider the single
anomaly and the multiple anomaly case.

####  5.2.1 Single anomaly

We consider a single anomaly interval located in the middle with three
different lengths, where the details can be found in Table 1 and the
coefficient matrices are presented in Figure 1. The non-zero entries in
$A^{(2)}-A^{(1)}$ are all equal to $\Delta=0.35$ and $A^{(2)}$ is made by
adding $\Delta$ to the first ten smallest positive entries of $A^{(1)}$.

Figure 1: The underlying coefficient matrices, $(A^{(1)},A^{(2)},A^{(1)})$, for the simulation setting in Section 5.2.1, where $A^{(2)}$ corresponds to an anomaly. | T | p | $[\eta_{1},\eta_{2}]$ | $\eta_{2}-\eta_{1}$ | $\Delta$ | $\|\boldsymbol{\Theta}\|_{0}$
---|---|---|---|---|---|---
case 1 | 500 | 10 | $[T(5/11),T(6/11)]$ | $45$ | 0.35 | 10
case 2 | 500 | 10 | $[T(7/15),T(8/15)]$ | $33$ | 0.35 | 10
Table 1: Simulation settings for two cases considered in Section 5.2.1, where $\Delta$ is the size of non-zero entries of $\boldsymbol{\Theta}$ and $\|\boldsymbol{\Theta}\|_{0}$ is the number of non-zero elements of $\boldsymbol{\Theta}$. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 100 | 82
(s = 1029) | Lasso | 100 | 94
deterministic | OLS | 100 | 85
(s = 1029) | Lasso | 100 | 95
deterministic | OLS | 100 | 85
(s = 540) | Lasso | 100 | 94
case 2 | random | OLS | 98 | 46
(s = 1029) | Lasso | 100 | 66
deterministic | OLS | 98 | 56
(s = 1029) | Lasso | 99 | 75
deterministic | OLS | 98 | 54
(s = 540) | Lasso | 99 | 72
Table 2: Empirical power ($\%$) from 100 simulation runs for two methods in all cases described in Section 5.2.1, where $s$ is the number of intervals examined. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 0.45 (0.31) | 10.50 (16.66)
(s = 1029) | Lasso | 0.40 (0.00) | 5.24 (10.34)
deterministic | OLS | 0.39 (0.25) | 7.18 (16.31)
(s = 1029) | Lasso | 0.35 (0.16) | 2.63 (9.96)
deterministic | OLS | 0.39 (0.33) | 7.15 (16.32)
(s = 540) | Lasso | 0.32 (0.22) | 3.06 (10.85)
case 2 | random | OLS | 3.42 (6.49) | 28.38 (21.04)
(s = 1029) | Lasso | 2.20 (0.91) | 19.08 (20.54)
deterministic | OLS | 1.32 (6.57) | 21.31 (23.20)
(s = 1029) | Lasso | 0.82 (4.67) | 12.46 (20.57)
deterministic | OLS | 1.27 (6.58) | 22.17 (23.34)
(s = 540) | Lasso | 0.76 (4.68) | 13.76 (21.34)
Table 3: The mean (standard deviation) of Hausdorff distance from 100
simulation runs for two methods in all cases described in Section 5.2.1, where
$s$ is the number of intervals examined.

As shown in Table 2, the lasso-based method tends to detect an anomaly more
often than the OLS-based approach in all cases regardless of the way of
choosing intervals to investigate and whether $A^{(1)}$ is known or unknown.
As expected, compared to the results when the true $A^{(1)}$ is known, both
OLS and lasso methods perform less well when $\hat{A}^{(1)}$ is used.
Comparing the randomly and the deterministically chosen segments with the size
($s$) equal to $1029$, for both the OLS and the lasso methods, the
deterministic way tends to give a slightly lower power when $A^{(1)}$ is known
but gives a similar or a slightly larger power when $A^{(1)}$ is estimated.
Note, when $A^{(1)}$ is estimated, the deterministically chosen intervals with
smaller sample size ($s=540$) shows a similar or a larger power than those
chosen randomly with a sample size ($s=1029$) for both methods, and the
difference becomes larger as the length of anomalous interval becomes shorter
(from case 1 to case 2 as presented in Table 1). Table 3 shows that the lasso
method also outperforms in terms of distance between the estimated and the
true anomaly and its variance.

####  5.2.2 Two collective anomalies

We now consider two collective anomalies, $[\eta_{1},\eta_{2}]$ and
$[\eta_{3},\eta_{4}]$, where the corresponding coefficient matrix is $A^{(2)}$
and $A^{(3)}$, respectively. $A^{(2)}$ and $A^{(3)}$ are obtained by adding
$\Delta_{1}$ and $\Delta_{2}$, respectively to the first five smallest
positive entries of $A^{(1)}$ and the true coefficient matrices are presented
in Figure 2. Two different cases are considered and the details are provided
in Table 4.

Figure 2: The underlying coefficient matrices, $(A^{(1)},A^{(2)},A^{(1)},A^{(3)},A^{(1)})$, for the simulation setting in Section 5.2.2, where $A^{(2)}$ and $A^{(3)}$ correspond to the anomalies. | T | p | $[\eta_{1},\eta_{2}]$ | $[\eta_{3},\eta_{4}]$ | $\eta_{2}-\eta_{1}$ | $\eta_{4}-\eta_{3}$ | $\Delta_{1}$ | $\Delta_{2}$ | $\|\boldsymbol{\Theta}\|_{0}$
---|---|---|---|---|---|---|---|---|---
case 1 | 500 | 10 | $[133,166]$ | $[333,366]$ | $33$ | $33$ | $0.6$ | $0.6$ | $5$
case 2 | 500 | 10 | $[33,66]$ | $[433,466]$ | $33$ | $33$ | $0.5$ | $0.5$ | $5$
Table 4: Simulation settings for two cases considered in Section 5.2.2, where $\Delta_{1}=|A^{(2)}-A^{(1)}|$, $\Delta_{2}=|A^{(3)}-A^{(1)}|$ and $\|\boldsymbol{\Theta}\|_{0}$ is the number of non-zero elements of $\boldsymbol{\Theta}$. |  |  |  $A^{(1)}$ is known |  |  $A^{(1)}$ is estimated
---|---|---|---|---|---
|  |  | 0 | 1 | 2 | 3 |  | 0 | 1 | 2 | 3
case 1 | random | OLS | 0 | 27 | 73 | 0 |  | 3 | 82 | 15 | 0
(s = 1944) | Lasso | 0 | 24 | 76 | 0 |  | 0 | 48 | 52 | 0
deterministic | OLS | 0 | 24 | 76 | 0 |  | 1 | 70 | 29 | 0
(s = 1944) | Lasso | 0 | 12 | 86 | 2 |  | 0 | 35 | 65 | 0
deterministic | OLS | 0 | 26 | 74 | 0 |  | 2 | 72 | 26 | 0
(s = 1029) | Lasso | 0 | 21 | 77 | 2 |  | 0 | 40 | 60 | 0
case 2 | random | OLS | 0 | 4 | 94 | 2 |  | 36 | 58 | 6 | 0
(s = 1944) | Lasso | 0 | 2 | 98 | 0 |  | 7 | 50 | 43 | 0
deterministic | OLS | 0 | 4 | 95 | 1 |  | 31 | 60 | 9 | 0
(s = 1944) | Lasso | 0 | 1 | 96 | 3 |  | 4 | 43 | 53 | 0
deterministic | OLS | 0 | 4 | 95 | 1 |  | 34 | 58 | 8 | 0
(s = 1029) | Lasso | 0 | 1 | 98 | 1 |  | 8 | 40 | 52 | 0
Table 5: Distribution of the number of detected anomalies for two methods in all cases described in Section 5.2.2 over 100 simulation runs, where $s$ is the number of intervals examined. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 3.01 (2.89) | 15.90 (12.21)
(s = 1944) | Lasso | 2.82 (2.83) | 7.38 (8.56)
deterministic | OLS | 2.46 (2.62) | 12.45 (13.04)
(s = 1944) | Lasso | 1.97 (2.68) | 4.41 (7.26)
deterministic | OLS | 2.59 (2.57) | 13.07 (13.06)
(s = 1029) | Lasso | 2.43 (2.66) | 5.11 (7.77)
case 2 | random | OLS | 3.93 (5.11) | 12.26 (2.83)
(s = 1944) | Lasso | 2.50 (1.88) | 8.38 (5.44)
deterministic | OLS | 1.92 (4.19) | 11.92 (3.76)
(s = 1944) | Lasso | 1.50 (3.80) | 6.66 (6.22)
deterministic | OLS | 1.98 (4.22) | 12.02 (3.58)
(s = 1029) | Lasso | 1.44 (3.68) | 6.79 (6.17)
Table 6: The mean (standard deviation) of Hausdorff distance from 100
simulation runs for two methods in all cases described in Section 5.2.2, where
$s$ is the number of intervals examined.

To detect multiple anomalies, the procedure presented in Algorithm 2 is
applied. From Table 5, we obtain similar interpretations to those from the
single anomaly case in Section 5.2.1. When the true $A^{(1)}$ is known, both
the OLS and the lasso methods give better results than when $A^{(1)}$ is
assumed to be unknown and estimated. For both methods and for both cases
($A^{(1)}$ is known and estimated), the deterministic settings ($s=1029$ and
$s=1944$) tend to return better results than the random setting with the
sample size of $s=1944$. In Table 6, the lasso method returns smaller mean and
standard deviation of the Hausdorff distance, regardless of the way of
choosing segments and whether $A^{(1)}$ is known or not.

###  5.3 Sparse $A^{(1)}$

| T | p | $[\eta_{1},\eta_{2}]$ | $\eta_{2}-\eta_{1}$ | $\|A^{(1)}\|_{\infty}$ | $\|A^{(2)}\|_{\infty}$ | $\|\boldsymbol{\Theta}\|_{0}$
---|---|---|---|---|---|---|---
case 1 | 500 | 20 | $[T(4/9),T(5/9)]$ | $55$ | $0.6$ | $0.05$ | $19$
case 2 | 500 | 20 | $[T(6/13),T(7/13)]$ | $39$ | $0.6$ | $0.05$ | $19$
Table 7: Simulation setting for two cases considered in Section 5.3, where $\|A^{(1)}\|_{\infty}$ and $\|A^{(2)}\|_{\infty}$ are the size of non-zero elements in $A^{(1)}$ and $A^{(2)}$, respectively and $\|\boldsymbol{\Theta}\|_{0}$ is the number of non-zero elements of $\boldsymbol{\Theta}$. Figure 3: The underlying coefficient matrices, $(A^{(1)},A^{(2)},A^{(1)})$, for the simulation setting in Section 5.3, where $A^{(2)}$ corresponds to an anomaly. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 100 | 90
(s = 499) | Lasso | 100 | 100
deterministic | OLS | 100 | 95
(s = 499) | Lasso | 100 | 100
case 2 | random | OLS | 100 | 30
(s = 499) | Lasso | 100 | 58
deterministic | OLS | 100 | 50
(s = 499) | Lasso | 100 | 92
Table 8: Empirical power ($\%$) from 100 simulation runs for two methods in all cases described in Section 5.3, where $s$ is the number of intervals examined. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 1.43 (0.53) | 6.00 (13.03)
(s = 499) | Lasso | 1.31 (0.45) | 1.28 (0.21)
deterministic | OLS | 0.43 (0.24) | 2.70 (9.71)
(s = 499) | Lasso | 0.34 (0.12) | 0.37 (0.13)
case 2 | random | OLS | 3.26 (1.19) | 32.80 (20.89)
(s = 499) | Lasso | 3.21 (1.18) | 20.12 (22.48)
deterministic | OLS | 0.43 (0.56) | 23.39 (23.12)
(s = 499) | Lasso | 0.30 (0.19) | 4.03 (12.56)
Table 9: The mean (standard deviation) of Hausdorff distance from 100
simulation runs for two methods in all cases described in Section 5.3, where
$s$ is the number of intervals examined.

We now consider the case when $A^{(1)}$ is sparse i.e. only a smaller number
of entries are non-zero. Similar to the settings used in Safikhani and
Shojaie, (2020), the 1-off diagonal values of the coefficient matrix are non-
zero as shown in Figure 14. The details of the simulation setting are given in
Table 7. Tables 8 and 9 show similar interpretations with those given in
Section 5.2.

##  6 Data analysis

###  6.1 Yellow cab demand in New York City

Figure 4: (Top) The differenced yellow taxi pickups recorded from March 11,
2019 to March 6, 2020 in Manhattan. (Middle) The 20 largest test statistics
with the corresponding interval. The blue horizontal dashed line indicates the
threshold. (Bottom) The portion of the top plot indicated with dashed green
vertical lines. Red vertical lines show the estimated anomaly, [Nov $3,2019$
$00:00:00$, Nov $3,2019$ $06:30:00$].

To demonstrate the usefulness of our method, we now turn to real data
applications. In our first example, we apply our method to the yellow taxi
trip data that is previously analysed by Safikhani and Shojaie, (2020). The
data can be downloaded from the New York City Taxi and Limousine Commission
(TLC) Database (https://www1.nyc.gov/site/tlc/about/tlc-trip-record-
data.page). This data consists of the number of yellow taxi pick-ups recorded
from 10 randomly selected zones in Manhattan, a borough in New York City. We
aggregate the number of yellow taxi pick-ups every 30 minutes from March 11,
2019 to March 6, 2020 which results in 17376 time points. We seek to detect
whether a collective anomaly exists after removing the first order
nonstationarity from the data. We use the differenced version of the time
series, using the first 4344 data points to estimate the underlying VAR
coefficient $A^{(1)}$ by applying a lasso penalty. The next 4344 data points
are used to obtain a threshold, where the threshold is selected as the 99$\%$
quantile of the test statistics from 100 deterministically chosen intervals.
Then we detect a single anomaly using the remaining 8687 data points.

Figure 5: Taxi demand (Top) and differenced Taxi demand (Bottom) for 2nd
(black), 4th (red) and 6th (green) zones in Manhattan recorded from October
24, 2019 to November 5, 2019.

The top plot in Figure 4 shows that two consecutive spikes are observed
between October 7 and November 6 in 2019, where the interval within green
vertical lines is enlarged in the bottom plot. From the middle plot, we see
that the largest test statistic is obtained for a small interval which
includes the spikes shown in the top plot. The bottom plot shows that the
spikes occur between 12am to 2am on November 3, 2019 and our method detects an
anomaly between 12am and 6:30am on November 3, 2019. From Figure 5, we see
that a sudden high demand occurred at the $4^{\text{th}}$ and $6^{\text{th}}$
zones located in Downtown Manhattan (also known as Lower Manhattan). This
anomaly seems to be related to traffic management for the 2019 New York City
Marathon which took placed on November 3, 2019 in New York City. We can
interpret that there was a sudden high demand in Downtown Manhattan where the
marathon route did not pass through, and this changes the relationship between
the 10 zones we investigate.

###  6.2 EEG Data

Figure 6: (Top) EEG data recorded at 18 different channels. Blue solid
vertical line is the time at which the neurologist thinks seizure starts and
the red dashed vertical line is the anomaly detected in the online setting.
(Bottom) The maximum test statistics at each time point obtained through
Algorithm 3 which stops when the anomaly is detected. The horizontal red line
presents the pre-specified threshold.

We now show how our method can be used in as an online changepoint detection
method. We demonstrate this on electroencephalogram (EEG) data collected from
an epileptic patient. Other ways of analysing this dataset can be found in
Ombao et al., (2001), Ombao et al., (2005) and Schröder and Ombao, (2019). The
data consists of brain electrical potentials recorded by placing electrodes on
18 locations on the scalp of a patient. The EEG signals are recorded during an
epileptic seizure, thus these exists a visible change in the data as shown in
Figure 6. The brain wave patterns are recorded over $500$ seconds with the
sampling rate 100 Hz (i.e. 100 points per second). As done in Safikhani and
Shojaie, (2020), to speed up computation, we use 2 observations per second
which reduces the number of time points to $T=1000$.

We separate the data into a training set of the size $T_{1}=600$ and a test
set of the size $T_{2}=400$. The first half of the training set is used to
estimate the underlying VAR coefficient $A^{(1)}$ by applying a lasso penalty
and the second half is used to have a threshold that is chosen as the 99$\%$
quantile of the test statistics computed from 327 deterministically chosen
intervals. Then we perform the single anomaly detection using a test set.

As mentioned in Section 2, here we show how our method can be applied to the
online framework. We refer the reader to Fisch et al., (2020) and Yu et al.,
(2021) for the recent works on online detection algorithm for change-points or
anomalies. In the online setting, we make sequential decisions about the
occurrence of an anomaly whenever each new observations is obtained. Our
algorithm for online anomaly detection is similar to Algorithm 2 of Yu et al.,
(2021). The detailed procedure is given in Algorithm 3 where we set
$t_{0}=10$. As shown in Figure 6, an anomaly is estimated at $t=119$ that has
the detection delay of $5$ time points compared to $t=114$ at which the
neurologist states that a seizure takes place.

INPUT: $\boldsymbol{X}$, $\lambda^{\text{thr}}$, $t_{0}$

$t\leftarrow t_{0}$

FLAG $\leftarrow 0$

while  _FLAG $=0$_ do

$t\leftarrow t+1$

$J\leftarrow\Big{\lfloor}\frac{\log t}{\log 2}\Big{\rfloor}$

$j\leftarrow 1$

while  _FLAG = 0 and $j\leq J$_ do

$s_{j}\leftarrow t-2^{j-1}$

$J\leftarrow[s_{j},t]$

FLAG
$\leftarrow\mathbb{1}\\{T^{\text{lasso}}(J)>\lambda^{\textsuperscript{thr}}\\}$

$j\leftarrow j+1$

end while

end while

OUTPUT : $t$.

Algorithm 3 Online anomaly detection

##  7 Discussion

Our lasso-based approach is motivated for data where we have substantially
more data about the current or normal behaviour of the time series than for
any anomaly or epidemic change. Thus it is natural to model the change as
sparse and thus a lasso-based test is more appropriate than a standard
likelihood-ratio or OLS-based test. We provide a numerical evidence that our
method outperforms existing competitors in detecting sparse change when
$A^{(1)}$ is either dense or sparse. Our method searches a set of local
segments to detect an anomalous interval, whereas the existing change
detection methodologies for the VAR model perform global optimisation. As
illustrated in real data example, the local optimisation aspect of our method
gives a flexibility to extend it to the online setting.

## Acknowledgements

The authors thank Hernando Ombao and Abolfazl Safikhani for kindly providing
access to the EEG data, and also thank Yi Yu for helpful conversations. The
authors gratefully acknowledge the financial support of EPSRC, grant
EP/N031938/1 (STATSCALE).

## References

  * Ansley and Kohn, (1986) Ansley, C. F. and Kohn, R. (1986).  A note on reparameterizing a vector autoregressive moving average model to enforce stationarity.  Journal of Statistical Computation and Simulation, 24:99–106.
  * Aue et al., (2009) Aue, A., Hörmann, S., Horváth, L., and Reimherr, M. (2009).  Break detection in the covariance structure of multivariate time series models.  The Annals of Statistics, 37:4046–4087.
  * Aue et al., (2006) Aue, A., Horváth, L., Hušková, M., and Kokoszka, P. (2006).  Change-point monitoring in linear models.  The Econometrics Journal, 9:373–403.
  * Bai et al., (2020) Bai, P., Safikhani, A., and Michailidis, G. (2020).  Multiple change points detection in low rank and sparse high dimensional vector autoregressive models.  IEEE Transactions on Signal Processing, 68:3074–3089.
  * Baltagi et al., (2017) Baltagi, B. H., Kao, C., and Wang, F. (2017).  Identification and estimation of a large factor model with structural instability.  Journal of Econometrics, 197:87–100.
  * Bańbura et al., (2010) Bańbura, M., Giannone, D., and Reichlin, L. (2010).  Large Bayesian vector auto regressions.  Journal of applied Econometrics, 25:71–92.
  * Baranowski et al., (2019) Baranowski, R., Chen, Y., and Fryzlewicz, P. (2019).  Narrowest-over-threshold detection of multiple change points and change-point-like features.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81:649–672.
  * Barigozzi et al., (2018) Barigozzi, M., Cho, H., and Fryzlewicz, P. (2018).  Simultaneous multiple change-point and factor analysis for high-dimensional time series.  Journal of Econometrics, 206:187–225.
  * Barigozzi and Hallin, (2017) Barigozzi, M. and Hallin, M. (2017).  A network analysis of the volatility of high dimensional financial series.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 66:581–605.
  * Basu et al., (2019) Basu, S., Li, X., and Michailidis, G. (2019).  Low rank and structured modeling of high-dimensional vector autoregressions.  IEEE Transactions on Signal Processing, 67:1207–1222.
  * Basu and Michailidis, (2015) Basu, S. and Michailidis, G. (2015).  Regularized estimation in sparse high-dimensional time series models.  The Annals of Statistics, 43:1535–1567.
  * Bernanke et al., (2005) Bernanke, B. S., Boivin, J., and Eliasz, P. (2005).  Measuring the effects of monetary policy: a factor-augmented vector autoregressive (favar) approach.  The Quarterly journal of economics, 120:387–422.
  * Breitung and Eickmeier, (2011) Breitung, J. and Eickmeier, S. (2011).  Testing for structural breaks in dynamic factor models.  Journal of Econometrics, 163:71–84.
  * Calomiris et al., (2008) Calomiris, C. W., Longhofer, S. D., and Miles, W. (2008).  The foreclosure-house price nexus: lessons from the 2007-2008 housing turmoil.  Technical report, National Bureau of Economic Research.
  * Chen and Gupta, (1997) Chen, J. and Gupta, A. K. (1997).  Testing and locating variance changepoints with application to stock prices.  Journal of the American Statistical association, 92:739–747.
  * Chen et al., (2014) Chen, L., Dolado, J. J., and Gonzalo, J. (2014).  Detecting big structural breaks in large factor models.  Journal of Econometrics, 180:30–48.
  * Chen and Tian, (2010) Chen, Z. and Tian, Z. (2010).  Modified procedures for change point monitoring in linear models.  Mathematics and computers in simulation, 81:62–75.
  * Cho, (2016) Cho, H. (2016).  Change-point detection in panel data via double cusum statistic.  Electronic Journal of Statistics, 10:2000–2038.
  * Cho and Fryzlewicz, (2015) Cho, H. and Fryzlewicz, P. (2015).  Multiple-change-point detection for high dimensional time series via sparsified binary segmentation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77:475–507.
  * Cribben and Yu, (2017) Cribben, I. and Yu, Y. (2017).  Estimating whole-brain dynamics by using spectral clustering.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 66:607–627.
  * Davis et al., (2016) Davis, R. A., Zang, P., and Zheng, T. (2016).  Sparse vector autoregressive modeling.  Journal of Computational and Graphical Statistics, 25:1077–1096.
  * De Mol et al., (2008) De Mol, C., Giannone, D., and Reichlin, L. (2008).  Forecasting using a large number of predictors: Is Bayesian shrinkage a valid alternative to principal components?  Journal of Econometrics, 146:318–328.
  * Dette and Gösmann, (2020) Dette, H. and Gösmann, J. (2020).  A likelihood ratio approach to sequential change point detection for a general class of parameters.  Journal of the American Statistical Association, 115:1361–1377.
  * Fisch et al., (2020) Fisch, A., Bardwell, L., and Eckley, I. A. (2020).  Real time anomaly detection and categorisation.  arXiv preprint arXiv:2009.06670.
  * Fisch et al., (2018) Fisch, A. T. M., Eckley, I. A., and Fearnhead, P. (2018).  A linear time method for the detection of point and collective anomalies.  arXiv preprint arXiv:1806.01947.
  * Fremdt, (2015) Fremdt, S. (2015).  Page’s sequential procedure for change-point detection in time series regression.  Statistics, 49:128–155.
  * Fryzlewicz, (2014) Fryzlewicz, P. (2014).  Wild binary segmentation for multiple change-point detection.  The Annals of Statistics, 42:2243–2281.
  * Horváth et al., (2004) Horváth, L., Hušková, M., Kokoszka, P., and Steinebach, J. (2004).  Monitoring changes in linear models.  Journal of Statistical Planning and Inference, 126:225–251.
  * Inclan and Tiao, (1994) Inclan, C. and Tiao, G. C. (1994).  Use of cumulative sums of squares for retrospective detection of changes of variance.  Journal of the American Statistical Association, 89:913–923.
  * Kim and Siegmund, (1989) Kim, H.-J. and Siegmund, D. (1989).  The likelihood ratio test for a change-point in simple linear regression.  Biometrika, 76:409–423.
  * Kirch et al., (2015) Kirch, C., Muhsal, B., and Ombao, H. (2015).  Detection of changes in multivariate time series with application to eeg data.  Journal of the American Statistical Association, 110:1197–1216.
  * Kovács et al., (2020) Kovács, S., Li, H., Bühlmann, P., and Munk, A. (2020).  Seeded binary segmentation: A general methodology for fast and optimal change point detection.  arXiv preprint arXiv:2002.06633.
  * Lai, (1995) Lai, T. L. (1995).  Sequential changepoint detection in quality control and dynamical systems.  Journal of the Royal Statistical Society: Series B (Methodological), 57(4):613–644.
  * Lin and Michailidis, (2017) Lin, J. and Michailidis, G. (2017).  Regularized estimation and testing for high-dimensional multi-block vector-autoregressive models.  The Journal of Machine Learning Research, 18:4188–4236.
  * Lütkepohl, (2005) Lütkepohl, H. (2005).  New introduction to multiple time series analysis.  Springer Science & Business Media.
  * Nicholson et al., (2020) Nicholson, W. B., Wilms, I., Bien, J., and Matteson, D. S. (2020).  High dimensional forecasting via interpretable vector autoregression.  Journal of Machine Learning Research, 21:1–52.
  * Ombao et al., (2005) Ombao, H., Von Sachs, R., and Guo, W. (2005).  Slex analysis of multivariate nonstationary time series.  Journal of the American Statistical Association, 100:519–531.
  * Ombao et al., (2001) Ombao, H. C., Raz, J. A., von Sachs, R., and Malow, B. A. (2001).  Automatic statistical analysis of bivariate nonstationary time series.  Journal of the American Statistical Association, 96:543–560.
  * Rapach et al., (2007) Rapach, D. E., Strauss, J. K., et al. (2007).  Forecasting real housing price growth in the eighth district states.  Federal Reserve Bank of St. Louis. Regional Economic Development, 3:33–42.
  * Safikhani and Shojaie, (2020) Safikhani, A. and Shojaie, A. (2020).  Joint structural break detection and parameter estimation in high-dimensional nonstationary var models.  Journal of the American Statistical Association, pages 1–14.
  * Schröder and Ombao, (2019) Schröder, A. L. and Ombao, H. (2019).  Fresped: Frequency-specific change-point detection in epileptic seizure multi-channel eeg data.  Journal of the American Statistical Association, 114:115–128.
  * Seth et al., (2015) Seth, A. K., Barrett, A. B., and Barnett, L. (2015).  Granger causality analysis in neuroscience and neuroimaging.  Journal of Neuroscience, 35:3293–3297.
  * Shojaie and Michailidis, (2010) Shojaie, A. and Michailidis, G. (2010).  Discovering graphical granger causality using the truncating lasso penalty.  Bioinformatics, 26:i517–i523.
  * Siegmund and Venkatraman, (1995) Siegmund, D. and Venkatraman, E. (1995).  Using the generalized likelihood ratio statistic for sequential detection of a change-point.  The Annals of Statistics, pages 255–271.
  * Sims, (1980) Sims, C. A. (1980).  Macroeconomics and reality.  Econometrica: journal of the Econometric Society, pages 1–48.
  * Smith, (2012) Smith, S. M. (2012).  The future of fmri connectivity.  Neuroimage, 62:1257–1266.
  * Song and Bickel, (2011) Song, S. and Bickel, P. J. (2011).  Large vector auto regressions.  arXiv preprint arXiv:1106.3915.
  * Stock and Watson, (2008) Stock, J. H. and Watson, M. (2008).  The evolution of national and regional factors in us housing construction.  Volatility and Time Series Econometrics.
  * Tveten et al., (2020) Tveten, M., Eckley, I. A., and Fearnhead, P. (2020).  Scalable changepoint and anomaly detection in cross-correlated data with an application to condition monitoring.  arXiv preprint arXiv:2010.06937.
  * Wang et al., (2017) Wang, D., Yu, Y., and Rinaldo, A. (2017).  Optimal covariance change point localization in high dimension.  arXiv preprint arXiv:1712.09912.
  * Wang et al., (2019) Wang, D., Yu, Y., Rinaldo, A., and Willett, R. (2019).  Localizing changes in high-dimensional vector autoregressive processes.  arXiv preprint arXiv:1909.06359.
  * Wang and Samworth, (2018) Wang, T. and Samworth, R. J. (2018).  High dimensional change point estimation via sparse projection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80:57–83.
  * Yao, (1993) Yao, Q. (1993).  Tests for change-points with epidemic alternatives.  Biometrika, 80(1):179–191.
  * Yau and Zhao, (2016) Yau, C. Y. and Zhao, Z. (2016).  Inference for multiple change points in time series via likelihood ratio scan statistics.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 895–916.
  * Yu et al., (2021) Yu, Y., Padilla, O. H. M., Wang, D., and Rinaldo, A. (2021).  Optimal network online change point localisation.  arXiv preprint arXiv:2101.05477.

Supplementary Material for Collective anomaly detection in High-dimensional
VAR Models

Hyeyoung Maeng, Idris Eckley and Paul Fearnhead
Lancaster University, United Kingdom

##  1 Technical proofs

We first give a preparatory lemma and then move onto the proofs of main
theorems and corollaries presented in Section 4.

###### Lemma 1

Let Assumptions 1-4 hold. For any $J\in\mathbb{J}_{T,p}(L)$ such that
$J\subseteq[\eta_{1},\eta_{2}]$ and any $|J|\succsim
c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{1}\log p$, with probability at
least $1-T^{-6}$, we have

${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\geq
c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}|J|\cdot\|{\boldsymbol{\Theta}}\|^{2}_{2}-c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}\log(p)\cdot\|{\boldsymbol{\Theta}}\|^{2}_{1}$

where
$c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{1},c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2},c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}>0$
are some constants depending on $\boldsymbol{\mathfrak{m}}$ and $\mathcal{M}$.

##### Proof of Lemma 1

The argument follows the proof of Lemma 13-(b) of Wang et al., (2019).

##### Proof of Theorem 1

By the construction of the candidate set $\mathbb{I}^{*}$, it is sufficient to
show that $T^{\text{lasso}}(J)\rightarrow 0$ under the null, where
$T^{\text{lasso}}(J)$ is as in (11). The KKT conditions for the lasso problem
in (11) is that any $\hat{\boldsymbol{\Theta}}$ is optimal if and only if
there exists a subgradient $\hat{s}$ such that

$\boldsymbol{X}_{J}^{\top}\Big{(}\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\Big{)}=\lambda\hat{s}_{J},$
(18)

where $\hat{s}_{J}=\partial\big{|}\hat{\boldsymbol{\Theta}}\big{|}_{1}$ is a
subgradient of the $l_{1}$ norm evaluated at $\hat{\boldsymbol{\Theta}}$ which
takes the form

$\hat{s}_{J}=\text{sgn}(\hat{\boldsymbol{\Theta}})\;\text{for}\;\hat{\boldsymbol{\Theta}}\neq
0,\quad|\hat{s}_{J}|\leq 1\;\text{otherwise}.$ (19)

As $\boldsymbol{Y}=\boldsymbol{E}$ under the null, (18) and (19) give a
condition on $\boldsymbol{X}$ and $\boldsymbol{E}$ to ensure that we estimate
$\boldsymbol{\Theta}=\boldsymbol{0}$ as follows: for any
$J\in\mathbb{J}_{T,p}(L)$ such that $J\cap[\eta_{1},\eta_{2}]=\emptyset$,

$\max_{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset}\Big{\|}\boldsymbol{X}_{J}^{\top}\boldsymbol{E}_{J}\Big{\|}_{\infty}\leq\lambda.$
(20)

We remind that $\mathcal{X}_{J}$ is the unvectorised covariates as follows:

$\mathcal{X}_{J}=\begin{pmatrix}0\\\ 0\\\ \vdots\\\ 0\\\
\boldsymbol{x}_{t+1}^{\prime}\\\ \vdots\\\
\boldsymbol{x}_{t+h-1}^{\prime}\end{pmatrix}_{(2h-1)\times p}$

and note that

$\displaystyle\max_{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset}\Bigg{\|}\frac{\boldsymbol{X}_{J}^{\top}\boldsymbol{E}_{J}}{|J|}\Bigg{\|}_{\infty}$
$\displaystyle=\max_{\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\},1\leq
i,j\leq
p}\Bigg{|}e^{\prime}_{i}\bigg{(}\frac{\mathcal{X}^{\prime}_{J}E_{J}}{|J|}\bigg{)}e_{j}\Bigg{|},$
$\displaystyle\leq\max_{\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\},1\leq
i,j\leq
p}\Bigg{|}e^{\prime}_{i}\bigg{(}\frac{\mathcal{X}^{\prime}_{J}E_{J}}{L}\bigg{)}e_{j}\Bigg{|},$
(21)

where $e_{i}\in\mathbb{R}^{p}$ with the $i$-th element equals to $1$ and zero
otherwise. Similar to the argument used in Proposition $2.4(b)$ of Basu and
Michailidis, (2015), for fixed $i,j,J$, there exist $k_{1},k_{2}>0$ such that
for all $\gamma>0$:

$P\Bigg{(}\bigg{|}e^{\prime}_{i}\bigg{(}\mathcal{X}^{\prime}_{J}E_{J}\bigg{)}e_{j}\bigg{|}>k_{1}L\gamma\Bigg{)}\leq
6\exp(-k_{2}L\min(\gamma,\gamma^{2})).$ (22)

As the number of intervals contained in $\mathbb{J}_{T,p}(L)$ is of the order
$O(T)$ when they are constructed through the seeded interval idea in Kovács et
al., (2020), we consider the union over $p^{2}\cdot T$ possible choices of
$i,j,J$ in (21). Then the result follows by setting
$\gamma=k_{3}\sqrt{\frac{2\log{p}+\log{T}}{L}}$ for a large enough $k_{3}>0$.
Therefore, with probability at least $1-C_{4}\exp(-C_{5}(2\log{p}+\log{T}))$,
we have

$\max_{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset}\Big{\|}\boldsymbol{X}_{J}^{\top}\boldsymbol{E}_{J}\Big{\|}_{\infty}\leq
C_{3}\sqrt{L\log({T\lor p})},$ (23)

where $C_{4}>0$ and $C_{5}>0$. Having the condition
$\lambda=C_{3}\sqrt{L\log({T\lor p})}$ with a large enough $C_{3}>0$ in (20),
we obtain $\hat{\boldsymbol{\Theta}}=\boldsymbol{0}$ with probability at least
$1-C_{4}\exp(-C_{5}(2\log{p}+\log{T}))$. Therefore, under the null, the
probability that $T^{\text{lasso}}(J)\rightarrow 0$ for
$J\in\mathbb{J}_{T,p}(L)$ such that $J\cap[\eta_{1},\eta_{2}]=\emptyset$ is at
least $1-C_{4}\exp(-C_{5}(2\log{p}+\log{T}))$ where $C_{4},C_{5}>0$.

We emphasise that (23) can be applied to any serially uncorrelated Gaussian
errors
$\boldsymbol{\varepsilon}_{t}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}N(\boldsymbol{0},\Sigma_{\varepsilon})$
as the constant $k_{1}$ in (22) presented in Proposition $2.4(b)$ of Basu and
Michailidis, (2015) has a form of

$k_{1}=2\pi\Lambda_{\max}(\Sigma_{\varepsilon})\bigg{(}1+\frac{1+\mu_{\max}(\mathcal{A})}{\mu_{\min}(\mathcal{A})}\bigg{)},$

where $\Lambda_{\max}(\Sigma_{\varepsilon})$ is the maximum eigenvalue of
$\Sigma_{\varepsilon}$,
$\mu_{\max}(\mathcal{A})=\max_{|z|=1}\Lambda_{\max}(\mathcal{A}^{*}(z)\mathcal{A}(z))$,
$\mu_{\min}(\mathcal{A})=\min_{|z|=1}\Lambda_{\min}(\mathcal{A}^{*}(z)\mathcal{A}(z))$
and $\mathcal{A}(z)=\mathit{I}_{p}-\mathbf{A}^{(1)}z$ for the VAR(1) model and
$\mathcal{A}(z)=\mathit{I}_{p}-\sum_{d=1}^{q}\mathbf{A}^{(1)}_{d}z^{d}$ for
the VAR(q) model. Therefore, even if $\Sigma_{\varepsilon}$ is not an identity
matrix, we can have (23) with a different constant $C_{3}$ which depends on
the maximum eigenvalue of $\Sigma_{\varepsilon}$.

##### Proof of Theorem 2

It is sufficient to prove that for any $J\in\mathbb{J}_{T,p}(L)$ such that
$J\subseteq[\eta_{1},\eta_{2}]$, with probability approaching to $1$ as
$T\rightarrow\infty$,

$\big{\|}\boldsymbol{Y}_{J}\big{\|}_{2}^{2}-\big{\|}\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\big{\|}_{2}^{2}>\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}.$
(24)

This is because the other part in equation (17),

$\Big{\\{}\big{\|}\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\big{\|}_{2}^{2}+\lambda\|\boldsymbol{\Theta}\|_{1}-\big{\|}\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\big{\|}_{2}^{2}-\lambda\|\hat{\boldsymbol{\Theta}}\|_{1}\Big{\\}},$
(25)

is always positive and the left-hand side of (24) dominates (25). We can
simplify (24) as

${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{E}_{J}>\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}.$
(26)

The left-hand side of (26) is a Gaussian variable that can be written as
$\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}+2\boldsymbol{\nu}_{J}^{\top}\boldsymbol{E}_{J}\sim
N(\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J},4\boldsymbol{\nu}_{J}^{\top}\Sigma_{\varepsilon}\boldsymbol{\nu}_{J})$,
where $\boldsymbol{\nu}_{J}=\boldsymbol{X}_{J}{\boldsymbol{\Theta}}$ and
$\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}\rightarrow\infty$. Then for
any $g(J)=o(\sqrt{\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}})$ that goes
to $\infty$, we have the following bound with probability approaching to $1$,

$\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}+2\boldsymbol{\nu}_{J}^{\top}\boldsymbol{E}_{J}\geq\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}-g(J)\sqrt{4\gamma\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}},$
(27)

where $\gamma$ is the maximum eigenvalue of $\Sigma_{\varepsilon}$. The right-
hand side of (27) is of order
$\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}$, thus we now show that
${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}$
is bounded by
$\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}$ with
probability tending to $1$. From Lemma 1, with probability approaching to $1$,
we have

$\displaystyle{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}$
$\displaystyle\geq
c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}|J|\cdot\|{\boldsymbol{\Theta}}\|^{2}_{2}-c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}\log(p)\cdot\|{\boldsymbol{\Theta}}\|^{2}_{1}$
$\displaystyle\geq
c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}L\cdot\|{\boldsymbol{\Theta}}\|^{2}_{2}-c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}\log(p)\cdot\|{\boldsymbol{\Theta}}\|^{2}_{1},$

where
$c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2},c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}>0$,
thus we now show

$c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}\|{\boldsymbol{\Theta}}\|^{2}_{2}>c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}\frac{\log(p)}{L}\cdot\|{\boldsymbol{\Theta}}\|^{2}_{1}+\frac{\lambda}{L}\|\boldsymbol{\Theta}\|_{1}+\frac{\lambda^{\textsuperscript{thr}}}{L},$
(28)

as $T,p\rightarrow\infty$. We can obtain (28) as $T,p\rightarrow\infty$ from
combining

$\displaystyle(a)\;c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}\|{\boldsymbol{\Theta}}\|^{2}_{2}>c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{3}\frac{\log(p)}{L}\cdot\|{\boldsymbol{\Theta}}\|^{2}_{1},$
$\displaystyle(b)\;c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}\|{\boldsymbol{\Theta}}\|^{2}_{2}>\frac{\lambda}{L}\|\boldsymbol{\Theta}\|_{1},$
$\displaystyle(c)\;c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}\|{\boldsymbol{\Theta}}\|^{2}_{2}>\frac{\lambda^{\textsuperscript{thr}}}{L},$

where (a) can be shown by using
$d_{0}\|\boldsymbol{\Theta}\|_{2}^{2}\geq\|\boldsymbol{\Theta}\|_{1}^{2}$ from
Assumption 3 and $\frac{\log p}{L}\rightarrow 0$ from Assumption 2. By using
$d_{0}\|\boldsymbol{\Theta}\|_{2}^{2}\geq\|\boldsymbol{\Theta}\|_{1}^{2}$, (b)
becomes
$c^{\boldsymbol{\mathfrak{m}},\mathcal{M}}_{2}\|{\boldsymbol{\Theta}}\|_{2}>\frac{\lambda}{L}\sqrt{d_{0}}$
that can be achieved from ${\frac{\lambda}{L}}=\sqrt{\frac{C_{3}\log(T\lor
p)}{L}}$ and
$\|\boldsymbol{\Theta}\|_{2}^{2}>C_{2}\cdot\frac{\log^{1+\xi}{(T\lor p)}}{L}$
in Assumption 4. Similarly (c) can be obtained from
$\frac{\lambda^{\textsuperscript{thr}}}{L}=O\Bigg{(}\sqrt{\frac{\log(T\lor
p)}{L}}\Bigg{)}$ and Assumption 4.

We now consider the case $\Sigma_{\varepsilon}$ is not an identity matrix. In
that case, (24) becomes

$\boldsymbol{Y}_{J}^{\top}\Sigma_{\varepsilon}^{-1}\boldsymbol{Y}_{J}-(\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}})^{\top}\Sigma_{\varepsilon}^{-1}(\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}})>\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}},$

thus (26) becomes

${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\Sigma_{\varepsilon}^{-1}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\Sigma_{\varepsilon}^{-1}\boldsymbol{E}_{J}>\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}},$
(29)

which holds by following the same argument used above with
$\boldsymbol{\nu}_{J}=\Sigma_{\varepsilon}^{-1/2}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}$
and different constants, as the left-hand side of (29) is a Gaussian random
variable bounded by a component that is of order
$\boldsymbol{\nu}_{J}^{\top}\boldsymbol{\nu}_{J}$.

Lastly, without repeating all the proofs, we argue that the theory we present
for the known $\Sigma_{\varepsilon}$ can be applied to the case when an
estimate of ${\Sigma_{\varepsilon}}$ is used. If $\hat{\Sigma}_{\varepsilon}$
is used instead of ${\Sigma_{\varepsilon}}$, the left-hand side of (29) can be
rewritten as

${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\Sigma_{\varepsilon}^{-1}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\Sigma_{\varepsilon}^{-1}\boldsymbol{E}_{J}+{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}(\hat{\Sigma_{\varepsilon}}^{-1}-\Sigma_{\varepsilon}^{-1})\boldsymbol{X}_{J}{\boldsymbol{\Theta}}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}(\hat{\Sigma_{\varepsilon}}^{-1}-\Sigma_{\varepsilon}^{-1})\boldsymbol{E}_{J},$
(30)

thus the test depends on the eigenvalues of the measure of the distance
between $\hat{\Sigma}^{-1}_{\varepsilon}$ and $\Sigma^{-1}_{\varepsilon}$. If
$\hat{\Sigma}^{-1}_{\varepsilon}$ converges to $\Sigma^{-1}_{\varepsilon}$ as
observation increases, the last two terms in (30) become under control, thus
we can obtain the same argument with extra constant terms.

##### Proof of Theorem 3

It is straightforward that the test statistic of the OLS method in (10) has a
$\chi^{2}_{p^{2}}$ distribution, where the degrees of freedom $p^{2}$ comes
from the difference in dimensionality of ${\Theta}_{0}$ and $\hat{\Theta}$.
Therefore, we get an asymptotic level $\alpha$ test if the null hypothesis is
rejected for $T(J)>\chi^{2}_{p^{2};(1-\alpha)}$, where
$\chi^{2}_{p^{2};(1-\alpha)}$ is the $(1-\alpha)$-quantile of chi-square
distribution with $p^{2}$ degrees of freedom. Using the threshold established
above, under the alternative, an upper bound on the power of the OLS method
can be obtained as

$\displaystyle P\bigg{(}T(J)>\chi^{2}_{p^{2};(1-\alpha)}\bigg{)}$
$\displaystyle=P\bigg{(}\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}+\Big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2}\Big{\\}}\geq\chi^{2}_{p^{2};(1-\alpha)}\bigg{)}$
$\displaystyle=P\bigg{(}\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}+\Big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2}\Big{\\}}-p^{2}\geq\chi^{2}_{p^{2};(1-\alpha)}-p^{2}\bigg{)}$
(31)
$\displaystyle\leq\frac{E(\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2})}{\chi^{2}_{p^{2};(1-\alpha)}-p^{2}}$
(32)
$\displaystyle\approx\frac{E(\|\boldsymbol{Y}_{J}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2})}{\frac{1}{2}\Big{(}z_{1-\alpha}+\sqrt{2p^{2}-1}\Big{)}^{2}-p^{2}},$
(33)

where $z_{1-\alpha}$ is the $(1-\alpha)$-quantile of Gaussian distribution.
The equality in (31) is obtained by subtracting
$E\Big{\\{}\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}{\boldsymbol{\Theta}}\|_{2}^{2}-\|\boldsymbol{Y}_{J}-\boldsymbol{X}_{J}\hat{\boldsymbol{\Theta}}\|_{2}^{2}\Big{\\}}=p^{2}$
from both sides, the inequality in (32) is obtained by using Markov’s
inequality and (33) is achieved as the quantile of chi-square distribution has
an approximation,
$\chi^{2}_{p^{2};(1-\alpha)}\approx\frac{1}{2}\Big{(}z_{1-\alpha}+\sqrt{2p^{2}-1}\Big{)}^{2}$.
Therefore, the upper bound on the power of the OLS method can be obtained as
in (15), which implies that
$E(\|\boldsymbol{Y}\|_{2}^{2}-\|\boldsymbol{Y}-\boldsymbol{X}{\boldsymbol{\Theta}}\|_{2}^{2})$
needs to be at least $O_{p}(p)$ to have power approaching to $1$.

##### Proof of Corollary 1

As we have
$\boldsymbol{Y}=\boldsymbol{E}+\text{vec}\big{(}\mathcal{X}^{(1)}(\boldsymbol{\theta}^{(1)}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime})\big{)}$
under the null rather than $\boldsymbol{Y}=\boldsymbol{E}$, the right-hand
side of the inequality in (21) can be represented as

$\displaystyle\max_{\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\},1\leq
i,j\leq
p}\Bigg{|}e^{\prime}_{i}\Bigg{(}\frac{{\mathcal{X}^{(2)}_{J}}^{\prime}\Big{(}E_{J}+\mathcal{X}_{J}^{(1)}\big{(}{\boldsymbol{\theta}^{(1)}}^{\prime}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}\Big{)}}{L}\Bigg{)}e_{j}\Bigg{|}$
$\displaystyle\leq\max_{\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\},1\leq
i,j\leq
p}\Bigg{|}e^{\prime}_{i}\Bigg{(}\frac{{\mathcal{X}^{(2)}_{J}}^{\prime}E_{J}}{L}\Bigg{)}e_{j}\Bigg{|}+\max_{\\{J:J\in\mathbb{J}_{T,p}(L),J\cap[\eta_{1},\eta_{2}]=\emptyset\\},1\leq
i,j\leq
p}\Bigg{|}e^{\prime}_{i}\Bigg{(}\frac{{\mathcal{X}^{(2)}_{J}}^{\prime}\mathcal{X}_{J}^{(1)}\big{(}{\boldsymbol{\theta}^{(1)}}^{\prime}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}}{L}\Bigg{)}e_{j}\Bigg{|},$
(34)

It is sufficient to show that both terms in (34) are less than or equal to
$C_{3}\lambda$ with probability approaching 1. The condition for the first
term is obtained from the proof of Theorem 1 and the one for the second term
is obtained from Assumption 5 and from the fact that
$\frac{\mathcal{X}^{(2)\prime}_{J}\mathcal{X}_{J}^{(1)}}{|J|}$ converges as
$T\rightarrow\infty$.

##### Proof of Corollary 2

It is sufficient to prove that (24) still holds where
$\boldsymbol{Y}_{J}=\text{vec}\big{(}\mathcal{Y}_{J}-\mathcal{X}^{(1)}_{J}{{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}$
is replaced by
$\boldsymbol{Y}_{J}^{\prime}=\text{vec}\big{(}\mathcal{Y}_{J}-\mathcal{X}^{(1)}_{J}{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}$.
The left-hand side of (24) can be simplified as

${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{E}_{J}+2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\text{vec}\Big{(}\mathcal{X}_{J}^{(1)}\big{(}{\boldsymbol{\theta}^{(1)}}^{\prime}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}\Big{)}>\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}.$
(35)

As shown in the proof of Theorem 2, it is sufficient to show that
${\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\boldsymbol{X}_{J}{\boldsymbol{\Theta}}$
is bounded by
$\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}$ with
probability tending to $1$ as the last component in left-hand side of (35),

$2{\boldsymbol{\Theta}}^{\top}\boldsymbol{X}_{J}^{\top}\text{vec}(\mathcal{X}_{J}^{(1)}\big{(}{\boldsymbol{\theta}^{(1)}}^{\prime}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)})=2\boldsymbol{\Theta}^{\top}\text{vec}\Big{(}{\mathcal{X}^{(2)}_{J}}^{\prime}\mathcal{X}_{J}^{(1)}\big{(}{\boldsymbol{\theta}^{(1)}}^{\prime}-{\hat{\boldsymbol{\theta}}^{(1)}}^{\prime}\big{)}\Big{)},$

is less than
$\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}$ with
probability approaching to $1$ from Assumption 5 and also from the fact that
$\frac{{\mathcal{X}^{(2)}_{J}}^{\prime}\mathcal{X}_{J}^{(1)}}{|J|}$ converges
as $T\rightarrow\infty$. Following the same logic presented in the proof
Theorem 2, it can be shown that the first component in left-hand side of (35)
is greater than
$\lambda\|\boldsymbol{\Theta}\|_{1}+\lambda^{\textsuperscript{thr}}$ with
probability approaching to $1$ which completes the proof.

##  2 Additional Simulation Results

In this section, additional simulation results are reported. As mentioned in
the main paper, our method is compared with the one proposed in Safikhani and
Shojaie, (2020) available from
https://github.com/abolfazlsafikhani/SBDetection. We first present a new
simulation scenario that is similar to the one used in Safikhani and Shojaie,
(2020), then give the additional results for those scenarios examined in the
main paper. Regarding the tuning parameters for Safikhani and Shojaie, (2020),
we follow the recommendation of their paper by using the default ones in
Section 2.1 where the simulation setting is a slightly modified version of
scenario 1 of Safikhani and Shojaie, (2020). However, the anomalies presented
in Section 2.2 are harder to detect as the size of change in coefficient
matrix is smaller and the noise has a larger variance compared to the one in
Section 2.1. Thus, to improve the performance of their method, we adjust
tuning parameters rather than using the default ones and the details can be
found in Section 2.2

###  2.1 Stronger signal-to-noise ratio and larger change size

####  2.1.1 Sparse $A^{(1)}$

Figure 7: The underlying coefficient matrices, $(A^{(1)},A^{(2)},A^{(1)})$, for the simulation setting in Section 2.1.1, where $A^{(2)}$ corresponds to an anomaly. T | p | $[\eta_{1},\eta_{2}]$ | $\eta_{2}-\eta_{1}$ | $\text{nzr}(A^{(1)})$ | $\text{nzr}(A^{(2)})$ | $\|\boldsymbol{\Theta}\|_{0}$ | $\Sigma_{\varepsilon}$
---|---|---|---|---|---|---|---
500 | 20 | $[T(1/3),T(2/3)]$ | $167$ | $-0.6$ | $0.75$ | $19$ | $0.01\textbf{{I}}_{p}$
Table 10: Simulation setting for Section 2.1.1, where $\text{nzr}(A^{(1)})$
and $\text{nzr}(A^{(2)})$ are the non-zero elements of $A^{(1)}$ and
$A^{(2)}$, respectively and $\|\boldsymbol{\Theta}\|_{0}$ is the number of
non-zero elements of $\boldsymbol{\Theta}$.

We borrow the simulation setting of scenario 1 used in Safikhani and Shojaie,
(2020). To make the single anomaly setting, we slightly modify the original
setting by changing the size of non-zero coefficients to $(-0.6,0.75,-0.6)$
for those intervals divided by an anomaly, whereas Safikhani and Shojaie,
(2020) consider the two change points with the corresponding size of non-zero
coefficients $(-0.6,0.75,-0.8)$. The details of the anomaly are given in Table
10, and we can see that the size of change is larger, the length of anomaly is
longer and the signal-to-noise ratio is larger than the simulation setting in
Section 5.3.

|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
$\Sigma_{\varepsilon}$ is known | random | OLS | 100 | 100
(s = 499) | Lasso | 100 | 100
deterministic | OLS | 100 | 100
(s = 499) | Lasso | 100 | 100
$\Sigma_{\varepsilon}$ is unknown | random | OLS | 100 | 100
(s = 469) | Lasso | 100 | 100
deterministic | OLS | 100 | 100
(s = 469) | Lasso | 100 | 100
Table 11: Empirical power ($\%$) from 100 simulation runs for the settings
described in Section 2.1.1, where $s$ is the number of intervals examined.

Comparing the simulation results of ours with Safikhani and Shojaie, (2020),
Tables 15-13 show that all methods detect one anomaly in all $100$ runs (this
is shown as “one” anomaly for the default and the lasso methods and “two”
change-points for the method of Safikhani and Shojaie, (2020)). In terms of
the localisation, the default and the lasso methods work better than Safikhani
and Shojaie, (2020) as they have smaller mean and sd of Hausdorff distance. We
emphasise that the anomaly presented in this section is easier to detect than
those used in the main paper in the sense that the size of change in
coefficient matrix is larger, the width of anomaly is longer and the noise
variance is smaller.

|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
$\Sigma_{\varepsilon}$ is known | random | OLS | 1.07 (0.10) | 0.80 (0.00)
(s = 499) | Lasso | 1.08 (0.10) | 0.80 (0.00)
deterministic | OLS | 0.40 (0.02) | 0.40 (0.03)
(s = 499) | Lasso | 0.40 (0.02) | 0.39 (0.04)
$\Sigma_{\varepsilon}$ is unknown | random | OLS | 1.14 (0.09) | 0.80 (0.00)
(s = 469) | Lasso | 1.03 (0.08) | 0.80 (0.00)
deterministic | OLS | 0.39 (0.03) | 0.40 (0.03)
(s = 469) | Lasso | 0.39 (0.03) | 0.39 (0.05)
Table 12: The mean (standard deviation) of Hausdorff distance from 100 simulation runs for the settings described in Section 2.1.1, where $s$ is the number of intervals examined. Empirical power ($\%$) |  | mean (sd) of Hausdorff distance
---|---|---
100 |  | 2.19 (1.55)
Table 13: Simulation results of Safikhani and Shojaie, (2020) under the
setting in Section 2.1.1.

####  2.1.2 Low rank + sparse $A^{(1)}$

Figure 8: The underlying coefficient matrices, $(A^{(1)},A^{(2)},A^{(1)})$, for the simulation setting in Section 2.1.2, where $A^{(2)}$ corresponds to an anomaly. T | p | $[\eta_{1},\eta_{2}]$ | $\eta_{2}-\eta_{1}$ | $\text{nzr}(A^{(1)}_{\text{sps}})$ | $\text{nzr}(A^{(2)}_{\text{sps}})$ | Range(low-rank) | $\Sigma_{\varepsilon}$
---|---|---|---|---|---|---|---
300 | 20 | $[T(1/3),T(2/3)]$ | $100$ | $-0.8682672$ | $0.8682672$ | $(-0.217,0.212)$ | $0.01\textbf{{I}}_{p}$
Table 14: Simulation setting for Section 2.1.2, where
$\text{nzr}(A^{(1)}_{\text{sps}})$ and $\text{nzr}(A^{(2)}_{\text{sps}})$ are
the non-zero elements in sparse part of $A^{(1)}$ and $A^{(2)}$, respectively
and Range(low-rank) is the range of non-zero elements in low-rank part.

In this section, we borrow the simulation setting of scenario $A.2$ used in
Bai et al., (2020) to compare their performance with ours. The underlying VAR
coefficient matrix has the low rank plus sparse structure and only the sparse
part (i.e. 1-off diagonal in this setting) undergoes change at anomaly. The
true coefficient matrices are presented in Figure 14. As all elements of
$A^{(1)}$ are non-zero, we estimate $A^{(1)}$ by using the ridge penalty when
$A^{(1)}$ is assumed to be unknown.

|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
$\Sigma_{\varepsilon}$ is known | random | OLS | 100 | 100
(s = 845) | Lasso | 100 | 100
deterministic | OLS | 100 | 100
(s = 845) | Lasso | 100 | 100
$\Sigma_{\varepsilon}$ is unknown | random | OLS | 100 | 100
(s = 941) | Lasso | 100 | 100
deterministic | OLS | 100 | 100
(s = 941) | Lasso | 100 | 100
Table 15: Empirical power ($\%$) from 100 simulation runs for the settings described in Section 2.1.2, where $s$ is the number of intervals examined. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
$\Sigma_{\varepsilon}$ is known | random | OLS | 0.67 (0.00) | 1.47 (0.52)
(s = 845) | Lasso | 0.67 (0.00) | 1.55 (0.66)
deterministic | OLS | 0.33 (0.00) | 0.33 (0.00)
(s = 845) | Lasso | 0.33 (0.00) | 0.34 (0.07)
$\Sigma_{\varepsilon}$ is unknown | random | OLS | 0.67 (0.00) | 1.74 (0.63)
(s = 941) | Lasso | 0.67 (0.00) | 1.31 (0.55)
deterministic | OLS | 0.61 (0.13) | 0.65 (0.08)
(s = 941) | Lasso | 0.38 (0.12) | 0.44 (0.15)
Bai et al., (2020) | 2.11 (1.59)
Table 16: The mean (standard deviation) of Hausdorff distance from 100
simulation runs for the settings described in Section 2.1.2, where $s$ is the
number of intervals examined. $\\#$ of estimated change points from 100 runs
---
2 | 3
80 | 20
Table 17: Simulation results of Bai et al., (2020) under the setting in
Section 2.1.2.

###  2.2 Simulation results when $\Sigma_{\varepsilon}$ is estimated

In this section, we repeat the simulation settings introduced in Section 5 for
the case when $\Sigma_{\varepsilon}$ is unknown. We use the maximum likelihood
estimator for $\hat{\Sigma}_{\varepsilon}$ and compare the performance of our
method with the change-point detection technique proposed by Safikhani and
Shojaie, (2020)). As mentioned earlier, the default tuning parameters
recommended by Safikhani and Shojaie, (2020) do not fit well in the simulation
settings presented in the following sections. Among three tuning parameters
$\lambda_{1}$, $\lambda_{2}$ and $\omega$, we adjust $\lambda_{1}$ and
$\omega$; a larger range is examined for finding the optimal $\lambda_{1}$ in
the initial break detection stage and a larger value of $\omega=4\log T\log p$
is used (instead of the default constant $1/1.75$) to allow smaller number of
break points to be selected in the screening stage.

####  2.2.1 Dense $A^{(1)}$, Single anomaly

|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 25 | 13
(s = 1969) | Lasso | 100 | 78
deterministic | OLS | 20 | 29
(s = 1969) | Lasso | 100 | 83
deterministic | OLS | 19 | 28
(s = 981) | Lasso | 100 | 82
case 2 | random | OLS | 7 | 6
(s = 1969) | Lasso | 99 | 34
deterministic | OLS | 10 | 14
(s = 1969) | Lasso | 100 | 34
deterministic | OLS | 10 | 13
(s = 981) | Lasso | 100 | 32
Table 18: Empirical power ($\%$) from 100 simulation runs for the settings described in Section 5.2.1, where $s$ is the number of intervals examined. | # of estimated change points
---|---
| 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 17
case 1 | 0 | 1 | 4 | 6 | 16 | 19 | 11 | 10 | 13 | 10 | 4 | 3 | 1 | 2
case 2 | 1 | 1 | 7 | 14 | 11 | 25 | 10 | 15 | 7 | 4 | 3 | 1 | 0 | 1
Table 19: Distribution of the number of estimated change-points by Safikhani and Shojaie, (2020) under the simulation setting described in Section 5.2.1 over 100 simulation runs. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 35.41 (18.24) | 43.00 (8.79)
(s = 1969) | Lasso | 2.67 (3.85) | 18.28 (17.75)
deterministic | OLS | 37.72 (16.39) | 35.73 (17.66)
(s = 1969) | Lasso | 1.84 (4.32) | 12.23 (18.11)
deterministic | OLS | 38.33 (15.92) | 36.62 (16.94)
(s = 981) | Lasso | 2.00 (4.33) | 13.95 (19.05)
case 2 | random | OLS | 43.86 (11.50) | 45.83 (6.16)
(s = 1969) | Lasso | 3.74 (6.07) | 36.70 (17.51)
deterministic | OLS | 42.58 (13.38) | 42.47 (13.23)
(s = 1969) | Lasso | 2.50 (5.15) | 36.80 (18.29)
deterministic | OLS | 42.80 (13.06) | 42.97 (12.47)
(s = 981) | Lasso | 2.28 (4.70) | 37.15 (18.01)
Table 20: The mean (standard deviation) of Hausdorff distance from 100 simulation runs for the settings described in Section 5.2.1, where $s$ is the number of intervals examined. | mean (sd) of Hausdorff distance
---|---
case 1 | 22.84 (4.64)
case 2 | 23.46 (4.45)
Table 21: The mean (standard deviation) of Hausdorff distance from 100
simulation runs obtained by Safikhani and Shojaie, (2020) under the simulation
setting described in Section 5.2.1.

####  2.2.2 Dense $A^{(1)}$, Two anomalies

| # of estimated change points
---|---
| $\mathbf{4}$ | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13
case 1 | 5 | 5 | 16 | 18 | $\mathbf{21}$ | 18 | 12 | 4 | 1 | 0
case 2 | 4 | 12 | 11 | 21 | $\mathbf{24}$ | 13 | 7 | 5 | 2 | 1
Table 22: Distribution of the number of estimated change points by Safikhani and Shojaie, (2020) under the simulation setting described in Section 5.2.2 over 100 simulation runs. |  |  |  $A^{(1)}$ is known |  |  $A^{(1)}$ is estimated
---|---|---|---|---|---
|  |  | 0 | 1 | 2 | 3 |  | 0 | 1 | 2 | 3
case 1 | random | OLS | 55 | 45 | 0 | 0 |  | 48 | 49 | 3 | 0
(s = 1969) | Lasso | 0 | 69 | 31 | 0 |  | 5 | 94 | 1 | 0
deterministic | OLS | 23 | 65 | 12 | 0 |  | 37 | 54 | 9 | 0
(s = 1969) | Lasso | 0 | 65 | 33 | 2 |  | 5 | 92 | 3 | 0
deterministic | OLS | 28 | 61 | 11 | 0 |  | 42 | 55 | 3 | 0
(s = 981) | Lasso | 0 | 66 | 34 | 0 |  | 5 | 91 | 4 | 0
case 2 | random | OLS | 72 | 28 | 0 | 0 |  | 86 | 14 | 0 | 0
(s = 1969) | Lasso | 0 | 5 | 93 | 2 |  | 38 | 62 | 0 | 0
deterministic | OLS | 71 | 23 | 6 | 0 |  | 75 | 23 | 2 | 0
(s = 1969) | Lasso | 0 | 5 | 95 | 0 |  | 36 | 64 | 0 | 0
deterministic | OLS | 81 | 16 | 3 | 0 |  | 82 | 17 | 1 | 0
(s = 981) | Lasso | 0 | 6 | 94 | 0 |  | 36 | 64 | 0 | 0
Table 23: Distribution of the number of detected anomalies for two methods in all cases described in Section 5.2.2 over 100 simulation runs, where $s$ is the number of intervals examined. | mean (sd) of Hausdorff distance
---|---
case 1 | 31.42 (3.70)
case 2 | 20.25 (7.08)
Table 24: The mean (standard deviation) of Hausdorff distance from 100 simulation runs obtained by Safikhani and Shojaie, (2020) under the simulation setting described in Section 5.2.2. |  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
---|---|---|---|---
case 1 | random | OLS | 29.74 (8.38) | 30.01 (7.85)
(s = 1969) | Lasso | 7.27 (4.38) | 19.02 (9.87)
deterministic | OLS | 28.51 (10.91) | 29.69 (9.50)
(s = 1969) | Lasso | 5.95 (3.33) | 17.81 (12.08)
deterministic | OLS | 29.02 (10.45) | 30.96 (7.39)
(s = 981) | Lasso | 5.72 (3.16) | 18.05 (12.16)
case 2 | random | OLS | 13.60 (0.18) | 13.54 (0.72)
(s = 1969) | Lasso | 3.42 (5.13) | 12.41 (1.68)
deterministic | OLS | 13.14 (2.77) | 13.44 (1.20)
(s = 1969) | Lasso | 2.35 (2.71) | 12.72 (1.64)
deterministic | OLS | 13.27 (1.76) | 13.41 (1.17)
(s = 981) | Lasso | 2.47 (2.89) | 12.62 (1.76)
Table 25: The mean (standard deviation) of Hausdorff distance from 100
simulation runs for the settings described in Section 5.2.2, where $s$ is the
number of intervals examined.

####  2.2.3 Sparse $A^{(1)}$, Single anomaly

|  |  |  $\\#$ anomaly detection from 100 runs
---|---|---|---
|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
case 1 | random | OLS | 2 | 3
(s = 469) | Lasso | 100 | 99
deterministic | OLS | 73 | 73
(s = 469) | Lasso | 100 | 100
case 2 | random | OLS | 2 | 1
(s = 469) | Lasso | 85 | 33
deterministic | OLS | 57 | 66
(s = 469) | Lasso | 100 | 52
Table 26: Empirical power ($\%$) from 100 simulation runs for the settings described in Section 5.3, where $s$ is the number of intervals examined. | # of estimated change points
---|---
| 0 | 1 | 2 | 3 | 4 | 5
case 1 | 1 | 34 | 51 | 12 | 1 | 1
case 2 | 1 | 29 | 56 | 12 | 1 | 1
Table 27: Distribution of the number of estimated change points by Safikhani and Shojaie, (2020) under the simulation setting described in Section 5.3 over 100 simulation runs. |  |  | mean (sd) of Hausdorff distance
---|---|---|---
|  |  |  $A^{(1)}$ is known |  $A^{(1)}$ is estimated
case 1 | random | OLS | 44.44 (2.55) | 43.92 (6.08)
(s = 469) | Lasso | 1.50 (0.67) | 1.75 (4.36)
deterministic | OLS | 20.72 (19.56) | 23.43 (19.20)
(s = 469) | Lasso | 0.32 (0.15) | 0.32 (0.15)
case 2 | random | OLS | 45.93 (3.47) | 46.35 (0.50)
(s = 469) | Lasso | 12.97 (14.19) | 31.43 (21.44)
deterministic | OLS | 29.41 (19.33) | 30.27 (17.77)
(s = 469) | Lasso | 0.46 (0.15) | 22.50 (23.08)
Table 28: The mean (standard deviation) of Hausdorff distance from 100 simulation runs for the settings described in Section 5.3, where $s$ is the number of intervals examined. | mean (sd) of Hausdorff distance
---|---
case 1 | 45.65 (1.84)
case 2 | 45.82 (2.46)
Table 29: The mean (standard deviation) of Hausdorff distance from 100
simulation runs obtained by Safikhani and Shojaie, (2020) under the simulation
setting described in Section 5.3.