Forecast reconciliation is a post-forecasting process that involves transforming a set of incoherent forecasts into coherent forecasts which satisfy a given set of linear constraints for a multivariate time series. In this paper we extend the current state-of-the-art cross-sectional probabilistic forecast reconciliation approach to encompass a cross-temporal framework, where temporal constraints are also applied. Our proposed methodology employs both parametric Gaussian and non-parametric bootstrap approaches to draw samples from an incoherent cross-temporal distribution.
To improve the estimation of the forecast error covariance matrix, we propose using multi-step residuals, especially in the time dimension where the usual one-step residuals fail.
To address high-dimensionality issues, we present four alternatives for the covariance matrix, where we exploit the two-fold nature (cross-sectional and temporal) of the cross-temporal structure, and introduce the idea of overlapping residuals. We evaluate the proposed methods through a detailed simulation study that investigates their theoretical and empirical properties. We further assess the effectiveness of the proposed cross-temporal reconciliation approach by applying it to two empirical forecasting experiments, using the Australian GDP and the Australian Tourism Demand datasets. For both applications, we show that the optimal cross-temporal reconciliation approaches significantly outperform the incoherent base forecasts in terms of the Continuous Ranked Probability Score and the Energy Score. Overall, our study expands and unifies the notation for cross-sectional, temporal and cross-temporal reconciliation, thus extending and deepening the probabilistic cross-temporal framework. The results highlight the potential of the proposed cross-temporal forecast reconciliation methods in improving the accuracy of probabilistic forecasting models.

Keywords: Coherent, GDP, Linear constraints, Multivariate time series, Temporal aggregation, Tourism flows

§ INTRODUCTION

Forecast reconciliation is a post-forecasting process intended to improve the quality of forecasts for a system of linearly constrained multiple time series <cit.>. There are many fields where forecast reconciliation is useful, such as when forecasting GDP and its components, electricity demand and power generation, demand in supply chains with product categories, tourist flows across geographic regions and travel purpose, and more. Moreover, effective decision-making depends on the support of accurate and coherent forecasts.

Classical reconciliation methods addressed the issue of incoherent forecasts in a cross-sectional hierarchy by forecasting only one level and using these to generate forecasts for the remaining series. The bottom-up approach <cit.> starts by generating forecasts at the most disaggregate level and summing these to arrive at the desired forecasts for aggregate levels. On the other hand, the top-down approach <cit.> forecasts the most aggregated level and then disaggregates it to lower levels <cit.>. The middle-out method <cit.> combines both approaches by selecting an intermediate level and applies top-down for lower levels and bottom-up for upper levels.

All of these approaches ignore useful information available at other levels <cit.>. Consequently, in the last decade, hierarchical forecasting has significantly evolved to include modern least squares-based reconciliation techniques in the cross-sectional framework <cit.>, later extended to temporal hierarchies <cit.>. Obtaining coherent forecasts across both the cross-sectional and temporal dimensions (known as cross-temporal coherence) has been limited to sequential approaches that address each dimension separately <cit.>. Recently, <cit.> suggested a unified reconciliation step that takes into account both the cross-sectional and temporal dimensions, instead of dealing with them separately, utilizing the entire cross-temporal hierarchy.

However, these cross-temporal works focus on point forecasting, and do not consider distributional or probabilistic forecasts <cit.>. In the cross-sectional framework, there have been some developments towards probabilistic forecasting including  <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. <cit.> made a significant contribution by formalizing cross-sectional probabilistic reconciliation using the geometric framework for point forecast reconciliation proposed by <cit.>. They show how a reconciled forecast can be constructed from an arbitrary base forecast when the density of the base forecast is available and when only a sample can be drawn. They also show that in the case of elliptical distributions, the correct predictive distribution can be recovered via linear reconciliation, regardless of the base forecast location and scale parameters, and derive conditions for this to hold in the special case of reconciliation via projection.

In this paper, we extend cross-sectional probabilistic reconciliation to the cross-temporal case, working on issues related to the two-fold nature of this framework. First, we revise and develop the notation proposed by <cit.> to generalize the work of <cit.>. This allows us to move from cross-temporal point reconciliation to a probabilistic setting through the generalization of definitions and theorems well-established in the cross-sectional framework. Second, we propose effective and practical solutions to draw a sample from the base forecast distribution according to either a parametric approach that assumes Gaussianity or a non-parametric approach that bootstraps the base model residuals. Third, we propose some solutions to specific problems that arise when combining the cross-sectional and temporal dimensions. We propose using multi-step residuals to estimate the relationships between different forecast horizons when we deal with temporal levels, since one-step residuals are not suitable for this purpose. To solve high-dimensionality issues we introduce the idea of overlapping residuals and consider alternative forms for constructing the covariance matrix. Fourth, we propose new shrinkage procedures for reconciliation that aim to identify a feasible cross-temporal structure. The methodological contributions described in this paper are implemented in the  package for R <cit.>.

The remainder of the paper is structured as follows. In <ref>, we provide a unified notation for the cross-sectional, temporal and cross-temporal point reconciliation. We generalize the cross-sectional definitions and theorems developed by <cit.> in <ref>, and propose both a parametric Gaussian and a non-parametric bootstrap approach to draw a sample from the base forecast distribution. In <ref>, we analyze the structure of the cross-temporal covariance matrix, proposing four alternative forms, and propose shrinkage approaches for reconciliation. In addition, we explore cross-temporal residuals (overlapping and multi-step) looking at their advantages and limitations. A simulation study is performed in <ref>, to better understand the properties of the methodology. Two empirical applications using the Australian GDP and the Australian Tourism Demand datasets are considered in Sections <ref> and <ref>, respectively[A complete set of results is available at the GitHub repository <https://github.com/danigiro/ctprob>.]. Finally, <ref> presents conclusions and a future research agenda on this and other related topics.

§ NOTATION AND DEFINITIONS

Let $\yvet_t = [y_{1,t},\dots,y_{i,t},\dots,y_{n,t}]'$ be an $n$-variate linearly constrained time series observed at the most temporally disaggregated level, with a seasonality of period $m$ (e.g., $m = 12$ for monthly data, $m = 4$ for quarterly data, $m = 24$ for hourly data). Suppose that the constraints are expressed by linear equations such that <cit.>

\begin{equation}
\label{eq:cs_con}
\Cvet_{cs}\yvet_t = \Zerovet_{(n_a \times 1)}, \qquad t = 1, \;\dots, \;T,
\end{equation}

where $\Cvet_{cs}$ is the $(n_a \times n)$ zero constraints cross-sectional matrix, that can be seen as the coefficient matrix of a linear system with $n_a$ equations and $n$ variables.

An example is a hierarchical time series where series at upper levels can be expressed by appropriately summing part or all of the series at the bottom level. <ref>(a) shows the two-level hierarchical structure for three linearly constrained time series such that $y_{T,t} = y_{X,t} + y_{Y,t}$, $\forall t = 1,...,T$. Now let $\yvet_t = [\uvet_t',~ \bvet_t']'$, where $\uvet_t = [y_{1,t}, \dots, y_{n_a,t}]'$ is the $n_a$-vector of upper levels time series and $\bvet_t = [y_{(n_a+1),t}, \dots, y_{n,t}]'$ is the $n_b$-vector of bottom level time series with $n = n_a+n_b$. The upper and lower level time series are connected by the cross-sectional aggregation matrix $\Avet_{cs}$ such that $\uvet_t = \Avet_{cs}\bvet_t$. Following <cit.>, we can always construct a zero-constraints cross-sectional matrix from the aggregation matrix, $\Cvet_{cs}=[\Ivet_{n_a} ~~ {-\Avet_{cs}}]$ . Finally, the cross-sectional structural matrix is given by $\Svet_{cs} = [\Avet' ~~ \Ivet_{n_b}]'$, providing the structural representation <cit.> $\yvet_t = \Svet_{cs} \bvet_t$. Considering the hierarchical example in <ref>(a), we have $\Avet_{cs} = \left[ 1 \quad 1 \right]$, $\Cvet_{cs} = \left[ 1 \quad -1\quad -1 \right]$ and $\Svet_{cs} = \left[\begin{array}{cc} 1 & 1 \\[-0.25cm] \multicolumn{2}{c}{\Ivet_2} \end{array}\right]$. In general there is no reason for $\uvet_t$ to be restricted to simple sums of $\bvet_t$; therefore $\Avet_{cs} \in \mathbb{R}^{n_a\times n_b}$ may contain any real values, and not only 0s and 1s.

Considering now the temporal framework, we denote as $\mathcal{K} = \{ k_p , k_{p-1}, \dots, k_2, k_1 \}$ the set of $p$ factors of $m$, in descending order, where $k_1= 1$ and $k_p= m$ <cit.>. Given a factor $k$ of $m$, and assuming that $T = N m$ (where $N$ is the length of the most temporally aggregated version of the series), we can construct a temporally aggregated version of the time series of a single variable $\{\yvet_{i,t}\}_{t = 1, \dots, T}$, through the non-overlapping sums of its $k$ successive values, which has a seasonal period equal to $M_k= \displaystyle\frac{m}{k}$: $x_{i,j}^{[k]} = \displaystyle\sum_{t=(j-1)k+1}^{jk} y_{i,t}$, where $j = 1,\dots, N_k$, $i = 1,\dots,n$, $N_k = \displaystyle\frac{T}{k}$ and $x_{i,j}^{[1]}=y_{i,t}$. Define $\tau$ as the observation index of the most aggregate level $k_p$. For a fixed temporal aggregation order $k \in \mathcal{K}$, we stack the observations in the column vector $\xvet_{i,\tau}^{[k]} = \left[x_{i,M_k(\tau-1)+1}^{[k]} \quad x_{i,M_k(\tau-1)+2}^{[k]} \quad \dots \quad x_{i,M_k\tau}^{[k]}\right]',$ and obtain the vector for all the temporal aggregation orders $\xvet_{i,\tau} = \left[x_{i,\tau}^{[k_p]} \quad \xvet_{i,\tau}^{[k_{p-1}]\prime} \quad \dots \quad \xvet_{i,\tau}^{[1]\prime} \right]'$, $\tau = 1,\dots,N$. The structural representation of the temporal hierarchy <cit.> is then $\xvet_{i,\tau} = \Svet_{te}\xvet_{i,\tau}^{[1]}$, where $\Svet_{te} = [\Avet_{te}' ~~ \Ivet_{m}]'$ is the $[(m+k^\ast) \times m]$ temporal structural matrix,

\Avet_{te} = \begin{bmatrix}
\multicolumn{3}{c}{\Unovet_{k_p}'}                       \\
\Ivet_{\frac{m}{k_{p-1}}} & \otimes & \Unovet_{k_{p-1}}' \\
& \vdots  &                    \\
\Ivet_{\frac{m}{k_{2}}}   & \otimes & \Unovet_{k_2}'
\end{bmatrix}

is the $(k^\ast \times m)$ temporal aggregation matrix with $k^\ast = \displaystyle\sum_{k \in \mathcal{K}\setminus\{k_1\}} M_k$, and $\otimes$ is the Kronecker product. For each series $x_{i,\tau}$, $i = 1,\dots,n$, we have also the zero-constrained representation

\begin{equation}
\label{eq:te_con}
\Cvet_{te}\xvet_{i,\tau} = \Zerovet_{[k^\ast \times (m+k^\ast)]}, \qquad \tau = 1,\dots,N, \qquad i = 1,\dots, n
\end{equation}

where $\Cvet_{te} = [\Ivet_{k^\ast} ~~ {-\Avet_{te}}]$ is the $[k^\ast \times (m+k^\ast)]$ zero constraints temporal matrix. <ref>(b) shows the hierarchical representation of a quarterly time series, for which $m = 4$, $\mathcal{K} = \{4,2,1\}$, and

\Avet_{te} = \begin{bmatrix}
1 & 1 & 1 & 1 \\
1 & 1 & 0 & 0 \\
0 & 0 & 1 & 1 \\
\end{bmatrix}, \quad \Cvet_{te} = \begin{bmatrix}
1 & 0 & 0 & -1 & -1 & -1 & -1 \\
0 & 1 & 0 & -1 & -1 & 0  & 0  \\
0 & 0 & 1 & 0  & 0  & -1 & -1 \\
\end{bmatrix} \quad \mathrm{and} \quad \Svet_{te} = \begin{bmatrix}
\Avet_{te} \\
\Ivet_4
\end{bmatrix}.

When we temporally aggregate each series, the cross-sectional constraints for the most temporally disaggregated series (<ref>) hold for all the temporal aggregation orders such that $\Cvet_{cs}\xvet^{[k]}_j = \Zerovet_{(n_a \times 1)}$, for $k \in \mathcal{K}$ and $j = 1, \dots, N_k$, where $\xvet_j^{[k]} = \left[\uvet_j^{[k]\prime}\quad \bvet_j^{[k]\prime}\right]'$ with $\uvet^{[k]}_j = \left[ x^{[k]}_{1,\;j}\quad \dots\quad x^{[k]}_{n_a,\;j}\right]'$ is the $n_a$-vector of upper time series and $\bvet^{[k]}_j = \left[x^{[k]}_{(n_a+1),\;j}\quad\dots\quad x^{[k]}_{n,\;j}\right]'$ is the $n_b$-vector of bottom time series in the temporal hierarchy.

To include both cross-sectional and temporal constraints at the same time in a unified framework, we stack the series into a $[n \times (m+k^\ast)]$ matrix $\Xvet_\tau$, whose rows and columns represent, respectively, the cross-sectional and the temporal dimension:

\begin{equation}
\label{eq:Xtau}
\Xvet_\tau = \begin{bmatrix}
\xvet_{1,\tau}' \\
\vdots          \\
\xvet_{n,\tau}'
\end{bmatrix} = \begin{bmatrix}
\Xvet_{\tau}^{[k_p]} & \dots & \Xvet_{\tau}^{[k_1]} \\ \end{bmatrix}
\quad \text{with} \quad \Xvet_{\tau}^{[k]} = \begin{bmatrix}
\Uvet_{\tau}^{[k]} \\
\Bvet_{\tau}^{[k]},
\end{bmatrix},
\end{equation}

where for any fixed $k$,
$\Uvet_{\tau}^{[k]}$ is the ($n_a\times N_k$) matrix grouping the upper time series, $\Bvet_{\tau}^{[k]}$ is the ($n_b\times N_k$) matrix grouping the bottom time series. Further, $\Cvet_{cs}\Xvet_\tau = \Zerovet_{\left[n_a \times (m+k^\ast)\right]}$ and $\Cvet_{te}\Xvet_\tau' = \Zerovet_{(k^\ast \times n)} $. We can consider the cross-temporal framework as a generalization of the cross-sectional and temporal frameworks, that simultaneously takes into account both types of constraints. The cross-sectional reconciliation approach proposed by <cit.> can be obtained by assuming $m = 1$, while the temporal one <cit.> is obtained when $n = 1$ (with $n_a = 0$ and $n_b = 1$).

<cit.> show that the cross-temporal constraints working on the complete set of observations corresponding to time period $\tau$ can be expressed in a zero-constrained representation through the full rank $\left[(n_am+nk^\ast)\times n(m+k^\ast)\right]$ zero constraints cross-temporal matrix $\Cvet_{ct}$ such that

\begin{equation}
\label{eq:Cct}
\Cvet_{ct} = \begin{bmatrix}
\Cvet_\ast \\
\Ivet_n \otimes \Cvet_{te}
\end{bmatrix} \quad \Longrightarrow \quad
\Cvet_{ct} \xvet_{\tau} = \Zerovet_{[(n_am+nk^\ast)\times1]} \quad \mathrm{for} \quad \tau = 1,\dots,N,
\end{equation}

where $\xvet_{\tau} = \mathrm{vec}(\Xvet_{\tau}') = [\xvet_{1, \tau}',~ 	\dots, ~ \xvet_{n, \tau}']'$, $\Cvet_\ast = [\Zerovet_{(n_a m\times nk^\ast)} ~~ \Ivet_m \otimes \Cvet_{cs}]\Pvet'$, and $\Pvet$ is the commutation matrix <cit.> such that $\Pvet \mathrm{vec}(\Yvet_{\tau}) = \mathrm{vec}(\Yvet_{\tau}')$. A structural representation can be considered as well: $\xvet_\tau = \Svet_{ct}\bvet^{[1]}_\tau = s(\bvet_{\tau}^{[1]})$, where

\begin{equation}
\label{eq:Sct}
\Svet_{ct} = \Svet_{cs} \otimes \Svet_{te}
\end{equation}

is the $\left[n(k^\ast+m)\times n_b m\right]$ cross-temporal summation matrix, $s: \mathbb{R}^{n_b m} \rightarrow \mathbb{R}^{n(m+k^\ast)}$ is the operator describing the pre-multiplication by $\Svet_{ct}$, and $\bvet^{[1]}_\tau = \mathrm{vec}(\Bvet^{[1]\prime}_{\tau})$. We observe that, in agreement with <cit.>, $\xvet_{\tau}$ lies in an $(n_b m)$-dimensional subspace $\mathfrak{s}_{ct}$ of $\mathbb{R}^{n(k^\ast+m)}$, which we refer to as the cross-temporal coherent subspace, spanned by the columns of $\Svet_{ct}$. In <ref>, we have represented $\Svet_{ct}$ for 3 linearly constrained quarterly time series, shown in <ref>.

§.§ Optimal point forecast reconciliation

Let $\widehat{\xvet}_{h} = \mathrm{vec}(\widehat{\Xvet}_{h}')$, $h = 1, \dots, H$, be the $h$-step ahead base forecasts (however obtained) with error covariance matrix given by $\bm{W}_h = \text{Var}(\widehat{\xvet}_h - \xvet)$, where $H$ is the forecast horizon for the most temporally aggregated time series. Denote

\widehat{\Xvet}_{h} = \begin{bmatrix}
\widehat{\xvet}_{1,h} \\
\vdots                \\
\widehat{\xvet}_{n,h}
\end{bmatrix} =\begin{bmatrix}
\widehat{\Uvet}_{h}^{[m]} & \dots & \widehat{\Uvet}_{h}^{[k]} & \dots & \widehat{\Uvet}_{h}^{[1]} \\[0.25cm]
\widehat{\Bvet}_{h}^{[m]} & \dots & \widehat{\Bvet}_{h}^{[k]} & \dots & \widehat{\Bvet}_{h}^{[1]} \\\end{bmatrix},

where $\widehat{\Uvet}_{h}^{[k]}$ is the ($n_a\times M_k$) matrix grouping the upper time series and $\widehat{\Bvet}_{h}^{[k]}$ is the ($n_b\times M_k$) matrix grouping the bottom time series for a given temporal aggregation order $k$. The matrix $\widehat{\Xvet}_{h}$, organized as ${\Xvet}_{\tau}$ in expression (<ref>), contains incoherent forecasts, such as $\Cvet_{ct} \widehat{\xvet}_{h} \neq \Zerovet_{[(n_am+nk^\ast)\times1]}$
with $h = 1, \dots, H$ and $\widehat{\xvet}_{h} = \mathrm{vec}(\widehat{\Xvet}_{h}')$. In this framework, the definition for forecast reconciliation in the cross-sectional framework given by <cit.> can be generalized as follows.

Forecast reconciliation aims to adjust the base forecast $\widehat{\xvet}_{h}$ by finding a mapping $\psi: \mathbb{R}^{n(m+k^\ast)} \rightarrow \mathfrak{s}$ such that $\widetilde{\xvet}_{h} = \psi\left(\widehat{\xvet}_{h}\right)$, where $\widetilde{\xvet}_{h} \in \mathfrak{s}$ is the vector of the reconciled forecasts.

For a given forecast horizon $h = 1,\dots, H$, the mapping $\psi$ may be defined as a projection onto $\mathfrak{s}$ given by <cit.>

\begin{equation}
\label{eq:Mvet}
\widetilde{\xvet}_{h} = \psi\left(\widehat{\xvet}_h\right) = \Mvet \widehat{\xvet_h},
\end{equation}

where $\Mvet = \Ivet_{n(m+ k^\ast)} - \Omegavet_{ct}\Cvet'_{ct}\left(\Cvet_{ct}\Omegavet_{ct}\Cvet'_{ct}\right)^{-1}\Cvet_{ct}$, for a positive definite matrix $\Omegavet_{ct}$, and $\widetilde{\xvet}_{h} = \mathrm{vec}(\widetilde{\Xvet}'_{h})$.
According to <cit.> showed that the minimum variance linear unbiased reconciled forecasts, satisfying the unbiased condition $\text{E}(\widetilde{\xvet}_h -\bm{\xvet}_h) = 0$, has solution (<ref>) when $\Omegavet_{ct} = \Wvet_h$.

Alternatively, the cross-temporal reconciled forecasts $\widetilde{\Xvet}_{h}$ may be found according to the structural approach proposed by <cit.> for the cross-sectional framework, yielding $\widetilde{\xvet}_h = \Svet_{ct}\Gvet \widehat{\xvet}_h$ for some matrix $\Gvet$. <cit.> showed that this leads to a solution equivalent to the cross-temporally reconciled forecasts in (<ref>), given by

\begin{equation}\label{eq:SGy}
\widetilde{\xvet}_{h} = \psi\left(\widehat{\xvet}_h \right) = \left(s \circ g \right)\left(\widehat{\xvet}_h\right)=\Svet_{ct}\Gvet \widehat{\xvet}_{h},
\end{equation}

where $\Gvet = (\Svet_{ct}' \Omegavet_{ct}^{-1}\Svet_{ct})^{-1} \Svet_{ct}'\Omegavet_{ct}^{-1}$,  and $\Mvet = \Svet_{ct} \Gvet$. In this case, $\psi$ is the composition of two transformations, say $s \circ g$, where $g: \mathbb{R}^{n(m+k^\ast)} \rightarrow \mathbb{R}^{n_b m}$ is a continuous function. In <ref> we report some cross-sectional, temporal and cross-temporal approximations for the covariance matrix to be used in (<ref>) and (<ref>).

§.§ Cross-temporal bottom-up forecast reconciliation

The classic bottom-up approach <cit.> simply consists in summing-up the base forecasts of the most disaggregated level in the hierarchy to obtain forecasts of the upper-level series. To reduce the computational cost involved in optimal cross-temporal reconciliation, we may be interested in applying a reconciliation along only one dimension (cross-sectional or temporal) and reconstructing the cross-temporal structure using a partly bottom-up approach <cit.>.

<ref> provides a visual representation of partly bottom-up in a two-step cross-temporal reconciliation approach. On the left (<ref>), we first compute the cross-sectionally reconciled forecasts at the highest frequency ($k = 1$), and then apply temporal bottom-up to obtain coherent cross-temporal forecasts. On the right (<ref>), we first compute temporally reconciled forecasts for the most disaggregated cross-sectional level, and then apply the cross-sectional bottom-up. We denote these two-step reconciliation approaches, respectively, as ct$(rec_{te},bu_{cs})$, and ct$(rec_{cs},bu_{te})$, where ‘$rec_{te}$’ and ‘$rec_{cs}$’ denote a forecast reconciliation approach in the temporal and cross-sectional dimensions and, ‘$bu_{cs}$’ and ‘$bu_{te}$’ denote using bottom-up in the cross-sectional and temporal dimensions, respectively. It is worth noting that the simple cross-temporal bottom-up approach corresponds to $\mathrm{ct}(bu_{cs}, bu_{te})=\mathrm{ct}(bu_{te}, bu_{cs})=\mathrm{ct}(bu)$.

$\widetilde{\Xvet}$ with ct$(rec_{cs}, bu_{te})$
[>=latex, line width=1pt,
matrix of nodes,
text width = 1.5cm,
text height = 0.65cm,
column sep=2pt,
row sep=7pt,
nodes in empty cells,
left delimiter=[,
right delimiter=],
ampersand replacement=&
[Matrix] (Mt)
$\widetilde{\Uvet}_{te(bu)}^{[m]}$ & …& $\widetilde{\Uvet}_{te(bu)}^{[k_2]}$ & $\widetilde{\Uvet}_{cs(rec)}^{[1]}$

$\widetilde{\Bvet}^{[m]}_{te(bu)}$ & …& $\widetilde{\Bvet}^{[k_2]}_{te(bu)}$ & $\widetilde{\Bvet}^{[1]}_{cs(rec)}$

[<-, opacity = 0] (Mt.north east)++(0.4,0) coordinate (temp) – (temp |- Mt.south) node [midway,label=[label distance=0.1cm,rotate=-90, xshift = 1.5mm, font=]Cross-sectional];
[<-] (Mt.south west)++(0,-0.15) coordinate (temp) – (temp -| Mt.east) node [midway,label=[label distance=0cm,xshift = 1.5mm, font=]below:Temporal];
rounded corners,
inner sep=0pt, fill = blue, fit=(Mt-1-4)(Mt-2-4)](Bt);
rounded corners,
inner sep=0pt, fill = red, fit=(Mt-1-1)(Mt-2-3)](At);

$\widetilde{\Xvet}$ with ct$(rec_{te}, bu_{cs})$
[>=latex, line width=1pt,
matrix of nodes,
text width = 1.5cm,
text height = 0.65cm,
column sep=2pt,
row sep=7pt,
nodes in empty cells,
left delimiter=[,
right delimiter=],
ampersand replacement=&
[Matrix] (Mcs)
$\widetilde{\Uvet}_{cs(bu)}^{[m]}$ & …& $\widetilde{\Uvet}_{cs(bu)}^{[k_2]}$ & $\widetilde{\Uvet}_{cs(bu)}^{[1]}$

$\widetilde{\Bvet}^{[m]}_{te(rec)}$ & …& $\widetilde{\Bvet}^{[k_2]}_{te(rec)}$ & $\widetilde{\Bvet}^{[1]}_{te(rec)}$

[<-] (Mcs.north east)++(0.4,0) coordinate (temp) – (temp |- Mcs.south) node [midway,label=[label distance=0.1cm,rotate=-90, xshift = 1.5mm, font=]Cross-sectional];
[<-, opacity = 0] (Mcs.south west)++(0,-0.15) coordinate (temp) – (temp -| Mcs.east) node [midway,label=[label distance=0cm,xshift = 1.5mm, font=]below:Temporal];
rounded corners,
inner sep=0pt, fill = blue, fit=(Mcs-2-1)(Mcs-2-4)](Bcs);
rounded corners,
inner sep=0pt, fill = red, fit=(Mcs-1-1)(Mcs-1-4)](Acs);

A visual representation of partly bottom up starting from (<ref>) cross-sectionally reconciled forecasts for the temporal order 1 ($\widetilde{\Uvet}^{[1]}$ and $\widetilde{\Bvet}^{[1]}$) followed by temporal bottom-up, and (<ref>) temporally reconciled forecasts of the cross-sectional bottom time series $(\widetilde{\Bvet}^{[k]}, \, k\in \mathcal{K})$ followed by cross-sectional bottom-up. The blue background indicates generating reconciled forecasts along one dimension, while the pink background indicates the forecasts obtained using bottom-up along the other.

§ PROBABILISTIC FORECAST RECONCILIATION

To introduce the idea of coherence and probabilistic forecast reconciliation, we adapt the notations and the formal definitions introduced in <cit.> and <cit.> for the cross-sectional probabilistic case. These definitions can also be generalized to the cross-temporal framework by following the approach developed by <cit.> for count data. However, in this paper we only focus on the continuous case.

Our aim is to extend these definitions to cross-temporal coherent probabilistic forecasts and cross-temporal probabilistic forecast reconciliation. Let $(\mathbb{R}^{n_b m}, \mathcal{F}_{\mathbb{R}^{n_b m}}, \nu)$ be a probability space for the bottom time series $\bvet_{\tau}^{[1]}$, where $\mathcal{F}_{\mathbb{R}^{n_b m}}$ is the Borel $\sigma$-algebra on $\mathbb{R}^{n_b m}$. Then a $\sigma$-algebra $\mathcal{F}_{\mathfrak{s}}$ can be constructed from the collection of sets $s(\mathcal{B})$ for all $\mathcal{B} \in \mathcal{F}_{\mathbb{R}^{n_b m}}$.

Given the probability space $(\mathbb{R}^{n_b m}, \mathcal{F}_{\mathbb{R}^{n_b m}}, \nu)$, we define the coherent probability space as the triple $(\mathfrak{s}, \mathcal{F}_{\mathfrak{s}}, \breve{\nu})$ satisfying the following property:
$\breve{\nu}(s(\mathcal{B}))=\nu(\mathcal{B})$, $\forall \mathcal{B} \in \mathcal{F}_{\mathbb{R}^{n_b m}}$.

Let $(\mathbb{R}^{n(m+k^\ast)}, \mathcal{F}_{\mathbb{R}^{n(m+k^\ast)}}, \hat{\nu})$ be a probability space referring to the incoherent probabilistic forecast ($\widehat{\xvet}_{h}$) for all the $n$ series in the system at any temporal aggregation order $k \in \mathcal{K}$.

The reconciled probability measure of $\hat{\nu}$ with respect to $\psi$ is a probability measure $\tilde{\nu}$ on $\mathfrak{s}$ with $\sigma$-algebra $\mathcal{F}_{\mathfrak{s}}$ satisfying

\begin{equation}\label{eq:pfr}
\tilde{\nu}(\mathcal{A})=\hat{\nu}(\psi^{-1}(\mathcal{A})), \quad \forall \mathcal{A} \in \mathcal{F}_{\mathfrak{s}},
\end{equation}

where $\psi^{-1}(\mathcal{A})=\{x \in \mathbb{R}^{n(m+k^\ast)}: \psi(x) \in \mathcal{A}\}$ denotes the pre-image of $\mathcal{A}$.

The map $\psi$ may be obtained as the composition $s \circ g$, as for the cross-temporal point reconciliation (<ref>). With the following result, we can separate the mechanism used to generate the base forecasts samples from the reconciliation phase.

Suppose that $(\widehat{\xvet}_1, \dots, \widehat{\xvet}_L)$ is a sample drawn from a (cross-temporal) incoherent probability measure $\widehat{\nu}$. Then $(\widetilde{\xvet}_1, \dots, \widetilde{\xvet}_L)$, where $\widetilde{\xvet}_\ell=\psi(\widehat{\xvet}_\ell)$ and $\ell= 1, \dots, L$, is a sample drawn from the (cross-temporal) reconciled probability measure $\widetilde{\nu}$ defined in (<ref>).

See Theorem 4.5 from <cit.> using Definition <ref>.

Theorem <ref> is the cross-temporal extension of Theorem 4.5 in <cit.>. It means that a sample from the reconciled distribution can be obtained by reconciling each member of a sample from the incoherent distribution.

§.§ Parametric framework: Gaussian reconciliation

It is possible to obtain a reconciled probabilistic forecast analytically for some parametric distributions, such as the multivariate normal <cit.>. In the cross-sectional framework, <cit.> show that, starting from an elliptical distribution for the base forecasts, the reconciled forecast distribution is also elliptical. Using the results shown in <ref>, we may extend[We assume $H =1$ and simplify the notation by removing the $h$ suffix without loss of generality] this results to the cross-temporal case. To obtain a reconciled forecast using the multivariate normal distribution, we start with a base forecast distributed as $\mathcal{N}(\widehat{\xvet}, \Sigmavet)$, where $\widehat{\xvet}_h$ is the mean vector and $\Sigmavet$ is the covariance matrix of the base forecasts. The reconciled forecast distribution is then given by $\mathcal{N}(\widetilde{\xvet}, \widetilde{\Omegavet})$, where

\begin{equation}\label{eq:meanvar}
\widetilde{\xvet} = \Mvet\widehat{\xvet} \quad \mbox{and} \quad \widetilde{\Omegavet} = \Mvet \Sigmavet \Mvet',
\end{equation}

where $\Mvet$ is the projection matrix defined in (<ref>).
Note that if we assume that $\Sigmavet = \Omegavet_{ct}$, then the covariance matrix in (<ref>) simplifies to $\widetilde{\Omegavet} = \Mvet \Omegavet_{ct}$. In the cross-temporal case, sensibly estimating the covariance matrix $\Sigmavet$ can be difficult because we need to simultaneously consider both the temporal and cross-sectional structures. This requires many parameters to be estimated, which can be challenging in practice. Additionally, naively using one-step residuals to estimate the cross-temporal correlation structure can lead to an inappropriate estimate of the covariance matrix[In particular, some temporal covariances are fixed to zero (see the online appendix for more details).]. These challenges will be explored in more depth in the following sections.

Focusing on the computational aspect[We use two R packages to sample from a the base forecast gaussian distribution:  <cit.> when the dimension are very big (<ref>) and  <cit.> otherwise  (Sections <ref> and <ref>).], we can take several steps to reduce the time required to obtain simulations from the reconciled forecast distribution. For example when dealing with a genuine hierarchical structure (that share the same top- and bottom-level variables, ), it is not necessary to simulate from a normal distribution with a defined covariance matrix for the entire structure. Instead, we can utilize the properties of elliptical distributions to simulate from the high frequency bottom time series and then obtain the complete simulation through the $\Svet_{ct}$ matrix. Furthermore, we do not need to calculate the reconciled mean and variance and generate a new sample if we already have a sample from the normal distribution of the base forecasts; we can simply apply the point forecast reconciliation (<ref>) as outlined in Theorem <ref>. The relationships between base and reconciled forecast distributions and their respective simulations through Theorem <ref> are depicted in <ref>.

§.§ Non-parametric framework: bootstrap reconciliation

Analytical expressions for the base and reconciled forecast distributions are sometimes challenging to express, and sometimes unrealistic parametric assumptions are used. We propose a procedure called cross-temporal joint (block) bootstrap (ctjb) to generate samples from the base forecast distributions that preserve cross-temporal relationships. This approach involves drawing samples of all series at the same time from the most temporally aggregated level, and using the most temporally aggregated level to determine the corresponding time indices for the other levels.

Let $\widehat{\Evet}^{[k]}$ be the ($n \times N_k$) matrix of the residuals for $k \in \mathcal{K}$. <ref> (on the left) provides a visualization of these matrices and how they are related to each other for the example in <ref>. It is assumed that the residuals cover four years ($N=4$): the green color corresponds to the first year, the blue to the second year, and so on. Further, let $\mathcal{M}_i$ be the model used to calculate the base forecasts and residuals for the $i^{th}$ series. In this work, we assume $\mathcal{M}_i$ to be a univariate model, however nothing prevents the use of multivariate models, perhaps for different temporal levels or for groups of time series.

Assuming $H = 1$, $\tau$ is a random draw with replacement from $1,\dots, N$ and the $\ell^{th}$ bootstrap incoherent sample is $\widehat{\xvet}_{i,\ell}^{[k]} = f_i(\mathcal{M}_i, \widehat{\evet}_{i}^{[k]})$, where $f_i(\cdot)$ depends on the fitted univariate model $\mathcal{M}_i$. That is, $\widehat{\xvet}_{i,l}^{[k]}$ is a sample path simulated for the $i^{th}$ series with error approximated by the corresponding block bootstrapped sample residual $\widehat{\evet}_{i}^{[k]}$, the $i^{th}$ row of

\widehat{\Evet}^{[k]}_{\tau} = \begin{bmatrix}
\widehat{e}^{[k]}_{1,M_k(\tau-1)+1} & \dots  & \widehat{e}^{[k]}_{1,M_k\tau}   \\
\vdots                              & \ddots & \vdots                          \\
\widehat{e}^{[k]}_{n,M_k(\tau-1)+1} & \dots  & \widehat{e}^{[k]}_{n,M_k\tau} \
\end{bmatrix}.

<ref> (on the right) shows $\widehat{\Evet}^{[k]}_{\tau}$, $k\in\{4,2,1\}$, for the quarterly cross-temporal hierarchy in <ref>.

One of the main advantages of the cross-temporal joint bootstrap is that it allows us to accurately account for the dependence between the different levels of temporal aggregation and not only the cross-sectional dependencies. By sampling residuals from the most temporally aggregated level and using it to determine the indices for the other levels, we can ensure that the bootstrap sample reflects the underlying data distribution. Additionally, the cross-temporal joint bootstrap is easy to implement in R <cit.> using the package  <cit.> for many forecasting models, making it a practical and efficient tool. Furthermore, this approach is easily scalable in order to utilize multiple computing power simultaneously for each individual series. This can be especially useful when dealing with large datasets or when trying to speed up the analysis process.

§ CROSS-TEMPORAL COVARIANCE MATRIX ESTIMATION

Representation of four types of covariance matrices that can be obtained from the cross-temporal hierarchical structure ($3$ time series and $m = 2$) for two different values of $\lambda\in\{0,1\}$, the shrinkage parameter. The cells that are not modified by shrinkage are colored black, those actively involved in the shrinkage phase are colored light blue, and those derived from and not estimated by the base forecasts errors are colored blue. Additionally, for $\lambda = 1$, the cells corresponding to a zero value are colored white.

As the covariance matrix $\Omegavet$ is unknown in practice, a natural estimate is the empirical sample covariance matrix of the base forecasts $\widehat{\Omegavet}$. In this section, our focus will be exclusively on the cross-temporal framework., this means that we have to estimate $r = n(k^\ast+m)[n(k^\ast+m)-1]/2$ different parameters. A possible solution to estimating many parameters when we have fewer observations than $r$, is to construct a shrinkage estimator <cit.>, using a convex combination of $\widehat{\Omegavet}$ and a diagonal target matrix $\widehat{\Omegavet}_D = \widehat{\Omegavet} \odot \Ivet_{n(k^\ast+m)}$, such that $\widehat{\Omegavet}_{G} = \lambda \widehat{\Omegavet}_D + (1-\lambda) \widehat{\Omegavet}$, where $\lambda \in [0,1]$ is the shrinkage intensity parameter that can be estimate using the unbiased estimator proposed by <cit.> (see ). The linear combination involving these two matrices is referred to as Global shrinkage (G), where all off-diagonal elements are shrunk towards zero. $\widehat{\Omegavet}_{G}$ corresponds to the matrix used by the reconciliation approach oct$(shr)$ shown in <ref>. However, shrinking all off-diagonal elements to zero, when we know that the covariance matrix has a cross-sectional and/or temporal structure, results in information loss. Therefore, we propose to estimate a smaller matrix, and to use the cross-sectional and/or temporal structure to obtain a better estimator for the covariance matrix of the entire system. Given that $\Svet_{ct} = \Svet_{cs} \otimes \Svet_{te}$, it is possible to express the actual covariance matrix in terms of three smaller matrices such that

\begin{equation}
\label{eq:OmSct}
\Omegavet = \Svet_{ct}\Omegavet_{\textit{hf-bts}}\Svet_{ct}' = \left(\Ivet_n \otimes \Svet_{te}\right)\Omegavet_{\textit{hf}}\left(\Ivet_n \otimes \Svet_{te}\right)' = \left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)\Omegavet_{bts}\left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)',
\end{equation}

where $\Omegavet_{\textit{hf-bts}}$ is the $(n_b m\times n_b m]$ covariance matrix for the bottom time series at temporal aggregation level $k = 1$ (highest frequency bottom time series), $\Omegavet_{\textit{hf}}$ is the $(nm\times nm)$ covariance matrix related to all the high frequency time series and $\Omegavet_{bts}$ is the $[n_b(k^\ast + m)\times n_b(k^\ast + m)]$ covariance matrix related to bottom time series at any temporal aggregation.

Therefore, we can apply the idea of “Stein-type shrinkage" <cit.> to $\Omegavet_{\textit{hf-bts}}$, $\Omegavet_{\textit{hf}}$ and $\Omegavet_{\textit{bts}}$ by using the corresponding empirical base forecasts residuals estimation. We obtain the following expressions (see the online appendix for details):

* High frequency Bottom time series shrinkage matrix (HB):

$\widehat{\Omegavet}_{HB} = \lambda \Svet_{ct}\widehat{\Omegavet}_{\textit{hf-bts}, D}\Svet_{ct}'+ (1-\lambda) \Svet_{ct}\widehat{\Omegavet}_{\textit{hf-bts}}\Svet_{ct}'$;

* High frequency shrinkage matrix (H):
$\widehat{\Omegavet}_{H}  = \lambda (\Ivet_{n} \otimes \Svet_{te})\widehat{\Omegavet}_{hf, D}(\Ivet_{n} \otimes \Svet_{te})' + (1-\lambda) (\Ivet_{n} \otimes \Svet_{te})\widehat{\Omegavet}_{\textit{hf}}(\Ivet_{n} \otimes \Svet_{te})'$;

* Bottom time series shrinkage matrix (B):
$\widehat{\Omegavet}_{B} = \lambda \left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)\widehat{\Omegavet}_{bts, D}\left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)' +  (1-\lambda) \left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)\widehat{\Omegavet}_{bts}\left(\Svet_{cs} \otimes \Ivet_{m+k^\ast}\right)'$,

where $\widehat{\Omegavet}_{l, D} = \Ivet_{n_b m}\odot\widehat{\Omegavet}_{j}$, $l = \{\textit{hf-bts}, \;\textit{hf}, \;\textit{bts}\}$, and $\lambda$ is the shrinkage parameter.

Method                # of different parameters                                 AR(2)     GDP     Tourism

G                              $r = \displaystyle\frac{n(k^\ast+m)[n(k^\ast+m)-1]}{2}$             36                 $221\,445$              $108\,052\,350$

HB     $r_{HB} = \displaystyle\frac{n_bm[n_bm-1]}{2}<r$           6
$(83\%)$                 $30\,876$
$(86\%)$               $6\,655\,776$

H                   $r_{HB}<\displaystyle\frac{nm[nm-1]}{2}<r$     15
$(58\%)$                 $72\,390$
$(67\%)$               $19\,848\,150$

B                  $r_{HB}<\displaystyle\frac{n_b(k^\ast+m)[n_b(k^\ast+m)-1]}{2}<r$     15
$(58\%)$                 $94\,395$
$(57\%)$               $36\,231\,328$

Number of different parameters that need to be estimated for the Monte Carlo simulation (AR(2), see <ref>), the Australian GDP (see <ref>) and the Australian Tourism Demand (see <ref>) forecasting experiments: the first one has $3$ time series (one upper and two bottom) with temporal aggregation $\mathcal{K} = \{2, 1\}$; the second one has $95$ quarterly ($m = 4$ and $k^\ast = 3$) time series ($62$ free and $33$ constraints, see ); the last one has a total of 525 monthly ($m = 12$ and $k^\ast = 16$) time series ($304$ bottom and $221$ upper). The percentage reductions in the number of parameters compared to the global approach are reported in parentheses.

<ref> gives some visual insights on the covariance matrices obtainable with $\lambda=0$ and $\lambda=1$, respectively, for a simple cross-temporal hierarchical structure with 3 time series and $\mathcal{K}=\{2,1\}$ (e.g, cross-temporal semi-annual, see the Monte Carlo simulation in <ref>). Another important aspect is the number of parameters to be estimated through the residuals of the base forecasts. In <ref> we report the number of different parameters for the Monte Carlo simulation (AR2) and for the two forecasting experiment: Australian GDP (see <ref>) and Australian Tourism Demand (see <ref>). In addition, we calculate also the percentage reductions in the number of parameters compared to the global approach, that is: $\% \text{ reduction} = 100(1-r_i/r_G)$ with $i \in \{\textit{HB}, H, B\}$. As we can see, G involves a considerable number of parameters compared to other procedures. We observe that HB leads to a decrease of around 85%, whereas the H and B approaches look as a compromise between the previous two. In general, as $m$ and $n$ increase (see the online appendix), using H requires the estimation of a smaller number of parameters.

In our simulations and forecasting experiments, we will be closely analyzing these different constructions with a dual purpose. First, we will use the full covariance matrix ($\lambda = 0$) of the base forecasts to obtain a base forecasts sample of the linearly constrained time series under Gaussianity assumption. Then, we will use the shrinkage versions as approximations of the covariance matrix to be used for reconciliation. This will allow us to better understand the properties and abilities of each parameterization.

§.§ Multi-step residuals

Model residuals may be used to estimate the covariance matrix in cross-temporal forecast reconciliation. In time series analysis, it is common to use residuals corresponding to one-step ahead forecasts, but because of the different temporal dimensions, we need residuals corresponding to different forecast horizons. Thus, we define multi-step residuals as $e_{i,h,j}^{[k]} = x_{i,j+h}^{[k]} - \widehat{x}_{i,j+h|j}^{[k]}$, where $i = 1,\dots,n$, $j = 1,\dots,N_k$ and $\widehat{x}_{i,j+h|t}^{[k]}$ is the $h$-step fitted value, calculated as the $h$-step-ahead forecast using data up to time $j$. In general, these residuals will be autocorrelated except when $h=1$.

Following <cit.>, we use a matrix organization of the residuals similar to the one for the base forecasts in <ref>. Specifically, let $N$ be the total number of observations for the most temporally aggregate time series. Then, the $N_k$-vectors of multi-step residuals for the temporal aggregation $k$ and the series $i$, $\evet_{i,h}^{[k]} =  \Big[e_{i,h,1}^{[k]} \quad e_{i,h,2}^{[k]} \quad \dots \quad e_{i,h,N_k}^{[k]}\Big]'$ with $h = 1,\dots, M_k$, can be organized in matrix form as

\Evet_i^{[k]} = \begin{bmatrix}
e_{i,1,1}^{[k]}                     & e_{i,2,2}^{[k]}                     & \dots & e_{i,M_k,M_k}^{[k]} \\
\vdots                            & \vdots                            &       & \vdots                  \\
e_{i,1,N_k - M_k + 1}^{[k]} & e_{i,2,N_k - M_k + 2}^{[k]} & \dots & e_{i,M_k,N_k}^{[k]}         \\
\end{bmatrix}.

Let $\Evet_i = \Big[\Evet_i^{[m]} \quad \Evet_i^{[k_p-1]} \quad \dots \quad \Evet_i^{[1]}  \Big]$. Then the $[N \times n(m+k^\ast)]$ cross-temporal residual matrix is given by $
\Evet = \Big[\Evet_1 \quad \Evet_2 \quad \dots \quad \Evet_n \Big]$.

§.§ Overlapping residuals

Another issue that arises in the case of cross-temporal reconciliation is the low number of available residuals, especially for the higher orders of temporal aggregation. A possible solution is to use residuals calculated using overlapping series by allowing the year to have a varying starting time. To better explain how to calculate overlapping residuals, assume we have a single series $\yvet = [y_1 \; y_2 \; y_3 \; \dots\; y_{T-1}\; y_{T}]'$. We can construct $k$ non overlapping series such that $\xvet^{[k], s} = \left\{x^{[k],s}_{j}\right\}_{j = 1}^{N_k-s}$ where $x^{[k],s}_{j} = \displaystyle\sum_{t = (j-1)k+s+1}^{jk-s} y_t$, with $s = 0, \dots, (k-1)$. For example, suppose we have a biannual series with $k = 2$ and $T = 6$, then we can construct two annual time series depending on which time is deemed the start of the year: $\xvet^{[2], 0} =  \Big[x_1^{[2], 0}\quad x_2^{[2], 0}\quad x_{3}^{[2], 0} \Big]' =\Big[y_1 + y_2\quad y_3 + y_4\quad y_5 + y_6\Big]'$ and $\xvet^{[2], 1} = \Big[x_1^{[2], 1}\quad x_2^{[2], 1} \Big]' = \Big[y_2 + y_3\quad  y_4 + y_5 \Big]'$. To calculate overlapping residuals, we propose the following steps:

* Fit a model to $\xvet^{[k], 0}$ (i.e., select an appropriate model and estimate the model parameters using the available data) and calculate the residuals.

* Apply the same model in step 1 to $\xvet^{[k], s}$ for $s = 1, \dots, k-1$, without re-estimating the parameters, and calculate the residuals.

The resulting residuals can be used to estimate the covariance matrix in cross-temporal forecast reconciliation. This increases the number of available residuals, particularly when working with higher frequency observations such as monthly or daily data.

It is important to note that this approach assumes that the model used in step 1 is appropriate for all the different series $\xvet^{[k], s}$. Some seasonal models will not be appropriate as the seasonal pattern will be shifted for different values of $s$. However, this will not affect seasonal ARIMA models as the seasonality is defined in terms of lags which are unaffected by the value of $s$.

§ MONTE CARLO SIMULATION

We study the effect of combining cross-sectional and temporal aggregations, using a simple hierarchy that allows us to effectively visualize the quantities involved, such as the covariance matrix. Additionally, the small size and nature of the data generating process make it possible to exactly calculate the true cross-temporal covariance structure, thus providing insights into the nature of the time series data involved in the forecast reconciliation process.

Consider a $2$-level hierarchical structure with three time series (one upper series, $A$, and two bottom series, $B$ and $C$) such that the cross-sectional aggregation matrix is $\Avet_{cs} = \left[ 1 \quad 1 \right]$ ($A = B+C$). The temporal structure we are considering is equivalent to using semi-annual data with $\mathcal{K} = \{2,1\}$ and $m = 2$. The assumed Data-Generating Processes (DPG) for the semi-annual bottom level series are two AR(2) given by

\begin{aligned}
y_{B,t} &= \phi_{B, 1} y_{B,t-1} + \phi_{B, 2} y_{B,t-2} + \varepsilon_{B, t}\\
y_{C,t} &= \phi_{C, 1} y_{C,t-1} + \phi_{C, 2} y_{C,t-2} + \varepsilon_{C, t}
\end{aligned}

with parameters[The $\phivet_B$ and $\phivet_C$ parameters are estimated from the “Lynx" and “Hare" time series contained in the  dataset of the  package for R <cit.>.] $\phivet_B = [\phi_{B,1}\; \phi_{B,2}]' = [1.34\; -0.74]'$ and $\phivet_C  = [\phi_{C,1}\; \phi_{C,2}]' = [0.95\;~-~0.42]'$. The error $\epsvet_t = \left[\varepsilon_{B, t}\quad \varepsilon_{C, t}\right]'$ driving the process is drawn from a multivariate normal distribution with standard deviations simulated from a uniform distribution between 0.5 and 2 and a fixed correlation of -0.8. The cross-sectional error covariance matrix is thus given by

\Omegavet_{cs} = \begin{bmatrix}
0.9 & 0   \\
0   & 1.8
\end{bmatrix} \begin{bmatrix}
1    & \rho \\
\rho & 1
\end{bmatrix} \begin{bmatrix}
0.9 & 0   \\
0   & 1.8
\end{bmatrix} = \begin{bmatrix}
\sigma_B^2  & \sigma_{BC} \\
\sigma_{BC} & \sigma_C^2
\end{bmatrix}.

To obtain the remaining series, the bottom series are then cross-temporally aggregated.

For the forecast experiment, the base forecasts are computing using AR models where the order is automatically determined by the algorithm proposed by <cit.> and implemented in the R package  <cit.>, thus allowing for possible mis-specification in the models. The training window length is 500 years, consisting of 1000 high frequency observations. The experiment is replicated 500 times, with a forecast horizon of 1 year.

Since the AR(2) models used as DPG for the bottom series $B$ and $C$ at the most disaggregated temporal level are known, we may compute the true covariance matrix for one-step ahead forecasts at the annual level $\Omegavet_{ct} = \Svet_{ct}\Omegavet_{\textit{hf-bts}}\Svet_{ct}'$, where

\Omegavet_{\textit{hf-bts}} = \begin{bmatrix}
\sigma^2_B            &                                                 &                      &                                        \\
\phi_{B,1}\sigma_B^2  & \sigma_B^2\left(1+\phi_{B,1}^2\right)           &                      &                                        \\
\sigma_{BC}           & \phi_{B,1}\sigma_{BC}                           & \sigma_C^2           &                                        \\
\phi_{C,1}\sigma_{BC} & \sigma_{BC}\left(1+\phi_{B,1}\phi_{C,1} \right) & \phi_{C,1}\sigma_C^2 & \sigma_C^2\left(1+\phi_{C,1}^2\right)\
\end{bmatrix}.

The detailed calculations can be found in the online appendix.
<ref> shows both the covariance matrix and the correlation matrix for fixed parameters.

Simulation experiment. True cross-temporal covariance (left) and correlation (right) error matrix of the reconciled forecasts with $\sigma_B = 0.9$, $\sigma_C = 1.8$, $\phivet_B = [1.34\; -0.74]'$, $\phivet_C = [0.95\; -0.42]'$ and $\rho = -0.8$.

To construct cross-temporal samples of the base forecasts, we use the Gaussian and bootstrap approaches discussed in Sections <ref> and <ref>, respectively. For the parametric approach we use multi-step residuals with the different covariance matrix structures analyzed in <ref>, while for the non-parametric approach, we use regular one-step residuals. We do not use overlapping residuals in our analysis as we have the advantage of generating a large number of observation. Ten different reconciliation approaches have been adopted (see <ref>): ct$(bu)$, ct$(shr_{cs}, bu_{te})$, ct$(wlsv_{te}, bu_{cs})$, oct$(wlsv)$, oct$(bdshr)$, oct$_h(shr)$, oct$_h(bshr)$, oct$_h(hshr)$ and oct$_h(hbshr)$.

Label     Description

ct$(bu)$     Simple cross-temporal bottom-up (<ref>).

ct$(\;\cdot\;, bu_{te})$     Partly bottom-up (<ref>) starting from cross-sectional reconciled forecasts using the $shr$ and $wls$ approaches (<ref>).

ct$(wlsv_{te}, bu_{cs})$     Partly bottom-up (<ref>) starting from temporally reconciled forecasts using the $wlsv$ approach (<ref>).

oct$(\;\cdot\;)$     Optimal cross-temporal reconciliation for the $ols$, $struc$, $wlsv$ and $bdshr$ approaches (see <ref>). One-step residuals were used with $wlsv$ and $bdshr$.

oct$_h(\;\cdot\;)$     Optimal cross-temporal reconciliation with multi-step residuals (see <ref>) for the shrinkage approaches presented in <ref>: $shr$ stands for Global shrinkage, $hshr$ for High frequency shrinkage, $bshr$ for bottom time series shrinkage, $hbshr$ for High frequency bottom time series shrinkage.

oct$_o(\;\cdot\;)$     Optimal cross-temporal reconciliation with overlapping residuals (see <ref>) for the $wlsv$ and $bdshr$ approaches (<ref>).

oct$_{oh}(\;\cdot\;)$     Optimal cross-temporal reconciliation with overlapping and multi-step residuals (see Section <ref> and <ref>) for the shrinkage approaches presented in <ref>: $shr$ stands for Global shrinkage, $hshr$ for High frequency shrinkage.

Cross-temporal reconciliation approaches for the Monte Carlo simulation (see <ref>), the Australian GDP (see <ref>) and the Australian Tourism Demand (see <ref>) forecasting experiments. All the reconciliation procedures are available in the R package  <cit.>.

§.§ Covariance matrix comparison and forecast accuracy scores

To compare the true covariance matrix $\Omegavet_{ct}$ with the estimated covariance matrix $\widehat{\Omegavet}$, we use the Frobenius norm to quantify the difference between two matrices: $\lVert \Dvet \rVert_F = \displaystyle\sqrt{\sum_{i = 1}^{n(k^\ast + m)}\sum_{j = 1}^{n(k^\ast + m)}|d_{i,j}|^2}$ where $\Dvet = \widehat{\Omegavet} - \Omegavet_{ct}$. The true covariance matrix, shown in <ref>, was compared to the estimated covariance matrices obtained using various reconciliation approaches and techniques for generating sample paths of the base forecasts. Thus, we should be able to determine which reconciliation approach and simulation technique produce an accurate estimate of the covariance matrix. Other types of matrix norms were also considered with similar results.

From <ref>, it appears that the reconciled covariance matrices are always closer to the true matrix than the base forecast matrix when using both the Gaussian and the bootstrap  approach. Overall, there are no major differences in the findings when using either one-step or multi-step residuals in cross-temporal forecast reconciliation. In fact, using approaches like oct$(bdshr)$, we obtain results that are consistent with approaches such as oct$_h(shr)$, where no temporal and/or cross-sectional correlation assumptions are imposed. It is worth noting that the $HB$ covariance matrix when used to calculate the base forecasts samples, is not changed by the reconciliation step (see the online appendix). In conclusion, our results suggest that using multi-step residuals or bootstrap techniques may help find a “good" estimate of the covariance matrix, which can be further improved by the reconciliation.

Simulation experiment. Frobenius norm between the true and the estimated covariance matrix for different reconciliation approaches and different techniques for simulating the base forecasts. Entries in bold represent the lowest value for each column, while the blue entry represent the global minimum. The reconciliation approaches are described in <ref>.

The accuracy of the probabilistic forecasts is evaluated using the Continuous Ranked Probability Score (CRPS, ), given by

\begin{equation}\label{eq:crps}
\operatorname{CRPS}(\widehat{P}_i, z_i)=\frac{1}{L} \sum_{l=1}^{L}\left|x_{i,l}-z_i\right|-\frac{1}{2 L^{2}} \sum_{l=1}^{L} \sum_{j=1}^{L}\left|x_{i,l}-x_{i,j}\right|, \quad i = 1,\dots,n,
\end{equation}

where $\widehat{P}_i(\omega)=\displaystyle\frac{1}{L} \sum_{l=1}^{L} \mathbf{1}\left(x_{i,l} \leq \omega\right)$, $\xvet_{1}, \xvet_{2}, \dots, \xvet_{L}\in \mathbb{R}^{n}$ is a collection of $L$ random draws from the predictive distribution and $\zvet \in \mathbb{R}^{n}$ is the observation vector. CRPS is an index that considers the single series and provides us a marginal evaluation of the approaches. In addition, we employ the Energy Score (ES, ), that is the CRPS extension to the multivariate case, to evaluate the forecasting accuracy for the whole system <cit.>,

\begin{equation}\label{eq:es}
\operatorname{ES}(\widehat{P}, \zvet)=\frac{1}{L} \sum_{l=1}^{L}\left\|\xvet_{l}-\zvet\right\|_{2}-\frac{1}{2(L-1)} \sum_{i=1}^{L-1}\left\|\xvet_{l}-\xvet_{l+1}\right\|_{2}
\end{equation}

where $	\lVert \cdot \rVert_2$ is the L$_2$ norm. In particular, we consider the geometric mean of the relative CRPS <cit.>, and the relative ES:

\begin{equation}\label{eq:skill}
\operatorname{AvgRelCRPS}_{j,s}^{[k]} = \left(\prod_{i = 1}^n \frac{CRPS^{[k]}_{i, j, s}}{CRPS^{[k]}_{i, 0, 0}}\right)^{\frac{1}{n}} \qquad \mathrm{and} \qquad \operatorname{RelES}_{j,s}^{[k]} = \frac{ES^{[k]}_{j, s}}{ES^{[k]}_{0, 0}},
\end{equation}

where $j$ denotes the reconciliation approach and $s$ indicates the approach used to simulate the base forecasts. As a reference approach ($s=0$ and $j=0$), we consider the base forecasts produce by the Bootstrap approach. If we consider all the temporal aggregation orders (i.e. $\forall k \in \{2,1\}$), the overall accuracy index is given by

\begin{equation}\label{eq:skillCRPS_all}
\operatorname{AvgRelCRPS}_{j,s} = \left(\prod_{\substack{i = 1, \dots, n \\ k \in \{1,2\}}}\frac{CRPS^{[k]}_{i, j, s}}{CRPS^{[k]}_{i, 0, 0}}\right)^{\frac{1}{n(k^\ast+m)}},
\end{equation}

\begin{equation}\label{eq:skillES_all}
\operatorname{AvgRelES}_{j,s}= \left(\prod_{k \in \{1,2\}}\frac{ES^{[k]}_{j, s}}{ES^{[k]}_{0, 0}}\right)^{\frac{1}{(k^\ast+m)}}.
\end{equation}

In Tables <ref> and <ref> are reported the results using these two scores where low values indicate better quality of the forecasts.

A limitation of this simulation setting is that we are using a high number of residuals, which may result in undervaluing the benefit from using the parameterization form of the covariance matrix such as $HB$, $H$, or $B$ (see <ref>). Additionally, shrinkage techniques often yield very similar results when we use the corresponding matrix with $\lambda = 0$ (full covariance matrix).

The good performance of the ct$(bu)$ approach can be explained by a good quality of the base forecasts at the bottom level for $k=1$, and therefore it is difficult for the other approaches to correctly adjust them using the somewhat less good forecasts of the higher temporal and cross-sectional levels. This also explains the good performance of ct$(shr_{cs}, bu_{te})$, which by definition only takes into account the information provided by the most temporally disaggregated base forecasts.
Looking at the optimal cross-temporal reconciliation approaches, it does not seem to be any advantage in using multi-step residuals to calculate the covariance matrix in the reconciliation step.

Simulation experiment. AvgRelCRPS defined in (<ref>) and (<ref>). Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

Simulation experiment. ES ratio indices defined in (<ref>) and (<ref>). Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

In conclusion, we found that simulating base forecasts from multi-step residuals allows us to estimate a covariance matrix close to the true one. Additionally, we observed that reconciliation could be used to further improve the accuracy of these estimates: accurate base forecasts for $k=1$ assist the good performance for bottom-up and optimal cross-temporal reconciliation approaches, such as oct$(wlsv)$ and oct$(bdshr)$, which perform well in terms of both CRPS and ES.

§ FORECASTING AUSTRALIAN GDP

The Australian Quarterly National Accounts (QNA) dataset has been widely used in the literature on forecast reconciliation <cit.>. Building on these results (see the online appendix), now we consider the cross-temporally reconciled probabilistic forecasts.

We use univariate ARIMA models[We use the  function from the R package  <cit.>.] to obtain quarterly base forecasts for the $n = 95$ QNA time series, spanning the period 1984:Q4 – 2018:Q1, defining GDP from both Income and Expenditure sides. We perform a rolling forecast experiment with an expanding window: the first training sample spans the period 1984:Q4 to 1994:Q3, and the last ends in 2017:Q1, for a total of 91 forecast origins. For the temporal aggregation dimension we aggregate the quarterly data to both semi-annual and annual. We obtain $4$-step, $2$-step and $1$-step ahead base forecasts respectively from the quarterly, semi-annual and annual frequencies, i.e., $\mathcal{K} = \{4,2,1\}$.

The base forecast samples in the Gaussian case are obtained using the sample covariance matrices with the Global (G) and High frequency (H) parameterization (<ref>), since it is not possible to identify a unique representation for the other cases[When simultaneously considering Income and Expenditure sides hierarchies, the result is a general linearly constrained time series, where bottom and upper time series are not uniquely defined, making unfeasible the cross-sectional bottom-up reconciliation approach <cit.>.]. We compare the results obtained using multi-step residuals with and without overlapping, in order to measure the benefit of obtaining overlapping residuals. In the non-parametric case, we use the cross-temporal joint bootstrap presented in <ref>. Finally, to reconcile the resulting (1000) base forecasts samples, we have applied the following techniques[In the online appendix, we show the results with shrunk covariance matrices. We also report the results obtained using one-step residuals in the reconciliation.] (see <ref>): ct$(shr_{cs}, bu_{te})$, ct$(wls_{cs}, bu_{te})$, oct$_o(wlsv)$, oct$_o(bdshr)$, oct$_{oh}(shr)$ and oct$_{oh}(hshr)$.

§.§ Results

Forecasting accuracy indices based on CRPS and ES are presented in Tables <ref> and <ref>, respectively. As a benchmark approach, we use the base forecasts calculated using the bootstrap method. For base forecasts, we find that using a parametric approach with the normal distribution performs better than the non-parametric bootstrap approach. This is likely due to the limited number of residuals available for bootstrapping, which does not allow for sufficient exploration of the data. Diagonal covariance matrices were found to be more effective than matrices that try to recover the correlation structure through shrinkage. Among all the procedures, ct$(wls_{cs},bu_{te})$ and oct$_o(wlsv)$ show the greatest relative gains. In contrast, oct$_{oh}(shr)$ and $oct_{oh}(hshr)$ do not show much improvement. Furthermore, the greatest improvements are observed for higher temporal aggregation levels.

AvgRelCRPS defined in (<ref>) and (<ref>) for the Australian QNA dataset. Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

ES ratio indices defined in (<ref>) and (<ref>) for the Australian QNA dataset. Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

MCB Nemenyi test for the Australian QNA dataset using CRPS at different temporal aggregation levels for the Gaussian (using overlapping and multi-step residuals, H) and the non-parametric bootstrap approaches. In each panel, the Friedman test p-value is reported in the lower right corner. The mean rank of each approach is shown to the right of its name. Statistical differences in performance are indicated if the intervals of two forecast reconciliation procedures do not overlap. Thus, approaches that do not overlap with the blue interval are considered significantly worse than the best, and vice-versa.

Additionally, we utilize the non-parametric Friedman test and the post hoc “Multiple Comparison with the Best" (MCB) Nemenyi test <cit.> to determine if the forecasting performances of the different techniques are significantly different from one another.
<ref> presents the MCB using the CRPS. We found that ct$(wls_{cs},bu_{te})$ and oct$_o(wlsv)$ are significantly better than the base forecasts at any level of aggregation.

Overall, we find that using overlapping residuals almost always leads to a greater improvement in terms of ES and CRPS, and this is also generally true for reconciliation. Forecasts at the most aggregated level (year) seem to benefit the most from reconciliation, and using one-step overlapping residuals appears to be sufficient to improve forecasts if the generation of the base forecasts sample paths take into account the multi-step structure.

§ FORECASTING AUSTRALIAN TOURISM DEMAND

The Australian Tourism Demand dataset <cit.> measures the number of nights Australians spent away from home. It includes 228 monthly observations of Visitor Nights (VN) from January 1998 to December 2016, and has a cross-sectional grouped structure based on a geographic hierarchy crossed by purpose of travel. The geographic hierarchy comprises seven states, 27 zones, and 76 regions, for a total of 111 nested geographic divisions. Six of these zones (see the online appendix) are each formed by a single region, resulting in a total of 105 unique nodes in the hierarchy. The purpose of travel comprises four categories: holiday, visiting friends and relatives, business, and other.
To avoid redundancies <cit.>, 24 nodes are not considered, resulting in an unbalanced hierarchy of 525 unique nodes instead of the theoretical 555 with duplicated nodes.
The dataset includes the 304 bottom series, which are aggregated into 221 upper time series. <ref> omits duplicated entries and updates the information in Table 7 from <cit.>. This data can be temporally aggregated into 2, 3, 4, 6, or 12 months ($\mathcal{K} = \{12,4,3,2,1\}$).

3cNumber of series

GD     PT     Tot.

Australia     1     4     5

States     7     28     35

Zones$^*$     21     84     105

Regions     76     304     380

Total     105                                      420       525

4l6cm* 6 Zones with only one Region are included in Regions.

Grouped time series for the Australian Tourism Demand dataset. GD: Geographic Division; PT: Purpose of Travel.

The forecasting experiment involves using a recursive training sample with an expanding window. The process begins by using the first 10 years of data, from January 1998 to December 2008, to generate forecasts for the entire following year (2009). Then, the training set is increased by one month. This process is repeated until the last training set is used (January 1998 to December 2015) with a total of 85 different forecast origins. For the temporal aggregation dimension we aggregate the monthly data up to annual data. We obtain $12$-step, $6$-step, $4$-step, $3$-step, $2$-step and $1$-step ahead base forecasts respectively from the monthly data and the aggregation over 2, 3, 4, 6, and 12 months. ETS models selected by minimizing the AICc criterion with the R package  are fitted to the log-transformed data, with the resulting base forecasts being back-transformed to produce non-negative forecasts, as described in <cit.>.

The (1000) base forecast samples are obtained using the Gaussian approach with sample[The results with shrunk covariance matrices are available in the online appendix.] covariance matrices (<ref>) using multi-step residuals[We do not include overlapping, as we are unable to correctly determine the residuals for the overlapping series using ETS models (see <ref>).] and the bootstrap approach (<ref>). For reconciliation, 11 different approaches have been adopted (see <ref>): ct$(bu)$, ct$(shr_{cs}, bu_{te})$, ct$(wlsv_{te}, bu_{cs})$, ct$(ols)$, oct$(struc)$, oct$(wlsv)$, oct$(bdshr)$, oct$_h(hbshr)$, oct$_h(bshr)$, oct$_h(hshr)$, and oct$_h(shr)$.

Negative forecasts may be produced during the reconciliation phase <cit.> thus generating unreasonable values (e.g., a negative forecast for tourism demand makes no sense). To overcome this limitation (see the online appendix), we applied the simple heuristic proposed by <cit.>. Following Theorem <ref>, we are thus able to obtain the reconciled sample respecting non-negativity constraints starting from an incoherent sample simulated from a Gaussian distribution. Finally, to evaluate the performance, we employed the Continuous Ranked Probability Score (CRPS), the Energy Score (ES), and the “Multiple Comparison with the Best" (MCB) Nemenyi test, introduced in Sections <ref> and <ref>.

§.§ Results

The CRPS and ES indices are shown, respectively, in Tables <ref> and <ref> for monthly, quarterly and annual forecasts[The complete set of results for all temporal aggregation levels is reported in the online appendix.]. These tables are divided by different temporal levels and each column uses a different approach to calculate the base forecasts, referred to as “base". The bootstrap method was used as a benchmark to calculate the accuracy indices.

AvgRelCRPS defined in (<ref>) and (<ref>) for the Australian Tourism Demand dataset.
Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

ES ratio indices defined in (<ref>) and (<ref>) for the Australian Tourism Demand dataset.
Approaches performing worse than the benchmark (bootstrap base forecasts, ctjb) are highlighted in red, the best for each column is marked in bold, and the overall lowest value is highlighted in blue. The reconciliation approaches are described in <ref>.

Base forecasts using a Gaussian approach are better in terms of both CRPS and ES compared to the bootstrap approach (the benchmark). Assumptions of truncated Gaussianity (Gaussian with negative values set to zero) may seem strict, but given the limited number of residuals, it can lead to improved forecasts in terms of CRPS and ES. Bootstrap forecasts suffer from the limited number of available residuals, leading in general to lower predictive accuracy. The Gaussian approach overcomes this limitation and provides better results. Regarding the different covariance matrix estimates for Gaussian base forecasts, there are no big differences. For this reason, using only the high frequency bottom time series can be useful to estimate fewer parameters and reduce the initial high dimensionality.

MCB Nemenyi test for the Australian Tourism Demand dataset using CRPS at different temporal aggregation orders for the Gaussian (multi-step residuals, HB) and the non-parametric bootstrap approaches. In each panel, the Friedman test p-value is reported in the lower right corner. The mean rank of each approach is shown to the right of its name. Statistical differences in performance are indicated if the intervals of two forecast reconciliation procedures do not overlap. Thus, approaches that do not overlap with the blue interval are considered significantly worse than the best, and vice-versa.

In the Gaussian case, bottom-up ct$(bu)$ and partly bottom-up techniques like ct$(shr_{cs}, bu_{te})$ and ct$(wlsv_{te}, bu_{cs})$ lead to better results than the benchmark (bootstrap base forecasts). However, it's not always guaranteed that the improvement is higher than the starting base forecasts (by comparing the value of each column). This is particularly true for higher levels of temporal aggregation (see online appendix for details). Overall, oct$(bdshr)$ in terms of CRPS is almost always the best. The shrinkage approach oct$_h(hshr)$ has good performance in the bootstrap case: it is competitive with oct$(bdshr)$ at lower temporal frequency ($k \in \{2,1\}$) and it is able to improve for $k\ge 3$. In terms of ES, oct$(bdshr)$ is still competitive, although it does not always show the best relative performance. In this case, approaches that attempt to explicitly take into account temporal and cross-sectional relationships, such as oct$_h(hbshr)$ and oct$_h(bshr)$, look better. It is also worth noting that techniques that don't make use of residuals like oct$(ols)$ and oct$(struc)$ prove to be competitive by consistently improving the base forecasts in terms of both CRPS and ES.

<ref> shows the MCB using the CRPS for the Gaussian approach using multi-step residuals (HB) and the non-parametric bootstrap approach. In general, partly bottom-up procedure improves with respect to base forecasts at monthly level, but optimal cross-temporal procedures are always better. In the bootstrap framework, we can identify a group of three procedures, oct$(bdshr)$, oct$(hshr)$ and oct$(struc)$ that are almost always in the group of the best approaches (denoted by the blue dot). In the Gaussian framework, oct$(wlsv)$, oct$(struc)$, and oct$(bdshr)$ are always significantly better than base forecasts and equivalent in terms of results for temporal aggregation orders greater than 2. For monthly series, oct$(bdshr)$ is always significantly better than all other approaches.

§ CONCLUSION

In this paper, we extend the probabilistic reconciliation setting developed by <cit.> for the cross-sectional case to the cross-temporal framework. Through appropriate notation, we show how theorems and definitions valid for the cross-sectional case can be reinterpreted and extended. The general notation proposed can help investigate extensions following different probabilistic approaches, such as those in <cit.>, <cit.> and <cit.>. We propose a Gaussian and a bootstrap approach to simulate the base forecasts able to take into account both cross-sectional and temporal relationships simultaneously, opening the way for further research into  cross-temporal probabilistic forecasting.

Moreover, we analyze the use of residuals, showing that one-step residuals fail to capture the temporal structure, and we propose multi-step residuals that can fully capture the full cross-temporal relationships. When dealing with covariance matrices (due to the high-dimensionality of the cross-temporal setting), we propose four alternative forms to reduce the number of parameters to be estimated, showing that the overlapping residuals may reduce the high-dimensionality burden by increasing the number of available residuals. We present a simple simulation to investigate using the different types of residuals in the cross-temporal setting arriving to some useful conclusions. These ideas are worth requiring further investigation in future works.

Finally, we perform empirical applications on two Australian data sets commonly used in forecast reconciliation research: GDP from Income and Expenditure sides and Australian Tourism Demand. We find that in both cases optimal cross-temporal reconciliation approaches significantly improve on base forecasts. We also compare these with partly bottom-up techniques that use uni-dimensional reconciliations (either cross-sectional or temporal) and find confirmation that simultaneously exploiting both dimensions in reconciliation gives better results, especially at higher levels of temporal aggregation. In conclusion, reconciliation approaches can play an important role to improve the accuracy of forecasts in a probabilistic framework while achieving the important attribute of producing coherent forecasts.

§ COVARIANCE APPROXIMATIONS

<ref> presents some approximations for the cross-sectional and the temporal covariance matrices. <cit.> consider the following approximations for the cross-temporal covariance matrix.

oct$(ols)$ - identity: $\Omegavet_{ct} = \Ivet_{n(k^*+m)}$.

oct$(struc)$ - structural: $\Omegavet_{ct} = \mathrm{diag}(\Svet_{ct} \mathbf{1}_{mn_b})$.

oct$(wlsv)$ - series variance scaling: $\Omegavet_{ct} = \widehat{\Omegavet}_{ct,wlsv}$, a straightforward extension of the series variance scaling matrix presented by <cit.>.

oct$(bdshr)$ - block-diagonal shrunk cross-covariance scaling: $\Omegavet_{ct} = \Pvet\widehat{\Wvet}^{BD}_{ct,shr}\Pvet'$, with $\widehat{\Wvet}^{BD}_{ct,shr}$ a block diagonal matrix where each $k-$block ($k = m,k_{p-1},\dots, 1$) is $\Ivet_{M_k} \otimes \widehat{\Wvet}^{[k]}_{shr}$, $\widehat{\Wvet}^{[k]}_{shr}$ is the shrunk estimate of the cross-sectional covariance matrix proposed by <cit.>, and $\Pvet$ is the commutation matrix such that $\Pvet \mathrm{vec}(\Yvet_{\tau}) = \mathrm{vec}(\Yvet_{\tau}')$.

oct$(shr)$ - MinT-shr:   $\Omegavet_{ct} = \hat{\lambda}\widehat{\Omegavet}_{ct,D} + (1-\hat{\lambda})\widehat{\Omegavet}_{ct}$, where $\hat{\lambda}$ is an estimated shrinkage coefficient (), $\widehat{\Omegavet}_{ct,D} = \Ivet_{n(k^\ast + m)} \odot \widehat{\Omegavet}_{ct}$ with $\odot$ denoting the Hadamard product, and $\widehat{\Omegavet}_{ct}$ is the covariance matrix of the cross-temporal one-step ahead in-sample forecast errors.

oct$(sam)$ - MinT-sam:  $\Omegavet_{ct} = \widehat{\Omegavet}_{ct}$.

Cross-sectional framework                                                         Temporal framework

identity            cs$(ols)$: $\Wvet = \Ivet_n$                                                               te$(ols)$: $\Omegavet = \Ivet_{k^\ast + m}$

structural          cs$(struc)$: $\Wvet = \mathrm{diag}(\Svet_{cs} \mathbf{1}_{nb})$                           te$(struc)$: $\Omegavet = \mathrm{diag}(\Svet_{te} \mathbf{1}_{m})$

series variance     cs$(wls)$: $\Wvet = \widehat{\Wvet}_D = \Ivet_n \odot \widehat{\Wvet}$                     te$(wlsv)$: $\Omegavet = \widehat{\Omegavet}_{wlsv}$

MinT-shr            cs$(shr)$: $\Wvet = \hat{\lambda}\widehat{\Wvet}_D + (1-\hat{\lambda})\widehat{\Wvet}$     te$(shr)$: $\Omegavet = \hat{\lambda}\widehat{\Omegavet}_D + (1-\hat{\lambda})\widehat{\Omegavet}$

MinT-sam            cs$(sam)$: $\Wvet = \widehat{\Wvet}$                                                       te$(sam)$: $\Omegavet = \widehat{\Omegavet}$

3p0.9Note: $\widehat{\Wvet}$ ($\widehat{\Omegavet}$) is the covariance matrix of the cross-sectional (temporal) one-step ahead in-sample forecast errors, $\widehat{\Omegavet}_{wlsv}$ is a diagonal matrix presented by <cit.>, and $\widehat{\Omegavet}_D = \Ivet_{k^\ast + m} \odot \widehat{\Omegavet}$, where $\odot$ denotes the Hadamard product.

Approximations for cross-sectional ($\Wvet$) and temporal ($\Omegavet$) covariance matrices.

§ ACKNOWLEDGMENTS

Tommaso Di Fonzo and Daniele Girolimetto acknowledge financial support from project PRIN2017 “HiDEA: Advanced Econometrics for High-frequency Data”, 2017RSMPZZ. Rob Hyndman acknowledges the support of the Australian Government through the Australian Research Council Industrial Transformation Training Centre in Optimisation Technologies, Integrated Methodologies, and Applications (OPTIMA), Project ID IC200100009.