Coherent Probabilistic Aggregate Queries on Long-horizon Forecasts

1 Introduction

Long-term forecasting is critical for decision support in several application domains including finance, governance, and climate studies. Interactive decision making tools need accurate probability distributions not just over the long-term base forecasts but also on dynamically chosen aggregates of the forecasts. An analyst may be interested in a variety of aggregates such as average, sum, or trend of values for various windows of granularity. For example, the base-level forecasts may be retail price of a commodity at a daily-level for the next two quarters, but the end user may also inspect average at weekly or monthly granularity, or change in price (inflation) from week to week. For the aggregates too, the user is interested in the associated uncertainty represented as a distribution that is coherent with the base distribution.

A time-series in our model is characterized by a sequence (x1, y1) . . . (xn, yn) where xt ∈ Rd denotes a vector of input features and yt is a real-value at a time t. We are given a prefix or history from such a series up to a present time T which we denote as HT . Our goal is to predict future values yT+1 . . . yT+R given input features x = xT+1 . . . xT+R and the history HT . Unlike in conventional time-series forecasting where the interest is short-term (say R between 1 and 10), in long-term forecasting we focus on large values of R (say between 100 and 1000). We wish to capture the uncertainty of the prediction by outputting a distribution over the predicted values. We use ˆ P(yT+1 . . . yT+R|HT , x) to denote the predicted distribution for future time horizon given history HT .

∗Contact author

Our goal is not just to provide accurate probabilistic forecasts for the the base-level values but also for aggregates of these values. For example, the base-level forecasts may be at the daily-level and we may wish to view monthly sales. In theory, we could marginalize the base-level forecast distribution ˆ P(yT+1 . . . yT+R|HT , x) to compute coherent distribution of any dynamic aggregate. In practice, we need to ensure that the marginalization can be efficiently computed. Furthermore, the accuracy of the predicted aggregate should be no worse than the accuracy of an independently trained model to directly predict the distribution of the aggregate. For example, a practitioner could train separate models for predicting monthly sales and daily sales for greater accuracy at the monthly-level. However, training separate models for each aggregate has two limitations: (1) independently predicted aggregate distributions may not be coherent with each other or with the base-level forecasts, and (2) the set of aggregate functions may be decided dynamically during data analysis where training a separate model is impractical. Thus, we seek solutions that explicitly marginalize a fixed distribution to ensure coherence of distributions of aggregate quantities.

Existing forecast models can be categorized as either autoregressive (AR) [Flunkert et al., 2017; Mukherjee et al., 2018] or non-auto regressive (NAR). We will show in this paper that both these models are inadequate for our goal of efficiently and dynamically defining the distribution of any aggregate quantity. In AR models the output distribution captures the full joint dependency among all predicted values by defining that the value at t depends on all values before it. During training we condition on true previous values but during inference we need to perform forward sampling where each predicted value is conditioned on sampled prior values. In long-term forecasting, such forward sampling suffers from drift caused by cascading errors [Deshpande et al., 2021]. A second major limitation of the AR model is that the joint distribution over the forecast variables is expressed via a blackbox neural network. Marginalizing the distribution to obtain distribution over dynamically chosen aggregates will involve repeated sampling steps, which may incur huge latency during interactive analysis.

This has led to greater interest in the non-auto-regressive (NAR) models which independently predict each variable in the forecast horizon conditioned on the history. NAR models have been shown to work better in practice both for medium-term and long-term forecasting [Wen et al., 2017; Deshpande and Sarawagi, 2019; Zhou et al., 2020]. Another major benefit of the NAR models is that during inference all values can be predicted in parallel, unlike in the AR method that entails sequential sampling. However, a limitation of the non-auto-regressive methods is that the dependency among the output variables is not modeled. The distribution over aggregates computed from independent base predictions may be inaccurate. For example, if an analyst marginalizes the data to inspect the distribution of the difference between adjacent values, the variance from the independent model might be over-estimated. Further, in the raw form a time-series often contains noise that hinders a forecasting model from capturing top-level patterns in the series over long prediction horizons. In applications that involve analysis over aggregates of forecast values, these might give rise to non-smooth or inconsistent aggregates.

In this paper we propose a novel method of preserving higher-level patterns in long-range forecasts of non-autoregressive forecast models while allowing accurate and efficient computation of the distribution of dynamically chosen aggregates. Our key idea is to train the model to independently output various high-order statistics over the forecast values, along with the base-level forecasts. Aggregated values provide a noise-removed signal of how the series evolves in time. We use such signals to guide inference to capture top-level trends in the data. We formulate that as a task of obtaining a revised consensus distribution that minimizes the KL-distance with each predicted aggregate and base forecasts. The consensus distribution is represented as a Gaussian distribution that can be easily marginalized. The parameters of the consensus can be obtained via an efficient numerical optimization problem. Unlike existing methods [Han et al., 2021], our consensus distribution not only provides coherence and accuracy for aggregate quantities, but also improves base-level forecasts.

1.1 Contributions

• We propose a probabilistic forecasting method that is more accurate than state of the art for long-term forecasting.
• Our method outputs an easy to marginalize joint distribution for computing accurate and coherent aggregate distributions over dynamically defined functions on forecast values.
• We propose an efficient inference algorithm for computing the consensus distribution.
• Empirical evaluation on benchmark datasets show that our method improves both base-level and aggregate forecasts compared to several existing baselines including a recent hierarchical consensus model.

2 Our Model

Our method seeks to retain the advantage of parallel training and inference that modern transformer-based NAR models offer while providing accurate distributions of forecasts at the base-level and dynamically chosen aggregates over them. We train multiple NAR models for forecasting different types of aggregate functions at different levels of granularity. During inference we obtain a consensus distribution, as a joint Gaussian distribution that can be marginalized in closed form for computing distributions over linear aggregates defined dynamically during analysis. The consensus distribution is obtained by minimizing the KL distance with the various base and aggregate-level forecasts by an efficient inference algorithm. We present the details of each of these steps next.

2.1 Aggregate Functions

We choose a set A of aggregate functions depending on likely use in downstream analysis. This set does not need to include all possible aggregate functions to be used during analysis, since our method generalizes well to new aggregates even when trained with two aggregate functions. Each aggregate function Ai is characterized with a window size Ki and a fixed real vector ai ∈ RKi that denotes the weight of aggregated values. For each Ai, we create an aggregated series by aggregating on disjoint windows of size Ki in the original series. The j-th value in the ith aggregated series is the aggregation of the values y(j-1)Ki+1 . . . yjKi and denoted zij . We will use the notation wi,j to refer to the window of indices [(j - 1)Ki, . . . , jKi] that the j-th value of the ith function aggregates. We denote i-th aggregate series as zi1 , . . . , z iTi where Ti = T Ki and calculate as follows:

Here are examples of two aggregations - (i) average and (ii) trend. Average is just:

Trend aggregate of the window captures how the series changes along time. A positive (negative) trend denotes that the value increases (decreases) as we move forward in the window. We compute trend as the slope of a linear fit on the values in each window as

We show examples of a base-level time series and one computed aggregate (average) in Figure 1. Our framework can work with any aggregate like the above two that can be expressed as a linear weighted sum of its arguments. For notational ease we refer to the base-level also as an aggregate function A0 with window size Ki = 1 and ai = [1].

2.2 Forecast Method

We train independent probabilistic forecasting models for the series corresponding to each aggregate Ai ∈ A which includes the original time-series as well. Our work is agnostic of the exact neural architecture used for the forecasts. We present details of how to train these models in Section 2.3.

During inference we use the known history HT and the trained models to get forecasts ˆ P(zij |HT , xj) for each of the variables zi1 , . . . , z iTi for each aggregate function Ai, which includes the base-level quantities.

Since models for all aggregates and for the original series are trained independently, these predicted distributions are not necessarily coherent. We infer a new consensus distribution Q(yT+1 . . . yT+R) that minimizes the KLdistance with each of the forecast distributions. We choose a tractable form for Q: a multivariate Gaussian with mean µ = [µT+1 . . . µT+R]T and covariance Σ. With this form we can compute the marginal distribution of aggregate variable zij as:

We use the notation Σwi,j to denote the sub-matrix of Σ that spans over the indices in wi,j . Likewise for µwi,j .

Using these we can write our objective of finding the parameters µ and Σ so as to minimize the KL-distance with each forecast distribution as:

Here, the weights αi are hyper-parameters which denote importance of KL-distance of i-th aggregate.

When the distribution of each forecast variable ˆ P (zij |·) is represented as a Gaussian with mean ˆµ(zij ) and variance ˆσ(zij ), their KL distance can be expressed in closed form 1 . Also, the optimization over the µ and Σ terms separate out into two independent optimization problems. For µ the objective is now:

This is a convex quadratic objective that can be solved in closed form as

where we use ˜aij ∈ RR to denote the padded version of the indices wi,j .

For Σ the objective is:

Algorithm 1 KLST Training and Inference Algorithm

• 1: Input: Training data {(x1, y1) . . . (xn, yn)}, aggregate functions A = {(ai, Ki)}, forecast horizon R.
• 2: for Ai in A do
• 3: /*Train i-th aggregate model using
• training data for i-th aggregate */
• 4: Ditrn = zi1 , . . . , z in/Ki .
• 5: θi ← TRAIN(Ditrn ) using Eqn 11
• 6: Get forecasts ˆ P(zij |θi ) for j = Ti + 1, . . . , Ti +Ri.
• 7: end for
• 8: Get coherent mean forecasts µ∗ using Eqn. 7.
• 9: Get coherent variance forecasts (V ∗ , σ∗ ) using Eqn. 9.
• 10: /* End of Training */
• 11: Inference: Given new aggregate Anew defined via ˜anew
• 12: µ∗new = µ∗T ˜anew .
• 13: Σ∗new = ∑Rj=1 (˜anew j )2σ∗ j 2 + (˜anew )T V ∗V ∗T (˜anew )

This objective cannot be solved in closed form but we follow a strategy to efficiently approximate it. We consider a simplification of Σ to be a low-rank matrix. Following [Salinas et al., 2019] we express the covariance Σ in a tractable format as follows:

where vT+r ∈ Rk is r-th row of V , and k is the chosen rank of the covariance matrix.

Using the above form of Σ, Equation 8 can be rewritten as:

With the above representation of the joint distribution Q(yT , . . . yT+R), the number of parameters needed to store the distribution is linear in R, and also we can compute the distribution of any linear aggregate efficiently. Algorithm 1 presents an overview of our method.

2.3 Training Parameters of Forecast Models

As mentioned earlier we train separate models for each aggregate series. For each aggregate Ai the forecast model can be viewed as a function F(HT , xR|θi ) of the history of known values HT = (x1, y1) . . . , (xT , yT ), and known input features xT+1 . . . xT+R of the next R values for which the model needs to forecast. We use θi to denote the parameters of the ith forecast model. The exact form of the forecast model F() that we used is described next, followed by a description of the training procedure.

Architecture of Forecast Model

Our forecasting model F (HT , xR|θi ) is based on Transformers like Informer [Zhou et al., 2020] but with subtle differences that provide gains even for the baseline predictions.

Our model first applies a convolution on the input z series and features x. The Transformer's multi-headed attention layers are on a concatenation of the convolution layer output.

where h ∈ RT×dmodel is the output of the transformer encoder and et denotes the positional encoding of time index t. In the decoder, we warm start to provide more context to the decoder apart from the encoder state h. In warm start, we select the last s values in the input series and use them in the decoder of the transformer.

where 0R is a vector of zeros of length R, hµ ∈ RR×dmodel , hσ ∈ RR×dmodel , µ ∈ RR×1 , and σ ∈ RR×1 .

Training Procedure

We start with base level series (x1, y1) . . . (xn, yn). First we split the series into training, validation, and test sets of lengths ltrn, lval, and ltest respectively as follows:

Note that when lval ≥ R or ltest ≥ R, we use rollingwindow setting as done in [Deshpande and Sarawagi, 2019; Rangapuram et al., 2018] on validation and test sets.

We create training batches by sampling chunks of length T + R from training set Dtrn as

Similarly, we process i-th aggregate series zi1 , . . . , z in/Ki . Validation and test sets across base-level and aggregated series are aligned. To aggregate features, we simply take average value of each feature in the aggregation window for all aggregates. For an aggregate model for i-th aggregate with window size Ki, the forecast horizon is Ri = R/Ki. Similarly, the size of the history also changes. We set Ti = B(T/Ki), so the aggregate model can effectively see the history B-times the size of base level model's history. We found setting B=2 works well in practice. It is also a good option for smaller datasets such as ETTH.

For each aggregate function Ai, we sample training chunks of size Ti +Ri and train the parameters θi of the Transformer model for the following training objective.

3 RelatedWork

Time series forecasting is an extensively researched problem. In deep learning, the problem got renewed interest [Benidis et al., 2020] after an RNN-based auto-regressive model [Flunkert et al., 2017; Mukherjee et al., 2018] established conclusive gains over conventional machine learning and statistical methods. While these early neural models were based on RNNs, more recently Transformer-based models have been found to provide faster and efficient predictions on time-series data [Li et al., 2019; Lim et al., 2020; Bansal et al., 2021; Zhou et al., 2020]. For long-term forecasting auto-regressive models are both slow and subject to cascading errors from previous predicted values for next-step predictions. Another option is non-auto-regressive models (NAR) which predict each future forecast independently in parallel [Wen et al., 2017; Deshpande and Sarawagi, 2019; Zhou et al., 2020]. However, since NAR models do not capture the joint distribution, they fail to provide accurate results for aggregate queries on forecasts. N-beats [Oreshkin et al., 2019] is another recent architecture for time-series forecasting based on layered residual connections, specifically designed for interpretability. Our technique is orthogonal to the underlying time-series forecasting model.

Multi-level Models Our approach of independently predicting models at different aggregation levels and then reconciling is related to the approach in [Taieb et al., 2017; Ben Taieb and Koo, 2019; Wickramasuriya et al., 2019]. Their hierarchies are over items and consider only sum-based aggregates whereas we create hierarchies along time and consider arbitrary linear aggregates. Also, our objective is to minimize the KL distance with the distribution at various levels whereas they focus on error minimization and their parameterization of the reconciliation is different from ours. [Han et al., 2021] is another multi-level model for item hierarchies that train models in a bottom up manner with higher level models regularized to be consistent with lower layers. They do not model inter-item correlation, can only efficiently handle fixed aggregates, and do not revise base-level forecasts. We obtain gains at the base-level too via the hierarchies as we show in our empirical comparison. [Rangapuram et al., 2021] propose an end-to-end method that produces coherent forecasts. However, their method relies on the sampling to obtain covariance of the joint distribution whereas we train the covariance through our inference method alone.

Beyond Base-level Forecasting Loss Our method of learning the slope as an aggregate in the forecast horizon is related to recent efforts at trying to learn the shape of the series. Recently, [LE GUEN and THOME, 2019; Guen and Thome, 2020] learn the shape of the series in forecast horizon using a differentiable approximation of Dynamic Time Warping (DTW) loss function. They also contain additional losses that minimize distortion and increase diversity. SRU

captures moving statistics of the history through a recurrent layer [Oliva et al., 2017]. However, they use the statistics only to capture dependence among observed variables.

Long-term Forecasting Model Recently [Zhou et al., 2020] proposed a Transformer architecture for long-range time-series forecasting. Their focus is to increase efficiency of the self-attention in Transformers for long series by enforcing sparsity. To the best of our knowledge, this is the best performing of existing methods for long-term forecasting. We compare empirically with this model and present significant gains on multiple datasets.

4 Experiments

We evaluate our method on four real-life datasets and contrast with state-of-the-art methods of long-term forecasting.

4.1 Datasets

Electricity This dataset contains national electricity load at Panama power system 2 . The dataset is collected at an hourly granularity from January 2015 to March 2020.

Solar This dataset is on hourly photo-voltaic production of 137 stations and was used in [Salinas et al., 2019].

ETT (Hourly and 15-minute) This dataset contain a timeseries of oil-temperature at an electrical transformer 3 collected from July 2016 to June 2017. The dataset is available at two granularities - 15 minutes (ETT) and one hour (ETTH). Both were used in [Zhou et al., 2020].

4.2 Base Models

Informer [Zhou et al., 2020] is a recent Transformer-based model specifically designed for long-term forecasting. They also use a non-auto-regressive model and predict each future value independently.

Trans-NAR refers to our Transformer architecture described in Section 2.3 but where we predict only a single model on only the base-level forecasts.

2 https://data.mendeley.com/datasets/byx7sztj59/1

3 https://github.com/zhouhaoyi/ETDataset

Trans-AR Trans-AR is an autoregressive version of the transformer architecture described in Section 2.3. At the t-th step during prediction, it performs convolution on the window ˆ zt-wy+1, . . . , ˆ zt-1 of predicted values. When t < wf , the convolution window overlaps with the history HT and observed values fromHT is used.

SHARQ SHARQ [Han et al., 2021] is a hierarchical forecasting method in which a regularization loss ensures that the coherency error between aggregate predictions and aggregation of base-level predictions is minimized. Since SHARQ is agnostic to the underlying forecasting model, we use our baseline Trans-NAR as the base model for SHARQ. For j-th value of an aggregate model, we write the regularization term of SHARQ as:

where ˆµ and ˆ qτ denote the predicted mean and aggregate model at j-th position. This model requires a fresh model to be trained for each aggregate unlike our method.

KLST is our method as described in Algorithm 1. The choice of aggregate functions and range of K values are based on validation set and chosen from aggregate functions = { Sum, Slope} and K={6,12}. For all datasets, we assign weights α={10,0.5} to Sum and Slope aggregates respectively. Since direct predictions on Sum aggregate are often more accurate, we assign higher weight to it.

4.3 EvaluationMetric

We evaluate probabilistic forecasts using the wellknown [Salinas et al., 2019] Continuous Ranked Probability Score (CRPS) metric. If we denote the predicted CDF for a target value yt by Ft(·), we can write the CRPS by integrating the quantile loss over all possible quantiles [Salinas et al., 2019]:

It reduces to Mean Absolute Error (MAE) when the forecast is a point,and is a proper scoring metric of the quality of an entire distribution including its mean and variance.

4.4 Qualitative Results

We present some anecdotes to provide insights on why KLST should enhance the accuracy of base-level forecasts. In Figure 1 (left column) we show that forecasts by KLST (blue) are closer to the ground truth (black) compared to the base-level model Trans-NAR (red). We explain these gains by showing the improved accuracy of the model that independently predicts the sum aggregate for a window of K = 6 in the right column of Figure 1. We observe that direct forecast values are much closer to the ground-truth than those from aggregating base-level forecasts. When these direct aggregates are used to revise the base-level forecasts to be consistent with these aggregates (as described in Section 2.2), we obtain better forecasts at the base-level.

4.5 Comparison of KLST with Baselines

Impact of KLST on Long-TermForecasts In Table 1, we compare KLST with the baselines described in Sec. 4.2. Under Sum aggregate and column K=1 (highlighted), we compare the CRPS of base level forecasts for Informer, Trans-AR, Trans-NAR, SHARQ and KLST. For base-level forecasts, KLST performs better than all baselines while Trans-NAR and SHARQ are share the second position. Since SHARQ does not update base level model parameters, base level forecasts of SHARQ are identical to Trans-NAR.

Results on Unseen Aggregates Post Inference An important aspect of time-series analytics is the effectiveness of base-level forecasts on various types of aggregates at various granularities. Training an aggregate model for each aggregate and at each granularity is rather expensive. Our approach is to train a few representative aggregate models using a fixed set of Ki values. Then, forecasts after consensus can be used to compute other aggregate and at other granularities.

In Table 1 (columns with K=4,8,12,24), we compare KLST with (i) all baselines for Sum aggregate and (ii) with Trans-NAR for Slope and Diff aggregates. The set of K values we use in Table 1 also contain the K for which no aggregate model is trained. Also, no aggregate model was trained corresponding to the Diff aggregate. For each aggregate and for each K, first we aggregate the forecasts of baselines and KLST and then evaluate the aggregates using CRPS. Since SHARQ does not update the bottom level forecasts, we have trained a separate aggregate model for each K and Sum aggregate using SHARQ regularization loss (Eqn. 12).

For Sum aggregate, KLST outperforms all baselines on ETT and ETTH dataset. SHARQ performs better than KLST on Solar and Electricity for K=24 and K=4,8 respectively. Since SHARQ requires training a separate model for each aggregate and K, we compare SHARQ only on Sum aggregate.

For Slope and Diff aggregates, KLST outperforms TransNAR on Solar and Electricity. Whereas for ETT and ETTH, KLST is better than Trans-NAR in 12 out of 20 cases.

5 Conclusion

In this paper we addressed the problem of long-range forecasting with the goal of providing accurate probabilistic predictions, not just at the base-level but also at dynamically chosen aggregates. We proposed a simple method of independently training forecast models at different aggregation levels, and then designing a consensus distribution to minimize distance with forecasts from each component distribution. We obtained significant gains in accuracy over the base-level and over new aggregate functions at new granularity levels defined during test time. Compared to existing hierarchical methods, we impact both base-level predictions and new aggregates defined dynamically.