Multivariate Probabilistic Time Series Forecasting with Correlated Errors
Authors: Vincent Zhihao Zheng
Lijun Sun
Venue: NeurIPS 2024 poster (2024)
Keywords: probabilistic forecasting
time series
error correlation
Primary area: probabilistic_methods
TL;DR: This paper proposes an efficient method to train multivariate probabilistic time series forecasting models with temporally-correlated errors.

Abstract: Accurately modeling the correlation structure of errors is critical for reliable uncertainty quantification in probabilistic time series forecasting. While recent deep learning models for multivariate time series have developed efficient parameterizations for time-varying contemporaneous covariance, but they often assume temporal independence of errors for simplicity. However, real-world data often exhibit significant error autocorrelation and cross-lag correlation due to factors such as missing covariates. In this paper, we introduce a plug-and-play method that learns the covariance structure of errors over multiple steps for autoregressive models with Gaussian-distributed errors. To ensure scalable inference and computational efficiency, we model the contemporaneous covariance using a low-rank-plus-diagonal parameterization and capture cross-covariance through a group of independent latent temporal processes. The learned covariance matrix is then used to calibrate predictions based on observed residuals. We evaluate our method on probabilistic models built on RNNs and Transformer architectures, and the results confirm the effectiveness of our approach in improving predictive accuracy and uncertainty quantification without significantly increasing the parameter size.

Multivariate Probabilistic Time Series Forecasting with Correlated Errors
Vincent Zhihao Zheng McGill University Montréal, QC, Canada
Lijun Sun" McGill University Montréal, QC, Canada 1ijun. www.mssnnmngggl..mmm
zhihao. thengtmaii.mcgilll
Abstract
Accurately modeling the correlation structure of errors is critical for reliable un- certainty quantification in probabilistic time series forecasting. While recent deep learning models for multivariate time series have developed efficient parameter- izations for time-varying contemporaneous covariance, but they often assume temporal independence of eπгогs for simplicity. However, real-world data often exhibit significant error autocorrelation and cross-lag correlation due to factors such as missing covariates. In this paper, we introduce a plug-and-plas method that learns the covariance structure of errors over multiple steps for autoregres- sive models with computational efficiency, we model the contemporaneous covariance using a low- aukkplus-diagonal parameterization and capture cross-covariance through group independent latent temporal processes. The learned covariance matrix is then used to calibrate predictions based observed residuals. We evaluate method on probabilistic models built on RNNs and Transformer architectures, and the results confirm the effectiveness of our approach in improving predictive accuracy and uncertainty quantification without significantly increasing the parameter size.
errors. To ensure scalable inference and
Introduction
Uncertainty quantification is crucial in time series forecasting, especially for applications that need more detailed insights than point forecasts. Probabilistic time series forecasting with deep learning (DL) has attracted attention for its ability to capture complex, nonlinear dependencies and provide the probability distribution of target variables 2 multivariate time series, autoregressive models are widely used for probabilistic forecasting B 51. modeling the joint one-step-aheat predictive distribution and generating learning, these models often assume that errors are independent over time. Typically, time series variables follow a Gaussian distribution z, = m(h)) + 7, where m(.) the mean function and 7) (0, Σ,) is stochastic error process with contemporaneous covariance matrix Σ, The assumption of time-independence implies Cov(n,, 7,) = 0, V's * t. This holds when the model account for all correlations between successive time steps through hidden states determined by previous values. However, real-world data often violate this assumption, as residuals exhibit substantial cross-correlation due to omission of important covariates and model misspecification.
callisttepahhel predictions
rolling manner. To enable scalable
Modeling error autocorrelation (or coosecoorrllttion) is key area of research in statistical time series models. A common approach is to assume that the error series follows a dependent temporal process, such as an autoregressive integrated moving average (ARIMA) model [6]. Deep learning models face similar challenges. Previous studies have attempted to incorporate
errors into the training process by modifying the loss function 7 g However, these methods, based on
Conresponding author.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

deterministic output, are not easily applicable to probabilistic forecasting models, particularly in multivariate setting. A notable innovation the batch training method introduced by Zheng et al. [9]. which trains univariate probabilistic forecasting model using generalized least squares (GLSS loss over batched errors. This approach parameterizes dynamic covariance matrix to capture error autocorrelation, which is then used to calibrate the predictive distribution of time series variables. While this method consistently improves probabilistic forecasting performance compared to naive training (i.e., without considering autocorrelated errors); however, it is only applicable to univariate models, such as DeepAR 10].
In this paper, we introduce an efficient method for learning error cross- cors-correlltion in multivariate probabilistic forecasting models. Our focus on deep learning models that are autoregressive with Gaussian-distribuuee distributed errors. Modeling cross- cross-correlation in multivariate models presents challenges due to increased dimensionality, as the covariance matrix scales with the number of time series N. To address this computational challenge, we propose characterizing error cross-correlation through set of independent latent temporal processes using low-rank parameterization of the covariance matrix. This approach prevents the computational cost from growing with the number of time series. Our method offers a general-purpose approach to multivariate probabilistic forecasting models, offering significantly improved predictive accuracy.
Contributions:
1. We introduce a plug-and-pp aag-anddplay method for training autoregressive multivariate probabilistic fore- casting models using a redesigned GLS loss. ($4)
2. We propose an efficient parameterization of the error covariance matrix across multiple steps,
enabling efficient computation of its inverse and determinant through matrix inversion and deter- minant lemmas. ($4.1)
The learned covariance matrix is used to fine-tune the predictive distribution based on observed residuals. ($4.2)
4. We demonstrate that the proposed method effectively captures error cross-correlation and improves prediction quality. Notably. these improvements are achieved through a statistical formulation without significantly increasing the size of model parameters. ($5)
2 Probabilistic Time Series Forecasting
Denote Z: [z,,r.... a E IR-V as the vector of time series variables at time step f, where N is the number of time series. Probabilistic time series forecasting can be formulated as estimating the joint conditional distribution p (ZT+I- I:T4 2T- - P+T++)) given the observed history {z: ) where ktystz = | 2t,] and x, are known time-dependent covariates (e.g.. time of day, day of week) for all future time steps. In essence, the problem involves predicting the time series values for Q future time steps using all available covariates and P steps of historical time series data:
II T 1+Q
p (ZT+17=-Q RT-P+1T: 1
(=T+1 ( 24-P--11:
(1)
which becomes autoregressive model that can be used for either one-step-ahead (Q 1) multistep-aheut forecasting in rolling manner. When performing multistep- forecasting, samples are drawn in the prediction range (1 2 1) and fed back for the next time step until the end of the desired prediction range. In neural networks, the conditioning information is commonly encoded into state vector h, Hence, Eq. (1 be expressed more concisely:
(ZT41:44 Xy X = II T+Q
p (zz hy),
(2)
where h is mapped to the parameters of a parametric distribution (c.g., multivariate Gaussian). Existing autoregressive models typically assume that the error at each time step is independent, meaning that z, follows a multivariate Gaussian distribution:
h, N' (μ(h;)) E(h,)),
where ua(-) and E(-) map h, to the mean and covariance parameters of multivariate Gaussian distribution. This formulation can be decomposed as 2, = 14 + n, with n; NN E,). The temporally independent error assumption corresponds to Cov(n,: m7) 0 for any time points and

---

Cov(nr
Cov(y- 1.
Cov(n- 2. m))
Cov(y -x 11)
0.6 0
0.6
0
0.6
0.4
0.4
0.4
10
10
10
10
0.2
0.2
0.2
0.2
15
15-
15
15
5
15
0.0
0
10
15
0.0
5
10
15
0.0
0
10
15
0.0
Figure 1: Contemporaneous' covariance matrix Cov 7.)) and sToss-covariance matrix Cov(N:-A- m.). A = 1, 2, calculated based on the one-step-ahead prediction residuals of GP- Var batch of time series from the m4 o4_hourly dataset. For visualization elarity, covariance clipped to the range [0, 0.6].
t where of t. Fig. provides an empirical example of the contemporaneous covariance matrix (T, 7.) and covarcoov matrix Cov(7)-A: 77), A 1,2, The results are calculated based on the prediction residuals of GPVar [3] on the m4 jo_uuurly dataset. While multivariate models primarily focus contemporaneous covariance, the residuals clearly exhibit temporal dependence, as Cov(1) 7i) + This non-zero cross-covariance suggests that residuals still contain valuable information, which can be leveraged to improve predictions.
Related Work
3.1 Probabilistic Time Series Forecasting
Probabilistic forecasting aims to model the probability distribution of target variables, unlike de- terministic forecasting, which produces only point estimates. There are two main approaches: parametric probability density functions (PDFs) and quantile functions (2]. For example, MQ-RNN M generates quantile forecasts using seguence-to-sqquence sequence (Seq2Seq) RNN architecture. In contrast, PDF-based approaches assume specific distribution (e.g., Gaussian, Poisson) and use neural networks generate the distribution parameters. DeepAR [10], for instance, uses an RNN to model hidden state transitions, while its multivariate version, GPVar [3] employs Gaussian copula to transform observations into Gaussian variables, assuming a joint multivariate Gaussian distribution.
Neural networks can also generate probabilistic model parameters. The deep state space model (SSM) [12] uses an RNN to generate SSM parameters. The normalizing Kalman filter (NKF) [33 combines normalizing flows (NFs) with the linear Gaussian state space model (LGM) model nonlinear dynamics and evaluate the PDF of observations. NKF uses RNNs to produce LGM parameters at each time step, then transforms the LGM output into observations using NFs. Wang al. 14 proposed the deep factor model, which includes deterministic global component parameterized by an RNN and random component from any classical probabilistic model (e.g.. Gaussian white noise) represent random effects. Some methods improve expressive conditioning for probabilistic forecasting by using Transformer instead RNNs to model latent state dynamics, thus breaking the Markovian assumption in RNNs [15.. Other approaches adopt more flexible distribution forms, including normalizing flows 4. diffusion models (5), and copulas [16 For a recent and comprehensive review, we refer readers to Benidis et al. [2]
3.2 Modeling Correlated Errors
Error correlation in time series has been extensively studied in econometrics and statistics [18, In multivariate time series, correlation structure characterized by contemporaneous covariance Var(n.) = Cov(n, n.) and cross cross-covariance Cov Cov(n:-A- 7T). Cross-covariance includes both the autocovariance of errors Cov(mj.t->> M).+) and the cross-lag covariance Cov(n.,-A 7T.{) between pairs of components in the multivariate series. Contemporaneous covariance captures the correlation among individual time series at a specific point in time. In the univariate setting, DeepAR [10) achieves probabilistic forecasting by modeling the contemporaneous covariance, assuming that emoor are independent over time. address autocorrelation, Sun al. [7] re-parameterized input and output of neural networks to model first-order emror autocorrelation, effectively capturing serially

---

correlated errors using an AR(1) process. This method improves the performance of one- mne-step-ahead neural forecasting models, allowing joint optimization of base and error regressors, but is limited to deterministic models. In spatial modeling, Saha et al. (20 introduced the RF-GLS model, which uses random forests to estimate nonlinear covariate effects and Gaussian processes (GP) to model spatial random effects. The RF-GLS model assumes that the error process follows an AR(p) process to accommodate autocorrelated errors. Zheng et al. [ proposed training a probabilistic forecasting model with GLS loss that explicitly models the time-varyigg autocorrelation of batched error terms, extending De to incorporate autocorrelated errors.
In the multivariate setting, most existing work focuses on modeling contemporaneous covariance, assuming that n, is independently distributed, which implies Cov(N): A mi
For example, GP [3] generalizes DeepAR [TO] to account for correlations between time series by viewing the distribution of time series variables as a Gaussian process. In Seq2Seq models, correlations can span across series and forecasting steps, as predictions for future time steps are generated simultaneously. Since predictions for future time steps are generated simultaneously, we refer to these correlations as contemporaneous correlations within the scope of this study. Choi et al. [21] introduced dynamic mixture of matrix Gaussian distributions to capture contemporaneous covariance of errors in Seq2Seq models. One exception that explicitly models error cross-correlation is (8, where the authors assume that the matrix- -aarii-vvriitte variate error term of a multivariate Seq2Seq model follows a matrix autoregressive (AR) process with seasonal lags. However, applying this technique to probabilistic forecasting models is not straightforward.
To the best of our knowledge, our work is the first to model cross-covariance in multivariate proba- bilistic time series forecasting. The closest related studies are by Zheng et al. g and Zheng et al. []. Zheng et al. [9] applies GLS loss in the temporal domain to model autocorrelated errors, but their approach is tailored for univariate time series. Zheng et al. [8) models cross-covariance in multivariate forecasting models, but their method is limited to deterministic models and requires predefined seasonal lags in the error autoregressive process, Our work extends p to the multivariate setting, enabling the modeling of the correlation structure of multivariate errors across multiple steps. In addition, we distinguish our approach from methods that directly model the distribution of time series variables, such as Copulas [16 17). where decomposition of the error term is provided.
4
Our Method
Our methodology builds upon the formulation outlined in Eq. 2 employing autoregressive model as its foundational framework. Using an RNN as an example, a probabilistic forecasting model consists of two components. Firstly, it incorporates a transition model Se to capture the dynamics of state transitions h, = fe (h,-1, Zt- x)). thus inherently having autoregressive properties. Second, it integrates distribution head, represented by 6, which maps h, to the parameters of the desired probability distribution. Following GPV [3, our approach employs the multivariate Gaussian distribution as the distribution head. The time series variable can be decomposed into deterministic mean component and a random error component 2, = H4 + 7h where n. ~ (O, E,). To efficiently model the covariance Σ, for large N, GPVar adopts E, L,L, diag (d:), where L, NXR (R < N) and d, R. Autoregressive models based on Gaussian likelihood typically assume that n, are independently distributed following a multivariate Gaussian distribution. The log-likelihood of the distribution serves as the loss function for optimizing the model:
parameterization
T 1 I E.l + n..T I..\mm) L = Σ 1=1 log (ze (b:)) Σ
(4)
1== 2
The parameters ((h:) are parameterized as (H) L,, de), where μe E RN represents the mean vector of the distribution. L, and d, correspond to the covariance factor and diagonal elements in the low- rank parameterization of the multivariate Gaussian distribution. We use shared mapping functions for all time series:
(bhh,,)) = -a(bee) = ww has
db(,..)) =d(bu,)) = =log(1 + exp(w] h),)))
(5)
4(b,,) = -i(har) =Wiba
4

---

Cross-correlation window
RxD
ExD
D
vec
hat
ver()
BxQ
A sample slice
time series
RxD
Ky-Deit
Figure 2: Graphic illustration of Eq. ( where the number time series in a batch, the rank of the covariance factor, D is the time window consider crosscoorrelaton, P and are the conditioning range and prediction range. Cross-correlation is modeled by introducing correlation in each of matrix Te-D+1:t-
where hir RR Wn RH Wa E and Wi Ri are parameters. Since the parameters of the mapping functions are shared across all time series, we can use random subset of time series to compute the Gaussian likelihood-basse loss in each optimization step, any subset of z will still follow a multivariate Gaussian distribution. In other words, we can train the model with a substantially reduced batch size N.
4.1 Training with Correlated Errors
We build upon the approach introduced in (9] to address cross-correlated errors in a multivariate context by introducing dime-dependent error terms 7,, into the GLS loss. In many existing deep probabilistic forecasting models, such as GPVar B1. a training batch typically consists of a sample slice of B time series spanning a temporal length of P Q, where P is the conditioning range and Q is the prediction range. The Gaussian likelihood is evaluated independently at each time step within the prediction range through one-step-ahead predictions. However, this approach overlooks the serial correlation of errors across consecutive time steps. To address this limitation, we propose modifying the likelihood function by introducing dynamic covariance that accommodates the To achieve this, we organize D smaller slices of time series with temporal length of + (i.e., Q 1), sorted by the prediction start time in sequential order, where D represents the time horizon over which we consider cross- cooss-oorrelatin.. The new batch structure effectively reconstructs the conventional training batch, covering the same time horizon when D Q. An example of the collection of target time series variables in batch covering cross-correlation horizon D is given by
temporal dependence of the error term, as illustrated in Fig. 2
Mt-D+1 + 7h-D+1
2, D42 14-D42 n: D42
(6)
24 14, + n
where for time point t', Myr, Le and de are the outputs of the model. The covariance parameterization in GPVar corresponds
ne Le +
where Fe (0, is low-dimensional latent variable, and E:r (O diag (dr)) is an additional error independent of Ty. We denote z = vec (Z:-D+1a) E DB as the collection of target time series variables in batch, where voc(.) an operator that stacks all the columns of matrix into vector. Similarly, we define 4a E RDB y DRDR ebot e DB d:st DB and L = blkdiag((LL)) ) ) E DB where L has block diagonal structure (see Fig. 2).
The batch-wise decomposition then expressed
hall 21
hal μ1,
+ L bat hal bat
+ Et
(8)
5

---

logp (4,
bata har, , m))
luge (4 "4,
Riz-Detic
Zi. 2o0e
Rop beg-nent
Pac
hgr-Deirt
dt
34 Bag-D-11e
PAc
Figure 3: Illustration of the training process. Following G time series dimensions are randomly sampled, and the base model (e.g., RNNss is unrolled for each dimension individually (e.g.. 1, 4, followed by 1, 3, 4 as depicted). The model parameters are shared across all time series dimensions. batch of time series variables z contains time series vectors 2, covering time steps from L- D to I. In contrast to [3. our approach explicitly models dependencies over the extended temporal window from t D 1 to during training.
The default GPVar model assumes the latent variable Tr is temporally independent, meaning Cov (r) 0, V's f l. However, this assumption cannot capture the potential cross-correlation in the errors. To address this, we introduce temporal dependencies in the latent variable within batch by assuming paat ~ N (0, C, @ IR). where C, is a dynamic x D correlation matrix. This approach assumes that the rows in the matrix rt-D+1:t = Et-D+1..... .r,] are independent and identically distributed, following N (0, C.). To efficiently capture dynamic patterss over time, we follow Zheng
et al. 9 and express C₁ as a dynamic weighted sum base kernel matrices: C, = EM Wm. Km where Wmmt 2 0 (with Σn Wm,t 1) represents the weights for each component. For simplicity, we model each component Km using a kernel matrix generated from kernel function, where the (i.j)-th entry is Ký exp(- (يجا with different lengthscales 4m
(SE)
(c.g., / = 3,2.3..... ). addition, we incorporate an identity matrix into the additive structure to account for the independent noise process. This parameterization ensures that C₁ is a positive definite symmetric matrix with unit diagonals, making it a valid correlation matrix. The weights for these components are derived from the hidden state h, at each time step through a small neural network, with the number of nodes in the output layer set to M (i.e., the number of components). A softmax layer is used ensure that these weights are summed up 1. Note that the parameters of this network will be learned simultaneously with those of the base model.
Marginalizing out r in Eq.
8 we have 2( ~ (nt hat Ihher) with covariance
(L) (C, @ (L)
but + diag(d))))
is straightforward to derive that for any i,j E {0, 1,.......D D 1} and * j. the proposed model creates cross-covariance Cov (n-.;, 7:-)) C L. between times t i and f - j. which non-stationary multivariate process through varying coregionalization (22 23), a key difference is that both the coregionalization coefficient matrix L, and the temporal correlation C, are generated by a deep neural network. In this sense, our model can better characterize the empirical cross-covariance matrices of the residuals (see empirical examples in Fig. 1). As fu! L. and d, are default outputs of the base probabilistic model, we can compute the overall likelihood (with overlapped data) as
no longer O. While this parameterization results
L
logp (z)
bal
but
(10)
Here, computing the log- logliike involves evaluating the inverse and the determinant of Inha with size DB x DB, for which naive implementation has prohibitive time complexity of O (D3B) B²).
However, our parameterization of >>t1 as E + ACAT, where E diag (d(tt) A L and Morrie---oorriooooooooe
C= C, & IR. allows us to leverage the lemma) and the companion matrix determinant lemma to simplify the computation:
identity (matrix inversion
EE + ACA
E- - E-
A(C + E
A) ATE
(11)
det (E + ACAT) = del (C-1 A A A) del (C) det (E).
Then, the likelihood calculation only requires computing the inverse and determinant of a DR x DR matrix, specifically C-1 + T E These computations can efficiently performed using Cholesky factorization. Detailed computations are provided in Appendix A.2

---

iteratively apply the process described in Modeling the latent process Te offers several advantages. Firstly, because Py has a much lower dimension than E,, modeling the cross-correlation of r, results in significantly smaller DR DR covariance matrix compared to the DB x DB covariance matrix of et. Secondly, since ry follows an isotropic Gaussian distribution, the covariance of r can be parameterized with a Kronecker structure C & I Rg- This greatly simplifies the task into learning D x D correlation matrix shared by all time series in batch. Lastly, similar to GPVar, we can still train the model in an end-to-end manner using subset of time series in each iteration to ensure computational efficiency (Fig. 3
4.2 Mult Rolling Prediction
Autoregressive models perform multistep-ahees forecasting in an iterative manner, where the model generates a sample at each time step during prediction, using it as input for the subsequent step, and continuing this process until the desired prediction range reached. Our approach enhances this process, similar to Zheng et al. [9], by offering additional calibration based on the learned correlation matrix C₁. Assuming observations are available up to time step f. the conditional distribution of 7c+-1 given errors in the past (D 1) steps, can be derived as
M Th 1.......... D-2 ~ N E.Eii Tin Jobs E++1- - E.I.E..
(12)
where 7obs vec ([):- (::=D22:.... 7h- n:]) (D-18) represents the set of residuals, accessible at forecasting step + 1. Here, Eote is a (D - 1)B x (D - B)B partition of 2/41 that captures the covariance of Mobs. and E, is a B x (D - 1)B partition of 2 representing the covariance hat between 1 and nobs. i.e., -
Eebs Σ.
For conciseness, we omit the time index / in Ec-1
Σ,
Eots. E. and Tohs- Since 1441 is a deterministic output from the base model, a sample of the target variables ž-+1 can be derived by first drawing a sample Mi+1 from Eq. (I2). then combining it with the predicted mean vector M4+1 as žt+1 = P + - It should be noted that we can still leverage the -orimaann-Moriio----er Morrison-Woodbury Woodbury identity when computing the inverse E
By taking the sample 7 observed residual,
Eq. (12 to derive a trajectory of {2,+g] 1 Repeating this procedure allows us to generate multiple samples, characterizing the predictive distribution each time step.
5
Experiments
5.1 Evaluation of Predictive Performance
Datasets. We use widely recognized time series benchmarking datasets from GluonTS 24] The prediction range (Q) for each dataset follows the configurations provided by GluonTS. We applied sequential split into training, validation, and testing sets for each dataset. Each dataset was standardized using the mean and standard deviation from the training set, and predictions were rescaled their original values for evaluation. Further details the datasets be found in Appendix $
Base probabilistic models. We integrated the proposed method into two distinct autoregressive models: the RNN-based GPVar 3 and the decoder-only Transformer (25]. These models are trained to generate distribution parameters as described in T Our approach can be applied to other autoregressive multivariate models with minimal adjustments, provided the final prediction follows a multivariate Gaussian distribution. The implementation is based on using PyTorch Forecasting (26). Both models use lagged time series values and additional features or covariates as inputs. Details on model training ($ (SA.3) hyperparameter tuning ($A.5). and the base model ($A.6) are provided in Appendix $A The code is available at https:/ggthuu.ccom//// wwwm//ottenii//myy/a... liarieaivvyyyypp__c_oee cs__corrllttedrrr
Dynamic correlation matrix. We introduce limited number of additional parameters to project the state vector h; into component weights Wm,t, which are used to generate the dynamic correlation matrix Ct. The number of base kernels (M) for generating C, and the associated lengthscale set
3.4 4 and two sets of (Im m=1 are treated as hyperparameters. We perform grid search over M = 2,233.
lengthscales- (0.5, 1.5,... } and (1.0, 2.0.... ... J. Models with the best validation loss are selected. These different lengthscales capture varying correlation decay rates, enabling the model to account
7

---

Table 1: CRPS accuracy comparison. "wlo" denotes methods without time-dependent errors, while "w/" indicates our method. Bold values show models with time-dependent errors performing better. Mean and standard deviation are obtained from 10 runs of each model.
GARCH
Transfomer
w/o
w/o
exchange_ rate
0.003320.0000
0.0435a0.0001
0.006 0004
0.01174000000
0.005550000002
0.0042==0.0002
solar
0.7663=0.000
0.8752=000015
0.7103+0.0065
0.6929:0.0039
0.4960=0.0033
0.4132=00.002
electricity
0,1264=0. 0006
0.2847=0.0015
0.0430+0. 0005
0.04033=00004
0.0494+0.000 00044
0.0638+0 0003
traffic
3,5241=0.0084
0,4459=0.00 0005
0.1095#0.000 0002
0.0649:0.000 0002
0.0717=00000 0002
0.0981=0.000
wikipedis
26,202550.0388
0.669910.0055
41745400.000
0.0743200.009
0.9841100.00)
0.0500=0.0005
m4 hourly
0,2352+0.0008
0.2758=0.0006
0.0613340.000 0004
0.0358:00.000
0.065 1+0
0.0616=00.0003
qi_quarteryy
NIA
0.394240.0030
0.3538:0.0017
0.444820.0027
0,437700,0028
pens03
0.0598=0.0002
0.3202+0.0007
0.0503+0.0001
0.0491=000002
0.0490&0.0001
0.0386=00.0001
uber_ _ouurll
0.0342+0.0006
0.0222#000004
0.0612+0.0003
0.051330.0005
avg. rel. impr.
13.79%
avg. rel. impr.
6.91%
for different temporal patterns. The viry-ng component weights enable dynamic adaptation to changing correlation structures over time.
Baselines. We evaluate the proposed method by comparing it with a baseline model trained without accounting for error cooss-coreelation (Eq. ( The baseline model represents a special case of our model with C, = I D- To ensure a straightforward and fair comparison, we align the cross-correlation range (D) with the prediction range (Q), ensuring identical data sampling processes for both methods. Additionally, we set P = Q following the default configuration in GluonTS. We also include VAR and GARCH as naive baseline models (see Appendix $A.4).
Metrics. We use the Continuous Ranked Probability Score (CRPS) 27 as the main metric:
1 CRPS(F, = EF 2 - Ep 1Z 2'1,
(13)
2
where F is the cumulative distribution function (CDF) of the predicted variable, 2 the observation. Z and Z' are independent copies of the prediction samples associated with the distribution F. To evaluate multivariate dependencies in the time series data, we compute CRPS..un by first summing both the forecast and ground-truth values across all time series and then calculating the CRPS over the resulting sums 3 16 17.. As CRPS may overlook model performance on individual dimensions 288 we also report additional metrics, e.g., the energy score 27 29), in Appendix $B.1
Training dynamics. Our approach incurs additional training costs per optimization step due to the more complex likelihood function. As shown in Appendix BB.3 the training time per epoch for models using our method is generally longer than that of baseline methods. However, our parameterization allows for scalability to large time series datasets by using small random subset of time series at each optimization step during training.
Benchmark results. The CRPS results are presented in Table I Our method achieves an average improvement of 13.79% for GPVar and 6.91% % for the Transformer model. is important to note that the degree of performance enhancement varies across different base models and datasets, influenced by factors such as the inherent data characteristics and the performance of different model architectures. The alignment between the actual correlation structure and our kernel assumption also plays a erucial role in the effectiveness of our method. Additionally, our approach demonstrates consistent improvements across five different metrics, with significant gains in multivariate metrics such the energy score (Appendix (B.2).
To provide further insights, we compare the residual autocorrelation and cross-lag correlation with and without applying our method in Appendix IB.5.I showing that our method effectively reduces cross-correlations in many scenarios. We use ACF plot comparisons to illustrate the reduction in autocorrelation and cross-correlatio plot comparisons to demonstrate the decrease in cross-lag correlation. The residuals generated by the model with our method exhibit weaker cross-correlations, which is particularly enhanced by the calibration process during prediction ($4.2).
Furthermore, Appendix $B.5.2 separates the accuracy improvement over forecast steps for each dataset. The performance improvement is shown to be related to both the absolute time across the
8

---

Component Weights
ACF: 17:00
Siaut: 17:00
0.6
1.00
0.03
0.75
0.4
0.02
160
uv/l-
Mx(f-225) 2.5)
0.2
0.01
0.25
320
17:00
0.0
01-12
01-13
01-14
01-15
0.00
10
160
320
0.00
Time (date)
Figure 4: (a) Component weights for generating C. for a batch of time series (B 8) from the m4 hourly dataset obtained by the GPVar model. Parameters wo, W1, W2 represent the component weights of the kernel matrices associated with lengthscales = 0.5, 1.5, 2.5, and w3 is the component weight of the identity matrix. Shaded areas distinguish different days; (b) The autocorrelation function (ACF) indicated by the correlation matrix C2 at 17:00. Given the rapid decay of the ACF, we only plot 12 lags to enhance visualization; (c) The corresponding covariance matrix of the associated target variables > 17:00. A zoom-in view of 3B 3B region is illustrated in the plot, where the diagonal blocks represent B x B covariance matrices Ev of Zy over three consecutive time steps. The off-diagonal blocks describe the cross-covariance Co zt), VA + 0. For visualization clarity, covariance values are clipped to the range [0, 0.03].
temporal span of the dataset (especially for time series with strong periodic patterns) and the relative time over the prediction horizon.
5.2 Model Interpretation
Our method captures error cross-correlation through the dynamic construction of a covariance matrix, achieved by combining kernel matrices with varying lengthscales in dynamically weighted sum. small lengthscale corresponds to short-range positive correlations, while a large lengthscale captures positive correlations over longer lags.
In Fig. we depict the dynamic component weights and the resulting autocorrelation function (the first of the correlation matrix C.) for batch of time series from the p4_hourly dataset spanning four-day window. We also provide the covariance matrix of hat using the correlation matrix and model outputs a specific time of day. The component weight w3, corresponding the identity matrix, dominates throughout the observation period. This suggests that the eггоr correlation is generally mild over time. This behavior is influenced by the Kronecker structure used to parameterize the covariance over the low-dimensional latent variables the same autocorrelation structure. Given the Kronecker structure, the model tends learn the mildest temporal correlation among the time series in a batch.
which assumes all latent processes share
Moreover, observe that the dynamic component weights adjust the correlation strengths. Specif- ically, when the weight assigned to the identity matrix (wa) increases, the error process tends to be more independent. In contrast, when the weights assigned the other kernel matrices (wo, W1. and w2) are larger, the егтоr process becomes more correlated, as the kernel matrices with different lengthscales combine formulate specific correlation structure. Fig. aa) demonstrates pronounced daily patterns in temporal correlation, particularly when errors exhibit increased correlation around 17:00 each day. The corresponding autocorrelation function shown in Fig. bb). Fig. 4c) illustrates the corresponding covariance matrix of the associated target variables within the conssscorrlltionn horizon. The diagonal blocks represent the contemporaneous covariance E, of z, each time step, while the off-diagonal blocks capture the cross-covariance Cov (2o(2z-0 A, 2z) for VA A O, effectively modeled by our approach. The zoomed-in view provides 3B 3B region that illustrates cross-covariance within two lags. We observe that the cross-covariance is most pronounced at lag 1, consistent with the observation in Fig. () that the component weight wo, assigned to the base kernel matrix with lengthscale [ 0.5, is more pronounced than w and W2.
Discussion
In this section, we discuss factors that influence the performance of our method. Specifically, we highlight the effectiveness of our model in long-term forecasting across various scenarios. We also

---

discuss the effect of scaling up to larger batch sizes during prediction. Additionally, we examine the impact of non- non-Gaussian errois on model performance.
Long-term forecasting. The advantage of modeling error correlation can vary in long-term fore- casting, especially in autoregressive predictions where errors accumulate and propagate over time. Using residuals from previous time steps to calibrate forecasts may be beneficial for non-stationary segments of the time series. However, for time series with strong periodic effects, the model may also rely on seasonal lags. As shown in Fig. 21 and Fig. 22 of the Appendix, the advantage of modeling error correlation can decrease in longer-term forecasts compared to shorter-term forecasts for some datasets with strong periodic effects (e.g.. the traffic ic dataset in Fig. 21). It is not necessarily true that the advantage diminishes for long-horizon predictions, as the effectiveness of our method depends on the quality of predictions during inference. In cases where the model provides accurate long-term forecasts, the benefit of modeling correlated errors may be less pronouneed.
Increasing the number of time series B in a batch leads to higher training costs. Because the model requires numerous iterations over the dataset for optimization, using a large B during training is not feasible. However, during prediction, the batch size can be increased to leverage more information. This may enhance both prediction accuracy and error calibration, provided sufficient memory is available. We demonstrate the effect of increasing batch size during prediction in Appendix |B.4 through additional experiments. Both models, with and without our method, show improvement from increased batch sizes during prediction, as reflected by decrease in CRPSsum-
Scalability.
Non-Gaussian errors. For the baseline model, assuming Gaussian errors may lead to model misspec- ification, resulting in more correlated residuals. To address this issue, we also trained the baseline models using the likelihood of multivariate f-distribution; the results the Appendix. Although using an alternative distribution can lead to better performance on some datasets without our method, we observed that our method effectively closes the performance gap when the i-distribution outperforms the Gaussian assumption. We chose the Gaussian distribution for its beneficial properties, including its marginalization rule and well-defined conditional distributions, both essential for statistically consistent model training and reliable inference. Thus, a more effective data
presented in Table 15
approach could involve first transforming the original observations into using a Gaussian Copula [3) followed by applying our method.
7 Conclusion and Broader Impacts
This paper presents a novel approach for addressing error cross-correlation in multivariate proba- bilistic time series forecasting, specifically for models with autoregressive properties and Gaussian distribution outputs. We construct a dynamic covariance matrix using small set of independent and identically distributed latent temporal processes. These latent processes effectively model temporal correlation and integrate seamlessly into the base model, where the contemporancous covariance is parameterized by a low-rank- raan--llusddaae plus-diagonal structure. This approach enables the modeling and prediction of time-varying covariance matrix for the target time series variables. The experimental results demonstrate its effectiveness in enhancing uncertainty quantification.
Our contributions are two-fold. First, our approach relaxes the time-indepene independent error assumption during the training process for probabilistic forecasting models, addressing the reality that residuals are typically time-dependent. Second, the learned cross-correlation improves multistep-ahead predictions by refining the distribution output at each forecasting step. These enhancements to existing models have broader implications for fields such as finance, healthcare, and energy, where improved forecasts and uncertainty quantification can lead to more informed decisions.
There are several avenues for future research. First, the Kronecker structure C, 3 I for the covariance matrix of the latent variable r may be too restrictive for multivariate time series problems. Exploring more flexible covariance structures, such as employing different C, matrices for each latent temporal process as the linear model of coregionalization (LMC, [22 23)), could be promising direction for further investigation. Second, the parameterization of C₂ could be expanded. Instead of using SE kernels, C, could be parameterized as fully learnable positive definite symmetric Toeplitz matrices. For example, an AR(p) process has a covariance structure in Toeplitz form, allowing for the modeling of negative correlations. This alternative approach could offer greater flexibility in capturing complex correlation patterns in multivariate time series data.
10

---

Acknowledgments and Disclosure of Funding
We acknowledge the support from the Natural Sciences and Engineering Research Council (NSERC) of Canada (Discovery Grant). Vincent Zhihao Zheng also acknowledges the support received from the FRONT B2X Doctoral Scholarship Program.
References
Tilmann Gneiting and Matthias Katzfuss. Probabilistic forecasting. Anmual Review of Statistics and lts Application, 1(1):125-55, 2014.
Konstantinos Benidis, Syama Sundar Rangapuram, Valentin Flunkert, Yuyang Wang. Danielle Maddix, Caner Turkmen, Jan Gasthaus, Michael Bohlke-Schneider, David Salinas, Lorenzo Stella, et al. Deep learning for time series forecasting: Tutorial and literature survey. ACM Compuring Surveys, 55(6):1-36, 2022.
[3] David Salinas, Michael Bohnke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank gaussian copula processes. Advances in Neural Information Processing Systems, 32, 2019.
[4] Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Multivariate probabilistic time series forecasting via conditioned normalizing flows. In International Conference on Learning Representations, 2021.
[5] Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf. Autoregressive denoising diffusion models for multivariate prohabilistic time series forecasting. In International Conference on Machine Learning, pages 8857-8868. PMLR. 2021.
Rob J Hyndman and George Athanasopoulos. Forecasting: Principles and Practice. OTexts, 2018.
[7] Fan- Fan-Keng Sun, Chris Lang, and Duane Boning. Adjusting for autocorrelated errors in neural networks for time series. Advances in Neural Information Processing Systems, 34:29806-29819, 2021.
[8] Vincent Zhihao Zheng. Seongjin Choi. and Lijun Sun. Enhancing deep traffic forecasting models with dynamic regression. arXiv preprint 2r30v:2330..00 3..66,060,, 2023.
[9] Vincent Zhihao Zheng, Scongjin Choi, and Lijun Sun. Better batch for deep probabilistic time series forecasting. In International Conference Artificial Intelligence and Statistics, pages 91-99. PMLR, 2024.
[10] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3):1181-119) 2020.
[11] Ruofeng Wen, kari Torkkola, Balakrishnan Narayanaswamy, and Dhruv Madeka. A multi-horizon quantile recurrent forecaster. arXiv preprint arXiv:1711.1105: 2017.
[12] Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski. Deep space models for time series forccasting. Advances in Neural Information Processing Systems, 31, 2018.
[13] Emmanuel de Bézenac, Syama Sundar Rangapuram, Konstantinos Benidis, Michael Bohike-Schneider, Richard Kurle, Lorenzo Stella, Hilaf Hasson, Patrick Gallinari, and Tim Januschowski. Normalizing kalman filters for multivariate time series analysis. Advances in Neural Information Processing Systems, 33:2995-3007, 2020.
[14] Yuyang Wang. Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski. Deep factors for forecasting. In International Conference on Machine Learning, pages 6607-6617. PMLR, 2019.
[15] Binh Tang and David S Matteson. Probabilistic transformer for time series analysis. Advances Neural Information Processing Systems, 34:23592 226608, 2021.
[16] Alexandre Drouin, Étienne Marcotte, and Nicolas Chapados. Tactis: Transformer-tttennions copulas for time series. In International Conference on Machine Learning, pages 5447-5493, PMLR, 2022.
[17] Arjun Ashok, Étienne Marcotte, Valentina Zantedeschi, Nicolas Chapados, and Alexandre Drouin. Tactis-2: Better, faster, simpler attentional copulas for multivariate time series. In International Conference on Learning Representations, 2024.

---

[18] Raquel Prado, Marco AR Ferreira, and Mike West. Time Series: Modeling, Computation, and Inference. CRC Press, 2021.
[19] James Douglas Hamilton. Time Series Analysis. Princeton University Press, 2020.
[20] Arkajyoti Saha, Sumanta Basu, and Abhirup Datta. Random forests spatially dependent data, Journal of the American Statistical Association, 1186541))65-6633, 2023.
[21] Seongjin Choi, Nicolas Saunier, Vincent Zhilao Zheng. Martin Trepanier, and Lijun Sun. Scalable dynamic mixture model with full covariance for probabilistic traffic forecasting. arXiv preprint wrXiv:2112.06653. 2022.
[22] Alan E. Gelfand, Alexandra M Schmidt, Sudipto Banerjee, and CF Sirmans. Nonstationary multivariate process modeling through spatially varying coregionalization. Test, 13:263-312, 2004.
[23] Rui Meng, Braden Soper, Herbert KH Lee, Vincent Liu, John D Greene, and Priyadip Ray. Nonstationary multivariate gaussian processes for electronic health records. Journal of Biomedical Informatics, 117: 103698, 2021.
[24] Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, et al. Gluonts: Probabilistic and neural time series modeling in python. The Journal of Machine Learning Research, 21 (1):4629-4634, 2020.
[25] Alec Radford, Karthik Narasimhan, Tim Salimans, llya Sutskever, et al. Improving language understanding by generative pre-training. 2018.
[26] Jan Beitner. Pytorch forecasting. http:://poorra--freccctinn..g readhhhhcs io 2020.
[27] Tilmann Gneiting and Adrian Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477) 2007.
[28] Alireza Koochali, Peter Schichtel, Andreas Dengel, and Sheraz Ahmed. Random noise vs. state-of-thc-art probabilistic forecasting methods: A case study on crps-sum discrimination ability. Applied Sciences, 12 (10):5104, 2022.
[29] Érienne Marcotte, Valentina Zantedeschi, Alexandre Drouin, and Nicolas Chapados. Regions reliability in the evaluation of multivariate probabilistic forecasts. In International Conference on Machine Learning. pages
PMLR, 2023.
[30] Dheeru Dua and Casey Graff. https://archive.t ics. www.mii.eed/////dde........ 2017.
Uci
machine
learning
repository.
URL
[31] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos. The m4 competition: 100,000 time series and 61 forecasting methods. International Journal of Forecasting, 36(1):54-74, 2020.
[32] Guokun Lai, Wei-Cheng Chang. Yiming Yang. and Hanxiao Liu. Modeling long-and short-term temporal patterns with deep neural networks. In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 95-104, 2018.
[33] Spyros Makridakis, Allan Andersen, Robert Carhone, Rohert Fildes, Michele Hihon, Rudolf Lewandowski, Joseph Newton, Emanuel Parzen, and Robert Winkler. The accuracy of extrapolation (time series) methods: Results of a forecasting competition. Journal of Forecasting, 1(2): 121-111, 1982.
[34] Chao Chen, Karl Petty, Alexander Skabardonis, Pravin Varaiya, and Zhanfeng Jia. Freeway performance measurement system: mining loop detector data. Transportation Research Record, 1748(1):961100) 2001.
[35] Caltrans. Caltrans performance measurement system. URL https//peems.dor...gov/ 2015.
[36] NYC Taxi Limousine Commission. Uber tle foil response, 2015.
[37] Jan Gasthaus, Konstantinos Benidis, Yuyang Wang, Syama Sundar Rangapuram, David Salinas, Valentin Flunkert, and Tim Januschowski. Probabilistic forecasting with spline quantile function mns. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1901-1910. PMLR, 2019.
[38] Helmut Lütkepohl. New Introduction to Multiple Time Series Analysis. Springer Science & Business Media, 2005.
[39] Roy Van der Weide. Go-garch: multivariate generalized orthogonal garch model. Journal of Applied Econometrics, 17(5):549 (1544449-664 2002.
12

---

[40] Robert Engle. Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models. Journal of Business & Economic Statistics, 20(3):339-350, 2002.
[41] Skipper Seabold and Josef Perktold. Statsmodels: econometric and statistical modeling with python. SciPy, 7:1. 2010.
[42] Prashant Srivastava. mgarch. https///pppi..rtr/rroooevee//g/p 2022.
[43] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in Neural Information Processing Systems, 30, 2017.
[44] Lei Bai, Lina Yao, Can Li, Xianzhi Wang, and Can Wang. Adaptive graph convolutional recurrent network for traffic forecasting. Advances in Neural Information Processing Systems, 33:17804-17815, 2020.
[45] Shun-Yar Shih, Fan- Fan-Keng Sun, and Hung-yi Lec. Temporal pattern attention for multivariate time series forecasting, Machine Learning, 108:1421-1441) 2019.
[46] Massimo Guidolin and Manuela Pedio. Essentials of Time Series for Financial Applications. Academic Press, 2018.
13

---

Appendix
Table of Contents
A Experimental Details
A.1 Datasets
A.2 Multivariate Likelihood with Autocorrelated Errors
A.3 Training Procedure
A.4 Naive Baseline Description
A.5 Hyperparameter Search
A.6 Base Model Description and Input Features
B Metrics and Additional Results
19
B.1 Metric Definition
B.1.1 Continuous Ranked Probability Score
B.1.2 Quantile Loss
20
B.1.3 Energy Score
20
B.1.4 Root Relative Mean Squared Error
20
B.2 Results on Other Metrics
20
B.3 Training Dynamics
20
B.4 Effect of the Number of Time Series during Prediction
20
B.5 Additional Model Interpretation
B.5.1 Comparison of Residual Correlation
2
B.5.2 Performance Breakdown at Each Forecast Step
B.6 Alternative Parametrization of C₁
29
B.6.1 Learnable Lengthscales
29
B.6.2 Using Autocorrelations an AR(p) process
30
B.7 Alternative Error Assumptions
B.8 Qualitative Results on Forecasting
A Experimental Details
A.1 Datasets
We performed experiments diverse set of real-world datasets obtained from GluonTS 244. These datasets include:
Hourly electricity consumption data collected from a total of 370 households over the period spanning from January 2012 June 2014.
electricity (30]:
m4 je__hurly [31]: Hourly time series data from various domains, covering microeconomics. macroeconomics, finance, industry, demographics, and various other fields, are sourced from the M4-competition.
exchange_raas rate [32]: Daily exchange rate information for eight different countries span-
ning the period from 1990 to 2016.
m1 __quurterlly 33 Quarterly time series data spanning seven different domains.
pems03 [34] Traffic flow records obtained from Caltrans District 3 and accessed through
the Caltrans Performance Measurement System (PeMS). The records are aggregated at
5-minute interval.
14

---

solar 32]
Hourly time series representing solar power production data in the state of Alabama for the year 2006.
Hourly traffic occupancy rates recorded by sensors installed in the San Francisco freeway system between January 2008 and June 2008.
traffic
35
uber_hourly [36]: Hourly time series of Uber pickups in New York City spanning from
February to July 2015.
wikipedia [37] Daily page views for 2,000 Wikipedia pages spanning from January 2012
to March 2014.
These datasets are widely employed for benchmarking time series forecasting models, following their default configurations in GluonTS, including granularity, prediction range (Q), and the number of rolling evaluations. For each dataset, we performed sequential split into training, validation, and testing sets, with the temporal length of the validation set matching that of the testing set. The temporal length of the testing set was determined by considering the prediction range and the required number of rolling evaluations. For instance, the testing horizon for the traffi. dataset is computed as 24 + - 1 30 time steps. As result, the model will predict 24 steps (Q) sequentially, with distinct consecutive prediction start timestamps, also known as 7 forecast instances. In our experiments, we align the conditioning range (P) with the prediction range (Q), maintaining consistency with the default setting in GluonTS. For simplicity, we set the autocorrelation horizon D. Each
(D) to also match the prediction range (Q). Essentially, in this paper, have P = dataset was standardized using the mean and standard deviation from the training set. Predictions rescaled to their original values for computing evaluation metrics. The statistics of all datasets are summarized in Table
Table 2: Dataset summary.
Dataset
Granularity
# of time series
# of time steps
Rolling evaluation
electricity
hourly
370
5,857
24
m4 hourly exchange rate
hourly
414
1,008
48
workday
6,101
m1_ quarterly
quarterly
281
48
pems03 solar
5min
358
26,208
12
24
hourly
137
7,033
24
7
traffic
hourly
963
4,025
24
uber_ uber_hourly
hourly
8,343
wikipedia
daily
2,000
792
30
A.2 Multivariate Likelihood with Correlated Errors
The probability density function of multivariate normal distribution with autocorrelated errors, as described in $4.1 defined in Eq. (14). For simplicity, we omit the subscript and superseript bat for all notations:
f(z) (2x) B/2
exp -
(14)
We use the negative log likelihood (NLL) of an observed z as the loss function for training our model. The NLL can be calculated the negative log of the probability density function in Eq. (14).
LXLL =
(z) 2 In (z nn) x- --) + Bln (2T))
(15)
where B is the number of time series in batch. The covariance matrix is parameterized Σ RR. DBDDR is the covariance factor, C DxD
L(C @ LT E. In this parameterization, is the autocorrelation matrix, E = diag(ds d). and d e RPB are the diagonal elements. The bottleneck in evaluating this NLL lies the calculation of the inverse and determinant of Σ. Therefore, we can simplify the calculation using the Sherman-Morrison- -Mrrrioon and the companion matrix determinant lemma:
identity (matrix inversion lemma)
Σ - = E LC IR) LT)
(16)
=E--_ E (00L(((O)
-T LE
15

---

Consequently, the Mahalanobis term in Eq. 15 becomes:
-1
エーー n E 1 1) E L (C (C RIIR L T.ELLEE E
L EE
=n -n-E-a E 'n- n E (eappoon)
L L'
(17)
=n E n L E E
(Loaph) E En'n)
=1 E n k k
where
Leap L E 7, Leep is the Cholesky factor of the capacitance matrix (C& L E L The computation of & can be efficiently resolved by solving the
linear system of equations Leapk L EE' n7. Since E is diagonal matrix, the only matrix inverse we need to calculate in Eq. (T7 is (C & Ix) _1 which can be further simplified as C & I Recall that C is a D x D autocorrelation matrix. Therefore, calculating its inverse is much easier than computing the inverse Σ, which is DB DB matrix. Moreover, the computational cost does not scale with the number of time series B in a batch.
The calculation of the determinant can also be greatly simplified with our parameterization:
In Σ = In| E + (C @ IR) L
II|[(C In) L E E In C &
Inl
In
E
(18)
DR
DB
= =2 = In [Leap::
+ 2R
Σ [Lcks +
Σ
In E 8,2
where Le the Cholesky factor of the autocorrelation matrix C.
A.3
Training Procedure
Compute used All models in the paper were trained in an Anaconda environment with access to one AMD Ryzen Threadripper PRO 5955WX CPU and four NVIDIA RTX A5000 GPUs (each with 24 GB of memory).
by using B 20 time series in sample slice and a batch size of 16. Because our data sampler selects one slice of time series as a batch instead of sampling 16 slices simultaneously, we set accumulate _aaeeeerrr____ee batcheebabheeeteeeeeeeeeeeeeeee_ to 16 achieve effective batch size of 16.
Batch size We adopt the approach of Ga 3
Training loop Each epoch involves training the model on up to 400 batches from the training set, followed by computing the NLL on the validation set. Training stops when any of the following conditions are met:
A total of 10,000 gradient updates have been performed during model training.
No improvement in the best NLL value on the validation set is observed for 10 consecutive epochs.
We select the version of the model that achieved the best NLL value the validation set.
A.4
Naive Baseline Description
In this paper, we employ VAR 38 (Vector Autoregression) and GARCH (39 (Generalized Au- toregressive Conditionally Heteroskedasticity as two naive baseline models. The VAR(p) model is defined as
Z, = c A .... ApZipp + E; (0 E.),
(19)
where A, is N N coefficient matrix, and is the intercept. We use VAR model of lag I (i.e., VAR(1) model) in the experiments. The parameters of Eq. 19 are estimated using ordinary least squares (OLS), following the procedure in (88.
The GARCH model describes the conditional covariance matrix of the error term a multivariate system. Suppose the model for the conditional mean is AR(1) model:
c A + Es,
(20)
16

---

where the enror term is modeled as
E: H
(21)
where H, is an N x conditional covariance matrix, and e is an N x standard normal vector, e; ~ (O, In the experiments, covariance matrix H, is defined as
the DCC-GARCH(L, 1) model [40], where the conditional
H, D.D,DD
(22)
where D, diag (h,)*² and h, contains the variances for each time series. R. is the conditional correlation matrix the DCC-GARCH model. The parameters the DCC-GARCH model are estimated with the logeiikel function:
N N (2 21n Da + In R. o R 01
(23)
In this paper, we implement the VAR model using statsmodels and the DCC-GARCH model using mgarch ch (2).
A.5 Hyperparameter Search
The hyperparameters and training configuration largely align with those used the GPVar paper B All DL models are trained using the Adam optimizer with 12 regularization set to 1c-8, and gradients are clipped at 10.0. For all methods, we limit the total number gradient updates 10,000 and decay the learning rate by a factor of 2 after 500 consecutive updates without improvement. Table lists the parameters that are tuned, as well the hyperparameters that are kept constant across all datasets and not subject to tuning.
Table 3: Hyperparameters values that are fixed or searched over a range during hyperparameter tuning.
Hyperparameter
Value or Range Searched
learning rate
[le-4, 1c-3, 1e-2]
LSTM cells / d_model of Transformer LSTM layers/ / Transformer decoder layers n heads (Transformer) rank
[10, 20, 40]
2
2
sampling dimension dropout
20
0.01
batch size
16
To tune the hyperparameters of each model, we conduct a grid search over nine parameters on each dataset. The best hyperparameters for each base model-dataset combination are selected based on the lowest validation loss. Once the optimal learning rate and hidden size are determined, we apply the same hyperparameters to models both with and without our method.
The number of base kernels (M) for generating C, and the associated lengthscale set {4m Jammm==. are two additional hyperparameters when applying our method. The optimal values of M and (Im I are selected in a similar manner via hyperparameter search. The values of M and {mm | explored
J.M M-MI during hyperparameter tuning are shown Table 4 There are possible combinations. For example, if we set MM = 3 and choose the initial lengthscale to be 1.0, the lengthscales for generating the component kernels will be {1.0, 2.0} since the last weight corresponds the identity matrix.
A.6 Base Model Description and Input Features
The input to the base models consists of lagged time series values and generic features that encode time and identify each time series. The number of lagged values used is determined by the time- frequency of each dataset. Specifically, we use lags [1, 24, 168] for hourly data; [1, 7, 14] for daily data; and [1, 2, 4, 12, 24, 48] for data with granularity of less than one hour. For all other datasets, only use the lag-1 values.
17

---

Table 4: Hyperparameters values of our method that are searched over range during hyperparameter tuning.
Hyperparameter
Value or Range Searched
(2. possible lengthscales (Im m I..., (1.0,
number of kernels
We use generic features to represent time. For datasets with a granularity of one hour less, include features for the hour of the day and the day of the week. For daily datasets, we use the day of the week feature. Additionally, each time series is distinguished by an identifier number. All features are encoded with a single value; for example, the hour of the day feature takes values in [0, 23J. These feature values are concatenated with the RNN or Transformer input each time step to generate the model input vector y:-
As illustrated in + our method requires state vector h, each time step generate the parameters for the predictive distribution and the dynamic weights for correlation matrix kernels. We use two different neural architectures for this purpose: RNN and Transformer, both of which preserve autoregressive properties. Specifically, we use an LSTM as our base model for the RNN and decoder-only Transformer (i.e., the GPT model [25] for the Transformer. Table 5 and Table summarize the number of parameters for the Gar and Transformer models across each dataset.
Table 5: Number of parameters of the GPVar model for each dataset.
covariate embedding
IBn
distribation proj
covariance pro) (our method)
exchange_
60
6,1k
252
84
solar
3.7k
26.6k
492
164
electricity
16.5%
29.6k
492
164
traffic
72.5k
34.6k
492
164
uikipedia
200k
5,7k
132
44
nd hourly
19.7k
10.2k
252
84
nt niquuarerly
164
peas03
26.4k
34.6k
492
164
uber_hourly
9.7%
28.3k
492
164
Table 6: Number of parameters of the Transformer model for each dataset.
target proj
covariate
covariate embedding
transformer
distribution proj
covariance (our method)
exchange_rate
160
400
26.5k
492
164
solar
400
3.7k
1.8k
132
electricity
160
2.4k
16.5k
26.5k
492
164
traffic
80
1.8k
72.5k
6.8k
252
vikipedia
nd hourly
80
1.2k
19.7k
6.8k
252
84
quartarly
pens03
70
870
26.4k
1.8k
132
ubar hourly
160
2k
9.7k
26.5k
492
164
LSTM, type of RNN architecture, designed model sequences and time series data. Unlike traditional RNNs, LSTMs can learn long-term dependencies, making them effective for tasks requiring context and memory over long sequences. A decoder-only Transformer primarily used for sequence generation tasks, such as text generation, language modeling, and machine translation. It is simplified version the original Transformer model introduced by Vaswani 43 consisting of only the decoder component. The LSTM model can be formulated
18

---

fr = a (W, Vr] + by),
i, a (W, [o + b,).
C, tanh (We y1 bc),
(24)
C, = f, @ C,-1 + i, @ C₁
o, a (W, y, bo),
h, = of @ tan)
where Yr is the input at each time step. The decoder-only Transformer can be formulated as
Q = Y.Wa,
K Y.WK, K;
V Y, WWW V.
M Mask(K),
QKT
(25)
Softmax
M V,
Vdl
H, LayerNorm (Y + Z);
FFN ReLU W + b1) W + b2,
H,- LayerNorm (H, FFN),
where H2 is the output containing state vectors for all time steps, and M is a square causal mask for the sequence to preserve autoregressive properties.
B Metrics and Additional Results
B.1 Metric Definition
In this paper, we repeated the evaluation process on the testing set ten times to compute the mean and standard deviation of all metrics. Metrics calculated in each independent evaluation based on the average results from all forecast instances in the testing set. For example, the CRPS reported for traffic is the average CRPS of seven forecast instances in its testing set. 100 prediction samples were drawn for all evaluation processes.
B.1.1 Continuous Ranked Probability Score
The Continuous Ranked Probability Score (CRPS) is defined as:
CRPS (F, EF 7
EplZ Z11,
(26)
where F is the cumulative distribution function (CDF) of the predicted variable, z is the observation, Z and Z' are independent copies of set of prediction samples associated with the distribution For a single forecast instance, we calculate the average CRPS across time series and over the prediction horizon:
E,, [CRPS (Fi,, 21,0),
(27)
where we use the empirical CDF to represent F,, when predicting 4,,, Since CRPS only compares
single ground-truth value to its predicted distribution, we also calculate the CRPS B 16
[7] to
assess multivariate dependencies in the time series data. CRPS.um computed by summing both the forecasted and ground-truth -tuth values across all time series and then calculating the CRPS over the resulting sums:
E, CRPS Fr,
Ziet
(28)
where the empirical F, is obtained by summing samples across time series.
19

---

B.1.2 Quantile Loss
The Quantile Loss (p-risk) is another metric used in 10 to evaluate the performance of probabilistic forecasting:
Lp (z, 29) ( =) ((1 (11-1) Iş+> plargz),
(29)
where I is binary indicator function that equals when the condition is met, 2 represents the predicted p-quantile, and 2 represents the ground truth value. The quantile loss serves metric to assess the accuracy of given quantile p from the predictive distribution. We summarize the quantile losses over the testing set across all time series segments by computing normalized summation of these losses: EEe Lp (zz,t, 26,1)
(Ei+.1) In this paper, evaluate the 0.5-risk and the
0.9-risk following Salinas et al. [10]
B.1.3
Energy Score
The Energy Score (ES) generalizes the CRPS to evaluate distributional forecasts of vector-valued random variable and is thus another multivariate metric used in this paper:
ES (, 2)
- z 2 Z
E.
Z Z' 11
(30)
2
Z~P Z~~P
where ||z|l2 is the Euclidean norm. In this paper, 3 = 1, following [17]. Since also want to aggregate over the prediction horizon, we calculate the Frobenius norm of the matrix || in practice.
B.1.4
Root Relative Mean Squared Error
The Root Relative Mean Squared Error (RRMSE) is a metric commonly used for point forecasts RRMSE is defined as:
44
32
45
EE=Ilzz 2
RRMSE
(31)
VE Iz, z112
where 2, is obtained by taking the mean of our prediction samples, and ž is the mean value of the entire forecast instance. We use this metric to evaluate the mean prediction performance our model.
B.2 Results on Other Forecasting Metrics
We present the results for CRPS (Table the 0.5-risk (Table the 0.9-risk (Table ES (Table 10 and RRMSE (Table An "N/A" entry in the tables indicates that the naive baseline models could not be properly fitted to this dataset. We observe consistent performance improvements in the base models using our method across different evaluation metrics. Notably, in the multivariate metrie ES, our method shows significant improvement, reducing the score by an average of 5.58% for GPVar and 3.21 % for the Transformer.
B.3 Training Dynamics
we compare the training dynamics of the base models trained with and without our method. Note that the likelihood losses of the two methods are not directly comparable, even for the same dataset and base model, due to differences in the likelihood structures. We observe that while our method introduces more complexity into the likelihood function, there no evidence that significantly prolongs model convergence. On the contrary, method can speed up convergence for some datasets in terms of the training steps used. also report the training time in Table [2
In Fig. 5 and Fig. 6
B.4 Effect of the Number of Time Series during Prediction
The number of time series does not impact training, the model is trained using a random subset of B time series at a time, independent of the total number time series N. However, during prediction, the batch size can be increased beyond the training batch size of B = 20 for multistep-ahead rolling
20

---

Table 7: Comparison of CRPS accuracy. "w/o" denotes methods without time- deeeddeenndent errors, while "w/" indicates our method. Boldface values indicate that models considering time-dependent errors have better performance. Mean and standard deviation are obtained from 10 runs of each model.
VAR
GARCH
GPVrr
Transformer
w/o
w/
w/o
w/
exchange_
0.0070=0.000:
0.0438+0.0001
0.0171=0.0004
0.0141=0.0003
0.0092=0.0002
0.0081=0.0001
solar
0.9566=0.mm 0022
0.9193+0.0011
0.7097+0.004)
0.752 I= 0027
0.5981=0.0021
0.5627=0.0015
electricity
0.1548=0.c0 0003
0.2778+0.0011
0,058620 0004
0.0568:0.000
0.0665=0.00
0.0775=0.0001
traffic
19.920810.0449
0.40633200000
0.1474=00.0001
0.129620.000
0.1260=0.0000
0.131820.0001
viki
334 66.6002165 4.93936
3,035 0048
0.371 0003
0.3705=0 0004
0,3737=0 0003
0,2937=0.00 0002
-4_bourly
0.2837=0.0004
0.356720.0004
0.117420.0000 0002
0.1237=0.0002
0.130620.0000
0.118920.0002
nt aiquuaterly
N/A
NUA
0.3942=0.0033
0.353820.001 0007
0.4448
0.436740.002
penada
0.1 144:0.0001
0.3533+0.0002
0.0828=0.0000
0.0815+0.0001
0.0826=0.000
0.0735=0.0000
uber_boud per_bourly
N/A
N/A
0,1488=0.0003
0. 1466=0000002
0.1576=0.000
0. 0.176200.003
rel. impr.
3.59%
rel. impr.
3.13%
Table 8: Comparison of 0.5-risk accuracy. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method. Boldface values indicate that models considering timedeepeneent errors have better performance. Mean and standard deviation are obtained from 10 runs of each model.
GARCH
GPVar
Transformer
w/o
w/
w/o
w/
exchange_rre
0.0049 0000
0.0256*0.000 1001
0.0090000000
0.0095 0004
0.0060=0.000
0.005620 0001
solar
0.6140.0.od
0.5621+0.0008
0.499820.0022
0.524620.0016
0.4233:0.0017
0.395820.0013
electricity
0.1113=0.0000
0.2014+0.000
0.0405=0.0000
0.0397=0.0000
0.0449=0.0000
0.0505=0.000
traffic
10.2654#0.002 0268
0.2722+0.000 0002
0.0933=0.0000
0.0859=0.0001
0.0803=0.0001
0.079420.000)
viki
(7150099002257 2573
0.722520.0006 0007
0.223120.00005
0.223620.0006
0.2030g0.00025
0.045320.0003
n4 hourly
0.1992=0.0003
0.2365+0.000 0005
0.0807+0.0001
0.084990.000 0001
0.0880=0.0002
0.0808=0.000
Luquurterlyy
N/A
N/A
0.2196=0.0022
0.1948=0.cmm 0.00
0.2328=0 0000
0.2327=0.0014
pene03
0.078420.0001
202800000022
0.0568=0.000
0.057420.000 0001
1.059=00.0000
0.050620.0000
aber_hourly
N/A
0.1035
0.1013:0.0002
0.123420.0002
avg. rel. impr.
2.75%
avg. rel. impr.
3.88%
Table 9: Comparison of 0.9-risk accuracy. "w/o" denotes methods without time- -eme-deeendent errors, while "w/" indicates our method. Boldface values indicate that models considering time-dependent errors have better performance. Mean and standard deviation are obtained from 10 runs of each model.
VAR
GARCH
GPVar
Transformer
w/
w/
exchange rrateeee_rttteeee
0.002
0.0070+0 0000
0.0042+0.0001
0.0057=0.0000
0.0030=000000
0.0023+0.0000
solar
0.467620 0016
0.4393+0.0005 0008
0.1617=0.00%
0.1597=0.0000
0.2744=0.0015
0.2710=0.0015
electricity
0.041420.0003
0.074420.0000
0.0211=0.0004
0.018520.000 00022
0.028100.0002
0.036620.0002
traffic
1.017044.0400 04n5
0.168910.0001
0.066620.000
0.0580=0.0001
0.0609=0.0001
0.069820.0000
viki
74.00 74.07556+.3777 3770
1. 1..9004000000 0044
0.2136+0.0002
0.2048=0.00 0001
0.2117=0.0006
0.1 J.1764-0 0003
hourly
0.1029=0.00 0003
0.1 1330000.0003
0.0452=0.00 0002
0.0463:0.00m 0000
0.0525=0.000 2.00
0.0475=0.0000
n1 quartarly
N/A
N/A
0.3049=0.0044
0.2787=0.0023
0.3784=0.0031
0.3621=0.0037
pena03
0.0399=0.000 00000
0.1783400001
0.0817=0.000
0.031 7=0
0.0304
0.0269=0.0000
uber uber_hourly
N/A
N/A
0.0533=0.0002
0.0528=0.0001
0.0562=0.0001
0.0638=0.002
avg. rel. impr
0.22%
avg. rel. impr.
0.91%
21

---

Table 10: Comparison of ES accuracy. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method. Boldface values indicate that models considering time-dependent have better performance. Mean and standard deviation are obtained from 10 runs of each model.
GPVa
Transformer
exchange solar (x (01)
1301
0,608560.0009
0.5674a0.
0.2615500007
0:1798:0.0079
0.1428:0.0126
1.742900.0043
1.7758H)/0015
1,002-00 0095
(
1.5307+000044)
14667H0000194
electricsty (X10")
0,3172
0.4754a0:0125
traffic 10")
3,3585,00010 (X )
1.419840.002
2,4008 0020
2.2408g000015
2.2240u0 0002
2 2.256620.001
970.004242.500
2.885700.0783
0,1155a0
0.1236e0.000
1075500.004
(L.ourly(x110
1.510910.0054
5.18/910.0008
2.2729400006,
2.3611a000060
2.387720.0098
2.341020.0081
al_quarterty
3.86/6e9.0147
4,214950000222
4.159620.0225
para03 (x 10")
1 1390
5.464200.000
1053550.0010
1,075660
1.0673+01.012
aber_hourly
0.90255000001
0.8773600027
0,9377200.0055
LA806200003
avg. impr
a/g impe
3.31%
Table 11: Comparison of RRMSE accuracy. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method. Boldface values indicate that models considering dependenn errors have better performance. Mean and standard deviation are obtained from 10 runs of each model.
VAR
GARCH
GPVar
Transformer
w/o
W/
w/b
wachange_rrat rate
00004700.0000
0.0955000000
0.069940.0012
0.0501+0.0010
0,0850HO DODS
10036500D000
solar
0:955kell.000
0.819550,DKES
03835440
0.511440.D02
0.776L20
wlectricity
27320000020
0,5564-0
0.1010H0.0015
0.1203+0
teaffic
0.6912e00017
1
4.9946601000
0.5385e00.0005
1.5061#0.00EE
0.5025e0.0001 DODS
0.500 DO03
0.651000
10
6.3386
hourly
061634000012
0065-
0.30220000008
nl quaureerry
NIA
NIA
19 1005a0..1220 1246
17.027740. 0845
20 20,233360.088 0830
20 20.2700 09
pwae03
0.372700,000
0.8834-0.001
0.2796+0.0003
2877+6.8925
0.2841 +0.000
2502+0
uber hourly
0.228210.000S
avg. impe
5.48%
impe
2.855
Table 12: Training cost comparison. "w/o" denotes methods without time-dependenn errors, while "w/" indicates our method.
GPVar
Transformer
w/o
w/
w/o
w/
sec./opoch
opochs
sec./epoch
epochs
sec./epoch
epochs
sec//opoch
epochs
exchange_ras
4.60
56
200,27
39
9.73
57
206.49
41
solar
.26
132
electricity
7.44
71
119.06
94
19.38
63
103.15
65
traffic
10.50
225,20
48
28.30
247.58
100
viki
12.45
164.81
33
28.16
351.06
48
m4 hourly
7.42
189.14
43
17.81
43
355.27
65
ml _iquurreryy
pans03
11.75
62
143.34
57
36.31
78
88.05
53
uber hourly
17.28
22

---

exchange rate ___ppsp
solar_ _npps
electricity_ __pps
80
train loss
-20
25
val loss
train loss
75
60
val loss
train loss
75
val loss
-30-
50
40
33
-70
70-
25
-0
1000
0
1000
2000
step
step
step
traffic__nige nips
wiki2000 nips
m4 hourly
-50
train loss
val loss
train loss
175-
val loss
train loss
80-
-40
val loss
-60-
66
190
60-
-70-
165-
20
180
1000
500
1000
step
step
step
m1_ quarterly
pems03 flow
uber tlc _ourly
110
92.5
45
train loss
40
train loss
train loss
val loss
100
90-
90.0
46
val
val
40-
90-
30
80
44
0
1000
1000
1000
step
step
step
train_loss: |lssolosssso w/o
train_loss: loss:
val loss: w/o
val val__osss Joss: w
Figure 5: Training loss/validation loss training time the GPVar model. "wlo" denotes methods without time- dene-dependent errors, while "w/" indicates our method.
exchange_re Bratiggeeeee_____ nips
solar_ _nnnssss
electricity_ nips
-00
25
} 80
train loss
val 188
train loss
75-
75
val loss
train loss
80-
val los6
30
75
-30-
50
50
-35
-25
70-
70
1000
0
2000
1000
step
step
step
traffic_nips
wikl2000 _nt00
m4 hourly
550
190
40
train loss
60
val loss
train loss
170-
val loss
train loss
val loss
-60
165
180
60
-70-
J.
-70
Loiiter
20
2000
1000
1000
step
step
step
m1 quarterly
pems03 flow
uber tlc hourly
100
92.5
train loss
75
train loss
90-
-0.00
val loss
train loss
42.5
val loss
90
50
val
46
40.0-
-5
80
A
87.5
44
500
0
2000
0
1000
step
step
step
loss: w/o
train_ Rain_
val. loss: w/o
val. loss: w
Figure 6: Training loss/validation loss vs training time of the Transformer model. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method.
23

---

predictions. This allows more information to be utilized, potentially improving both predictions and error calibration, provided that memory capacity permits. We conducted an additional experiment to demonstrate the effect of increasing the batch size during inference (Fig.) 7
solar nips
electricity nips
traffic nips
0.70
wio
0.043
D.10
0.042
CHIPS
CRaPS...
w/o
CRPS
who
0.65
0.041
w/
0.08
w/
0.60
0.039
20
60
100
140
180
20
60
100
140
180
20
60
100
140
180
of time series
of time series
1f of time series
wiki2000 nips
m4 hourly
0.175
0.06
CItPS am
0.125
w/o
CRPS
0.05
w/o
w!
w/
0.100
0.075
20
60
100
140 180
20
60
100 140
180
time series
of time series
Figure 7: The influence of the number of time series in a batch on the performance of inference. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method. We only show some datasets here because the remaining datasets have fewer than B = 20 time series in the testing set.
B.5 Additional Model Interpretation
In this section, we provide further insights into how our method improves the base model. We illustrate these improvements by comparing the cross-correlations of the residuals from models with and without our method. Additionally, we demonstrate the performance of our method over the prediction horizon in multistep-ahead forecasting.
B.5.1 Comparison of Residual Correlation
Recall that our method models both the autocovariance of errors Cov Th.)) and the cross-lag covariance Cov(n,,-4: ().+)) between all pairs of components in the multivariate series. With the calibration process introduced in $4.2
method is expected to reduce error cross-correlations, including autocorrelation and cross-lag correlation. Here, we compare the empirical ACF of the residuals 7 of single time series i, as well the empirical cross-correlations of 7) across multiple time series.
We begin by comparing the ACF of the one-step-ahead prediction residuals with and without our method. The comparisons are provided for the following datasets: solar (Fig. electricity
(Fig. 99. traffic (Fig. 10), wiki (Fig. [T). m4 m4_hourly (Fig. [2]. pems03 (Fig. 13) and
uber hourhh (Fig. 14). We observe that autocorrelation of the residuals reduced after applying our method.
Next, we compare the cross-correlations of the one-step-ahead prediction residuals with and without our method. The comparisons are provided for the following datasets: electricity (Fig. (5). traffic (Fig. 16). wiki (Fig. 17p. m4_ hourly (Fig. T8). pems03 (Fig. 19). and uber_hourl
(Fig. 20). We also observe that the cross-correlations of the residuals are reduced after applying our method.
B.5.2 Performance Breakdown at Each Forecast Step
To investigate our performance gain each forecast step, we calculate the CRPS for each forecast step. The results are shown Fig. 2T for GPVar and Fig. [2 for the Transformer. Note that the CRPS reported in this section may have different scales compared to previous sections because
24

---

Time Series
Series
ACF
ACF
-1
lags
lags
Time Series
Time Series D
ACF
ACF
lags
lags
wio
Figure 8: ACF comparison of the one-step-ahead prediction residuals with and without our method.
The results depict the prediction outcomes generated by GPVar for four time series in the solar dataset.
Time Series
Time Series B
AGF
10
15
0
10 lags
15
lags
Time Series
Time Series D
10 lags
15
-1
0
10
15
lags
w/o
w/
Figure ACF comparison of the one-step- ahead prediction residuals with and without our method. The results depict the prediction outcomes generated by GPVar for four time series in the electricity dataset.
Time Series
Time Series B
ACF
10
15
0
10
15
lags
lags
Time Series
Time Series
ACF
ACF
10 lags
15
15
lags
wio
Figure 10: ACF comparison the one-step-ahead prediction residuals with and without our method.
The results depict the prediction outcomes generated by GPVar for four time series in the traffic dataset.
25

---

Time Series
Time Series B
ACF
ACF
15
lags
lags
Time Series
Time Series D
ACF
ACF
10 lags
15
lags
wio
Figure 11: ACF comparison of the one-step-ahead prediction residuals with and without method.
The results depict the prediction outcomes generated by GPVar for four time series in the wiki dataset.
Time Series
Time Series
ACF
10
15
-1
15
lags
lags
Series
Series
ACF
AOF
-1
10
15
-1
15
lags
lags
wio
Figure 12: ACF comparison the one-step-ahead prediction residuals with and without our method.
The results depict the prediction outcomes generated by GPVar for four time series in the m4 14_huurly dataset.
Time Series
Time Series B
ACF
ACF
10
15
0
15
lags
lags
Time Series
Time Series D
ACF
ACF
10 lags
15
0
10
15
wio
w/
Figure 13: ACF comparison the one-step-ahead prediction residuals with and without our method.
The results depict the prediction outcomes generated by GPVar for four time series in the pems03 dataset.
26

---

Time Series
Time Series B
ACF
ACF
10
15
0
10 lags
15
lags
Time Series
Time Series
ACF
0
10 lags
15
-1
lags 15
wio
Figure 14: ACF comparison the one-step-ahead prediction residuals with and without our method. The results depict the prediction outcomes generated by GPVar for four time series in the uber hourly dataset.
wild
Coviw to
Coviy >
Covla ach)
Covla
0.10
0.10
0.10
0.10
6-
0-
6-
0.00
0.00
0.00
0.00
15
0..05 15
16
16
0.10
0.10
0.10
0.10
10
15
w/
Cavfa in ())
Cavia o m.)
Covin a.1))
Covia (.)
0.10
0.10
0.10
0.10
5-
0.05
0.05
0.05 5
0.06
0.00
10
0,00 10
0.00
10
0.00
0.05
15
0.05 15
0.05 15
0.05
0.10
0.10
0.10
0.10
15
Figure 15: Cross-correlation comparison the one-step-ahead prediction residuals with and without method. The results depict the prediction outcomes generated by GPVar for four time series in the electricity dataset.
%)
Covin -.
Covin aa
- 1.
0.10
0.10
0.10
0.10
5
0.05
0.05
5
0.05
0.00
10
0.00 10
00.00
4.00
0.05 15
0.05 15
00.0 15
0005
-0.10
-0.10
10
-0.10
1's
-0.10
W/
Coviy w)
Coviy ())
Covio
Cowio (.
0-
0.10
0.10
5-
0.05
0.05
0.05
0.05
10
10
10-
10-
-0.05
-0.05 15+
00.05
0.05
15
0.10
-0.10
10
15
Figure 16: Cross-coreelation comparison of the one-step-ahead prediction residuals with and without our method. The results depict the prediction outcomes generated by GPVar for four time series in the traffic dataset.
27

---

wio
Covia 10
CovEn 2: ad
Cav(s 3: m0
Covin 90
0
0.5
0.5
0.5
0.5
5
10.0 10
+0.0
0.0
10-
0.0
15
15
15
15
-0.5
-0.5
0.5
0
10
15
5
10
15
10
15
5
10
15
Cov(a 1.)
Cov(n b.)
Cov(s s. th)
Cov(i w)
0.5
0.5
0.5
0.5
5-
5+
-0.0 10-
-0.0
-..0
10-
0.0
15-
15
15-
-0.5
-0.5
-0.5
-0.5
0
10
15
5
10
15
10
15
0
5
10
15
Figure 17: Cross-correlation comparison of the one-step-ahead prediction residuals with and without our method. The results depict the prediction outcomes generated by GPVar for four time series in the wiki dataset.
Covia 1. m)
Covia 2. a)
Cow a md
Covla 6.
0
0.1
0.1
5
0.0
0.0
0.0
10
0.0
15
-5,1 15
-0,115
-011
10
15
10
15
10
15
5
10
15
Covia
Covéa m.)
Cowfa md
Covla- (.
0.1
0.1
0.1
0.1
5
5
10
10
10-
10
10. 15
15
-0.1 16
-0.1
15
15
Figure 18: Cross-correlation comparison of the one-step-ahead prediction residuals with and without method. The results depict the prediction outcomes generated by GPVar for four time series in the m4 hourly dataset.
Covin 1. m)
Covia a. a)
Cow a m)
Covlx 6.3))
0
0.5
0.5
0
0.5
0.5
5-
5-
0.0
0.0
0.0
0.0
15
15
-0.5
-0.5
-0.5
-0.5
10
15
5
10
15
0
10
15
5
10
15
Covia
Covéa
())
Cow m)
Covlat -c
0.5
5
5-
10
10
10-
10-
15-
16
16-
15
5
15
Figure 19: Cross-correlation comparison the one-step-aheet prediction residuals with and without method. The results depict the prediction outcomes generated by GPVar for four time series in the pems03 dataset.
28

---

wio
Govia
Gavéa 8)
Cowly ad
Cavia (.1))
0.5
0.5
0.5
5
5-
10
4.00
10
0.0
10-
+0.0
10-
0.0
15-
15
15-
0..5
0.55
-0.5
-0.5
10
15
5
10
15
Wi
Covin 1.
Covia 2.00
Cowle 1))
Cowlay (.70)
0-
4.5
0.5
0-
0.5
0.5
5-
0.0
0.0
100
0.0
15
15
-0.5
-0.5
-0.5
005
10
15
5
10
15
0
10
15
5
10
15
Figure 20: Cross-correlation comparison of the one-step-ahead prediction residuals with and without our method. The results depict the prediction outcomes generated by GPVar for four time series in the uber hourly dataset.
they are not normalized. In multistep-ahead forecasting, since the predicted values used as inputs for subsequent predictions within the prediction range, the residuals accumulate the effects of inaccuracies from previous steps. Therefore, the performance improvement depends not only on our modeling of eиror correlations but also on the properties of the residuals. These properties can be influenced by the absolute and relative time of the forecast and the seasonality of the data. For data without strong scasonality, residuals tend to be larger when predicting further ahead, making error accumulation more apparent. Conversely, for data with strong scasonality, the impact of error accumulation can vary. We observe that, in most scenarios, CRPS.um is reduced at the early forecasting stages. As predictions extend further into the future, some datasets (e.g., traffic c in Fig. (1)) show decreased improvement, likely due to seasonality effects. Conversely, other datasets (e.g., wiki in Fig. 21) exhibit larger improvements further into the future, possibly because the residuals accumulate over the steps.
electr city
traffic réow
3.100
2.0/5
3.025
3,000
15
28
21
0.0
12
20
356p
584p
5t4p
Stap
0300
0.15
W
0.20
0.40
0.050
0.05
E.02
0.025
10
24
30
36
20
12
15
18
Figure 21: Step-wise CRPS accuracy of GPVar. "w/o" denotes methods without time-dependent errors, while "w/" indicates our method.
B.6 Alternative Parametrization of C₁
B.6.1 Learnable Lengthscales
In this paper, the lengthscales are fixed when generating the correlation matrix C,, and the flexibility of C; comes from dynamically generating the component weights of the kernel matrices. Making these lengthscales learnable parameters to find the optimal set of (lm 1M-1 is another approach we can explore to increase modeling flexibility. Based on the best model identified in Table we experiment with treating the lengthscales learnable parameters, jointly optimized with the base
29

---

exchange rate, nips
solar nips
wik12000 nipp
crps sum
0.008
W/
1.5
w)
0.20
W
0.15
0.006
1..0
2 0.10
0.004
0.5
12
16
20
24
0.04
12
15
18
21
24
0.00
0
4
12
16
20
24
28
Step
Stop
Step
quarterly
pams00
hourly
weo
who
0.20
w/o
0.08-
crps sum
sum
0.15
0.04
0.10
0.02
0.05
0.3
0.00
0
2
3
5
6
2
6
10
0
3
6
9
12
15
18
21
Step
Step
Step
Figure 22: Step-wise CRPS,um accuracy of Transformer. "w/o" denotes methods without time- dependent errors, while "w/" indicates our method.
model. The results are shown in Table 13 lengthscales.
We do not observe significant improvement from learnable
Table 13: Comparison CRPS accuracy. "w/o" denotes methods without time-dependent error,, while "w/" indicates method. "w/ (1)" indicates the lengthscales are learnable parameters. Boldface values indicate that models considering time-dependent errors have better performance. Mean and standard deviation obtained from 10 runs of each model.
w/o
w/ (i)
w/o
(L)
exchange_rate rate
0.0068#0.0004
0.01
0.0045=0.0001
0.0055#0 0002
0.00422:0000 2002
0. 0.0072#0. 0002
solar
0.710320.0065
0.6929=0.000
0.7727 0040
0.496010.0034
0,4132:+0.027
0.41 384000022
electricity
0.0430=0.00 0.000=000005
0.0403=00.000
00.050
0.0494#0.0004
0.0638+0. 0003
0 0.0858+0.00 0000
traffic
0.1095=0.0002
0.0649=0.00022
0.1297=0.00
0.0717+0.0002
0.0981 0002
0.0950=0
viki
0.174520.0008
0.074320.0000
0.483920.0021
0.08412100013
0.0500g0. 0005
0,0472:10,0004
=4, hourly
0.0613:0.0004
0.035%:00. 0002
0.0527=0.0003
0.0651 0004
0.061640.00D 0000
0.0355+0.00M 000033
mI quartarrly
0.3942=0.0073
0.3538=0.0017
0.3534:0.0017
0.444840
0,436740.0022
0.3709H0.0120
pena03
0.0509+0.000)
0.0491=0.0000
0.045620.0001
0.049000 0001
0.0386+0 0001
0.033020.0001
hourly
0.034220.0000
0.0222:0.000x
0.021820.0003
0.063240.0000
0,0513:0.0005
0.096920.0005
B.6.2 Using Autocorrelations of an AR(p) process
One could parameterize Ct as fully learnable, positive definite symmetric Toeplitz matrices. For instance, an AR(p) process has autocorrelation matrix with Toeplitz structure, allowing the modeling of negative correlations. This alternative approach may offer more flcxibility in capturing complex correlation patterns in multivariate time series data. The autocorrelations of (R(p) process can obtained by solving a set of equations known as the Yule- Walker equations 46 For example, if consider AR(2) process and let Pk be the autocorrelation at lag A:
=
Φ1 2 22-2-2 et-
(32)
where Φ and Φ are the coefficients. We have Ps = 1 by definition and:
Pn = 01P0 + 0221,
(33)
Ps = 61 + Φ2Pk-22 k 2.
Since Po I, we can solve for P1:
Φ1
P1
(34)
1 02
30

---

and for any & 2 2, we can solve Pk: iteratively by:
Pk = 11 22
(35)
The collection (po: Pl:.... Pk: .........1111 PD-1) } forms the first row or column a Toeplitz matrix and can be used to parameterize C₁. We perform hyperparameter search to find the best AR order based on the validation loss. As shown in Table [4] while the correlation matrix C₁ parameterized by an AR process shows promise in modeling both positive and negative correlations, it does not empirically provide an overall improvement compared to the kernel method used this paper. This may be because pross-correlations in time series are predominantly positive. However, the AR method does show significant improvements on certain datasets where the kernel method does not perform well. For example, the AR method greatly improves GPVar exchange rate and the Transformer electricity.
Table 14: Comparison of CRPS accuracy. "wwo" denotes methods without time- deme-ddeendent errors, while "w/" indicates our method. "w/ (AR)" indicates C. is parameterized by an AR process. Boldface values indicate that models considering time-dependent errors have better performance. Mean and standard deviation are obtained from 10 runs of each model.
GPVar
Transformer
(AR)
exchange_rrae
0.006R+D.0004
0.0117=0.0004
0.0051=00.0002
0.0055#0.00 00002
0.0042+00000 00022
0.008840.000 0004
solar
0.7103=0 0066
0.6929=0.000
0.5923#0.00M 0012
0.496040.000 0004
0,413240 0027
0.3362+0.002
electricity
0.043020.000 0.05
0.0403=0.0000
0,043320 0007
0.049440 0004
0.0638#0 0003
0.0252 00002
traffic
0.109520.0002
0.0649::0.0002
0.07172200000
0.0878820.0K
viki
0.1745=0.0000
0.0743=0.001P
0.2375=0.0013
0.0841#0.0013
0,0500 00005
0.051 0008
hourly
0.061.3+0.0004
0.0358=0.0002
0.0298=0.0002
0.0651+0.0004
0.0616=00.0003
0.0680+0.0000 000?
quaarterly
0. 0.3942=0.000 0030
0.3538=0.0011
0.1692=0.00/29
0.4448+#0 0027
0.4367=0.0023
0,4348+0..022
pana03
0.0509#0.0001
0.0491=0.0002
0.0787#0.0002
0.0490+0.000 0001
0.0386+0.000
0.0656++0000 0001
uber_bourly
0.0842+0.0006
0.0222=0.0004
0.0575#0.0004
0.063240.0003
0,0513+0.0005
0.0770+0.0007
B.7 Alternative Error Assumptions
A more suitable likelihood function can regularize the training process, potentially reducing residual correlations. For example, assuming the errors follow a multivariate t-distribution improves the robustness of the model outliers. Additionally, stronger base model can help produce residuals that are more independent. Based on these considerations, we designed our approach to adapt dynamically varying levels of error correlation. The weighted correlation matrix assigns greater weight to the identity matrix when the errors exhibit lower correlation.
We also trained the baseline models using the likelihood of the multivariate t-distribution, and the results are shown in Table 15 While using an alternative distribution can lead to better performance on certain datasets when our method is not applied, we observed that method effectively closes the performance gap in cases where the multivariate Gaussian assumption is outperformed by the t-distribution.
An important feature of our method the ability to use a subset of time series in each training batch for model optimization, which enhances scalability. For the multivariate (-distribution, the distribution of these subsets of z, should have the same degrees of freedom as the full distribution of zy However, since the degrees of freedom are treated an additional output of the model in each training batch, they are not guaranteed to be consistent across batches. While this is not problematic for deep learning, it violates the marginalization property of the t-distribution from statistical standpoint.
We chose Gaussian noise for its beneficial properties, including its marginalization rule and well- defined conditional distribution, both essential for statistically consistent model training and reliable inference. To address model misspecification, more effective approach could involve first trans- forming the original observations into Caussian-distributee data using Gaussian Copula 3 and then applying our method.
31

---

Table 15: CRPS accuracy comparison. "w/o" denotes methods without tine-dependen errors, while "w/" indicates our method. Bold values show models with time-dependent errors performing better. Mean and standard deviation are obtained from 10 runs of each model. "N/A" indicates that the model could not be properly fitted..
GPVar
Gaussian (wio)
Gaussian (w/)
B-distribution (w/o)
Gaussian (w/o)
Gaussian (w/
E-distribution (w/o)
exchangge_ rate
0.0068:0.000 0004
0.01 1910 0004
0.01599000000
10055=000002
0.0042 0000
0.0101=0.0003
solar
0.710320.000
0.6929s0.0039
0.4960)20 0034
0.4132+0.
N/A
electricity
0,0430=0.0005
0,0403:000000
0,0467=0.0004
0,0494+0.0008 0004
0,0638+0.000 00003
0.0466200
traffic
0.1095a0.000
0,064910.0002
0.U679
0.071720.0002
0.0981 1000
N/
vikipedia
0.1745=0.0008
0.0743+00.0009
0.07300000000
0.0841=0.0013
0.0500:0.0000 0005
0.197920.0005
nd hourly
0.061:++0.0000
0.0358=0.0002
0.0365=0.0007
0.0651 +0,0004
0.0616:20 0003
0.0665=0.0000
ni qi__uurterly
0.3942=0.0000
0.3538+0.0017
0.3550#0,0064
0,4448=0.0022
0.4367
0,4466:0 0044
peas03
0.0503=0.0000
0.0491#0.0002
0,0679+0.0002
0.0490+0.0001
0.03862t..00
1.052920.0002
uber hourly
0.0222+0.0004
0.0666200.0010
0.0632a0.0003
0.03:40:0.00004
B.8 Qualitative Results on Forecasting
In this section, we provide qualitative analysis of the actual prediction performance by visualizing the predictions.
TS:1
TS:
2
TS: 3
TS: 4
1.8-
0.18
0.90
Value
0.85
Value
1.1
1.0-
0.80-
1.6
0.75
0.15
-20
20
Time Step
Time Step
Time Step
Time Step
TS:
TS: 6
TS:
TS:
1.2
1.2
0.9-
Volue
1.1-
Velue
0.8-
1.0-
0.8
0.9
0.9
-20
-20
20
-26
20
Time Step
Time Step
Time Step
Time Step
TS:
12
1.2
0.90
1.2-
0.90
Valua
1.1-
0.85-
1.1
0.85
0.80
0.80
0.9-
0.75-
0.9
0.75-
-20
-20
20
-25
-20
Time Step
Step
Time Step
observed
prodicted
Figure 23: Visualization of forecasting results on exchange_ran paat using GPVar with our method.
32

---

TS: 1
TS:
TS:
TS:4
100
200
AA
150
Value
100-
Value
100-
Value
100
Value
50-
50-
0-
20
-20
-20
20
Time Step
Time Step
Time Step
Time Step
TS:5 5
TS:
TS:7
TS:
150
150
400-
Value
100
Value
100
Value
Valuee
200
50-
20
-20
20
Time Step
Time Step
Time Step
Time Step
TS: 9
TS:10 10
TS: 11
TS: 12
155
Value
100-
Value
100
Value
Value
50
50
Time Step
Time Step
Time Step
Time Step
observed
predicted
Figure 24: Visualization of forecasting results on solar using GP with our method.
TS:
TS:
TS:
400
Value
Value
Value
300
40-
200
20
20
-20
20
-20
20
-20
20
-20
20
Time Step
Time Step
Time Step
Time Step
TS: 5
TS:
TS:7 7
TS:
4000
2000
Valuee
Value
2000
Value
1000
2000
1000
0
-20
-20
Time Step
Time Step
Time Step
Time Step
TS:9 9
TS: 10
TS:
TS: 12
400
Value
40
Value
1000
200-
40
500
0-
-20
-20
-20
-20
20
Stap
Time Step
Step
Time Step
observed
prodicted
Figure 25: Visualization of forecasting results on electricity using GPVar with our method.
33

---

TS:
TS:4
0.3
0.16-
Value
Value
0.10-
Value
0.05
Value
0.2-
0.06
0.0
0.00-
0.00
0.0
-20
20
-20
20
-20
20
-20
20
Time Slep
Time Slep
Time Step
Time Slep
TS:
TS: 6
TS:
TS:8
0.3-
0.4
Valbe
0.1-
Value
0.1
Value
Velue
0.1
0.0-
0.0
0.0
-20
20
-20
20
-20
0
-20
20
Time Step
Time Step
Time Step
Time Step
TS:
TS: 10
TS: 11
TS: 12
0.15
0.10
0.2-
Value
Velue
Value
0.1
Velue
0.10
0.05
0.05
0.00
0.0
0.0
0.00
20
-20
20
Time Step
Time Step
Time Step
Time Step
observed
predicted
Figure 26: Visualization of forecasting results on traffic using GPVar with our method.
TS:
TS:2
TS:3
TS:
200000
8000
7500
40000
Vallee
Vallee
velue
5000
Velue
100000
20000
2000
2500
0-
0-
0-
-20
20
-25
25
-25
25
-20
Step
Step
Step
Step
16000
7500
10000
40000
Value
Value
10000
Value
Value
5000
20000
2500
0-
25
25
-20
Time Step
Time Step
Time Step
Time Step
TS:9
TS:
TS:11
7S: 12
10000
20000
Value
Value
Value
Value
5000
5000
0
-20
20
-25
25
-25
25
-20
20
Time Step
Time Step
Time Step
Time Step
observed
precicted
Figure 27: Visualization of forecasting results on wiki using GPVar with our method.
34

---

TS: 1
TS:2
TS:3
TS:
25
100
Valux
Value
Value
20
50
Time Slep
Time Step
Time Step
Time Step
TS:
TS:6
TS: 7
TS:
1000
25
Value
200
Value
500
Value
20-
Value
15-
15
Step
Step
Step
Step
M
300
200
Value
Value
200
Value
Value
100
100
20
100
-50
-50
50
-00
50
-50
Time
Time Step
Time
Time Step
observed
predicted
Figure 28: Visualization of forecasting results on m4 14_uurrly using GPVar with our method.
TS:
TS:4
800
100
100
Value
400
Value
Value
50
50
200-
0-
10
-10
10
10
-10
10
Step
Step
Time Step
Step
TS:
1000
300
Value
Value
Value
200
500
200
100
100
0-
-10
10
10
10
Step
Step
Step
Step
TS:
200
250
300
Value
Value
150
Value
200
250
150
100
200
100
-10
-10
10
-10
10
10
Time
Time
Time Step
Time
observed
prodicted
Figure 29: Visualization of forecasting results on pema03 using GPVar with our method.
35

---

TS:1
TS: 2
TS:
TS: 4
60
500
Value
10
Valu
100
0-
Time Step
Time Step
Time Step
Time Step
TS:
TS: 6
TS:
TS:
2.5
20
Value
MLAN
Volue
Value
W
Value
0.0
10-
-2.5
-20
-20
20
-20
20
Time Step
Time Step
Step
Time Step
10-
Value
20
Value
10-
Value
5-
Q
-20
20
-20
20
-20
20
-20
20
Time Step
Time
Time
Time Step
observed
predicted
Figure 30: Visualization of forecasting results on uber hourly using GPVar with our method.
36

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and/or introduction have clearly stated the claims made, including the contributions made the paper and important assumptions and limitations.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We have discussed the limitations of the work in the "Conclusion" section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
37

---

Justification:
We do not have theoretical result in this study.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibiliyy
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have fully disclosed the information needed to reproduce the main experimental results of the paper in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
38

---

Answer: [No]
Justification: The code will released after the paper accepted. However, we have provided a sufficient amount of experimental details in the Appendix.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
haapai//guiddeee/eeSuuiii bonoeeooo for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to to reproduce the results. See the NeurIPS code and data submission guidelines (https
ssicnlolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: The paper specified all the training and details necessary to understand the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We ran all of our experiments for times calculate the standard deviation.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
39

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
The paper has indicated the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: The research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We have discussed societal impacts in the last section of this paper.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
40

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
The creators or original owners of assets (e.g., code, data, models), used in the paper, have been properly credited. The license and terms of use have been explicitly mentioned and properly respected.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URI.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g.. website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
41

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/moodl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
42


Reviews:
[{"review_id": "ZL0AiqWDNn", "reviewer": "Reviewer_uzY6", "summary": "This paper introduces a method for multi-variate time series forecasting where the residual errors are not assumed to be independent and modeled to be correlated using a Gaussian process model. Standard time-series approaches assume that the errors are temporally independent, however, this assumption does not hold in many datasets due to model misspecification. It is shown experimentally that the residuals have a non-zero correlation between different time steps.\n\nThe proposed methods accounts for the correlation between the time steps using a Gaussian process approach. The correlation matrix of the GP is designed to be both temporal and spatial. Spatial correlations occur due to a low rank assumption of the correlation matrix in the spatial dimension for a particular time step. Temporal correlations occur due to a non-identical correlation matrix for the latent variables across different time steps.\n\nExperimental results on benchmark datasets show improved performance with temporal correlations than without. An extensive analysis and interpretation of the results are provided in the appendix.", "questions": "- How is the model actually learned? A full algorithm listing all the steps would be very useful.\n- Are the base models and residual GP model trained independently? Or are they trained end-to-end?\n- How is inference actually done? Let's say we have historical data up to step t, and would like to predict step t+1. The correlation between the residual errors and information about the residual errors till step t helps predict a more accurate residual for step t. Decreasing correlation over longer time ranges means that the advantage of modeling correlated errors will be small for long horizon predictions. Why don't we observe this phenomena in Figures 19 and 20?", "limitations": "Please comment on the usefulness of the method for various time ranges. It seems that the method may not be an advantage when predicting over longer time ranges.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- The presented approach is well motivated as experimentally shown in the extensive analysis. Residual errors are indeed correlated and hence this is an interesting and significant problem in the time series community.\n- The proposed method is technically sound. Gaussian processes are an ideal approach to learn correlations between dimensions and hence a natural choice to account for the temporal correlations.\n- Experimental results show improved performance on standard benchmark datasets when modeling the correlations. A full set of experiments and extensive analysis shows potential of impact.", "weaknesses": "- The proposed method may not be a completely novel approach to time series forecasting. Gaussian processes have been used for time series predictions extensively and a simple way to account for temporal correlations is the following:\n  - Learn a deterministic base model on the data, compute the residual and fit a Gaussian process on the residuals.\n  - The Gaussian process kernel is a product of two kernels, one along the time dimension and the other across time series.\n  - It is unclear how this simple baseline will perform compared to the proposed method which is complex.\n- Some strong baselines are missing from the experiments. Just to give some examples:\n  - Rangapuram, Syama Sundar, et al. \"Deep state space models for time series forecasting.\" Advances in neural information processing systems 31 (2018).\n  - Salinas, David, et al. \"DeepAR: Probabilistic forecasting with autoregressive recurrent networks.\" International journal of forecasting 36.3 (2020): 1181-1191", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's encouraging comments and insightful suggestions.\n\n# Weaknesses\n\n> The proposed method may not be a completely novel approach to time series forecasting. \n\n***Response:*** Thank you for your valuable suggestions. Our model is not a two-stage model that first fits the mean of the data and then fits the residuals. Instead, our method employs an end-to-end learning approach that dynamically outputs the hyperparameters of the Gaussian Process (GP) on the errors. This approach avoids the need for constant re-training of the GP to obtain hyperparameters over the temporal horizon in the context of a two-stage model.\n\n> Some strong baselines are missing from the experiments.\n\n***Response:*** Thank you for your feedback. We did not use Deep SSM and DeepAR as baselines because they are univariate models that do not account for inter-series relationships. Therefore, we included only multivariate baselines in this paper, following [1].\n\n[1] Salinas, David, et al. High-dimensional multivariate forecasting with low-rank gaussian copula processes. In NeurIPS, 2019.\n\n# Questions\n> How is the model actually learned? \n\n***Response:*** Thank you for your comments. Our model is trained end-to-end, similar to any deep learning model. Therefore, we did not provide a list of algorithm steps, as it is not a multi-stage training process. The base forecasting model provides a hidden state $\\mathbf{h}\\_{t}$ at each time step. In the baseline models, $\\mathbf{h}\\_{t}$ is projected onto the distribution parameters $(\\boldsymbol{\\mu}\\_{t}, \\boldsymbol{L}\\_{t}, \\mathbf{d}\\_{t})$ using a simple neural network. With our method, $\\mathbf{h}\\_{t}$ is also projected onto the weight parameters $w\\_{m,t}$ of the sum of kernels in the dynamic correlation matrix $\\boldsymbol{C}\\_t$ using an additional simple neural network. The training process remains the same, with the addition of trainable parameters in the neural networks for projecting $w\\_{m,t}$ and a modified loss function (Eq. (10)) that accounts for temporal correlation.\n\n> Are the base models and residual GP model trained independently? Or are they trained end-to-end?\n\n***Response:*** As clarified in our response to the previous question, our model is trained in an end-to-end manner.\n\n> How is inference actually done? \n\n***Response:***\nThe inference process is performed by applying Eq.(12) **recursively** until the desired prediction range is reached. For example, when predicting $\\mathbf{z}\\_{t+1}=\\boldsymbol{\\mu}\\_{t+1}+\\boldsymbol{\\eta}\\_{t+1}$:\n\n$$\n\\boldsymbol{\\eta}\\_{t+1} \\mid \\boldsymbol{\\eta}\\_{t}, \\boldsymbol{\\eta}\\_{t-1},\\ldots,\\boldsymbol{\\eta}\\_{t-D+2}  \\sim \\mathcal{N} ( \\boldsymbol{\\Sigma}\\_{\\star}\\boldsymbol{\\Sigma}\\_{\\text{obs}}^{-1}\\boldsymbol{\\eta}\\_{\\text{obs}}, \\boldsymbol{\\Sigma}\\_{t+1}- \\boldsymbol{\\Sigma}\\_{\\star}\\boldsymbol{\\Sigma}\\_{\\text{obs}}^{-1} \\boldsymbol{\\Sigma}\\_{\\star}^{\\top})\n$$\n\nwhere $\\boldsymbol{\\eta}\\_{\\text{obs}}=\\operatorname{vec}\\left(\\left[\\boldsymbol{\\eta}\\_{t-D+2},\\ldots,\\boldsymbol{\\eta}\\_{t-1},\\boldsymbol{\\eta}\\_t\\right]\\right)$ represents the observed residuals, accessible at forecasting step $t+1$. $\\boldsymbol{\\Sigma}\\_{\\text{obs}}$ is the partition of $\\boldsymbol{\\Sigma}\\_{t+1}^{\\text{bat}}$ that gives the covariance of $\\boldsymbol{\\eta}\\_{\\text{obs}}$, and $\\boldsymbol{\\Sigma}\\_{\\star}$ is the partition of $\\boldsymbol{\\Sigma}\\_{t+1}^{\\text{bat}}$ that gives the covariance of $\\boldsymbol{\\eta}\\_{t+1}$ and $\\boldsymbol{\\eta}\\_{\\text{obs}}$.\n\n**By \"recursively\", we mean that the samples of $\\boldsymbol{\\eta}\\_{t+1}$ can be used as observed residuals for the next prediction step.** Therefore, the decreasing correlation indicated by $\\boldsymbol{C}\\_t$ does not imply that as we move further from time $t+1$, $\\boldsymbol{\\eta}_{\\text{obs}}$ becomes less useful. On the contrary, $\\boldsymbol{\\eta}\\_{\\text{obs}}$ will be continuously updated as we progress. Thus, it is not necessarily the case that the advantage of modeling correlated errors will diminish for long-horizon predictions.\n\nHowever, the advantage of modeling error correlation can indeed vary over long-term forecasting. Generally, errors will accumulate and propagate to future time steps when performing autoregressive predictions. Using residuals from previous time steps can be advantageous, especially for non-stationary segments of the time series. However, this may not be the case for time series with strong periodic effects, where the model can leverage information from the seasonal lags of the data. As shown in Figures 19 and 20, the advantage of modeling error correlation decreases over time for some datasets, while for others it does not.\n\n# Limitations\n> Please comment on the usefulness of the method for various time ranges. \n\n***Response:*** As clarified in an earlier response, it is not necessarily the case that the advantage of modeling correlated errors will diminish for long-horizon predictions. The usefulness of our method in longer-range forecasting can vary based on the quality of predictions generated during the inference process. Generally, errors will accumulate and propagate to future time steps when performing autoregressive predictions. Using residuals from previous time steps can be advantageous, especially for non-stationary segments of the time series. However, this may not be the case for time series with strong periodic effects, where the model can leverage information from the seasonal lags of the data. In cases where the model provides good predictions for long-term forecasting, the advantage of our method can be less obvious.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tzrGZUFcgY", "reviewer": "Reviewer_RQB7", "summary": "The paper proposes a method to improve multivariate probability forecasting by accounting for potential temporal dependencies of the residuals. The paper introduces a dynamic covariance using a small number of latent temporal processes. The method is evaluated on standard dataset of multivariate time-series on two forecasting models, with the CRPS as a metric.", "questions": "How does the number of time-series influence your method? In particular, it could be interesting to measure the decorrelation of the residuals with your method (with an appropriate metric that averages over all the time-series), with methods without accounting for it. \n\nMinor comments. \n- l18 \"vital\", this is a strong word\n- equation (1) is not accurate. Indeed, to simplify, let's take $P=1, Q=2$, and $x = 0$ (in order to remove it). To simplify further, let's assume $z$ takes discrete values. Writing $A=\\{z_T=z^0_T\\}, B=\\{z_{T+1}=z^0_{T+1}\\}, C=\\{z_{T+2}=z^0_{T+2}\\}$, equation (1) implies the equation $\\mathbb{P}(B\\cap C | A) = \\mathbb{P}(B|A)\\mathbb{P}(C|B)$ where $\\mathbb{P}$ is the probability measure. This identity is wrong. Indeed, one just need to take $A\\cap B\\cap C = \\emptyset$ but $A\\cap B \\ne \\emptyset$ and $B\\cap C \\ne \\emptyset$.", "limitations": "The limitations of the method should be explored.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The method proposed in the paper is original and clearly motivated. \nThe paper shows improvements of the CRPS accuracy over a number of time-series for two models.", "weaknesses": "The paper partly relies on the premise that incorporating time-dependency in the residuals is crucial for providing better uncertainty estimates. Demonstrating that the new residuals generated by the paper's method are less auto-correlated than those without the method is therefore crucial. However, the paper lacks a comprehensive metric for this evaluation. The theoretical introduction could be enhanced for mathematical clarity and efficiency.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's encouraging comments and insightful suggestions.\n# Weaknesses\n> The paper partly relies on the premise that incorporating time-dependency in the residuals is crucial for providing better uncertainty estimates. Demonstrating that the new residuals generated by the paper's method are less auto-correlated than those without the method is therefore crucial. However, the paper lacks a comprehensive metric for this evaluation. The theoretical introduction could be enhanced for mathematical clarity and efficiency.\n\n***Response:*** Thank you for the suggestions. To the best of our knowledge, there is no widely used metric that can comprehensively evaluate both autocorrelation and cross-lag correlation in our case. Therefore, we chose to use the comparison of ACF plots (Figures 6-12 in the Appendix) to demonstrate the effect of decreased autocorrelation and the comparison of cross-correlation plots (Figures 13-18 in the Appendix) to show the effect of decreased cross-lag correlation. We will enhance the theoretical introduction in the next version of the paper to improve the flow.\n\n# Questions\n> How does the number of time-series influence your method? In particular, it could be interesting to measure the decorrelation of the residuals with your method (with an appropriate metric that averages over all the time-series), with methods without accounting for it.\n\n***Response:*** The number of time series does not influence the training since the model is trained based on a random subset of $B$ time series at a time, regardless of the total number of time series $N$. However, during inference, one can increase the number of time series in a batch to leverage the information from more time series, as long as memory allows. We have provided an additional experiment showing the effect of increasing the number of time series in a batch when performing inference (see Fig. 1 in the PDF for the global rebuttal).\n\nTo the best of our knowledge, there is no widely used metric that can comprehensively evaluate both autocorrelation and cross-lag correlation in our case. Therefore, we chose to use the comparison of ACF plots (Figures 6-12 in the Appendix) to demonstrate the effect of decreased autocorrelation and the comparison of cross-correlation plots (Figures 13-18 in the Appendix) to show the effect of decreased cross-lag correlation.\n\n> l18 \"vital\", this is a strong word.\n\n***Response:*** Thank you for pointing this out. We will revise the narrative in the next version of the paper.\n\n> equation (1) is not accurate.\n\n***Response:*** Thank you for your observation. The example given by the reviewer is indeed correct in a general context. However, in Eq.(1), the factorization is based on the joint distribution of Gaussian variables. Since the PDF of a Gaussian is always positive, the factorization in Eq.(1) holds true in this context.\n\n# Limitations\nThe limitations of the method should be explored.\n\n***Response:*** Methodology-wise, the limitations of our method are mainly two-fold. Firstly, there is a potential model misspecification issue by assuming Gaussian noise. A potential solution is to first transform the original observations $\\mathbf{z}\\_{t}$ into Gaussian-distributed data $\\mathbf{x}\\_{t}$ using a Gaussian Copula, and then apply our method. The second limitation comes from the parameterization of the covariance matrix. Specifically, the Kronecker structure $\\boldsymbol{C}\\_t \\otimes \\mathbf{I}\\_{R}$ for the covariance matrix of the latent variable $\\boldsymbol{r}\\_t^{\\text{bat}}$ may be too restrictive. This structure assumes that the rows (i.e., latent temporal processes) in the matrix $\\left[\\mathbf{r}\\_{t-D+1},\\ldots,\\mathbf{r}\\_{t}\\right]$ are identically distributed following $\\mathcal{N}\\left(\\boldsymbol{0},\\boldsymbol{C}\\_t\\right)$. Additionally, the parameterization of $\\boldsymbol{C}\\_t$ could be expanded. Instead of using SE kernels, $\\boldsymbol{C}\\_t$ could be parameterized as fully learnable positive definite symmetric Toeplitz matrices, allowing for the modeling of negative correlations. Computation-wise, our method can lead to increased training costs due to the larger $DB \\times DB$ covariance matrix used in the likelihood calculation.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PL7k7BWNxt", "reviewer": "Reviewer_Vi7g", "summary": "The paper proposes an approach for multivariate time series forecasting, modelling the correlation between the multivariate errors in close time steps.\nThe authors show that indeed multivariate errors in close time instants are correlated.\nThey generalise the work of Zheng et al (Aistats 2024) where the idea of modelling the correlation between errors in different time step has been used only for probabilistic univariate models.\nExperimentally, they consider two different deep  architectures suitable for multivariate time series forecasting.\nIn both cases they report better performance compared to the traditional approach which does not model the correlation between the multivariate errors in different time steps.", "questions": "* You show a better CRPS for the models trained with correlation compared to those without correlation. Could you elaborate if the improvement of CRPS is due to better point forecasts, better variance of the predictive distribution, or both?\n\n* line 169-171: could you clarify (also in the final paper) the sentence: Since the parameters of mapping functions are shared for all time series, we can view z_t as a Gaussian process assessed at points hi_t?", "limitations": "None in particular.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Generalizing  the model of the temporal correlation of error from the univariate to the multivariate case is a challenging problem.\nThe proposed approach, based on modelling the correlation of the error as the sum of kernels with different length scales,  can be applied to any deep network model suitable for multivariate forecasting.\nThe approach is tested with two different neural architectures and positive results are reported in both cases.\nThe paper is well written.", "weaknesses": "None in particular.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's encouraging comments and insightful suggestions.\n# Questions\n> Could you elaborate if the improvement of CRPS is due to better point forecasts, better variance of the predictive distribution, or both?\n\n***Response:*** In this paper, we use $\\operatorname{CRPS}\\_{\\text{sum}}$ as the main metric, while energy score, quantile loss, and RRMSE are included as additional metrics in the Appendix. The $\\operatorname{CRPS}\\_{\\text{sum}}$ and energy score primarily assess the **overall quality** of the predictive distribution. While quantile loss evaluates the accuracy of specific quantiles, denoted by $\\rho$, from the predictive distribution. Specifically, we assess the $0.5$-risk and $0.9$-risk, with $0.5$-risk being equivalent to the MAE of the mean prediction, which reflects the quality of point forecasts. The $0.9$-risk offers a snapshot of the quality of a specific quantile, which is influenced by the variance/covariance of the predictive distribution. The comparison of $0.5$-risk and $0.9$-risk using GPVar as an example is shown in Table 1. We observe that the results vary across different datasets. In some cases, both $0.5$-risk and $0.9$-risk are improved. In other instances, our method improves either the point forecasts or the higher quantile risk. The outcomes can depend on the characteristics (e.g., noise level) of each model and dataset. For GPVar, our method improves $0.5$-risk on 5 out of 9 datasets and $0.9$-risk on 6 out of 9 datasets. **It is important to note that $0.5$-risk and $0.9$-risk provide only two snapshots of the quality of the predictive distribution and are thus not as comprehensive as $\\operatorname{CRPS}\\_{\\text{sum}}$**.\n\n**Table 1. Comparison of $0.5$-risk and $0.9$-risk using GPVar. \"w/o\" denotes methods without time-dependent errors, while \"w/\" indicates our method. Mean and standard deviation are obtained from 10 runs of each model.**\n|                           |               | $0.5$-risk    |                       |                        | $0.9$-risk               |                       |\n|---------------------------|---------------|---------------|-----------------------|------------------------|--------------------------|-----------------------|\n|                           | w/o           | w/            | rel. impr.       | w/o                    | w/                       | rel. impr.       |\n| $\\mathtt{exchange\\\\_rate}$ | 0.0109±0.0003 | 0.0095±0.0004 | 12.84% | 0.0042±0.0001 | 0.0057±0.0001 | -35.71% |\n| $\\mathtt{solar}$          | 0.4998±0.0025 | 0.5246±0.0016 | -4.96% | 0.1617±0.0004 | 0.1597±0.0003 | 1.24%   |\n| $\\mathtt{electricity}$    | 0.0405±0.0003 | 0.0397±0.0002 | 1.98%  | 0.0211±0.0004 | 0.0185±0.0002 | 12.32%  |\n| $\\mathtt{traffic}$        | 0.0933±0.0001 | 0.0859±0.0001 | 7.93%  | 0.0666±0.0001 | 0.0580±0.0001 | 12.91%  |\n| $\\mathtt{wikipedia}$      | 0.2231±0.0005 | 0.2236±0.0006 | -0.22% | 0.2136±0.0002 | 0.2048±0.0001 | 4.12%   |\n| $\\mathtt{m4\\\\_hourly}$     | 0.0807±0.0001 | 0.0849±0.0001 | -5.20% | 0.0452±0.0002 | 0.0463±0.0001 | -2.43%  |\n| $\\mathtt{m1\\\\_quarterly}$  | 0.2196±0.0023 | 0.1948±0.0005 | 11.29% | 0.3049±0.0044 | 0.2787±0.0027 | 8.59%   |\n| $\\mathtt{pems03}$         | 0.0568±0.0000 | 0.0574±0.0001 | -1.06% | 0.0317±0.0000 | 0.0317±0.0001 | 0.00%   |\n| $\\mathtt{uber\\\\_hourly}$   | 0.1035±0.0002 | 0.1013±0.0002 | 2.13%  | 0.0533±0.0002 | 0.0528±0.0001 | 0.94%   |\n\n> line 169-171: could you clarify (also in the final paper) the sentence: Since the parameters of mapping functions are shared for all time series, we can view z_t as a Gaussian process assessed at points hi_t?\n\n***Response:*** We apologize for any confusion caused. In this paper, the time series variables $\\mathbf{z}\\_{t}$ jointly follow a multivariate Gaussian distribution. Since any subset of $\\mathbf{z}\\_{t}$ will also follow a multivariate Gaussian distribution, we can use a training batch comprising $B$ random time series from a total of $N$ time series, on which the likelihood is computed.\n\nFollowing [1], we assume $\\left.\\mathbf{z}\\_{t} \\mid \\mathbf{h}\\_{t}\\right.\\sim\\mathcal{N}\\left( \\boldsymbol{\\mu}(\\mathbf{h}\\_{t}), \\boldsymbol{\\Sigma}(\\mathbf{h}\\_{t})\\right)$, where $\\boldsymbol{\\Sigma}\\_t=\\boldsymbol{L}\\_{t}\\boldsymbol{L}\\_{t}^\\top+\\text{diag}{(\\mathbf{d}\\_{t})}$, and $\\boldsymbol{\\mu}\\_{t}$, $\\boldsymbol{L}\\_{t}$, and $\\mathbf{d}\\_{t}$ are conditioned on $\\mathbf{h}\\_{t}$:\n\n$$\n\\boldsymbol{\\Sigma}\\left(\\mathbf{h}\\_t\\right)=\\left[\\begin{array}{ccc}d_1\\left(\\mathbf{h}\\_{1, t}\\right) & & 0 \\\\\\ & \\ddots & \\\\\\ 0 & & d_N\\left(\\mathbf{h}\\_{N, t}\\right)\\end{array}\\right]+\\left[\\begin{array}{c}\\mathbf{l}_1\\left(\\mathbf{h}\\_{1, t}\\right) \\\\\\ \\cdots \\\\\\ \\mathbf{l}_N\\left(\\mathbf{h}\\_{N, t}\\right)\\end{array}\\right]\\left[\\begin{array}{c}\\mathbf{l}_1\\left(\\mathbf{h}\\_{1, t}\\right) \\\\\\ \\cdots \\\\\\ \\mathbf{l}_N\\left(\\mathbf{h}\\_{N, t}\\right)\\end{array}\\right]^\\top=\\boldsymbol{D}_t+\\boldsymbol{L}\\_{t} \\boldsymbol{L}\\_{t}^\\top.\n$$\n\nwhere $\\mu_i$, $d_i$, and $\\mathbf{l}_{i}$ are the mapping functions used to transform the hidden state $\\mathbf{h}\\_{i,t}$ of time series $i$ at time $t$ to the corresponding distribution parameters. By shared mapping functions, we mean that we use a set of global mapping functions $\\tilde{\\mu}(\\cdot)$, $\\tilde{d}(\\cdot)$, and $\\tilde{\\mathbf{l}}(\\cdot)$ for all time series, instead of having a separate function for each time series. The shared mapping functions enable us to take any subset of time series and calculate the mean and covariance.\n\n[1] Salinas, David, et al. High-dimensional multivariate forecasting with low-rank gaussian copula processes. In NeurIPS, 2019.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TlSUJhedYk", "reviewer": "Reviewer_szpP", "summary": "Briefly summarise the paper and its contributions. This is not the place to critique the paper; the authors should generally agree with a well-written summary.\n\nThe paper proposes an extension to multivariate forecasting model that takes into account the error auto-correlation. The issue of remaining auto-correlation residual is first illustrated clearly in an example with a multivariate forecasting method and the author then introduces their method. The method consists in modelling the cross-covariance through a learned weighted sum of kernels matrices that are projected to obtain the final estimated correlation. This approach leverages previous work on low-rank parametrisation of multivariate forecasting and the authors manages to find an approach that remains tractable despite the complexity of the task. Experiments are conducted on various real-world datasets where their method is shown to improve several forecasting metrics such as CRPS-SUM, RMSE or energy scores.", "questions": "I am separating my questions between the most important ones that would help to raise my score and minor comments. \n\n## Most important\n\n**Lack of clarity whether the method will be competitive with better likelihoods.** Currently, the method only studies a Gaussian noise. However, this noise is a poor fit for the dataset studied and as such is inflating the residuals. This leaves open the question of the method works just because the initial residuals are very poor (and adds a smoothing effect which helps their estimation) and not so much because of the reason given by the paper (that the approach models the autocorrelation aspect).\n\nI can see two ways to remove this potential cofounding factor, one is to include experiments with paranormal distribution which would be more robust from [2], the other is to consider student-t distribution which should also be tractable.\n\n**Details of Hyperparameter Optimization.** The authors mentions:\n```\nThe best hyperparameters for a base model—dataset combination are selected based on the best validation loss. Once the optimal learning rate and hidden size are determined, we apply the same hyperparameters for models both with and without applying our approach \n```\nHow does it work exactly? Which model are you using to determine and fix the hyperparameters? (Selecting the best hyperparameters for your method may favour their transformation)\n\n**Matching budget.** The method proposed by the author is roughly 10x slower on average (and up to 30x slower for exchange-rate from table 10). This point is important and currently a bit buried down, for instance the authors mentions: “* However, there is no evidence that our method slows down convergence. On the contrary, it accelerates convergence for 5 out of 9 datasets for GPVar.” However, the wallclock time are clearly much slower I think this is fine for the paper overall but it should be stated clearly as a limitation (the method proposed by the authors would likely be worse than 10 ensemble of the baseline for instance, I don’t expect this particular point to appear in text or experiment but I want to explain why this is important).\n\n**Citation of multivariate evaluation metrics.** The authors uses CRPS-sum as their main error metrics (in addition to many other metrics including energy scores in their appendix), I would highly recommend citing this paper which discuss the limitation for this metric when assessing multivariate methods (not mine):\n\n*Regions of Reliability in the Evaluation of Multivariate Probabilistic Forecasts*\nÉtienne Marcotte, Valentina Zantedeschi, Alexandre Drouin, Nicolas Chapados\n*Proceedings of the 40th International Conference on Machine Learning*, PMLR 202:23958-24004, 2023.\n\nTo be clear, I do not think this is a large problem given that the authors also conducted many high quality qualitative analysis of their results.\n\n## Details\n- It would be very useful to clearly spell out the dimension/numbers of projection parameters in a table, possibly in the appendix\n- I highly encourage the authors to release their code (the intention is mentioned in the supp, the work will be hard to reproduce w/o given that it is quite technical, and the work will be more impactful as other will be able to directly compare to their setup)\n- a small point regarding \"real-world data frequently exhibit significant error autocorrelation and cross-lag correlation due to factors such as missing covariates\", I agree that missing covariates is indeed a potential use-case but it feels very niche to me (none of the dataset have covariates for instance) compared to model mispecification", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* the paper is well written, even though the approach is very technically involved, the authors did a great job at explaining it well \n* the paper is well motivated: insighful data-analysis are made to illustrate the presence of autocorrelated residual and motivate their approach\n* large coverage of experiments and empirical evaluations: the experiments are done on a large set of datasets covering many metrics and are also well assessed qualitively", "weaknesses": "* potential gaps in experiments:\n  * matching budget: the method is currently much more expensive that the baseline considered\n  * some details of Hyperparameter Optimization are currently missing (but will be added at the camera ready according to the rebuttal)\n* the paper lacks clarity on whether the method will be competitive with better likelihoods. Currently the method may work well only because of non-competitive likelihood is used (Gaussian) as opposed to Student-t or Non-paranormal which would get better residuals", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's encouraging comments and insightful suggestions.\n\n# Questions\n## Most important\n> Lack of clarity whether the method will be competitive with better likelihoods. \n\n***Response:*** Thank you for your valuable feedback. The reviewer is correct in noting that a better likelihood function could potentially regularize the training process, leading to less correlated residuals. For instance, assuming the errors follow a multivariate $t$-distribution can enhance the model's robustness to outliers. Additionally, a more powerful base model can also produce more independent residuals. Given these factors, we designed our approach to dynamically adapt to varying levels of error correlation. Our weighted correlation matrix can assign a higher weight to the identity matrix when errors are less correlated.\n\n**We tested training the baseline models using the likelihood of the multivariate $t$-distribution, and the results are presented in Table 1 (see PDF for the global rebuttal).** It is true that using an alternative distribution can achieve better performance in certain datasets when our method is not implemented. However, we also observed that our method successfully remedies the performance gap for those datasets where the Gaussian assumption is outperformed by the $t$-distribution.\n\n*Note: An important aspect of our method is the ability to use a subset of time series in a training batch for model optimization, which improves scalability. For the multivariate $t$-distribution, the distribution of the subsets of $\\boldsymbol{z}\\_t$ should have the same degrees of freedom as the distribution of $\\boldsymbol{z}\\_t$. Since the degrees of freedom becomes an additional output of the model in each training batch, they are not guaranteed to be the same. While this is not a problem for deep learning, it violates the marginalization rule of the $t$-distribution from a statistical perspective.*\n\nTherefore, we chose Gaussian noise for its numerous beneficial properties, including its marginalization rule and well-defined conditional distribution, which are essential for statistically consistent model training and reliable inference. To address model misspecification, a more effective approach could involve first transforming the original observations $\\mathbf{z}\\_{t}$ into Gaussian-distributed data $\\mathbf{x}\\_{t}$ using a Gaussian Copula [1], and then applying our method.\n\n> Details of Hyperparameter Optimization. \n\n***Response:*** Thank you for bringing this to our attention. When optimizing the hyperparameters of our model, we first tune each base model (e.g., GPVar) on each dataset individually and select the optimal hyperparameters based on the loss observed on the validation set. For example, there will be an optimal learning rate and hidden size for the GPVar model on the $\\mathtt{traffic}$ dataset. We perform this tuning using the baseline models (i.e., without implementing our method). To ensure fair comparison, we use the same set of hyperparameters for the models implemented with our method.\n\n> Matching budget. \n\n***Response:*** The reviewer rightly points out the increased training cost introduced by our method. In Table 10, the clock time per epoch for training with our method is significantly higher compared to the baseline methods. This increase is attributed to the larger $DB \\times DB$ covariance matrix used in the likelihood calculation. We will clearly state this limitation in the next version of the paper.\n\n> Citation of multivariate evaluation metrics. \n\n***Response:*** The reviewer is correct about the limitations of $\\operatorname{CRPS}\\_{\\text{sum}}$ in assessing multivariate probabilistic forecasting, despite its wide usage in existing works [1-5]. $\\operatorname{CRPS}\\_{\\text{sum}}$ has been reported to overlook the performance of the model on each dimension [6]. In [7], the authors suggest that neither $\\operatorname{CRPS}$ nor the energy score is a good surrogate for the NLL score in small-sample regimes, which are common in practice. However, comparing NLL is infeasible in our case, as sampling must be conducted to calibrate predictions autoregressively over the prediction range. Therefore, we include various metrics for assessment in this paper. We will discuss the limitations of $\\operatorname{CRPS}\\_{\\text{sum}}$ in the next version of the paper.\n\n[1] Salinas, David, et al. High-dimensional multivariate forecasting with low-rank gaussian copula processes. In NeurIPS, 2019.\n\n[2] Rasul, Kashif, et al. Multivariate probabilistic time series forecasting via conditioned normalizing flows. In ICLR, 2021.\n\n[3] Rasul, Kashif, et al. Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. In ICML, 2021.\n\n[4] Drouin, Alexandre, et al. Tactis. In ICML, 2022.\n\n[5] Ashok, Arjun, et al. Tactis-2. In ICLR, 2024.\n\n[6] Koochali, Alireza, et al. Random noise vs. state-of-the-art probabilistic forecasting methods: A case study on CRPS-Sum discrimination ability. Applied Sciences 12.10 (2022): 5104.\n\n[7] Marcotte, Étienne, et al. Regions of reliability in the evaluation of multivariate probabilistic forecasts. In ICML, 2023.\n\n## Details\n> It would be very useful ..\n\n***Response:*** Thank you for the suggestion. We have provided more details about the number of projection parameters (see Table 2 and Table 3 in the PDF for the global rebuttal.). It can be seen that we introduce only a negligible amount of new parameters to the model.\n\n> I highly encourage ...\n\n***Response:*** Yes, we will release the code as soon as the work is published.\n\n> a small point regarding ...\n\n***Response:*** We indeed have covariates for each dataset, which include generic features that encode time and time series identification, as detailed in Appendix A.6. We agree that a better likelihood, a stronger base forecasting model, and comprehensive covariates all contribute to less correlated residuals.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
Dear AC and Reviwers,

We would like to express our gratitude for your thorough and insightful reviews. We have carefully considered each of your comments and suggestions. Below, we provide a summary of our responses to the main points raised by the reviewers.

1. Clarifications and Misunderstandings
    - Learning and inference process: We have clarified that our model is trained end-to-end and performs inference in an autoregressive manner.
    - Clarification of viewing $\mathbf{z}_{t}$ as a Gaussian process.
2. Comparison with Related Work
    - Is the method competitive with better likelihood? We have added experiments showing that our method is competitive with a likelihood based on the $t$-distribution (Table 1 in pdf).
3. Experimental Results
    - How does the number of time series influence our method? We have added experiments showing the effect of increasing the number of time series in a batch when performing inference (Figure 1 in pdf).
    - How can the decorrelation of the residuals be measured? We have shown the decorrelation effect in our Appendix using ACF plots and cross-correlation plots.
    - Does the improvement in CRPS come from the point forecast or better variance/covariance? We have provided quantile loss as two snapshots from the predictive distribution to assess the improvement (Table 1 in the response to Reviewer Vi7g). 
    - Experiment details: We have provided information about the number of parameters in each component of our models (Table 2&3 in pdf).

We believe that these revisions have significantly strengthened the paper and addressed the reviewers' concerns. We appreciate the reviewers’ time and effort in providing constructive feedback, which has been invaluable in improving our work.

Thank you for your consideration.

Sincerely,

The authors

Meta-review: The reviews are generally strong and agreement is high. I would strongly
encourage the authors to incorporate the feedback and suggestions into their
revision, as well as their clarifications. Please prioritize the addition of a figure 
showing all components the discussion of error correlations over longer forecast 
horizons. Please also be sure to modify the paper along the lines of the author
rebuttal, especially to reviewer szpP.
Decision: Accept (poster)