# Foundations of Sequence-to-Sequence Modeling for Time Series

Vitaly Kuznetsov<EMAIL_ADDRESS>
Google Research, New York, NY
Zelda Mariet<EMAIL_ADDRESS>
Massachusetts Institute of Technology, Cambridge, MA Authors are in
alphabetical order.

###### Abstract

The availability of large amounts of time series data, paired with the
performance of deep-learning algorithms on a broad class of problems, has
recently led to significant interest in the use of sequence-to-sequence models
for time series forecasting. We provide the first theoretical analysis of this
time series forecasting framework. We include a comparison of sequence-to-
sequence modeling to classical time series models, and as such our theory can
serve as a quantitative guide for practitioners choosing between different
modeling methodologies.

##  1 Introduction

Time series analysis is a critical component of real-world applications such
as climate modeling, web traffic prediction, neuroscience, as well as
economics and finance. We focus on the fundamental question of time series
forecasting. Specifically, we study the task of forecasting the next $\ell$
steps of an $m$-dimensional time series ${\bf Y}$, where $m$ is assumed to be
very large. For example, in climate modeling, $m$ may correspond to the number
of locations at which we collect historical observations, and more generally
to the number of sources which provide us with time series.

Often, the simplest way to tackle this problem is to approach it as $m$
separate tasks, where for each of the $m$ dimensions we build a model to
forecast the univariate time series corresponding to that dimension. Auto-
regressive and state-space models [8, 4, 6, 5, 12], as well as non-parametric
approaches such as RNNs [3], are often used in this setting. To account for
correlations between different time series, these models have also been
generalized to the multivariate case [23, 24, 35, 13, 14, 1, 2, 35, 2, 38, 31,
41, 25]. In both univariate and multivariate settings, an observation at time
$t$ is treated as a single sample point, and the model tries to capture
relations between observations at times $t$ and $t+1$. Therefore, we refer to
these models as _local_.

In contrast, an alternative methodology based on treating $m$ univariate time
series as $m$ samples drawn from some unknown distribution has also gained
popularity in recent years. In this setting, each of the $m$ dimensions of
${\bf Y}$ is treated as a separate example and a single model is learned from
these $m$ observations. Given $m$ time series of length $T$, this model learns
to map past vectors of length $T-\ell$ to corresponding future vectors of
length $\ell$. LSTMs and RNNs [15] are a popular choice of model class for
this setup [9, 11, 42, 22, 21, 43].111Sequence-to-sequence models are also
among the winning solutions in the recent time series forecasting competition:
https://www.kaggle.com/c/web-traffic-time-series-forecasting. Consequently, we
refer to this framework as _sequence-to-sequence modeling_.

While there has been progress in understanding the generalization ability of
local models [40, 27, 28, 29, 17, 18, 19, 20, 44], to the best of our
knowledge the generalization properties of sequence-to-sequence modeling have
not yet been studied, raising the following natural questions:

  * –

What is the generalization ability of sequence-to-sequence models and how is
it affected by the statistical properties of the underlying stochastic
processes (e.g. non-stationarity, correlations)?

  * –

When is sequence-to-sequence modeling preferable to local modeling, and vice
versa?

We provide the first generalization guarantees for time series forecasting
with sequence-to-sequence models. Our results are expressed in terms of
simple, intuitive measures of non-stationarity and correlation strength
between different time series and hence explicitly depend on the key
components of the learning problem.

Table 1: Summary of local, sequence-to-sequence, and hybrid models. Learning Model | Training set | Hypothesis | Example
---|---|---|---
UniVar. local | ${\bf Z}_{i}=\\{(Y_{t-1-p}^{t-1}(i),Y_{t}(i)):p\leq t\leq T\\}$ | $h_{i}:\mathcal{Y}^{p}\to\mathcal{Y}$ | ARIMA
MultiVar. local | ${\bf Z}=\\{(Y_{t-1-p}^{t-1}(\cdot),Y_{t}(\cdot)):p\leq t\leq T\\}$ | $h:\mathcal{Y}^{m\times p}\to\mathcal{Y}^{m}$ | VARMA
Seq-to-seq | ${\bf Z}=\\{(Y_{1}^{T-1}(i),Y_{T}(i)):1\leq i\leq m\\}$ | $h:\mathcal{Y}^{T}\to\mathcal{Y}$ | Neural nets
Hybrid | ${\bf Z}=\\{(Y_{t-p}^{t-1}(i),Y_{t}(i)):1\leq i\leq m,p\leq t\leq T\\}$ | $h:\mathcal{Y}^{p}\to\mathcal{Y}$ | Neural nets

test inputtest targetinputtarget (a) The local model trains each
$h_{\text{loc},i}$ on time series $Y(i)$ split into multiple (partly
overlapping) examples.

test inputtest targetinputtarget (b) The sequence-to-sequence trains
$h_{\text{s2s}}$ on $m$ time series split into (input, target) pairs.

Figure 1: Local and sequence-to-sequence splits of a one dimensional time
series into training and test examples.

We compare our generalization bounds to guarantees for local models and
identify regimes under which one methodology is superior to the other.
Therefore, our theory may also serve as a quantitative guide for a
practitioner choosing the right modeling approach.

The rest of the paper is organized as follows: in Section 2, we formally
define sequence-to-sequence and local modeling. In Section 3, we define the
key tools that we require for our analysis. Generalization bounds for
sequence-to-sequence models are given in Section 4. We compare sequence-to-
sequence and local models in Section 5. Section 6 concludes this paper with a
study of a setup that is a hybrid of the local and sequence-to-sequence
models.

##  2 Sequence-to-sequence modeling

We begin by providing a formal definition of sequence-to-sequence modeling.
The learner receives a multi-dimensional time series ${\bf
Y}\in\mathcal{Y}^{m\times T}$ which we view as $m$ time series of same length
$T$. We denote by $Y_{t}(i)$ the value of the $i$-th time series at time $t$
and write $Y_{a}^{b}(i)$ to denote the sequence
$(Y_{a}(i),Y_{a+1}(i),\ldots,Y_{b}(i))$. Similarly, we let
$Y_{t}(\cdot)=(Y_{t}(1),\ldots,Y_{t}(m))$ and
$Y_{a}^{b}(\cdot)=(Y_{a}(\cdot),\ldots,Y_{b}(\cdot))$. In particular, ${\bf
Y}\equiv Y_{1}^{T}(\cdot)$. In addition, the sequence $Y_{1}^{T-1}(\cdot)$ is
of a particular importance in our analysis and we denote it by
$\mathbf{Y}^{\prime}$.

The goal of the learner is to predict $Y_{T+1}(\cdot)$.222We are often
interested in long term forecasting, i.e. predicting $Y_{T+1}^{T+\ell}(\cdot)$
for $\ell\geq 1$. For simplicity, we only consider the case of $\ell=1$.
However, all our results extend immediately to $\ell\geq 1$. We further assume
that our input ${\bf Y}$ is partitioned into a training set of $m$ examples
${\bf Z}=\\{Z_{1},\ldots,Z_{m}\\}$, where each
$Z_{i}=(Y_{1}^{T-1}(i),Y_{T}(i))\in\mathcal{Y}^{T}$. The learner’s objective
is to select a hypothesis $h:\mathcal{Y}^{T}\to\mathcal{Y}$ from a given
hypothesis set $\mathcal{H}$ that achieves a small generalization error:

$\mathcal{L}(h\mid{\bf
Y})=\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}_{\mathcal{D}}\left[L(h(Y_{1}^{T}(i)),Y_{T+1}(i))\mid{\bf
Y}\right],$

where $L\colon\mathcal{Y}\times\mathcal{Y}\rightarrow[0,M]$ is a
bounded333Most of the results in this paper can be straightforwardly extended
to unbounded case assuming $Y$ is sub-Gaussian. loss function and
$\mathcal{D}$ is the distribution of $Y_{T+1}$ conditioned on the past ${\bf
Y}$.

In other words, the learner seeks a hypothesis $h$ that maps sequences of past
$Y_{1}(i),\ldots,Y_{T}(i)$ values to sequences of future values
$Y_{T+1}(i),\ldots,Y_{T+\ell}(i)$, justifying our choice of “sequence-to-
sequence” terminology.444In practice, each $Z_{i}$ may start at a different,
arbitrary time $t_{i}$, and may furthermore include some additional features
$X_{i}$, i.e. $Z_{i}=(Y_{t_{i}}^{T-1}(i),X_{i},Y_{T}(i))$. Our results can be
extended to this case as well using an appropriate choice of the hypothesis
set. Incidentally, the machine translation problem studied by Sutskever et al.
[36] under the same name represents a special case of our problem when
sequences (sentences) are independent and data is stationary. In fact, LSTM-
based approaches used in aforementioned translation problem are also common
for time series forecasting [9, 11, 42, 22, 21, 43]. However, feed-forward NNs
have also been successfully applied in this framework [34] and in practice,
our definition allows for any set of functions $\mathcal{H}$ that map input
sequences to output sequences. For instance, we can train a feed-forward NN to
map $Y_{1}^{T-1}(i)$ to $Y_{T}(i)$ and at inference time use $Y_{2}^{T}(i)$ as
input to obtain a forecast for $Y_{T}(i)$.555As another example, a runner-up
in the Kaggle forecasting competition (https://www.kaggle.com/c/web-traffic-
time-series-forecasting) used a combination of boosted decision trees and
feed-forward networks, and as such employs the sequence-to-sequence approach.

We contrast sequence-to-sequence modeling to _local_ modeling, which consists
of splitting each time series $Y(i)$ into a training set ${\bf
Z}_{i}=\\{Z_{i,1},\ldots,Z_{i,T}\\}$, where
$Z_{i,t}=(Y_{t-p}^{t-1}(i),Y_{t}(i))$ for some $p\in\mathbb{N}$, then learning
a separate hypothesis $h_{i}$ for each ${\bf Z}_{i}$. Each $h_{i}$ models
relations between observations that are close in time, which is why we refer
to this framework as local modeling. As in sequence-to-sequence modeling, the
goal of a local learner is to achieve a small generalization error for
$h_{\text{loc}}=(h_{1},\ldots,h_{m})$, given by:

$\mathcal{L}(h_{\text{loc}}\mid{\bf
Y})=\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}_{\mathcal{D}}\left[L(h_{i}(Y_{T-p}^{T}(i)),Y_{T+1}(i))\mid{\bf
Y}\right].$

In order to model correlations between different time series, it is also
common to split ${\bf Y}$ into one single set of multivariate examples ${\bf
Z}=\\{Z_{1},\ldots,Z_{T}\\}$, where
$Z_{i}=(Y_{t-p}^{t-1}(\cdot),Y_{t}(\cdot))$, and to learn a single hypothesis
$h$ that maps $\mathcal{Y}^{m\times p}\to\mathcal{Y}^{m}$. As mentioned
earlier, we consider this approach a variant of local modeling, since $h$ in
this case again models relations between observations that are close in time.

Finally, _hybrid_ or _local sequence-to-sequence_ models, which interpolate
between local and sequence-to-sequence approaches, have also been considered
in the literature [43]. In this setting, each local example is split across
the temporal dimension into smaller examples of length $p$, which are then
used to train a single sequence-to-sequence model $h$. We discuss bounds for
this specific case in Section 6.

Our work focuses on the statistical properties of the sequence-to-sequence
model. We provide the first generalization bounds for sequence-to-sequence and
hybrid models, and compare these to similar bounds for local models, allowing
us to identify regimes in which one methodology is more likely to succeed.

Aside from learning guarantees, there are other important considerations that
may lead a practitioner to choose one approach over others. For instance, the
local approach is trivially parallelizable; on the other hand, when additional
features $X_{i}$ are available, sequence-to-sequence modeling provides an
elegant solution to the _cold start problem_ in which at test time we are
required to make predictions on time series for which no historical data is
available.

##  3 Correlations and non-stationarity

In the standard supervised learning scenario, it is common to assume that
training and test data are drawn i.i.d. from some unknown distribution.
However, this assumption does not hold for time series, where observations at
different times as well as across different series may be correlated.
Furthermore, the data-generating distribution may also evolve over time.

These phenomena present a significant challenge to providing guarantees in
time series forecasting. To quantify non-stationarity and correlations, we
introduce the notions of mixing coefficients and discrepancy, which are
defined below.

The final ingredient we need to analyze sequence-to-sequence learning is the
Rademacher complexity $\mathfrak{R}_{m}(\mathcal{F})$ of a family of functions
$\mathcal{F}$ on a sample of size $m$, which has been previously used to
characterize learning in the i.i.d. setting [16, 30]. In App. A, we include a
brief discussion of its properties.

###  3.1 Expected mixing coefficients

To measure the strength of dependency between time series, we extend the
notion of $\beta$-mixing coefficients [7] to _expected $\beta$-mixing
coefficients_, which are a more appropriate measure of correlation in
sequence-to-sequence modeling.

######  Definition 1 (Expected $\beta_{\text{{s2s}}}$ coefficients).

Let $i,j\in[m]\triangleq\\{1,\ldots,m\\}$. We define

$\displaystyle\beta_{\text{{s2s}}}(i,j)=$ $\displaystyle~{}\mathbb{E}_{{\bf
Y}^{\prime}}\Big{[}\|P(Y_{T}(i)|{\bf Y}^{\prime})P(Y_{T}(j)|{\bf
Y}^{\prime})-P(Y_{T}(i),Y_{T}(j)|{\bf Y}^{\prime})\|_{TV}\Big{]},$

where $TV$ denotes the total variations norm. For a subset $I\subseteq[m]$, we
define

$\beta_{\text{{s2s}}}(I)=\sup_{i,j\in C}\beta_{\text{{s2s}}}(i,j).$

The coefficient $\beta_{\text{{s2s}}}(i,j)$ captures how close $Y_{T+1}(i)$
and $Y_{T+1}(j)$ are to being independent, given ${\bf Y}^{\prime}$ (and
averaged over all realizations of ${\bf Y}^{\prime}$). We further study these
coefficients in Section 4, where we derive explicit upper bounds on expected
$\beta_{\text{{s2s}}}$-mixing coefficients for various standard classes of
stochastic processes, including spatio-temporal and hierarchical time series.

We also define the following related notion of $\bar{\beta}$-coefficients.

######  Definition 2 (Unconditional $\bar{\beta}$ coefficients).

Let $i,j\in[m]\triangleq\\{1,\ldots,m\\}$. We define

$\displaystyle\bar{\beta}(i,j)\\!=$
$\displaystyle\|\Pr(Y_{1}^{T}(i),Y_{1}^{T}(j))-\Pr(Y_{1}^{T}(i))\Pr(Y_{1}^{T}(j))\|_{TV}$
$\displaystyle\bar{\beta}^{\prime}(i,j)\\!=$
$\displaystyle\|\Pr(Y_{1}^{T-1}(i),Y_{1}^{T-1}(j))-\Pr(Y_{1}^{T-1}(i))\Pr(Y_{1}^{T-1}(j))\|_{TV}$

and as before, for a subset $I$ of $[m]$, write $\bar{\beta}(I)=\sup_{i,j\in
I}\bar{\beta}(i,j)$ (and similarly for $\bar{\beta}^{\prime}$).

Note that $\beta_{\text{{s2s}}}$ coefficients measure the strength of
dependence between time series conditioned on the history observed so far,
while $\bar{\beta}$ coefficients measure the (unconditional) strength of
dependence between time series. The following result relates these two
notions.

######  Lemma 1.

For $\bar{\beta}$ (and $\bar{\beta}^{\prime}$ similarly), we have the
following upper bound:

$\displaystyle\bar{\beta}(i,j)\leq$
$\displaystyle\beta_{\text{{s2s}}}(i,j)+\mathbb{E}_{{\bf
Y}^{\prime}}\Big{[}\text{Cov}\Big{(}\Pr(Y_{T}(i)\mid{\bf
Y}^{\prime}),\Pr(Y_{T}(j)\mid{\bf Y}^{\prime})\Big{)}\Big{]}$

The proof of this result (as well as all other proofs in this paper) is
deferred to the supplementary material.

Finally, we require the notion of _tangent collections_ , within which time
series are independent.

######  Definition 3 (Tangent collection).

Given a collection of time series $C=\\{Y(1),\ldots,Y(c)\\}$, we define the
tangent collection $\widetilde{C}$ as
$\\{\widetilde{Y}(1),\ldots,\widetilde{Y}(c)\\}$ such that $\widetilde{Y}(i)$
is drawn according to the marginal $\Pr(Y(i))$ and such that
$\widetilde{Y}(i)$ and $\widetilde{Y}(i^{\prime})$ are independent for $i\neq
i^{\prime}$.

The notion of tangent collections, combined with mixing coefficients, allows
us to reduce the analysis of correlated time series in $C$ to the analysis of
independent time series in $\widetilde{C}$ (see Prop. 6 in the appendix).

###  3.2 Discrepancy

Various notions of discrepancy have been previously used to measure the non-
stationarity of the underlying stochastic processes with respect to the
hypothesis set $\mathcal{H}$ and loss function $L$ in the analysis of local
models [18, 44]. In this work, we introduce a notion of discrepancy
specifically tailored to sequence-to-sequence modeling scenario, taking into
account both the hypothesis set and the loss function.

######  Definition 4 (Discrepancy).

Let $\mathcal{D}$ be the distribution of $Y_{T+1}$ conditioned on ${\bf Y}$
and let $\mathcal{D}^{\prime}$ be the distribution of $Y_{T}$ conditioned on
$\mathbf{Y}^{\prime}$. We define the discrepancy $\Delta$ as
$\Delta=\sup_{h\in\mathcal{H}}\left|\mathcal{L}(h\mid{\bf
Y})-\mathcal{L}(h\mid\mathbf{Y}^{\prime})\right|$ where $\mathcal{L}(h\mid{\bf
Y}^{\prime})=\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}_{\mathcal{D}^{\prime}}\left[L(h(Y_{1}^{T-1}(i)),Y_{T}(i))\mid{\bf
Y}\right]$.

The discrepancy forms a pseudo-metric on the space of probability
distributions and can be completed to a Wasserstein metric (by extending
$\mathcal{H}$ to all Lipschitz functions). This also immediately implies that
the discrepancy can be further upper bounded by the $l_{1}$-distance and by
relative entropy between conditional distributions of $Y_{T}$ and $Y_{T+1}$
(via Pinsker’s inequality). However, unlike these other divergences, the
discrepancy takes into account both the hypothesis set and the loss function,
making it a finer measure of non-stationarity.

However, the most important property of the discrepancy is that it can be
upper bounded by the related notion of _symmetric discrepancy_ , which can be
estimated from data.

######  Definition 5 (Symmetric discrepancy).

We define the symmetric discrepancy $\Delta_{s}$ as

$\displaystyle\Delta_{s}=$
$\displaystyle\frac{1}{m}\sup_{h,h^{\prime}\in\mathcal{H}}\Big{|}\sum\nolimits_{i=1}^{m}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))\Big{|}.$

######  Proposition 1.

Let $\mathcal{H}$ be a hypothesis space and let $L$ be a bounded loss function
which respects the triangle inequality. Let $h\in\mathcal{H}$ be any
hypothesis. Then, $\Delta\leq\Delta_{s}+\mathcal{L}(h\mid{\bf
Y})+\mathcal{L}(h\mid{\bf Y}^{\prime}).$

We do not require test labels to evaluate $\Delta_{s}$. Since $\Delta_{s}$
only depends on the observed data, $\Delta_{s}$ can be computed directly from
samples, making it a useful tool to assess the non-stationarity of the
learning problem.

Another useful property of $\Delta_{s}$ is that, for certain classes of
stochastic processes, we can provide a direct analysis of this quantity.

######  Proposition 2.

Let $I_{1},\cdots,I_{k}$ be a partition of $\\{1,\ldots,m\\}$,
$C_{1},\ldots,C_{k}$ be the corresponding partition of ${\bf Y}$ and
$C^{\prime}_{1},\ldots,C^{\prime}_{k}$ be the corresponding partition of ${\bf
Y}^{\prime}$. Write $c=\min_{j}|C_{j}|$, and define the _expected discrepancy_

$\displaystyle\Delta_{e}=$
$\displaystyle\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\mathbb{E}_{Y}[L(h(Y_{1}^{T}),h^{\prime}(Y_{1}^{T}))]-\mathbb{E}_{Y}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))]\Big{]}.$

Then, writing $\mathfrak{R}$ the Rademacher complexity (see Appendix A) we
have with probability $1-\delta$,

$\displaystyle\Delta_{s}\leq$
$\displaystyle\Delta_{e}+{\max\Big{(}\max_{j}\mathfrak{R}_{|C_{j}|}(\widetilde{C}^{\prime}_{j}),\max_{j}\mathfrak{R}_{|C_{j}|}(\widetilde{C}_{j})\Big{)}}+\sqrt{\frac{1}{2c}\log\frac{2k}{\delta-\sum_{j}(|I_{j}|-1)[\bar{\beta}(I_{j})+\bar{\beta}^{\prime}(I_{j})]}}.$

The expected discrepancy $\Delta_{e}$ can be computed analytically for many
classes of stochastic processes. For example, for stationary processes, we can
show that it is negligible. Similarly, for covariance-stationary666Recall that
a process $X_{1},X_{2},\ldots$ is called stationary if for any $l,k,m$, the
distributions of $(X_{k},\ldots,X_{k+l})$ and $(X_{k+m},\ldots,X_{k+m+l})$ are
the same. Covariance stationarity is a weaker condition that requires that
$\mathbb{E}[X_{k}]$ be independent of $k$ and that
$\mathbb{E}[X_{k}X_{m}]=f(k-m)$ for some $f$. processes with linear hypothesis
sets and the squared loss function, the discrepancy is once again negligible.
These examples justify our use of the discrepancy as a natural measure of non-
stationarity. In particular, the covariance-stationary example highlights that
the discrepancy takes into account not only the distribution of the stochastic
processes but also $\mathcal{H}$ and $L$.

######  Proposition 3.

If $Y(i)$ is stationary for all $1\leq i\leq m$, and $\mathcal{H}$ is a
hypothesis space such that $h\in\mathcal{H}:\mathcal{Y}^{T-1}\to\mathcal{Y}$
(i.e. the hypotheses only consider the last $T-1$ values of $Y$), then
$\Delta_{e}=0$.

######  Proposition 4.

If ${\bf Y}$ is covariance stationary for all $1\leq i\leq m$, $L$ is the
squared loss, and $\mathcal{H}$ is a linear hypothesis space $\\{x\to w\cdot
x\mid\|w\|\in\mathbb{R}^{p}\leq\Lambda\\}$, $\Delta_{e}=0$.

Another insightful example is the case when $\mathcal{H}=\\{h\\}$: then,
$\Delta=0$ even if ${\bf Y}$ is non-stationary, which illustrates that
learning is trivial for trivial hypothesis sets, even in non-stationary
settings.

The final example that we consider in this section is the case of non-
stationary periodic time series. Remarkably, we show that the discrepancy is
still negligible in this case provided that we observe all periods with equal
probability.

######  Proposition 5.

If the $Y(i)$ are periodic with period $p$ and the observed starting time of
each $Y(i)$ is distributed uniformly at random in $[p]$, then $\Delta_{e}=0$.

##  4 Generalization bounds

We now present our generalization bounds for time series prediction with
sequence-to-sequence models. We write $\mathcal{F}=\\{L\circ
h:~{}h\in\mathcal{H}\\}$, where $f=L\circ h$ is the loss of hypothesis $h$
given by $f(h,Z_{i})=L(h(Y_{1}^{T-1}(i)),Y_{T})$. To obtain bounds on the
generalization error $\mathcal{L}(h\mid{\bf Y})$, we study the gap between
$\mathcal{L}(h\mid{\bf Y})$ and the empirical error $\widehat{\mathcal{L}}(h)$
of a hypothesis $h$, where

$\widehat{\mathcal{L}}(h)=\frac{1}{m}\sum_{i=1}^{m}f\big{(}h,Z_{i}\big{)}.$

That is, we aim to give a high probability bound on the supremum of the
empirical process $\Phi({\bf Y})=\sup_{h}[\mathcal{L}(h\mid{\bf
Y})-\widehat{\mathcal{L}}(h)]$. We take the following high-level approach: we
first partition the training set ${\bf Z}$ into $k$ collections
$C_{1},\ldots,C_{k}$ such that within each collection, correlations between
different time series are as weak as possible. We then analyze each collection
$C_{j}$ by comparing the generalization error of sequence-to-sequence learning
on $C_{j}$ to the sequence-to-sequence generalization error on the _tangent_
collection $\widetilde{C}_{j}$.

######  Theorem 4.1.

Let $C_{1},\ldots,C_{k}$ form a partition of the training input ${\bf Z}$ and
let $I_{j}$ denote the set of indices of time series that belong to $C_{j}$.
Assume that the loss function $L$ is bounded by 1. Then, we have for any
$\delta>\sum_{j}(|I_{j}|-1)\beta(I_{j})$, with probability $1-\delta$,

$\displaystyle\Phi$ $\displaystyle({\bf
Y})\leq\max_{j}\left[\widehat{}\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})\right]+\Delta+\frac{1}{\sqrt{2\min_{j}|I_{j}|}}\sqrt{\log\left(\frac{k}{\delta-\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})}\right)}.$

Theorem 4.1 illustrates the trade-offs that are involved in sequence-to-
sequence learning for time series forecasting. As
$\sum_{j}(|I_{j}|-1)\beta(I_{j})$ is a function of $m$, we expect it to
decrease as $m$ grows (_i.e._ more time series we have), allowing for smaller
$\delta$ as $m$ increases.

Assuming that the $C_{j}$ are of the same size, if $\mathcal{H}$ is a
collection of neural networks of bounded depth and width then
$\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})=\mathcal{O}\Big{(}\sqrt{{kT}/{m}}\Big{)}$
(see Appendix A). Therefore,

${\mathcal{L}(h\mid{\bf
Y})\leq\widehat{\mathcal{L}}(h)+\Delta}+\mathcal{O}\Big{(}\sqrt{\tfrac{kT}{m}}\Big{)}$

with high probability uniformly over $h\in\mathcal{H}$, provided that
$\tfrac{m}{k}\sum_{j=1}^{k}\beta_{\text{{s2s}}}(I_{j})=o(1)$. This shows that
extremely high-dimensional ($m\gg 1$) time series are beneficial for sequence-
to-sequence models, whereas series with a long histories $T\gg m$ will
generally not benefit from sequence-to-sequence learning. Note also that
correlations in data reduce the effective sample size from $m$ to $m/k$.

Furtermore, Theorem 4.1 indicates that balancing the complexity of the model
(e.g. depth and width of a neural net) with the fit it provides to the data is
critical for controlling both the discrepancy and Rademacher complexity terms.
We further illustrate this bound with several examples below.

###  4.1 Independent time series

We begin by considering the case where all dimensions of ${\bf Y}$ are
independent. Although this may seem a restrictive assumption, it arises in a
variety of applications: in neuroscience, different dimensions may represent
brain scans of different patients; in reinforcement learning, they may
correspond to different trajectories of a robotic arm.

######  Theorem 4.2.

Let $\mathcal{H}$ be a given hypothesis space with associated function family
$\mathcal{F}$ corresponding to a loss function $L$ bounded by 1. Suppose that
all dimensions of ${\bf Y}$ are independent and let $I_{1}=[m]$; then
$\beta(I_{1})=0$ and so for any $\delta>0$, with probability at least
$1-\delta$ and for any $h\in\mathcal{H}$:

$\mathcal{L}(h|{\bf
Y})\leq\widehat{\mathcal{L}}(h)+2\mathfrak{R}_{m}(\mathcal{F})+\Delta+\sqrt{\frac{\log(1/\delta)}{m}}.$

Theorem 4.2 shows that when time series are independent, learning is not
affected by correlations in the samples and can only be obstructed by the non-
stationarity of the problem, captured via $\Delta$.

Note that when examples are drawn i.i.d., we have $\Delta=0$ in Theorem 4.2:
we recover the standard standard generalization results for regression
problems.

###  4.2 Correlated time series

We now consider several concrete examples of high-dimensional time series in
which different dimensions may be correlated. This setting is common in a
variety of applications including stock market indicators, traffic conditions,
climate observations at different locations, and energy demand.

Suppose that each $Y(i)$ is generated by the auto-regressive (AR) processes
with correlated noise

$y_{t+1}(i)=\Theta_{i}(y^{t}_{0}(i))+\varepsilon_{t+1}(i)$ (4.1)

where the $w_{i}\in\mathbb{R}^{p}$ are unknown parameters and the noise
vectors $\epsilon_{t}\in\mathbb{R}^{m}$ are drawn from a Gaussian distribution
$\mathcal{N}(0,\Sigma)$ where, crucially, $\Sigma$ is not diagonal. The
following lemma is key to our analysis.

######  Lemma 2.

Two AR processes $Y(i),Y(j)$ generated by (4.1) such that
$\sigma=\text{Cov}(Y(i),Y(j))\leq\sigma_{0}<1$ verify
$\beta_{\text{{s2s}}}(i,j)=\max\left(\frac{3}{2(1-\sigma_{0}^{2})},\frac{1}{1-2\sigma_{0}}\right)\sigma=\mathcal{O}(\sigma)$.

#### Hierarchical time series.

As our first example, we consider the case of hierarchical time series that
arises in many real-world applications [39, 37]. Consider the problem of
energy demand forecasting: frequently, one observes a sequence of energy
demands at a variety of levels: single household, local neighborhood, city,
region and country. This imposes a natural hierarchical structure on these
time series.

ptGlobalptFranceptParispt1st ar.pt2nd ar.pt…fptUSAptNew Yorkpt5th Ave.pt8th
Ave.pt…f Figure 2: Hierarchical dependence. Collections for $d=1$ are given by
$C_{1}$ which contains the red leaves and $C_{2}$ which contains the blue
right-side leaves.

Formally, we consider the following hierarchical scenario: a binary tree of
total depth $D$, where time series are generated at each of the leaves. At
each leaf, $Y(i)$ is given by the AR process (4.1) where we impose
$\Sigma_{ij}=(\tfrac{1}{m})^{d(i,j)}$ given $d(i,j)$ the length of the
shortest path from either leaf to the closest common ancestor between $i$ and
$j$. Hence, as $d(i,j)$ increases, $Y(i)$ and $Y(j)$ grow more independent.

For the bound of Theorem 4.1 to be non-trivial, we require a partition
$C_{1},\ldots,C_{k}$ of ${\bf Z}$ such that within a given $C_{j}$ the time
series are close to being independent. One such construction is the following:
fix a depth $d\leq D$ and construct $C_{1},\ldots,C_{2^{d}}$ such that each
$C_{i}$ contains exactly one time series from each sub-tree of depth $D-d$;
hence, $|C_{i}|=2^{D-d}$. Lemma 2 shows that for each $C_{i}$, we have
$\beta(C_{i})=\mathcal{O}(m^{d-D})$. For example, setting
$d=\tfrac{D}{2}=\tfrac{\log m}{2}$, it follows that for any $\delta>0$, with
probability $1-\delta$,

$\displaystyle\mathcal{L}(h|{\bf Y})\leq$
$\displaystyle\widehat{\mathcal{L}}(h)+\max_{j}\left[\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})\right]+\Delta+\frac{1}{\sqrt{2}\sqrt[4]{m}}\sqrt{\log\left(\frac{\sqrt{m}}{\delta-\frac{m}{\mathcal{O}(\sqrt{m^{\log
m}})}}\right)}.$

Furthermore, suppose the model is a linear AR process given by
$y_{t+1}(i)=w_{i}\cdot(y^{t}_{t-p}(i))+\varepsilon_{t+1}(i)$. Then, the
underlying stochastic process is weakly stationary and by Prop. 3 our bound
reduces to: $\mathcal{L}(h|{\bf
Y})\leq\widehat{\mathcal{L}}(h)+\max_{j}\left[\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})\right]+\mathcal{O}\Big{(}\frac{\sqrt{\log
m}}{m^{1/4}}\Big{)}$. By Proposition 5, similar results holds when
$\Theta_{i}$ is periodic.

#### Spatio-temporal processes.

Another common task is spatio-temporal forecasting, in which historical
observation are available at different locations. These observations may
represent temperature at different locations, as in the case of climate
modeling [26, 10], or car traffic at different locations [22].

It is natural to expect correlations between time series to decay as the
geographical distance between them increases. As a simplified example,
consider that the sphere $\mathbb{S}^{3}$ is subdivided according to a
geodesic grid and a time series is drawn from the center of each patch
according to (4.1), also with $\Sigma_{ij}=m^{-d(i,j)}$ but this time with
$d(i,j)$ equal to the (geodesic) distance between the center of two cell
centers. We choose subsets $C_{i}$ with the goal of minimizing the strength of
dependencies between time series within each subsets. Assuming we divide the
sphere into $\sqrt{m}$ collections size approximately $c=\sqrt{m}$ such that
the minimal distance between two points in a set is $d_{0}$, we obtain

$\displaystyle\mathcal{L}(h\mid{\bf Y})\leq$
$\displaystyle\widehat{\mathcal{L}}(h)+\max_{j}\left[\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})\right]+\Delta+\frac{1}{\sqrt{2}\sqrt[4]{m}}\sqrt{\log\left(\frac{\sqrt{m}}{\delta-\mathcal{O}(m^{1-d_{0}})}\right)}.$

As in the case of hierarchical time series, Proposition 3 or Proposition 5 can
be used to remove the dependence on $\Delta$ for certain families of
stochastic processes.

##  5 Comparison to local models

This section provides comparison of learning guarantees for sequence-to-
sequence models with those of local models. In particular, we will compare our
bounds on the generalization gap $\Phi({\bf Y})$ for sequence-to-sequence
models and local models, where the gap is given by

$\displaystyle\Phi_{\text{loc}}({\bf Y})=\sup_{(h_{1},\ldots,h_{m})\in
H^{m}}\Big{[}\mathcal{L}(h_{\text{loc}}\mid{\bf
Y})-\widehat{\mathcal{L}}(h_{\text{loc}})\Big{]}$ (5.1)

where $\widehat{\mathcal{L}}(h_{\text{loc}})$ is the average empirical error
of $h_{i}$ on the sample ${\bf Z}_{i}$, defined as
$\widehat{\mathcal{L}}(h_{\text{loc}})=\frac{1}{mT}\sum\nolimits_{i=1}^{m}\sum\nolimits_{t=1}^{T}f(h_{i},Z_{t,i})$
where $f(h_{i},Z_{t,i})=L(h_{i}(Y_{t-p}^{t-1}(i)),Y_{t}(i))$.

To give a high probability bound for this setting, we take advantage of
existing results for the single local model $h_{i}$ [18]. These results are
given in terms of a slightly different notion of discrepancy $\Delta$, defined
by

$\displaystyle\Delta({\bf Z}_{i})=$
$\displaystyle\sup_{h\in\mathcal{H}}\Bigg{[}\mathbb{E}\Big{[}L(h(Y_{t-p+1}^{T}),Y_{T+1})\mid
Y_{1}^{T}\Big{]}-\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\Big{[}L(h(Y_{t-p}^{t-1}),Y_{t})\mid
Y_{1}^{t-1}\Big{]}\Bigg{]}.$

Another required ingredient to state these results is the expected sequential
covering number $\mathbb{E}_{v\sim
T(\mathbb{P})}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)]$ [18]. For many
hypothesis sets, the log of the sequential covering number admits upper bounds
similar to those presented earlier for the Rademacher complexity. We provide
some examples below and refer the interested reader to [33] for a details.

######  Theorem 5.1.

For $\delta>0$ and $\alpha>0$, with probability at least $1-\delta$, for any
$(h_{1},\ldots,h_{m})$, and any $\alpha>0$,

$\displaystyle\Phi_{\text{loc}}({\bf Y})\leq$
$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\Delta({\bf
Z}_{i})+2\alpha+\sqrt{\frac{2}{T}\log\frac{m\,\max_{i}\mathbb{E}_{v\sim T({\bf
Z}_{i})}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)]}{\delta}}.$

Choosing $\alpha=1/\sqrt{T}$, we can show that, for standard local models such
as the linear hypothesis space $\\{x\to
w^{\top}x,w\in\mathbb{R}^{p},\|w\|_{2}\leq\Lambda\\}$, we have

$\displaystyle\sqrt{\frac{1}{T}\log\frac{2m\,\mathbb{E}_{v\sim
T(Z)}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)]}{\delta}}=\mathcal{O}\Big{(}\sqrt{\frac{\log
m}{T}}\Big{)}.$

In this case, it follows that $\Phi_{\text{loc}}({\bf
Y})\leq\frac{1}{m}\sum_{i=1}^{m}\Delta({\bf
Z}_{i})+\mathcal{O}\Big{(}\sqrt{\frac{\log m}{T}}\Big{)}$. where the last term
in this bound should be compared with corresponding (non-discrepancy) terms in
the bound of Theorem. 4.1, which, as discussed above, scales as
$\mathcal{O}(\sqrt{T/m})$ for a variety of different hypothesis sets.

Hence, when we have access to relatively few time series compared to their
length ($m\ll T$), learning to predict each time series as its own independent
problem will with high probability lead to a better generalization bound. On
the other hand, in extremely high-dimensional settings when we have
significantly more time series than time steps ($m\gg T$), sequence-to-
sequence learning will (with high probability) provide superior performance.
We also expect the performance of sequence-to-sequence models to deteriorate
as the correlation between time series increases.

A direct comparison of bounds in Theorem 4.1 and Theorem 5.1 is complicated by
the fact that discrepancies that appear in these results are different. In
fact, it is possible to design examples where
$\tfrac{1}{m}\sum_{i=1}^{m}\Delta({\bf Z}_{i})$ is constant and $\Delta$ is
negligible, and vice-versa.

Consider a tent function $g_{b}$ such that $g_{b}(s)=2bs/T$ for $s\in[0,T/2]$
and $g_{b}(s)=-2bs/T+2b$ for $s\in[T/2,T]$. Let $f_{b}$ be its periodic
extension to the real line, and define $\mathcal{S}=\\{f_{b}\colon
b\in[0,1]\\}$. Suppose that we sample uniformly $b\in[0,1]$ and
$s\in\\{0,T/2\\}$ $m$ times, and observe time series
$f_{b_{i}}(s_{i}),\ldots,f_{b_{i}}(s_{i}+T)$. Then, as we have shown in
Proposition 5, $\Delta$ is negligible for sequence-to-sequence models.
However, unless the model class is trivial, it can be shown that $\Delta({\bf
Z}_{i})$ is bounded away from zero for all $i$.

Conversely, suppose we sample uniformly $b\in[0,1]$ $m$ times and observe time
series $f_{b_{i}}(0),\ldots,f_{b_{i}}(T/2+1)$. Consider a set of local models
that learn an offset from the previous point $\\{h\colon x\mapsto
x+c,c\in[0,1]\\}$. It can be shown that in this case $\Delta({\bf Z}_{i})=0$,
whereas $\Delta$ is bounded away from zero for any non-trivial class of
sequence-to-sequence models.

From a practical perspective, we can simply use $\Delta_{s}$ and empirical
estimates of $\Delta({\bf Z}_{i})$ to decide whether to choose sequence-to-
sequence or local models.

We conclude this section with an observation that similar results to Theorem
4.1 can be proved for multivariate local models with the only difference that
the sample complexity of the problem scales as $O(\sqrt{m/T})$, and hence
these models are even more prone to the curse of dimensionality.

##  6 Hybrid models

In this section, we discuss models that interpolate between local and
sequence-to-sequence models. This hybrid approach trains a single model $h$ on
the union of _local_ training sets ${\bf Z}_{1},\ldots,{\bf Z}_{m}$ used to
train $m$ models in the local approach. The bounds that we state here require
the following extension of the discrepancy to $\Delta_{t}$, defined as

$\displaystyle\Delta_{t}=$
$\displaystyle\frac{1}{m}\sup_{h\in\mathcal{H}}\Big{|}\sum_{i=1}^{m}\mathbb{E}_{\mathcal{D}}[L(h(Y_{t-p-1}^{t-1}(i)),Y_{t}(i))|Y_{1}^{t-1}]-\mathbb{E}_{\mathcal{D}^{\prime}}[L(h(Y_{T-p}^{T}(i)),Y_{T+1}(i))|{\bf
Y}]\Big{|}$

Many of the properties that were discussed for the discrepancy $\Delta$ carry
over to $\Delta_{t}$ as well. The empirical error in this case is the same as
for the local models:

$\widehat{\mathcal{L}}(h)=\frac{1}{mT}\sum_{i=1}^{m}\sum_{t=1}^{T}f(h,Z_{t,i}).$

Observe that one straightforward way to obtain a bound for hybrid models is to
apply Theorem 5.1 with $(h,\ldots,h)\in\mathcal{H}^{m}$. Alternatively, we can
apply Theorem 4.1 at every time point $t=1,\ldots,T$.

Combining these results via union bound leads to the following learning
guarantee for hybrid models.

######  Theorem 6.1.

Let $C_{1},\ldots,C_{k}$ form a partition of the training input ${\bf Z}$ and
let $I_{j}$ denote the set of indices of time series that belong to $C_{j}$.
Assume that the loss function $L$ is bounded by 1. Then, for any $\delta>0$,
with probability $1-\delta$, for any $h\in\mathcal{H}$ and any $\alpha>0$

${\mathcal{L}(h\mid{\bf Y})\leq\widehat{\mathcal{L}}(h)+\min(B_{1},B_{2})},$

where

$\displaystyle B_{1}=$
$\displaystyle\frac{1}{T}\sum\nolimits_{t=1}^{T}\Delta_{t}+\max_{j}\widehat{}\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})+\frac{1}{\sqrt{2\min_{j}|I_{j}|}}\sqrt{\log\left(\frac{2Tk}{\delta-2\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})}\right)}$
$\displaystyle B_{2}=$
$\displaystyle\frac{1}{m}\sum\nolimits_{i=1}^{m}\Delta({\bf
Z}_{i})+2\alpha+\sqrt{\frac{2}{T}\log\frac{2m\,\max_{i}\mathbb{E}_{v\sim
T({\bf Z}_{i})}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)]}{\delta}}.$

Using the same arguments for the complexity terms as in the case of sequence-
to-sequence and local models, this result shows that hybrid models are
successful with high probability when $m\gg T$ or correlations between time
series are strong, as well as when $T\gg m$.

Potential costs for this model are hidden in the new discrepancy term
$\frac{1}{T}\sum_{t=1}^{T}\Delta_{t}$. This term leads to different bounds
depending on the particular non-stationarity in the given problem. As before
this trade-off can be accessed empirically using the data-dependent version of
discrepancy.

Note that the above bound does not imply that hybrid models are superior to
local models: using $m$ hypotheses $h_{1},\ldots,h_{m}$ can help us achieve a
better trade-off between $\widehat{\mathcal{L}}(h)$ and $B_{2}$, and vice
versa.

##  7 Conclusion

We formally introduce sequence-to-sequence learning for time series, a
framework in which a model learns to map past sequences of length $T$ to their
next values. We provide the first generalization bounds for sequence-to-
sequence modeling. Our results are stated in terms of new notions of
discrepancy and expected mixing coefficients. We study these new notions for
several different families of stochastic processes including stationary,
weakly stationary, periodic, hierarchical and spatio-temporal time series.

Furthermore, we show that our discrepancy can be computed from data, making it
a useful tool for practitioners to empirically assess the non-stationarity of
their problem. In particular, the discrepancy can be used to determine whether
the sequence-to-sequence methodology is likely to succeed based on the
inherent non-stationarity of the problem.

Furthermore, compared to the local framework for time series forecasting, in
which independent models for each one-dimensional time series are learned, our
analysis shows that the sample complexity of sequence-to-sequence models
scales as $\mathcal{O}(\sqrt{T/m})$, providing superior guarantees when the
number $m$ of time series is significantly greater than the length $T$ of each
series, provided that different series are weakly correlated.

Conversely, we show that the sample complexity of local models scales as
$\mathcal{O}(\sqrt{\log(m)/T})$, and should be preferred when $m\ll T$ or when
time series are strongly correlated. We also study hybrid models for which
learning guarantees are favorable both when $m\gg T$ and $T\gg m$, but which
have a more complex trade-off in terms of discrepancy.

As a final note, the analysis we have carried through is easily extended to
show similar results for the sequence-to-sequence scenario when the test data
includes _new_ series not observed during training, as is often the case in a
variety of applications.

## References

  * Banbura et al. [2010] Marta Banbura, Domenico Giannone, and Lucrezia Reichlin.  Large Bayesian vector auto regressions.  _Journal of Applied Econometrics_ , 25(1):71–92, 2010.
  * Basu and Michailidis [2015] Sumanta Basu and George Michailidis.  Regularized estimation in sparse high-dimensional time series models.  _Ann. Statist._ , 43(4):1535–1567, 2015.
  * Bianchi et al. [2017] Filippo Maria Bianchi, Enrico Maiorino, Michael C. Kampffmeyer, Antonello Rizzi, and Robert Jenssen.  _Recurrent Neural Networks for Short-Term Load Forecasting - An Overview and Comparative Analysis_.  Springer Briefs in Computer Science. Springer, 2017.
  * Bollerslev [1986] Tim Bollerslev.  Generalized autoregressive conditional heteroskedasticity.  _Journal of Econometrics_ , 31(3):307 – 327, 1986\.
  * Box and Jenkins [1990] George Edward Pelham Box and Gwilym Jenkins.  _Time Series Analysis, Forecasting and Control_.  Holden-Day, Incorporated, 1990.
  * Brockwell and Davis [1986] Peter J Brockwell and Richard A Davis.  _Time Series: Theory and Methods_.  Springer-Verlag New York, Inc., 1986.
  * Doukhan [1994] P. Doukhan.  _Mixing: Properties and Examples_.  Lecture notes in statistics. Springer, 1994.
  * Engle [1982] Robert Engle.  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  _Econometrica_ , 50(4):987–1007, 1982.
  * Flunkert et al. [2017] Valentin Flunkert, David Salinas, and Jan Gasthaus.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _Arxiv:1704.04110_ , 2017.
  * Ghafarianzadeh and Monteleoni [2013] Mahsa Ghafarianzadeh and Claire Monteleoni.  Climate prediction via matrix completion.  In _Late-Breaking Developments in the Field of Artificial Intelligence_ , volume WS-13-17 of _AAAI Workshops_. AAAI, 2013.
  * Goel et al. [2017] Hardik Goel, Igor Melnyk, and Arindam Banerjee.  R2N2: Residual recurrent neural networks for multivariate time series forecasting.  _arXiv:1709.03159_ , 2017.
  * Hamilton [1994] James Douglas Hamilton.  _Time series analysis_.  Princeton Univ. Press, 1994.
  * Han et al. [2015a] Fang Han, Huanran Lu, and Han Liu.  A direct estimation of high dimensional stationary vector autoregressions.  _Journal of Machine Learning Research_ , 16:3115–3150, 2015a.
  * Han et al. [2015b] Fang Han, Sheng Xu, and Han Liu.  Rate optimal estimation of high dimensional time series.  Technical report, Technical Report, Johns Hopkins University, 2015b.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural Comput._ , 9(8):1735–1780, 1997.  ISSN 0899-7667.
  * Koltchinskii and Panchenko [2002] V. Koltchinskii and D. Panchenko.  Empirical margin distributions and bounding the generalization error of combined classifiers.  _Ann. Statist._ , 30(1):1–50, 2002.
  * Kuznetsov and Mohri [2014] Vitaly Kuznetsov and Mehryar Mohri.  _Generalization Bounds for Time Series Prediction with Non-stationary Processes_ , pages 260–274.  Springer International Publishing, 2014.
  * Kuznetsov and Mohri [2015] Vitaly Kuznetsov and Mehryar Mohri.  Learning theory and algorithms for forecasting non-stationary time series.  In _Advances in Neural Information Processing Systems 28_ , pages 541–549. Curran Associates, Inc., 2015.
  * Kuznetsov and Mohri [2016] Vitaly Kuznetsov and Mehryar Mohri.  Time series prediction and online learning.  In _29th Annual Conference on Learning Theory_ , volume 49 of _Proceedings of Machine Learning Research_ , pages 1190–1213. PMLR, 2016\.
  * Kuznetsov and Mohri [2017] Vitaly Kuznetsov and Mehryar Mohri.  Discriminative state space models.  In _NIPS, Long Beach, CA, USA_ , 2017.
  * Laptev et al. [2017] Nikolay Laptev, Jason Yosinski, Li Erran Li, and Slawek Smyl.  Time-series extreme event forecasting with neural networks at uber.  In _ICML Workshop_ , 2017.
  * Li et al. [2017] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _arXiv:1707.01926_ , 2017.
  * Lütkepohl [2006] Helmut Lütkepohl.  Chapter 6 Forecasting with VARMA models.  In _Handbook of Economic Forecasting_ , volume 1, pages 287 – 325\. Elsevier, 2006.
  * Lütkepohl [2007] Helmut Lütkepohl.  _New Introduction to Multiple Time Series Analysis_.  Springer Publishing Company, Incorporated, 2007.
  * Lv et al. [2015] Yisheng Lv, Yanjie Duan, Wenwen Kang, Zhengxi Li, and Fei-Yue Wang.  Traffic flow prediction with big data: A deep learning approach.  _IEEE Transactions on Intelligent Transportation Systems_ , 16:865–873, 2015.
  * McQuade and Monteleoni [2012] Scott McQuade and Claire Monteleoni.  Global climate model tracking using geospatial neighborhoods.  In _Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence, July 22-26, 2012, Toronto, Ontario, Canada._ AAAI Press, 2012.
  * Meir and Hellerstein [2000] Ron Meir and Lisa Hellerstein.  Nonparametric time series prediction through adaptive model selection.  In _Machine Learning_ , pages 5–34, 2000.
  * Mohri and Rostamizadeh [2009] Mehryar Mohri and Afshin Rostamizadeh.  Rademacher complexity bounds for non-i.i.d. processes.  In _Advances in Neural Information Processing Systems 21_ , pages 1097–1104. Curran Associates, Inc., 2009.
  * Mohri and Rostamizadeh [2010] Mehryar Mohri and Afshin Rostamizadeh.  Stability bounds for stationary $\phi$-mixing and $\beta$-mixing processes.  _J. Mach. Learn. Res._ , 11:789–814, 2010.  ISSN 1532-4435.
  * Mohri et al. [2012] Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  _Foundations of Machine Learning_.  The MIT Press, 2012.
  * Negahban and Wainwright [2011] Sahand Negahban and Martin J. Wainwright.  Estimation of (near) low-rank matrices with noise and high-dimensional scaling.  _Ann. Statist._ , 39(2):1069–1097, 2011.
  * Neyshabur et al. [2015] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  Norm-based capacity control in neural networks.  In _Proceedings of The 28th Conference on Learning Theory_ , volume 40 of _Proceedings of Machine Learning Research_ , pages 1376–1401. PMLR, 2015.
  * Rakhlin et al. [2015] Alexander Rakhlin, Karthik Sridharan, and Ambuj Tewari.  Sequential complexities and uniform martingale laws of large numbers.  _Probability Theory and Related Fields_ , 161(1-2):111–153, 2015.
  * Romeu et al. [2013] Pablo Romeu, Francisco Zamora-Martínez, Paloma Botella-Rocamora, and Juan Pardo.  _Time-Series Forecasting of Indoor Temperature Using Pre-trained Deep Neural Networks_ , pages 451–458.  Springer Berlin Heidelberg, 2013.
  * Song and Bickel [2011] Song Song and Peter J. Bickel.  Large vector auto regressions.  _arXiv:1106.3915_ , 2011.
  * Sutskever et al. [2014] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In _Advances in Neural Information Processing Systems 27_ , pages 3104–3112. Curran Associates, Inc., 2014.
  * Taieb et al. [2017] Souhaib Ben Taieb, James W. Taylor, and Rob J. Hyndman.  Coherent probabilistic forecasts for hierarchical time series.  In _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 3348–3357. PMLR, 2017.
  * W. Sun [2015] D. Malioutov W. Sun.  Time series forecasting with shared seasonality patterns using non-negative matrix factorization.  In _The 29th Annual Conference on Neural Information Processing Systems (NIPS)_. Time Series Workshop, 2015.
  * Wickramasuriya et al. [2015] Shanika L Wickramasuriya, George Athanasopoulos, and Rob J Hyndman.  Forecasting hierarchical and grouped time series through trace minimization.  Technical Report 15/15, Monash University, Department of Econometrics and Business Statistics, 2015.
  * Yu [1994] Bin Yu.  Rates of convergence for empirical processes of stationary mixing sequences.  _Ann. Probab._ , 22(1):94–116, 1994.
  * Yu et al. [2016] Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon.  Temporal regularized matrix factorization for high-dimensional time series prediction.  In _Advances in Neural Information Processing Systems 29_ , pages 847–855. Curran Associates, Inc., 2016.
  * Yu et al. [2017] Rose Yu, Stephan Zheng, Anima Anandkumar, and Yisong Yue.  Long-term forecasting using tensor-train RNNs.  _Arxiv:1711.00073_ , 2017.
  * Zhu and Laptev [2017] Lingxue Zhu and Nikolay Laptev.  Deep and confident prediction for time series at Uber.  _arXiv:1709.01907_ , 2017.
  * Zimin and Lampert [2017] Alexander Zimin and Christoph H. Lampert.  Learning theory for conditional risk minimization.  In _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, AISTATS 2017_ , pages 213–222, 2017\.

##  Appendix A Rademacher complexity

######  Definition 6 (Rademacher complexity).

Given a family of functions $\mathcal{F}$ and a training set ${\bf
Z}=\\{Z_{1},\ldots,Z_{m}\\}$, _the Rademacher complexity of $\mathcal{F}$
conditioned on ${\bf Y}^{\prime}$_ is given by

$\widehat{\mathfrak{R}}_{\bf Z}(\mathcal{F})=\mathbb{E}_{{\bf
Z},\sigma}\left[\max_{f\in\mathcal{F}}\frac{1}{m}\sum\nolimits_{i=1}^{m}\sigma_{i}f(Z_{i})\Big{|}{\bf
Y}^{\prime}\right]$

where $\sigma_{1},\ldots,\sigma_{m}$ are i.i.d. random variables uniform on
$\\{-1,+1\\}$. _The Rademacher complexity_ of $\mathcal{F}$ for sample size
$m$ is given by

$\mathfrak{R}_{m}(\mathcal{F})=\mathbb{E}_{{\bf
Y}^{\prime}}\left[\widehat{\mathfrak{R}}_{\bf Z}(\mathcal{F})\right].$

The Rademacher complexity has been studied for a variety of function classes.
For instance, for the linear hypothesis space $\mathcal{H}=\\{x\rightarrow
w^{\top}x,\|w\|_{2}\leq\Lambda\\}$, $\widehat{}\mathfrak{R}_{\bf Z}$ can be
upper bounded by $\widehat{}\mathfrak{R}_{\bf
Z}(\mathcal{H})\leq\frac{\Lambda}{\sqrt{m}}\max_{i}\|Z_{i}\|_{2}$. As another
example, the hypothesis class of ReLu feed-forward neural networks with $d$
layers and weight matrices $W_{k}$ such that
$\prod_{k=1}^{d}\|W\|_{F}\leq\gamma$ verifies $\widehat{}\mathfrak{R}_{\bf
Z}(\mathcal{H})\leq\frac{2^{d-1/2}\gamma}{\sqrt{m}}\max_{i}\|Z_{i}\|_{2}$
[32].

##  Appendix B Discrepancy analysis

######  Proposition 1.

Let $\mathcal{H}$ be a hypothesis space and let $L$ be a bounded loss function
which respects the triangle inequality. Let $h^{\prime}\in\mathcal{H}$. Then,

$\Delta\leq\Delta_{s}+\mathcal{L}(h\mid{\bf Y})+\mathcal{L}(h\mid{\bf
Y}^{\prime})$

###### Proof.

Let $h,h^{\prime}\in\mathcal{H}$. For ease of notation, we write

$\displaystyle\Delta_{s}(h,h^{\prime},{\bf Y}^{\prime})=$
$\displaystyle\frac{1}{m}\sum_{i}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-\frac{1}{m}\sum_{i}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i))).$

Applying the triangle inequality to $L$,

$\displaystyle\mathcal{L}(h\mid{\bf Y})=$
$\displaystyle\frac{1}{m}\sum_{i}\mathbb{E}[L(h(Y_{1}^{T}(i)),Y_{T+1}(i))\mid{\bf
Y}]$ $\displaystyle\leq$
$\displaystyle\frac{1}{m}\sum_{i}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))+\frac{1}{m}\sum_{i}\mathbb{E}[L(h^{\prime}(Y_{1}^{T}(i)),Y_{T+1}(i))\mid{\bf
Y}]$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{i}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))+\mathcal{L}(h^{\prime}\mid{\bf
Y}).$

Then, by definition of $\Delta_{s}(h,h^{\prime},{\bf Y}^{\prime})$, we have

$\displaystyle\mathcal{L}(h\mid{\bf Y})\leq$
$\displaystyle\frac{1}{m}\sum_{i}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-\frac{1}{m}\sum_{i}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))$
$\displaystyle+\frac{1}{m}\sum_{i}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))+\mathcal{L}(h^{\prime}\mid{\bf
Y})$ $\displaystyle\leq$ $\displaystyle\Delta_{s}(h,h^{\prime},{\bf
Y}^{\prime})+\mathcal{L}(h^{\prime}\mid{\bf
Y})+\frac{1}{m}\sum_{i}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i))).$

By an application of the triangle inequality to $L$,

$\displaystyle\mathcal{L}(h,\mathcal{D})\leq$
$\displaystyle\Delta_{s}(h,h^{\prime},{\bf
Y}^{\prime})+\mathcal{L}(h^{\prime}\mid{\bf
Y})+\frac{1}{m}\sum_{i}\mathbb{E}[L(h(Y_{1}^{T-1}(i)),Y_{T}(i))\mid{\bf
Y}^{\prime}]$
$\displaystyle+\frac{1}{m}\sum_{i}\mathbb{E}[L(h^{\prime}(Y_{1}^{T-1}(i)),Y_{T}(i))\mid{\bf
Y}^{\prime}]$ $\displaystyle=$ $\displaystyle\Delta_{s}(h,h^{\prime},{\bf
Y}^{\prime})+\mathcal{L}(h^{\prime}\mid{\bf Y})+\mathcal{L}(h\mid{\bf
Y}^{\prime})+\mathcal{L}(h^{\prime}\mid{\bf Y}^{\prime}).$

Finally, we obtain

$\displaystyle\mathcal{L}(h\mid{\bf Y})-\mathcal{L}(h\mid{\bf
Y}^{\prime})\leq$ $\displaystyle\Delta_{s}(h,h^{\prime},{\bf
Y}^{\prime})+\mathcal{L}(h^{\prime}\mid{\bf Y})+\mathcal{L}(h^{\prime}\mid{\bf
Y}^{\prime})$

and the result announced in the theorem follows by taking the supremum over
$\mathcal{H}$ on both sides. ∎

######  Proposition 2.

Let $I_{1},\cdots,I_{k}$ be a partition of $\\{1,\ldots,m\\}$, and
$C_{1},\ldots,C_{k}$ be the corresponding partition of ${\bf Y}$. Write
$c=\min_{j}|C_{j}|$. Then we have with probability $1-\delta$,

$\displaystyle\Delta_{s}\leq$
$\displaystyle\Delta_{e}+{\max\Big{(}\max_{j}\mathfrak{R}_{|C_{j}|}(\widetilde{C}^{\prime}_{j}),\max_{j}\mathfrak{R}_{|C_{j}|}(\widetilde{I}_{j})\Big{)}}+\sqrt{\frac{1}{2c}\log\frac{2k}{\delta-\sum_{j}(|I_{j}|-1)[\bar{\beta}(I_{j})+\bar{\beta}^{\prime}(I_{j})]}}.$

###### Proof.

By definition of $\Delta_{s}$,

$\displaystyle\Delta_{s}=$
$\displaystyle\sup_{h,h^{\prime}\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\Big{[}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))\Big{]}$
$\displaystyle\leq$
$\displaystyle\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\frac{1}{m}\sum_{i=1}^{m}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-\mathbb{E}_{Y}[L(h(Y_{1}^{T}),h^{\prime}(Y_{1}^{T}))]\Big{]}$
$\displaystyle+\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\mathbb{E}_{Y}[L(h(Y_{1}^{T}),h^{\prime}(Y_{1}^{T}))]-\mathbb{E}_{Y}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))]\Big{]}$
$\displaystyle+\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\mathbb{E}_{Y}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))]-\frac{1}{m}\sum_{i=1}^{m}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))\Big{]}$

by sub-additivity of the supremum. Now, define

$\displaystyle\phi({\bf Y})\triangleq$
$\displaystyle\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\frac{1}{m}\sum_{i=1}^{m}L(h(Y_{1}^{T}(i)),h^{\prime}(Y_{1}^{T}(i)))-\mathbb{E}_{Y}[L(h(Y_{1}^{T}),h^{\prime}(Y_{1}^{T}))]\Big{]}$
$\displaystyle\psi({\bf Y}^{\prime})\triangleq$
$\displaystyle\sup_{h,h^{\prime}\in\mathcal{H}}\Big{[}\mathbb{E}_{Y}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))-\frac{1}{m}\sum_{i=1}^{m}L(h(Y_{1}^{T-1}(i)),h^{\prime}(Y_{1}^{T-1}(i)))]\Big{]}.$

By definition of $\Delta_{e}$, we have from the previous inequality

$\Delta_{s}\leq\Delta_{e}+\phi({\bf Y}_{1}^{T})+\psi({\bf Y}_{1}^{T-1}).$

We now proceed to give a high-probability bound for $\phi$; the same reasoning
will yield a bound for $\psi$. By sub-additivity of the max,

$\displaystyle\phi({\bf Y})\leq$
$\displaystyle\sum_{j}\frac{|C_{j}|}{m}\sup_{h\in\mathcal{H}}\Big{[}\mathbb{E}_{Y}[f(h,Y_{1}^{T})]-\frac{1}{|C_{j}|}\sum_{Y\in
C_{j}}f(h,Y_{1}^{T})\Big{]}$ $\displaystyle\leq$
$\displaystyle\sum_{j}\frac{|C_{j}|}{m}\phi(C_{j})$

and so by union bound, for $\epsilon>0$

$\Pr(\phi({\bf Y})>\epsilon)\leq\sum_{j}\Pr(\phi(C_{j})>\epsilon).$

Let $\epsilon>\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})]$ and set
$\epsilon_{j}=\epsilon-\mathbb{E}[\phi(\widetilde{C}_{j})]$.

Define for time series $Y(i),Y(j)$ the mixing coefficient

$\bar{\beta}(i,j)=\|\Pr(Y_{1}^{T}(i),Y_{1}^{T}(j))-\Pr(Y_{1}^{T}(i))\Pr(Y_{1}^{T}(j))\|_{TV}$

where we also extend the usual notation to $\bar{\beta}(C_{j})$.

$\displaystyle\Pr\left(\phi(C_{j})>\epsilon\right)=$
$\displaystyle\Pr\Big{(}\phi(C_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})]>\epsilon_{j}\Big{)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}$
$\displaystyle\Pr\Big{(}\phi(\widetilde{C}_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})]>\epsilon_{j}\Big{)}+{(|I_{j}|-1)\bar{\beta}(I_{j})}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}$ $\displaystyle
e^{-2c\epsilon_{j}^{2}}+(|I_{j}|-1)\bar{\beta}(I_{j}),$

where (a) follows by applying Prop. 6 to the indicator function of the event
$\Pr(\phi(C_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})]\geq\epsilon)$, and (b) is
a direct application of McDiarmid’s inequality to
$\phi(\widetilde{C}_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})]$.

Hence, by summing over $j$ we obtain

$\displaystyle\Pr\left(\phi({\bf Y})>\epsilon\right)\leq$ $\displaystyle
ke^{-2\min_{j}|C_{j}|(\epsilon-\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})])^{2}}+\sum_{j}(|I_{j}|-1)\bar{\beta}(I_{j})$

and similarly

$\displaystyle\Pr\left(\psi({\bf Y}^{\prime})>\epsilon\right)\leq$
$\displaystyle
ke^{-2\min_{j}|C^{\prime}_{j}|(\epsilon-\max_{j}\mathbb{E}[\psi(\widetilde{C}^{\prime}_{j})])^{2}}+\sum_{j}(|I_{j}|-1)\bar{\beta}^{\prime}(I_{j}),$

which finally yields

$\displaystyle\Pr$
$\displaystyle(\Delta_{s}-\Delta_{e}>\epsilon)\leq\Pr(\phi({\bf
Y})>\epsilon)+\Pr(\psi({\bf Y}^{\prime})>\epsilon)$ $\displaystyle\leq
2k\exp(-2c(\epsilon-\max(\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})],\max_{j}\mathbb{E}[\psi(\widetilde{C}^{\prime}_{j})]))^{2})+\sum_{j}(|I_{j}|-1)[\bar{\beta}(I_{j})+\bar{\beta}^{\prime}(I^{\prime}_{j})],$

where we recall that we write $c=\min_{j}|C_{j}|$. We invert the previous
equation by setting

$\displaystyle\epsilon=$
$\displaystyle\max(\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})],\max_{j}\mathbb{E}[\psi(\widetilde{C}^{\prime}_{j})])+\sqrt{\frac{1}{2c}\log\frac{2k}{\delta-\sum_{j}(|I_{j}|-1)[\beta(I_{j})+\beta(I^{\prime}_{j})]}},$

yielding with probability $1-\delta$,

$\displaystyle\Delta_{s}\leq$
$\displaystyle\Delta_{e}+\max(\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})],\max_{j}\mathbb{E}[\psi(\widetilde{C}^{\prime}_{j})])+\sqrt{\frac{1}{2c}\log\frac{2k}{\delta-\sum_{j}(|I_{j}|-1)[\bar{\beta}(I_{j})+\bar{\beta}^{\prime}(I_{j})]}}.$

We now bound $\mathbb{E}[\phi(\widetilde{C}_{j})]$ by
$\mathfrak{R}_{|C_{j}|}(\widetilde{C}_{j})$. A similar argument yields the
bound for $\psi$. By definition, we have

$\displaystyle\mathbb{E}[\phi(\widetilde{C}_{j})]$
$\displaystyle=\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\frac{1}{|C_{j}|}\sum_{Z\in\widetilde{C}_{j}}f(h,Y_{1}^{T}(i))-\mathbb{E}_{Y}[f(h,Y_{1}^{T})]\Big{]}$
$\displaystyle=\frac{1}{|C_{j}|}\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\sum_{Z\in\widetilde{C}_{j}}\underbrace{f(h,Y_{1}^{T}(i))\\!-\\!\mathbb{E}_{Y}[f(h,Y_{1}^{T})}_{g(h,Y_{1}^{T}(i))}]\Big{]}$
$\displaystyle=\frac{1}{|C_{j}|}\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\sum_{Z\in\widetilde{C}_{j}}g(h,Y_{1}^{T}(i))\Big{]}$

Standard symmetrization arguments as those used for the proof of the famous
result by [16], which hold also when data is drawn independently but not
identically at random, yield

$\mathbb{E}[\phi(\widetilde{C}_{j})]\leq\mathfrak{R}_{|C_{j}|}(\widetilde{C}_{j}).$

The same argument yields for $\psi$

$\mathbb{E}[\psi(\widetilde{C}^{\prime}_{j})]\leq\mathfrak{R}_{|C_{j}|}(\widetilde{C}^{\prime}_{j}).$

To conclude our proof, it only remains to prove the bound

$\displaystyle\bar{\beta}(i,j)\leq$
$\displaystyle\beta_{\text{{s2s}}}(i,j)+\mathbb{E}_{{\bf
Y}^{\prime}}\Big{[}\text{Cov}\Big{(}\Pr(Y_{T}(i)\mid{\bf
Y}^{\prime}),\Pr(Y_{T}(j)\mid{\bf Y}^{\prime})\Big{)}\Big{]}$

Let $Y(i),Y(j)$ be two time series, and write
$X_{i}=\mathbb{E}[\Pr(Y_{1}^{T}(i))\mid{\bf Y}^{\prime}]$. Then the following
bound holds

$\displaystyle\bar{\beta}(i,j)=$
$\displaystyle\|\Pr(Y_{1}^{T}(i),Y_{1}^{T}(j))-\Pr(Y_{1}^{T}(i))\Pr(Y_{1}^{T}(j))\|_{TV}$
$\displaystyle=$
$\displaystyle\|\mathbb{E}[\Pr(Y_{1}^{T}(i),Y_{1}^{T}(j))\mid{\bf
Y}^{\prime}]-\mathbb{E}[X_{i}]\mathbb{E}[X_{j}]\|_{TV}$ $\displaystyle=$
$\displaystyle\|\mathbb{E}[\Pr(Y_{1}^{T}(i),Y_{1}^{T}(j))\mid{\bf
Y}_{1}^{T-1}]-\mathbb{E}[X_{i},X_{j}]-\mathbb{E}[\text{Cov}(X_{i},X_{j})]\|_{TV}$
$\displaystyle\leq$ $\displaystyle\beta_{\text{{s2s}}}(i,j)+\mathbb{E}_{{\bf
Y}^{\prime}}[\text{Cov}(X_{i},X_{j})],$

which is the desired inequality. ∎

We now show two useful lemmas for various specific cases of time series and
hypothesis spaces.

######  Proposition 3.

If $Y(i)$ is stationary for all $1\leq i\leq m$, and $\mathcal{H}$ is a
hypothesis space such that $h\in\mathcal{H}:\mathcal{Y}^{T-1}\to\mathcal{Y}$
(i.e. the hypotheses only consider the last $T-1$ values of $Y$), then
$\Delta_{e}=0$.

###### Proof.

Let $h,h^{\prime}\in\mathcal{H}$. For stationary $Y(i)$, we have
$\Pr(Y_{1}^{T}(i))=\Pr(Y_{2}^{T}(i))$, and so

$\mathbb{E}[L(h(Y_{2}^{T}),h^{\prime}(Y_{2}^{T}))]-\mathbb{E}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))]=0$

and so taking the supremum over $h,h^{\prime}$ yields the desired result. ∎

######  Proposition 4.

If $Y(i)$ is covariance stationary for all $1\leq i\leq m$, $L$ is the squared
loss, and $\mathcal{H}$ is a linear hypothesis space $\\{x\to w\cdot
x\mid\|w\|\in\mathbb{R}^{p}\leq\Lambda\\}$, then $\Delta_{e}=0$.

###### Proof.

Recall that a time series $Y$ is covariance stationary if
$\mathbb{E}_{Y}[Y_{t}]$ does not depend on $t$ and
$\mathbb{E}_{Y}[Y_{t}Y_{s}]=f(t-s)$ for some function $f$.

Let now $(h,h^{\prime})\in\mathcal{H}\equiv(w,w^{\prime})\in\mathbb{R}^{p}$.
We write $\Sigma=\Sigma_{2}^{T}(Y)=\Sigma_{1}^{T}(Y)$ the covariance matrix of
$Y$ where the equality follows from covariance stationarity. Without loss of
generality, we consider $p=T-1$. Then,

$\displaystyle\mathbb{E}$
$\displaystyle[L(h(Y_{2}^{T}),h^{\prime}(Y_{2}^{T}))]-\mathbb{E}[L(h(Y_{1}^{T-1}),h^{\prime}(Y_{1}^{T-1}))]$
$\displaystyle=\mathbb{E}[((w-w^{\prime})^{\top}\Sigma_{2}^{T}(Y)(w-w^{\prime})]-\mathbb{E}[((w-w^{\prime})^{\top}\Sigma_{1}^{T-1}(Y)(w-w^{\prime})]$
$\displaystyle=0.$

Taking the supremum over $h,h^{\prime}$ yields the desired result. ∎

######  Proposition 5.

If the $Y(i)$ are periodic of period $p$ and the observed starting time of
each $Y(i)$ is distributed uniformly at random in $[p]$, then $\Delta_{e}=0$.

###### Proof.

This proof is similar to the stationary case: indeed, we can write
$\Pr(Y_{1}^{T-1}(i))=\frac{1}{p}\Pr(Y(i))$ due to the uniform distribution on
starting times. Then, by the same reasoning, we have also

$\Pr(Y_{2}^{T}(i))=\frac{1}{p}\Pr(Y(i))=\Pr(Y_{1}^{T-1}(i)),$

from which the result follows. ∎

##  Appendix C Generalization bounds

######  Proposition 6.

_Yu [40, Corollary 2.7]_. Let $f$ be a real-valued Borel measurable function
such that $0\leq f\leq 1$. Then, we have the following guarantee:

$\left|\mathbb{E}[f(\widetilde{C})]-\mathbb{E}[f(C)]\right|\leq(|C|-1)\beta,$

where $\beta$ is the total variation distance between joint distributions of
$C$ and $\widetilde{C}$.

######  Theorem 4.1.

Let $\mathcal{H}$ be a hypothesis space, and $h\in\mathcal{H}$. Let
$C_{1},\ldots,C_{k}$ form a partition of the training input ${\bf Y}_{1}^{T}$,
and consider that the loss function $L$ is bounded by 1. Then, we have for
$\delta>0$, with probability $1-\delta$,

$\displaystyle\Phi_{\text{s2s}}(h)\leq$
$\displaystyle\Delta+\max_{j}\left[\mathfrak{R}_{|C_{j}|}(\widetilde{C}_{j}\mid{\bf
Y})\right]+\frac{1}{\sqrt{2\min_{j}|I_{j}|}}\sqrt{\log\left(\frac{k}{\delta-\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})}\right)}.$

For ease of notation, we write

$\displaystyle\phi({\bf Y})=$
$\displaystyle\sup_{h\in\mathcal{H}}\mathcal{L}(h\mid{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h,{\bf Y})$ $\displaystyle=$
$\displaystyle\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}[f(h,Y_{1}^{T}(i))\mid{\bf
Y}^{\prime}]-\frac{1}{m}\sum_{i=1}^{m}f(h,Y_{1}^{T}(i)).$

We begin by proving the following lemma.

######  Lemma 3.

Let $\bar{}{\bf Y}$ be equal to ${\bf Y}$ on all time series except for the
last, where we have $\bar{Y}(m)=Y(m)$ at all times except for time $t=T$. Then

$\left|\phi({\bf Y})-\phi(\bar{}{\bf Y})\right|\leq\frac{1}{m}$

###### Proof.

Fix $h^{*}\in\mathcal{H}$. Then,

$\displaystyle\mathcal{L}$ $\displaystyle(h^{*}\mid{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h^{*},{\bf
Y})-\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}(h\mid\bar{}{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h,\bar{{\bf Y}})\Big{]}$
$\displaystyle\leq\mathcal{L}(h^{*}\mid{\bf
Y}^{\prime})-\widehat{L}(h^{*},{\bf
Y})-\Big{[}\mathcal{L}(h^{*},\mid\bar{}{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h^{*},\bar{{\bf Y}})\Big{]}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\widehat{\mathcal{L}}(h^{*},\bar{}{\bf
Y})-\widehat{\mathcal{L}}(h^{*},{\bf Y})$
$\displaystyle\leq\frac{1}{m}\Big{[}f(h^{*},\bar{Y}_{1}^{T}(m))-f(h^{*},Y_{1}^{T}(m))\Big{]}\leq\frac{1}{m}.$

where (a) follows from the fact that ${\bf Y}^{\prime}=\bar{}{\bf Y}^{\prime}$
and the last inequality follows from the fact that $f$ is bounded by 1.

By taking the supremum over $h^{*}$, the previous calculations show that
$\phi({\bf Y})-\phi({\bar{}{\bf Y}})\leq 1/m$; by symmetry, we obtain
$\phi({\bar{}{\bf Y}})-\phi({{\bf Y}})\leq 1/m$ which proves the lemma. ∎

We now prove the main theorem.

###### Proof.

Observe that the following bounds holds

$\displaystyle\Phi_{\text{s2s}}({\bf Y})=$ $\displaystyle\mathcal{L}(h\mid{\bf
Y})-\widehat{\mathcal{L}}(h,{\bf Y})$ $\displaystyle\leq$
$\displaystyle\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}(h\mid{\bf
Y})-\mathcal{L}(h\mid{\bf
Y}^{\prime})\Big{]}+{\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}(h\mid{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h,{\bf Y})\Big{]}}.$

and so

$\Phi_{\text{s2s}}({\bf
Y})-\Delta\leq\underbrace{\sup_{h\in\mathcal{H}}\mathcal{L}(h,\mid{\bf
Y}^{\prime})-\widehat{\mathcal{L}}(h,{\bf Y})}_{\phi({\bf Y})}.$

Define $M=\max_{j}\mathbb{E}[\phi(\widetilde{C}_{j})\mid\widetilde{}{\bf
Y}^{\prime}]$. Then,

$\displaystyle\Pr$ $\displaystyle\Big{(}\Phi_{\text{s2s}}({\bf
Y})-\Delta-M>\epsilon\mid{\bf Y}^{\prime}\Big{)}\leq\Pr(\phi({\bf
Y})-M>\epsilon\mid{\bf Y}^{\prime}).$ (C.1)

By sub-additivity of the supremum, we have

$\phi({\bf
Y})-M\leq\sum_{j}\frac{|C_{j}|}{m}{\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}(h\mid{\bf
Y})-\widehat{\mathcal{L}}(h,C_{j})-M\Big{]}}$

and so by union bound,

$\Pr(\phi({\bf Y})-M\geq\epsilon\mid{\bf
Y}^{\prime})\leq\sum_{j}\Pr(\phi(C_{j})-M\geq\epsilon\mid{\bf Y}^{\prime}).$

By definition of $M$,

$\displaystyle\Pr\Big{(}\phi(C_{j})-M\geq\epsilon\mid{\bf Y}^{\prime}\Big{)}$
$\displaystyle\leq\Pr(\phi(C_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})\mid\widetilde{}{\bf
Y}^{\prime}]\geq\epsilon\mid{\bf Y}^{\prime})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\Pr(\phi(\widetilde{C}_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})\mid\widetilde{}{\bf
Y}^{\prime}]\geq\epsilon\mid{\bf
Y}^{\prime})+(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j}\mid{\bf Y}^{\prime})$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}e^{-2|C_{j}|\epsilon^{2}}+(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j}\mid{\bf
Y}^{\prime}).$

where (a) follows by applying Prop. 6 to the indicator function of the event
$\Pr(\phi(C_{j})-\mathbb{E}[\phi(\widetilde{C}_{j})\mid\widetilde{}{\bf
Y}^{\prime}]\geq\epsilon)$, and (b) is a direct application of McDiarmid’s
inequality, following Lemma 3. The notation
$\beta_{\text{{s2s}}}(I_{j}\mid{\bf Y}^{\prime})$ indicates the total
variation distance between the joint distributions of $C_{j}$ and
$\widetilde{C}_{j}$ conditioned on ${\bf Y}^{\prime}$. In particular, we have
$\mathbb{E}_{{\bf Y}^{\prime}}\beta_{\text{{s2s}}}(C_{j}\mid{\bf
Y}^{\prime})=\beta_{\text{{s2s}}}(C_{j})$.

Finally, taking the expectation of the previous term over all possible ${\bf
Y}^{\prime}$ values and summing over $j$, we obtain

$\displaystyle\Pr$ $\displaystyle(\mathcal{L}(h\mid{\bf
Y})-\widehat{\mathcal{L}}(h,{\bf
Y})-\mathbb{E}_{\widetilde{C}_{j^{\prime}}}[\phi(\widetilde{C}_{j^{\prime}})\mid\widetilde{}{\bf
Y}]\geq\epsilon)\leq\sum_{j}e^{-2|C_{j}|\epsilon^{2}}+\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j}).$

Combining this bound with (C.1), we obtain

$\displaystyle\Pr\Big{(}\Phi_{\text{s2s}}({\bf Y})-\Delta-M>\epsilon\Big{)}$
$\displaystyle\leq\sum_{j}e^{-2|C_{j}|\epsilon^{2}}+\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})$
$\displaystyle\leq
ke^{-2\min_{j}|C_{j}|\epsilon^{2}}+\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})$

We invert the previous equation by choosing
$\delta>\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})$ and setting

$\epsilon=\sqrt{\frac{\log\frac{k}{\delta-\sum_{j}(|I_{j}|-1)\beta_{\text{{s2s}}}(I_{j})}}{2\min_{j}|I_{j}|}},$

which yields that with probability $1-\delta$, we have

$\Phi_{\text{s2s}}(Z)\leq
M+\Delta+\sqrt{\frac{\log\left(\frac{k}{\delta-\sum_{j}(|I_{j}|-1)\beta(I_{j})}\right)}{2\min_{j}|I_{j}|}}.$

To conclude our proof, it remains to show that

$M\leq\mathfrak{R}_{|C_{j}|}(\widetilde{C}_{j}\mid\widetilde{}{\bf
Y}^{\prime}).$
$\displaystyle\mathbb{E}[\phi(\widetilde{C}_{j})\mid\widetilde{}{\bf
Y}^{\prime}]=$
$\displaystyle\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\mathcal{L}(h\mid\widetilde{}{\bf
Y}^{\prime})-\frac{1}{|C_{j}|}\sum_{i=1}^{m}f(h,\widetilde{Y}_{1}^{T}(i))\mid\widetilde{}{\bf
Y}^{\prime}\Big{]}$ $\displaystyle\ =$
$\displaystyle\frac{1}{|C_{j}|}\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\sum_{\widetilde{Y}_{1}^{T}\in\widetilde{C}_{j}}\mathbb{E}[f(h,\widetilde{Y}_{1}^{T})\mid\widetilde{}{\bf
Y}^{\prime}]-f(h,\widetilde{Y}_{1}^{T}(i))\mid\widetilde{}{\bf
Y}^{\prime}\Big{]}$ $\displaystyle\leq$
$\displaystyle\frac{1}{|C_{j}|}\mathbb{E}\Big{[}\sup_{h\in\mathcal{H}}\sum_{\widetilde{Y}_{1}^{T}\in\widetilde{C}_{j}}g(h,\widetilde{Y}_{1}^{T}(i))\mid\widetilde{}{\bf
Y}^{\prime}\Big{]}$

where we’ve defined

$g(h,\widetilde{Y}_{1}^{T}(i))\triangleq\mathbb{E}[f(h,\widetilde{Y}_{1}^{T}(i))\mid\widetilde{}{\bf
Y}^{\prime}]-f(h,\widetilde{Y}_{1}^{T}(i)).$

Similar arguments to those used at the end of Appendix B yield the desired
result, which concludes the proof of Theorem 4.1. ∎

##  Appendix D Generalization bounds for local models

######  Theorem 5.1.

Let $h=(h_{1},\ldots,h_{m})$ where each $h_{i}$ is a hypothesis learned via a
local method to predict the univariate time series $Z_{i}$. For $\delta>0$ and
any $\alpha>0$, we have w.p. with $1-\delta$

$\displaystyle\Phi_{\text{loc}}({\bf Z})\leq$
$\displaystyle\frac{1}{m}\sum_{i}\Delta(Y(i))+2\alpha+\sqrt{\frac{2}{T}\log\frac{m\max_{i}(\mathbb{E}_{v\sim
T(Y(i))}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)])}{\delta}}$

###### Proof.

Write

$\displaystyle\Phi(Y_{1}^{T}(i))=$
$\displaystyle\sup_{h\in\mathcal{H}}\mathbb{E}[f(h,Y_{1}^{T+1})\mid
Y_{1}^{T}]-\frac{1}{T}\sum_{t=1}^{T}f(h,Y_{t}^{t+T}(i)).$

By [18, Theorem 1], we have that for $\epsilon>0$, and $1\leq i\leq m$,

$\displaystyle\Pr(\Phi(Y_{1}^{T}(i)-\Delta(Y(i))>\epsilon)\leq$
$\displaystyle\mathbb{E}_{v\sim
T(p)}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)]\exp{\Big{(}-\frac{T(\epsilon-2\alpha)^{2}}{2}\Big{)}}.$

By union bound,

$\displaystyle\Pr(\frac{1}{m}\sum_{i}\Phi(Y_{1}^{T}(i))-\Delta(Y(i))>\epsilon)$
$\displaystyle\leq m\max_{i}(\mathbb{E}_{v\sim
T(Y(i))}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)])\exp{\Big{(}-\frac{T(\epsilon-2\alpha)^{2}}{2}\Big{)}}$

We invert the previous equation by letting

$\epsilon=2\alpha+\sqrt{\frac{2}{T}\log\frac{m\max_{i}(\mathbb{E}_{v\sim
T(Y(i))}[\mathcal{N}_{1}(\alpha,\mathcal{F},v)])}{\delta}}.$

which yields the desired result. ∎

##  Appendix E Analysis of expected mixing coefficients

######  Lemma 2.

Two AR processes $Y(i),Y(j)$ generated by (4.1) such that
$\sigma=\text{Cov}(Y(i),Y(j))\leq\sigma_{0}<1$ verify
$\beta_{\text{{s2s}}}(i,j)=\max\left(\frac{3}{2(1-\sigma_{0}^{2})},\frac{1}{1-2\sigma_{0}}\right)\sigma$.

###### Proof.

For simplicity, we write $U=Y(i)$ and $V=Y(j)$.

Write

$\displaystyle\beta=$ $\displaystyle\|P(U_{T}|{\bf Y}^{\prime})P(V_{T}|{\bf
Y}^{\prime})-P(U_{T},V_{T}|{\bf Y}^{\prime})\|_{TV}$ $\displaystyle=$
$\displaystyle{\sup_{u,v}\left|P(U_{T}\\!=\\!u)P(V_{T}\\!=\\!v)-P(U_{T}=u,V_{T}=v)\right|}$
$\displaystyle=$ $\displaystyle\sup_{u,v}\Big{|}{P(U_{T}=u\mid
U_{0}^{T-1})P(V_{T}=v\mid v_{0}^{T-1})}-{P(U_{T}=u,V_{T}=v\mid
v_{0}^{T-1},u_{0}^{T-1})}\Big{|}$ $\displaystyle=$
$\displaystyle\sup_{u,v}\Big{|}\Big{[}{P(u,v\mid
U_{0}^{T-1},V_{0}^{T-1})}+f(\sigma,\delta,\epsilon)\Big{]}-{P(u,v\mid
U_{0}^{T-1},V_{0}^{T-1})}\Big{|}$

where we’ve written $\delta=u-\Theta_{i}(U_{0}^{T-1})$ (and $\epsilon$
similarly for $v$), and we’ve defined

$\displaystyle f(\sigma,\delta,\epsilon)$
$\displaystyle={P(u|U_{0}^{T-1})P(v|V_{0}^{T-1})}-{P(u,v|U_{0}^{T-1},V_{0}^{T-1})}$
$\displaystyle=e^{-\frac{1}{2}(\delta^{2}+\epsilon^{2})}-\frac{1}{1-\sigma^{2}}e^{-\frac{1}{2}\frac{1}{1-\sigma^{2}}(\delta^{2}+\epsilon^{2}-2\sigma\epsilon\delta)}.$

Assuming we can bound $f(\sigma,\delta,\epsilon)$ by a function $g(\sigma)$
independent of $\delta,\epsilon$, we can then derive a bound on $\beta$.

Let $x=\sqrt{\delta^{2}+\epsilon^{2}}$ be a measure of how far the AR process
noises lie from their mean $\mu=0$. Using the inequality

$|\delta\epsilon|\leq\delta^{2}+\epsilon^{2},$

we proceed to bound $|f(\sigma,\delta,\epsilon)|$ by bounding $f$ and $-f$.

$\displaystyle f(\sigma,\delta,\epsilon)$ $\displaystyle\leq
e^{-\frac{1}{2}(\delta^{2}+\epsilon^{2})}-e^{-\frac{1}{2}\frac{1}{1-\sigma^{2}}(\delta^{2}+\epsilon^{2}+2\sigma|\delta\epsilon|)}$
$\displaystyle\leq
e^{-\frac{1}{2}x^{2}}-e^{-\frac{1}{2}\frac{1}{1-\sigma^{2}}(1+2\sigma)x^{2}}$
$\displaystyle\leq
e^{-\frac{1}{2}x^{2}}\Big{(}1-e^{-\frac{1}{2}\frac{2\sigma+\sigma^{2}}{1-\sigma^{2}}x^{2}}\Big{)}$

Using the inequality $1-x\leq e^{-x}$, it then follows that

$\displaystyle f(\sigma,\delta,\epsilon)$ $\displaystyle\leq
e^{-\frac{1}{2}x^{2}}(1-(1-\tfrac{1}{2}\tfrac{2\sigma+\sigma^{2}}{1-\sigma^{2}}x^{2}))$
$\displaystyle\leq\tfrac{1}{2}\tfrac{3}{1-\sigma^{2}}\sigma
x^{2}e^{-\frac{1}{2}x^{2}}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\tfrac{3}{e(1-\sigma^{2})}\sigma$
(E.1)

where inequality $(a)$ follows from the fact that $y\to ye^{-y}$ is bounded by
$1/e$.

Similarly, we now bound $-f$:

$\displaystyle-f(\sigma,\delta,\epsilon)$
$\displaystyle\leq\frac{1}{1-\sigma^{2}}e^{-\frac{1}{2}\frac{1}{1-\sigma^{2}}(\delta^{2}+\epsilon^{2}-2\sigma|\epsilon\delta|)}-e^{-\frac{1}{2}(\delta^{2}+\epsilon^{2})}$
$\displaystyle\leq\frac{1}{1-\sigma^{2}}e^{-\frac{1}{2}\frac{1-2\sigma}{1-\sigma^{2}}x^{2}}-e^{-\frac{1}{2}x^{2}}$
$\displaystyle\leq\frac{1}{1-\sigma^{2}}e^{-\frac{1}{2}(1-2\sigma)x^{2}}-e^{-\frac{1}{2}x^{2}}.$

One shows easily that this last function reaches its maximum for
$x_{0}^{2}=\frac{1}{\sigma}\log(\frac{1-\sigma^{2}}{1-2\sigma})$, at which
point it verifies

$-f(\sigma,x_{0})=\frac{2\sigma}{1-2\sigma}e^{-\frac{1}{2\sigma}\log(\frac{1-\sigma^{2}}{1-2\sigma})}\leq\frac{2\sigma}{1-2\sigma}$
(E.2)

Putting (E.1) and (E.2) together, we obtain

$\displaystyle|f(\sigma,\delta,\epsilon)|$
$\displaystyle\leq\sigma\max\left(\frac{3}{e(1-\sigma^{2})},\frac{1}{1-2\sigma}\right)$
$\displaystyle\leq\max\left(\frac{3}{2(1-\sigma_{0}^{2})},\frac{1}{1-2\sigma_{0}}\right)\sigma$

Taking the expectation over all possible realizations of ${\bf Y}^{\prime}$
yields the desired result. ∎

###### Proof.

Recall that ${\bf Y}$ contains $m^{\prime}=mT$ examples, which we denote
$Y_{t-p}^{t}(i)$ for $1\leq i\leq m$ and $1\leq t\leq T$ (when $t-p<0$, we
truncate the time series approprietly). We define

$\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y})=\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}[L(h(Y_{T-p+1}^{T}(i)),Y_{T+1}(i))\mid{\bf
Y}]$ $\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y}^{\prime})=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}[L(h(Y_{t-p}^{t-1}(i)),Y_{t}(i))\mid{\bf
Y}^{\prime}]$
$\widehat{}\mathcal{L}_{\text{hyb}}(h)=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{T}\sum_{t=1}^{T}L(h(Y_{t-p}^{t-1}(i)),Y_{t}(i))$

where we note that here ${\bf Y}^{\prime}$ indicates each of the $mT$ training
samples excluding their last time point.

Observe that the following chain of inequalities holds:

$\displaystyle\Phi_{\text{hyb}}({\bf Y})=$
$\displaystyle\sup_{h\in\mathcal{H}}\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y})-\widehat{}\mathcal{L}_{\text{hyb}}(h)$ $\displaystyle\leq$
$\displaystyle{\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y})-\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y}^{\prime})\Big{]}}+{\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y}^{\prime})-\widehat{}\mathcal{L}_{\text{hyb}}(h,{\bf Y})\Big{]}}$
$\displaystyle\leq$
$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y})-\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}_{\mathcal{D}^{\prime}}[L(h(Y_{t-p}^{t-1}(i)),Y_{t}(i))\mid{\bf
Y}^{\prime}]\Big{]}$
$\displaystyle+\sup_{h\in\mathcal{H}}\Big{[}\mathcal{L}_{\text{hyb}}(h\mid{\bf
Y}^{\prime})-\widehat{}\mathcal{L}_{\text{hyb}}(h,{\bf Y})\Big{]}.$

and so

$\Phi_{\text{hyb}}({\bf
Y})-\frac{1}{T}\sum_{t}\Delta_{t}\leq\underbrace{\sup_{h\in\mathcal{H}}\mathcal{L}_{\text{hyb}}(h,\mid{{\bf
Y}^{\prime}})-\widehat{}\mathcal{L}_{\text{hyb}}(h,{\bf Y})}_{\phi({\bf Y})}.$

Then, following the exact same reasoning as above for $\Phi_{\text{s2s}}$
shows that for $\delta>0$, we have with probability $1-\delta/2$

$\displaystyle\Phi_{\text{hyb}}({\bf Y})$
$\displaystyle\leq\underbrace{\max_{j}\widehat{}\mathfrak{R}_{\widetilde{C}_{j}}(\mathcal{F})+\frac{1}{T}\sum_{t}\Delta_{t}+\sqrt{\frac{\log\left(\frac{2k}{\delta-\sum_{j}(|I_{j}|-1)\beta(I_{j})}\right)}{2\min_{j}|I_{j}|}}}_{B_{1}}$

However, upper bounding $\Phi_{\text{hyb}}$ can also be approached using the
same techniques as Kuznetsov and Mohri [18], which we now describe. Let
$\alpha>0$. For a given $h$, computing $\mathcal{L}_{\text{hyb}}(h,{\bf Y})$
is similar in expectation to running $h$ on each of the $m$ time series,
yielding for each time series $Y_{T-p+1}^{T}(i)$ the bound

$\displaystyle\mathbb{E}$
$\displaystyle[L(h(Y_{T-p+1}^{T}(i)),Y_{T+1}(i))\mid{\bf Y}]$
$\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}L(h(Y_{t-p}^{t-1}(i)),Y_{t}(i))+\Delta({\bf
Y}_{i})+2\alpha+\sqrt{\frac{2}{T}\log\frac{\max_{i}\mathbb{E}_{v\sim T({\bf
Y}_{i})}[\mathcal{N}_{i}(\alpha,\mathcal{F},v)]}{\delta}}$

and so by union bound, as above, we obtain with probability $1-\delta/2$

$\displaystyle\Phi_{\text{hyb}}({\bf Y})\leq$
$\displaystyle\frac{1}{m}\sum\Delta({\bf
Y}_{i})+2\alpha+\sqrt{\frac{2}{T}\log\frac{2m\max_{i}\mathbb{E}_{v\sim T({\bf
Y}_{i})}[\mathcal{N}_{i}(\alpha,\mathcal{F},v)]}{\delta}}$ $\displaystyle\leq$
$\displaystyle B_{2}$

We conclude by a final union bound on the event $\\{\Phi_{\text{hyb}}({\bf
Y})\geq B_{1}\cup\Phi_{\text{hyb}}({\bf Y})\geq B_{2}\\}$, we obtain with
probability $1-\delta$,

$\Phi_{\text{hyb}}({\bf Y})\leq\min(B_{1},B_{2})$

∎