lhf.Rmd

---
title: "Fast forecast reconciliation using linear models"
author:
- familyname: Ashouri
  othernames: Mahsa
  address: Institute of Service Science, National Tsing Hua University, Taiwan
  email: mahsa.ashouri@iss.nthu.edu.tw
  correspondingauthor: true
- familyname: Hyndman
  othernames: Rob J
  address: Monash University, Clayton VIC 3800, Australia
  email: rob.hyndman@monash.edu
- familyname: Shmueli
  othernames: Galit
  address: Institute of Service Science, National Tsing Hua University, Taiwan
  email: galit.shmueli@iss.nthu.edu.tw
abstract: "Forecasting hierarchical or grouped time series using a reconciliation approach involves two steps: computing base forecasts and reconciling the forecasts. Base forecasts can be computed by popular time series forecasting methods such as Exponential Smoothing (ETS) and Autoregressive Integrated Moving Average (ARIMA) models. The reconciliation step is a linear process that adjusts the base forecasts to ensure they are coherent. However using ETS or ARIMA for base forecasts can be computationally challenging when there are a large number of series to forecast, as each model must be numerically optimized for each series. We propose a linear model that avoids this computational problem and handles the forecasting and reconciliation in a single step. The proposed method is very flexible in incorporating external data. We illustrate our approach using a dataset on monthly Australian domestic tourism, as well as a simulated dataset. We compare our approach to reconciliation using ETS and ARIMA, and show that our approach is much faster while providing similar levels of forecast accuracy."
keywords: "hierarchical forecasting, grouped forecasting, reconciling forecast, linear regression"
wpnumber: 29/19
jelcodes: C10,C14,C22
blind: false
cover: true
toc: false
bibliography: references.bib
biblio-style: authoryear-comp
output:
  MonashEBSTemplates::workingpaper:
    fig_caption: yes
    fig_height: 5
    fig_width: 8
    includes:
      in_header: preamble.tex
    keep_tex: yes
    number_sections: yes
    citation_package: biblatex
---

```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  cache = FALSE,
  message = FALSE,
  warning = FALSE
)
library(tidyverse)
library(knitr)
library(kableExtra)
library(patchwork)
library(gridExtra)
```

# Introduction

Modern data collection tools have dramatically increased the amount of available time series data [@januschowski2013forecasting]. For example, the internet of things (IoT) and point-of-sale (POS) scanning produce huge volumes of time series in a short period of time. Naturally, there is an interest in forecasting these time series, yet forecasting large collections of time series is computationally challenging. Although computers have become faster and cheaper, computational intensity remains an important issue for many businesses and non-for-profit organizations that employ forecasting [@makridakis2018statistical].

In many cases, time series can be structured and disaggregated based on hierarchies or groups such as geographic location, product type, gender, etc.  Such hierarchical or grouped time series arise in a variety of applications where forecasts at different levels lead to different types of actions. For example, electricity demand and smart grid management can benefit from forecasts at the level of individual households, transformers, feeders, substations, balancing areas, cities, regions, and countries [@taieb2020hierarchical].  In tourism, forecasting at the national and various sub-national levels are fundamental to managers, government, and policy-makers for effective planning [@makoni2021hierarchical]. 

An example of a hierarchical time series is the monthly number of Australian domestic tourists, which can be disaggregated into different states and then into different zones. Figure \@ref(fig:hierarchicalexample) shows a schematic of such a hierarchical time series structure with three levels. The top level is the total series, formed by aggregating all the bottom-level series. In the middle level, series are aggregations of their own child series; for instance, series A is the aggregation of AW and AX. Finally, the bottom-level series includes the most disaggregated series. In our example, A might represent the Northern Territory (NT) state, which can be disaggregated into northern coast NT and central NT.

```{r hierarchicalexample, out.width = "200px", out.height= "170px,trim=0 0 190 0,clip=true", fig.align="center", fig.cap="An example of a two level hierarchical structure."}
knitr::include_graphics("Paper-Figures/hierarchical_example.jpg")
```

Grouped time series involve more complicated aggregation structures compared to strictly hierarchical time series. To continue our Australian tourism example, suppose we have two grouping factors for each tourist that are not nested: purpose of travel (Business/Holiday) and sex (Male/Female). The disaggregated series for each combination of purpose of travel and sex can be combined to form purpose of travel subtotals, or sex subtotals. These subtotals can be combined to give the overall total. Both subtotals are of interest.

We can think of such structures as hierarchical time series without a unique hierarchy. A schematic of this grouped time series structure is shown in Figure \@ref(fig:groupexample) with two grouping factors, each of two levels (A/B and C/D). The series in this structure can be split first into groups A and B and then subdivided further into C and D (left side), or split first into C and D and then subdivided into A and B (right side). The final disaggregation is identical in both cases, but the middle level aggregates are different.

```{r groupexample, out.width = "330px", out.height= "180px", fig.align="center", fig.cap="An example of a two level grouped structure."}
G1 <- grid::rasterGrob(as.raster(png::readPNG("Paper-Figures/Group_1.png")), interpolate = FALSE)
G2 <- grid::rasterGrob(as.raster(png::readPNG("Paper-Figures/Group_2.png")), interpolate = FALSE)
gridExtra::grid.arrange(G1, G2, ncol = 2)
```

Hierarchy and grouping structures offer useful information for improving forecasts of individual series within those structures, due to the correlation between different hierarchy/grouping levels. Useful information may be extracted from the aggregate series that would otherwise be potentially lost at the individual series [@syntetos2016supply]. Moreover, if we just forecast each series individually, we are ignoring the hierarchical or grouping structure, and the forecasts will not be 'coherent'. That is, forecasts for lower level series will not necessarily add up to forecasts for higher level series. This means forecasts will not add up in a way that is consistent with the aggregation structure of the time series collection [@fpp2].

There are several available methods that consider the hierarchical structure information when forecasting time series. These include the top-down [@gross1990disaggregation;@fliedner2001hierarchical], bottom-up [@kahn1998revisiting], middle-out and optimal combination [@hyndman2011optimal] approaches. In the top-down approach, we first forecast the total series and then disaggregate the forecast to form lower level series forecasts based on a set of historical or forecasted proportions [for details see @athanasopoulos2009hierarchical]. In the bottom-up approach, the forecasts in each level of the hierarchy can be computed by aggregating the bottom-level series forecasts. However, we may not get good upper-level forecasts because the most disaggregated series can be noisy and so their forecasts are often inaccurate. In the middle-out approach, the process can be started from one of the middle levels and other forecasts can be computed using aggregation for upper levels and disaggregation for lower levels. Finally, optimal combination uses all the forecasts for all of the series in the entire structure, and then uses an optimization process to reconcile the resulting forecasts. The advantage of the optimal combination method, compared with the other methods, is that it considers all information in the hierarchy, including any correlations among the series.

We use the same notation for both hierarchical and grouped time series, [following @fpp2] for both hierarchical and grouped time series. We denote the total series at time \(t\) by \(y_t\), and the series at node \(Z\) (subaggregation level \(Z\)) and time \(t\) by \(y_{Z,t}\). For describing the relationships between series, we use an \(N\times M\) matrix, called the `summing matrix', denoted by \(\bm{S}\), in which \(N\) is the overall number of nodes and \(M\) is the number of bottom-level nodes. For example, in Figure \@ref(fig:hierarchicalexample), \(N = 7\) and \(M = 4\), while in Figure \@ref(fig:groupexample), \(N=9\) and \(M=4\). Then we can write \(\bm{y}_t=\bm{S}\bm{b}_t\), where \(\bm{y}_t\) is a vector of all the level nodes at time \(t\) and \(\bm{b}_t\) is the vector of all the bottom-level nodes at time \(t\). For the example shown in Figure \ref{fig:groupexample}, the equation can be written as follows:
\begin{equation}\label{eq:Smatrixexample}
  \begin{pmatrix}
    y_{t}\\y_{A,t}\\y_{B,t}\\y_{C,t}\\y_{D,t}\\y_{AC,t}\\y_{AD,t}\\y_{BC,t}\\y_{BD,t}
  \end{pmatrix} =
  \begin{pmatrix}
    1&1&1&1\\1&1&0&0\\0&0&1&1\\1&0&1&0\\0&1&0&1\\1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\\
  \end{pmatrix}
  \begin{pmatrix}
    y_{AC,t}\\y_{AD,t}\\y_{BC,t}\\y_{BD,t}\\
  \end{pmatrix}.
\end{equation}

In the optimal combination method, reconciled forecasts can be computed using the following generalized least squares equation [@mint2018]
\begin{equation}\label{eq:mint}
  \tilde{\bm{y}}_{t+h}=\bm{S}(\bm{S}'\bm{W}_{h}^{-1}\bm{S})^{-1}\bm{S}'\bm{W}_{h}^{-1}\hat{\bm{y}}_{t+h},
\end{equation}
where \(\hat{\bm{y}}_{t+h}\) represents a vector of \(h\)-step-ahead base forecasts for all levels of the hierarchy, and \(\bm{W}_{h}\) is the covariance matrix of forecast errors for the \(h\)-step-ahead base forecasts.

Several possible methods for estimating \(\bm{W}_{h}\) are available. @mint2018 discuss a simple 'structural scaling' approximation whereby \(\bm{W}_{h} = k_{h} \bm{\Lambda}\) with \(k_{h}\) being a positive constant, \(\bm{\Lambda} = \text{diag}(\bm{S}\bm{1})\), and \(\bm{1}\) being a unit vector of dimension \(M\) (the number of bottom-level series). Note that \(\bm{\Lambda}\) simply contains the row sums of the summing matrix \(\bm{S}\), and that \(k_{h}\) will cancel out in \@ref(eq:mint).^[For simplicity of exposition, we used the structural scaling (wls\_struct) summing matrix for reconciliation in all the results.  In [Appendix A.1](#appendixA.1), we compare two alternatives: Tables \ref{tab:Tourismdatadifrecrolling} and \ref{tab:Tourismdatadifrecfix} display the results for estimates of \(\bm{W}_h\) using a shrinkage estimator (\texttt{mint\_shrink}) and variance scaling (\texttt{wls\_var}) for the Australian tourism example.] Thus
\begin{equation}\label{eq:mint2}
  \tilde{\bm{y}}_{t+h}=\bm{S}(\bm{S}'\bm{\Lambda}^{-1}\bm{S})^{-1}\bm{S}'\bm{\Lambda}^{-1}\hat{\bm{y}}_{t+h}.
\end{equation}

The most computationally challenging part of the optimal combination method is producing all the base forecasts that make up \(\hat{\bm{y}}_{t+h}\). In many applications, there may be thousands or even millions of individual series, and each of them must be forecast independently. The most popular time series forecasting methods, such as ETS and ARIMA models [@fpp2] , involve non-linear optimization routines to estimate the parameters via maximum likelihood estimation. Usually, multiple models are fitted for each series, and the best is selected by minimizing Akaike's Information Citerion [@akaike1998information]. This computational challenges increases with the number of lower-level series as well as in the number of aggregations of interest.

We therefore propose a new approach to compute the base forecasts that is both computationally fast while maintaining an acceptable forecasting accuracy level. Our proposed approach avoids the computational challenge of using ETS or ARIMA that require numerical optimization for each series. It is useful in terms of incorporating external data, handling missing values, and model selection. And finally, our approach handles the forecasting and reconciliation in a single step.

# Proposed approach: Linear model \label{sec:proposedapproach1}

<!-- Our proposed approach is based on using OLS-LR models for computing base forecasts. We use $\bm{X}_{Z}$ to denote the matrix of $k$ predictors corresponding to the series at node $Z$. Then we can write
\begin{equation}\label{eq:linearmodel}
  \bm{y}_Z = \bm{X}_Z \bm{\beta}_Z+\bm{\varepsilon}_Z
\end{equation}
where $\bm{y}_Z = \{y_{Z,1},\dots,y_{Z,T}\}, \bm{\beta}_Z$ is a vector of coefficients and $\bm{\varepsilon}_Z$ is an error term with mean zero and variance matrix $\sigma^2_Z\bm{I}$. Then using standard regression results \citep{SeberLee}, the OLS estimate of $\bm{\beta}_Z$ is given by
\begin{equation}\label{eq:linearcoefficientstwosteps}
  \hat{\bm{\beta}}_Z = (\bm{X}_{Z}'\bm{X}_Z)^{-1}\bm{X}_Z'\bm{y}_Z,
\end{equation}
and the base forecasts at horizon $h$ can be written as
\begin{equation}
  \hat{y}_{Z,T+h} = \bm{x}_{Z,T+h}'(\bm{X}_{Z}'\bm{X}_Z)^{-1}\bm{X}_Z'\bm{y}_Z,
\end{equation}
with corresponding variance
$$
\hat\sigma^2\left[1 + \bm{x}_{Z,T+h}'(\bm{X}_{Z}'\bm{X}_Z)^{-1}\bm{x}_{Z,T+h}\right],
$$
where $\bm{x}_{Z,T+h}$ denotes the $k$-vector of predictors for time period $T+h$ and
$$
\hat\sigma^2 = \frac{1}{T-k-1}(\bm{y}_Z - \bm{X}_Z\hat{\bm{\beta}}_Z)'(\bm{y}_Z - \bm{X}_Z\hat{\bm{\beta}}_Z).
$$

\todo[inline]{Not complete - Matrix notation: we need to write this for all $\hat{y}_Z$.} -->

Our proposed approach is based on using linear regression models for computing base forecasts. Suppose we have a linear model that we use for forecasting, and we wish to apply it to \(N\) different series which have some aggregation constraints. We have observations \(y_{t,i}\) from times \(t=1,\dots,T\) and series \(i=1,\dots,N\). Then

\begin{equation}
\label{eq:basicequation}
  y_{t,i} = \bm{\beta}_{i}' \bm{x}_{t,i} + \varepsilon_{t,i}
\end{equation}

where \(\bm{x}_{t,i}= (1, x_{t,1,i},\dots,x_{t,p,i})\) is a \((p+1)\)-vector of regression variables, \({\bm{\beta}}_i = (\beta_{0,i}, \beta_{1,i}, \beta_{2,i}, \dots, \beta_{p,i})\) is a \((p+1)\)-vector of coefficients, and \({\varepsilon}_{t,i}\) is the error\footnote{If different predictors are chosen for different series, we can still choose a constant \(p\) that includes all the predictors, and then set the relevant coefficients in the \(X\) matrix to \(0\).}. This equation for all the observations in matrix form can be written as follows: 
\begin{align}\label{eq:linearmodel}
\begin{split}
  \begin{pmatrix}
  \bm{y}_1\\
  \bm{y}_2\\
  \bm{y}_3 \\
  \vdots\\
  \bm{y}_N
  \end{pmatrix}&=
  \begin{pmatrix}
  \bm{X}_1 & 0        & 0        & \dots  & 0\\
  0        & \bm{X}_2 & 0        & \dots  & 0\\
  0        & 0        & \bm{X}_3 & \ddots & \vdots \\
  \vdots   & \vdots   & \ddots   & \ddots & 0\\
  0        & 0        & \dots    & 0      & \bm{X}_N
  \end{pmatrix}
  \begin{pmatrix}
  \bm{\beta}_1\\
  \bm{\beta}_2\\
  \bm{\beta}_3\\
  \vdots\\
  \bm{\beta}_N
  \end{pmatrix}+
  \begin{pmatrix}
  \bm{\varepsilon}_1\\
  \bm{\varepsilon}_2\\
  \bm{\varepsilon}_3\\
  \vdots \\
  \bm{\varepsilon}_N
  \end{pmatrix}
  \\[2ex]
  \bm{Y} &= \bm{XB} + \bm{E},
\end{split}
\end{align}
where \(\bm{y}_i = (y_{1,i}, y_{2,i}, \dots, y_{T,i})\) is a \(T\)-vector, \({\bm{\beta}}_i = (\beta_{0,i}, \beta_{1,i}, \beta_{2,i}, \dots, \beta_{p,i})\) is a \((p+1)\)-vector, \({\bm{\varepsilon}}_i = (\varepsilon_{1,i}, \varepsilon_{2,i}, \dots, \varepsilon_{T,i})\) is a \(T\)-vector, and \(\bm{X}_i\) is the \(T\times (p+1)\)-matrix
\begin{equation}\label{eq:Xmatrixdefinition}
  \bm{X}_i = \begin{pmatrix}
  1 & x_{1,i,1} & x_{1,i,2} & \dots & x_{1,i,p}\\
  1 & x_{2,i,1} & x_{2,i,2} & \dots & x_{2,i,p}\\
  \vdots & \vdots & \vdots & & \vdots \\
  1 & x_{T,i,1} & x_{T,i,2} & \dots & x_{T,i,p}
\end{pmatrix}.
\end{equation}

Equation \@ref(eq:linearmodel) can be written as \(\bm{Y} = \bm{X} \bm{B} + \bm{E}\), with parameter estimates given by \(\hat{\bm{B}} = (\bm{X}'\bm{X})^{-1} \bm{X}'\bm{Y}\). Then the base forecasts are obtained using
\begin{equation}\label{eq:baseforecasts}
  \hat{\bm{y}}_{t+h} = \bm{X}_{t+h}^* \hat{\bm{B}},
\end{equation}
where \(\hat{\bm{y}}_{t+h}\) is an \(N\)-vector of forecasts, \(\hat{\bm{B}}\) comprises \(N\) stacked \((p+1)\)-vectors of estimated coefficients, and \(\bm{X}_{t+h}^*\) is the \(N\times N(p+1)\) matrix
\pagebreak[3]
\begin{equation}
  \bm{X}_{t+h}^* =
  \begin{pmatrix}
  \bm{x}_{t+h,1}' & 0               & 0               & \dots  & 0\\
  0               & \bm{x}_{t+h,2}' & 0               & \dots  & 0\\
  0               & 0               & \bm{x}_{t+h,3}' & \ddots & \vdots \\
  \vdots          & \vdots          & \ddots          & \ddots & 0\\
  0               & 0               & \dots           & 0      & \bm{x}_{t+h,N}'
  \end{pmatrix}.
\end{equation}
Note that we use \(\bm{X}^*_{t}\) to distinguish this matrix, which combines \(\bm{x}_{t,i}\) across all series for one time, from \(\bm{X}_i\) which combines \(\bm{x}_{t,i}\) across all times for one series.

Finally, we can combine the two linear equations for computing base forecasts and reconciled forecasts (Equations \@ref(eq:mint2) and \@ref(eq:baseforecasts)) to obtain the reconciled forecasts with a single equation:
\begin{equation}\label{eq:singlestep}
    \tilde{\bm{y}}_{t+h} = \bm{S}(\bm{S}'\bm{\Lambda}\bm{S})^{-1}\bm{S}'\bm{\Lambda}
                            (\bm{X}_{t+h}^* \hat{\bm{B}})
                         = \bm{S}(\bm{S}'\bm{\Lambda}\bm{S})^{-1}\bm{S}'\bm{\Lambda}
                            \bm{X}_{t+h}^* (\bm{X}'\bm{X})^{-1} \bm{X}'\bm{Y}.
\end{equation}


## Simplified formulation for a fixed set of predictors ($\bf {X}$) \label{sec:proposedapproach2}

If we have the same set of predictor variables, \(\bm{X}\), for all the series, we can write Equations \@ref(eq:linearmodel) to \@ref(eq:singlestep) more easily using multivariate regression equations, and we can obtain all the reconciled forecasts for all the series in one equation. In that case, Equation \@ref(eq:linearmodel) can be rearranged as follows:
\begin{equation}\label{eq:linearmodelsameX}
  \begin{pmatrix}
  y_{11} & \dots & y_{1N}\\
  y_{21} & \dots & y_{2N}\\
  \vdots &       & \vdots\\
  y_{T1} & \dots & y_{TN}
  \end{pmatrix} =
  \begin{pmatrix}
  1      & X_{11} & \dots & X_{1p}\\
  1      & X_{21} & \dots & X_{2p}\\
  \vdots & \vdots &       & \vdots\\
  1      & X_{T1} & \dots & X_{Tp}
  \end{pmatrix}
  \begin{pmatrix}
  \beta_{01} & \dots & \beta_{0N}\\
  \beta_{11} & \dots & \beta_{1N}\\
  \vdots     &       & \vdots\\
  \beta_{p1} & \dots & \beta_{pN}
  \end{pmatrix} \\
  +
  \begin{pmatrix}
  \varepsilon_{11} & \dots & \varepsilon_{1N}\\
  \varepsilon_{21} & \dots & \varepsilon_{2N}\\
  \vdots           &       & \vdots\\
  \varepsilon_{T1} & \dots & \varepsilon_{TN}
  \end{pmatrix},
\end{equation}
where \(\bm{Y}\), \(\bm{X}\), \(\bm{B}\) and \(\bm{E}\) are now matrices of size \(T\times N\), \(T\times (p+1)\), \((p+1)\times N\) and \(T \times N\), respectively. Equations  \@ref(eq:baseforecasts) to \@ref(eq:singlestep) can be written accordingly using Equation \@ref(eq:linearmodelsameX), and here \(\bm{X}^*_{t+h,i} = \bm{X}^*_{t+h}\), where \(\bm{X}^*_{t+h}\) is an \(h\times (p+1)\) matrix. This simpler formulation translates into a computational advantage, as the \(\bm X\) matrix is smaller, thereby easing matrix multiplication operations.

<!--We can represent Equation \@ref(eq:linearmodelsameX) as follows
\begin{equation}\label{eq:linearmodelmatrixsameX}
\bm{Y}= \bm{X} \bm{B} + \bm{E},
\end{equation}
where $\bm{Y}$, $\bm{X}$, $\bm{B}$ and $\bm{E}$ are now matrices of size $T\times N$, $T\times (p+1)$, $(p+1)\times N$ and $T \times N$, respectively. Using Equation \@ref(eq:linearmodelmatrixsameX), parameter estimates are given by
\begin{equation}
\hat{\bm{B}} = (\bm{X}'\bm{X})^{-1} \bm{X}'\bm{Y}.
\end{equation}
Then the base forecasts can be written as
\begin{equation}
\begin{pmatrix}
 \hat{y}_{t+1,1} & \hat{y}_{t+1,2} & \dots & \hat{y}_{t+1,N}\\
 \hat{y}_{t+2,1} & \hat{y}_{t+2,2} & \dots & \hat{y}_{t+2,N}\\
 \vdots & \vdots & & \vdots\\
 \hat{y}_{t+h,1} & \hat{y}_{t+h,2} & \dots & \hat{y}_{t+h,N}\\
 \end{pmatrix} =
 \begin{pmatrix}
 1 & X_{t+1,1} & X_{t+1,2} & \dots & X_{t+1,p}\\
 1 & X_{t+2,1} & X_{t+2,2} & \dots & X_{t+2,p}\\
 \vdots & \vdots & & \vdots\\
 1 & X_{t+h,1} & X_{t+h,2} & \dots & X_{t+h,p}
 \end{pmatrix}
 \begin{pmatrix}
 \hat\beta_{01} & \hat\beta_{02} & \dots & \hat\beta_{0N}\\
 \hat\beta_{11} & \hat\beta_{12} & \dots & \hat\beta_{1N}\\
 \vdots & \vdots & & \vdots\\
 \hat\beta_{p1} & \hat\beta_{p2} & \dots & \hat\beta_{pN}
 \end{pmatrix},
\end{equation}
or in shorter form,
\begin{equation}\label{eq:baseforecastssameX}
\hat{\bm{Y}} = \bm{X}^* \hat{\bm{B}}.
\end{equation}
Then the reconciliationfor all the forecasted points, $h$, and all the series, $N$, can be written as
\begin{equation}\label{eq:recforecastssameX}
\tilde{\bm{Y}} = (\bm{S}'\bm{W}\bm{S})^{-1}\bm{S}'\bm{W} \hat{\bm{Y}},
\end{equation}
where again $\bm{W}$ is an $N\times N$ diagonal matrix with $(i,i)$th element $\sigma_i^2$, and $\bm{S}$ is the summation matrix containing the aggregation constraints. Also in the above equation, $\tilde{\bm{Y}}$ includes all the reconciled forecasts for all the bottom level series.

At the end, we can combine the two linear equations for computing base forecasts and reconciliation (Equations \@ref(eq:baseforecastssameX) and \@ref(eq:recforecastssameX)) to obtain the reconciled forecasts with a single equation
\begin{equation}\label{eq:singlestepsameX}
\tilde{\bm{Y}} = (\bm{S}'\bm{W}\bm{S})^{-1}\bm{S}'\bm{W}
                        (\bm{X}^* \hat{\bm{B}})
                        = (\bm{S}'\bm{W}\bm{S})^{-1}\bm{S}'\bm{W}
                        (\bm{X}^* (\bm{X}'\bm{X})^{-1} \bm{X}'\bm{Y})
\end{equation}
!-->

## OLS predictors

As an example of the \(\bm{X}_t\) matrix in Equation \@ref(eq:linearmodel), we can refer to the set of predictors proposed in @ashouri2018 for modeling trend, seasonality, and autocorrelation by using lagged values (\(y_{t-1}\), \(y_{t-2}\), \dots), trend variables, and seasonal dummy variables:
\begin{equation}\label{eq:linearmodelexample}
   y_t = \alpha_0 + \alpha_1 t + \beta_1 s_{1,t} + \cdots + \beta_{m-1} s_{m-1,t} + \gamma_1 y_{t-1} + \cdots + \gamma_p y_{t-p} + \delta z_t + \varepsilon_t.
\end{equation}
Here, \(s_{j,t}\) is a dummy variable taking value 1 if time \(t\) is in season \(j\) (\(j=1, 2, \dots, m\)), \(y_{t-k}\) is the \(k\)th lagged value for \(y_t\), and \(z_t\) is some external information at time \(t\). The seasonal period \(m\) depends on the problem; for instance, if we have daily data with day-of-week seasonality, then \(m=7\).

When a single OLS model is fitted to a collection of time series (e.g.~several bottom-level series), then trend and seasonality predictors are the same for all series and we can use the simpler multivariate regression models in Equation \@ref(eq:linearmodelsameX). However, this formulation is inappropriate when including lags, which differ for each series, and/or series-specific external series, in which case we use the formulation in Equation \@ref(eq:linearmodel).

<!--Because of using lags and external series as predictors in Equation \@ref(eq:linearmodelexample), we do not have same set of predictors for all the series, $y_t$. However, if we just use trend and seasonality dummies as the predictors, then the simpler equations, Equation \@ref(eq:linearmodelsameX),  can be written using multivariate regression models.!-->

When there are many options for choosing predictors, such as many seasonal dummy variables, lags, or high order trend terms, we can consider using  model selection criteria such as Akaike's Information Criterion (AIC) or model selection approaches such as leave-one-out cross-validation (LOOCV) to select the best set of predictors in terms of forecast accuracy. In practice, LOOCV can be computationally heavy except in the special case of linear model [@christensenplane] aand therefore using linear models provide a viable solution. Also, when the number of seasons \(m\) is large (e.g. in hourly data), Fourier terms can result in fewer predictors than dummy variables. The number of Fourier terms can also be determined using the same AIC or LOOCV approach  [@fpp2].

<!--While OLS is popular in practice for forecasting time series, it is often frowned upon due to its independence assumption. This can cause issues for parametric inference but is less of a problem for forecasting. In fact it often performs sufficiently well for forecasting as can be seen by its popular use in practice. Further, the use of autoregressive terms in the above model should model most of the autocorrelation in the data.-->

## Computational considerations \label{sec:computationalconsiderations}

There are two ways for implementing the above equations and} computing forecasts. First, we can create the matrices \(\bm{Y}\), \(\bm{X}\) and \(\bm{E}\), and then directly use the above equations to obtain the forecasts by taking advantage of sparse matrix routines [@bunch2014sparse], such as R's Matrix package (@maechler20062nd) to obtain the forecasts. Alternatively, we could use separate regression models to compute the coefficients for each linear model individually. Although the matrix \(\bm{X}'\bm{X}\) which we need to invert is sparse and block diagonal, we find that it is still faster to use the second approach involving separate regression models.

## Prediction intervals

For obtaining prediction intervals, we need to compute the variance of reconciled forecasts as follows [@mint2018]:
\begin{equation}\label{eq:variance}
    \text{Var}(\tilde{\bm{y}}_{t+h})
        = \bm{S}\bm{G}{\bm{\Sigma}_{t+h}} \bm{G}'\bm{S}',
\end{equation}
where \(\bm{G} = (\bm{S}'\bm{\Lambda}\bm{S})^{-1}\bm{S}'\bm{\Lambda}\) and \({\bm{\Sigma}_{t+h}}\) denotes the variance of the base forecasts given by the usual linear model formula [@fpp2]
\begin{equation}\label{eq:recvariance}
  \bm{\Sigma}_{t+h} = \sigma^2\left[1 + \bm{X}_{t+h}^*(\bm{X}'\bm{X})^{-1}(\bm{X}_{t+h}^*)'\right].
\end{equation}
where \(\sigma^2\) is the variance of the base model residuals. Assuming normally distributed errors, we can easily obtain any required prediction intervals corresponding to elements of \(\tilde{\bm{y}}_{t+h}\) using the diagonals of \eqref{eq:variance}.

# Applications

In this section we illustrate our approach using a real dataset and a simulated dataset^[All methods were run on a Linux server with Intel Xeon Silver 4108 (1.80GHz / 8-Cores / 11MB Cache)*2 and 8GB DDR4 2666 DIMM ECC Registered Memory. R version 3.6.1. All the displayed computation times are for the unreconciled point forecasts.]. The real data study includes forecasting monthly Australian domestic tourism. This dataset contains 304 series with both hierarchical and grouped structure with strong seasonality. In the simulation studies, we simulate series based on the monthly Australian domestic tourism data and systematically modify the forecasting horizon, noise level, hierarchy levels, and number of series. In  [Appendix B](#appendixB),  we use another real dataset involving daily Wikipedia pageviews. In all these cases, we compare the forecasting accuracy of ETS, ARIMA^[ For running ETS and ARIMA, we applied \texttt{ETS} and \texttt{ARIMA} functions from the \texttt{fable} package [@o2019fable]. The two sets of functions were run independently and not immediately one after the other.] and the proposed linear OLS forecasting model, with and without the reconciliation step. In these applications, we use the weighted reconciliation approach from Equation \@ref(eq:mint2). For comparing these methods, we use the average of Root Mean Square Errors (RMSEs) across all series and also display box plots for forecast errors along with the raw forecast errors. A forecast error is defined as the difference between an observed value and its forecast. To aid visibility, we suppress plotting the outliers.

<!--The two real datasets differ in terms of structure, size and behavior.The tourism data contains 304 series with both hierarchical and grouped structure, while the Wikipedia pageviews dataset contains 913 series with grouped structure. The tourism dataset has strong seasonality while the Wikipedia data are noiser.-->

We apply two methods for generating forecasts that align with two different practical forecasting scenarios. The first approach is \emph{rolling origin} forecasting, where we generate one-step-ahead forecasts (\(\tilde{\bm{y}}_{t+1}\) where \(t\) changes). This mimics the scenario where data are refreshed every time period. In the second \emph{fixed origin} method, forecasts are generated at a fixed time \(t\) for \(h\) steps ahead: \(\tilde{\bm{y}}_{t+1}, \tilde{\bm{y}}_{t+2},\dots, \tilde{\bm{y}}_{t+h}\) (we replace lagged values of \(y\) by their forecasts if they occur at periods after the forecast origin). See algorithms \@ref(alg:algorithmrolling) and \@ref(alg:algorithmfixed) for details of the rolling and fixed origin approaches).

\begin{algorithm}[h]
  \begin{algorithmic}[1]
  \caption{Hierarchical and grouped time series rolling origin OLS forecast reconciliation}\label{alg:algorithmrolling}

  \STATE $data \gets matrix(y_{1:T, 1:N}, x_{1:T, 1:N, 1:p})$
    \FOR{$i \in \{1,\dots,N\}$}
      \FOR{$k \in \{1,\dots,h\}$}
       \STATE $data \gets data_{(1:(T-h) +(k-1)),i}$
       \STATE $newdata \gets data_{((T-h)+k),i}$
       \STATE $fit \gets lm(y_{t,i} \sim x_{t,i,1} + x_{t,i,2} + \dots + x_{t,i,p}, data = data)$
       \STATE $\hat{y}_{t+k,i} \gets predict(fit, newdata = newdata)$
      \ENDFOR
    \ENDFOR
    \STATE ${\bm{S}} \gets summing~matrix(groups_{data})$
    \STATE $\Lambda \gets diag(\bm{S1})$
    \STATE $adjust \gets \bm{S}(\bm{S'}\Lambda\bm{S})^{-1}\bm{S}\Lambda$
    \FOR{$i \in {1, \dots,h}$}
    \STATE $\tilde{y}_{(T-h)+i,1:N} \gets adjust \times \hat{y}_{(T-h)+i, 1:N}$
    \ENDFOR
    \STATE \bf{return} $\tilde{y}_{((T-h):T,1:N}$
  \end{algorithmic}
\end{algorithm}

\begin{algorithm}[h]
  \begin{algorithmic}[1]
  \caption{Hierarchical and grouped time series fixed origin OLS forecast reconciliation}\label{alg:algorithmfixed}
  \STATE $data \gets matrix(y_{1:T, 1:N}, x_{1:T, 1:N, 1:p})$
    \FOR{$i \in \{1,\dots,N\}$}
    \STATE $data \gets data_{(1:(T-h),i}$
    \STATE $fit \gets lm(y_{t,i} \sim x_{t,i,1} + x_{t,i,2} + \dots + x_{t,i,p}, data = data)$
    \STATE $newdata \gets data_{((T-h)+1),i}$ \FOR{$k \in \{1,\dots,h\}$}
        \STATE $\hat{y}_{t+k,i} \gets predict(fit, newdata = newdata)$
        \STATE $newdata \gets \hat {data}_{((T-h)+k),i}$
      \ENDFOR
    \ENDFOR
    \STATE ${\bm{S}} \gets summing~matrix(groups_{data})$
    \STATE $\Lambda \gets diag(\bm{S1})$
    \STATE $adjust \gets \bm{S}(\bm{S'}\Lambda\bm{S})^{-1}\bm{S}\Lambda$
    \FOR{$i \in {1, \dots,h}$}
    \STATE $\tilde{y}_{(T-h)+i,1:N} \gets adjust \times \hat{y}_{(T-h)+i, 1:N}$
    \ENDFOR
    \STATE \bf{return} $\tilde{y}_{((T-h):T,1:N}$
  \end{algorithmic}
\end{algorithm}

<!--
We apply two methods for generating forecasts, which differ in how they handle unobserved lagged values as inputs. The first approach is ex post in that it uses actual values, even when they are future to the forecast origin. These values are known to us because they are in the test set. We call these *rolling origin* forecasts. In the second ex ante method, we replace lagged values of $y$ by their forecasts if they occur at periods after the forecast origin. We call these *fixed origin* forecasts.
-->

## Australian domestic tourism

```{r Readtourism, results='hide'}

OLS.rolling.tourism <- readr::read_csv("Australia_data/fc.rolling.tourism.OLS.csv") %>%
  mutate(ForecastInterval = '1Step')
OLS.fix.tourism <- readr::read_csv("Australia_data/fc.fix.tourism.OLS.csv") %>%
  mutate(ForecastInterval = '24Step')
OLS.tourism <- bind_rows(OLS.rolling.tourism, OLS.fix.tourism)

OLS.tourism <- bind_rows(OLS.tourism %>%
                                   filter(Series == 'Total') %>%
                                   mutate (Level = 'Total'),
                         OLS.tourism %>% filter(grepl('State/', Series)) %>%
                                   mutate (Level = 'State'), 
                         OLS.tourism %>% filter(grepl('Zone/', Series)) %>%
                                   mutate (Level = 'Zone'), 
                         OLS.tourism %>% filter(grepl('Region/', Series)) %>%
                                   mutate (Level = 'Region'), 
                         OLS.tourism %>% filter(Series %in% c('Purpose/Bus','Purpose/Hol', 
                                                                              'Purpose/Oth','Purpose/Vis')) %>%
                                   mutate (Level = 'Purpose'), 
                         OLS.tourism %>% filter(grepl('State x Purpose/', Series)) %>%
                                   mutate (Level = 'State x Purpose'), 
                         OLS.tourism %>% filter(grepl('Zone x Purpose/', Series)) %>%
                                   mutate (Level = 'Zone x Purpose'), 
                         OLS.tourism %>% filter( !grepl('State', Series) & !grepl('Zone', Series) & 
                                                                   !grepl('Region', Series) & !grepl('Purpose', Series) & 
                                                                   !grepl('Total', Series)) %>%
                                   mutate (Level = 'Region x Purpose'))

error.OLS.tourism <-bind_rows( dplyr::select ( OLS.tourism, error = error.rec, Level, ForecastInterval) %>% 
                                         mutate(Rec = 'rec') %>%
                                         mutate(Method = 'OLS') , dplyr::select ( OLS.tourism, error = error.unrec, Level, ForecastInterval) %>% 
                                         mutate(Rec = 'unrec') %>%
                                         mutate(Method = 'OLS'))


ets.arima.rolling.tourism <- readr::read_csv("Australia_data/fc.rolling.tourism.ets.arima.csv") %>%
  select(-h, -new_index) %>%
  mutate(ForecastInterval = '1Step')
ets.arima.fix.tourism <- readr::read_csv("Australia_data/fc.fix.tourism.ets.arima.csv") %>%
  mutate(ForecastInterval = '24Step')

ets.arima.tourism <- bind_rows(ets.arima.rolling.tourism, ets.arima.fix.tourism)

ets.arima.tourism <- bind_rows( ets.arima.tourism %>%
                                          filter(State == '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                          mutate (Level = 'Total'), 
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                          mutate (Level = 'State'), 
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone != '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                          mutate (Level = 'Zone') ,
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone != '<aggregated>', Region != '<aggregated>', Purpose == '<aggregated>') %>%
                                          mutate (Level = 'Region') ,
                                ets.arima.tourism %>%
                                          filter(State == '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose != '<aggregated>') %>%
                                          mutate (Level = 'Purpose') ,
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose != '<aggregated>') %>%
                                          mutate (Level = 'State x Purpose') ,
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone != '<aggregated>', Region == '<aggregated>',  Purpose != '<aggregated>') %>%
                                          mutate (Level = 'Zone x Purpose') ,
                                ets.arima.tourism %>%
                                          filter(State != '<aggregated>', Zone != '<aggregated>', Region != '<aggregated>', Purpose != '<aggregated>') %>%
                                          mutate (Level = 'Region x Purpose') )

ets.arima.tourism <- bind_rows(ets.arima.tourism %>% 
                                         filter(.model %in% c('arima', 'ets')) %>%
                                         mutate(Rec = 'unrec'), 
                               ets.arima.tourism %>% 
                                         filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                                         mutate(Rec = 'rec'))

ets.arima.tourism <- bind_rows(ets.arima.tourism %>% 
                                         filter(.model %in% c('arima', 'arima_adjusted')) %>%
                                         mutate(Method = 'ARIMA'), 
                               ets.arima.tourism %>% 
                                         filter(.model %in% c('ets', 'ets_adjusted')) %>%
                                         mutate(Method = 'ETS'))



error.tourism <- bind_rows(
  dplyr::select ( ets.arima.tourism, error , Level, Rec, Method, ForecastInterval), error.OLS.tourism) %>%
  mutate( facet = factor(Level,
                         levels = c("Total", "State", "Zone", "Region", "Purpose", "State x Purpose", "Zone x Purpose", "Region x Purpose")
  )
  )

error.tourism <- bind_rows( error.tourism %>%
                              filter(Level == 'Total') %>%
                              mutate (y.max = 6000 , y.min = -2000),
                            error.tourism %>%
                              filter(Level == 'State') %>%
                              mutate (y.max = 1500 , y.min = -500),
                            error.tourism %>%
                              filter(Level == 'Zone') %>%
                              mutate (y.max = 400 , y.min = -400),
                            error.tourism %>%
                              filter(Level == 'Region') %>%
                              mutate (y.max = 200 , y.min = -100),
                            error.tourism %>%
                              filter(Level == 'Purpose') %>%
                              mutate (y.max = 2000 , y.min = -1000),
                            error.tourism %>%
                              filter(Level == 'State x Purpose') %>%
                              mutate (y.max = 400 , y.min = -200),
                            error.tourism %>%
                              filter(Level == 'Zone x Purpose') %>%
                              mutate (y.max = 100, y.min = -100),
                            error.tourism %>%
                              filter(Level == 'Region x Purpose') %>%
                              mutate (y.max = 60, y.min = -60))

#### different reconciliation approaches
fc.rolling.mint.shrink.ets.arima <- readr::read_csv("Australia_data/fc.rolling.mint.shrink.ets.arima.csv") %>%
  mutate(Rec = 'mint_shrink', ForecastInterval = '1Step') %>%
  filter(.model %in% c('ets_adjusted', 'arima_adjusted'))
fc.rolling.wls.var.ets.arima <- readr::read_csv("Australia_data/fc.rolling.wls.var.ets.arima.csv") %>%
  mutate(Rec = 'wls_var', ForecastInterval = '1Step') %>%
  filter(.model %in% c('ets_adjusted', 'arima_adjusted'))
fc.rolling.OLS.mint.shrink.wls.var <- readr::read_csv("Australia_data/fc.rolling.OLS.mint.shrink.wls.var.csv") %>%
  mutate(ForecastInterval = '1Step')
fc.fix.ets.arima.mint.shrink <- readr::read_csv("Australia_data/fc.fix.ets.arima.mint.shrink.csv") %>%
  mutate(Rec = 'mint_shrink', ForecastInterval = '24Step') %>%
  filter(.model %in% c('ets_adjusted', 'arima_adjusted'))
fc.fix.ets.arima.wls.var <-  readr::read_csv("Australia_data/fc.fix.ets.arima.wls.var.csv") %>%
  mutate(Rec = 'wls_var', ForecastInterval = '24Step') %>%
  filter(.model %in% c('ets_adjusted', 'arima_adjusted'))
fc.fix.OLS.mint.shrink.wls.var <- readr::read_csv("Australia_data/fc.fix.OLS.mint.shrink.wls.var.csv") %>%
  mutate(ForecastInterval = '24Step')

ets.arima.dif.rec <- bind_rows(fc.rolling.mint.shrink.ets.arima, fc.rolling.wls.var.ets.arima, 
                               fc.fix.ets.arima.mint.shrink, fc.fix.ets.arima.wls.var)
fc.OLS.mint.shrink.wls.var <- bind_rows(fc.rolling.OLS.mint.shrink.wls.var, fc.fix.OLS.mint.shrink.wls.var)

fc.OLS.mint.shrink.wls.var <- bind_rows(fc.OLS.mint.shrink.wls.var %>%
                           filter(Series == 'Total') %>%
                           mutate (Level = 'Total'),
                           fc.OLS.mint.shrink.wls.var %>% filter(grepl('State/', Series)) %>%
                           mutate (Level = 'State'), 
                           fc.OLS.mint.shrink.wls.var %>% filter(grepl('Zone/', Series)) %>%
                           mutate (Level = 'Zone'), 
                           fc.OLS.mint.shrink.wls.var %>% filter(grepl('Region/', Series)) %>%
                           mutate (Level = 'Region'), 
                           fc.OLS.mint.shrink.wls.var %>% filter(Series %in% c('Purpose/Bus','Purpose/Hol', 
                                                              'Purpose/Oth','Purpose/Vis')) %>%
                           mutate (Level = 'Purpose'), 
                           fc.OLS.mint.shrink.wls.var %>% filter(grepl('State x Purpose/', Series)) %>%
                           mutate (Level = 'State x Purpose'), 
                           fc.OLS.mint.shrink.wls.var %>% filter(grepl('Zone x Purpose/', Series)) %>%
                           mutate (Level = 'Zone x Purpose'), 
                           fc.OLS.mint.shrink.wls.var %>% filter( !grepl('State', Series) & !grepl('Zone', Series) & 
                                                   !grepl('Region', Series) & !grepl('Purpose', Series) & 
                                                   !grepl('Total', Series)) %>%
                           mutate (Level = 'Region x Purpose'))

error.OLS.dif.rec <- dplyr::select ( fc.OLS.mint.shrink.wls.var, error , Level, ForecastInterval, Rec) %>% 
                                 mutate(Method = 'OLS')


ets.arima.dif.rec <- bind_rows( ets.arima.dif.rec %>%
                                  filter(State == '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                  mutate (Level = 'Total'), 
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                  mutate (Level = 'State'), 
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone != '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                                  mutate (Level = 'Zone') ,
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone != '<aggregated>', Region != '<aggregated>', Purpose == '<aggregated>') %>%
                                  mutate (Level = 'Region') ,
                                ets.arima.dif.rec %>%
                                  filter(State == '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose != '<aggregated>') %>%
                                  mutate (Level = 'Purpose') ,
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone == '<aggregated>', Region == '<aggregated>', Purpose != '<aggregated>') %>%
                                  mutate (Level = 'State x Purpose') ,
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone != '<aggregated>', Region == '<aggregated>',  Purpose != '<aggregated>') %>%
                                  mutate (Level = 'Zone x Purpose') ,
                                ets.arima.dif.rec %>%
                                  filter(State != '<aggregated>', Zone != '<aggregated>', Region != '<aggregated>', Purpose != '<aggregated>') %>%
                                  mutate (Level = 'Region x Purpose') )


ets.arima.dif.rec <- bind_rows(ets.arima.dif.rec %>% 
                                 filter(.model %in% c('arima_adjusted')) %>%
                                 mutate(Method = 'ARIMA'), 
                               ets.arima.dif.rec %>% 
                                 filter(.model %in% c('ets_adjusted')) %>%
                                 mutate(Method = 'ETS'))

error.tourism.struct <- error.tourism %>%
  select(-Rec) %>%
  mutate(Rec = 'wls_struct')

error.dif.rec <- bind_rows(
  dplyr::select ( ets.arima.dif.rec, error , Level, Rec, Method, ForecastInterval), error.OLS.dif.rec, error.tourism.struct) %>%
  mutate( facet = factor(Level,
                         levels = c("Total", "State", "Zone", "Region", "Purpose", "State x Purpose", "Zone x Purpose", "Region x Purpose")
  )
  )


# Read csv files for tourism results
rmse <- error.tourism %>%
  group_by(Rec, Method, facet, ForecastInterval) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, ForecastInterval, facet, ETS, ARIMA, OLS) %>%
  mutate(
    facet = str_replace(facet, "level([0-9])", "facet \\1")
  )

rmse.dif.rec <- error.dif.rec %>%
  group_by(Rec, Method, facet, ForecastInterval) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, ForecastInterval, facet, ETS, ARIMA, OLS) %>%
  mutate(
    facet = str_replace(facet, "level([0-9])", "facet \\1")
  )

## forecast results
forecast.tourism.OLS <- bind_rows(OLS.tourism %>% filter(Series == 'Total') , OLS.tourism %>% filter(Series == 'AAAVis'))
forecast.tourism.OLS <- forecast.tourism.OLS %>% select(-X1, -error.rec, -error.unrec, -Level) 

arima.unrec <- bind_rows( ets.arima.tourism %>% filter(.model == 'arima' , State == '<aggregated>', 
                                                               Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                            select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`),
                          ets.arima.tourism %>% filter(.model == 'arima' , State == 'A', 
                                                               Zone == 'AA', Region == 'AAA', Purpose == 'Vis') %>%
                            select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`))

arima.rec <- bind_rows( ets.arima.tourism %>% filter(.model == 'arima_adjusted' , State == '<aggregated>', 
                                                             Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                          select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`),
                        ets.arima.tourism %>% filter(.model == 'arima_adjusted' , State == 'A', 
                                                             Zone == 'AA', Region == 'AAA', Purpose == 'Vis') %>%
                          select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`))

ets.unrec <- bind_rows( ets.arima.tourism %>% filter(.model == 'ets' , State == '<aggregated>', 
                                                             Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                          select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`), 
                        ets.arima.tourism %>% filter(.model == 'ets' , State == 'A', 
                                                             Zone == 'AA', Region == 'AAA', Purpose == 'Vis') %>%
                          select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`))

ets.rec <- bind_rows( ets.arima.tourism %>% filter(.model == 'ets_adjusted' , State == '<aggregated>', 
                                                           Zone == '<aggregated>', Region == '<aggregated>', Purpose == '<aggregated>') %>%
                        select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`), 
                      ets.arima.tourism %>% filter(.model == 'ets_adjusted' , State == 'A', 
                                                           Zone == 'AA', Region == 'AAA', Purpose == 'Vis') %>%
                        select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`))
forecast.tourism.data <- bind_cols (forecast.tourism.OLS, arima.unrec, arima.rec, ets.unrec, ets.rec)

forecast.tourism <- forecast.tourism.data %>%
  select(-OLS.lower.rec, -OLS.upper.rec, -OLS.lower.unrec, -OLS.upper.unrec,
         -ARIMA.lower.rec, -ARIMA.upper.rec, -ARIMA.lower.unrec, -ARIMA.upper.unrec,
         -ETS.lower.rec, -ETS.upper.rec, -ETS.lower.unrec, -ETS.upper.unrec) %>%
  gather(-Series, -ForecastInterval, -date, key = "Method", value = "Count") %>%
  mutate(
    Rec = str_extract(Method, "[a-z]*$"),
    Rec = if_else(Rec == "ctual", "Actual", Rec),
    Model = str_extract(Method, "^[A-Z]*"),
    Model = if_else(Model == "A", "Actual", Model)
  )

```

This dataset has 19 years of monthly visitor nights in Australia by Australian tourists, a measure used as an indicator of tourism activity [@mint2018]. The data were collected by computer-assisted telephone interviews with 120,000 Australians aged 15 and over [@researchAustralia2005]. The dataset includes 304 time series, each of length 228 observations. The hierarchy and grouping structure for this dataset is made using geographic and purpose of travel information.

```{r Australiageographicaldivision, results='asis'}
tibble(
  series = c(
    "Total", "1",
    "State", "2", "3", "4", "5", "6", "7", "8",
    "Zone", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
    "Region", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54"
  ),
  name = c(
    "", "Australia",
    "", "NSW", "VIC", "QLD", "SA", "WA", "TAS", "NT",
    "", "Metro NSW", "Nth Coast NSW", "Sth Coast NSW", "Sth NSW", "Nth NSW", "ACT", "Metro VIC", "West Coast VIC", "East Coast VIC", "Nth East VIC", "Nth West VIC", "Metro QLD", "Central Coast QLD", "Nth Coast QLD", "Inland QLD", "Metro SA", "Sth Coast SA", "Inland SA", "West Coast SA", "West Coast WA", "Nth WA", "Sth WA", "Sth TAS", "Nth East TAS", "Nth West TAS", "Nth Coast NT", "Central NT",
    "", "Sydney", "Central Coast", "Hunter", "North Coast NSW", "Northern Rivers Tropical NSW", "South Coast", "Snowy Mountains", "Capital Country", "The Murray", "Riverina", "Central NSW", "New England North West", "Outback NSW", "Blue Mountains", "Canberra", "Melbourne", "Peninsula", "Geelong", "Western"
  ),
  label = c(
    "", "Total",
    "", "A", "B", "C", "D", "E", "F", "G",
    "", "AA", "AB", "AC", "AD", "AE", "AF", "BA", "BB", "BC", "BD", "BE", "CA", "CB", "CC", "CD", "DA", "DB", "DC", "DD", "EA", "EB", "EC", "FA", "FB", "FC", "GA", "GB",
    "", "AAA", "AAB", "ABA", "ABB", "ABC", "ACA", "ADA", "ADB", "ADC", "ADD", "AEA", "AEB", "AEC", "AED", "AFA", "BAA", "BAB", "BAC", "BBA"
  ),
  series2 = c("Region", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "111"),
  name2 = c("", "Lakes", "Gippsland", "Phillip Island", "General Murray", "Goulburn", "High Country", "Melbourne East", "Upper Yarra", "Murray East", "Wimmera+Mallee", "Western Grampians", "Bendigo Loddon", "Macedon", "Spa Country", "Ballarat", "Central Highlands", "Gold Coast", "Brisbane", "Sunshine Coast", "Central Queensland", "Bundaberg", "Fraser Coast", "Mackay", "Whitsundays", "Northern", "Tropical North Queensland", "Darling Downs", "Outback", "Adelaide", "Barossa", "Adelaide Hills", "Limestone Coast", "Fleurieu Peninsula", "Kangaroo Island", "Murraylands", "Riverland", "Clare Valley", "Flinders Range and Outback", "Eyre Peninsula", "Yorke Peninsula", "Australia's Coral Coast", "Experience Perth", "Australia's SouthWest", "Australia's North West", "Australia's Golden Outback", "Hobart and the South", "East Coast", "Launceston, Tamar and the North", "North West", "Wilderness West", "Darwin", "Kakadu Arnhem", "Katherine Daly", "Barkly", "Lasseter", "Alice Springs", "MacDonnell"),
  label2 = c("", "BCA", "BCB", "BCC", "BDA", "BDB", "BDC", "BDD", "BDE", "BDF", "BEA", "BEB", "BEC", "BED", "BEE", "BEF", "BEG", "CAA", "CAB", "CAC", "CBA", "CBB", "CBC", "CBD", "CCA", "CCB", "CCC", "CDA", "CDB", "DAA", "DAB", "DAC", "DBA", "DBB", "DBC", "DCA", "DCB", "DCC", "DCD", "DDA", "DDB", "EAA", "EAB", "EAC", "EBA", "ECA", "FAA", "FBA", "FBB", "FCA", "FCB", "GAA", "GAB", "GAC", "GBA", "GBB", "GBC", "GBD")
) %>%
  kable(
    align = c("r", "l", "l", "r", "l", "l"),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    longtable = TRUE,
    caption = "Australian geographic hierarchical structure.",
    col.names = c("Series", "Name", "Label", "Series", "Name", "Label")
  ) %>%
  kable_styling(position = "center", font_size = 9)
```

```{r Australiahierarchystructure, out.width = "450px", out.height= "150px", fig.align="center", fig.cap="Australian geographic hierarchical structure."}
knitr::include_graphics("Paper-Figures/Australian_hierarchy_structure.jpg")
```

```{r Australiahierarchystructuremap, out.width = "450px", out.height= "360px", fig.align="center", fig.cap="Australia map, showing the eight states."}
knitr::include_graphics("Paper-Figures/ausTurRegionsBW.pdf")
```

This dataset uses the three levels of geographic divisions in Australia. In the first level, Australia is divided into seven 'States' including New South Wales (NSW), Victoria (VIC), Queensland (QLD), South Australia (SA), Western Australia (WA), Tasmania (TAS) and Northern Territory (NT). In the second and third levels, it is divided into 27 'Zones' and 76 'Regions'. This geographical structure is detailed in Table \ref{tab:Australiageographicaldivision}, shown schematically in Figure \ref{fig:Australiahierarchystructure} , and as a map chart in Figure \ref{fig:Australiahierarchystructuremap}^[Reproduced from \url{www.tra.gov.au/tra/2016/Tourism_Region_Profiles/Region_profiles/index.html}].

In addition to geography, we have four purposes of travel: Holiday (Hol), Visiting friends and relatives (Vis), Business (Bus) and Other (Oth). This results in \(76\times4 = 304\) series at the most disaggregate level. Based on the geographic hierarchy and purpose grouping, we end up with 8 aggregation levels with 555 series in total, as shown in Table \ref{tab:Australiageographicalpurposedivision}.

```{r Australiageographicalpurposedivision, results='asis'}
groups <- tribble(
  ~Aggregation_Level, ~Series,
  "Australia", 1,
  "State", 7,
  "Zone", 27,
  "Region", 76,
  "Purpose", 4,
  "State x Purpose", NA,
  "Zone x Purpose", NA,
  "Region x Purpose", NA,
  "Total", NA
)
groups[6:8, 2] <- groups[2:4, 2] * c(groups[5, 2])
groups[9, 2] <- sum(groups[1:8, 2])
groups %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Number of Australian domestic tourism series at each aggregation level."
  ) %>%
  kable_styling(
    position = "center",
    latex_options = "hold_position"
  ) %>%
  row_spec(row = 8, hline_after = TRUE)
```

We compare the forecast error distributions for all these aggregation levels (Figures \@ref(fig:boxplotrollingtourism) and \@ref(fig:boxplottourism), as well as the average RMSE (Tables \@ref(tab:Tourismdataresulrolling) and \@ref(tab:TourismdataresultRMSE)). We use the same predictors in the OLS predictor matrix for the rolling and fixed origin approaches: a quadratic trend, 11 dummy variables, and lags 1 and 12. This set is intended to capture trend and monthly seasonality. In addition, before running the model, we partition the data into training and test sets, with the last 24 months (2 years) as our test set, and the rest as our training set.

```{r Tourismdataresulrolling, results='asis', dependson="Readtourism"}
bind_cols(
  rmse %>% filter(ForecastInterval == "1Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  rmse %>% filter(ForecastInterval == "1Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Mean(RMSE) on 24 months test set for ETS, ARIMA and OLS with and without reconciliation - Rolling origin.",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 3, "Reconciled" = 3)) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```

```{r TourismdataresultRMSE, results='asis', dependson="Readtourism"}
bind_cols(
  filter(rmse, ForecastInterval == "24Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  filter(rmse, ForecastInterval == "24Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Mean(RMSE) on 24 months test set for ETS, ARIMA and OLS with and without reconciliation - Fixed origin.",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 3, "Reconciled" = 3)) %>%
  kable_styling(position = "center")
```

```{r boxplotrollingtourism, fig.align="center", fig.cap="Box plots of rolling origin forecast errors from reconciled and unreconciled ETS, ARIMA and OLS methods at each hierarchical level.", out.width="100%"}
## 1-step-ahead
boxplot.stat <- function(x) {
  coef <- 1.5
  n <- sum(!is.na(x))
  # calculate quantiles
  stats <- quantile(x, probs = c(0.0, 0.25, 0.5, 0.75, 1.0))
  names(stats) <- c("ymin", "lower", "middle", "upper", "ymax")
  iqr <- diff(stats[c(2, 4)])
  # set whiskers
  outliers <- x < (stats[2] - coef * iqr) | x > (stats[4] + coef * iqr)
  if (any(outliers)) {
    stats[c(1, 5)] <- range(c(stats[2:4], x[!outliers]), na.rm = TRUE)
  }
  return(stats)
}
error.tourism %>%
  filter(ForecastInterval == "1Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  ggplot(aes(x = id, y = error, fill = id)) +
  stat_summary(fun.data = boxplot.stat, geom = "boxplot", alpha = 0.5) +
  # geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 4, scales = "free_y") +
  geom_blank(aes(y = y.max)) +
  geom_blank(aes(y = y.min)) +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

```{r boxplottourism, fig.align="center", fig.cap="Box plots of fixed origin forecast errors for reconciled and unreconciled ETS, ARIMA and OLS methods at each hierarchical level.", out.width="100%"}
## 24-step-ahead
error.tourism %>%
  filter(ForecastInterval == "24Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  ggplot(aes(x = id, y = error, fill = id)) +
  stat_summary(fun.data = boxplot.stat, geom = "boxplot", alpha = 0.5) +
  # geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 4, scales = "free_y") +
  geom_blank(aes(y = y.max)) +
  geom_blank(aes(y = y.min)) +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

In Figures \@ref(fig:boxplotrollingtourism) and \@ref(fig:boxplottourism) we display the forecast error box plots for reconciled and unreconciled forecasts using all three methods, for the rolling origin and fixed origin scenarios, respectively. The box plots allow comparing the error distributions across the different models. We see that all methods behave similarly, and without bias, for all levels except Total and Purpose. For Total and Purpose, all methods are biased in the same direction (under-forecasting). In both cases, for OLS reconciliation has the least effect.

Note that because higher-level series have higher counts, their errors are larger in magnitude^[[Appendix A.2](#appendixA.2) shows the box plots with scaled errors, to better compare errors across all the hierarchy levels. Scaled errors are computed by subtracting the mean and dividing by the standard deviation.], Finally, we see that (as expected) for rolling origin 1-step-ahead forecasts, the error densities are closer and more tightly distributed around zero than the fixed origin multi-step-ahead forecasts.

We also compare forecast accuracy by examining time plots. Figures \@ref(fig:forecstrolling24tourismtotal) and \@ref(fig:forecstrolling24tourism) show the rolling and fixed origin forecast results for the total series and one of the bottom-level series, 'AAAVis' (Sydney - Business). Comparing the reconciled (dashed lines) and unreconciled (dotted lines) forecasts, we see that the reconciliation step improves the forecasts in this series. We also see the OLS model's forecast accuracy is similar to that of the other methods. 

```{r forecstrolling24tourismtotal, fig.align="center", fig.cap="Comparing ETS, ARIMA and OLS forecasts (reconciled and unreconciled) for `Total' series. (Top: rolling origin, Bottom: fixed origin).", out.width="100%"}
### total series
ylim <- forecast.tourism %>%
  filter(Series == "Total") %>%
  pull(Count) %>%
  range()
g1 <- forecast.tourism %>%
  filter(Series == "Total", ForecastInterval == "1Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
  scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(linetype = FALSE)
g2 <- forecast.tourism %>%
  filter(Series == "Total", ForecastInterval == "24Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
  scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(color = FALSE)
(g1 / g2)
```

```{r forecstrolling24tourism, fig.align="center", fig.cap="Comparing ETS, ARIMA and OLS forecasts (reconciled and unreconciled) for `AAAVis' bottom-level series. (Top: rolling origin, Bottom: fixed origin).", out.width="100%"}
### one of the bottom level series
ylim <- forecast.tourism %>%
  filter(Series == "AAAVis") %>%
  pull(Count) %>%
  range()
p1 <- forecast.tourism %>%
  filter(Series == "AAAVis", ForecastInterval == "1Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
    scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(linetype = FALSE)
p2 <- forecast.tourism %>%
  filter(Series == "AAAVis", ForecastInterval == "24Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
    scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(color = FALSE)
p1 / p2
```

Next, we expand the comparison to include prediction intervals. We display the prediction intervals (based on reconciled forecasts) for the total series (Figure \@ref(fig:predinttotal)) and for the bottom-level series `AAAVis' (Sydney - Visiting) (Figure \@ref(fig:predintAAAVis)).We again see that OLS is similar to the two other methods^[We compare the reconciled and unreconciled prediction intervals in [Appendix A.3](#appendixA.3).].

```{r predinttotal, fig.align="center", fig.cap="Comparing ETS, ARIMA and OLS reconciled forecasts and prediction intervals for `Total' series. (Top: rolling origin, Bottom: fixed origin).", out.width="100%"}
### total series
p1 <- forecast.tourism.data %>%
   filter(Series == "Total", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = .5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw() +
    theme(legend.position="none")

p2 <- forecast.tourism.data %>%
   filter(Series == "Total", ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
    geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()
p1 / p2
```

```{r predintAAAVis, fig.align="center", fig.cap="Comparing ETS, ARIMA and OLS reconciled forecasts and prediction intervals for `AAAVis' bottom-level series. (Top: rolling origin, Bottom: fixed origin).", out.width="100%"}
### one series example series
p1 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
    geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw() +
    theme(legend.position="none")

p2 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()

p1 / p2
```
Table \@ref(tab:Tourismdatacomputationtime) compares the computation time of the three methods for rolling and fixed origin forecasting. We see that the OLS model is much faster compared to the other methods. Also, since reconciliation is a linear process, in all methods it is very fast and does not affect computation time significantly.

Since we are using a linear model, we can easily include exogenous variables which can often be helpful in improving forecast accuracy. In this application, we tried including an 'Easter' dummy variable indicating the timing of Easter. However, we found its effect on forecast accuracy was minimal, and therefore omit it in the model reported here.

Finally, Table \@ref(tab:Tourismdatacomputationtimeappendix) shows that, as mentioned in Section \@ref(sec:computationalconsiderations), computation is faster using separate regression models compared to the matrix approach (even using sparse matrix algebra).

```{r Tourismdatacomputationtime, results='asis'}
tibble(
  model = c("ETS", "ARIMA", "OLS"),
  #rolling_unreconciled = c(, , 48.40),
  rolling_reconciled = c(14648, 30346, 48),
  #fixed_unreconciled = c(, , 17.42),
  fixed_reconciled = c(618, 1085, 18)
) %>%
  kable(
    align = c("l", rep("r", 2)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Computation time (seconds) for OLS using the matrix approach and separate regression models, with reconciliation, for rolling and fixed origin, on a 24 months test set.",
    col.names = c("", "Rolling origin", "Fixed origin"),
    digits = 1,
  ) %>%
  #column_spec(1:2, width = "3cm") %>%
  #add_header_above(c("", "Rolling origin" = 3, "Fixed origin" = 3)) %>%
  kable_styling(position = "center", full_width = FALSE)
```

```{r Readloopsvsmatrix, results='hide', eval=FALSE}
forecast.tourism.loops.matrix <- readr::read_csv("Paper-Figures/results_Tourism/loopsvsmatrix.csv")
```

```{r loopsvsmatrix, fig.align="center", fig.cap="Comparison of the forecasts obtained using a matrix approach and separate regression models to reconcile forecasts for rolling and fixed origin tourism demand (bottom level series only).", out.width="100%", eval=FALSE}
### qqplot - rec
lvsm.1 <-
  forecast.tourism.loops.matrix %>%
  filter(Level == "level7", ForecastInterval == "1Step") %>%
  ggplot(aes(x = OLS.unrec, y = OLS.unrec.matrix)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ylim(0, 2600) +
  xlab("Separate regression models") +
  ylab("Matrix approach") +
  ggtitle("Rolling origin") +
  theme_bw() +
  guides(linetype = FALSE)

lvsm.24 <-
  forecast.tourism.loops.matrix %>%
  filter(Level == "level7", ForecastInterval == "24Step") %>%
  ggplot(aes(x = OLS.unrec, y = OLS.unrec.matrix)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ylim(0, 2600) +
  xlab("Separate regression models") +
  ylab("Matrix approach") +
  ggtitle("Fixed origin") +
  theme_bw() +
  guides(linetype = FALSE)

lvsm.1 | lvsm.24
```

```{r Tourismdatacomputationtimeappendix, results='asis'}
tibble(
  method = c("Matrix approach", "Separate models"),
  #rolling_unreconciled = c(202.06, 48.40),
  rolling_reconciled = c(210, 48),
  #fixed_unreconciled = c(87.73, 16.66),
  fixed_reconciled = c(106, 18)
) %>%
  kable(
    align = c("l", rep("r", 4)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Computation time (seconds) for OLS using the matrix approach and separate regression models, with reconciliation, for rolling and fixed origin, on a 24 months test set.",
    col.names = c("", "Rolling origin", "Fixed origin"),
    digits = 1
  ) %>%
  # column_spec(1:3, width = "3cm") %>%
  # add_header_above(c("", "Rolling origin" = 1, "Fixed origin" = 1)) %>%
  kable_styling(position = "center", full_width = FALSE)
```

<!--Now since we are using a linear model for forecasting, we can easily include information about the timing of Easter to check its effect on forecasting results. We also add this information on ARIMA models and compare with the OLS forecasting model. In Tables \@ref(tab:easterroolingRMSE) and \@ref(tab:easterRMSE), we display the average RMSE of ARIMA and OLS including the easter information, ARIMAX and OLSX, across different levels with and without reconciliation. These tables are for 1-step-ahead and 24-step-ahead forecasts. Figure \@ref(fig:forecstrolling24tourism) shows the 1-step-ahead and 24-step-ahead forecast results for one of the bottom level series, AAAVis (Sydney/Visiting). In these plots we have both reconciled (solid lines) and unreconciled (dashed lines) forecasts and we see that the reconciliation step improves the forecasts in this series. We also see that the OLS model forecast accuracy is similar to the other two methods. These results show that adding this external data does not change the forecasting results significantly. However in different cases, trying different external data can be helpful in improving the forecasting results.

```{r easterroolingRMSE, results='asis', eval = FALSE}
bind_cols(
  filter(easter.rmse, ForecastInterval == "1Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  filter(easter.rmse, ForecastInterval == "1Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -Level)
) %>%
  select(Level, ARIMA, OLS, ARIMAX, OLSX, ARIMA1, OLS1, ARIMAX1, OLSX1) %>%
  kable(
    booktabs = TRUE, format = "latex", digits = 1, linesep = "",
    caption = "Mean(RMSE) for ETS, ARIMA and OLS with and without reconciliation - 1-step-ahead - Tourism dataset",
    col.names = c("Level", rep(c("ARIMA", "OLS", "ARIMAX", "OLSX"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 4, "Reconciled" = 4)) %>%
  kable_styling(position = "center")
```

```{r easterRMSE, results='asis', eval = FALSE}
bind_cols(
  filter(easter.rmse, ForecastInterval == "24Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  filter(easter.rmse, ForecastInterval == "24Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -Level)
) %>%
  select(Level, ARIMA, OLS, ARIMAX, OLSX, ARIMA1, OLS1, ARIMAX1, OLSX1) %>%
  kable(
    booktabs = TRUE, format = "latex", digits = 1, linesep = "",
    caption = "Mean(RMSE) for ETS, ARIMA and OLS with and without reconciliation - 24-step-ahead - Tourism dataset",
    col.names = c("Level", rep(c("ARIMA", "OLS", "ARIMAX", "OLSX"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 4, "Reconciled" = 4)) %>%
  kable_styling(position = "center")
```

\todo[inline]{I can't see much point including these Easter results.}
-->
In summary, for this dataset the OLS method performs similar to ETS and ARIMA in terms of forecast accuracy while being computationally much cheaper. The reason is that the OLS model only involves linear operations whereas ETS and ARIMA require nonlinear optimization methods.

## Australian domestic tourism simulation study

We provide results from two simulation studies based on the Australian domestic tourism dataset, to evaluate the sensitivity of our results to several factors: forecast horizon, noise level, and hierarchy/grouping structure.

\paragraph{Effect of Forecast Horizon and Noise Level:}

In the first study, we simulate bottom-level series similar to the real bottom-level series of the tourism data, with the same number of series and the same length. We then generate forecasts for four forecast horizons (12, 24, 36 and 48 months) with four different noise levels (standard deviation = 0.01, 0.1, 0.5 and 1)^[Since the levels of the series are different, we first scale the simulated series (subtracting the mean and dividing by standard deviation), add the white noise series, and then rescale the series.].  For the OLS approach we use the same set of predictors as in Section \ref{australian-domestic-tourism}.

```{r Readsimtourism, results='hide'}
fc.fix.noise0.01.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.fix.noise0.01.FH.csv')
fc.fix.noise0.1.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.fix.noise0.1.FH.csv')
fc.fix.noise0.5.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.fix.noise0.5.FH.csv')
fc.fix.noise1.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.fix.noise1.FH.csv')
fc.rolling.noise0.01.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.rolling.noise0.01.FH.csv') 
fc.rolling.noise0.1.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.rolling.noise0.1.FH.csv')
fc.rolling.noise0.5.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.rolling.noise0.5.FH.csv')
fc.rolling.noise1.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.rolling.noise1.FH.csv')

fc.fix.noise.FH <- bind_rows(fc.fix.noise0.01.FH, fc.fix.noise0.1.FH, fc.fix.noise0.5.FH, fc.fix.noise1.FH) %>%
  mutate(ForecastInterval = '24Step')
fc.rolling.noise.FH <- bind_rows(fc.rolling.noise0.01.FH, fc.rolling.noise0.1.FH, fc.rolling.noise0.5.FH, fc.rolling.noise1.FH) %>%
  mutate(ForecastInterval = '1Step')%>%
  select(-order)
fc.noise.FH <- bind_rows(fc.fix.noise.FH, fc.rolling.noise.FH)
fc.noise.FH <- bind_rows(fc.noise.FH %>% 
            filter(.model %in% c('arima', 'ets')) %>%
            mutate(Rec = 'unrec'), 
            fc.noise.FH %>% 
            filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
            mutate(Rec = 'rec'))
fc.noise.FH <- bind_rows(fc.noise.FH %>% 
                                     filter(.model %in% c('arima', 'arima_adjusted')) %>%
                                     mutate(Method = 'ARIMA'), 
                             fc.noise.FH %>% 
                                     filter(.model %in% c('ets', 'ets_adjusted')) %>%
                                     mutate(Method = 'ETS'))
fc.OLS.fix.noise.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.OLS.fix.noise.FH.csv')
fc.rolling.OLS.noise.FH <- readr::read_csv('Australia_data/sim/noiseFH/fc.rolling.OLS.noise.FH.csv')
error.OLS.fix.noise.FH <-bind_rows( dplyr::select ( fc.OLS.fix.noise.FH, error = error.rec, noise, FH) %>% 
                                         mutate(Rec = 'rec') %>%
                                         mutate(Method = 'OLS') %>%
                                         mutate(ForecastInterval = '24Step'), 
                                    dplyr::select ( fc.OLS.fix.noise.FH, error = error.unrec, noise, FH) %>% 
                                      mutate(Rec = 'unrec') %>%
                                      mutate(Method = 'OLS') %>%
                                      mutate(ForecastInterval = '24Step'))
error.OLS.rolling.noise.FH <-bind_rows( dplyr::select ( fc.rolling.OLS.noise.FH, error = error.rec, noise, FH) %>% 
                                      mutate(Rec = 'rec') %>%
                                      mutate(Method = 'OLS') %>%
                                      mutate(ForecastInterval = '1Step'), 
                                    dplyr::select ( fc.rolling.OLS.noise.FH, error = error.unrec, noise, FH) %>% 
                                      mutate(Rec = 'unrec') %>%
                                      mutate(Method = 'OLS') %>%
                                      mutate(ForecastInterval = '1Step'))
error.noise.FH <- bind_rows(
  dplyr::select ( fc.noise.FH, error ,noise, FH,  Rec, Method, ForecastInterval), 
  dplyr::select (error.OLS.rolling.noise.FH, error ,noise, FH,  Rec, Method, ForecastInterval), 
   dplyr::select (error.OLS.fix.noise.FH, error ,noise, FH,  Rec, Method, ForecastInterval))

forecast.OLS.noise.FH <- bind_rows(fc.OLS.fix.noise.FH %>% filter(Series == 'AAAHol', FH =='48') %>% mutate(ForecastInterval = '24Step'), 
                                   fc.rolling.OLS.noise.FH %>% filter(Series == 'AAAHol', FH =='48') %>% mutate(ForecastInterval = '1Step'))
forecast.OLS.noise.FH <- forecast.OLS.noise.FH %>% select(-OLS.var.rec, -error.rec, -error.unrec) 

forecast.OLS.noise.FH <- bind_rows(forecast.OLS.noise.FH %>% filter(noise == '0.01', ForecastInterval == '24Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '0.1', ForecastInterval == '24Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '0.5', ForecastInterval == '24Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '1', ForecastInterval == '24Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '0.01', ForecastInterval == '1Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '0.1', ForecastInterval == '1Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '0.5', ForecastInterval == '1Step'),
                                   forecast.OLS.noise.FH %>% filter(noise == '1', ForecastInterval == '1Step'))

forecast.arima.rec.noise.FH <- fc.noise.FH %>% filter(Method == 'ARIMA', Rec == 'rec', State == 'A', 
                                                               Zone == 'AA', Region == 'AAA', Purpose == 'Hol', FH =='48') %>%
                            select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`)
forecast.arima.unrec.noise.FH <- fc.noise.FH %>% filter(Method == 'ARIMA', Rec == 'unrec', State == 'A', 
                                                      Zone == 'AA', Region == 'AAA', Purpose == 'Hol', FH =='48') %>%
  select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`)
forecast.ets.rec.noise.FH <- fc.noise.FH %>% filter(Method == 'ETS', Rec == 'rec', State == 'A', 
                                                      Zone == 'AA', Region == 'AAA', Purpose == 'Hol', FH =='48') %>%
  select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`)
forecast.ets.unrec.noise.FH <- fc.noise.FH %>% filter(Method == 'ETS', Rec == 'unrec', State == 'A', 
                                                        Zone == 'AA', Region == 'AAA', Purpose == 'Hol', FH =='48') %>%
  select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`)
forecast.noise.FH <- bind_cols (forecast.OLS.noise.FH, forecast.arima.rec.noise.FH, forecast.arima.unrec.noise.FH, 
                                forecast.ets.rec.noise.FH, forecast.ets.unrec.noise.FH)

fc.fix.sim.ets.arima.levelNS8 <-  readr::read_csv('Australia_data/sim/levelNS/fc.fix.sim.ets.arima.levelNS8.csv')
fc.fix.sim.ets.arima.levelNS10 <- readr::read_csv('Australia_data/sim/levelNS/fc.fix.sim.ets.arima.levelNS10.csv')
fc.fix.sim.ets.arima.levelNS12 <- readr::read_csv('Australia_data/sim/levelNS/fc.fix.sim.ets.arima.levelNS12.csv')
fc.fix.sim.ets.arima.levelNS18 <- readr::read_csv('Australia_data/sim/levelNS/fc.fix.sim.ets.arima.levelNS18.csv')
fc.rolling.sim.ets.arima.levelNS8 <- readr::read_csv('Australia_data/sim/levelNS/fc.rolling.sim.ets.arima.levelNS8.csv')
fc.rolling.sim.ets.arima.levelNS10 <- readr::read_csv('Australia_data/sim/levelNS/fc.rolling.sim.ets.arima.levelNS10.csv')
fc.rolling.sim.ets.arima.levelNS12 <- readr::read_csv('Australia_data/sim/levelNS/fc.rolling.sim.ets.arima.levelNS12.csv')
fc.rolling.sim.ets.arima.levelNS18 <- readr::read_csv('Australia_data/sim/levelNS/fc.rolling.sim.ets.arima.levelNS18.csv')
fc.OLS.fix.levelNS <- readr::read_csv('Australia_data/sim/levelNS/OLS.fix.levelNS.csv') %>%
 mutate(cat = paste(NS,Level,  sep = '-'))
fc.OLS.rolling.levelNS <- readr::read_csv('Australia_data/sim/levelNS/OLS.rolling.levelNS.csv') %>%
 mutate(cat = paste(NS,Level,  sep = '-'))

fc.fix.sim.ets.arima.levelNS8 <- fc.fix.sim.ets.arima.levelNS8 %>%
  mutate(ForecastInterval = '24Step')
fc.rolling.sim.ets.arima.levelNS8 <- fc.rolling.sim.ets.arima.levelNS8 %>%
  mutate(ForecastInterval = '1Step')%>%
  select(-new_index)
fc.NS.level8 <- bind_rows(fc.fix.sim.ets.arima.levelNS8, fc.rolling.sim.ets.arima.levelNS8)

fc.fix.sim.ets.arima.levelNS10 <- fc.fix.sim.ets.arima.levelNS10 %>%
  mutate(ForecastInterval = '24Step')
fc.rolling.sim.ets.arima.levelNS10 <- fc.rolling.sim.ets.arima.levelNS10 %>%
  mutate(ForecastInterval = '1Step')%>%
  select(-new_index)
fc.NS.level10 <- bind_rows(fc.fix.sim.ets.arima.levelNS10, fc.rolling.sim.ets.arima.levelNS10)

fc.fix.sim.ets.arima.levelNS12 <- fc.fix.sim.ets.arima.levelNS12 %>%
  mutate(ForecastInterval = '24Step')
fc.rolling.sim.ets.arima.levelNS12 <- fc.rolling.sim.ets.arima.levelNS12 %>%
  mutate(ForecastInterval = '1Step')%>%
  select(-new_index)
fc.NS.level12 <- bind_rows(fc.fix.sim.ets.arima.levelNS12, fc.rolling.sim.ets.arima.levelNS12)

fc.fix.sim.ets.arima.levelNS18 <- fc.fix.sim.ets.arima.levelNS18 %>%
  mutate(ForecastInterval = '24Step')
fc.rolling.sim.ets.arima.levelNS18 <- fc.rolling.sim.ets.arima.levelNS18 %>%
  mutate(ForecastInterval = '1Step')%>%
  select(-new_index)
fc.NS.level18 <- bind_rows(fc.fix.sim.ets.arima.levelNS18, fc.rolling.sim.ets.arima.levelNS18)

fc.NS.level8 <- bind_rows(fc.NS.level8 %>% 
                           filter(.model %in% c('arima', 'ets')) %>%
                           mutate(Rec = 'unrec'), 
                          fc.NS.level8 %>% 
                           filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                           mutate(Rec = 'rec'))
fc.NS.level8 <- bind_rows(fc.NS.level8 %>% 
                           filter(.model %in% c('arima', 'arima_adjusted')) %>%
                           mutate(Method = 'ARIMA'), 
                          fc.NS.level8 %>% 
                           filter(.model %in% c('ets', 'ets_adjusted')) %>%
                           mutate(Method = 'ETS'))
fc.NS.level10 <- bind_rows(fc.NS.level10 %>% 
                            filter(.model %in% c('arima', 'ets')) %>%
                            mutate(Rec = 'unrec'), 
                          fc.NS.level10 %>% 
                            filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                            mutate(Rec = 'rec'))
fc.NS.level10 <- bind_rows(fc.NS.level10 %>% 
                            filter(.model %in% c('arima', 'arima_adjusted')) %>%
                            mutate(Method = 'ARIMA'), 
                          fc.NS.level10 %>% 
                            filter(.model %in% c('ets', 'ets_adjusted')) %>%
                            mutate(Method = 'ETS'))
fc.NS.level12 <- bind_rows(fc.NS.level12 %>% 
                            filter(.model %in% c('arima', 'ets')) %>%
                            mutate(Rec = 'unrec'), 
                          fc.NS.level12 %>% 
                            filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                            mutate(Rec = 'rec'))
fc.NS.level12 <- bind_rows(fc.NS.level12 %>% 
                            filter(.model %in% c('arima', 'arima_adjusted')) %>%
                            mutate(Method = 'ARIMA'), 
                          fc.NS.level12 %>% 
                            filter(.model %in% c('ets', 'ets_adjusted')) %>%
                            mutate(Method = 'ETS'))
fc.NS.level18 <- bind_rows(fc.NS.level18 %>% 
                            filter(.model %in% c('arima', 'ets')) %>%
                            mutate(Rec = 'unrec'), 
                          fc.NS.level18 %>% 
                            filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                            mutate(Rec = 'rec'))
fc.NS.level18 <- bind_rows(fc.NS.level18 %>% 
                            filter(.model %in% c('arima', 'arima_adjusted')) %>%
                            mutate(Method = 'ARIMA'), 
                          fc.NS.level18 %>% 
                            filter(.model %in% c('ets', 'ets_adjusted')) %>%
                            mutate(Method = 'ETS'))

OLS.fix.levelNS <-bind_rows( dplyr::select ( fc.OLS.fix.levelNS, error = error.rec, NS, Level) %>% 
                                      mutate(Rec = 'rec') %>%
                                      mutate(Method = 'OLS') %>%
                                      mutate(ForecastInterval = '24Step'), 
                                    dplyr::select ( fc.OLS.fix.levelNS, error = error.unrec, NS, Level) %>% 
                                      mutate(Rec = 'unrec') %>%
                                      mutate(Method = 'OLS') %>%
                                      mutate(ForecastInterval = '24Step'))
OLS.rolling.levelNS <-bind_rows( dplyr::select ( fc.OLS.rolling.levelNS, error = error.rec, NS, Level) %>% 
                                          mutate(Rec = 'rec') %>%
                                          mutate(Method = 'OLS') %>%
                                          mutate(ForecastInterval = '1Step'), 
                                        dplyr::select ( fc.OLS.rolling.levelNS, error = error.unrec, NS, Level) %>% 
                                          mutate(Rec = 'unrec') %>%
                                          mutate(Method = 'OLS') %>%
                                          mutate(ForecastInterval = '1Step'))
error.NS.Level <- bind_rows(
  dplyr::select (fc.NS.level8, error ,NS, Level,  Rec, Method, ForecastInterval), 
  dplyr::select (fc.NS.level10, error ,NS, Level,  Rec, Method, ForecastInterval), 
  dplyr::select (fc.NS.level12, error ,NS, Level,  Rec, Method, ForecastInterval), 
  dplyr::select (fc.NS.level18, error ,NS, Level,  Rec, Method, ForecastInterval), 
  dplyr::select (OLS.fix.levelNS, error ,NS, Level,  Rec, Method, ForecastInterval), 
  dplyr::select (OLS.rolling.levelNS, error ,NS, Level,  Rec, Method, ForecastInterval))

rmse.rolling.noise.FH <- error.noise.FH %>%
  filter(ForecastInterval == '1Step', Rec == "rec") %>%
  group_by(Rec, Method, noise, FH) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, noise, FH, ETS, ARIMA, OLS)

rmse.fix.noise.FH <- error.noise.FH %>%
  filter(ForecastInterval == '24Step', Rec == "rec") %>%
  group_by(Rec, Method, noise, FH) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, noise, FH, ETS, ARIMA, OLS)

com.time.noise.FH <- readr::read_csv("Australia_data/sim/noiseFH/computation.noise.FH.csv")

## changing number of levels and series
rmse.sim.level.NS.rolling <- error.NS.Level %>%
  filter(ForecastInterval == '1Step',  Rec == "rec") %>%
  group_by(Rec, Method, NS, Level) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, NS, Level, ETS, ARIMA, OLS)

rmse.sim.level.NS.fix <- error.NS.Level %>%
  filter(ForecastInterval == '24Step',  Rec == "rec") %>%
  group_by(Rec, Method, NS, Level) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, NS, Level, ETS, ARIMA, OLS)

## forecast
forecast.OLS.levelNS.fix <- fc.OLS.fix.levelNS %>% filter(Series %in% c('AAAHol', 'AAAAHol1', 'AAAAAHol1', 'AAAAAHol1G1'), 
                              cat %in% c('304-8', '608-10', '1520-12', '3040-18')) %>%
  mutate(ForecastInterval = '24Step')

forecast.OLS.levelNS.rolling <- fc.OLS.rolling.levelNS %>% filter(Series %in% c('AAAHol', 'AAAAHol1', 'AAAAAHol1', 'AAAAAHol1G1'), 
                                  cat %in% c('304-8', '608-10', '1520-12', '3040-18')) %>%
  mutate(date = forecast.OLS.levelNS.fix$date) %>%
 mutate(ForecastInterval = '1Step')


forecast.OLS.levelNS <- bind_rows(forecast.OLS.levelNS.fix, forecast.OLS.levelNS.rolling) %>%
                               select(-OLS.var.rec, -error.rec, -error.unrec) 

## level = 8
fc.NS.level8.arima.rec <- fc.NS.level8 %>% filter(Method == 'ARIMA', Rec == 'rec', State == 'A', 
                                                      Zone == 'AA', Region == 'AAA', Purpose == 'Hol', NS == '304') %>%
  select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level8.arima.unrec <- fc.NS.level8 %>% filter(Method == 'ARIMA', Rec == 'unrec', State == 'A', 
                                                  Zone == 'AA', Region == 'AAA', Purpose == 'Hol', NS == '304') %>%
  select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`, ForecastInterval)
fc.NS.level8.ets.rec <- fc.NS.level8 %>% filter(Method == 'ETS', Rec == 'rec', State == 'A', 
                                                  Zone == 'AA', Region == 'AAA', Purpose == 'Hol', NS == '304') %>%
  select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level8.ets.unrec <- fc.NS.level8 %>% filter(Method == 'ETS', Rec == 'unrec', State == 'A', 
                                                  Zone == 'AA', Region == 'AAA', Purpose == 'Hol', NS == '304') %>%
  select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`, ForecastInterval)
## level = 10
fc.NS.level10.arima.rec <- fc.NS.level10 %>% filter(Method == 'ARIMA', Rec == 'rec', level1 == 'A', State == 'AA', 
                                                  Zone == 'AAA', Region == 'AAAA', Purpose == 'Hol1', NS == '608') %>%
  select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level10.arima.unrec <- fc.NS.level10 %>% filter(Method == 'ARIMA', Rec == 'unrec', level1 == 'A', State == 'AA', 
                                                      Zone == 'AAA', Region == 'AAAA', Purpose == 'Hol1', NS == '608') %>%
  select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`, ForecastInterval)
fc.NS.level10.ets.rec <- fc.NS.level10 %>% filter(Method == 'ETS', Rec == 'rec', level1 == 'A', State == 'AA', 
                                                  Zone == 'AAA', Region == 'AAAA', Purpose == 'Hol1', NS == '608') %>%
  select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level10.ets.unrec <- fc.NS.level10 %>% filter(Method == 'ETS', Rec == 'unrec', level1 == 'A', State == 'AA', 
                                                    Zone == 'AAA', Region == 'AAAA', Purpose == 'Hol1', NS == '608') %>%
  select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`, ForecastInterval)

## level = 12
fc.NS.level12.arima.rec <- fc.NS.level12 %>% filter(Method == 'ARIMA', Rec == 'rec',level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                    Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1', NS == '1520') %>%
  select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level12.arima.unrec <- fc.NS.level12 %>% filter(Method == 'ARIMA', Rec == 'unrec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                      Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1', NS == '1520') %>%
  select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`, ForecastInterval)
fc.NS.level12.ets.rec <- fc.NS.level12 %>% filter(Method == 'ETS', Rec == 'rec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                  Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1', NS == '1520') %>%
  select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level12.ets.unrec <- fc.NS.level12 %>% filter(Method == 'ETS', Rec == 'unrec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                    Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1', NS == '1520') %>%
  select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`, ForecastInterval)
## level = 18
fc.NS.level18.arima.rec <- fc.NS.level18 %>% filter(Method == 'ARIMA', Rec == 'rec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                    Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1' , level3 == 'G1', NS == '3040') %>%
  select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level18.arima.unrec <- fc.NS.level18 %>% filter(Method == 'ARIMA', Rec == 'unrec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                      Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1' , level3 == 'G1', NS == '3040') %>%
  select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`, ForecastInterval)
fc.NS.level18.ets.rec <- fc.NS.level18 %>% filter(Method == 'ETS', Rec == 'rec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                  Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1' , level3 == 'G1', NS == '3040') %>%
  select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`, ForecastInterval)
fc.NS.level18.ets.unrec <- fc.NS.level18 %>% filter(Method == 'ETS', Rec == 'unrec', level1 == 'A', level2 == 'AA', State == 'AAA', 
                                                    Zone == 'AAAA', Region == 'AAAAA', Purpose == 'Hol1' , level3 == 'G1', NS == '3040') %>%
  select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`, ForecastInterval)

fc.NS.level.arima.rec <- bind_rows(bind_rows(fc.NS.level8.arima.rec , fc.NS.level10.arima.rec, fc.NS.level12.arima.rec, fc.NS.level18.arima.rec) %>%
                               filter(ForecastInterval == '24Step'), 
                               bind_rows(fc.NS.level8.arima.rec , fc.NS.level10.arima.rec, fc.NS.level12.arima.rec, fc.NS.level18.arima.rec) %>%
                                 filter(ForecastInterval == '1Step')) %>%
                              select(-ForecastInterval)
fc.NS.level.arima.unrec <- bind_rows(bind_rows(fc.NS.level8.arima.unrec , fc.NS.level10.arima.unrec, fc.NS.level12.arima.unrec, fc.NS.level18.arima.unrec) %>%
                                         filter(ForecastInterval == '24Step'), 
                                       bind_rows(fc.NS.level8.arima.unrec , fc.NS.level10.arima.unrec, fc.NS.level12.arima.unrec, fc.NS.level18.arima.unrec) %>%
                                         filter(ForecastInterval == '1Step'))%>%
                                      select(-ForecastInterval)
fc.NS.level.ets.rec <- bind_rows(bind_rows(fc.NS.level8.ets.rec , fc.NS.level10.ets.rec, fc.NS.level12.ets.rec, fc.NS.level18.ets.rec) %>%
                                     filter(ForecastInterval == '24Step'), 
                                   bind_rows(fc.NS.level8.ets.rec , fc.NS.level10.ets.rec, fc.NS.level12.ets.rec, fc.NS.level18.ets.rec) %>%
                                     filter(ForecastInterval == '1Step'))%>%
                                   select(-ForecastInterval)
fc.NS.level.ets.unrec <- bind_rows(bind_rows(fc.NS.level8.ets.unrec , fc.NS.level10.ets.unrec, fc.NS.level12.ets.unrec, fc.NS.level18.ets.unrec) %>%
                                       filter(ForecastInterval == '24Step'), 
                                     bind_rows(fc.NS.level8.ets.unrec , fc.NS.level10.ets.unrec, fc.NS.level12.ets.unrec, fc.NS.level18.ets.unrec) %>%
                                       filter(ForecastInterval == '1Step'))%>%
                                     select(-ForecastInterval)

forecast.NS.Level <- bind_cols (forecast.OLS.levelNS, fc.NS.level.arima.rec, fc.NS.level.arima.unrec, 
                                fc.NS.level.ets.rec, fc.NS.level.ets.unrec)


com.time.level.NS <- readr::read_csv("Australia_data/sim/levelNS/computation.level.NS.csv")
```

Table \@ref(tab:TourismdatasimrollingnoiseFH) and \ref{tab:TourismdatasimfixnoiseFH} display the average RMSE for 12 to 48 month-ahead forecasts with different noise levels. Results are shown for reconciled forecasts for both rolling and fixed origin approaches. We see that, as expected, increasing the forecast horizon and/or noise level increases the average RMSE in all the three methods. We also see that the proposed OLS approach performs similarly to ETS and ARIMA. 

```{r TourismdatasimrollingnoiseFH, results='asis', dependson="Readsimtouris"}
kable(rmse.rolling.noise.FH,
  align = c("l", rep("r", 4)),
  format = "latex",
  booktabs = TRUE,
  linesep = "",
  digits = c(0, 2, 0, 1, 1, 1),
  caption = "Mean RMSE of rolling origin forecasts for simulated data (304 bottom-level series and 8 levels of hierarchy), for different error levels,  methods (ETS, ARIMA, OLS), forecast horizons, with reconciliation.",
  col.names = c("Reconciliation", "Error", "Forecast horizon", "ETS", "ARIMA", "OLS")
) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```

```{r TourismdatasimfixnoiseFH, results='asis', dependson="Readsimtouris"}
kable(rmse.fix.noise.FH,
  align = c("l", rep("r", 6)),
  format = "latex",
  booktabs = TRUE,
  linesep = "",
  digits = c(0, 2, 0, 1, 1, 1),
  caption = "Mean RMSE of fixed origin forecasts for simulated data (304 bottom-level series and 8 levels of hierarchy), for different error levels, methods (ETS, ARIMA, OLS), forecast horizons, with reconciliation.",
  col.names = c("Reconciliation", "Error", "Forecast horizon", "ETS", "ARIMA", "OLS")
) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```

This similarity between OLS and the other methods in terms of forecast accuracy can also be seen when looking at time plots of individual series.} For example, Figures \@ref(fig:TourismdatasimPIrollingnoiseFH) and \@ref(fig:TourismdatasimPIfixnoiseFH) display one of the bottom-level series with 12 to 48 months ahead reconciled forecasts and prediction intervals, while changing the noise levels. In these two figures we see that the OLS prediction intervals are almost identical to those from ARIMA and ETS.

```{r TourismdatasimPIrollingnoiseFH, results='asis', dependson="Readsimtouris", fig.align="center", fig.cap="Comparing reconciled 'rolling origin' forecasts and prediction intervals for a sample bottom-level series across different error levels (different panels).", out.width="100%"}
forecast.noise.FH %>%
  filter(ForecastInterval == '1Step') %>%
  ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  geom_vline(xintercept = 12, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 24, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 36, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 48, linetype="dotted",
             color = "gray", size=1) +
  xlab("Horizon") +
  ylab("Count") +
  facet_wrap(. ~ noise, scales = "free_y", ncol = 2) +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red")) +
      scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  scale_x_discrete(limits = c(12,24,36,48), labels = c("12" = "h = 12", "24" = "h = 24" ,"36" = "h = 36","48" = "h = 48"))+
  theme_bw()
```



```{r TourismdatasimPIfixnoiseFH, results='asis', dependson="Readsimtouris", fig.align="center", fig.cap="Comparing reconciled 'fixed origin' forecasts and prediction intervals for a sample bottom-level series across different error levels (different panels).", out.width="100%"}
forecast.noise.FH %>%
  filter(ForecastInterval == '24Step') %>%
  ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = "ARIMA.rec")) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  geom_vline(xintercept = 12, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 24, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 36, linetype="dotted",
             color = "gray", size=1) +
  geom_vline(xintercept = 48, linetype="dotted",
             color = "gray", size=1) +
  xlab("Horizon") +
  ylab("Count") +
  facet_wrap(. ~ noise, scales = "free_y", ncol = 2) +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red")) +
      scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  scale_x_discrete(limits = c(12,24,36,48), labels = c("12" = "h = 12", "24" = "h = 24","36" = "h = 36","48" = "h = 48"))+
  theme_bw()
```

Finally, Table \@ref(tab:TourismsimcomputationtimenoiseFH) compares the computation time (seconds) for the different methods^[The table reports the computation time only for the 0.5 noise level since changing the noise level does not affect  computation time; Also, since the reconciliation step takes less than a second, computation times are practically identical for reconciled and unreconciled forecasts.]. We see that increasing the forecast horizon from 12 to 48 months increases computation time almost linearly, and that the computation time for ARIMA and ETS is significantly longer than OLS.

<!--
```{r TourismsimcomputationtimerollingnoiseFH, fig.align="center", fig.cap="Computation time (seconds) for rolling origin reconciled forecasts using ETS, ARIMA and OLS, by horizon and by different error values (panels), for 304 bottom-level series and 8 levels of hierarchy.", out.width="100%"}
com.time.noise.FH %>%
  filter(ForecastInterval == "rolling") %>%
  ggplot() +
  geom_line(aes(x = ForecastHorizon, y = ComputationTime, color = Method, linetype = Method)) +
  scale_color_manual(name = "Method", values = c(ETS = "green", ARIMA = "blue", OLS = "red")) +
  scale_linetype_manual( values = c(ETS = "solid", ARIMA = "dashed", OLS = "dotted")) +
  facet_wrap(. ~ Error, nrow = 2) +
  scale_x_continuous(breaks = c(12, 24, 36, 48)) +
  xlab("Horizon") +
  ylab("Computation time") +
  theme_light() +
  theme(
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold")
  )
```

```{r TourismsimcomputationtimefixednoiseFH, fig.align="center", fig.cap=" Computation time (seconds) for fixed origin reconciled forecasts using ETS, ARIMA and OLS, by horizon and by different error values (panels), for 304 bottom-level series and 8 levels of hierarchy.", out.width="100%"}
com.time.noise.FH %>%
  filter(ForecastInterval == "fixed") %>%
  ggplot() +
  geom_line(aes(x = ForecastHorizon, y = ComputationTime, color = Method, linetype = Method)) +
  scale_color_manual(name = "Method", values = c(ETS = "green", ARIMA = "blue", OLS = "red")) +
    scale_linetype_manual( values = c(ETS = "solid", ARIMA = "dashed", OLS = "dotted")) +
  facet_wrap(. ~ Error, nrow = 2) +
  scale_x_continuous(breaks = c(12, 24, 36, 48)) +
  xlab("Horizon") +
  ylab("Computation time") +
  theme_light() +
  theme(
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold")
  )
```
-->

```{r TourismsimcomputationtimenoiseFH, results='asis'}
tibble(
  ForecastH = c("12", "24", "36", "48"),
  rolling_ETS = c( 5188, 10429, 15724, 21093),
  rolling_ARIMA = c(3379, 7029, 10884, 14631),
  rolling_OLS = c(19, 37, 53, 69),
  fixed_ETS = c( 443, 449, 454, 460),
  fixed_ARIMA = c(293, 291, 305, 308),
  fixed_OLS = c(8, 14, 18, 24)
) %>%
  kable(
    align = c("l", rep("r", 2)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Computation time (seconds) for rolling and fixed origin reconciled forecasts using ETS, ARIMA and OLS, by horizon and for 0.5 error value, for 304 bottom-level series and 8 levels of hierarchy.",
    col.names = c( "Forecast Horizon", rep(c("ETS", "ARIMA", "OLS"), 2)),
    digits = 0,
  ) %>%
  #column_spec(1:2, width = "3cm") %>%
  add_header_above(c("", "Rolling origin" = 3, "Fixed origin" = 3)) %>%
  kable_styling(position = "center", full_width = FALSE)
```

\paragraph{Effect of Hierarchy/Grouping Structure:}
In the second simulation study, we fix the forecast horizon at \(h=24\) and the noise level at 0.5, and then create four different aggregation levels (8, 10, 12 and 18). The aggregation levels were obtained based on the original hierarchy structure of the Australian tourism example (Table \@ref(tab:Australiageographicalpurposedivision)). The new structures include the original structure and three new structures: The original structure plus one additional three-level hierarchy factor, the original structure plus two five-level additional hierarchy factors, and the original structure plus two hierarchy factors and one grouping factor. Table \@ref(tab:simlevel1) displays these four structures and the number of series in each level. 

```{r simlevel1, results='asis'}
tibble(
 Aggregation_Level = c("Australia", "Level 1", "Level 2", "State", "Zone", "Region", "Purpose", "Group 1", "Level 1 x Purpose", "Level 2 x Purpose", "State x Purpose", "Zone x Purpose", "Level 1 x Group 1", "Level 2 x Group 1", "State x Group 1", "Zone x Group 1", "Purpose x Group 1", "Bottom level", "Total"),
  Original_structure = c(1, "-", "-", 7, 27, 76, 4, "-", "-", "-", 28, 108, "-", "-", "-", "-", "-", 304, 555),
  Structure_1 = c(1, 3, "-", 7, 27, 76, 4, "-", 12, "-", 28, 108, "-", "-", "-", "-", "-", 304, 570),
   Structure_2 = c(1, 3, 5, 7, 27, 76, 4, "-", 12, 20, 28, 108, "-", "-", "-", "-", "-", 304, 595),
   Structure_3 = c(1, 3, 5, 7, 27, 76, 4, 5, 12, 20, 28, 108, 5, 6, 7, 27, 20, 304, 665)
) %>%
  kable(
    align = c("l", rep("r", 4)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Four simulated hierarchy/grouping structures.",
    col.names = c("Aggregation Level", "Original structure", "Structure 1", "Structure 2", "Structure 3")
  ) %>%
  # column_spec(1:3, width = "3cm") %>%
  add_header_above(c("", "Number of series" = 4)) %>%
  kable_styling(
    position = "center",
    full_width = FALSE
  )
```

<!---
```{r simlevel1test, results='asis'}
groups <- tribble(
  ~Division, ~Series,
  "Australia", 1,
  "Level 1", 3,
  "State", 7,
  "Zone", 27,
  "Region", 76,
  "Purpose", 4,
  "Level 1 x Purpose", NA,
  "State x Purpose", NA,
  "Zone x Purpose", NA,
  "Region x Purpose", NA,
  "Total", NA
)
groups[7:10, 2] <- groups[2:5, 2] * c(groups[6, 2])
groups[11, 2] <- sum(groups[1:10, 2])
groups %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Number of simulated Australian domestic tourism series at each aggregation level - adding one hierarchy variable (Level 1)"
  ) %>%
  kable_styling(
    position = "center",
    latex_options = "hold_position"
  ) %>%
  row_spec(row = 10, hline_after = TRUE)
```

```{r simlevel1level2, results='asis'}
groups <- tribble(
  ~Division, ~Series,
  "Australia", 1,
  "Level 1", 3,
  "Level 2", 5,
  "State", 7,
  "Zone", 27,
  "Region", 76,
  "Purpose", 4,
  "Level 1 x Purpose", NA,
  "Level 2 x Purpose", NA,
  "State x Purpose", NA,
  "Zone x Purpose", NA,
  "Region x Purpose", NA,
  "Total", NA
)
groups[8:12, 2] <- groups[2:6, 2] * c(groups[7, 2])
groups[13, 2] <- sum(groups[1:12, 2])
groups %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Number of simulated Australian domestic tourism series at each aggregation level - adding two hierarchy variables (Level 1 and Level 2)"
  ) %>%
  kable_styling(
    position = "center",
    latex_options = "hold_position"
  ) %>%
  row_spec(row = 12, hline_after = TRUE)
```

```{r simlevel1level2group1test, results='asis'}
groups <- tribble(
  ~Division, ~Series,
  "Australia", 1,
  "Level 1", 3,
  "Level 2", 5,
  "State", 7,
  "Zone", 27,
  "Region", 76,
  "Purpose", 4,
  "Group 1", 5,
  "Level 1 x Purpose", NA,
  "Level 2 x Purpose", NA,
  "State x Purpose", NA,
  "Zone x Purpose", NA,
  "Level 1 x Group 1", 5,
  "Level 2 x Group 1", 6,
  "State x Group 1", 7,
  "Zone x Group 1", 27,
  "Purpose x Group 1", NA,
  "Bottom level", 304,
  "Total", NA
)
groups[9:12, 2] <- groups[2:5, 2] * c(groups[7, 2])
groups[17, 2] <- groups[7, 2] * c(groups[8, 2])
groups[19, 2] <- sum(groups[1:18, 2])
groups %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Number of simulated Australian domestic tourism series at each aggregation level - adding two hierarchy and one grouping variables (Level 1, Level 2 and Group 1)"
  ) %>%
  kable_styling(
    position = "center",
    latex_options = "hold_position"
  ) %>%
  row_spec(row = 18, hline_after = TRUE)
```
-->
We also simulated four sets of bottom-level series with different numbers of series (304, 608, 1520 and 3040). In order to add series, we changed the number of `Purpose' categories (grouping factor) in the Australian domestic tourism example. Table \@ref(tab:TourismsimlevelNS) displays the total number of series based on 304, 608, 1520 and 3040 bottom-level series with 8, 10, 12 and 18 hierarchy levels.

```{r TourismsimlevelNS, results='asis'}
tibble(
  bottom = c(304, 608, 1520, 3040),
  levels_8 = c(555, 999, 2331, 4551),
  levels_10 = c(570, 1026, 2394, 4674),
  levels_12 = c(595, 1071, 2499, 2643),
  levels_18 = c(665, 1161, 2649, 5129)
) %>%
  kable(
    align = c("l", rep("r", 4)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Total number of series in hierarchy structure based on the number of bottom-level series and aggregation levels.",
    col.names = c("# Bottom-level series", "8", "10", "12", "18")
  ) %>%
  # column_spec(1:3, width = "3cm") %>%
  add_header_above(c("", "No. of aggregation levels" = 4)) %>%
  kable_styling(
    position = "center",
    full_width = FALSE
  )
```

Tables  \@ref(tab:TourismdatasimrollinglevelNS) and \@ref(tab:TourismdatasimfixlevelNS) display the average RMSEs by number of bottom-level series and number of aggregation levels for the different forecasting approaches. We see that increasing the number of bottom-level series and/or the number of aggregation levels increases the mean RMSE, and this similarly affects all forecasting methods, in both rolling origin and fixed origin scenarios. 

We verify this further by examining time plots and prediction intervals for one of the bottom-level series. Figures \@ref(fig:TourismdatasimPIrollinglevelNS) and \@ref(fig:TourismdatasimPIfixlevelNS) compare four aggregation levels (in different panels). We see that all three approaches indeed have similar prediction intervals across aggregation levels, for both the rolling origin and fixed origin forecasting scenarios.


```{r TourismdatasimrollinglevelNS, results='asis', dependson="Readsimtouris"}
rmse.sim.level.NS.rolling %>%
  kable(
    align = c("l", rep("r", 6)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    digits = c(0, 0, 0, 1, 1, 1),
    caption = "Mean RMSE by number of hierarchy levels, number of bottom-level series, method, with reconciliation. Simulated series has error value 0.5. Forecasting uses rolling origin for a 24-month horizon.",
    col.names = c("Reconciliation", "Series", "Levels", "ETS", "ARIMA", "OLS")
  ) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```


```{r TourismdatasimfixlevelNS, results='asis', dependson="Readsimtouris"}
rmse.sim.level.NS.fix %>%
  kable(
    align = c("l", rep("r", 6)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    digits = c(0, 0, 0, 1, 1, 1),
    caption = "Mean RMSE by number of hierarchy levels, number of bottom-level series, method, with reconciliation. Simulated series has error value 0.5. Forecasting uses fixed origin for a 24-month horizon.",
    col.names = c("Reconciliation", "Series", "Levels", "ETS", "ARIMA", "OLS")
  ) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```


```{r TourismdatasimPIrollinglevelNS, results='asis', dependson="Readsimtouris", fig.align="center", fig.cap="Comparing reconciled rolling origin forecasts and prediction intervals for a sample bottom-level series, for different number of bottom-level series and hierarchy levels (different panels). Simulated series has error value 0.5 and 24 months test set.", out.width="100%"}

forecast.NS.Level$SL_f = factor(forecast.NS.Level$cat, levels=c('304-8','608-10','1520-12', '3040-18'))

forecast.NS.Level %>%
  filter(ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
    geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  facet_wrap(. ~ SL_f, scales = "free_y", ncol = 2) +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red")) +
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()
```



```{r TourismdatasimPIfixlevelNS, results='asis', dependson="Readsimtouris", fig.align="center", fig.cap="Comparing reconciled fixed origin forecasts and prediction intervals for a sample bottom-level series, for different number of bottom-level series and hierarchy levels (different panels). Simulated series has error value 0.5 and 24 months test set.", out.width="100%"}


forecast.NS.Level$SL_f = factor(forecast.NS.Level$cat, levels=c('304-8','608-10','1520-12', '3040-18'))

forecast.NS.Level %>%
  filter(ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
    geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  facet_wrap(. ~ SL_f, scales = "free_y", ncol = 2) +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green",  "red")) +
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()
```

Finally, we compare the computation time of the different methods across the different number of aggregation levels and number of bottom-level series (Table \@ref(tab:TourismsimcomputationtimelevelNS)). We see that increasing the number of bottom-level series increases computation time linearly and increasing the number of aggregation levels only slightly increases computation time. For both rolling and fixed origin scenarios, we see a substantial difference in computation time between ARIMA and ETS compared with our OLS approach^[The reconciliation step computation time varies from less than one second to around three seconds, as a function of the number of series. Computation times displayed in the table are for unreconciled forecasts.].


<!--
 ```{r TourismsimcomputationtimerollinglevelNS, fig.align="center", fig.cap="Computation time (seconds) for rolling origin reconciled forecasts using ETS, ARIMA and OLS, by number of bottom-level series (x-axis), and by levels of hierarchy (panels). Simulated series has error value 0.5 and 24 months test set.", out.width="100%"}
 com.time.level.NS %>%
   filter(ForecastInterval == "rolling") %>%
   ggplot() +
   geom_line(aes(x = NS, y = ComputationTime, color = Method, linetype = Method)) +
   scale_color_manual(name = "Method", values = c(ETS = "green", ARIMA = "blue", OLS = "red")) +
   scale_linetype_manual( values = c(ETS = "solid", ARIMA = "dashed", OLS = "dotted")) +
   facet_wrap(. ~ Level, nrow = 2) +
   scale_x_continuous(breaks = c(304, 608, 1520, 3040)) +
   xlab("Number of series") +
   ylab("Computation time") +
   theme_light() +
   theme(
     axis.text = element_text(size = 10),
     strip.text = element_text(size = 12),
     axis.title = element_text(size = 12, face = "bold")
   )
```

 ```{r TourismsimcomputationtimefixedlevelNS, fig.align="center", fig.cap="Computation time (seconds) for fixed origin reconciled forecasts using ETS, ARIMA and OLS, by number of bottom-level series (x-axis), and by levels of hierarchy (panels). Simulated series has error value 0.5 and 24 months test set.", out.width="100%"}
 com.time.level.NS %>%
   filter(ForecastInterval == "fixed") %>%
   ggplot() +
   geom_line(aes(x = NS, y = ComputationTime, color = Method, linetype = Method)) +
   scale_color_manual(name = "Method", values = c(ETS = "green", ARIMA = "blue", OLS = "red")) +
   scale_linetype_manual( values = c(ETS = "solid", ARIMA = "dashed", OLS = "dotted")) +
   facet_wrap(. ~ Level, nrow = 2) +
   scale_x_continuous(breaks = c(304, 608, 1520, 3040)) +
   xlab("Number of series") +
   ylab("Computation time") +
   theme_light() +
   theme(
     axis.text = element_text(size = 10),
     strip.text = element_text(size = 12),
     axis.title = element_text(size = 12, face = "bold")
   )
```
-->

```{r TourismsimcomputationtimelevelNS, results='asis'}
tibble(
  Series = c("304", "", "", "", "608", "", "", "", "1520", "", "", "", "3040", "", "", ""),
  Levels = c("8", "10", "12", "18", "8", "10", "12", "18", "8", "10", "12", "18", "8", "10", "12", "18"),
  rolling_ETS = c(10490, 10468, 11618, 13737, 15049, 15538, 16322, 17891, 35323, 36202, 37584, 41458, 73406, 72505, 80636, 86112),
  rolling_ARIMA = c(7079, 7092, 7559, 8321, 10074, 10403, 10996, 12133, 24661, 25564, 26770, 28975, 54080, 58184, 61646, 64088),
  rolling_OLS = c(36, 36, 38, 44, 63, 64, 67, 73, 200, 210, 194, 226, 271, 297, 309, 326),
  fixed_ETS = c(19, 20, 20, 21, 35, 35, 36, 41, 79, 81, 85, 88, 159, 160, 165, 212),
  fixed_ARIMA = c(290, 299, 323, 356, 719, 760, 794, 836, 1935, 2037, 2186, 2414, 7691, 8081, 9247, 10437),
  fixed_OLS = c(19, 20, 20, 21, 35, 35, 36, 41, 79, 81, 85, 88, 159, 160, 165, 212)
) %>%
  kable(
    align = c("l", rep("r", 2)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Computation time (seconds) for rolling and fixed origin reconciled forecasts using ETS, ARIMA and OLS, by number of bottom-level series (x-axis), and by levels of hierarchy (panels). Simulated series has error value 0.5 and 24 months test set.",
    col.names = c("Series", "Levels", rep(c("ETS", "ARIMA", "OLS"), 2)),
    digits = 0,
  ) %>%
  #column_spec(1:2, width = "3cm") %>%
  add_header_above(c("", "", "Rolling origin" = 3, "Fixed origin" = 3)) %>%
  kable_styling(position = "center", full_width = FALSE)
```

In summary, these two simulation studies show that our result -- that OLS performs much faster but with close accuracy compared with ETS and ARIMA -- is robust to changes in forecast horizon, noise levels, and hierarchy/grouping structure. We see that computation time is not affected by the noise level while it linearly increases in the forecast horizon. Also, changing the number of aggregation levels in the hierarchy/grouping structure only slightly increases the computation time while increasing the number of bottom-level series increases computation time linearly.

\FloatBarrier

# Conclusion

We have proposed a linear model approach to fast forecasting of hierarchical or grouped time series, with accuracy that nearly matches that of forecast methods such as ETS and ARIMA. This is especially useful in large collections of time series, as is typical in hierarchical and grouped structures. Although ETS and ARIMA are advantageous in terms of forecasting power and accuracy, they can be computationally heavy when facing large collections of time series in the hierarchy.

An important feature of our proposed OLS model is its ability to easily include external information such as holiday dummies or other external series.We also note that OLS has the additional practical advantage of handling missing data while ETS and ARIMA require imputation. Another advantage of the OLS approach is its model selection flexibility, where computationally cheap criteria (e.g. AIC) and approaches (e.g. leave-one-out-cross-validation, LOOCV) can be used for choosing a best subset of predictors. We leave the study of missing data and model selection for future research.


@pennings2017 proposed another approach for forecasting hierarchical time series using state space models. Although their approach is flexible in handling outliers, missing data and external features, it is less flexible across different kinds of datasets and is computationally much more demanding.

Another advantage of our approach is that it can be computed in a single matrix equation (Equation \eqref{eq:singlestep}). This makes it extremely fast and easy to implement, and enables standard results to be derived with minimal effort (e.g., prediction intervals).

In this paper, we computed prediction intervals based on a normality assumption for all the real and simulated examples. One future direction for potentially improving prediction intervals is applying bootstrap-based methods which only assume uncorrelated forecast errors.


# Acknowledgements {-}

Ashouri and Shmueli were partially funded by Ministry of Science and Technology (MOST), Taiwan [Grant 106-2420-H-007-019]. Hyndman's research is supported by the Australian Center of Excellence in Mathematical and Statistical Frontiers.

# Supplementary material {-}
The datasets and R codes used in this paper are publicly available  at \url{github.com/robjhyndman/linear-hierarchical-forecasting}. We have also created a Binder interface for our Australian tourism and Wikipedia examples which allows easily running our code in a browser. The demo folder is reachable from \url{github.com/mahsaashouri/AUS-Wiki-Binder}.

\clearpage

# (APPENDIX) Appendix {-}

# Further analyses for tourism demand forecasting

## Different reconciliation matrices

We present reconciliation results using two other reconciliation matrices introduced in [@mint2018]: shrinkage estimator (\texttt{mint\_shrink}) and variance scaling (\texttt{wls\_var}). These results are presented for both rolling origin (Table \@ref(tab:Tourismdatadifrecrolling)) and fixed origin forecasts (Table \@ref(tab:Tourismdatadifrecfix)). We see that in this case the shrinkage estimator provides the best reconciled forecasts.

```{r Tourismdatadifrecrolling, results='asis', dependson="Readtourism"}
bind_cols(
  filter(rmse.dif.rec, ForecastInterval == "1Step", Rec == "mint_shrink") %>% select(-Rec, -ForecastInterval),
  filter(rmse.dif.rec, ForecastInterval == "1Step", Rec == "wls_var") %>% select(-Rec, -ForecastInterval, -facet),
  filter(rmse.dif.rec, ForecastInterval == "1Step", Rec == "wls_struct") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Comparing Mean(RMSE) of three   different reconciliation matrices for rolling origin forecasts on a 24 months test set.",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 3))
  ) %>%
  add_header_above(c("", "mint_shrink" = 3, "wls_var" = 3, 'wls_struct' = 3)) %>%
  kable_styling(position = "center")
```


```{r Tourismdatadifrecfix, results='asis', dependson="Readtourism"}
bind_cols(
  filter(rmse.dif.rec, ForecastInterval == "24Step", Rec == "mint_shrink") %>% select(-Rec, -ForecastInterval),
  filter(rmse.dif.rec, ForecastInterval == "24Step", Rec == "wls_var") %>% select(-Rec, -ForecastInterval, -facet),
  filter(rmse.dif.rec, ForecastInterval == "24Step", Rec == "wls_struct") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Comparing Mean(RMSE) of three   different reconciliation matrices for fixed origin forecasts on a 24 months test set.",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 3))
  ) %>%
  add_header_above(c("", "mint_shrink" = 3, "wls_var" = 3, 'wls_struct' = 3)) %>%
  kable_styling(position = "center")
```

## Scaled forecast errors boxplots

For further comparison of the forecast errors for the tourism example across different aggregation levels, we plot the \emph{scaled errors} for rolling forward forecasts (Figures \@ref(fig:boxplotrollingtourismappendix)) and fixed origin forecasts (Figure \@ref(fig:boxplottourismappendix)). Scaled errors are computed by subtracting the mean and dividing by the standard deviation. We see that the scaled error magnitudes are similar at different hierarchy/grouping levels.

```{r boxplotrollingtourismappendix, fig.align="center", fig.cap="Box plots of scaled forecast errors from reconciled and unreconciled ETS, ARIMA, and OLS rolling origin forecasts for tourism demand data. Panels are different hierarchy levels.", out.width="88%"}
## 1-step-ahead
error.tourism.scale <- error.tourism %>%
  mutate(Series = rep(1:6660, each = 24))
error.tourism.scale %>%
  filter(ForecastInterval == "1Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  group_by(id, Series) %>%
  mutate(error.scale = (error - mean(error)) / sd(error)) %>%
  ungroup() %>%
  ggplot(aes(x = id, y = error.scale, fill = id)) +
  geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 4, scales = "fixed") +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

```{r boxplottourismappendix, fig.align="center", fig.cap="Box plots of scaled forecast errors from reconciled and unreconciled ETS, ARIMA, and OLS fixed origin forecasts for tourism demand data. Panels are different hierarchy levels.", out.width="88%"}
## 24-step-ahead
error.tourism.scale %>%
  filter(ForecastInterval == "24Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  group_by(id, Series) %>%
  mutate(error.scale = (error - mean(error)) / sd(error)) %>%
  ungroup() %>%
  ggplot(aes(x = id, y = error.scale, fill = id)) +
  geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 4, scales = "fixed") +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

## Reconciled and unreconciled prediction intervals comparison

Figures \@ref(fig:boxplotrollingtourismappendix) and \@ref(fig:boxplottourismappendix) display boxplots of the scaled forecast errors for the tourism example. These plots are displayed for both rolling forward and multiple-step-ahead forecasts, respectively. We see that scaled error magnitudes are similar at different hierarchy/grouping levels.

```{r predintAAAVisrecunrec, fig.align="center", fig.cap="Comparing forecasts and prediction intervals of 'AAAVis' bottom-level series across methods. Left: rolling origin. Right: fixed origin.", out.width="100%", out.lenght="100%"}
### one series example series
p1 <- forecast.tourism.data  %>%
   filter(Series == "AAAVis", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = ARIMA.upper.unrec, ymin = ARIMA.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightblue", colour = "lightblue", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "ARIMA.unrec", size = 'ARIMA.unrec')) +
   geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("Rolling origin 1-step forecasts") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.unrec", "ARIMA.rec"),
                      values = c("black", "gray",  "blue"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.unrec", 'ARIMA.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw() +
    theme(legend.position="none")

p2 <- forecast.tourism.data  %>%
   filter(Series == "AAAVis", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = ETS.upper.unrec, ymin = ETS.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "ETS.unrec", size = 'ETS.unrec')) +
   geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "ETS.unrec", 'ETS.rec'),
                      values = c("black", "gray",  "green"))+
   scale_size_manual(breaks = c("Actual", "ETS.unrec", 'ETS.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw() +
    theme(legend.position="none")

p3 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = OLS.upper.unrec, ymin = OLS.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "OLS.unrec", size = 'OLS.unrec')) +
   geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "OLS.unrec", 'OLS.rec'),
                      values = c("black", "gray",  "red"))+
  scale_size_manual(breaks = c("Actual", "OLS.unrec", 'OLS.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw() +
    theme(legend.position="none")

p4 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = ARIMA.upper.unrec, ymin = ARIMA.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightblue", colour = "lightblue", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "ARIMA.unrec", size = 'ARIMA.unrec')) +
   geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("Fixed origin multi-step forecasts") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.unrec", "ARIMA.rec"),
                      values = c("black", "gray",  "blue"))+
scale_size_manual(breaks = c("Actual", "ARIMA.unrec", 'ARIMA.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw()

p5 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = ETS.upper.unrec, ymin = ETS.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "ETS.unrec", size = 'ETS.unrec')) +
   geom_line(aes(y = ETS.rec, colour = "ETS.rec", size= 'ETS.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "ETS.unrec", 'ETS.rec'),
                      values = c("black", "gray",  "green"))+
  scale_size_manual(breaks = c("Actual", "ETS.unrec", 'ETS.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw()

p6 <- forecast.tourism.data %>%
   filter(Series == "AAAVis", ForecastInterval == "24Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
   geom_ribbon(aes(x = date, ymax = OLS.upper.unrec, ymin = OLS.lower.unrec), fill = "gray83", colour = "gray83", alpha = .2, size = 0.5) +
   geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
   geom_line(aes(y = ETS.rec, colour = "OLS.unrec", size = 'OLS.unrec')) +
   geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
   geom_line() +
   xlab("Horizon") +
   ylab("Count") +
   ggtitle("") +
   scale_colour_manual("Method",
                      breaks = c("Actual", "OLS.unrec", 'OLS.rec'),
                      values = c("black", "gray",  "red"))+
   scale_size_manual(breaks = c("Actual", "OLS.unrec", 'OLS.rec'),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
   theme_bw()


grid.arrange(p1, p4, p2, p5, p3, p6,  nrow=3, widths=c(2.1/5,2.9/5))
```



# Wikipedia pageviews: Grouped structure


We illustrate our method on another real dataset. The Wikipedia dataset includes 12 months of daily data (2016-06-01 to 2017-06-29) on Wikipedia pageviews for the most popular social networks articles [@ashouri2018]. This dataset is noisier than the Australian monthly tourism data, making forecasting more challenging. The data has a grouped structure with the following attributes (see Table \@ref(tab:wikipediagroupingstructure)):

  * *Agent*: Spider, User;
  * *Access*: Desktop, Mobile app, Mobile web;
  * *Language*: en (English), de (German), es (Spanish), zh (Chinese); and
  * *Purpose*: Blogging related, Business, Gaming, General purpose, Life style, Photo sharing, Reunion, Travel, Video.

Figure \@ref(fig:wikigroupstructure) hows one possible hierarchy for this dataset, but the order of the hierarchy can be switched.

```{r wikigroupstructure, echo=FALSE, out.width = "500px", out.height= "250px", fig.align="center", fig.cap="One possible hierarchical structure for the Wikipedia pageviews dataset."}
knitr::include_graphics("Paper-Figures/Wiki_group_structure.jpg")
```

```{r Readwikipedia, results='hide'}
fc.rolling.wiki.OLS <- readr::read_csv("Wikipedia_data/fc.rolling.wiki.OLS.csv") %>%
  mutate(ForecastInterval = '1Step')
fc.fix.wiki.OLS <- readr::read_csv("Wikipedia_data/fc.fix.wiki.OLS.csv") %>%
  mutate(ForecastInterval = '28Step')

fc.wiki.OLS <- bind_rows(fc.rolling.wiki.OLS, fc.fix.wiki.OLS)

fc.wiki.OLS <- bind_rows(fc.wiki.OLS %>%
                                   filter(Series == 'Total') %>%
                                   mutate (Level = 'Total'),
                         fc.wiki.OLS %>% filter(Series %in% c('Access/desktop','Access/mobilea', 
                                                                              'Access/mobilea')) %>%
                                   mutate (Level = 'Access'), 
                         fc.wiki.OLS %>% filter(Series %in% c('Agent/sp','Agent/us')) %>%
                                   mutate (Level = 'Agent'), 
                         fc.wiki.OLS %>% filter(Series %in% c('Language/de','Language/es', 
                                                                              'Language/en', 'Language/zh')) %>%
                                   mutate (Level = 'Language'), 
                         fc.wiki.OLS %>% filter(Series %in% c('Purpose/Blo', 'Purpose/Bus', 'Purpose/Gam',
                                                                              'Purpose/Gen', 'Purpose/Lif', 'Purpose/Pho', 
                                                                              'Purpose/Reu', 'Purpose/Tra', 'Purpose/Vid')) %>%
                                   mutate (Level = 'Purpose'), 
                         fc.wiki.OLS %>% filter(grepl('Access x Purpose/', Series)) %>%
                                   mutate (Level = 'Access x Purpose'), 
                         fc.wiki.OLS %>% filter(grepl('Agent x Purpose/', Series)) %>%
                                   mutate (Level = 'Agent x Purpose'), 
                         fc.wiki.OLS %>% filter(grepl('Language x Purpose/', Series)) %>%
                                   mutate (Level = 'Language x Purpose'),
                         fc.wiki.OLS %>% filter(grepl('Access x Agent/', Series)) %>%
                                   mutate (Level = 'Access x Agent'),
                         fc.wiki.OLS %>% filter(grepl('Access x Language/', Series)) %>%
                                   mutate (Level = 'Access x Language'),
                         fc.wiki.OLS %>% filter(grepl('Agent x Language/', Series)) %>%
                                   mutate (Level = 'Agent x Language'),
                         fc.wiki.OLS %>% filter( !grepl('Access', Series) & !grepl('Agent', Series) & 
                                                                   !grepl('Language', Series) & !grepl('Purpose', Series) & 
                                                                   !grepl('Total', Series)) %>%
                                   mutate (Level = 'Bottom level'))

error.wiki.OLS <-bind_rows( dplyr::select ( fc.wiki.OLS, error = error.rec, Level, ForecastInterval) %>% 
                                         mutate(Rec = 'rec') %>%
                                         mutate(Method = 'OLS') , dplyr::select ( fc.wiki.OLS, error = error.unrec, Level, ForecastInterval) %>% 
                                         mutate(Rec = 'unrec') %>%
                                         mutate(Method = 'OLS'))

ets.arima.rolling.wiki <- readr::read_csv("Wikipedia_data/fc.rolling.wiki.ets.arima.csv") %>%
  select(-new_index) %>%
  mutate(ForecastInterval = '1Step')

ets.arima.fix.wiki <- readr::read_csv("Wikipedia_data/fc.fix.wiki.ets.arima.csv") %>%
  mutate(ForecastInterval = '28Step')

ets.arima.wiki <- bind_rows(ets.arima.rolling.wiki, ets.arima.fix.wiki)

ets.arima.wiki <- bind_rows( ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent == '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Total'), 
                             ets.arima.wiki %>%
                                       filter(Access != '<aggregated>', Agent == '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Access'), 
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent != '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Agent') ,
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent == '<aggregated>', Language != '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Language') ,
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent == '<aggregated>', Language == '<aggregated>', Purpose != '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Purpose') ,
                             ets.arima.wiki %>%
                                       filter(Access != '<aggregated>', Agent != '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Access x Agent') ,
                             ets.arima.wiki %>%
                                       filter(Access != '<aggregated>', Agent == '<aggregated>', Language != '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Access x Language') ,
                             ets.arima.wiki %>%
                                       filter(Access != '<aggregated>', Agent == '<aggregated>', Language == '<aggregated>', Purpose != '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Access x Purpose') ,
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent != '<aggregated>', Language != '<aggregated>', Purpose == '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Agent x Language') ,
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent != '<aggregated>', Language == '<aggregated>', Purpose != '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Agent x Purpose') ,
                             ets.arima.wiki %>%
                                       filter(Access == '<aggregated>', Agent == '<aggregated>', Language != '<aggregated>', Purpose != '<aggregated>'
                                              , Article == '<aggregated>') %>%
                                       mutate (Level = 'Language x Purpose') ,
                             ets.arima.wiki %>%
                                       filter(Access != '<aggregated>', Agent != '<aggregated>', Language != '<aggregated>', Purpose != '<aggregated>'
                                              , Article != '<aggregated>') %>%
                                       mutate (Level = 'Bottom level') )

ets.arima.wiki <- bind_rows(ets.arima.wiki %>% 
                                      filter(.model %in% c('arima', 'ets')) %>%
                                      mutate(Rec = 'unrec'), 
                                    ets.arima.wiki %>% 
                                      filter(.model %in% c('arima_adjusted', 'ets_adjusted')) %>%
                                      mutate(Rec = 'rec'))

ets.arima.wiki <- bind_rows(ets.arima.wiki %>% 
                                      filter(.model %in% c('arima', 'arima_adjusted')) %>%
                                      mutate(Method = 'ARIMA'), 
                            ets.arima.wiki %>% 
                                      filter(.model %in% c('ets', 'ets_adjusted')) %>%
                                      mutate(Method = 'ETS'))


error.wiki <- bind_rows(
  dplyr::select ( ets.arima.wiki, error , Level, Rec, Method, ForecastInterval), error.wiki.OLS) %>%
  filter(Level %in% c("Total", "Access", "Agent", "Language", "Purpose", "Bottom level")) %>%
  mutate(
    facet = factor(Level, levels = c("Total", "Access", "Agent", "Language", "Purpose", "Bottom level")))


# Read csv files for wiki results
rmse <- error.wiki %>%
  group_by(Rec, Method, facet, ForecastInterval) %>%
  summarise(
    rmse = sqrt(mean(error^2))
  ) %>%
  spread(value = rmse, key = Method) %>%
  ungroup() %>%
  select(Rec, ForecastInterval, facet, ETS, ARIMA, OLS) %>%
  mutate(
    facet = str_replace(facet, "level([0-9])", "facet \\1")
  )


## forecast results
m <- c(746, 450, 420, 576, 596, 631, 611, 514, 411, 452, 582, 927, 1106, 1313, 926, 422, 442, 619, 630, 605, 585, 606, 414, 405, 623, 614, 644, 653) 
forecast.wiki.OLS <- bind_rows(fc.wiki.OLS %>% filter(Series == 'Total'), fc.wiki.OLS %>% filter(Series == 'desktopusenPho') %>% filter(Actual%in% m))
forecast.wiki.OLS <- forecast.wiki.OLS %>% select(-error.rec, -error.unrec, -Level) 

arima.unrec <- bind_rows( ets.arima.wiki %>% filter(.model == 'arima' , Access == '<aggregated>', 
                                                               Agent == '<aggregated>', Language == '<aggregated>', 
                                                            Purpose == '<aggregated>', 
                                                               Article == '<aggregated>') %>%
                            select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`),
                          ets.arima.wiki %>% filter(.model == 'arima' , Access == 'desktop', 
                                                               Agent == 'us', Language == 'en', Purpose == 'Pho', Article == '21') %>%
                            select(ARIMA.unrec = .mean, ARIMA.lower.unrec = `95%_lower`, ARIMA.upper.unrec = `95%_upper`))

arima.rec <- bind_rows( ets.arima.wiki %>% filter(.model == 'arima_adjusted' , Access == '<aggregated>', 
                                                          Agent == '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>',
                                                          Article == '<aggregated>') %>%
                          select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`),
                        ets.arima.wiki %>% filter(.model == 'arima_adjusted' , Access == 'desktop', 
                                                          Agent == 'us', Language == 'en', Purpose == 'Pho', Article == '21') %>%
                          select(ARIMA.rec = .mean, ARIMA.lower.rec = `95%_lower`, ARIMA.upper.rec = `95%_upper`))

ets.unrec <- bind_rows( ets.arima.wiki %>% filter(.model == 'ets' , Access == '<aggregated>', 
                                                          Agent == '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>',
                                                          Article == '<aggregated>') %>%
                          select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`),
                        ets.arima.wiki %>% filter(.model == 'ets' , Access == 'desktop', 
                                                          Agent == 'us', Language == 'en', Purpose == 'Pho', Article == '21') %>%
                          select(ETS.unrec = .mean, ETS.lower.unrec = `95%_lower`, ETS.upper.unrec = `95%_upper`))

ets.rec <- bind_rows( ets.arima.wiki %>% filter(.model == 'ets_adjusted' , Access == '<aggregated>', 
                                                        Agent == '<aggregated>', Language == '<aggregated>', Purpose == '<aggregated>',
                                                        Article == '<aggregated>') %>%
                        select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`),
                      ets.arima.wiki %>% filter(.model == 'ets_adjusted' , Access == 'desktop', 
                                                        Agent == 'us', Language == 'en', Purpose == 'Pho', Article == '21') %>%
                        select(ETS.rec = .mean, ETS.lower.rec = `95%_lower`, ETS.upper.rec = `95%_upper`))
forecast.wiki.data <- bind_cols (forecast.wiki.OLS, arima.unrec, arima.rec, ets.unrec, ets.rec)

forecast.wiki <- forecast.wiki.data %>%
  select(-OLS.lower.rec, -OLS.upper.rec, -OLS.lower.unrec, -OLS.upper.unrec,
         -ARIMA.lower.rec, -ARIMA.upper.rec, -ARIMA.lower.unrec, -ARIMA.upper.unrec,
         -ETS.lower.rec, -ETS.upper.rec, -ETS.lower.unrec, -ETS.upper.unrec) %>%
  gather(-Series, -ForecastInterval, -date, key = "Method", value = "Count") %>%
  mutate(
    Rec = str_extract(Method, "[a-z]*$"),
    Rec = if_else(Rec == "ctual", "Actual", Rec),
    Model = str_extract(Method, "^[A-Z]*"),
    Model = if_else(Model == "A", "Actual", Model)
  )

```

```{r wikipediagroupingstructure, echo=FALSE, message = FALSE, results='asis'}
tibble(
  group1 = c("Total", "", "Access", "", "", "Agent", "", "", "", "Language", "", "", ""),
  series1 = c("", "1. Social Network", "", "2. Desktop", "3. Mobile app", "", "4.  Mobile web", "5. Spider", "6. User", "", "7. en (English)", "8. de (German)", "9. es (Spanish)"),
  group2 = c("Language", "", "Purpose", "", "", "", "", "", "", "", "", "", ""),
  series2 = c("", "10. zh (Chinese)", "", "11. Blogging related", "12. Business", "13. Gaming", "14. General purpose", "15. Life style", "16. Photo sharing", "17. Reunion", "18. Travel", "19. Video", "")
) %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Social networking Wikipedia article grouping structure",
    col.names = c("Grouping", "Series", "Grouping", "Series")
  ) %>%
  kable_styling(position = "center")
```

We consider the main aggregation factors and their two-way combinations^[First, while we present results for single groups, we applied all the two-way combinations in the reconciliation step. Second, there are four more 3-way aggregation combinations that we do not include: Agent \(\times\) Access \(\times\) Language, Agent \(\times\) Access \(\times\) Purpose, Agent \(\times\) Language \(\times\) Purpose, and Access \(\times\) Language \(\times\) Purpose. Including these four additional aggregations might slightly improve the results but for simplicity, we excluded them.]. The final dataset includes 913 time series, each with length 394. Table \@ref(tab:wikidivision) shows the group structure's different levels and the number of series in each level.

```{r wikidivision, results='asis'}
groups <- tribble(
  ~Aggregation_Level, ~Series,
  "Total pageviews", 1,
  "Access", 3,
  "Agent", 2,
  "Language", 4,
  "Purpose", 9,
  "Access x Agent", 5,
  "Access x Language", 12,
  "Access x Purpose", 27,
  "Agent x Language", 8,
  "Agent x Purpose", 18,
  "Language x Purpose", 33,
  "Bottom level", 913,
  "Total", NA
)
groups[13, 2] <- sum(groups[1:12, 2])
groups %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    caption = "Number of Wikipedia pageviews series at each aggregation level."
  ) %>%
  kable_styling(
    position = "center",
    latex_options = "hold_position"
  ) %>%
  row_spec(row = 12, hline_after = TRUE)
```

<!--\todo[inline]{What about all the cross-products such as the aggregates of Access*Agent, or Language*Purpose, or Agent*Access*Language, etc?}-->

For this daily dataset, in the OLS forecasting model we include in the predictor matrix the following terms: a quadratic trend, 6 seasonal dummies, and lags 1 and 7 for rolling and fixed origin models. We partitioned the data into training and test sets. We used the last 28 days for our test set and the rest for the training set.

Tables \@ref(tab:wikipediadataresulrolling) and \@ref(tab:wikipediadataresultRMSE) show the RMSE results. Although these time series are noisier, we still get acceptable results for the OLS forecasting model compared with ETS and ARIMA. In this case, we get similar results with and without the reconciliation step.

```{r wikipediadataresulrolling, results='asis', dependson="Readwikipedia"}
bind_cols(
  filter(rmse, ForecastInterval == "1Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  filter(rmse, ForecastInterval == "1Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Mean RMSE for ETS, ARIMA and OLS with and without reconciliation - Rolling origin - Wikipedia dataset",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 3, "Reconciled" = 3)) %>%
  kable_styling(position = "center", latex_options = "hold_position")
```

```{r wikipediadataresultRMSE, results='asis', dependson="Readwikipedia"}
bind_cols(
  filter(rmse, ForecastInterval == "28Step", Rec == "unrec") %>% select(-Rec, -ForecastInterval),
  filter(rmse, ForecastInterval == "28Step", Rec == "rec") %>% select(-Rec, -ForecastInterval, -facet)
) %>%
  kable(
    booktabs = TRUE,
    format = "latex",
    digits = 0,
    linesep = "",
    caption = "Mean RMSE for ETS, ARIMA and OLS with and without reconciliation - Fixed origin - Wikipedia dataset",
    col.names = c("Level", rep(c("ETS", "ARIMA", "OLS"), 2))
  ) %>%
  add_header_above(c("", "Unreconciled" = 3, "Reconciled" = 3)) %>%
  kable_styling(position = "center")
```

```{r boxplotrollingwiki, fig.align="center", fig.cap="Box plots of forecast errors for reconciled and unreconciled ETS, ARIMA and OLS methods at each hierarchy level for rolling origin forecasts of Wikipedia pageviews.", out.width="100%"}
## 1-step-ahead
error.wiki %>%
  filter(ForecastInterval == "1Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  ggplot(aes(x = id, y = error, fill = id)) +
  stat_summary(fun.data = boxplot.stat, geom = "boxplot", alpha = 0.5) +
  # geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 3, scales = "free_y") +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

```{r boxplotwiki, fig.align="center", fig.cap="Box plots of forecast errors for reconciled and unreconciled ETS, ARIMA and OLS methods at each hierarchy level for fixed origin forecasts of Wikipedia pageviews.", out.width="100%"}
## 28-step-ahead
error.wiki %>%
  filter(ForecastInterval == "28Step") %>%
  mutate(id = factor(paste(Method, Rec, sep = "."),
    levels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec"),
    labels = c("ETS.rec", "ETS.unrec", "ARIMA.rec", "ARIMA.unrec", "OLS.rec", "OLS.unrec")
  )) %>%
  ggplot(aes(x = id, y = error, fill = id)) +
  stat_summary(fun.data = boxplot.stat, geom = "boxplot", alpha = 0.5) +
  # geom_boxplot(alpha = 0.5) +
  xlab("Method") +
  ylab("Forecast errors") +
  facet_wrap(~facet, ncol = 3, scales = "free_y") +
  guides(fill = guide_legend(nrow = 1, bycol = TRUE)) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "ETS.rec" = "green",
    "ETS.unrec" = "lightgreen",
    "ARIMA.rec" = "blue",
    "ARIMA.unrec" = "lightblue",
    "OLS.rec" = "pink4",
    "OLS.unrec" = "pink"
  )) +
  theme_light() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "none"
  )
```

```{r forecstrolling28wikitotal, fig.align="center", fig.cap="Comparing reconciled and unreconciled ETS, ARIMA and OLS rolling and fixed origin forecasts for Wikipedia pageviews 'Total' Series.", out.width="100%"}
### one of the bottom level series
ylim <- forecast.wiki %>%
  filter(Series == "Total") %>%
  pull(Count) %>%
  range()
h1 <- forecast.wiki %>%
  filter(Series == "Total", ForecastInterval == "1Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
    scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      OLSX = "yellow",
      ARIMAX = "orchid1",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(linetype = FALSE)
h2 <- forecast.wiki %>%
  filter(Series == "Total", ForecastInterval == "28Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
    scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      OLSX = "yellow",
      ARIMAX = "orchid1",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(color = FALSE)
h1 / h2
```

```{r forecstrolling28wiki, fig.align="center", fig.cap="Comparing reconciled and unreconciled ETS, ARIMA and OLS rolling and fixed origin forecasts for Wikipedia pageviews 'desktopusenPho21' bottom-level series", out.width="100%"}
### one of the bottom level series
ylim <- forecast.wiki %>%
  filter(Series == "desktopusenPho") %>%
  pull(Count) %>%
  range()
f1 <- forecast.wiki %>%
  filter(Series == "desktopusenPho", ForecastInterval == "1Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle(" Rolling origin 1-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
    scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      OLSX = "yellow",
      ARIMAX = "orchid1",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(linetype = FALSE)

f2 <- forecast.wiki %>%
  filter(Series == "desktopusenPho", ForecastInterval == "28Step") %>%
  # rename(Reconciled = Rec) %>%
  ggplot(aes(x = date, y = Count, colour = Model, linetype = Rec, size = Model)) +
  geom_line() +
  ylim(ylim) +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_linetype_manual(name = "Reconciled", values = c(Actual = "solid", rec = "dashed", unrec = "dotted")) +
      scale_size_manual(values = c(Actual = 0.8, ETS = 0.5, ARIMA = 0.5, OLS = 0.5), guide = "none") +
  scale_color_manual(
    name = "Series",
    values = c(
      ETS = "green",
      ARIMA = "blue",
      OLS = "red",
      OLSX = "yellow",
      ARIMAX = "orchid1",
      Actual = "black"
    )
  ) +
  theme_bw() +
  guides(color = FALSE)
f1 / f2
```


```{r predinttotalwiki, fig.align="center", fig.cap="The actual test set for the 'Total series' compared to the forecasts from reconciled and unreconciled OLS methods with prediction interval for rolling and fixed origin Wikipedia pageviews.", out.width="100%"}
### total series
p1 <- forecast.wiki.data %>%
   filter(Series == "Total", ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, size = 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green" , "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()+
    theme(legend.position="none")

p2 <- forecast.wiki.data %>%
   filter(Series == "Total", ForecastInterval == "28Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  #geom_ribbon(aes(x = date, ymax = ETS.rec.upper, ymin = ETS.rec.lower), fill = "lightgreen", colour = "lightgreen", alpha = .2) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  #geom_line(aes(y = ETS.rec, colour = "ETS.rec")) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",   "OLS.rec"),
                      values = c("black", "blue",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "OLS.rec"),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
  theme_bw()

p1 / p2
```

```{r predintdesktopus, fig.align="center", fig.cap="The actual test set for the 'desktopusenPho21' bottom level series compared to the forecasts from reconciled and unreconciled OLS methods with prediction interval for rolling and fixed origin  Wikipedia pageviews.", out.width="100%"}
## one series example series
p1 <- forecast.wiki.data %>%
   filter(Series == "desktopusenPho", Actual %in% m , ForecastInterval == "1Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  geom_ribbon(aes(x = date, ymax = ETS.upper.rec, ymin = ETS.lower.rec), fill = "lightgreen", colour = "lightgreen", alpha = .2, szie= 0.5) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  geom_line(aes(y = ETS.rec, colour = "ETS.rec", size = 'ETS.rec')) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Rolling origin 1-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
                      values = c("black", "blue", "green" ,"red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec",  "ETS.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5,  0.5), guide = "none") +
  theme_bw()+
    theme(legend.position="none")

p2 <- forecast.wiki.data %>%
   filter(Series == "desktopusenPho", Actual %in% m, ForecastInterval == "28Step") %>%
   ggplot(aes(x = date, y = Actual, colour = "Actual", size = 'Actual')) +
  geom_ribbon(aes(x = date, ymax = ARIMA.upper.rec, ymin = ARIMA.lower.rec), fill = "lightskyblue", colour = "lightskyblue", alpha = .2, size = 0.5) +
  #geom_ribbon(aes(x = date, ymax = ETS.rec.upper, ymin = ETS.rec.lower), fill = "lightgreen", colour = "lightgreen", alpha = .2) +
 geom_ribbon(aes(x = date, ymax = OLS.upper.rec, ymin = OLS.lower.rec), fill = "pink", colour = "pink", alpha = .2, size = 0.5)  +
  geom_line(aes(y = ARIMA.rec, colour = "ARIMA.rec", size = 'ARIMA.rec')) +
  #geom_line(aes(y = ETS.rec, colour = "ETS.rec")) +
  geom_line(aes(y = OLS.rec, colour = "OLS.rec", size = 'OLS.rec')) +
  geom_line() +
  xlab("Horizon") +
  ylab("Count") +
  ggtitle("Fixed origin multi-step forecasts") +
  scale_colour_manual("Method",
                      breaks = c("Actual", "ARIMA.rec",  "OLS.rec"),
                      values = c("black", "blue",  "red"))+
   scale_size_manual(breaks = c("Actual", "ARIMA.rec", "OLS.rec"),
        values = c( 0.8, 0.5,  0.5), guide = "none") +
  theme_bw()

p1 / p2
```


Figures \@ref(fig:boxplotrollingwiki) and \@ref(fig:boxplotwiki) display the forecast error box plot. These plots are for rolling and fixed origin forecasts over 28 days in each level of grouping. We see that the error distribution is similar in all levels across the different methods. The only exception is the Total series, where ETS performs significantly better than ARIMA and OLS. We also note that the reconciliation is less effective. As in the tourism example, in higher levels series have higher counts and therefore their error magnitudes are larger.

In Figures \@ref(fig:forecstrolling28wikitotal) and \@ref(fig:forecstrolling28wiki), we display results for the total and one of the bottom-level series, 'desktopusenPho21' (desktop-user-english-photo sharing). The plot shows rolling and fixed origin forecast results over the 28 day test set for ETS, ARIMA and OLS, with (dashed lines) and without (dotted lines) applying reconciliation. We see that the OLS forecasting model performs close to the other two methods, and reconciliation improves the forecasts. In Figures \@ref(fig:predinttotalwiki) and \@ref(fig:predintdesktopus), we present the reconciled prediction intervals for the total and 'desktopusenPho21' series. In the fixed origin results, the ETS prediction intervals were much wider and hence have been omitted to allow for easier visualization of the results.

Table \@ref(tab:wikipediadatacomputationtime) presents the computation times for all three methods. ETS and ARIMA are clearly much more computationally heavy compared with OLS. As in the Australian tourism dataset, running reconciliation does not have much effect on computation time.


```{r wikipediadatacomputationtime, echo=FALSE, results='asis',message = FALSE}
tibble(
  model = c("ETS", "ARIMA", "OLS"),
  #rolling_unreconciled = c(, , 116.27),
  rolling_reconciled = c(22592, 20016, 116),
  #fixed_unreconciled = c(, , 61.33),
  fixed_reconciled = c(882, 1682, 61)
) %>%
  kable(
    align = c("l", rep("r",2)),
    format = "latex",
    booktabs = TRUE,
    linesep = "",
    digits = 2,
    caption = "Computation time (seconds) for ETS, ARIMA and OLS with reconciliation - Rolling and fixed origin forecasts - Wikipedia dataset",
    col.names = c("", "Rolling origin", "Fixed origin")
  ) %>%
  # column_spec(1:3, width = "3cm") %>%
  # add_header_above(c("", "Rolling origin" = 2, "Fixed origin" = 2)) %>%
  # add_header_above(c("", "Computation time (secs)" = 2)) %>%
  kable_styling(position = "center", full_width = F)
```


\clearpage