SAMoSSA:  Multivariate Singular Spectrum Analysis with Stochastic Autoregressive Noise
Authors: Abdullah Omar Alomar
Munther A. Dahleh
Sean Mann
Devavrat Shah
Venue: NeurIPS 2023 poster (2023)
Keywords: Time series
System Identification
Singular Spectrum Analysis

Abstract: The well-established practice of time series analysis involves estimating deterministic, non-stationary trend and seasonality components followed by learning the residual stochastic, stationary components. Recently, it has been shown that one can learn the deterministic non-stationary components accurately using  multivariate Singular Spectrum Analysis (mSSA) in the absence of a correlated stationary component; meanwhile, in the absence of deterministic non-stationary components, the Autoregressive (AR) stationary component can also be learnt readily, e.g. via Ordinary Least Squares (OLS). However, a theoretical underpinning of multi-stage learning algorithms involving both deterministic and stationary components has been absent in the literature despite its pervasiveness. We resolve this open question by establishing desirable theoretical guarantees for a natural two-stage algorithm, where mSSA is first applied to estimate the non-stationary components despite the presence of a correlated stationary AR component, which is subsequently learned from the residual time series. We provide a finite-sample forecasting consistency bound for the proposed algorithm, SAMoSSA, which is data-driven and thus requires minimal parameter tuning. To establish theoretical guarantees, we overcome three hurdles: (i) we characterize the spectra of Page matrices of stable AR processes, thus extending the analysis of mSSA; (ii) we extend the analysis of AR process identification in the presence of arbitrary bounded perturbations; (iii) we characterize the out-of-sample or forecasting error, as opposed to solely considering model identification. Through representative empirical studies, we validate the superior performance of SAMoSSA compared to existing baselines. Notably, SAMoSSA's ability to account for AR noise structure yields improvements ranging from 5% to 37% across various benchmark datasets.

SAMoSSA: Multivariate Singular Spectrum Analysis with Stochastic Autoregressive Noise
Abdullah Alomar MIT aalomar@mit.co edu
Munther Dahleh MIT dahleh@mit. edu
Sean Mann MIT seanmann@mit.cm edu
Devavrat Shah MIT devavrat@mit.en edu
Abstract
The well-established practice of time series analysis involves estimating determin- istic, non-stationary trend and seasonality components followed by learning the residual stochastic, stationary components. Recently, it has been shown that one can learn the deterministic non-stationary components accurately using multivari- ate Singular Spectrum Analysis (mSSA) the absence of a correlated stationary component; meanwhile, in the absence of deterministic non-stationary components, the Autoregressive (AR) stationary component can also be learnt readily, e-g. via Ordinary Least Squares (OLS). However, a theoretical underpinning of multi-stage learning algorithms involving both deterministic and stationary components has been absent in the literature despite its pervasiveness. We resolve this open question by establishing desirable theoretical guarantees for a natural two-stage algorithm, where mSSA is first applied to estimate the non-stationary components despite the presence a correlated stationary AR component, which is subsequently learned from the residual time series. We provide finite-sample forecasting consistency bound for the proposed algorithm, SAMoSSA, which is data-driven and thus re- quires minimal parameter tuning. To establish theoretical guarantees, we overcome three hurdles: (i) we characterize the spectra of Page matrices of stable AR pro- cesses, thus extending the analysis of mSSA: (ii) we extend the analysis of AR process identification in the presence of arbitrary bounded perturbations; (iii) we characterize the out-of-sample or forecasting error, as opposed to solely consid- ering model identification. Through representative empirical studies, we validate the superior performance of SAMOSSA compared existing baselines. Notably, SAMoSSA's 's ability to account for AR noise structure yields improvements ranging from 5% to 37% % across various benchmark datasets.
Introduction
Background. Multivariate time series have often been modeled as mixtures of stationary stochastic processes (c.g. AR process) and deterministic non-stationary components (e.g. polynomial and harmonics). To handle such mixtures, classical time series forecasting algorithms first attempt to estimate and then remove the non-stationary components. For example, before fitting an Autoregres- sive Moving-averin and then removed from the time series. Once the non-stationary components have been eliminated" the ARMA model is learned. This estimation procedure often relies on domain knowledge and/or fine-tuning, and theoretical analysis for such multiple-stage algorithms is limited the literature.
(ARMA) model, polynomial trend and seasonal components must be estimated
Prior work presented a solution for estimating non-stationary deterministic components without domain knowledge or fine-tuni tne-tuning using mSSA 3 This framework systematically models wide
Note that Autoregressive Integrated Moving Moving-average average (ARIMA) and Seasonal ARIMA models use differenc- ing in conjunction with unit root tests to remove stochastic non- stationary components, and are not suited for the model we consider.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

class of (deterministic) non-stationary multivariate time series as linear recurrent formulae (LRF). encompassing wide class spatio-temporal factor models that includes harmonics, exponentials, and polynomials. However, mSSA, both algorithmically and theoretically, does not handle additive correlated stationary noise, an important noise structure time series analysis. Indeed, all theoretical results of mSSA are established in the noiseless setting [16] or under the assumption of independent and identically distributed (i.idd)) noise [4.
On the other hand, every stable stationary process can be approximated as a finite order AR process 8 26 The classical OLS procedure has been shown to accurately learn finite order AR processes 17 However, the feasibility of identifying AR processes in the presence of non-stationary deterministic component has not been addressed.
In summary, despite the pervasive practice of first estimating non-stationary deterministic compo- nents and then learning the stationary residual component, neither elegant unified algorithm nor associated theoretical analyses have been put forward in the literature.
A step towards resolving these challenges is answer the following questions: (i) Can mSSA consistently estimate non-stationary deterministic components in the presence of correlated stationary noise? (ii) Can the AR model be accurately identified using the residual time series, after removing the non-stationary deterministic components estimated by mSSA, potentially with errors? (iii) Can the out-of oft-of-sample sample forecasting error of such multi-stage algorithm be analyzed?
In this paper, we resolve all three questions in the affirmative: we present SAMoSSA, a two-stage procedure, where in the first stage, we apply mSSA on the observations to extract the non-stationary components; and in the second stage, the stationary AR process is learned using the residuals.
Setup. We consider the discrete time setting where we observe a multivariate time series (t)) = [ya (t),... ..... yN (t)] RN at each time index & E [T] {1..... T} where T 2 N For each n E [V]. and each timestep t E [T]. the observations take the form
yn (t)= (t) + (1).
(1)
where fn =z + R denotes the non-stationary deterministic component, and v.. (4) is a stationary AR noise process of order pr. (AR(p.)). Specifically, each (1) is governed by
(t) ani (t - 1) + 17 (t),
where on (l) refers the per-step noise, modeled as mean-zero i.i.d. random variables, and Ozy RVi i E [Pn] are the parameters of the n-th AR(p,) process.
Goal. For each n [N], our objective is threefold. The first is estimating fn (t( from the noisy observations ya (t) for all t E (T]. The second is identifying the AR process' parameters Ond V [p.]. The third is out-of-sample forecasting of yn (t) for >>> T.
Contributions. The main contributions of this work is SAMoSSA, elegant two-stage algorithm, which manages to learn both non-stationary deterministic and stationary stochastic components of the underlying time series. A detailed summary of the contributions as follows.
(a) Estimating non-stationary component with mSSA under AR noise. The prior theoretical results for mSSA are established in the deterministic setting [16] or under the assumption of i.i.d. noise 4 4.3.. In Theorem 4.1 we establish that the mean squared estimation error scales as ~ 1/ 1/VNT under AR noise the same rate was achieved by [3] under i..dd noise. Key to this result is establishing spectral properties of the "Page" matrix of AR processes, which may be of interest in its own right (see Lemma
(b) Estimating AR model parameters under bounded perturbation. We bound the estimation error for the OLS estimator of AR model parameters under arbitrary and bounded observation noise, which could be of independent interest. Our results build upon the recent work of [17] which derives similar results but without any observation noise. In that sense, ours can be viewed robust generalization of [17]
if I' < N. we divide the N time series into [N/1] sets of time series where this condition will hold.

---

(c) Out-of-samee finite- consistency of SAMOSSA. We provide finite- sample analysis for the forecasting error of SAMoSSA- such analysis of two-stage procedures in this setup nascent despite its ubiquitous use in practice. Particularly, we establish in Theorem 4.4 that the out-of- sample forecasting error for SAMoSSA for the T time-steps ahead scales as probability.
with high
VNT
(d) Empirical results. We demonstrate superior performance of SAMoSSA using both real-world and synthetic datasets. We show that accounting for the AR noise structure, as implemented in SAMOSSA, consistently enhances the forecasting performance compared to the baseline presented in [3. Specifically, these enhancements range from 5% to 37% across benchmark datasets.
Related Work. Time series analysis is a well-developed field. We focus on two pertinent topics.
mSSA. Singular Spectrum Analysis (SSA), and its multivariate extension mSSA, are well-studied methods for time series analysis which have been used heavily in wide array of problems including
imputation 4 3 forecasting [18 2). and change point detection [22 6 Refer to [66 15] for a
good overview of the SSA literature. The classical SSA method consists of the following steps: (1) construct a Hankel matrix of the time series of interest; (2) apply Singular Value Decomposition (SVD) on the constructed Hankel matrix, (3) group the singular triplets to separate the different components of the time series; (4) learn a linear model for each component to forecast. mSSA is an extension of SSA which handles multiple time series simultaneously and attempts exploit the shared structure between them [8]. The only difference between mSSA and SSA is in the first step, where Hankel matrices of individual series are "stacked" together to create single stacked Hankel matrix. Despite its empirical success, mSSA's 's classical analysis has mostly focused on identifying which time series have low-rank Hankel representation, and defining sufficient asymptotic conditions for signal extraction, i.e., when the various time series components are separable. All of the classical analysis mostly focus on the deterministic case, where no observation noise is present.
Recently, a variant of mSSA was introduced for the tasks of forecasting and imputation BI. This variant, which we extend in this paper, uses the Page matrix representation instead of the Hankel matrix, which enables the authors to establish finite-sample bounds on the imputation and forecasting error. However, their work assumes the observation noise to be i.i.d., and does not accommodate correlated noise structure, which is often assumed in the time series literature [26]. Our work extends the analysis in Agarwal et al. [ to observations under AR noise structure. We also extend the analysis of the forecasting error by studying how well we can learn (and forecast) the AR noise process with perturbed observations.
Estimating AR parameters. AR processes are ubiquitous and of interest in many fields, including time series analysis, control theory, and machine learning. In these fields, it is often the goal to estimate the parameters of an AR process from a sample trajectory. Estimation is often carried out through OLS, which is asymptotically optimal [II]. The asymptotic analysis of the OLS estimator is well established, and recently, given the newly developed statistical tools of high dimensional probability, cf.
well. For example, several works established results for the finite-time identification for general first order vector AR systems (20 19 27 24 21.. Further, González al. (17] provides finite-time bound on the deviation of the OLS estimate for general AR(p) processes. To accommodate our setting, we extend the results of 17 in two ways. First, we extend the analysis accommodate any sub-gaussian noise instead of assuming gaussianity; Second, and more importantly, we extend it to handle arbitrary bounded observation errors in the sampled trajectory, which can be of independent interest.
B1 29). various recent works have tackled its finite-tmme analysis
2 Model
Deterministic Non-Stationary Component. We adopt the spatio-temporal factor model of 3 3 described next. The spatio-temporal factor model holds when two assumptions are satisfied: the first assumption concerns the spatial structure, i.e., the structure across the N time series fi, J'v the second assumption pertains to the "temporal" structure. Before we describe the assumptions, we first define a key time series representation: the Page matrix.
Definition 2.1 (Page Matrix). Given time series f 2 + R. and an initial time index to > 0, the Page matrix representation over the T entries [to)...., foo,f((ot -1) with parameter < L S T is given by the matrix [J, L,, T, 1o) E [LX(T/L) with Z(J, I,, (o))) = f(0 +i + (j - 1) x I.) for i E [L], j E [[]
3

---

In words, to construct the L x [T/L] Page matrix of the entries f(fo)..... f(o + T 1), partition them into [ segments of L contiguous entries and concatenate these segments column-wise (see Figure D. Now we introduce the main assumptions for .1.....JN:
Assumption 2.1 (Spatial structure). For each n E [N], fa(t) -Et Uny Wy (t) for some Uur E IR and w, Z R.
Assumption 2.1 posits that each time series (.)() An E [M] can be described as a linear combination of R "fundamental" time series. The second assumption relates to the temporal structure of the fundamental time series w.)() Vr E [RC which we describe next.
Assumption 2.2 (Temporal structure). For each r E [R] and for any T 1, 1 1 L - T, 4, >
rank(Z(w; w., L, (,T) S G.
Assumption 2.2 posits that the Page matrix of each "fundamental" time series w, has finite rank. While this imposed temporal structure may seem restrictive first, it has been shown that many standard functions that model time series dynamics satisfy this property [4] 3]. These include any finite sum of products of harmonics, low-degree polynomials, and exponential functions (refer to Proposition 2.1 in (1).
Stochastic Stationary Component. We adopt the following two assumptions about the AR processes Ta Vn [AV]. We first assume that these AR processes are stationary (see Definition H.4 in Appendix H
Assumption 2.3 (Stationarity and distinct roots). Iy (t) is a stationary AR(pa) process Vn [N]. That is, let Ans E C, E [pm] denote the roots of (n(2) == 2P - On Then, [xae < Vi E [P.]. Further, Vn E [M] the roots of (2) are distinct.
Further, assume the per-step noise nn (t) are i.i.d. sub-gaussian random variables (see Definition H.1 in Appendix H).
Assumption 2.4 Sub-gaussian noise). For n E [N], t E [T), (t())) are zero-mean i.i.d. sub- sub-gaussian random variables with variance o².
Model Implications. In this section, state two important implications of the model stated above. The first is that the stacked Page matrix defined as
7 T, (o) [f,,, t,T,to)
L,
Z(fn, L.fn,,TT,tol L,T,to)),
low-rank. Precisely, recall the following Proposition stated in 3].
Proposition 2.1 (Proposition in (3) Let Assumptions 2.1 and 2.2 hold. Then for any L S LVT with any T Z 1, to > O, the rank of the Page matrix ZJJ L, T to for n E [N] is at most R x G.
Further, the rank of the stacked Page matrix (L, T, (o) is R
Throughout, will use the shorthand Z, = Z₁ T, and Zy(()) 7 (L, (o) The second implication is that there exists a linear relationship between the last row of 2 and its top L rows. The following proposition establishes this relationship. First, Let [Z,]æ. denote the L-th row of Z and Z', E IR LL- (N denote the sub-matrix that consist of the top L rows of Zj.
Proposition 2.2 (Proposition 3 in 3 ). Let Assumptions 2.7 and hold. Then, for L > RG, there exists 8 E TR L-1 such that [T = Z 3 Further |8 1o S RG.
3 Algorithm
The proposed algorithm provides two main funetionalities. The first one is decomposing the observa- tions yn (t) into estimate of the non-stationary and stationary components for < T. The second is forecasting yo (1) for / > T, which involves learning a forecasting model for both (a( and x (1).
Univariate Case. For ease of exposition, we will first describe the algorithm for the univariate case (N = 1). The algorithm has the following parameters: 1 < L < VT, and < < L (refer
We will use the same shorthand for the stacked page matrices of and x: namely, Zy and Z...

---

to Appendix B.2 for how to choose these parameters). For clarity, we will assume, without loss of generality, that L is chosen such that T/L is an integer In the first step of the algorithm, we transform the observations yı (t), E (T) into the I. x T/ Page matrix Z(nn I,, 1)). We will use the shorthand Z = Zyy ; L,, T, henceforth.
Decomposition. We compute the SVD Zyn = Ei= Syuge, where $1 2 32 $2... 2 SL 2 denote its ordered singular values, and uy E B², Us
; E RR/L denote its left and right singular vectors, respectively, for & E [4]. Then, we obtain 2 by retaining the top singular components of y (i.e., by applying Hard Singular Value Thresholding (HSVT) with threshold k). That is, 2, E/ Spuer.
We denote by (t) and क (t) the estimates of (t) and x (t) respectively. We read off the estimates (t) directly from the matrix 2,, using the entry that corresponds to t E [7]. More precisely, for t E [T), (t) equals the entry of Z in row (t - mod L) + 1 and column [t/L]. while (1) - y1 (t) (t).
Forecasting. To forecast y1 for t > T. we produce a forecast for both (1()) and (1).. Both forecasts are performed through linear models (3 for (t) and â2 for 3n (t).) For fi(,, we first learn linear model 3 defined as
T/L Σ (w((Lm) TTT F,)), m==1
3 argminseRt-1
where F,,. [L 1) + 1),.... Rl 1)] for
[T/ We then use 3 and the lagged observations to produce (t). That is fi(t) Yi (t 1), where Yi 1)
[vn - 1),
Z)].
For (t), we first estimate the parameters for the AR(pi) process. Specifically, define the P1- dimensional vectors (1) == (1), Eo (t Pr 1) Then, define the OLS estimate as,
a = argmin
Σ
(1 -
Then, let ((') y1 (t) (((e) for any 1' < t. We then use ā and the pi lagged entries of ž (.) to produce (t). That is (t) = a_ Xi 1), where Xi - 1) = [ - 1).... T (t - p1 p1)]- Finally, produce the forecast 01 (t) fi(t) + (t).
Multivariate Case. The key change for the case N is the use of the stacked Page matrix Zy. which the column-wise concatenation of the Page matrices induced by individual time series. Specifically, consider the stacked Page matrix Zy E RLXNT/L defined as L,
Zy - 7 L,T,1 Z(22,
Z(yn. L.T.
Decomposition. The procedure for learning the non-stationary component f1...., ......fx fN for t < T is similar to that of the univariate case. Specifically, we perform HSVT with threshold & on Z, to produce Zy. Then, we read off the estimate of (t) from Zy. Specifically, for t E [T', and E [N], let 2, refer to sub-matrix of 2, induced by selecting only its [(n - 1) (T/L) + 1,..., T LL]
'tteewwise, one can apply this algorithm to the two ranges {1,..., L x [T/L]} and {(T mod L) +
1. 1,...TT), .T).
In the forecasting algorithm, the estimates (L - 1))mm-1))),,.-..... x. 1)] are obtained by applying
HSVT sub-matris of Zy, which consists of its first L - rows. This is done establish the theoretical
results as it helps avoid dependencies in the noise between (Lm) and F... for m E [T]
'Note that we assume knowledge of the true parameter p1 (and p.. in the multivariate case).
"Note the subtle difference between 7, (t) and (t) precisely, for >> T. Ty (t) is an estimate of (t)
before observing y1 (t) (i.e., forecast), whereas Z (t) estimate of (t) after observing (t).
5

---

III. Learn
X60
X,0
ND1
Jun
7:20
A00
7600
3,03
2007
Fd--
10- 11
A0 10
Juk 10
301
2ylo
3000
1201
Stacked Page
Matrix transformation
Leern
NO x,C 1)
500
MVD
9,09
9,01
9p7-p1
2,01
$,00
2177)
9,0-
R,00
7,07
2,
II. HSVT
3,50
singulie
V. Forecast for t>T
L-X
6
7.20
-1- -1
2,00
2.09 Pg00
2g17)
Jun
2-- 2,
Figure visual depiction of SAMOSSA. The algorithm five steps are (i) transform the time series into its stacked Page matrix representation; (ii) decompose time series into non-stationary and stationary components (iii) estimate 3; (iv) estimates â, Vn [N]; (v) produce the forecast Zr (t) for > T.
columns. Then for E (7), (t) equals the entry of f, row (t - mod L) and column [t/L]. We then produce an estimate for n, (f) as &, (t) (t) ft()).
Forecasting. To forecast (t) for t > T. as in the univariate case, we learn a linear model 3 defined as
NXT/L Σ m=1
argmin
(ym 3TT
(6)
(2(().... 1/2 (T)..... .......................... (T) E where Um is the m-th component of [yn(L), 31 (2 L)...., (T)),, NXTLL and F., R corresponds to the vector formed by the entries of the first L rows in the moth column of Z, for m [N x T/L]. We then use 3 to produce (1) = 3 Y, 1), where yn (t-1). That is (t-1) = Vyn (t-11) 1)) ... yn (t- L)
again (t-1) is the vector of the L lags of yn
Then, for each [N], estimate agui Vi E [pm], the parameters for the n-th AR(p,) process. Let Gn denote the OLS estimate defined as
an argmin @RPP
Σ (
Then, produce forecast for Tn. (1) a X1 1), where again Xa(1 1) (L ,)....,ūr, Pw)]. Finally, produce the forecast for yn as (t) (t) + En (t). For a visual depiction of the algorithm, refer to Figure
Results
In this section, we provide sample high probability bounds on the following quantities:
1. Estimation error of non-stationary component. First, we give an upper bound for the estimation error of each one of ft),...... (t) for & [Z]. Specifically, we upper bound the following metric
1 T
EstErr( N, T,m) = + Σ
(8)
1=1
TTo be precise, 2, here the truncated SVD of sub-matrix of Zy which consist of its first L rows.

---

where fi(t) n [M]. E [T] are the estimates produced by the algorithm we proposed in Section
2. Identification of AR parameters. Second, analyze the accuracy of the estimates an Vn [N] produced by the proposed algorithm. Specifically, we upper bound the following metrics
a -anll2
Vn E [N],
where a (1 . Omp...
Out-of-sample forecasting error. Finally, we provide an upper bound for the forecasting error of y1 (t)...... (t( for E {I' + 1.... 2T'}. Specifically, we upper bound the following metric
2T Σ Σ
(N,Er(N,T)) = NT 1=1-1
(1) - E (v)(() vn - ..... M. (1))²
(10)
where V [N], E +T+ } are the forecasts produced by the algorithm we propose in Section 3 Before stating the main results, we state key additional assumptions.
Assumption 4.1 (Balanced spectra). Let (ta) denote the NT LxNT stacked Page matrix associated for their T consecutive entries starting from to- Let k (ttk()) and let (z (to)) denote the k-th singular value for (f((t) Then, for any to > 0,
with all N time series f...)..... (-)
2 (f)) is such that a (Z)()) 2 n/NT/ / V for some absolute constant 7 0.
This assumption holds whenever the the non-zero singular values are "well-balanced" a standard as- sumption in the matrix/tensor estimation literature 3 6. Note that this assumption, as stated, ensures balanced spectra for the stacked Page matrix of any set of T consecutive entries of fi f-()).
Finally, we will impose an additional necessary restriction on the complexity of the N time series [i!),... fx for T (as is done in 3 () Let Z', denote the (L-1)x (NT/L) LL) matrix formed using the top rows of Z,. Further, for any lo E [T + 1], let Z (l)) denote the (I.-1) x (NT/L) matrix formed using the top L - rows of 2 (fo). For any matrix M, let colspan(M)) denote the subspace spanned by the its columns. We assume the following property.
Assumption 4.2 (Subspace inclusion). For any to E [T' + 1]. colspan (t)())) C colspan (Z₃))
This assumption necessary as it requires the stacked Page matrix of the out-of-sample time series Z'y (fo) be only "rich" that of the stacked Page matrix of the "in-sample" time series Z'y.
4.1 Main Results
First, recall that R is defined in Assumption 2.1 G in Assumption 2.2 k is in Proposition 2.1 while y is defined in Assumption 4.1. Further, recall that AnS for Pnn are the roots of the characteristic polynomial of the n-th AR process, as defined in Assumption 2.3 Throughout, let and C absolute constants, fma == max,, [n)) := max, Pus Omax max, |aall2: and C(Jmax, 7) denote constant that depends only (polynomially) on model parameters fiax and Last but not least, define the key quantity
CXT
x.)
where X, max;,n Dass and ex is constant that depends only on A, Note that o2 a key quantity that we use to bound important spectral properties of AR processes, as Lemma A.2 establishes.
4.1.1 Estimation Error of Non-Stationary Component
Theorem 4.1 (Estimation error /). Assume access to the observations y1 (I), ... (N(()) fort E [T]
as defined in ( Let assumptions 2.1 2.2 2.3 2.4 and 4.1 hold. Let L /NT and le = k, then, for any n [N), with probability of at least
o2 VNT
EstErr T, T, S C(Smax: 7)
(11)
See Appendix oor an explicit expression of es.

---

Theorem 4.4. Let the conditions of Theorem and Assumption 4.2 hold. Then, with probability of least c This theorem implies that the mean squared error of estimating f scales Ô probability. The proof of Theorem 4.1 in Appendix C
with high
VNT
4.1.2 Identification of AR Processes
Theorem 4.2 (AR Identification). Let the conditions of Theorem 4.1 hold. Let an be as defined
in 7 and an = a111 anpal as defined in (2) Then, for sufficiently large T such that
Co²
Xmax (T5) and for any E [N] where EstErr(N, T,, x² Ami we have
Ammin Ammin(())
pot max
with probability of least
Tí..
Cp
(1108 log
( (Tmaax(W) Amia
(s, T, log(?) o² min o Amin (V))
an
anl12
(12)
Amin (v)
Recall that on is estimated using En (+), perturbed version of the true AR process In (-). This perturbation comes from the estimation error of (.,(),, which is consequence of EstErr(N, Hence, not surprising that EstErr(NN, I,n) shows up the upper bound. Indeed, the upper bound we provide here consists of two terms: the first, which scales Ô (+) is the bound one would get with access the true AR process x. (.).), shown in 17.. The second term characterizes the contribution of the perturbation caused by the estimation error, and it scales as O (EstErr(N, T, n)),
which we show is Ô
with high probability in Theorem 4.1 Also, note that the theorem VNT
holds when the estimation error is sufficiently small (EstErr(AV (,, T,n) S o² (²)m)n())) This condition ensures that the perturbation arising from the estimation error remains below the lower bound for the minimum eigenvalue of the sample covariance matrix associated with the autoregressive processes. Note that Amax (D), Amin (V) and Amax (T) are quantities that relate to the parameters of the A AR(p) processes and their "controllability Gramian"
Up
detail in Appendix A.I The proof of Theorem
4.2 is in Appendix [
4.1.3
Out-of-sample Forecasting Error
We first establish upper bound on the error of our estimate 3.
Theorem 4.3 (Model Identification (3*)). Let the conditions of Theorem 4.1 hold. Let be defined as in (6) and 3 as defined in Proposition 2.2 Then, with probability of at least 1 UN7y
G² Rd log(NT) max (1.8
13 13 C(Jmax: 7)
(13)
VNT
This shows that the error of estimating 3* (in squared euclidean norm) scales as δ
VNT
Tm
p log (T).
GR
log(
ForErr(N,T) < CO ROP'' 2
T
min (a2 (o2,o4)
VNT
where e. is an absolute constant, and Ò denotes a constant that depends only (polynomially) on fmax; 71 Amiul (V) Amax (T), and Omnax-
This theorem establishes that the forecasting enror for the next T time steps scales as Ô ( /NT with high probability. Note that the Ô (+) term function of T only, the error incurred when estimate each Ogg- Recall that we estimate an separately for [N] using the available (perturbed) T observation of each process. The Ô VNT hand is a consequence of the error incurred when learning and forecasting f1...,,a, which is done collectively across the N time series. Finally, Theorem 4.4 implies that when = the
consequence of
term on the other
forecasting error scales Ô (+). The proof of Theorems 4.3 and 4.4 are in Appendix and respectively.
10The Ö(-) notation is analogous the standard O(-) while ignoring log dependencies.

---

Table Performance of algorithms on various datasets, measured by mean n².
Traffic
Electricity
Exchange
Synthetic
SAMoSSA
0.776
0.829
0.731
0.476
mSSA
0.747
0.605
0.674
0.366
ARIMA
0.723
<-10
0.756
0.305
Prophet
0,462
0.197
<-10
-0,445
DeepAR
0.824
0.764
0.579
0.323
LSTM
0.821
-1.261
-1,825
0.381
Experiments
In this section, support theoretical results through several experiments using synthetic and real-world data. In particular, we draw the following conclusions:
1. Our results align with numerical simulations concerning the estimation of non-stationary compo- nents under AR stationary noise and the accuracy of AR parameter estimation (Section 5.T]. Modeling and learning the autoregressive process in SAMoSSA led to consistent improvement over mSSA The improvements range from 5% to 37% across standard datasets (Section 5.2)
2.
5.1 Model Estimation
Setup. We generate a synthetic multivariate time series (N 10) that is mixture of both harmonics and an AR noise process. The AR process is stationary, and its parameter is chosen such that X m Aryi is one of three values {0.3, 0.6, 0.95]. Refer to Appendix B. for more details about the generating process. We then evaluate the estimation error EstErr(N,T, T, and the AR parameter estimation error ||a ã1 llz as we increase T from 200 to 500000.
shows Results. Figure 2a visualizes the mean squared estimation error of f1, while Figure sho the AR parameter estimation error. The solid lines in both figures indicate the mean across ten trials, whereas the shaded areas cover the minimum and maximum error across the ten trials. We find that, as the theory suggests, the estimation error for both the non-stationary component and the AR parameter decay to zero as NT increases. We also see that the estimation error is inversely proportional to (1 )*), which is being reflected in the AR parameter estimation error well.
10
EstErt( T,1)
0.3
A
10-1
0.6
A
0.95
X 0.95
10
0.5
10
0
102
10
10
10a
VNT
/ NT
(a)
(b)
Figure 2: The error in SAMoSSA's estimation of the non-stationary components and the AR parameters decays to zero as NT increases as the theorem suggests.
5.2
Forecasting
We showcase SAMoSSA's forecasting performance relative to standard algorithms. Notably, we compare SAMoSSA to the mSSA variant in to highlight the value of learning the autoregressive process.
Setup. The forecasting ability of SAMoSSA was compared against (i) mSSA B. (ii) ARIMA, very popular classical time series prediction algorithm, (iii) Prophet (28), (iv) DeepAR 23 and (v) LSTM [14jon three real-life datasets and one synthetic dataset containing both deterministic trend/seasonality and stationary noise process (see details in Appendix B}. Each dataset was split into train, validation, and test sets (see Appendix B. Each dataset has multiple time series,

---

aggregate the performance of each algorithm using the mean R² score. We use the R score since it invariant to scaling and gives reasonable baseline: a negative value indicates performance inferior to simply predicting the mean.
Results. In Table 1 we report the mean R for each method on each dataset. We highlight that, across all datasets, SAMoSSA consistently performs the best or is otherwise very competitive with the best. The results also underscore the significance of modeling and learning the autoregressive process in real-world datasets. Notably, learning the autoregressive model in SAMoSSA consistently led to an improvement over mSSA, with the increase in R values ranging from 5% to 37%. We note that this improvement is due to the fact that mSSA, as described in [3], overlooks any potential structure in the stochastic processes 21 z1(-)... zN (+) and assumes .ii.dmmen--a. zero noise process. While in SAMoSSA, we attempt to capture the structure of z (-)....... Zx(( through the learned AR process.
6 Discussion and Limitations
We presented SAMoSSA, a two-stage procedure that effectively handles mixtures of deterministic non- stationary and stationary AR processes with minimal model assumptions. We analyze SAMoSSA's ability to estimate non-stationary components under stationary AR noise, the error rate of AR system identification via OLS under observation errors, and a finite-sample forecast error analysis.
We note that our results can be readily adapted to accommodate (i) approximate low-rank settings (as [3]); and (ii) scenarios with incomplete data. We do not discuss these
in the model by Agarwal settings to focus on our core contributions, but they represent valuable directions for future studies.
Our analysis reveals some limitations, providing avenues for future research. One limitation of our model is that it only considers stationary stochastic processes. Consequently, processes with non-stationary stochastic trends, such a random walk, are not incorporated. Investigating the inclusion of such models and their interplay with the SSA literature is a worthy direction for future work. Second, our model assume non-interaction between the N stationary processes 21,... ........... NN. Yet, it might be plausible to posit the existence of interactions among them. possibly through a vector AR model (VAR). Examining this setting represents another compelling direction for future work.
References
[1] Y. Abbasi-yadkori, Pál, and C. Szepesvári. Improved algorithms for linear stochastic bandits. J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
[2] A. Agarwal, A. Alomar, and D. Shah. tspdb: Time series predict db. In NeurIPS 2020 Competition and Demonstration Track. pages 27-56. PMLR, 2021.
[3] A. Agarwal, A. Alomar, and D. Shah. On multivariate singular spectrum analysis and its variants. ACM SIGMETRICS Performance Evaluation Review, 50(1):79 80, 2022.
[4] A. Agarwal, M. J. Amjad, D. Shah, and D. Shen. Model agnostic time series analysis via matrix estimation. Proceedings of the ACM on Measurement and Analysis of Computing Systems, 2(3):1-39, 2018.
[5] A. Agarwal, D. Shah, and D. Shen. On principal component regression in a high-dimensional error-in-variables setting. arXiv preprint (xXiv:2010.14449) 2020.
[6] A. Alanqary, A. Alomar, and D. Shah. Change point detection via multivariate singular spectrum analysis. Advances in Neural Information Processing Systems, 34:23218 -22200, 2021.
[7] A. Alexandrov, K. Benidis, M. Bohlke-Schneider, V. Flunkert, J. Gasthaus, T. Januschowski, D. C. Maddix, S. Rangapuram, D. Salinas. J. Schulz, al. Gluonts: Probabilistic time series
models in python. arXiv preprint arXiv:1906.05266 2019.
[8] T. W. Anderson. The statistical analysis of time series. John Wiley & Sons, 2011.
[9] Chollet. keras. httpps//ggtaab..cc///callleee//e February 2020.
2015. Online; accessed 25
10

---

[10] C. Davis and W. M. Kahan. The rotation of eigenvectors by perturbation. iii. SIAM Journal on Numerical Analysis, 7(1):1-46, 1970.
[11] Durbin. Estimation of parameters in time-series regression models. Journal of the royal
statistical society; Series B (Methodologiccl)) 22():139-133, 1960.
[12] Facebook. Prophet. httpes///gtaul.ccco//c/aabbppppp 2020.
[13] M. Gavish and D. L. Donoho. The optimal hard threshold for singular values is 4/v3. IEEE Transactions Information Theory, (0(8):5040-5053 2014.
[14] A. Gers, Schmidhuber, and Cummins. Learning forget: Continual prediction with Istm. 1999.
[15] N. Golyandina, A. Korobeynikov, and A. Zhigljavsky. Singular spectrum analysis with R. Springer, 2018
[16] N. Golyandina, V. Nekrutkin, and A. A. Zhigljavsky. Analysis of time series structure: SSA and related techniques. Chapman and Hal//CRC, 2001.
[17] R. A. González and C. R. Rojas. A finite-sample deviation bound for stable autoregressive processes. In Learning for Dynamics and Control, pages 191-200. PMLR, 2020.
[18] H. Hassani and R. Mahmoudvand. Multivariate singular spectrum analysis: A general view and new vector forecasting approach. International Journal of Energy and Statistics, 1(01):55-83, 2013.
[19] Y. Jedra and A. Proutiere. Sample complexity lower bounds for linear system identification. In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 2676-2681. IERE, 2019.
[20] Y. Jedra and A. Proutiere. Finite-time identification of stable linear systems optimality of the least-squares estimator. In 2020 59th IEEE Conference on Decision and Control (CDC), pages 996-1001. IEEE, 2020.
[21] Y. Jedra and A. Proutiere. Finite-time identification of linear systems: Fundamental limits and optimal algorithms. IEEE Transactions on Automatic Control, 2022.
[22] Microsoft.
nimbusml package, ssa class, 2020. https///ocs.miicrooff..oom//
[23] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 2019.
[24] T. Sarkar and A. Rakhlin. Near optimal finite time identification of arbitrary linear dynamical systems. In International Conference on Machine Learning, pages 5610-5618. PMLR, 2019.
[25] S. Seabold and Perktold. statsmodels: Econometric and statistical modeling with python. In 9th Python in Science Conference, 2010.
[26] R. H. Shumway, D. S. Stoffer, and D. S. Stoffer. Time series analysis and its applications, volume 3. Springer, 2000.
[27] M. Simehowitz, H. Mania, S. Tu, M. I. Jordan, and B. Recht. Learning without mixing: Towards sharp analysis of linear system identification, 2018.
[28] S. J. Taylor and B. Letham. Forecasting at scale. The American Statistician, 72(1):37-45, 2018.
[29] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint arXiv: 2010.
[30] R. Vershynin. High dimensiiins probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
[31] M. Wainwright. High-dimensional statistics: A non-asymptotic viewpoint, volume 48.
Cambridge university press, 2019.

---

[32] P.-A. Wedin. Perturbation bounds in connection with singular value decomposition. BIT Numerical Mathematics, 12(1):99-111 1972.
[33] P.-A. Wedin. Perturbation theory for pseudo-inverses. BIT Numerical Mathematics, 13:217-232.
1973.
12

---

A Concentration inequalities for AR processes
In this section, we present important lemmas that central to the main results. Before stating the two lemmas, we start with key important quantities of AR Processes.
A.1 Controllability Gramian of AR Processes
In this section, we explain the notation of key quantities used in the main results. First, note that the process in (2 can be written in matrix-vectie -vector form by defining transition matrix and input vector
T ml:m(pp- -1)
anpp
A,
Bn 0p--1
pp.-1
0p.-1
where -1 the column vector of Pn zeroes. Then, with (t) (tn,t))...,; In Pa +1);
x X, + BRdne-
With this setup, we define two important matrices,
U = Σ B, 1=
(14)
where ,, is known as the controllability Gramian of the n-th AR(p,) process In. For a positive semi-definite matrix M, let (M) 2 0 and Amin (M) 0 be its maximum and minimum eigenvalues, respectively. We define Amax ()) = max, Amax ()))) and Amin (V) = min, Amiu n
We also define Amax (T) and Amin (T') analogously for r.
A.2
Key Lemmas
Now, we present two important lemmas that are central to the main results. The first lemma states that a vector of consecutive entries of stationary autoregressive process is a sub-gaussian vector with variance proxy ozz
Lemma A.1. Let assumptions 2.3 and 2.4 hold, then, for any E [N]. t 2 1,K 2 1, the vector (t). + K)] sub- nub--aussian vector with variance proxy of
This lemma, along with the well known bound on the norm of sub-gaussian vectors, is key to the bounds presented in Section 4.1 Further, this Lemma. along with E-net argument, allows us to bound the operator norm of the stacked Page matrix of the AR processes 21 (f),...., IN (t). Specifically, we establish next that its operator norm grows as O(VNT/L) with high probability.
Lemma A.2. For any T > 1, 1 L S VNT, and to > 0, let Z(x2, T,t,)
(to) = [ZII,,
Z(In, to))
be the NT L stacked Page matrix of the AR processes x1 Assumptions 2.3 and 2.4 hold. Then, each row (ta)); is a sub-gaussian vector with variance proxy o. Further, With probability exp (-eNT/L)
as defined 2 Let
12 (10)1l2 200 V NT/L.
Next, we prove Lemma A. and Lemma A.2 Before stating the proofs, we start with the important helper Lemma A.3
Lemma A.3 (MA(oc) representation of an AR(p) process). Let xxt) be a stationary AR(p) process defined as
p xxt)) = -)altt (yt).
Assume that the roots (Aı,..., Xp) of its characteristic polynomial g(z) 20 Ei-1 0:20-1 distinct. Then, xo() has the MA(20) representation
20
x(t)= Σ
Ban(t k).
2 -1 (1-
with 3, = Ei- a,2 a; =====1sspp.. 1#1
13

---

Proof. First, note that that we can write xxt) in lag operator form
1
P
J(L) x(1)
Σ
a
x(() ntt).
(15)
Here, L is a lag operator such that (t( srll 1). Note that since we have assumed the AR process to be stationary, all the roots of g(2) lie inside the unit circle [26]. Using the roots of g(2). 135) can be written as,
II - x() n(t).
Further, using the partial fraction decomposition we can write the process as
Zt
I i= X; (t))
P
Σ
1 X₂L (t)
where,
1 A3
II
0,
A
1515p,
Recall that we can expand each fraction as
a2
Σ
L
A
K=U
Thus, finally, we can write r(t) as,
xxt) = Σ
Bun(t k),
k=0
with 3g = Dd aa
A.3 Proof of Lemma A.1
Proof. First, let (t, K) [x, (t), ... F. (T K)] be the vector of interest. Let Ex(t, K)
(tm(t)...... ...).................. Tn (T + K)]. then using Lemma A.3 we can write X,, (t, K) as
(t, = Σ Et (t-k,
k=
Note that Ban Ft kc K)) clearly sub-gaussian vector with variance proxy 02 Further note that the sum of two sub- sub-gaussian random vectors, not necessarily independent, is also sub-gan -aobsaussian random vector. To see this, suppose x and y are both sub-gaussian vectors with variance proxies o² and og. Then for any E SK
v (x vT vTy
subG ((o, + o,))))
Thus, X, (1, K), for any n, is also sub- aub-gaussian random vector with variance proxy
oc:
oc:
2
Σ
1sntlor
S 0²
Σ
max |3ntk]
14

---

Now we are interested in bounding Σk=o k== max, 15 To do so, first recall that Buk =
Pre Eis
ane Ans
where And C, i E [pm] denote the roots of (n(2) and a₁
Hisss tx (1-
Let
Ana
A. = maxi,n Anal and CA.
-E
PM II TisSSP. j#i
Ani Ani
Then, we can bound Brik as
Pu
18nk S X
Σ
ass
==1
S Cx Xx
Thus, with cy = max, CAn
C2
max 1Bukk
(16)
1- X
A=
630²
Hence, (I, K) a sub-gaussian random vector with variance proxy of = (1-A.)
A.4 Proof of Lemmal A.2
Proof. Let = NT/L. We will find upper bound the operator norm by following steps: (i) we will represent each z. (t) Vn E [N] as an MA(co) process; then, (ii) we will bound the operator norm of the stacked Page matrix of this MA(oc) processes.
Step 1: MA representation for AR(p).
As direct application of Lemma A.3 write xn (t) as,
(x = Σ
Basiia (t k:),
k=
with 3nk El=1 Ane Oni IIssis.. iti
Anz Anne
Step Bound the Page matrix for the MA()c) Page matrix.
First, let's consider, without loss of generality, the case when to = 1. That is, let's consider Z₁ Zzz(1). (1). To get the desired bound, we will first obtain high probability bound on the Euclidean norm of Z v for any vector E S = {v E RA: ||vl2 = 1}. Then, we will bound the operator norm using this and an s-net argument.
Note that by step (1), we have
Fa (1) = Σ Baknn (t k))
k=0
Now define E; as the analogous stacked Page matrix for 7), (t). (t) for t E (i....,++7" 1}. Observe that E; are simply identically distributed (but not independent) matrices i.. sub-gaussian random variables. Then we can write
Z,
E k=0
B
(17)
where Bg E R8x0 is a diagonal matrix with the i-th diagonal [B] = Br Then
|
Σ
B
2
15

---

It is easy to verify that E v is a sub- sub-gaassiin random vector (see Definition H.3 with variance proxy o², because for any u 8 we have (using E shorthand for E;)
L
T Eu
Bay Eij
1= f=1
4
subG
ΣΣ i=1
of o2
sub
Using the above, B; v is smbG (Bnk|20²) Since Z is of sub-gaussian random vectors, it is also a bub--gassiinn random vector with variance proxy of as we showed in the proof of Lemma A.1
Using Lemma H.2 we get, with probability least 1 6,
Z v12 40 20 log
Equivalently, for 2 0, holds that
t²
P P[I 40, v 1) S exp
40
Now we have shown that for any fixed S the quantity |Z |Zalll grows as O ( (q) with high probability. What remains is to apply the union bound over an e-net to extend this to the maximum over St i.c. the operator norm of Z..
Let Σ be a maximal 1/2-net of S i.e. a set of points S spaced at least 1/2 from each other such that other point be added without violating this property. Now consider v, such that [Zz 1- Since Σ is maximal, there is some v E Σ within 1/2 of v.. Since
[ v=
2 (v v. v2)112 < IZzlla ||v
vx S Zz /2,
and
(v- v.))l2 2 1 v. 112 -
IZ
we have that
| v/l
Z /2,
i.e taking the maximum over E is equivalent to doing over - up factor of two. Using the fact that Σ < Cd for some absolute constant C, the fact that < q, and using a sufficiently large absolute constant > 0 have
(|| 1rz 200 S Σ P(
ca
VEE
< 20 exp (-d')
exp (pc"q))
where c" > 0 and > 0 are absolute constants. Recalling that NT/L concludes the proof for the operator norm.
Finally, we show that each row [Zz(4)); the stacked page matrix a sub-gaussian vector with variance proxy 02. Consider, without loss of generality, the first row [Zz]1.- Using (17). have for any v E S
v Z
Σ v EII k=
note that v [E subG 22) and hence, using the same argument for the
sum of dependent sub-gaussian random variables used in the proof above, have v T [Zz]1- subG(o?)..
16

---

B Experiment details
In Appendix B.I we detail the datasets used. Appendix B.2 we explain the implementations and hyperparameter choices for our method and benchmark algonttmms. In Appendix [B.] we describe our hyperparameter tuning strategy.
B.1 Datasets
In the model estimation experiments, generate synthetic multivariate time series (N 10) that a mixture of both harmonics and an AR(2) process. The harmonics are generated as follow: we first generate R fundamental time series of the form (t) sin (wat + xxx) for & E [R]. where we and Φk are uniformly randomly sampled on the ranges [2Tr/ / 100, 2 /50], [0, 2T] respectively. i= Ct.gg(t)) (t) for A E [N], where the mixture
We then sample the N 10 series using fa(t) = - coefficients tjj are independent standard normal random variables. The AR process is stationary, and its parameter is chosen such that X maxn,1 Avail is set to one of {0.3, 0.6, 0.95). The noise of the process has variance o² = 0.2.
In the forecasting experiments, we evaluate method and benchmark algorithms on three datasets and one synthetic dataset. The preprocessing steps and setup of each dataset are described below.
Traffic Dataset. This public dataset obtained from the UCI repository shows the occupancy rate of traffic lanes in San Francisco. The data is sampled every 15 minutes but to be consistent with previous work, we aggregate the data into hourly data and use the first 10248 time-poine for training, the next 48 points for validation, and another 48 points for testing in the forecasting experiments. Specifically, in our testing period, we do h-hour ahead forecasts for the next 48 hours.
Electricity Dataset. This is public dataset obtained from the UCI repository which shows the 15-minutes electricity load of 370 households. We aggregate the data into hourly intervals and use the first 25824 time-poine for training, the next 48 points for validation, and another 48 points for testing in the forecasting experiments. Specifically, in our testing period, we do I-hour ahead forecasts for the next 48 hours.
Exchange Dataset. This is a dataset containing the daily exchange rates of eight foreign currencies, including those of Australia, the UK, Canada, Switzerland, China, Japan, New Zealand and Singapore, between 1990 and 2016. We standardize each time series to have zero mean and unit variance since the range of typical exchange rates varies greatly across the currencies. We use the first 7528 time- points for training, the next 30 for validation, and the final 30 for testing. All forecasts are made 1 --day ahead.
Synthetic Dataset. We generate R fundamental time series of the form 9N (t) sin (wyt Φx) ) m for k E [AR], where Whe dx, and my are uniformly randomly sampled on the ranges [2T/100, 2T 100, [0, 2m.. and [-5 x 10- ,,5 10-4] respectively. We then sample N mixtures of the form (t)
ER= Ck:gk(t) (t) for E [N]. where the mixture coefficients Okt are independent standard normals. We then inject AR(1) autoregressive noise to each series to produce yu (t) (t) + 2x(t) for k E [N]. where (t) (a(()))eEN [E[N] are independent processes, each with fixed autoregressive coefficient a = -0.5
and noise variance a² = 1. We fix the coefficient to ensure that each time series has significantly autocorrelated noise that can be exploited to improve predictions.
Table 2: Dataset and training/validationnnn split details.
Dataset
No. time series
Training period
Validation period
Test period
Traffic
963
to 10248
10249 10296
10296 to 10344
Electricity
370
to 25824
25825 to 25872
25872 to 25919
Exchange
to 7528
7529 to 7558
7559 to 7588
Synthetic
25
1 to 10000
10001 to 10025
10026 to 10050
17

---

B.2 Algorithms
In this section, we discuss the implementations used and hyperparameter tuning details for each algorithm that was evaluated.
SAMoSSA & mSSA. Since mSSA is special case of SAMoSSA where the second stage fitting autoregressive models - is skipped, we use one in-house implementation to evaluate both algorithms. The relevant hyperparameters
follows:
1. The number of retained singular values, k. This hyperparameter is relevant for both algorithms. We select k using one of three methods: (i) the data-driven procedure suggested in [33 which computes a threshold based on the shape and median singular value of the page matrix; (ii) picking k as the minimum number of singular values required to capture 90% of the spectral energy of the page matrix; (iii) picking fixed low rank, specifically = 5.
2. The shape parameter of the stacked Page matrix. This is the ratio between the number of columns and rows of the stacked Page matrix used to carry out the matrix estimation procedure. While our theoretical analysis fixes this value to 1 (recall that we set L = /NT), in practice we have found that a slightly wider Page matrix can improve performance. This parameter is chosen from 31)
3. The number of autoregressive lag coefficients, p. This hyperparameter is only relevant for SAMoSSA. While each time series does not necessarily need to use the same value, we enforce this to be the case for computational efficiency. We choose p E {0, 1,2,3) 3], i.e., the model is allowed not fit residual model if it brings no apparent benefit, corresponding to the case = 0.
Prophet. We use Prophet's Python library with the parameters selected using a grid search of the following parameters as suggested in [12]
1. Changepoint prior scale. This parameter determines how much the trend changes the detected trend changepoints. We choose this parameter from {0.001, 0.05, 0.2].
2. Seasonality prior scale. This parameter controls the magnitude of the seasonality. We choose this parameter from {0.01, 10).
3. Seasonality Mode. We choose between an "additive" and "multiplicative" seasonality term.
The parameters for each time series are chosen independently.
ARIMA. We used the ARIMA implementation of the Python library statsmodels [25]. The hyperparameters are grid searched independently for each time series, and are follows:
Autoregressive order. This is the number of autoregressive lag coefficients fitted, chosen from {1,2, 3}.
2. Differencing order: This denotes the number of times the time series is differenced before model is fitted, and we choose between O (no differencing) and 1.
3. Moving average order. This the number of moving average lag coefficients, chosen from (1,2,3)
LSTM. We use Keras implementation 9 We perform a grid search the number of layers (2, 4).
DeepAR. We the implementation provided by the GluonTS package We the default
parameters.
B.3 Parameter Selection
We tune hyperparameters for all algorithms evaluated using rolling cross validation. During the validation step, each model is fitted only once on the training set; then, repeatedly makes one-step- ahead forecasts and is afterwards provided with the realized value. in order to use the updated data to make the next forecast. During test time, the model is fitted on the training and validation sets, and similarly repeatedly makes one-step-ahead forecasts with the fitted parameters. For the models that choose one set of hyperrparameters for all N time series, the mean R² score over all time series on the validation set predictions is the metric of choice.
18

---

C Proof of Theorem 4.1
While the analysis generally follows the argument for the imputation error in Bl we adapt it for our different assumptions and metric of interest. Metric. First, recall that our metric of interest is
In this section, provide the proof for Theorem 4.1
T T Σ (1.00 fn()))
EstErr(N, n,
(18)
Recall that 2 is the stacked Page matrix of fI,..., JN. and Zy is the stacked Page matrix of fi.....fr: To bound this error, we first establish a deterministic bound through a general lemma for Hard Singular Value Thresholding (HSVT).
C.1 Deterministic Bound
We state the following result,
general version of which is stated in 3
Lemma C.1 (.3), For k 2 1, let Y = M + E E Raxp with (an((M) k. Let UR and U En V denote the top singular components of the SVD of M and Y respectively. Then, the
HSVT estimate M Ú Σx V is such that for all j E [4].
MM, M, M < 2
EIS
(|
2 + [m²
+2 Vo T E.J
(M)
Proof. First, note that
M, M.. V V.J. V.V[M.. + (V.V - -M,).
where M, T a subspace orthogonal to V TT and hence we have by the Pythagorean theorem,
by definition. Note that the vector (v M.. - -M,) is in the span of
M
M
2 2
V. V.
V.VM V...M MT
V.
V MM
(19)
The first term can be futher decomposed as
V V.V.E 22
V V.
Vv Vv
E
+
V& E
<2
V.V.
E
+2
VAVA
(20)
Next, we bound the first term on the right hand side of (0 To that end, by Wedin sin Theorem (see 10 32 and recalling rank(M) k
E
VAVE VAVN
(M)
Then it follows that
V.V. E
2
EE
E 1
+ V2 V ET
V2
(21)
(M)2
Then, we turn to bounding the second term from
19
V.VAM - MJ.
V.VI V&Va M
V T
V.V.
M
||
M
(22)
ox (M
19

---

Using 19
21
and
22
we have,
M,
M
E ok M
(1=771²
MM 2) +2 VaV& E
(23)
which completes the proof.
Now, Let's apply Lemma C.I to our setting. Let zr be the stacked page matrix for 21 (f)...................... In (t).
Using Lemma C.1 and setting Y = ZT MM = z E = z and reusing U₁ En V to represent the top k singular components of the SVD Z,, results in the following deterministic bound for HSVT with rank set to k.
Z 0h (2z)²
| +2 U.U T [2mlo]
(24)
C.2 Deterministic To High-Probability
Next, convert the bound in 24 to bound in expectation (as well one in high-probability) for [2, 11112,00 In particular, establish
Theorem C.1. Let assumptions 2.7
2.2
2.3
and
2.4 hold. Let On [I-X.)) where X
Ani Ani 711
maxi, DAnil and = max, E
Pa
I1 TisisP itf
Then, the HSVT estimate Z
i=
with parameter l such that
(Z)) (02 + fmax) + (o
CNT(2
max jE
L
[1), - [2)
L
(25)
with probability of at least 1 (NT)
Proof. We start by identifying certain high probability events. Subsequently, using these events and we will conclude the proof.
High Probability Events. Let q = /T the number of columns in the stacked page matrix 2, and let C >> 0 be some positive absolute constant, Define
E1 /Z < CoLva).
Eg l0 2: Zz 12,00 Co.va).
Eg
max jEly
U T [Z.] S Cof Cooklog(o)
Lemma C.2. For some positive constant 01 and Eg,
0 and large enough definitions of Er, Ez,
P F1 > 2c 2019
PP E2) > -19
PP E3 > (NT))
Proof. We bound the probability of events above below.
Bounding E,, This is an immediate consequence of Lemma A.2
Bounding Eg. Recall that we assume L < q. Observe that for any matrix A E R4xg,
|| ||
Al|2,50 S || Thus using the argument to bound E,, concludes the proof.
20

---

Bounding E3. Let ui...., Ug be the orthonormal basis for Ug. Then, consider for j E [g].
UR [z]
[Zn]
Z
i=1
i=1
i=1
where Zi u: [Z j.
As shown in Lemma A.1 [Zn] is sub-gaussian vector with variance proxy we have
o.. Using Lemma H.4
k
Σ
2 11 - cexp 16ko
Therefore, for choice of = log()) with large enough constant C' > 368, we have,
Σ
Co log(1))
C1
923
1=1
Recalling that L S q, and taking a union bound over all j E [9]. we have that
Di I (qL)
The following is an immediate corollary of the above stated bounds.
Corollary C.1. Let E == Ei E2 Eg. Then,
C
P(E)) S (NT)IT*
(26)
where C₁ an absolute positive constant.
Thus, under event E, and using 24 we have, with probability 1 (NV)
Cofy (noz + a
[ [2)) s
max jE(g]
(Z,)²
Co log(4)
Recall that q NT/ then,
+ - 5E(g
CO [Z,,((,)) (0" + (2 + Co looklog(NT/L
max
(27)
I
This completes the proof of Theorem let P, = [M-1)T + 1,..., af
Finally, now we are ready bound EstErr(N, T). First, be the set of columns in the stacked Page matrix Z, that belongs to the n-th time series ft). Note that
C.1
1 (t.))) - fr(())² - Σ
2 /zyl.//
jen.
CNT
(02 + f log(N/L))
T JER.
(Z,)
CNT(N
(02 + Plax)) + (og
Lon (2
Choosing L = VNT. and using oz (Z,)) N N2TT
we get.
(,()))))((()3)---)) C1 2 (og(NT)
T
VNT
Finally and setting & = RG concludes the proof.
21

---

D Proof of Theorem 4.2
In this section, we prove Theorem 4.2 which bounds the estimation error l|ā - anlll Vn [N]. First, recall the OLS estimate 0 [nn,1 2.. Gn.pl defined as
ā argmin
Σ -
(28)
where = ... (l p + 1)]. Alternatively, the OLS estimate can be defined as,
a = X xn) X,, Yn,
29)
where X, E 1R(T-P))xP is the row-wise concatenation of X (t) fort e (p.p+1..... T 1}. and
Pn = (p + 1)..... (T) Further, consider the OLS estimate with access to the true AR series (.).). Precisely, let an defined as
a = (x xn) X Y..
(30)
where X, E R(''-)xp is the row-wise concatenation of (t) = [=n((),...,, In p + 1)] for LE (p,p+ 1,.., 1}. and Y, = [z,, (p +1),+ II,...,,,,(T)) (T)]-
In what follows, as we do the analysis for one process In, we drop the subscript n. We note that the analysis holds for all n [N]. First, note that,
|| - aa22 < 2 f|ā - al12 + f|a al12
(31)
To bound ||ā arliz we use Corollary E.2 which states that with probability 1 -mT.
p²
log T'max VV)
aā- 1 < C
(32)
TAmin V
Amiu (V)
M denote the Moore-Penrose inverse of matrix M. Further, let 5(t) r(t) a(t) f(t) f(), and A X X.i.e. A is the row-wise
Next, to bound ||a | - ālp let M = (M
concatenation of (1) = ((()) ...., 5(( p 1)] for 1 {p.p 1,...,?
T 1}. Then, with
Ay = [6 + 1)....,de ((T)) we have
X x x'rl|
1r/12
/|2
||a al|2
X Ay + +xxt X Y
2
X + (xl - x )r/
X
| l12 + x X
1/, ||x|||
+2 max {||x
2
x 11²
|| 11ą
11
| 112 |1
(33)
22"
Where Lemma H.6 is used in the last inequality. Note that,
x |-=
inf
| = V Amii (xTx x)
(34)
2
vERr:llo|,
inf
X
DERP:||ul1
inf DEB::||u||===
|1(X A A)ull2
V Amin (X T X) - | 1
(35)
22

---

Using 33
34
35
and the theorem conditions, we get
2 A |
|ā 2 S
+
Amin
X X)
-111112
min
Amin (X X),
(v Vimin (X X)
||ا||
(36)
Using Corollary E.I we get that with probability 7
Amin ( (X (xx)= 2 0T - Amin VV).
(37)
Now, note that | Et=p+ p+1 (F(1) /())) and
| 1 < S
Σ E=
(f(1) S()))
(38)
- Recall that (ft)- f(t)² = TEstErr(N T, n), which we will refer to herein as EstErr(N, T)
as we drop the dependence n. Therefore, we have,
A 2 S pl EstErr((N
(39)
A 11 S TEstEr((N,TT)
(40)
Finally, using Lemma H.4
with probability at least Tir it holds for an absolute constant C' > 4,
|| s Co T
(41)
First, note that, using 37 and 39
we have
11 (T Amiu (V))
Amia (X X
PTEttEm((NNT
V2
2 Co CoVCTTi T>mmn (V),
(42)
for sufficiently small EstErr(N, T) such that EstErr(N, < 6p
e'Amin(
Now using 36).
37
44
and (42) we have,
CoNVDEStEr)(N, log(?)
l| 01l2 5 (onfoo2mmnn (1) Amiu (V)
Therefore,
||ā -ll2 < c D (T) o² Amin (V) x2
EstErr(N, T) max {
1
(43)
admin
finally, using 31
and 43 we get with probability
Tii
l| all2 s c.P Plog(T T o2 Amin (V) o²
1
Cp² T'mmn (V)
T Thmaxe
EstErr(N, TT) max
+
log
o² Amin (11
Amin V/)
(44)
which completes the proof of Theorem 4.2
23

---

E Identification of AR processes
In this section, state key theorems for the identification of AR processes without the perturbation we consider in this paper. The theorems in this section and their proofs follow closely the ones in [77, we modify them accommodate general sub-gaussian noise and state them fully for completeness.
First, consider the stationary AR(p) process defined by the transition matrix
A
a 1:p pp-1
Ap
(45)
Op-1
where 0p-1 is the column vector 1 zeroes, and input vector
B 1
(46)
0p-1
Then, with ((1) = (x((), .... llt p + 1)], the AR process is described as
(t) = Ax + Byp-
(47)
With being sub-gaussian random variable with variance proxy o² With this setup, we define two important matrices,
r
A²
(48)
[=0
ABBB (A
(49)
where is known as the controllability Gramian. Theorem E.1. Consider the AR(p) process described in (47 {nd } is i.i.d. mean-zero sub- gaussian noise with variance o². Given 0 < E < 1. define the following quantities:
V2 = o² o² - p) (T - -()
Vu = x² (T p) ( + el')
8(c, T') == cexp
P
To(c) = 752 log(p)²
max (1, [|al]]]
o²²²
where and are constants that depends only on the process coefficients a1,..., Op. Then, for T> Tole) and valid values of e such that Vz > we have
P (V, = X x = V.) 2 1 - 8(e,T).
Remark E.1. To get Ve O, must choose sufficientll small such that ( eI') > 0. Using Weyl's 's inequality, we have
Amin (T cT) 2 Amin (V) edmaax
(50)
This suggest that < Ve 2 OT - p) Amin I. Note that this ensures < 1/2, since e < Aman and > W as can be seen from
(1) Amnaze (T) yields Ve x For example, a choice of =
2Am [T' will yield
Annin (W)
Amaax (T)
80
2 (T-W)==
Σ
27 - AB A
Σ |1(AT) (AT)))²
Σ
)211
|
(A )212
0
0.
24

---

Proof. First, write
T-1 x E X,X,
Expanding the summand using the process dynamics yields
X (AX, Bm (AX, + B
=AX,X, A + 1 + BX1BXX
AT + - BB
Now, defining Vr X (/(I-pp) we can sum the expression above - p- 1, 1,... 210 get
T-2
V AV T-p A(X X A +
Σ (741AX, + + A + (B. t=p-1
This is called a Lyapunov equation, and due to the stability of the process it has solution V7 Et=o A Er(AT)² The key step is argue that Ex converges R² BB" with high probability, which will allow us show that Vr converges to ouy To do this, we will split the terms in Er into three groups, handling them separately. For some c > 0, we will bound the probability of the following events:
X--17--1117 -o & = {o Et-p1 (1et1 BB - -B² S ST
& = {p (A(Xp-1XX)
-p)/33
Es = {o (Etep T M BX, AT) 3 c0² -p)/3
where (-) denotes the spectral radius of matrix.
Lemma E.1. With the setup above, it holds that
602 VT
P ((E) cccep cexp
303p(11a112 +
for some absolute constant 0.
Proof. Note that
p (A(X X - -
S cT (022TT-P))) p)/3
pAA(Xp-188 p- X ()A) S co²2 p3/3
||Ax--1 11/ ||AX 1/2 5 60 2 - 3
602 IT pp)
1Xp-1 122 + |
3 ||
(T p) 11Xp-1 12 ||
|||A|1
so, considering the fact that X,-1 and Xx are identically distributed by stationarity, suffices to upper bound the probability of
(T)
8'= -'=11 11p-11111 6 ||ll
a²g² First, recall that Xp-1 is a sub-uuussiaa vector with variance proxy (1-X-> = =2 (see Lemma A.2
Thus, using Lemma H.2 for
25

---

co² -DT-DD
co² - p)
P 1
S cexp
6 61|A|12
960 |A
co²vT
S cexp
30 ||l|
where in the last inequality, we used the assumption T To(c) 752p². Finally, note that
||A|| I Tr(A) = ||al2 + p- p(laa|2 2 + 1). Recalling that P (E1) > - 2P (8") concludes proof.
Lemma E.2. Continuing with the setup above, it holds that
-
((2)) cexp
48
Proof. Since the spectral radius of BB T is unity, &g is equivalent to the complement of the event
E' ='= (T (D- -top-1 ( - -0²) > /000/33
Since 7 sub-gaussian with variance proxy 0², n o² sub-exponential with parameter 16g2. Thus, Lemma H.5 applies, and therefore,
-EVT p
PP (E') S cexp
(51)
48
Lemma F.3. With the same setup, we have
eT²VT
P ((E)) - cexp
ccecp
75po
(2 ||(1112
as long as
75 7520"p² log(p)2
(2 laa1))
T Tole)
ote²
Proof. Recall the definition
Eg E=== (EE-- B S eo² -p)/3]
(XTTXXTAT)
For the sake of variational analysis of the spectral radius, consider any vector q [q1 qT R such that ||1|| = For notational simplicity. Then
T ET-p-
1 AX + B 9 S 2|
M+MT
Next, we aim to bound ||Mall First note that only the first column of M nonzero which evaluates
ET- TX,
1
M. M1+11
(52)
Etsp1 st-pt1
Therefore, have Mq q1 M. which implies that for any unit vector g.
q TM + M g s
M2 1
26

---

Therefore, we will bound P (Eg) by
P
P (EE)
P
Mai 2 e p)2
P
60 -T
-p) p)
IP
|Mi111
(53)
6vp
i=1
First
note that
each
term
P (Ma 2
50' 6/8
can
be
upper
bounded by
21P EE 1m 2T
(T-p)) 6
where 21 and
24 ei-1 for i
where &, is the
i-th standard basis vector. Using Lemma H.7 and choosing Ze 2, X and W, Th, have that
P ({ETEt me Xi-1 2 x n E (2) x2-1)² ===)) s exp
x²
2
and thus
x²
P Etp 1/ Xe-1 2 x) s exp
Σ (1] Xi-11) x=-1)>>>11
+
2
n²
Σ || | 112 > A/ 112113
I cxp
+P
2022
What is left is to bound P
Et= x > p
Using Lemma H.2
we get,
T-2
I cexp
(54)
p ||=11²
2
16pTo? 1
1=1
Choosing x = = -p ()/ and A = | 233/2 yields
(-T pp)
(- c(T cT
ca2
P Ets n2 X 2
(T 6v8
S exp
Sexp
72p (T72p111)
+ cexp
||=112 T3/2
16po
Which implies, that for T > 50p, we have,
PETT 2 X so (T-))
EVT
1
}).
I cexp
min
(55)
p
75 1=
1602
for some absolute constants c, C' 0. Now, note that
min
min
75 || 12
16o2
75
||:
02
20
o²
min
75 0
|| 2 7
0²
05 1
2
1 1
1
02 2
+ ||
27

---

Thus, using 55 and 533 0, we have
3O- - p)
P(Eg)
P
P
M1,12
(56)
6vD
i=1
60 (T p)
2TP
Σ
2
T X 1-
(57)
6vP
i=1
2T2VT
S cpeepp
(58)
75po?
2 [|all)]
Using the condition T (e) concludes the proof.
Now we argue that the sum of the probabilities deduced from lemmas expressed in the simple form of 8(e, I') as defined above. First, from lemma
E.
E.2
and E.3
can be
E
we have
a² BvT
IP (EE) S cexp
3(11a12 12 1) o2
From lemma F.2 and T > To(c)= > 50p we get
-cVT p
-CVT
((E) S cexp
cexp
48
49
Lastly, lemma
E
says, for 2 To(e).
P P(E)) < ccxp
75po2 (2 + l1a11(2)
Thus we can write
cT²/T po?
P ( (UUUUU) < cexp
- 75 2 [|al))
= =((c
Where in the statement we set ē subadditivity of the spectral radii of symmetric matrices, we get with probability at least 1 - -((,,T) that
75 (2(+11a)) .. Using lemmas E..
and
E.3 and the
& EnEEE p -2 S ca²
dBBB el) 1 ET 1 (BB +el) - (v -
=l') x x 3 (T +el)..
Finally, we establish the following corollary. which is an immediate consequence of Theorem
Corollary E.1. Consider the AR(p) process x, described in (A.I) Let X (T)))xx such that
[X]s 214j-1- Then, for sufficiently large T such that log(T) probability 1 7m.
8 Amnin (aiii(()) ). we have that with
/T
22p
Amax Am)a(T))
1
Amin [
x) (x
> - 0
2 p)mmin (V),
(59)
2
Amax X
x) <
0²
2T >mmmmxx (
(60)
where e, C are constants and 7
d² o
vb (the controllability Gramian) and T are defined
75(2+
in (14)
28

---

Proof. Note that according to Theorem
E.1 Weyl's inequality (see Lemma H.3 and using c
11p (T)
TVT
have, with probability 1
TD
Amin (X x)
2 (T - p) (Amin (V) Edmax (T))
(61)
2 0T - p)) Amin (V) - 11 Plogg)
öVT Mmxx(T)
(62)
1 0 02((T p) (Amin(()))
(63)
2
Also,
Amax X
T
X) S 0² (T -p) (Amiul ( (V) + edmax (T))
(64)
3
Σ 2 0² ( pp) (Amin (V)).
(65)
What is left check the condition on T and c. Note that Theorem E.1 requires T 7 log(p)" m || 12 which with our choice of implies T p. To see that the choice of ofce
is valid, note that (T)
2/T < Amin( which validates the choice of e.
(I).
29

---

Theorem E.2. Consider the AR(p) process described in 477 where {mi) is an i.i.d. mean-zero sub-gaussian noise with variance o²_ Let a == [a1,.... 0p and a be its least squares estimator: If T N To(e), then for any unit vector w Sp.
det (VuVe" + 1p) 1/2
P ww (ā aa) 20 w T VE 1/2
log
28(e, T),
o(c,T)
where Ve, Vu To(e), and ò(e, T) are as described Theorem E.1
Proof. Write 11pt1...... 7/1 I (a(T) a:) = w [X X
Using
we have
T
X
T
E
1/2
(x
1/2
w
Now will need to introduce lemma from
11:
Lemma E.4. Let {Fr] be a filtration. Let {neee= be a real-valued stochastic process such that ne is Fe-measurable and is conditionally sub-gaussian with variance proxy o² > 0. Further, let (X,)e be an Rd- -valued stochastic process such that Xe is Fi-1- Consider any positive definite matrix V IRdxd For any 2 1, define
= Σ X,X, S Sc=End no
Then for any > O, with probability at least - for all 1. 2 1,
||alle S 20² dett /det(V)) (V:) det (V) 11/2
We will apply this result to bound the second term above. Choosing V Vv, Ee=1 = x and X E, we have with probability least 1 1-88 -(e,
det (XTX Vi) 1/2 det (Vi)-1/2
|X FllexT++V-- S $120
20 log
ō(c,T)
Now further condition on the event Theorem E. which occur together with the condition
above with probability at least 20(e, 7). Recall, then, that X X + Ve I 2X X, and so
(X x 2 xx X) This means that
|xTX)
F.
|| (X((XTx)-1
S 2
(x
det (XTX + Vo) 1/2 det (V2)-1/2
< 20
log
d(e,T)
det (XT XV, + 1,) 1/2
20
log
8(c,T)
det (V.VV- + Ip)' 1/2
T,)
< 2x log
ö(e, T')
Finally, it is clear that under the conditions in Theorem E. we have (x 5
V 1/2
2'
30

---

Corollary E.2. Lete= = 12 1l
where è is as defined in Theorem
E.1
Then, for T such that
EVT
(7) VT
Amin(() 24p
Amaa(l))
we have
| [ all2 > CM
log
(Tmax
c
Ddmin (Vo)
Xminn (V)
T
Proof. First, note that with the choice of c 12 (Tg
have,
EvT
8(e, T) = cexp( loooo (2)
(66)
T12
Hence, - log(((e,,
C Co)g(T)) Further,
w
V/
1/2
V
1/2
2
2
Amax (V² (V
1
Vmmin(s)
/ (Aminl (V) cdmaxe (T))
(67)
-p) (V))
where the last inequality holds with probability -(e(, T') using Theorem E. Finally, consider
log dot VV 1 + Ip) 1/2)
det (Vi + Vt) det (Ve)
log
log 2P 2 - p)d det (M) det (V))
(68)
Note that, again by our assumed condition e,
1
det (V)) 2 det (T)
2
0 2p - p) Amin (
(69)
2
Using 69 and (68). we have,
2 det (I)
log deet (VaV, + 11)
< log
Amin (V))
4Amax VT)
- p log
(70)
Amin (V))
From 70
67
and (66). we have,
det (VuV,- + 1,) 1/2 8(c,T)
20 w V 1/2
log
Cp
4Amax (V)
log
+ log(T).
(71)
(T p)
A)min
Amiu Amin(()) (v)
31

---

Using 71 and Theorem E.2 we get that with probability 1
712
2c
Cp
(4)max ( (V)) Amin ( (1)
la, (T,(T) asl $ V
(T P)Amiu (4) V
log
+ log(T))
(72)
For any E [p). This implies.
Cp²
44 (V) Xamin (V) + log(T))
a 0lI2 S with probability Tir which completes the proof.
(T( - ossl3 V (T Pmmmi ( V
log
(73)
32

---

Proof of Theorem 4.3
F.1
Setup
Recall 8* as defined in Proposition 2.2 and its estimate 3 defined as
NxTyL
B argmin
(ym 3 (mm-3"PP.) PT.)
(74)
To establish this bound, first define (and recall) the following notations,
Recall that Z, (yn, 1,) Z(y2, L,T,1) matrix of yı (t).... (t) V E [T].
Z(yn L, T, 1)] is the stacked Page
Let Zg be the stacked Page matrix of 2 21 Zn (t)() Vt [T].
Let Zy be the stacked Page matrix of fi ...t,..................... f (t) Vt e [T). That is, Z, Z, Zg-
Let z E [2(L11x be the sub-matrix obtained by dropping the L-th row from the
stacked Page matrix Zy Define Z', and Z analogously.
Let UEVT denote the SVD of z
V and U+ matrices of orthonormal basis vectors that span the null space of Z' and z' respectively.
Let ÜΣVT denote the top k singular components of the SVD of Z, while U Σ (V T denote the remaining L - k - components such that Z, = -.VV + (VL)T.
Let Z be the HSVT estimate of Z'y with parameter k. That is 2,- UΣVT.
*
F.2
Deterministie Bound
First, given the notation above, the solution for 74), can be written as,
B
Z2
Z
US
V IZyl.
(75)
Then, let's consider the following expansion for || -0-
| -
100-0-)10--))3--000000-0----
3²)
|
(76)
where the last equality follow from the fact that (Ū4)TB= = ((0))Tū) E V Zylle. 0. Now,
consider the first term in 76). and note that 8" TU
| 1 (Ü ) 3 2 U (U
TUU 3 +
u - 8/11
(Ú-(Ü"T--''0
8-
- 18
<
- 118"
(77)
33

---

Next, by Wedin sin θ Theorem (see 10 32
we have:
12 - Z f12
Ö uu 110 -
118
ox(lz))
11zl12
Tn (Z') 11
(78)
That is,
[ la
12 S
(Z,)) 11
(79)
What is left is bounding
|U (3 - 12. To that end, first consider
||
< 2 22, T 3 (2,) + 21(2)) - -12)) 1
s 2| (2) 3
111 + 21 22,11,000 1138 1|?
(80)
Also, consider
[122) (3 10 = (8 ) -- (8- 8r
0 20, (3 - >1
(81)
From
8 and 88 we get,
| - 18- S
2
(1(2,) - (Z,) 8 12 + 1 - 21,12,000 02,1,,00111"111 118
(82)
σk 12,)2
First, recall that [Zy] s.. = | + Zz] and that by definition, (Z²) = [zz)... Then the term (2)) 1 - (2,) ) 1/113 can be bounded as follows. First, consider
[1 - 1
[[(2)) - 14 1[ +2
[.,].
S [(2)) 8- - - -1- +2 (2,) - (2z) 8*) 1Zol..
3-
((2)) 8- - -14 11 - || 12 +2 (22) - 3
[2m]t.
(2) 2;) T 3" -IZzL
-11Z.11-1 +2 2 (2,)
(x)) T [Za|z
[Zul..
2
[ zi) T
8r + (2)) (2,) 8t)
(83)
Note that T3 VER 3 VEŪTÜΣ'V [zyle VVTIZ and (3)) VEÜT B*, thus,
(Z 3 - (2,) 8)) T
[Zo]z. InyII. VV - 3d UEV ZolL
Zunl
[Z.lI. V 8 -U [m,
(14 V [Z. V 8r 3'T -U)
T [Zm] I..
(o 4 A 8" 2, V
VT Z 1z. +
|2 V
< 0 Z - 8 TZ
V11 T Zz.l
|
+ I V
< 11 - 2/ 12,00011111111 /1
Zz 12 V
+
Z I.
(84)
34

---

Finally, from
82 and 83
we get
2
1
Ü
11)1
21|21,
3 00 1
ox(27))
+ 1z
1112,,001|/*111
[ V
4
21121,
2,112,0 | 111
[Z
(Z,))
2 can be further decomposed as
The term [z 1. V
Zz
(VV
VVT || +2 [z VV 1|
Using Wedin sin θ Theorem,
VV VV S 12
(Z,)
Therefore,
1
11
<
8
||Z/, 2/ 12 11
IZ 11 ok m
|| 12++
[Z. V
oN (2
(85)
From
76
79
and
85
we have
I 12
|18-112 + || 1. 112
|
1
T (Z2)
σx (2,))2
2 2| 11, 18 + [|
(86)
T TR (2,))
High probability bound
Let = /T/L be the number of columns the stacked page matrix Z/,. and be some positive absolute constant. Define
E
|| 1l2 Convaj.
Ey [I[Ze|L.1l S Coovaj.
E = {1Z) -2,11,00 5 Co
(0 sZax) + Cozklog(())
(Z,)2
E = Z..... 1/, I Co log(g)
Lemma F.1. For some positive constant ej > 0 and and Ey
> 0 large enough in definitions of E2, Eyg E¹3
PPE)) 2 2c 019
P(E)) > 1 2e -019
IPE)) 2 C1
(IT)
I(E) > (NT)
35

---

Proof. We bound the probability of events above below.
Bounding E and E. Note that Z is sub-matrix of Zg. and its operator norm is bounded by that of Zg. Hence, as the probability is an immediate consequence of Lemma A.2 The second event probability follow from the fact that [|Zz]e-Ilz I Zz
Bounding E'G. This is immediate from Corollary C.1 and its consequence in (2
Bounding E First, recall that [Z L. is q- q-dimensional sub-gaussian vector with variance proxy of.
That is, for any unit vector v Sy-1, [ subG(o Now, recall that V E Tge* is set of orthonormal vectors V,, [k]. Let Z, [Z,]t. V,, Clearly, Z; ~ subG(o Therefore, [a,r.. V is a k-dimensional vector with dependent sub-gaussian entries with variance proxy oz Using Lemma H.4
we have
PP (I[Zz]) V122 cexp
16ko
Therefore, for choice of = Co, v k (o)((g) with large enough constant C 2 19, we have,
P(I1 P(llZ. V Co. Co.V/llo(()) < 9²2
Recalling that q > L, and q = NT L, concludes the proof.
Let E' = Ey Egg
Fgg 0 E.
Then, under event F., and using (6). we have, with probability
[TT)T
Co
||3
3
11188111
σk (Z2)
ox
Co
g²
(0z + m + k kog(o)) max max(1|18111 (118111
Tz (22)2
(Z,)
Choosing L /NT, and k RG, recalling Assumption 4.1 and using Weyl's lemma (H.3 to bound 1ox (2,)) ok (Z' )2 we get,
Cg NT
- 3011
VNT 10 RG log (NT)
( RG (0 faan)) + (Glog(NTT)
||9*||7
Cy S2ne o RRR lNT max |11." 1),
VNT
which concludes the proof of Theorem 4.3
36

---

G Proof of Theorem 4.4
In this section, we prove an upper bound on the forecasting error,
N 27
ForErr(N,
(Tn(1) - Bn(())²
NT
I=T+1
where (t)
= 18 Im (t) (1), ,.0)...,,..... Un (t
1)]
fn(t) + In((), where Zn(()
In (t) Fn (1). It 1)]. First, recall that for T
(t) = (t) + (t) 3 (t 1) + a xn (t 1),
where (t - 1) = (tn(t (I. - 1)), . (t 1)) and Xt 1) [T p),..., 1)]. Then, we have,
N 27
N 27
2
Σ
Σ
(u()))
Σ
Σ ++1
(fa(t) - f.)(1)) + (3.(() Fn (t))
NT
NT
-=11=
Now, we will bound each error term separately.
G.1 Bounding Forecasting Error of
To bound the forecasting error of f, we use the following Lemma.
Lemma G.1. Let the conditions of Theorem 4.1 hold. Then, with probability 1 (NT)
N
Σ (m(0)~((1)) < C(fmax, G max (1 I7 log (NT)
2T
Σ
NT
VNT
n=li=T+1
where is absolute constant, and C(Jmax: 7) denotes constant that depends only (polynomially) on model parameters fmax," "..
Proof. The proof of this lemma is similar that in (3) but adapt it for different settings. First, note that, where we assume that NT/LL is an integer,
NT Σ Σ (I.(0)
2
N T/L-1
f mL Sa +mL+ mL++0)
1=T+1
In what follows, we provide upper bound for the inner sum NT EN ET/Lo (I +mL
mL
0) - (T +m+++++)))) for all é E [Z]. Next, we show, without loss of generality, the case for l = 1.
To establish this bound, we first define (and recall) the following notations,
Recall that Z, = [Z(y, L,T) Z(y2, L,T) Z(ya, L,] the stacked Page matrix of (t)....... (t) for t [T].
Similarly, let Zyoout denote the stacked Page matrix of the out-of-sample observations y1 (t), ...)................... (t) for t E -T- L 2, .,.... 2T L 1}, i... Zy,out E I...NT/L
Let Zz be the stacked Page matrix of z (t), ..... xn (t) for E [T', and ou be the stacked Page matrix of 01 (t)()..,,,()) for E {T + 2, ..... 2T
1}.
Similarly, let Zg and oout be defined similarly but for (I).
Let Z, E ()(T)))(NT/LL) be the sub-matrix obtained by dropping the L-th row from the stacked Page matrix Zy Define Z Z'. f.,outt L'z, and Z'z out analogously.
37

---

Let UEVT and Oout Dout Vout" denote the SVD of Z' and 2 fout respectively.
let V and U be matrices of orthonormal basis vectors that span the null space of Z and Z'; respectively. Define Vous" and U similarly for the matrix Z f,ot
Let UEVT denote the top k singular components of the SVD of Z while U- (V-)7 denote
the remaining I, - k components such that Z, = UEVT -- VU'E(VE (V-)T.
Similarly, Let Unus Dont V out denote the top k singular components of the SVD of Z yoou while Ü. D out Four (ou) denote the remaining k components such that Z U Outt Dow VT out + U Dtt
y,out
out out (Vdh))
Let 2, UEV and f Uout Dont V out
Recall that we are interested in high-probability bound for the following out-of-sample prediction errоr:
L
T/L-1
+L fn(T mL + +1)
Σ
(87)
NT
m=
Recall that Fi(() TT,
Yr (t 1) and hence we can rewrite 88 as
T/L ΣV
2 = Zyyoout 3- (2 [2.ooul
fn m
(88)
Z noutt T zoo T 3000 1/
(89)
Next, we derive deterministic upper found for the expression above.
Deterministic Bound. Through adding and subtracting 2, Joue and triangle inequality, we have
(Zyoouc) T3 (Zi) 3+ 1=11- < ||(Zyonni 2B,oont) 13/ 1 1.. (o)) 3r (90)
Next, we proceed to bound each of the two terms on the right hand side.
First term: || iz
2 fout))
(Z) (Zyonu
2 T
1 - out! (Ū ont)) /1
<
Σ out 1 (
out ) 313
(91)
Note that 1 out ll2 equals the (A: + 1)-th singular value of Z "oo Using Weyl's inequality (see Lemma H.3 and recalling that Zy = L'iioou + Zz and that the rank Z' y.uu is k, we can bound the (1 + 1) singular value Z oo by the operator norm of Z' Zront That is,
Σ Doul 112 < Iz out
(92)
Next, we bound the term || ((+)))) 311.
38

---

|I )ut) TB12 = || out Ü out 13 12
U out Ü out) T 3" U out (U T (3- B
- 2 U out U out) T 3 + 2113 3"1
(93)
Further expanding the first term,
Ü out
Ü ont)
8
U out (U (oue) 8
U out (U (Wott*) T) Uout (Uout) T
U (TnelTuntt (Tont)
U Uout" ' (Uout ()T)
Uout (Uout)
T (Uoht) Uout (U )7)
13
U out (Udut)
- U out - (U "Uout" 4)T
118
Ü ÜT out
Uout Uout T ||
Where in the first equality we use the fact that 3 Uout(U (out) T 8r, i.e., 3 lives the column
space of 7 yout (Assumption 4.2)
Ü Tout Tout U Uout Uout" T
Next, by Wedin sin θ Theorem (see
10
32 we bound
as follows:
|| 1Zzzouellz ||35
U ÜT
U Oout
U out ||
|| -
(94)
out
0x (Z
Using this definition, 91). and 94) we have
1 1133113 Th (oue))
oot 2 ( TB| S 2 2 out 1
+ 113 3
(95)
Second term:
2 j.or (Z/, )T To bound the second term, we follow similar proof to
that shown in 5).
7 ooue 3 (Zj
1/0 = foout 3 2 jout 8r Ž J.our -(Z' TT T/31
2 7 oout (8 301 + ( f.out - (Z, g/
< 2 2 loout [.o0( - 188 )13 + 2| Z 2// Joou 22,,oull 112,0011330 1|?
(96)
Next, we bound the term
2 fout (ā- 11.
2 [.onn (3 18²))1 < 11 ( jout I's,ont + Z (.ou)) TT 1B- 33
(97)
< 2/ 2 f.out (ou) T( 3²' + 21|2Z jouu T [8 - 8*)
(98)
< 2 1', oooout 112113 3r + 2 2Z' f.uut [@ 1-1")1
(99)
Next, we bound IZ' oott [ 11 Recall that U spans the column space of Z your. Thus
Z'youut = yoot U therefore,
/uu (8 - -31 < Zz f. out 1 (Ø - -3)
IZyoout T TU (8- - 8²)
(100)
39

---

Now recall that we denote the top l left singular vectors of z by U, then consider
1000 -0 1000
22||UU* 1813 10 33-
(101)
Recall that, as we prove in
85 in Appendix
|
| (8 - 3-11
IZ' 1 na 18 ox (Zz)) )2 1L + || v11
ox 12 )
(102)
Using 101p. 1102 and Wedin sin θ Theorem, we obtain,
12,12
||3 1112+ 8 81 1-11
| (8-3") - 2
Th (2,))
Th (2())
8
IZ' 1 118 11 vI2
(103)
Tx (2,))2
Therefore, using 96
100). and 103 we have
1 f, out
(Z) (Z T 3/2
$ 2 2 T
four (3 - )1 1 + 22|Z' foot 2 fou 11. |
< 4 Z oort 1121 11
1 ||
|13- 3- + 811[2,]r.112
4| 1 our 113
(Z))
(x(2))
2|Z 2111 118 11 v11
| yout
T% (Z,))
(104)
Combining. Incorporating the two bounds in
95 and 104 yields,
Z out 2 $133112
(Z) (oo T -3--zz,oo (z))
S e
+ 113 -3113
Tk Z fome))
| 12
12 11 || [Z,, L. 12
+ cl I'' pulllll
1311
0x (Z2))
ox )4
+ IZ four l13
2/|Z/, 7, 11,350115 || + 7 |, v
0 (2)))
+ c (.out
(105)
For some absolute constant c > 0.
High probability bound. With our choice I. VNT. Let == VNT be the number of columns
in the stacked page matrix Zy and subsequently a + to be the number of columns the stacked
page matrix out Let C' > 0 be some positive absolute constant, and C(fmax: 7) be constant that
40

---

depends only on model parameters fmax and 7. Define
E = | 12 S Coavaj.
(106)
E2 : 1 112 S Convaj.
(107)
G2 max 119)11
Eg
| - 8r) I =((mmax:7) 7)
(108)
9
Co392
E1
1 11100000
ox(21)² (oz + fmax) + Co log(a)}
(109)
Cozy²
Es 12' %,ouellB... S
(oz + fmax) + Cook
(110)
ox (2,))
E6 = [Z Vl22 S Co2 v /loog())
(111)
E, = [|z,]1-122 - Jo.vaj.
(112)
First, with e > 0 is an absolute constant, recall from Lemma P(E)) 2
F.1
2e
-19,
E(E)) 2
C1
(NT))
P(E6) 2
Ci
(NTT)
TP [7) 2 2e 000
Further, note that out is sub-matrix of out, hence its operator norm is bounded by that of Zu,out. Therefore, as an immediate consequence of Lemma A.2
P(F2) 2 019
The probability of event Eg is bounded by Theorem 4.3 as
PP(E)) 2 (NT)
C1
Finally, as an immediate results of Corollary and its consequence in 27 we have
TP(E))
C1
(NT)
Now consider the event E = En Eg Eg n E. Es n Eg n Er. Then, under event E, and using 105 we have, with probability (NTH
ZZ (oo TB -(Zy,oou) 8- s C(fmax 7)0
11911 The (2ooo)))
+ (
max 1} SGG2R2lo((0)
o2
+ C(fmaxi Dlll''',.n
max{ 12, 11
024
+ c| 2
S.ouell
0x (2,)4
Co,92
+
1188*13
(Z,) (0 + Plnx)) + Co2k log(g)
Tk (2,))
Z /.out 112
+
2 log(4)
Nk (Z2)2
02.92
+ (2,)
(02 + f ||6* 11
(113)
00 3 k log [3) R.
(114)
41

---

Recalling VNT, k = RG, and Assumption 4.1 we get,
( outt) 3 (ot 3 2 C(fmax: ( 7) max 13 log(NTT.
(115)
Therefore, with probability 1 - 1-c/(NT)"
N TVL-1
Σ Σ
+ + 1) [T mm++1)
(116)
NT
n=l m=0
max{11 (og((TT)
I C(mmax 1)
(117)
VNT
Further, note that the same argument can be used bound the sum NT Σn= E (nu(T
f. for all & E [Z]. Thus, with p 1-cL/ > 1- c c/(T))
we have,
2T Σ
(fa(t) 100--f1h (a(0)) s 7 Go max 112 13 log(NT))
(118)
NT T=
VNT
G.2 Bounding Forecasting Error of
Now, bound the error of forecasting 2, (t). First, recall (t) yn((t (t(,, and let (t) (t)(t),.... In (t p + 1)] and (t) =
2, (t -p--pp--- p 1)]. Now, let Xn, X, and A,, all
in ((T=))x,, be the row-wise concatenations of (1), X. (1), and (1) = =xx(() (() for E [T p,T L..... 2T + p p- 1}.
Then, we have, where ān (t)= = E | (t) an (I), 1)]
2T
T
(in(1) =1)) =
T x
(119)
7=
X 1
i X + +ā
X
Ox 12
(120)
2
2
2
à X a x n
a x
T.
a
T
(121)
2
2
| | x n.
a x
+ || a x
(122)
A
||اm 1 + 1 aml1)) + +|
16
-an 1
(123)
First note that
Σ
22 < Σ
A
2
(124)
N
27-1
<p p Σ
Σ
(In(t) - (1))
(125)
1=11=T-p41
#=1
N
2T-1
=
Σ
Σ
( ?
(126)
1=1t=7-p41
42

---

Therefore,
27-1
Σ n= 1=++1
[n(1) (En())=a(0))
3 4px (1as + l anll))
E Σ ==-=
(In(1)
(n(())
Pa(1))
(127)
N an.1
x a..
(128)
High probability bound. For the high probability bound, consider the event E Es n Eg Es where
N
GRp of (Ng( V
Cp²
T L'maa (1L) Amin 41
E =
a - an112 C(Jmax 7)
max
1, 0² xm2mm (W)
log
Amin
0²
VNT
T TXmin (V)
n=l
(129)
2T-1
Eg
Σ
Σ
(ratt)-Fa(t)? s [((fmss.))R8&'0 (maxf]ls*1[]
13 (oTT))
1-p+1
(130)
N B. n
{
|x
00 p) Xmmaxx (V)).
(131)
2
n=1
itT using Theorems 4.2 4.1 and the the union bound. Further, note that (NT) To using both Lemma G.1 and Theorem 4.1 Further, note that according to Corollary E.1 with probability 1 Tír, we have
Note that (P(E)) 2 1
P(E2)
> 1
C
3 0T p)>max( X , 112 Amaax X x S
(132)
Therefore, under E, and recalling N
T, have with probability 1 - 718
C
2T+L
Σ Σ T+h
(n(() sn((()
(133)
NT
S C(fmax: G) R Amax (V) o max [1811711 Amain (v)
2po² log
T (Txmax (V) Amin (1)
(134)
GRa min {1, o2 Amin (V)
o (NT(
(135)
VNT
where O'mx = max, llanll2 therefore, with probability
TIo. we have
C
N
27
(1)()) vn(()) 2
Σ
(136)
NT
1=T+1
Amaax ( (V)
ppo2 T
T T>xax (V)
C(fmax: 3) Ra
Amin 4 0f max 11.
log
(137)
Amin v
GR GRIlaali min o2 Amin (V)) 02 VNT
σ log(N))
(138)
which completes the proof.
43

---

H Helper Lemmas and Definitions
A zero-mean random variable X E IR is said to be sub-gaussian with variance proxy o² (demoted as X ~ subG(o²)) if its moment generating function satisfies
Definition H.1 (Sub-gaussian random variable 30 31)..
02
IB EleexpsxX) v exp
ss R.
2
A zero-mean random variable X E IR said to be sub-exponential with parameter (denoted as X subE(v)) if its moment generating function satisfies
Definition H.2 (Sub-exponential random variable 30
31
12
1|s S 1 V
E[exp(xXX] S exp
Lemma H.1 (31 30 6). Let ~ subG(o??) and Z NO., ozz), For i E [N), let Y; ~ subE(vi).
Then:
1. sub)((o)))
2. 7 ~ subF(ag)
3. X2 -[[X] subE(160
4. Em Yi subE v4)
Y are independent:
5. Σ Yi subF( 22)
Definition H.3 (Sub-gause gaussgansiinn vector). random vector X E TRid is sub sub-gaussian random vector with parameter o² if
of x subG(o?), Ve e gd-1
where Sd-1 = fir RR; ||]|= 1].
Lemma H.2 (Norm of sub gaussian vector). Ler E Rd be sub gaussian random vector with parameter o². Then, with probability at least 1 - 8 for d E (0.
| 40vo+ 20 2 log (1)
This implies that for any t 20,
| S t
with probability at least exp(- 16d
26 ). process x(l) E. is said be stationary if its expectation does not depend on t and its autocovariance function depend only on the time difference. That is,
Definition H.4 (Stationary process
([(x(t)-H) [E[+(t)] 4
Vt
11) = 7())
Vt and any j.
Lemma H.3 (Weyl's inequality). Given A, B E Rm let σ; and G; be the i-th singular values of A and B, respectively, in decreasing order and repeated by multiplicities. Then for all E [m A ni),
1o, न. < A|A Bl22-
Lemma H.4. Let X [Xi...., Xa) where each X,, [n]] is sub- sub-gaussian (not necessarily independent) random variable with variance proxy o² and E, x?] 7. Then with probubility at least 1 cexp "Iooa"
|X112 <t
(139)
44

---

Proof. From Lemma H.l we can see that | = Ei=1 x; - -ny is sub-exponential with
parameter 16no². Let d Ei= i= x
ny, then E[exp(sd]] exp(st') I exp
u23² st'
1ss| S -
P(d> t')
(140)
v
Choosing =2 n7, v 16n0 and 16n
we get,
P P(d> 12 - -m)) = P (1x11 12) $ exp
exp
-
mn
(141)
16nd2
exp
exp
exp
(142)
160
16n
t
cexp
(143)
16no
where the last inequality follows from the fact the second moment of x is bounded by its variance proxy up to a multiplicative constant.
Lemma H.5. Ler Xi...., Xn be ii.d sub-exponential random variables with parameter v. Ler
S Σi= X,. Then
W/
P ([S| $ nt) cexp
(144)
Proof. Recall that the sum of n i.i.d. sub-exponential random variables is also x sub-exponential random variable with parameter v (see Lemma H.I...H.II.... Thus, using the definition of sub-exponential random variable, we have
[lexp(AS)]
/nw222
1 Vnv"
(S > nt) S
I exp
Ant
1ss S
(145)
exp(Anf)
2
Choosing À- Vnoi get,
1 exp
nu
PSS > nt) S exp
(146)
VTV
-lyn
< 2 exp
(147)
Applying the same inequality to S completes the proof.
Lemma H.6.
133 Let E mx and B = A E. Then
B All max {|
E
where M (M M) M) '' "M denote the Moore-Penross inverse of M.
Lemma H.7. [Lemma 4.2 of [271] Let (Filizo be filtration, and (Ze Ze and (W, W.hh2 be real- valued processes udapted to F. and Fetl respectively. Moreover, assume Wr|Fr mean zero and o2-sub-gaun gaussian. Then, for any positive real numbers or, 3 we have
PP (
8))
a²
Σ
Z.W, 2 Lann
Σ 1=1
2 <
< eexp
(148)
2003
45


Reviews:
[{"review_id": "0CSLwbKx0S", "reviewer": "Reviewer_rVQz", "summary": "The paper discusses a two-stage algorithm for time series analysis, which involves estimating deterministic, non-stationary trend and seasonality components, followed by learning the residual stochastic, stationary components. The first stage involves using multivariate Singular Spectrum Analysis (mSSA) to estimate the non-stationary components, even in the presence of a correlated stationary Autoregressive (AR) component. The AR component is then learned from the residual time series in the second stage. The authors provide a finite-sample forecasting consistency bound for SAMoSSA, which is data-driven and requires minimal parameter tuning. The paper also presents empirical studies that validate the superior performance of SAMoSSA compared to existing baselines. Notably, SAMoSSA's ability to account for AR noise structure yields improvements ranging from 5% to 37% across various benchmark datasets. The authors also provide a detailed explanation of the model and assumptions used in their analysis, including the spatial and temporal structure of the time series and the properties of the AR processes. They then describe the algorithm for the univariate case and explain how it decomposes the observations into estimates of the non-stationary and stationary components.", "questions": "The authors should provide simpler explanations or visual aids alongside the more complex mathematical definitions and proofs to make the paper more accessible. More extensive empirical evaluation should be carried out with a broader array of datasets to prove the efficacy of the method. The impact of the violation of some assumptions should be discussed, indicating potential limitations of the proposed method. The authors should provide more discussion on how this method could be generalized or modified to tackle different types of time series.", "limitations": "The authors discuss but do not highlight the potential limitations of this work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.\tThis paper assesses the importance and novelty of the proposed SAMoSSA methodology and discusses the paper's originality in the use of mSSA and AR processes to jointly learn the deterministic non-stationary and stationary stochastic components of time series.\n2.\tThe paper addresses a gap in the literature, namely the theoretical underpinning of multi-stage learning algorithms involving deterministic and stationary components.\n3.\tThe results demonstrate significant performance improvements over existing baseline approaches.\n", "weaknesses": "1.\tThe complexity of the paper could be a potential barrier, making it difficult for people without theoretical interests in the area of time series analysis to understand. The paper could provide clearer intuitions and simplifications alongside the complex definitions and theorems to make it more approachable.\n2.\tThe empirical evaluation could be slightly limited. The variety of benchmark baselines used to demonstrate the efficacy of the method could be extended. The analysis of the results should be discussed more carefully, providing insight into why SAMoSSA outperforms other models. It is unclear whether the assumptions made in the paper hold in some of the typical practical scenarios like the electricity and health datasets.\n3.\tThe paper should discuss in more depth where these assumptions come from and under what conditions they might not hold. It should also explore how the method could perform if some of the assumptions are violated. The method proposed seems tailored for a specific kind of problem and it is unclear how generalizable it is. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We express our gratitude to the reviewer for their positive and insightful feedback. In what follows, we address the specific questions and comments raised.\n\n> The authors should provide simpler explanations or visual aids alongside the more complex mathematical definitions and proofs to make the paper more accessible. \n\nThank you for the feedback. While we have included some discussion and aimed to provide intuition for the various theorems/definitions, we will revise them to make them more approachable. Can the reviewer kindly help us identify particular definitions/theorems that were particularly complex and inaccessible?  Thank you again!. \n\n> More extensive empirical evaluation should be carried out. \n\nAs the reviewer correctly points out,  the theoretical characterization of multi-stage learning algorithms in this setup is missing in the literature, and this is the main focus and deliverable of our paper.  That said, based on the collective feedback from reviewers, we added a new baseline to our experiments (DeepAR, see attached pdf), and we will attempt to add TRMF [1] in our revised manuscript. Interestingly, DeepAR only outperforms SAMoSSA in the traffic dataset. \n\t\n\n> The impact of the violation of some assumptions should be discussed, indicating potential limitations of the proposed method. \n\nWe thank the reviewer for raising this important point. To make our contributions more practical, we will add a discussion of a diagnostic test (developed in [2]) that can verify whether the main assumptions we make are valid. This will hopefully indicate the limits of the proposed method and allow practitioners to know when to use or not use the method. The test verifies whether mSSA is likely to succeed, given the observed data.  In particular, the test measures the effective rank of the Page matrix associated with the multivariate time series with parameter $L$ ~ $\\sqrt{NT}$. If the effective rank does not scale much slower than L then mSSA is unlikely to be effective.. We will add a reference and discussion of this test in the revised version. \n\n\n\n> The authors should provide more discussion on how this method could be generalized or modified to tackle different types of time series.\n\nIn the limitation section, we identify two areas of future work that can generalize the current method. The first is extending the model to include non-stationary stochastic processes (e.g., integrated processes). The second is considering a VAR-like model where there exists some dependence between the stationary processes $x_1$,..., $x_N$. \n\n\nWe again thank the reviewer for their comments, constructive feedback and suggestions. We hope the reviewer will take these clarifications into account in their revised scores.\n\n**References**\n\n[1] Yu, Hsiang-Fu, Nikhil Rao, and Inderjit S. Dhillon. \"Temporal regularized matrix factorization for high-dimensional time series prediction.\" Advances in neural information processing systems 29 (2016).\n\n[2] Agarwal, Anish, Abdullah Alomar, and Devavrat Shah. \"On multivariate singular spectrum analysis and its variants.\" ACM SIGMETRICS Performance Evaluation Review 50.1 (2022): 79-80.\n"}, {"review_id": "Zi2h1psP1o", "reviewer": "Reviewer_Xzdt", "summary": "This paper extends previous work on multivariate Singular Spectrum Analysis (mSSA) to observations with autoregressive (AR) noise. The method constructs a sliding window representation of the target univariate or multivariate time series called the Page matrix and learns the deterministic non-stationary component (with a linear model) from the singular value decomposition of this representation. The model then produces forecasts by first estimating the parameters of a linear model for the deterministic component followed by the AR parameter estimates of the AR model in the second step. The authors provide theoretical bounds for the error estimates under AR noise of this method, the AR model parameters under perturbation, and the out-of-sample forecasting error. Finally, the authors compare the performance on their method with other time series models (including the mSSA method that is the basis for the presented algorithm. ", "questions": "Why is there such a large improvement from SAMoSSA over mSSA for Electricity? I would appreciate any insight on why the proposed method performs much better here. ", "limitations": "Limitations are briefly discussed (and the discussion is sufficient in my opinion). ", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Originality\n\nFrom an algorithmic point of view, the proposed method is a straightforward extension of the previously developed mSSA model to incorporate autoregressive noise and use a two-step model to learn it’s parameters (and the parameters of the deterministic part). Of course, I think the main contribution lies in the theoretical error bound results but I am not familiar enough with the theory to judge the originality on these and will leave this to other reviewers.\n\n\nQuality\n\nAgain, I am not an expert on the underlying theory but so I cannot judge on the correctness of the theoretical results. The presented method is an interesting approach to time series forecasting and the theoretical guarantees are appealing. The empirical evaluation on the synthetic results confirm the established error bounds on that synthetic example and the paper also contains a (albeit brief) empirical evaluation with different baselines. \n\nClarity\n\nThe paper and its contribution are clearly written. However, the paper is very much written from the perspective Singular Spectrum Analysis point-of-view. I would have appreciated more background and discussion of related approaches in classical time series literature (which I will go into detail in Weaknesses). \n\nSignificance\n\nTime series methods that perform well in practice and also have provable error bounds are not very common in the recent ML literature that has focused mostly on deep learning models that lack these guarantees. As such, the method presented in this paper is significant, particularly for high-stakes applications where provable error bounds are required. \n", "weaknesses": "The paper is brief on related work, especially on related approaches in classical time series forecasting (STL decomposition, simple exponential smoothing (SES) with drift). I think it would be useful to the reader to contrast the presented method to these approaches. Another class of related models are matrix factorization models that are not discussed either. Again, I think it would be useful to the reader to discuss this class of methods. \n\nI understand that the main focus on the paper is on the derivation of the theoretical results and I appreciate the effort here. However, the quantitative evaluation is very brief and does not give much insight into the performance of the method relative to either closely related models (other than mSSA) or modern deep learning methods. \n\nI would appreciate if the authors would at least compare their method that is conceptually or algorithmically similar. To this end, I would propose to compare against SES with drift and at least one matrix factorization method for forecasting (TRMF for example: https://dl.acm.org/doi/abs/10.5555/3157096.3157191). A comparison with deep learning baselines such as DeepAR or Fedformer (https://arxiv.org/abs/2201.12740) would also be interesting to understand how the method compares to currently used deep learning models. Simple baselines are also missing (seasonal naive for Traffic/Electricity; the naive method proposed in Bergmeir et al., KDD 2022 (https://link.springer.com/article/10.1007/s10618-022-00894-5) for Exchange).\n\nI would also appreciate more insight on the non-stationary deterministic part learned by SAMoSSA. I’m actually wondering on the added benefit of that part on the Electricity and Traffic datasets, which are quite stationary (but this might be up to debate).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their detailed and constructive feedback. In the following sections, we address each of the questions and comments they have raised.\n\n> However, the quantitative evaluation is very brief and does not give much insight into the performance of the method relative to either closely related models (other than mSSA) or modern deep learning methods.\n\nAs the reviewer correctly points out, the main focus of this paper is to address the missing theoretical characterization of commonly used multi-stage learning algorithms.  That said, based on the collective feedback from reviewers, we added a new baseline to our experiments (DeepAR, see attached pdf), and we will attempt to add TRMF in our revised manuscript.  Interestingly, DeepAR only outperforms SAMoSSA in the traffic dataset. \n\n> Why is there such a large improvement from SAMoSSA over mSSA for Electricity? I would appreciate any insight on why the proposed method performs much better here. I would also appreciate more insight on the non-stationary deterministic part learned by SAMoSSA. I’m actually wondering on the added benefit of that part on the Electricity and Traffic datasets, which are quite stationary (but this might be up to debate).\n\nWe thank the reviewer for this great question, as it gives us the chance to clarify our results and highlight the advantage of our proposed method. To explain the improvement, let us first describe how the two forecast estimates (mSSA vs. SAMoSSA) are different. In mSSA, assuming the univariate case, the forecast estimate is\n$$\\hat{y}(t+1) = \\hat{f}(t+1)  =  \\sum_i \\hat{\\beta}_i y(t+1-i)  $$\n\nWhereas in SAMoSSA, the forecast estimate is (using the same $\\hat{\\beta}$ above)\n$$ \\hat{y}(t+1) = \\hat{f}(t+1) +    \\hat{x}(t+1)  =  \\sum_i \\hat{\\beta}_i y(t+1-i)   + \\sum_i \\alpha_i  \\tilde{x}_1(t+1-i) $$  \n\nThat is, mSSA as described in [3] overlooks any potential structure in the stochastic process $x(\\cdot)$ as it assumes i.i.d mean-zero noise process, while in SAMoSSA the structure of $x(\\cdot)$ is captured through the learned AR process. Given that, the difference in performance would be attributed to the structure in the (estimated) stochastic process $\\hat{x}_(\\cdot)$. That is, if there is an AR structure in $ \\hat{x}(\\cdot) = y(\\cdot) - \\hat{f}(\\cdot)$, then we expect SAMoSSA to perform better. \n\n\nInterestingly, we indeed find this to be the case -- in the electricity dataset, the partial autocorrelation coefficient of  $\\hat{x}(\\cdot)$ at lag $1$ is significant (we focus on lag $=1$ for brevity and since we find it the highest coefficient on average). In particular, Figure 1 in the attached PDF, shows that on average (across the different univariate time series in the electricity dataset), the partial autocorrelation coefficient at lag $1$ equals $0.2$. \n\nWe see a similar but weaker partial autocorrelation in the traffic datasets, with an average partial autocorrelation coefficient (at lag 1) of $0.1$. This could also explain why the improvement in the traffic dataset is relatively smaller.\n\n\n> The paper is brief on related work, especially on related approaches in classical time series forecasting (STL decomposition, simple exponential smoothing (SES) with drift). I think it would be useful to the reader to contrast the presented method to these approaches. Another class of related models are matrix factorization models that are not discussed either. Again, I think it would be useful to the reader to discuss this class of methods.\n\nWe thank the reviewer for their great suggestions. We agree that the related work should discuss both topics and we will add them in the revised version. \n\n\nWe again thank the reviewer for their comments and constructive feedback. We hope the reviewer will take these clarifications into account in their revised scores.\n"}, {"review_id": "x3Jw2zrXud", "reviewer": "Reviewer_s3mh", "summary": "This is a comprehensive work on a new extended variant of multivariate Singular Spectrum Analysis (mSSA), which manages to handle time series of deterministic trend/seasonality with AR stationary components, with rigorous theoretical guarantees. The algorithm is a natural extension of the variant of mSSA using the Page matrix. It further provides guarantees on the bound of estimation error of non-stationary components, the bound of AR parameter estimation error, and the bound of out-of-sample forecasting error. And its superior performance is shown in experiments on four datasets.\n", "questions": "Besides the questions in the \"Weaknesses\" part, there could be more deserving the authors' attentions.\n\n1. Is it correct that we consider such a simplified multivariate time series modelling problem that each channel/variable (i.e. y_n) is driven an independent AR process (x_n)? Together with Assumption 2.1, we actually can model each channel signal y_n separately. This also allows the stack of Page matrix for an extension from the univariate to the multivariate.\n\n2. Following the previous question, the analysis for SAMOSSA for multivariate case has nothing essentially different from the univariate case, right? Maybe we missed some technical challenges.\n\n3. In the proof of out-of-sample forecasting error bound, what is the particular issue that our problem has compared to the analysis in [3]?\n\n4. For time-series forecasting, there has been plenty of SOTA methods using deep learning (like Informer, PatchTST, etc.). We are curious about the comparative results. Indeed yours may be less accurate than DL SOTA, while yours has additional explainability. The performance gap might tell the cost for explainable modelling.\n", "limitations": "We do not see any potential negative societal impact of this paper.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The extension of the algorithm naturally follows by adding the second step of AR estimation. It is impressive that this paper manages to show all these three bounds on estimation/forecasting errors. What is surprising that, for such a two-stage approach, the bound of out-of-sample forecasting error could be attained.\n", "weaknesses": "The numerical experiments could be richer to show the performance boundary of the proposed method, considering it is dealing with a general signal as the non-stationary plus the AR components. The types of non-stationary components may affects the performance; in particular, instead of the complicated one, the simple deterministic components that could be also be captured by (close to marginal stable) AR parameters may raise issues. And the signal-noise-ratio, the intensity of AR process compared to the deterministic, etc. may also be tested in experiments. All these simulation may help us to further understand the strengths and weaknesses of the proposed method. The reason why the reviewer expects more numerical testing comes from the concerns on the power of truncated SVD of Page matrix for approximately splitting the deterministic/non-stationary and the AR components. If this step fails completely, the whole would fail. We like to know its performance boundary and reliability. \n\nAnother concern may be taken into account by area chair. This paper is more typical in the fields of econometrics or statistics. The reviewer is not sure if it fits NeurIPS well.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive feedback and valuable questions. Below we address the specific questions and comments they raise.\n\n> The numerical experiments could be richer, if this step fails completely, the whole would fail. We like to know its performance boundary and reliability.\n\nWe agree that it is important to characterize the performance of the algorithm under different settings. However, we maintain that our analytical approach provides a deeper understanding of the algorithm's behavior compared to adding more numerical experiments, which can only cover a limited number of scenarios.\n\nFor example, although it can be conveyed more explicitly, the analysis we provide helps characterize the performance in terms of the signal-to-noise ratio. In particular, we capture the noise effect through the dependence on the term $\\sigma_x$ in the various bounds; on the other hand, the signal (for which on can use the sum of the first $k=RG$ singular value of the stacked Page matrix of $f$ as a reasonable estimate) is assumed to obey a lower bound as the balanced spectra assumption states. We will add a discussion of how the various bounds can be stated as a function of the signal-to-noise ratio as suggested. Thank you for the great suggestion.\n\n> This paper is more typical in the fields of econometrics or statistics. The reviewer is not sure if it fits NeurIPS well.\n\nWhile we acknowledge the reviewer's concerns, we respectfully disagree that our paper does not fit NeurIPS. Time series papers, such as the examples cited in references [1-4], have been prominently featured at NeurIPS, demonstrating a range of technical flavors relevant to our work. [1] in particular is a paper of a very similar flavor (about mSSA for change point detection and is of similar technical depth). \n\n> Is it correct that we consider such a simplified multivariate time series modelling problem that each channel/variable (i.e. y_n) is driven an independent AR process (x_n)? Together with Assumption 2.1, we actually can model each channel signal y_n separately. This also allows the stack of Page matrix for an extension from the univariate to the multivariate.\nFollowing the previous question, the analysis for SAMOSSA for multivariate case has nothing essentially different from the univariate case, right? Maybe we missed some technical challenges.\n\nThe analysis of SAMoSSA encompasses studying both the deterministic and stochastic components of the time series.\n\nFor the deterministic part, the analysis of the multivariate case has its own challenges compared to that of the univariate case (i.e., SSA) and it characterizes how one can exploit assumption 2.1 to learn across different time series which results in the better scaling (w.r.t. N). \n\nFor the stochastic part, given that we assume that each channel is driven by its own stochastic process, the analysis of the AR parameter identification for the multivariate case is the same as the univariate case.  Considering a VAR-like model (or others), where there exists some dependence between the stationary processes $x_1$,..., $x_N$, is indeed an interesting direction of future work. \n\n> In the proof of out-of-sample forecasting error bound, what is the particular issue that our problem has compared to the analysis in [3]?\n\nThe difference between our setting and that of [3] is the noise settings, which raises the following difficulties when establishing the results of thm 4.4:\n1. **Spectral properties of AR noise.** given that the noise is not i.i.d., both reconstructing and forecasting the deterministic part ($f_n$) required establishing an upper bound on the operator norm of the Page matrix of AR processes; a result that may be of interest in its own right (see Lemma A.2). \n2. **Analysis of parameter identification of AR processes under arbitrary noise.** To establish the out-of-sample forecasting results for the observations $y_n$, we had to establish finite sample bounds for forecasting both $f_n$ and $x_n$ (whereas in the i.i.d. setting, $f_n$ was the only goal). This requires learning the AR parameters (system identification) under arbitrary and bounded observation noise (due to the reconstruction error).\n3. **Dependence between learned parameters and future noise terms.** The forecasting analysis in [3] benefits from the i.i.d. assumption as learned parameters remain independent of future “noise” terms $x_n$. For example,  the fact that $\\hat{\\beta}$ is independent of future noise terms $x_n(t), t > T$ is key when proving the forecasting error theorem in [3]. \n\n> We are curious about the comparative results. Indeed yours may be less accurate than DL SOTA, while yours has additional explainability. The performance gap might tell the cost for explainable modelling.\n\nBased on the collective feedback from reviewers, we added a new baseline to our experiments (DeepAR, see attached pdf), and we will attempt to add TRMF in our revised manuscript. Interestingly, DeepAR only outperforms SAMoSSA in the traffic dataset. \n\n\nWe again thank the reviewer for their comments and constructive feedback. We hope the reviewer will take these clarifications into account in their revised scores.\n \n**References**\n\n[1] Alanqary, Arwa, Abdullah Alomar, and Devavrat Shah. \"Change point detection via multivariate singular spectrum analysis.\" Advances in Neural Information Processing Systems 34 (2021): 23218-23230.\n\n[2] Mi, Xuelong, et al. \"BILCO: An Efficient Algorithm for Joint Alignment of Time Series.\" Advances in Neural Information Processing Systems 35 (2022): 36270-36281.\n\n[3] Liu, Yong, et al. \"Non-stationary transformers: Exploring the stationarity in time series forecasting.\" Advances in Neural Information Processing Systems 35 (2022): 9881-9893.\n\n[4] Wang, Zhiyuan, et al. \"Learning latent seasonal-trend representations for time series forecasting.\" Advances in Neural Information Processing Systems 35 (2022): 38775-38787.\n"}, {"review_id": "hpQE9sKJNM", "reviewer": "Reviewer_KKGZ", "summary": "This paper proposes a two-stage approach based on multivariate Singular Spectrum Analysis (mSSA) to  estimate the non-stationary components in a time series in the presence of a correlated stationary AR noise, which is subsequently estimated from the residual time series. Theoretical results on the performance of the algorithm in this novel setting, are established along with a finite-sample forecasting consistency bound. Empirical results demonstrate significant improvements in forecasting performance due to identification of the AR noise structure, across various benchmark datasets.", "questions": "The abbreviation SAMoSSA is used in the abstract and introduction without being explicitly defined (I can guess what it is through the title but not sure why the order is reversed). can you please state this explicitly?\n\nAssumption 2.1: Please provide some intuitive reasoning behind this assumption and what you mean by `fundamental' time series. \n\nIn Assumption 2.2, L<=T. Why does it become L<=sqrt(T) in line 168 on p.4?\n\nSection 3: The paper mentions that the proposed algorithm builds on the work of A. Agarwal, et al. [3]. It is not clear how the algorithm in Section 3 is different from what exists in the literature -it would be helpful to state this clearly. ", "limitations": "Yes, limitations of the approach are included in the final section. ", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 4, "strengths": "mSSA allows estimation of non-stationary deterministic components without domain knowledge or fine-tuning, however, cannot handle additive correlated stationary noise. The paper deals with the important problem of estimating non-stationary deterministic components in the presence of correlated stationary noise, through a unified approach via mSSA. Theoretical results for the mSSA are established beyond the i.i.d noise setting (specifically with AR noise). The paper also provides theoretical results on out-of-sample forecasting error for the proposed two-step algorithm. Overall, I found the contribution to be sufficiently novel and of high quality; addressing a crucial gap in the literature. ", "weaknesses": "Overall, the paper is well-written but is difficult to follow in some parts, and is lacking in conveying the key ideas. For example, presentation of the algorithm (Section 3): The details of the algorithm are presented well however, the main idea behind the approach is missing. A few sentences to convey the big picture behind the steps would have been very helpful.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nWe thank the reviewer for their positive and constructive feedback. Below we address the specific questions and comments they raise.\n\n> ... A few sentences to convey the big picture behind the steps would have been very helpful in the algorithm section.\n\nThank you, we will revise the algorithm description in the revised manuscript to convey the key ideas. \n\n> The abbreviation SAMoSSA is used in the abstract and introduction without being explicitly defined (I can guess what it is through the title but not sure why the order is reversed). can you please state this explicitly?\n\nThank you for this question. We think of SAMoSSA more like a pseudo-acronym, and not really an abbreviation. As you may probably surmised, \"SAMoSSA\" is a combination of the terms \"Stochastic\", \"Autoregressive\", and \"Multivariate Singular Spectrum Analysis\" (MSSA).   \n\n> Assumption 2.1: Please provide some intuitive reasoning behind this assumption and what you mean by `fundamental' time series.\n\nAnother way to state this assumption is: $f_n \\forall n \\in [N]$ is such that if you construct a matrix $M$ where $M_{ij}  = f_i(j)$, then rank($M$) $\\leq R$. That is, you can factorize the matrix M into “channel/time series” factors and temporal factors. These R temporal factors (in $\\mathbb{R}^{R\\times T}$) are what we call fundamental time series.   This assumption, while worded differently, is quite common in high dimensional time series analysis (e.g., see [1]).\n\n> In Assumption 2.2, L<=T. Why does it become L<=sqrt(T) in line 168 on p.4?\n\nThe $L\\leq\\sqrt{T}$ condition in the algorithm is used to simplify the exposition of the analysis but is not crucial for the results. More precisely, it's included to guarantee that the number of rows in the Page matrix is fewer than the number of columns.  For consistency, we will fix the wording of Assumption 2.2 to be $L\\leq\\sqrt{T}$.\n\n\n> The paper mentions that the proposed algorithm builds on the work of A. Agarwal, et al. [3]. It is not clear how the algorithm in Section 3 is different from what exists in the literature -it would be helpful to state this clearly.\n\nWe will revise the algorithm section to make this clearer. In brief, the first three steps in our algorithm (see figure 1) are the same as [3], but steps 4 and 5 are different. In particular, SAMoSSA, in addition to imputing and forecasting $f(\\cdot)$, also learns the AR process $x(\\cdot)$, and exploits that learned structure to forecast $x(t)$ for $t>T$.  This additional step, as illustrated in the empirical experiments, helps increase the forecast accuracy. \n\nWe would like to note though that a few technical challenges in the analysis arise with our model and algorithm compared to that of [3]. Specifically, \n- Given that we assume the noise process to be an AR process, most of the analysis presented in [3] breaks down. Specifically,  we now have to (1) establish spectral properties of the Page matrix of the correlated AR noise; a result that may be of interest in its own right (see Lemma A.2);  and (2) establish a forecasting error bound when the learned parameters (i.e. $\\hat{\\beta}$) and the stochastic terms (x_i’s) are correlated. \n- Further, as we attempt to learn the AR structure in the stochastic component, we had to improve over current results in system identification to accommodate for the case of arbitrary and bounded observation noise (due to the two-stage nature of the algorithm). In that sense, ours can be viewed as a robust generalization of previous work in AR system identification. \n\nWe again thank the reviewer for their comments and constructive feedback. \n\n**References**\n\n[1] Yu, Hsiang-Fu, Nikhil Rao, and Inderjit S. Dhillon. \"Temporal regularized matrix factorization for high-dimensional time series prediction.\" Advances in neural information processing systems 29 (2016).\n"}, {"review_id": "W1g4muUlHl", "reviewer": "Reviewer_ki9S", "summary": "This paper proposed SAMoSSA, a two-stage procedure that effectively handles mixtures of deterministic nonstationary and stationary AR processes with minimal model assumptions. The authors analyze SAMoSSA’s ability to estimate non-stationary components under stationary AR noise, the error rate of AR system identiﬁcation via OLS under observation errors, and a ﬁnite-sample forecast error analysis.", "questions": "1. The biggest difference between this paper and [3] lies in the different noise setting, i.e., this paper replaces the i.i.d noise with stationary AR noise. However, the theoretical derivation in this article is almost the same as it in [3]. What technical difficulties did the authors encounter during the proof process? \n\n2. How limited are  Assumptions 2.1 and 4.2 in practice?\n\n3. In experiment, more baselines are needed. \n\n4. How to verify that the dataset (especially the real dataset) satisfies the assumptions?\n\n", "limitations": "Not applicable", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. This paper is well-written and easy to follow.\n2. The theoretical results are solid. ", "weaknesses": "1. The theoretical contributions are somewhat incremental.\n2. The experiment are insufficient.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive and constructive feedback. Below we address the specific questions and comments they raise.\n\n> The biggest difference between this paper and [3] lies in the different noise setting, i.e., this paper replaces the i.i.d noise with stationary AR noise. However, the theoretical derivation in this article is almost the same as it in [3]. What technical difficulties did the authors encounter during the proof process?\n\nThe different noise settings raised the following three difficulties\n- **Establishing spectral properties of AR noise.** Given that the noise is not i.i.d., both reconstructing and forecasting the deterministic part ($f_n(t)$) using mSSA required establishing an upper bound on the operator norm of the “Page” matrix of AR processes; a result that may be of interest in its own right (see Lemma A.2). \n- **Analysis of parameter identification of AR processes under arbitrary noise.** To establish the out-of-sample forecasting results for the observations $y_n$, we had to establish a finite sample bound for forecasting both $f_n$ and $x_n$ (whereas in the i.i.d. noise setting, $f_n$ was the only goal). This requires learning the AR parameters (system identification) under  arbitrary and bounded observation noise.\n- **Dependence between learned parameters and future noise terms.** The forecasting analysis in [3] benefits from the iid assumption. For example,  the fact that $\\hat{\\beta}$ is independent of future noise terms $x_n(t), t > T$ is key to cancel a few of the terms when proving the forecasting error theorem. \n\n> How limited are Assumptions 2.1 and 4.2 in practice?\nHow to verify that the dataset (especially the real dataset) satisfies the assumptions?\n\t\nAssumption 2.1 implies that if you construct a matrix M where $M_{ij}  = f_i(j)$, then rank($M$) $\\leq R$. We would argue that this assumption is actually common in high dimensional time series analysis (see [1]). Indeed, TRMF, and also mSSA, has been shown to perform very well in practice on a variety of high dimensional multivariate time series which suggests that the assumption is not limiting in practice. \n\nTo verify whether the spatio-temporal model we assume (assumptions 2.1 and 2.2) hold in practice, one can use the diagnostic test for the Spatio-Temporal Model is furnished in [2]. This test verifies whether mSSA is likely to succeed based on the spectral properties of the observed matrix.  Specifically, the test measures the effective rank of the Page matrix associated with the multivariate time series with parameter $L \\sim \\sqrt{NT}$. If the effective rank does not scale much slower than $L$ then mSSA is unlikely to be effective. We will add a reference and discussion of this test in the revised version. \n\nAssumption 4.2 will hold in practice if the underlying deterministic part of the time series (i.e., $f_i(\\cdot)$) does not change, assuming that we have observed enough samples to capture the time series complexity. We believe this assumption is both necessary and natural for our settings; Note that typically, to establish generalization error, modern statistical estimators assume the data-generating process to be i.i.d. Herein, and in line with the work in [3], we make a less restrictive assumption and rely on a purely linear algebraic condition that is natural with the proposed model. \n\n\n> In experiment, more baselines are needed.\nAs other reviewers correctly point out,  the theoretical characterization of multi-stage learning algorithms in this setup is missing in the literature, and this is the main focus and deliverable of our paper.  That said, based on the collective feedback from reviewers, we added a new baseline to our experiments (DeepAR, see attached pdf), and we will attempt to add TRMF [1] in our revised manuscript. Interestingly, DeepAR only outperforms SAMoSSA in the traffic dataset. \n\t\nWe again thank the reviewer for their comments and constructive feedback. We hope the reviewer will take these clarifications into account in their revised scores.\n\n**References**\n\n[1] Yu, Hsiang-Fu, Nikhil Rao, and Inderjit S. Dhillon. \"Temporal regularized matrix factorization for high-dimensional time series prediction.\" Advances in neural information processing systems 29 (2016).\n\n[2] Agarwal, Anish, Abdullah Alomar, and Devavrat Shah. \"On multivariate singular spectrum analysis and its variants.\" ACM SIGMETRICS Performance Evaluation Review 50.1 (2022): 79-80.\n\n[3] Agarwal, Anish, Devavrat Shah, and Dennis Shen. \"On Model Identification and Out-of-Sample Prediction of Principal Component Regression: Applications to Synthetic Controls.\" arXiv preprint arXiv:2010.14449 (2020).\n"}, {"review_id": "ZqJajqknuS", "reviewer": "Reviewer_xDu7", "summary": "The authors propose SAMoSSA, an algorithm that combines deterministic trend estimation via mSSA with estimation of an autoregressive component of a time series.  They provide error rates for trend estimation, estimation of the AR coefficients, as well as the prediction error.   In addition, they consider real and simulated examples and show that the proposed method offers very competitive performance.", "questions": "The more general ARMA model is first discussed in the introduction.  What challenges do estimating the MA components pose for the analysis?", "limitations": "The authors do acknowledge some limitations with the studied framework, including assuming stationarity of the stochastic component.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- While the time series literature is vast and well-developed, it does appear that theoretical results that analyze the effect of trend estimation on AR estimates and the prediction error were quite limited.  Detrending a time series and then modeling using ARMA or VAR is a common routine so this is a practically important problem.  \n- The derived error bounds appear to be novel.  \n- The paper is well-written.  Assumptions and implications of theorems are for the most part explained well.", "weaknesses": "- The authors could do a more thorough literature review when it comes to inference for trends.  For example, there is a statistics literature on testing for trends. See Chen and Wu (2019): Testing for Trends in High-Dimensional Time Series and references therein.  \n- While it is a nice result, the implications for real-world data are questionable due to the conditions on the trend as well as the assumed autoregressive structure.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive feedback and valuable questions. Below we address the specific questions and comments they raise.\n\n> The more general ARMA model is first discussed in the introduction. What challenges do estimating the MA components pose for the analysis?\n\nThank you for this great question. If an ARMA process is stable, then it can be represented as an infinite MA process whose coefficients are absolutely summable. Hence, we believe that the Non-Stationary deterministic component can be learned (with a similar bound to the one in Theorem 4.1). However, more careful analysis is needed to determine how well we can learn the ARMA parameters and how well we can forecast. This is indeed an interesting venue for future work. We will add a discussion about ARMA in the revision. \n\n\n> While it is a nice result, the implications for real-world data are questionable due to the conditions on the trend as well as the assumed autoregressive structure.\n\nWe believe that the model we provide is rich. In terms of the deterministic part, many standard functions that model time series dynamics satisfy assumptions 2.1-2.2 either exactly or approximately.  Indeed,  as [3] shows, these functions include any finite sum of products of harmonics, low-degree polynomials, and exponential functions. \n\nIn terms of the stochastic part, it is worth noting that, in theory,  stationary stochastic processes can be approximated by AR processes suggesting that AR processes are also quite rich. However, there are valid concerns about the limitations we set in our model. In future work, we aim to also explore when the stochastic part is:\n1. an ARMA structure. \n2. an integrated process (e.g., ARIMA), as we discuss in the limitation\n\nThese are indeed great extensions to further increase the richness of the model. That being said, the model we’re considering now is indeed rich. The improvements of SAMoSSA over mSSA and other baselines, as well as the analysis we provide in the rebuttal, strongly suggest the richness and practicality of our model. \n\nWe again thank the reviewer for their comments and constructive feedback. We hope the reviewer will take these clarifications into account in their revised scores.\n\n"}]

Author rebuttal:
We thank all the reviewers for constructive feedback. Here's a succinct highlight of our paper's key contributions, beyond our individual responses to reviewers:

The main contribution of this paper is to showcase the effectiveness of a simple multi-stage algorithm in time series forecasting. For which we do the following:
1. To establish its effectiveness in theory, we had to extend prior work in the following ways:
- We extend the analysis of mSSA in [3] to accommodate the case of correlated autoregressive noise, a prevalent noise structure in time series analysis. In particular, we establish that one can effectively estimate and forecast the *deterministic* part of the signal under correlated autoregressive noise. 
- We extend the analysis of AR parameter estimation in [14] to accommodate the case when the process is observed under arbitrary bounded noise. Our results can be thought of as a robust generalization of [14] which derives similar results but without any observation noise. Thus, we establish that one can effectively forecast the *stochastic* part of the signal. 

That is, we overcame two key challenges in analyzing the proposed multi-stage algorithm, significantly building upon recent prior work..

2. We showcase that empirically, through limited but representative datasets and baselines, that SAMoSSA outperforms mSSA. This showcases the effectiveness of the multi-stage process we propose and that the model we proposed is indeed reasonable (see Figure 1 in the rebuttal for further evidence). We believe that further extensions (e.g., considering VAR-like process for the stochastic component) can help achieve even better performances -- which we will consider in future work.




Meta-review: This paper is about multivariate time series analysis. The authors consider N-dimensional time series y=(y_n)_{n=1}^N where each coordinate is modelled as a sum of a non-stationary deterministic component (f_n) and an AR(p_n) process x_n, as written in (1)-(2) with Assumption 2.1-2.4, 4.1-4.2. The goal is to estimate f_n, the parameters of x_n and to forecast the observation y_n. To address these challenges they design a novel two-stage algorithm based on multivariate singular spectrum analysis (mSSA). The resulting SAMoSSA technique is accompanied with finite-sample l_2 estimation error guarantee of the non-stationary component (Theorem 4.1) and the AR parameters (Theorem 4.2), and that of the forecasting (Theorem 4.4). Numerical results on both toy and 3 real-world datasets compared to existing baselines underpin the efficiency of the method.

Time series are omnipresent in data science. Developing principled techniques to analyze them is of fundamental interest to the machine learning community. The authors provide new valuable tools in this context, in a nicely written and well-organized manuscript.
Decision: Accept (poster)