Reconciling Hierarchical Forecasts via Bayes' Rule

TL;DR

This paper introduces a Bayesian approach for hierarchical forecast reconciliation, deriving closed-form solutions under Gaussian assumptions and comparing with existing methods like MinT and Kalman filter.

Contribution

It proposes a novel Bayesian framework for hierarchical forecast reconciliation, providing new algorithms and theoretical insights.

Findings

Closed-form Bayesian updating formulas under Gaussian assumptions

Algorithms that differ based on independence assumptions

Experimental comparison with MinT and Kalman filter methods

Abstract

We present a novel approach for reconciling hierarchical forecasts, based on Bayes rule. We define a prior distribution for the bottom time series of the hierarchy, based on the bottom base forecasts. Then we update their distribution via Bayes rule, based on the base forecasts for the upper time series. Under the Gaussian assumption, we derive the updating in closed-form. We derive two algorithms, which differ as for the assumed independencies. We discuss their relation with the MinT reconciliation algorithm and with the Kalman filter, and we compare them experimentally.

Tables

Table 1. Table 1 : Mean energy score, averaged over 1000 simulations and over h ℎ h =1,2,3,4. In each row we highlight the lower result.

$T$		method			BU		pMinT		LG
50		arima			9.7		9.5		9.4
50		ets			9.9		9.8		9.7
100		arima			9.1		9.0		9.0
100		ets			9.5		9.4		9.3
1000		arima			8.8		8.7		8.8
1000		ets			9.5		9.3		9.4

Table 2. Table 2 : Averaged energy scores. Each cell is the average over 200 reconciliations (50 different training/test splits × \times h ℎ h =1,2,3,4) . The rows corresponding to the best-performing values of k h subscript 𝑘 ℎ k_{h} are highlighted.

dset			$k_h$			method			BU			pMinT			LG
infantgts			1			arima			334.1			346.9			348.5
			1			ets			334.0			320.0			334.7
			$h$			arima			327.2			335.1			331.0
			$h$			ets			328.2			313.7			318.7
tourism			$1$			arima			2,737.6			2,412.0			2,547.4
			$1$			ets			2,496.0			2,403.7			2,520.1
			$h$			arima			2,785.3			2,380.3			2,448.2
			$h$			ets			2,527.1			2,353.6			2,410.3

Table 3. Table 3 : Posterior probability of the Bayesian signed rank test.

					Posterior probabilities
dset	$k_h$	method		pMinT <BU		pMinT <LG		LG <BU
infantgts	1	arima		0.24		0.02		0.47
		ets		1.00		0.85		1
	h	arima		0.89		0.00		1.00
		ets		1.00		0.00		1
tourism	1	arima		1.00		1.00		1.00
		ets		1.00		1.00		0.25
	h	arima		1.00		1.00		1.00
		ets		1.00		1.00		1.00

Equations

$${\mathbf{y}}_t=[y_{\text{Total}},y_{R_1},y_{R_2},y_{R_{11}},y_{R_{12}},y_{R_{21}},y_{R_{22}}{]}^T=[{\mathbf{u}}_t^T,{\mathbf{b}}_t^T{]}^T,$$

$${\mathbf{u}}_t=[y_{\text{Total}},y_{R_1},y_{R_2}{]}^T$$

$${\mathbf{b}}_t=[y_{R_{11}},y_{R_{12}},y_{R_{21}},y_{R_{22}}{]}^T.$$

$${\mathbf{y}}_t=\mathbf{S}{\mathbf{b}}_t.$$

$$\mathbf{S}=\left[\begin{array}{ccccccc}1& & 1& & 1& & 1\\ 1& & 1& & 0& & 0\\ 0& & 0& & 1& & 1\\ 1& & 0& & 0& & 0\\ 0& & 1& & 0& & 0\\ 0& & 0& & 1& & 0\\ 0& & 0& & 0& & 1\end{array}\right]=\left[\begin{array}{ccccccc}& & & & & & \\ & & & \text{makebox}(0.0,0.0)[]\text{A}& & & \\ & & & & & & \\ \text{cdashline}2-6& & & & & & \\ & & & \text{makebox}(0.0,0.0)[]\text{I}& & & \\ & & & & & & \\ & & & & & & \end{array}\right],$$

$${\tilde{\mathbf{b}}}_{t+h}={\mathbf{P}}_h{\hat{\mathbf{y}}}_{t+h}.$$

$${\mathbf{P}}_h=({\mathbf{S}}^T{\mathbf{W}}_h^{-1}\mathbf{S}{)}^{-1}{\mathbf{S}}^T{\mathbf{W}}_h^{-1}$$

$$p({\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b})$$

$$p({\mathbf{Y}}_{t+h}\mid {\mathcal{I}}_{t,b})$$

$${\hat{\mathbf{U}}}_{t+h}=\mathbf{A}{\mathbf{B}}_{t+h}+{\mathbf{\epsilon }}_{t+h}^u,$$

$${\mathbf{\epsilon }}_{t+h}^u\sim N\left(0,{\hat{\mathbf{\Sigma }}}_{U,h}\right),$$

$$p({\hat{\mathbf{U}}}_{t+h}\mid {\mathbf{B}}_{t+h})$$

$$p({\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{U}}}_{t+h})$$

$$=\frac{p({\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b})p({\hat{\mathbf{U}}}_{t+h}|{\mathbf{B}}_{t+h})}{p({\hat{\mathbf{U}}}_{t+h}\mid {\mathcal{I}}_{t,b})}\propto$$

$$\propto p({\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b})p({\hat{\mathbf{U}}}_{t+h}\mid {\mathbf{B}}_{t+h})=$$

$$=p({\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b})p(\mathbf{A}{\mathbf{B}}_{t+h}+{\mathbf{\epsilon }}_{t+h}^u\mid {\mathbf{B}}_{t+h})$$

$$\mathbb{E}[{\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b}]={\hat{\mathbf{b}}}_{t+h},$$

$$\mathbb{E}[{\hat{\mathbf{U}}}_{t+h}\mid {\mathcal{I}}_{t,b}]=\mathbf{A}{\hat{\mathbf{b}}}_{t+h}.$$

$$Cov({\mathbf{B}}_{\mathrm{t}+\mathrm{h}},{\hat{\mathbf{U}}}_{\mathrm{t}+\mathrm{h}}|{\mathcal{I}}_{\mathrm{t},\mathrm{b}})$$

$$=Cov({\mathbf{B}}_{\mathrm{t}+\mathrm{h}},{\mathbf{B}}_{\mathrm{t}+\mathrm{h}}|{\mathcal{I}}_{\mathrm{t},\mathrm{b}}){\mathbf{A}}^T+Cov({\mathbf{B}}_{\mathrm{t}+\mathrm{h}},{\mathbf{\epsilon }}_{\mathrm{t}+\mathrm{h}}^{\mathrm{u}}|{\mathcal{I}}_{\mathrm{t},\mathrm{b}})$$

$$=k_h({\hat{\mathbf{\Sigma }}}_{B,1}{\mathbf{A}}^T+{\mathbf{M}}_1)\in {\mathbb{R}}^{n\times (m-n)}$$

$$\mathrm{Cov}({\hat{\mathbf{U}}}_{t+h}\mid {\mathcal{I}}_{t,b})$$

$$=k_h\mathbf{A}{\hat{\mathbf{\Sigma }}}_{B,1}{\mathbf{A}}^T+k_h{\hat{\mathbf{\Sigma }}}_{U,1}+Cov({\mathbf{AB}}_{\mathrm{t}+\mathrm{h}},{\mathbf{\epsilon }}_{\mathrm{t}+\mathrm{h}}^{\mathrm{u}})+Cov({\mathbf{\epsilon }}_{\mathrm{t}+\mathrm{h}}^{\mathrm{u}},{\mathbf{AB}}_{\mathrm{t}+\mathrm{h}})$$

$$=k_h(\mathbf{A}{\hat{\mathbf{\Sigma }}}_{B,1}{\mathbf{A}}^T+{\hat{\mathbf{\Sigma }}}_{U,1}+{\mathbf{AM}}_1+{\mathbf{M}}_1^T{\mathbf{A}}^T).$$

$$\left(\begin{array}{c}{\mathbf{B}}_{t+h}\\ {\hat{\mathbf{U}}}_{t+h}\end{array}\right)\mid {\mathcal{I}}_{t,b}$$

$$\mathbf{G}$$

$$=\left({\hat{\mathbf{\Sigma }}}_{B,1}{\mathbf{A}}^T+{\mathbf{M}}_1\right){\left(\mathbf{A}{\hat{\mathbf{\Sigma }}}_{B,1}{\mathbf{A}}^T+{\hat{\mathbf{\Sigma }}}_{U,1}+{\mathbf{AM}}_1+{\mathbf{M}}_1^T{\mathbf{A}}^T\right)}^{-1},$$

$${\tilde{\mathbf{b}}}_{t+h}=\mathbb{E}\left[{\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{u}}}_{t+h}\right]={\hat{\mathbf{b}}}_{t+h}+\mathbf{G}({\hat{\mathbf{u}}}_{t+h}-\mathbf{A}{\hat{\mathbf{b}}}_{t+h})$$

$$\mathrm{Var}\left[{\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{u}}}_{t+h}\right]=k_h\left({\hat{\mathbf{\Sigma }}}_{B,1}-\mathbf{G}(\mathbf{A}{\hat{\mathbf{\Sigma }}}_{B,1}+{\mathbf{M}}_1^T)\right)$$

$$\mathbb{E}[{\mathbf{Y}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{u}}}_{t+h}]={\tilde{\mathbf{y}}}_{t+h}=\mathbf{S}{\tilde{\mathbf{b}}}_{t+h}$$

$$\mathrm{Var}\left[{\mathbf{Y}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{u}}}_{t+h}\right]=\mathbf{S}\mathrm{Var}\left[{\mathbf{B}}_{t+h}\mid {\mathcal{I}}_{t,b},{\hat{\mathbf{u}}}_{t+h}\right]{\mathbf{S}}^T.$$

$$\mathbb{E}[\mathrm{Cov}(B_{t+1}^i,B_{t+1}^j\mid {\mathbf{B}}_{1:t})]=\mathrm{Cov}(E_{t+1}^i,E_{t+1}^j)i,j=1,\dots ,n,$$

$$\mathrm{Cov}(B_{t+1}^i,B_{t+1}^j\mid {\mathbf{B}}_{1:t})$$

$$=\mathbb{E}\left[\left(B_{t+1}^i-{\hat{B}}_{t+1}^i\right)\left(B_{t+1}^j-{\hat{B}}_{t+1}^j\right)\mid {\mathbf{B}}_{1:t}\right]$$

$$\mathbb{E}[\mathrm{Cov}(B_{t+1}^i,B_{t+1}^j\mid {\mathbf{B}}_{1:t})]$$

$$=\mathbb{E}\left[\left(B_{t+1}^i-{\hat{B}}_{t+1}^i\right)\left(B_{t+1}^j-{\hat{B}}_{t+1}^j\right)\right]$$

$$=\mathrm{Cov}\left(B_{t+1}^i-{\hat{B}}_{t+1}^i,B_{t+1}^j-{\hat{B}}_{t+1}^j\right)=\mathrm{Cov}\left(E_{t+1}^i,E_{t+1}^j\right)$$

$${\mathbf{\epsilon }}_{t+1}^u=\mathbf{A}{\mathbf{B}}_{t+1}-{\hat{\mathbf{U}}}_{t+1},$$

$$\left(\begin{array}{c}{B}_1\\ B_2\end{array}\right)$$

Full text

1 Introduction

Often time series are organized into a hierarchy. For example, the total visitors of a country can be divided into regions and the visitors of each region can be further divided into sub-regions. The most disaggregated time series of the hierarchy are referred to as bottom time series, while the remaining time series are referred to as upper time series.

Forecasts of hierarchical time series should be coherent; for instance, the sum of the forecasts of the different regions should equal the forecast for the total. The forecasts are incoherent if they do not satisfy such constraints. A simple way for generating coherent forecasts is bottom-up: one takes the forecasts for the bottom time series and sums them up according to the summing constraints in order to produce forecasts for the entire hierarchy. Yet this approach does not consider the forecasts produced for the upper time series, which contain useful information. For instance, upper time series are smoother and allow to better estimate of the seasonal patterns and the effect of the covariates.

Thus, modern reconciliation methods [9, 18] proceed in two steps. First, base forecasts are computed by fitting an independent model to each time series. Then, the base forecasts are adjusted to become coherent; this step is called reconciliation. The forecasts for the entire hierarchy are then obtained by summing up the reconciled bottom time series. Reconciled forecasts are generally more accurate than the base forecasts, as they benefit from information coming from multiple time series. The state-of-the art reconciliation algorithm is MinT [18], which minimizes the mean squared error of the reconciled forecasts by solving a generalized least squares problem; its point forecast, besides being coherent, are generally more accurate than the base forecast.

Hierarchical probabilistic forecasting is however still an open area of research. The algorithm by [17] constructs a coherent forecast in a bottom-up fashion, modelling via copulas the joint distribution of the bottom time series, while [13] proposes a top-down approach, where the top time series is forecasted and then disaggregated. Both algorithms are based on numerical procedures which have no closed-form solution; hence they are not easily interpretable. In [6], a geometric interpretation of the reconciliation process is provided. It is moreover shown that the log score is not proper with respect to incoherent probabilistic forecasts. As an alternative, the energy score can be used for comparing reconciled to unreconciled probabilistic hierarchical forecasts. In [1] multivariate Gaussian predictive densities and bootstrap densities are experimentally compared for hierarchical probabilistic forecasting.

We address probabilistic reconciliation using Bayes’ rule. We define the prior beliefs about the bottom time series, based on the base forecasts for the bottom time series. We then update them incorporating the information contained in the forecasts for the upper time series. Under the Gaussian assumption, we compute the update in closed form, obtaining the posterior distribution about the bottom time series and then about the entire hierarchy. Our reconciled forecasts minimize the mean squared error; indeed, we prove that they match the point predictions of MinT, whose optimality has been proven in a frequentist way. Our algorithm provides the joint predictive distribution for the hierarchy; thus we call it pMinT, which stands for probabilistic MinT. We also provide a variant of pMinT, obtained by making an additional independence assumption; we call it LG, as it is related to the linear-Gaussian model [3, Chap.8.1.4].

We show a link between the reconciliation problem and the Kalman filter, opening the possibility of borrowing from the literature of the Kalman filter for future research. We then compare the algorithms on synthetic and real data sets, eventually concluding that pMinT yields more accurate probabilistic forecasts than both bottom-up and LG.

The paper is organized as follows: we introduce the reconciliation problem in Section 2, we discuss the algorithms in Sec. 3, we discuss the reconciliation of a simple hierarchy in Section 4 and we present the experiments in Section 5.

2 Time series reconciliation

Fig. 1 shows a hierarchy. We could interpret it as the visitors of a country, which are disaggregated first by region (R1, R2) and then by sub-regions (R11, R12, R21, R22). The most disaggregated time series (bottom time series) are shaded. The hierarchy contains $m$ time series, of which $n$ are bottom time series.

We denote by uppercase letters the random variables and by lowercase letters their observations. The vector of observations available at time $t$ for the entire hierarchy is $\mathbf{y}_{t}\in\mathbb{R}^{m}$; they are observations from the set of random variables $\mathbf{Y}_{t}$. Vector $\mathbf{y}_{t}$ can be broken down in two parts, namely $\mathbf{y}_{t}=[\mathbf{u}_{t}^{T},\mathbf{b}_{t}^{T}]^{T}$; $\mathbf{b}_{t}\in\mathbb{R}^{n}$ contains the observations of the bottom time series while $\mathbf{u}_{t}\in\mathbb{R}^{m-n}$ contains the observations of the upper time series. At time $t$, the observations available for the hierarchy of Fig. 1 are thus:

[TABLE]

The structure of the hierarchy is represented by the summing matrix $\mathbf{S}\in\mathbb{R}^{m\times n}$ such that:

[TABLE]

The $\mathbf{S}$ matrix of hierarchy in Fig.1 is:

[TABLE]

where the sub-matrix $\mathbf{A}\in\mathbb{R}^{(m-n)\times n}$ encodes which bottom time series should be summed up in order to obtain each upper time series.

We denote by $\widehat{\mathbf{y}}_{t+h}\in\mathbb{R}^{m}$ the base forecasts issued at time $t$ about of $y$ and referring to $h$ steps ahead. We separate base forecasts for bottom time series ($\widehat{\mathbf{b}}_{t+h}\in\mathbb{R}^{n}$) and upper time series ($\widehat{\mathbf{u}}_{t+h}\in\mathbb{R}^{m-n}$), namely $\widehat{\mathbf{y}}_{t+h}=[\widehat{\mathbf{u}}_{t+h}^{T},\widehat{\mathbf{b}}_{t+h}^{T}]^{T}$. The variances of the error of the base forecasts will be used later. If forecasts for different time horizons are needed (e.g., $h$=1,2,3,..), the reconciliation is performed independently for each $h$. In the following we generically assume to reconcile the forecasts for $h$ steps ahead.

The MinT reconciliation

Most reconciliation algorithms [9], including MinT [18], assume the reconciled bottom forecasts ($\widetilde{\mathbf{b}}_{t+h}$) to be a linear combination of the base forecasts ($\widehat{\mathbf{y}}_{t+h}$) available for the whole hierarchy, i.e. their objective is to find a matrix $\mathbf{P}_{h}\in\mathbb{R}^{n\times m}$ such that:

[TABLE]

Let us denote by $\widehat{\mathbf{E}}_{t+h}=\mathbf{Y}_{t+h}-\widehat{\mathbf{y}}_{t+h}\in\mathbb{R}^{m}$ the vector of the errors of the base forecast $h$-steps ahead and by $\mathbf{W}_{h}=\mathbb{E}[\widehat{\mathbf{E}}_{t+h}\ \widehat{\mathbf{E}}^{T}_{t+h}\mid\mathcal{I}_{t}]$ their covariance matrix, where $\mathcal{I}_{t}$ denotes all the information available up to time $t$. In [18] it is proven that the reconciliation matrix given by:

[TABLE]

is optimal, in the sense that it minimizes the trace of the reconciliation errors’ covariance matrix. The reconciled forecasts for the whole hierarchy are obtained by summing the reconciled bottom forecasts, and they are proven to minimize the mean squared error over the entire hierarchy.

2.0.1 Estimation of $\mathbf{W}_{h}$

Estimating $\mathbf{W}_{h}$ differently for each $h$ is an open problem. For the case $h$=1 the estimation is simpler. The variance of the forecasts equals the variance of the residuals (i.e., the errors on 1-step predictions made on the training data) and cross-covariances are estimated as the covariance of the residuals. The best estimates are obtained [18] by shrinking the full covariance matrix towards a diagonal matrix, using the method of [15].

The case $h$ >1 is instead problematic. The variance of the forecasts are obtained by increasing the 1-step variance through analytical formulas, which differ from the variance of the $h$-steps ahead residuals. Morevoer, the covariances in $\mathbf{W}_{h}$ have to be numerically estimated by looking at the $h$-steps residuals. However, the number of $h$-steps residuals decreases with $h$, making the estimate more noisy.

As a workaround, [18] assumes $\mathbf{W}_{h}=k_{h}\mathbf{W}_{1}$, where $k_{h}>0$ is an unknown constant which depends on $h$ while $\mathbf{W}_{1}$ is the covariance matrix of the one-step ahead errors. The underlying assumption is thus that all terms within the variance/covariance matrix of the errors grow in the same way with $h$. The advantage of this approach is that $k_{h}$ cancels out when computing the reconciled forecasts, as it can be seen by setting $\mathbf{W}_{h}=k_{h}\mathbf{W}_{1}$ in Eq.(4). Yet, $k_{h}$ appears in the expression of the variance of the reconciled forecasts. In the following we refer to the assumption $\mathbf{W}_{h}=k_{h}\mathbf{W}_{1}$ as “the $k_{h}$ assumption”.

3 Probabilistic Reconciliation

We address the reconciliation problem by merging the probabilistic information contained in the base forecasts for the bottom and the upper time series. We perform the fusion using Bayes’ rule.

We first define the prior about the bottom time series. We have observed the time series up to time $t$ and we are interested in the reconciled forecasts for time $t+h$. We denote by $\mathbf{B}_{t+h}$ the vector of the bottom time series at time $t+h$; this is thus a vector of random variables and $B^{i}_{t+h}$ represents its $i$th element. We moreover denote by $\widehat{\mathbf{b}}_{t+h}$ the vector of base forecasts for the bottom time series for time $t+h$, and by $\mathbf{b}_{t+h}$ the actual observation of the bottom random variables at time $t+h$. Finally, $\mathcal{I}_{t,b}$ is the information available up to time $t$ regarding the bottom time series, i.e. the past values of the bottom time series: $\mathcal{I}_{t,b}=\{\mathbf{b}_{1},\ldots,\mathbf{b}_{t}\}$.

In the following we adopt the $k_{h}$ assumption for all the covariance matrices assuming moreover that, for a given $h$, the value of $k_{h}$ is shared among all the involved covariance matrices. As we will show later, this is equivalent to the $k_{h}$ assumption made by MinT. Let us hence denote the covariance matrix of the forecast $h$-steps ahead by $\widehat{\boldsymbol{\Sigma}}_{B,h}=k_{h}\widehat{\boldsymbol{\Sigma}}_{B,1}$. Assuming the bottom time series to be jointly Gaussian we have:

[TABLE]

Probabilistic bottom-up

If we have no information about the upper time series, we can build a joint predictive distribution for the entire hierarchy by summing the bottom forecast via matrix $\mathbf{S}$:

[TABLE]

which is a probabilistic bottom-up reconciliation. Note that, in this case, $k_{h}$ appears only in the expression of the variance.

Updating

If the forecasts $\widehat{\mathbf{U}}_{t+h}$ about the upper time series are available, then we can use them in order to update our prior. We assume:

[TABLE]

where $\widehat{\boldsymbol{\Sigma}}_{U,h}=k_{h}\widehat{\boldsymbol{\Sigma}}_{U,1}$ is the covariance of the noise. We thus treat $\widehat{\mathbf{U}}_{t+h}$ as a set of different sums of the future values of the bottom time series, corrupted by noise. Hence:

[TABLE]

The posterior distribution of the bottom time series is given by Bayes’ rule:

[TABLE]

3.1 Computing Bayes’ rule

The posterior of Eq. (9) can be computed in closed form by assuming the vector $(\mathbf{B}_{t+h},\widehat{\mathbf{U}}_{t+h})$ to be jointly Gaussian distributed. The linear-Gaussian (LG) model [3, Chap.8.1.4]. computes analytically the updating by further assuming $\boldsymbol{\varepsilon}_{t+h}^{u}$ to be independent from $\mathbf{B}_{t+h}$. Yet this independence might not always hold in our case. Consider for instance a special event driving upwards most time series. As a result we would observe both high values of $\mathbf{B}_{t+h}$ and negative values of $\boldsymbol{\varepsilon}_{t+h}^{u}$, due to the underestimation of the upper time series. This would result in a correlation between $\mathbf{B}_{t+h}$ and $\boldsymbol{\varepsilon}_{t+h}^{u}$.

We thus generalize the LG model by accounting for such correlation. We will later compare experimentally the results obtained adopting the LG model and its generalized version. We denote $\operatorname{Cov(\mathbf{B}_{t+1},\boldsymbol{\varepsilon}_{t+1}^{u}\mid\mathcal{I}_{t,b}})=\mathbf{M}_{1}\in\mathbb{R}^{n\times(m-n)}$ and we assume $\operatorname{Cov(\mathbf{B}_{t+h},\boldsymbol{\varepsilon}_{t+h}^{u}\mid\mathcal{I}_{t,b}})=k_{h}\mathbf{M}_{1}$.

Our first step for computing Bayes’ rule is to express the joint distribution $p(\mathbf{B}_{t+h},\widehat{\mathbf{U}}_{t+h}\mid\mathcal{I}_{t,b})$. Since $\widehat{\mathbf{U}}_{t+h}=\mathbf{A}\mathbf{B}_{t+h}+\boldsymbol{\varepsilon}_{t+h}^{u}$, the expected values are:

[TABLE]

We now derive the different blocks of the covariance matrix. The cross-covariance between $\mathbf{B}_{t+h}$ and $\widehat{\mathbf{U}}_{t+h}$ is:

[TABLE]

where $k_{h}>0$ is the multiplicative constant (Sec.2.0.1) that yields the variance of the forecasts $h$-steps ahead, given the covariances of the forecasts 1-step ahead.

The covariance of the upper forecasts is:

[TABLE]

Hence the joint prior (i.e., before observing $\widehat{\mathbf{U}}_{t+h}$) is:

[TABLE]

Now we receive the forecast $\widehat{\mathbf{u}}_{t+h}$ for the upper time series (recall that $\widehat{\mathbf{U}}_{t+h}$ denotes the random variables while $\widehat{\mathbf{u}}_{t+h}$ denotes observations). We obtain the posterior distribution for the bottom time series $P(\mathbf{B}_{t+h}|\widehat{\mathbf{u}}_{t+h})$ by applying the standard formulas for the conditional distribution of a MVN distribution [12, Sec.4.3.1]. To have a shorter notation, let us define:

[TABLE]

where $k_{h}$ disappears from the expression of $\mathbf{G}$.

The reconciled bottom time series have then the following mean and variance:

[TABLE]

Thus the adjustment applied to the base forecasts is proportional to ($\widehat{\mathbf{u}}_{t+h}-\mathbf{A}\widehat{\mathbf{b}}_{t+h}$), i.e. the difference between the prior mean and the uncertain observation (i.e., the forecasts) of the upper time series. The term ($\widehat{\mathbf{u}}_{t+h}-\mathbf{A}\widehat{\mathbf{b}}_{t+h}$) is called the incoherence of the base forecasts in [18]. The mean of the reconciled bottom time series does not depend on $k_{h}$, while the variance does.

The reconciled point forecast and the covariance for the entire hierarchy are:

[TABLE]

3.2 Related works and optimality of the reconciliation

Bayes’ rule is a well-known tool for information fusion [5, Sec. 2], and we apply it for the first time for forecast reconciliation. We will later prove that the posterior mean of our approach yields the same point predictions of MinT.

Yet, our algorithm additionally provides the predictive distribution for the entire hierarchy; we thus call it pMinT, where p stands for probabilistic. We also contribute a novel reconciliation approach based on the linear-Gaussian (LG) model, which is obtained by setting $\mathbf{M}=\mathbf{0}$ in the definition of $\mathbf{G}$. Both pMinT and LG are thus probabilistic reconciliation algorithms.

We point out for the first time a link between the reconciliation problem and the Kalman filter, whose state-update equation can be derived from the linear-Gaussian model [14]. In particular, Equation (11) has the same structure of the state-update of a Kalman filter. According to the definition of $\mathbf{G}$ of Eq.(10), the LG reconciliation corresponds to the standard Kalman filter [16, Chap.5], while the pMinT reconciliation corresponds to a generalized Kalman filter which assumes correlation between the noise of the state and the noise of the output [16, Chap.7.1]. Thus future research could explore the literature of the Kalman filter in order to borrow ideas for the reconciliation problem.

The optimality of our approach can be informally proven by considering that it yields the posterior mean of the reconciled forecasts, which is the minimizer of the quadratic loss under the Gaussian assumption [12, Chap. 5.7]. Moreover, it has the same equation of the state-update step of the Kalman filter, which provably minimizes the mean squared error of the estimates without any distributional assumption [16, Chap.5]. In Sec. (4.1) we will moreover show that our point predictions correspond to those of MinT, which have been proven to be the minimizer of the mean-squared error.

3.3 The covariance matrices $\widehat{\boldsymbol{\Sigma}}_{B,1}$ and $\widehat{\boldsymbol{\Sigma}}_{U,1}$

The element $(i,j)$ of $\widehat{\boldsymbol{\Sigma}}_{B,1}$ is the covariance $\operatorname{Cov}(B^{i}_{t+1},B^{j}_{t+1}\mid\mathcal{I}_{t,b})=\operatorname{Cov}(B^{i}_{t+1},B^{j}_{t+1}\mid\mathbf{B}_{1:t}=\mathbf{b}_{1:t})$, where $\mathbf{B}_{1:t}=\mathbf{b}_{1:t}$ denotes a realization $\mathbf{b}_{1:t}$ of $\mathbf{B}_{1:t}$. Yet we only have one observation of $B^{i}_{t+1},B^{j}_{t+1}$ conditional on $\mathcal{I}_{t,b}$, which prevents estimating the covariance. We can overcome the problem with the following result, which shows that we can approximate $\widehat{\boldsymbol{\Sigma}}_{B,1}$ by computing the covariance of the residuals.

Let us consider the vectors of bottom time series $\mathbf{B}_{1},\ldots,\mathbf{B}_{t}$ and the conditional expectation $\widehat{\mathbf{B}}_{t+1}=\mathbb{E}[\mathbf{B}_{t+1}\mid\mathbf{B}_{1},\ldots,\mathbf{B}_{t}]=\mathbb{E}[\mathbf{B}_{t+1}\mid\mathbf{B}_{1:t}]$. Note that $\widehat{\mathbf{B}}_{t+1}$ is a random vector as we have not yet observed $\mathbf{B}_{i}$, $i=1,\ldots,t$. In this section we show

[TABLE]

where $E^{i}_{t+1}:=B^{i}_{t+1}-\widehat{B}^{i}_{t+1}$, $i=1,\ldots,n$ denotes the residual of the model fitted on the $i$-th time series, for the forecast horizon $t+1$. If we observe $\mathcal{I}_{t,b}$, then we can approximate $\operatorname{Cov}(B^{i}_{t+1},B^{j}_{t+1}\mid\mathcal{I}_{t,b})$ with $\operatorname{Cov}(B^{i}_{t+1}-\widehat{b}^{i}_{t+1},B^{j}_{t+1}-\widehat{b}^{j}_{t+1})$, the covariance of the residuals of the models fitted on the bottom time series.

Consider now the conditional covariance on the left side of eq. (15), we have

[TABLE]

By taking the expectation on both sides we obtain

[TABLE]

We thus estimate the covariance matrix $\widehat{\boldsymbol{\Sigma}}_{B,1}$ using the covariance of the residuals of the models fitted on the bottom time series.

3.3.1 Computation of $\widehat{\boldsymbol{\Sigma}}_{U,1}$

According to Eq. (7),

[TABLE]

whose variances and covariances can be readily computed from the residuals.

4 Reconciliation of a simple hierarchy

We now illustrate how the base forecasts interact during the reconciliation of a simple hierarchy. We consider a hierarchy constituted by two bottom time series ($B_{1}$ and $B_{2}$) and an upper time series $U$.

The base forecast for the bottom time series are the point forecasts $\widehat{b}_{1}$ and $\widehat{b}_{2}$ with variances $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$. The prior beliefs about $B_{1}$ and $B_{2}$ are:

[TABLE]

where for simplicity we remove the forecast horizon ($t+h$) from the notation.

The summing matrix is:

[TABLE]

We start considering the simpler case of reconciliation via the LG algorithm. The matrix $\mathbf{G}$ is:

[TABLE]

since $\mathbf{A}=[1\,\,1]$, $\widehat{\boldsymbol{\Sigma}}_{U,1}=\sigma_{u}^{2}$ and moreover:

[TABLE]

Note that $\mathbf{G}$ does not depend on $h$, as also shown in Eq. (10).

The reconciled bottom forecasts are:

[TABLE]

where $\widehat{u}$ is the base forecast for $U$.

The reconciled bottom forecast can be written as:

[TABLE]

Eq. (17) shows that the adjustment applied to the base forecasts depends on $\sigma_{u}^{2}$. If $\sigma_{u}^{2}$ is large the adjustment is small, since the upper forecast is not informative. If on the contrary $\sigma_{u}^{2}=0$, the sum of the reconciled bottom forecasts is forced to match $\hat{u}$, i.e., $\tilde{b}_{1}+\tilde{b}_{2}=\hat{u}$ (this can be shown by re-working Eq. (17)).

We now show that the reconciled bottom forecast are a linear combination of the base forecasts. Let us define:

[TABLE]

After some algebra we obtain:

[TABLE]

Thus $\widetilde{b}_{1}$ is a weighted average of two estimates: $\widehat{b}_{1}$ and $(\widehat{u}-\widehat{b}_{2})$; the weight of $\widehat{b}_{1}$ decreases with $\sigma_{1}^{2}$ and increases with ($\sigma_{2}^{2}+\sigma_{u}^{2})$.

The reconciliation carried out by pMinT is similar to what already discussed, once we adopt $g_{1}^{*}$ and $g_{2}^{*}$ in place of $g_{1}$ and $g_{2}$:

[TABLE]

Thus pMinT accounts also for the cross-covariances $\sigma_{u,1},\sigma_{u,2}$ between the bottom time series and of noise affecting the forecasts for the upper time series.

4.1 Relationship with MinT

Our reconciled bottom time series can be written as:

[TABLE]

The matrix $\mathbf{P}_{h}$ of pMinT is thus:

[TABLE]

Proposition 1

The matrices $\mathbf{P}_{h}$ of MinT and pMinT are equivalent. The proof is given in the supplementary material.

5 Experiments

In this paper we take a probabilistic point of view; we thus assess the reconciled predictive distributions rather than the point forecasts. Our metric is the energy score (ES), a scoring rule for multivariate distributions [7]. The ES is the multivariate generalization of the continuous ranked probability score (CRPS), which is obtained by integrating the Brier score over the predictive distribution of the forecast [7]. Let $\mathbf{y}$ be the actual multivariate observation, and let us assume that we have $k$ samples $\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{k}$, from the multivariate predictive distribution $F$. The energy score is:

[TABLE]

The energy score is a loss function: the lower, the better. We consider three methods for probabilistic reconciliation: probabilistic bottom-up (BU), LG and pMinT. We did not find any package implementing the algorithms of [17, 13]; thus we did not include them in our comparison. We estimate all the covariance matrices via the shrinkage estimator [15].

Base forecasts

We consider two time series models to compute the base forecasts. The first is ets, which fits different exponential smoothing variants and eventually performs model selection via AICc. The second method is auto.arima. It first decides how to differentiate the time series to make it stationary; then, it looks for the best arma model which fits the stationary time series, performing model selection via AICc. Both approaches performed well in forecasting competitions; they are available from the forecast package [10] for R. In all simulations, we compute forecasts up to $h$=4. As no particular pattern exists in the relative performance of the methods as $h$ varies, we present the results averaged over $h$=1,2,3,4.

Setting $k_{h}$

We are unaware of previous studies on how to set $k_{h}$. In this paper we compare two heuristics, acknowledging that this remains an open problem. The two options are $k_{h}$=$h$ and $k_{h}$=$1$. The choice $k_{h}$=$h$ is based on the following approximation. The variance of $\widehat{y}_{t+1}$ around $y_{t+1}$ is $\sigma^{2}$; assuming the independence of the errors, the variance of $\widehat{y}_{t+1}$ around $y_{t+h}$ is $h\sigma^{2}$. The approximation lies in the fact that we are modeling the variance of $\widehat{y}_{t+1}$ (not $\widehat{y}_{t+h}$) around $y_{t+h}$.

Instead, the option $k_{h}=1$ keeps the variance fixed with $h$. This represents the short-term behavior of models which contain only seasonal terms and no autoregressive terms. For instance, when dealing with a monthly time series, the variance of such models is constant up to $h$=12.

Code

The code of our experiments is available at: https://github.com/iamthejao/BayesianReconciliation.

5.1 Synthetic data

We generate synthetic data sets using the hierarchy above, previously considered in the experiments of [18]. We simulate the four bottom time series as AR(1) processes, drawing their parameters uniformly from the stationary region. The noises of the bottom time series at each time instant are correlated, multivariate Gaussian distributed, with mean $\mu=[0,0,0,0]^{T}$and covariance:

[TABLE]

Thus $\Sigma$ enforces a stronger correlation between time series which have the same parents. At each time instant $t$ we add the noise $\eta_{t}\sim N(0,10)$ to the time series AA and BA and the noise $(-\eta_{t})$ to the time series AB and BB. In this way we simulate noisy bottom time series ($\eta_{t}$ and $-\eta_{t}$ cancel out when dealing with the upper time series) which can be encountered in real cases when several disaggregations are applied to the total time series. We consider the following length $T$ of the time series: $\{50;100;1000\}$. For each value of $T$ we perform 1000 simulations. The averaged energy scores are given in Tab. 1; in each cell we report the lower energy score between the case $k_{h}$=$h$ and $k_{h}$=$1$.

Since the time series are stationary, they basically fluctuate around their mean. In this case the magnitude of the incoherence is generally limited, allowing also the bottom-up reconciliation to be competitive. The ES of both pMinT and LG is on average 1.5% smaller than that of BU. We also note an advantage of LG over pMinT for small $T$, and instead the reverse for large $T$; this might be the effect of the additional covariances estimated by pMinT (see $\sigma_{u,1}$ and $\sigma_{u,2}$ in Sec.4).

5.2 Experiments with real data sets

We consider two hierarchical time series: infantgts and tourism. Both are grouped time series, which is a generalization of hierarchical time series. In particular the time series of a given level are always sums of some bottom time series, but they are not necessarily sums of time series of the adjacent lower level.

Infantgts

The infantgts is available within the hts [8] package for R. It contains infant mortality counts in Australia, disaggregated by sex and by eight different states. Each time series contains 71 yearly observations, covering the period 1933-2003. The bottom level contains 16 time series (8 states x 2 genders). The second level contains 2 time series: the counts of males and females, aggregated over the states. The third level sums males and females in each state, yielding 8 time series (one for each state). The fourth level is the total.

Tourism

The tourism data set regards the number of nights spent by Australians away from home. It is available in raw format from https://robjhyndman.com/publications/MinT/. The time series cover the period 1998–2016 with monthly frequency. There are 304 bottom time series, referring to 76 regions and 4 purposes. The first level sums over the purposes, yielding 76 time series (one for each region); such values are further aggregated into macro-zones (27 time series) and states (7 time series). Other levels of the hierarchy aggregate the bottom time series of the same zone (yielding 108 time series: 27 zones x 4 purposes), which are then further aggregated into 28 time series (7 states x 4 purposes) and then 4 time series (4 purposes). The last level is the total. Overall the hierarchy contains 555 time series.

We repeat 50 times the following procedure: split the time series into training and test, using a different split point; compute the base forecasts up to $h$=4; reconcile the forecasts. The reconciliation is independently computed for each $h$. Each value of Tab. 2 is thus the average over 200 experiments (50 training/test splits $\times$ $h$=1,2,3,4). On infantgts, all the three reconciliation methods perform better with $k_{h}=h$, probably because the variance of the fitted time series models steadily increases with $h$. On the contrary, on tourism all reconciliation algorithms perform better with $k_{h}$=1; in this case most models only contain the seasonal part. The rows referring to the best values of $k_{h}$ are highlighted in Tab. 2.

We call setup the combination of a data set and a forecasting method, such as <infangts,arima >. The pMinT algorithm yields the lowest energy score in most setups; in the next section we check whether the differences between methods are significant.

Statistical analysis

For each setup we perform a significance tests for each pair of algorithms (pMinT vs BU, LG vs BU and pMinT vs LG), using the Bayesian signed-rank test [2], which returns the posterior probability of a method having lower median energy score than another (Tab. 3). Such posterior probabilities are numerically equivalent to (1 - p-value), where p-value is the p-value of the one-sided frequentist signed-rank test.

In most setups (Tab. 3) high posterior probabilities (implying low p-values) support the hypothesis of the pMinT having lower energy score than both BU and LG; moreover they also support the hypothesis of LG having lower energy score than BU.

Meta-analysis

We now perform a meta-analysis for for each pair of algorithms across the different setups, adopting the Poisson-binomial approach [11, 4]. Consider for instance pMinT vs BU. We model each setup as a Bernoulli trial, whose possible outcomes are the victory of pMinT or BU. The probability of pMinT winning is taken for each setup from Tab. 3 (the probability of BU winning is just its complement to 1). We then repeat 10,000 simulations, in which we draw the outcome of each setup according to the probabilities of Tab. 3.

We now report the probability of each method outperforming another method in more than half the setups, based on out of 10,000 simulations. Both pMinT and LG wins in more than half the setup with probability 1 against BU. Moreover, there is 0.85 probability of pMinT winning in more than half of the setups against LG. We thus recommend pMinT as a general default method for probabilistic reconciliation.

6 Conclusions

We have derived two algorithms (pMinT and LG) based on Bayes’ rule for probabilistic reconciliation. We have also shown a didactic example which clarifies how base forecast and their variances interact during the reconciliation, In general pMinT yields better predictive distributions and thus we recommend it as a default. The LG method can be anyway an interesting alternative when dealing with small sample sizes. Future research could borrow ideas from the extensive literature of the Kalman filter, based on the link we pointed out between reconciliation and Kalman filter

Acknowledgements

We acknowledge support from grant n. 407540_167199 / 1 from Swiss NSF (NRP 75 Big Data).