End-to-End Modeling of Hierarchical Time Series Using Autoregressive Transformer and Conditional Normalizing Flow-based Reconciliation

The hierarchical time series can be expressed as a tree structure (see Figure 2) with linear aggregation constraints, represented by aggregation matrix $S\in\mathbb{R}^{n\times m}$ ($m$ is number of bottom-level nodes, $n$ is total number of nodes). Each node in the hierarchy represents one time series, to be predicted over a time horizon.

Given a time horizon $t\in\{1,2,...,T\}$, we use $y_{t,i}\in\mathbb{R}$ to denote the values of a multivariate hierarchical time series, where $i\in\{1,2,...,n\}$ is the index of the individual univariate time series. Here we assume that the index $i$ of the individual time series abides by the level-order traversal of the hierarchical tree, going from left to right at each level, and we use $x_{t,i}$ to denote time-varying covariate vectors associated with each univariate time series $i$ at time step $t$.

In our tree hierarchy, the time series of leaf nodes are called the bottom-level series $b_{t}\in\mathbb{R}^{m}$, and those of the remaining nodes are termed the upper-level series $u_{t}\in\mathbb{R}^{r}$. Obviously, the total number of nodes $n=r+m$, and $y_{t}=[u_{t},b_{t}]^{T}\in\mathbb{R}^{n}$ contains observations at time t for all levels , which satisfy (by an aggregation matrix $S\in\{0,1\}^{n\times m}$):

for every time step $t$. For example, in Figure 2, the total number of series in the hierarchy is $n=r+m=3+4$, i.e., the number of series at the bottom-level is $m=4$ and the number of upper-level series is $r=3$. For every time step $t$, $u_{t}=[y_{A,t},y_{B,t},y_{C,t}]\in\mathbb{R}^{3}$ and $b_{t}=[y_{D,t},y_{E,t},y_{F,t},y_{G,t}]\in\mathbb{R}^{4}$, the $y_{A,t}=y_{B,t}+y_{C,t}=y_{D,t}+y_{E,t}+y_{F,t}+y_{G,t}$, and the aggregation matrix $S\in\{0,1\}^{7\times 4}$.

Recently, the encoder-decoder transformer structure has been highly successful in advancing the research on multivariate time series forecasting, enabled by its multi-head self-attention mechanism to capture both long- and short-term dependencies in time series data (Zhou et al., 2021; Li et al., 2019). Meanwhile, combined with classical autoregressive models for time series, we can extend the transformer structure to an autoregressive deep learning model using causal masking, which preserves the autoregressive property by utilizing a mask that reflects the causal direction of the progressing time, i.e., masking out data from future (Katharopoulos et al., 2020).

Specifically, given $\mathcal{D}$, defined as a batch of time series $\bm{Y}=[\bm{y}_{1},...,{\bm{y}_{T}}]\in\mathbb{R}^{T\times D}$, the transformer takes in the above sequence $\bm{Y}$, and then transforms this into dimension $l$ to obtain the query matrix $\bm{Q}$, the key matrix $\bm{K}$ and the value matrix $\bm{V}$ as follows:

where $\bm{Q}\in\mathbb{R}^{dl\times dh}$, $\bm{K}\in\mathbb{R}^{dl\times dh}$, and $\bm{V}\in\mathbb{R}^{dl\times dh}$ (we denote that $dh$ is the length of time steps), and the $\bm{W}_{l}^{Q}$, $\bm{W}_{l}^{K}$, $\bm{V}_{l}^{V}$ are learnable parameters. After these linear transformation, the scaled dot-product attention computes the sequence of vector outputs via:

where the mask matrix $M$ can be applied to filter out right-ward attention (or future information leakage) by setting its upper-triangular elements to $-\infty$ and normalization factor $d_{K}$ is the dimension of the $W_{h}^{K}$ matrix. Finally, all outputs $S$ are concatenated and linearly projected again into the next layer.

Transformer is commonly used in an encoder-decoder architecture, where some warm-up time series as context are passed through the encoder to train the decoder and autoregressively obtain predictions.

In this work, we employ the multivariate autoregressive transformer architecture, to model each time series simultaneously via globally shared parameters. By doing so, we achieve the information fusion of all levels in hierarchy to generate the base forecasts. In fact, in our case, base forecasts from the base model do not directly correspond to un-reconciled forecasts for the base series, but rather represent predictions of an unobserved latent states, which can be used for subsequent density estimations.

Normalizing flows(NF) (Papamakarios et al., 2019), which learn a distribution by transforming the data to samples from a tractable distribution where both sampling and density estimation can be efficient and exact, have been proven to be powerful density approximators. The change of variables formula (Equation 5 below) empowers the computation of exact likelihood (Kobyzev et al., 2020), which is in contrast to other powerful density estimator such as Variational Autoencoders or Generative Adversarial Networks. Impressive estimation results, especially in the field of nonlinear high-dimensional data generation, have lead to great popularity of flow-based deep models (Kingma and Dhariwal, 2018).

NF are invertible neural networks that typically transform isotropic Gaussians to characterize a more complex data distribution. They map from $\mathbb{R}^{D}$ to $\mathbb{R}^{D}$ such that densities $p_{Y}$ on the input space $Y\in\mathbb{R}^{D}$ are transformed into some tractable distribution $p_{Z}$ (e.g., an isotropic Gaussian) on space $Z\in\mathbb{R}^{D}$. This mapping function, $f:Y\rightarrow Z$ and inverse mapping function, $f^{-1}:Z\rightarrow Y$ are composed of a sequence of bijections or invertible functions, and we can express the target distribution densities $p_{Y}(\bm{y})$ by

where $\partial f(y)/\partial y$ is the Jacobian of $f$ at $y$. NF have the property that the inverse $y=f^{-1}(\bm{z})$ is easy to evaluate and computing the Jacobian determinant takes $O(D)$ time.

For mapping functions $f$, the bijection introduced by RealNVP architecture (the coupling layer) (Dinh et al., 2016) satisfies the above properties, leaving part of its inputs unchanged, while transforming the other part via functions of the un-transformed variables (with superscript denoting the coordinate indices)

where $\odot$ is an element wise product, $s()$ is a scaling and $t()$ is a translation function from $\mathbb{R}^{D}\mapsto\mathbb{R}^{D-d}$, using neural networks. To model a nonlinear density map $f(x)$, a number of coupling layers, mapping $\mathcal{Y}\mapsto\mathcal{Y_{1}}\mapsto...\mapsto\mathcal{Y_{K-1}}\mapsto\mathcal{Y_{K}}\mapsto\mathcal{Z}$, are composed together with unchanged dimensions.

A schematic overview of our hierarchical end-to-end architecture can be found in Figure 3.

During training, our loss function is directly computed on the coherent forecast samples. Specifically, given $\mathcal{D}$, defined as a batch of time series $Y:=\{y_{1},y_{2},...,y_{T}\}$, and the associated covariates $X:=\{x_{1},x_{2},...,x_{T}\}$, we can maximize the likelihood given by Equation 9 via SGD using Adam, i.e.

where the globally shared parameters ($\theta$, $\phi$) and $\psi$ are from the transformer module and the CNF module, respectively.

Please note that we can easily obtain $\theta$, $\phi$, and $\psi$ for other loss functions such as quantile loss, CRPS (continuous ranked probability score) or any other metrics preferred in the forecasting community, as long as we can compute the sufficient statistics from the Monte Carlo samples $\{\tilde{y_{t}}\}$ via the empirical distribution function.

During inferencing, we can predict using the autoregressive transformer in a step-by-step fashion over the time horizon. Specifically, we first generate the base forecasts $\hat{y}_{t}$ using autoregressive transformer for one time step using the covariate vector $x_{1:t-1}$ and the observed value $y_{1:t-1}$. Then we can incorporate base forecasts $\hat{\bm{y}}_{t}$ from all levels into CNF as additional condition $\bm{h}_{t}$ in the latent space to model conditional joint distribution $p_{Y}(\tilde{\bm{y}}_{t}|\bm{h}_{t})$. Lastly, we directly obtain a set of Monte Carlo samples from the distribution $p_{Y}(\tilde{\bm{y}}_{t}|\bm{h}_{t})$ as the reconciled bottom forecasts, which are multiplied by the aggregation matrix $S$ to obtain the coherent probabilistic forecasts $\tilde{\bm{y}}_{t}$.

Note that we can repeat the above procedure for the $h$-period-ahead forecasting to obtain a set of coherent forecasts $\{\tilde{y_{T}},\tilde{y_{T+1}},...,\tilde{y_{T+h}}\}$.

We conduct the performance comparison against state-of-the-art reconciliation algorithms below:

• •

  ARIMA-NaiveBU and ETS-NaiveBU: NaiveBU generates univariate point forecasts for the bottom-level time series independently and then sums them according to the hierarchical constraint to obtain point forecasts for the aggregate series. Specifically, we use ARIMA and ETS with auto-tuning for base forecast in NaiveBU with the R package hts (Hyndman and Wickramasuriya, 2020).

• •

  ARIMA-MinT-shr, ARIMA-MinT-ol, ETS-MinT-shr, ETS-MinT-ols : For MinT, we use the covariance matrix with shrinkage operator (MinT-shr) and the diagonal covariance matrix corresponding to ordinary least squares weights (MinT-ols) (Wickramasuriya et al., 2019). In addition, we also use ARIMA and ETS for auto-tuning of base forecasts in Mint with the R package hts (Hyndman and Wickramasuriya, 2020).

• •

  ARIMA-ERM and ETS-ERM: ERM is from Regularized Regression for Hierarchical Forecasting (Ben Taieb and Koo, 2019) in related work discussed in the previous section and we also use ARIMA and ETS for auto-tuning of base forecasts in ERM with the R package hts (Hyndman and Wickramasuriya, 2020).

• •

  N-BEATS-NaiveBU: N-BEATS (Oreshkin et al., 2019) is an interpretable time Series forecasting model, which is a deep neural architecture based on backward and forward residual links and a very deep stack of fully-connected layers. Also We use the above NaiveBU method to ensure coherency. The source code is available at https://github.com/philipperemy/n-beats.

• •

  Transformer-NaiveBU: Transformer (Vaswani et al., 2017) is a recently popular method modeling sequences based on self-attention mechanism. We use the native transformer implementation from pytroch. We use the above NaiveBU method to ensure coherency.

• •

  Informer-NaiveBU: Informer (Zhou et al., 2021) is a probsparse self-attention mechanism based model to enhance the prediction capacity in the Long sequence time-series forecasting problem, which validates the Transformer-like model’s potential value to capture individual long-range dependency between long sequence time-series outputs and inputs. Also We use the above NaiveBU method to ensure coherency. The source code is available at https://github. com/zhouhaoyi/Informer2020.

• •

  DeepAR-NaiveBU: DeepAR (Salinas et al., 2020)is a univariate probabilistic forecasting model based rnn, which is widely applied in real-wrold industrial application. Also We use the above NaiveBU method to ensure coherency. The code is from the Gluonts implementatin on https://github.com/awslabs/gluon-ts.

• •

  DeepVAR-lowrank-Copula: DeepVAR is a multivariate, nonlinear generalization of classical autoregressive model based rnn (Salinas et al., 2019). We use the vanilla model with no reconciliation, parameterizing a low-rank plus diagonal convariance via Copula process. We set the rank of low-rank convariance to be 5. The code is from the Gluonts implementatin on https://github.com/awslabs/gluon-ts.

• •

  Hier-E2E: In Hier-E2E (Rangapuram et al., 2021), we use the DeepVAR(Salinas et al., 2019) with the diagonal covariance matrix to obtain the base forecast, and use pre-calculated reconciliation matrix $M$ to project into the coherent subspace. The code is from the authors’ implementation on https://github.com/rshyamsundar/gluonts-hierarchical-ICML-2021.

All experiments run on the Linux server(Ubuntu 16.04) with the Intel(R) Xeon(R) Silver 4214 2.20GHz CPU, 16GB memory, and the signle Nvidia V-100 GPU.

For evaluation, we generate 200 samples from our probabilistic model to create an empirical predictive distribution. We run our method 5 times and report the mean and standard deviation of CRPS scores. Table 1 shows the performances of different approaches on the four hierarchical datasets, and we can observe that our proposed approach (Hier-Transformer-CNF) achieves all best results, with significant improvements of accuracy on most datasets.

Furthermore the results of experiments proven that effect of traditional classical statistical methods is far inferior to the deep learning model based transformer or rnn. By harnessing the power of deep model, We not only achieve best results but also overcome the shortcomings of statistical models, which rely on any assumption such as unbiased estimates or Gaussian distribution.

To further demonstrate the effectiveness of the designs described in Section 4, we conduct ablation studies using the variant of our model and the related models. The results of each aggregation level are shown in Figure 7.

DeepVAR-lowrank-copula is an RNN-based model using the lowrank Gaussian coupla (Salinas et al., 2019) and Autoregressive Transformer is based the encoder-decoder transformer architecture (Katharopoulos et al., 2020), enabled by its multi-head self-attention mechanism to capture both long- and short-term dependencies in time series data. Hier-E2E is the end-to-end model combing DeepVAR(Gaussian distribution) and projection matrix (Rangapuram et al., 2021)

It is obvious that transformer-based models outperform RNN-based models, especially in the upper-level data. The transformer allows the access to any part of the historic data regardless of temporal distance and is thus capable of generating better conditioning for NF head, evidenced by our experiments.

The non-Gaussian data distribution and nonlinear correlations depicted in Figure 7 attest the need for more expressive density estimators, i.e., CNF, which makes our approach significantly outperforming other methods especially in the bottom-level data. Moreover, the end-to-end reconciliation method also beats direct prediction model in the forecasting at all levels.

In summary, significant improvements at all levels prove the the effectiveness of the our approach.