longtable volume+number+eid [cbx@textcite]nameyeardelim

A novel forecast linear augmented projection (FLAP) method is introduced,
which reduces the forecast error variance of any unbiased multivariate
forecast without introducing bias. The method first constructs new component
series which are linear combinations of the original series. Forecasts are
then generated for both the original and component series. Finally, the full
vector of forecasts is projected onto a linear subspace where the constraints
implied by the combination weights hold. It is proven that the trace of the
forecast error variance is non-increasing with the number of components, and
mild conditions are established for which it is strictly decreasing. It is
also shown that the proposed method achieves maximum forecast error variance
reduction among linear projection methods. The theoretical results are
validated through simulations and two empirical applications based on
Australian tourism and FRED-MD data. Notably, using FLAP with Principal
Component Analysis (PCA) to construct the new series leads to substantial
forecast error variance reduction.

## 1 Introduction

We will show that any multivariate forecasting method can be improved through
a post-processing framework involving linear combinations of the original
series. We call this FLAP: Forecast Linear Augmented Projection. This has
immediate implications for multivariate forecasting in all disciplines
including macroeconomics [macroreview] and finance [financereview].

Our FLAP post-processing framework: (i) augments the data by constructing new
series that are linear combinations of the original series; (ii) forecasts
both the original and new series; and (iii) recovers a new set of forecasts
for the original series via projections. We prove that the method reduces the
forecast error variance of the original series in a way that is agnostic both
with respect to the weights of the linear combinations used at step (i) and
with respect to the model used to generate forecasts at step (ii). The model
is inspired by the forecast reconciliation literature [AthEtAl2023a] whereby
forecasts are adjusted to cohere with known linear constraints. In contrast to
that literature, the FLAP method focuses on multivariate forecasting where
such constraints are not present. Indeed, the method need not only be applied
to forecasting problems, but multivariate predictions in general.

It may appear puzzling that forecast accuracy can be improved, not by
introducing any new information, but by simply taking linear combinations of
existing time series. To give an intuition into how this puzzle can be
resolved, we consider a toy example of two series $y_{1}$ and $y_{2}$ that are
of concern to the forecaster, and two linear combinations or _components_ of
the series $c_{1}=0.5z_{1}+0.5z_{2}$ and $c_{2}=0.5z_{1}-0.5z_{2}$. Denote by
$\hat{y}_{1}$, $\hat{y}_{2}$, $\hat{c}_{1}$, $\hat{c}_{2}$ any forecasts of
these original series and components, which we collectively refer to as _base
forecasts_. The base forecasts may be generated by univariate methods,
multivariate methods, or even based on expert judgement. When considering
$y_{1}$, there is both a _direct_ forecast $\hat{y}_{1}$ and _indirect_
forecast $\hat{c}_{1}+\hat{c}_{2}$, similarly for $y_{2}$ the direct forecast
is $\hat{y}_{2}$ and the indirect forecast $\hat{c}_{1}-\hat{c}_{2}$.111This
argument, as well as the terminology direct and indirect forecasts, is
inspired by [HolEtAl2021] who discuss this in the forecast reconciliation
setting. With the exception of some pathological cases, the direct and
indirect forecasts for the same variable will not, in general, be equal.
Therefore forecast accuracy can be improved by combining direct and indirect
forecasts, something implicitly achieved by the proposed FLAP method. The
puzzle is thus resolved; while no new information is created at the data
augmentation step, there is new information embedded into forecasts of the
augmented series, which can be leveraged via model combination (see
[WanEtAl2023] for a review of forecast combination). Something obscured by the
simple toy example is the way our FLAP method differs from the usual forecast
combination methods, in particular our combinations are potentially non-convex
since they are obtained via projections, in a way that we now elaborate upon.

More formally and more generally, the FLAP method considers a vector of
original series $\mathbf{y}\in\mathbb{R}^{m}$ and a vector of components
$\mathbf{c}\in\mathbb{R}^{p}$. While
$\left(\mathbf{y}^{\prime},\mathbf{c}^{\prime}\right)^{\prime}$ is a
$p+m$-vector, the construction of components as linear combinations of the
original series implies that
$\left(\mathbf{y}^{\prime},\mathbf{c}^{\prime}\right)^{\prime}$ lies on a
linear subspace of at most dimension $m$. The corresponding vector of
forecasts
$\left(\hat{\mathbf{y}}^{\prime},\hat{\mathbf{c}}^{\prime}\right)^{\prime}$
will, in general, have support on $\mathbb{R}^{m+p}$. FLAP projects
$\left(\hat{\mathbf{y}}^{\prime},\hat{\mathbf{c}}^{\prime}\right)^{\prime}$
onto the $m$-dimensional linear subspace on which
$\left(\mathbf{y}^{\prime},\mathbf{c}^{\prime}\right)^{\prime}$ has support.
The setup of this problem bears similarities to the well known problem of
forecast reconciliation where [PanEtAl2021] provide similar geometric
intuition, while [WicEtAl2019], [AthEtAl2017], and [DiGir2023] have all shown
that reconciliation can reduce forecast error variance theoretically and
empirically. However, we note that these papers establish that reconciliation
improves forecast accuracy for the hierarchy _as a whole_. In the general
multivariate setting that we consider, this would imply improvements in
forecast accuracy for $\mathbf{y}$ and $\mathbf{c}$ taken together. This poses
a problem if improvements in forecast accuracy for $\mathbf{c}$ could be
offset by a deterioration in forecast accuracy for $\mathbf{y}$, since the
former are not of interest in and of themselves. A key insight we make in this
paper, is that reductions in forecast error variance accrue even for a subset
of variables in the hierarchy. It is this contribution that allows us to
propose a method that goes beyond the case where time series adhere to linear
constraints, and that instead applies to the more general setting.

While the theoretical results apply for any linear combinations of the
original series, in practice we propose to augment the data with principal
components (PCs). When doing so, the FLAP method bears a resemblance to
Dynamic Factor Models (DFMs), specifically those common in macroeconomic
forecasting [StoWat2002a, StoWat2002, StoWat2012], their extensions in the
machine learning literature [DeEtAl2019], as well as the factor augmented VAR
[BerEtAl2005]. The factor models assume that the multivariate time series
possesses common components and the dynamics of the observed series are
governed by the dynamics of these unobserved factors, often assumed to follow
some parametric model. In contrast, FLAP is a post-forecasting step; indeed,
forecasts can even be made using a DFM and then further improved by
implementing the FLAP method, something we demonstrate in Section 3 and
Section 4.

In the sense that forecast accuracy can be improved without any new
information, FLAP has parallels with bootstrap aggregation or “bagging”
[Bre1996, BerEtAl2016]. Bagging can reduce prediction variance without
increasing bias [HasEtAl2003], by mitigating model uncertainty [PetEtAl2018],
and does so without introducing any new data, but rather resampled versions of
the existing data. Our FLAP method also reduces forecast error variance
without introducing new data, but using linear combinations of the existing
data, rather than bootstrapping. The FLAP method improves forecast accuracy
through forecast combination and data augmentation. Other methods that also do
this include the theta method [AssNik2000], temporal aggregation [KouEtAl2014,
AthEtAl2017], forecasting with sub-seasonal series [LiEtAl2022c, FOSS,] and
forecast combination with multiple starting points [DisPet2015]; a review of
all these methods can be found in [PetSpi2021] who refer to them as using “the
wisdom of data”. Our FLAP method is distinct in that it aims to exploit
information in the data with a focus on linear combinations of multivariate
series.

The remainder of the paper is structured as follows. In Section 2, we propose
the FLAP method, and highlight its theoretical properties and associated
estimation methods. In Section 3, we present a simulation example
demonstrating its performance and discuss the implications for sources of
uncertainty. Section 4 examines the performance of FLAP in two empirical
applications: forecasting Australian domestic tourism and forecasting
macroeconomic variables in the FRED-MD data set. Section 5 concludes with some
thoughts on future research directions. The methods introduced in this paper
are implemented in the flap package for R [flap]. This paper is fully
reproducible with code and documentation provided at
https://github.com/FinYang/paper-forecast-projection.

## 2 Forecast Linear Augmented Projection (FLAP)

### 2.1 Method and theoretical properties

In the following, all vectors and matrices are denoted in bold font. We use
$\bm{I}_{n}$ to denote the $n\times n$ identity matrix, and $\bm{O}_{n\times
k}$ to denote the $n\times k$ zero matrix.

Let $\bm{y}_{t}\in\mathbf{R}^{m}$ be a vector of $m$ observed time series we
are interested in forecasting. The FLAP method involves three steps:

  1. 1.

_Form components_. Form $\bm{c}_{t}=\bm{\Phi}\bm{y}_{t}\in\mathbf{R}^{p}$, a
vector of $p$ linear combinations of $\bm{y}_{t}$ at time $t$, where
$\bm{\Phi}\in\mathbf{R}^{p\times m}$. We call $\bm{c}_{t}$ the components of
$\bm{y}_{t}$ and the component weights $\bm{\Phi}$ are known in the sense that
they are chosen by the user of FLAP. Let
$\bm{z}_{t}=\big{[}\bm{y}_{t}^{\prime},\bm{c}^{\prime}_{t}\big{]}^{\prime}$ be
the concatenation of series $\bm{y}_{t}$ and components $\bm{c}_{t}$. We note
that $\bm{z}_{t}$ will be constrained in the sense that
$\bm{C}\bm{z}_{t}=\bm{c}_{t}-\bm{\Phi}\bm{y}_{t}=\bm{0}$ for any $t$ where
$\bm{C}=\big{[}-\bm{\Phi}~{}~{}~{}\bm{I}_{p}\big{]}$ is referred to as the
constraint matrix.

  2. 2.

_Generate forecasts_. Denote as $\hat{\bm{z}}_{t+h}$ the $h$-step-ahead base
forecast of $\bm{z}_{t}$. The method used to generate forecasts is again
selected by the user. This can be univariate or multivariate. In the setting
where $\bm{z}_{t}$ are not time series but cross sectional data, any
prediction method can be used. In general, the constraints that hold for
$\bm{z}_{t}$ will not hold for $\hat{\bm{z}}_{t+h}$, i.e
$\bm{C}\hat{\bm{z}}_{t+h}\neq\bm{0}$

  3. 3.

_Project the base forecasts_. Let $\tilde{\bm{z}}_{t+h}$ be a set of projected
forecasts such that,

${\tilde{\bm{z}}_{t+h}=\bm{M}\hat{\bm{z}}_{t+h}}$ (1)

with projection matrix

${\bm{M}=\bm{I}_{m+p}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C},}$
(2)

where $\operatorname{Var}(\bm{z}_{t+h}-\hat{\bm{z}}_{t+h})=\bm{W}_{h}$ is the
forecast error covariance matrix. For the proofs of this section, we will
assume that $\bm{W}_{h}$ is known, in practice a plug-in estimate can be used
that will be discussed in Section 2.5.

In practice we are interested in forecasts of $\bm{y}_{t}$ and not the full
vector $\bm{z}_{t}$. We now introduce some notation to handle this issue.
Define the selection matrix
$\bm{J}_{n,k}=\big{[}\bm{I}_{n}~{}~{}~{}\bm{O}_{n\times k}\big{]}$, so that
$\bm{J}_{n,k}\bm{A}$ selects the first $n$ rows of a matrix $\bm{A}$. Let
$\hat{\bm{y}}_{t+h}$ and $\tilde{\bm{y}}_{t+h}$ denote the first $m$ elements
of $\hat{\bm{z}}_{t+h}$ and $\tilde{\bm{z}}_{t+h}$, comprising the base and
projected forecasts of $\bm{y}_{t}$ respectively. Similarly, let
$\hat{\bm{c}}_{t+h}$ and $\tilde{\bm{c}}_{t+h}$ denote the last $p$ elements
of $\hat{\bm{z}}_{t+h}$ and $\tilde{\bm{z}}_{t+h}$, comprising the base and
projected forecasts of $\bm{c}_{t}$ respectively. Then the projected forecast
of $\bm{y}_{t}$ can be found by

${\tilde{\bm{y}}_{t+h}=\bm{J}\tilde{\bm{z}}_{t+h}=\bm{J}\bm{M}\hat{\bm{z}}_{t+h},}$
(3)

where $\bm{J}=\bm{J}_{m,p}$.

We now present some theoretical results regarding the FLAP method, with proofs
provided in the Appendix. Theorem 2.1 establishes that the forecasts produced
by the FLAP method, $\tilde{\bm{y}}_{t+h}$, dominate the base forecasts,
$\hat{\bm{y}}_{t+h}$, in the sense that the difference between their forecast
error variances is always positive semi-definite. Theorem 2.2 establishes that
the trace of the covariance of the forecast errors of $\tilde{\bm{y}}_{t+h}$
is non-increasing with the number of components $p$. The conditions needed to
make the trace strictly decreasing are discussed in Theorem 2.3. In Theorem
2.4 we prove that the projection in Equation 2 achieves the minimum forecast
error variance amongst the class of all projections. Finally, while the
theoretical results imply that components could in principle continue to be
added to improve forecasts, in practice, larger values of $p$ will make
estimation of the plug-in covariance matrix $\bm{W}_{h}$ unreliable. We
explore this issue in a simulation setting and empirically in Section 3 and
Section 4, respectively.

### 2.2 Positive Semi-Definiteness of Error Variance Reduction

We first provide some intermediate results.

######  Lemma 2.1.

Matrix $\bm{M}$ is a projection onto the space where the constraint
$\bm{C}\bm{z}_{t}=\bm{0}$ is satisfied.

Proof in Appendix, page 2.2.

#### Proof of Lemma 2.1

We have

$\displaystyle\bm{M}\bm{M}$
$\displaystyle=\bm{I}_{m+p}-2\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}$
$\displaystyle\mbox{}\hskip
28.45274pt+\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}$
$\displaystyle=\bm{I}_{m+p}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}$
$\displaystyle=\bm{M},$

so $\bm{M}$ is a projection matrix. For any $\bm{z}$ such that
$\bm{M}\bm{z}=\bm{y}$ for some $\bm{y}$, we have

$\bm{C}\bm{y}=\bm{C}\bm{M}\bm{z}=\bm{C}\bm{z}-\bm{C}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{z}=\bm{0}.$

Thus, $\bm{M}$ projects any vector onto the space where the constraint
$\bm{C}\bm{y}=0$ is satisfied.

Based on the properties of projections, we have the following corollaries.

######  Corollary 2.1.

  1. 1.

The projected forecast $\tilde{\bm{z}}_{t+h}$ satisfies the constraint
$\bm{C}\tilde{\bm{z}}_{t+h}=\bm{0}$.

  2. 2.

For $\bm{z}_{t+h}$ that already satisfies the constraint, the projection does
not change its value, i.e., $\bm{M}\bm{z}_{t+h}=\bm{z}_{t+h}$
[Rao1974-JRSSSBSM, Lemma 2.4].

  3. 3.

If the base forecasts are unbiased such that
$\operatorname{E}(\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t})$,
then the projected forecasts are also unbiased, i.e.,
$\operatorname{E}(\tilde{\bm{z}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t})$.

Proof in Appendix, page 2.2.

#### Proof of Corollary 2.1

Items 1 and 2 are trivial application of Lemma 2.1. To prove 3, we have

$\operatorname{E}(\tilde{\bm{z}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{M}\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\bm{M}\operatorname{E}(\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\bm{M}\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{M}\bm{z}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t}).$

The assumption of unbiasedness of the base forecasts is not unreasonable in
practice, and where it does not hold, a bias correction method can be applied.
Note this is not a requirement on model specification. We do not assume the
model producing the base forecast is correctly specified like in the DFM
literature [StoWat2002, e.g.,]. In fact, the power of FLAP manifests when the
models are misspecified, as discussed in Section 3.

######  Lemma 2.2.

The forecast error covariance matrix of the component-augmented projected
$h$-step-ahead forecasts $\tilde{\bm{z}}_{t+h}$ is

$\operatorname{Var}(\bm{z}_{t+h}-\tilde{\bm{z}}_{t+h})=\bm{M}\bm{W}_{h}\bm{M}^{\prime}=\bm{M}\bm{W}_{h},$

and the forecast error covariance matrix of the projected $h$-step-ahead
forecasts $\tilde{\bm{y}}_{t+h}$ is

$\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h})=\bm{J}\bm{M}\bm{W}_{h}\bm{J}^{\prime}.$

Proof in Appendix, page 2.2.

#### Proof of Lemma 2.2

$\operatorname{Var}(\tilde{\bm{z}}_{t+h}-\bm{z}_{t+h})=\operatorname{Var}(\bm{M}\hat{\bm{z}}_{t+h}-\bm{M}\bm{z}_{t+h})\\\
=\bm{M}\operatorname{Var}(\hat{\bm{z}}_{t+h}-\bm{z}_{t+h})\bm{M}^{\prime}\\\
=\bm{M}\bm{W}_{h}\bm{M}^{\prime}.$

If we simplify it further, we have

$\displaystyle\bm{M}\bm{W}_{h}\bm{M}^{\prime}$
$\displaystyle=(\bm{I}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C})\bm{W}_{h}(\bm{I}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C})^{\prime}$
$\displaystyle=\bm{W}_{h}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}$
$\displaystyle\mbox{}\hskip
28.45274pt+\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}$
$\displaystyle=\bm{W}_{h}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}$
$\displaystyle=\bm{M}\bm{W}_{h}.$

To get $\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h})$, we just need
to recognise that it is the first $m\times m$ leading principal submatrix of
$\operatorname{Var}(\tilde{\bm{z}}_{t+h}-\bm{z}_{t+h})$.

Lemma 2.2 is a well-known result in the forecast reconciliation literature
[DiGir2023, e.g.,].

######  Theorem 2.1 (Positive Semi-Definiteness of Error Variance Reduction).

The difference between the forecast error variances of the base and projected
component-augmented forecasts,

$\displaystyle\operatorname{Var}(\bm{z}_{t+h}-\hat{\bm{z}}_{t+h})-\operatorname{Var}(\bm{z}_{t+h}-\tilde{\bm{z}}_{t+h})$
$\displaystyle=\bm{W}_{h}-\bm{M}\bm{W}_{h}$
$\displaystyle=\bm{W}_{h}-(\bm{I}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C})\bm{W}_{h}$
$\displaystyle=\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h},$

is positive semi-definite. The difference between the forecast error variances
of the base and projected forecasts of the original series,

$\operatorname{Var}(\bm{y}_{t+h}-\hat{\bm{y}}_{t+h})-\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h})=\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime},$

is therefore also positive semi-definite.

Proof in Appendix, page 2.2.

#### Proof of Theorem 2.1

Trivially,
$\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}$
and
$\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime}$
are positive semi-definite. Note that
$\operatorname{Var}(\hat{\bm{y}}_{t+h}-\bm{z}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{z}_{t+h})$
is the leading principal submatrix of
$\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}$,
and the leading principal submatrix of a positive semi-definite matrix is
positive semi-definite.

Theorem 2.1 is why FLAP works. Note when
$\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime}$
is semi-definite, each of its diagonal elements, the reduction in forecast
error variance for each series, is non-nagative. It implies that the forecast
error variance can be reduced by simply forecasting the components (the
artificially constructed linear combinations of the original data), and
mapping the forecasts using matrix $\bm{M}$. For the improvement to be zero,
the trace must be zero. This implies that
$\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime}=\bm{O}_{m\times
m}$, as this is a positive semi-definite matrix, something rarely observed in
practice. See Theorem 2.3 for more discussion.

The following example illustrates the mechanism of the reduction in the
forecast error variance.

######  Example 2.1.

Suppose $\bm{z}_{t}$ comprises $m$ original series $\bm{y}_{t}$ and $p$
components $\bm{c}_{t}$. Let $\hat{\bm{z}}_{t+h}$ and $\tilde{\bm{z}}_{t+h}$
be $h$-step-ahead base and projected forecasts of $\bm{z}_{t}$. Assume that
their corresponding forecast errors are uncorrelated with unit variance such
that $\bm{W}_{h}=\bm{I}_{m+p}$. Then

$\displaystyle\operatorname{Var}(\bm{z}_{t+h}-\hat{\bm{z}}_{t+h})-\operatorname{Var}(\tilde{\bm{z}}_{t+h}-\bm{z}_{t+h})$
$\displaystyle=\bm{C}^{\prime}(\bm{C}\bm{C}^{\prime})^{-1}\bm{C}$
$\displaystyle=\begin{bmatrix}-\bm{\Phi}^{\prime}\\\
\bm{I}_{p}\end{bmatrix}(\bm{\Phi}\bm{\Phi}^{\prime}+\bm{I}_{p})^{-1}\begin{bmatrix}-\bm{\Phi}&\bm{I}_{p}\end{bmatrix},$

where

$\bm{C}=\begin{bmatrix}-\bm{\Phi}&\bm{I}_{p}\end{bmatrix}.$

Let $\bm{\Phi}$ consist of $p\leq m$ orthogonal unit vectors, for example,
those obtained from Principal Component Analysis [Jol2002, PCA,]. In this case
$\bm{\Phi}\bm{\Phi}^{\prime}=\bm{I}_{p}$ and

$\operatorname{Var}(\bm{z}_{t+h}-\hat{\bm{z}}_{t+h})-\operatorname{Var}(\bm{z}_{t+h}-\tilde{\bm{z}}_{t+h})=\frac{1}{2}\begin{bmatrix}\bm{\Phi}^{\prime}\bm{\Phi}&-\bm{\Phi}^{\prime}\\\
-\bm{\Phi}&\bm{I}_{p}\end{bmatrix}.$

Focusing on the forecast error variance reduction of the forecasts of the
original series $\bm{y}_{t}$, i.e.,
$\operatorname{tr}(\operatorname{Var}(\bm{y}_{t+h}-\hat{\bm{y}}_{t+h})-\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h}))=\frac{1}{2}\operatorname{tr}(\bm{\Phi}^{\prime}\bm{\Phi})$.

  * •

When $p<m$, since $\bm{\Phi}^{\prime}\bm{\Phi}$ is idempotent,
$\operatorname{tr}(\bm{\Phi}^{\prime}\bm{\Phi})=\operatorname{rank}(\bm{\Phi}^{\prime}\bm{\Phi})=p$.
Hence, focusing on the original series the reduction in the total forecast
error variance of the FLAP forecasts relative to the base forecasts, is $p/2$.

  * •

When $p=m$, where all principal components are used,
$\bm{\Phi}^{\prime}\bm{\Phi}=\bm{I}_{m}$. This implies a total reduction in
the error variance of $m/2$, and that for each of the $m$ individual series
the error variance is halved.

If we keep increasing the number of components beyond $m$, the result in
Theorem 2.1 still holds, although $\bm{\Phi}$ can no longer contain orthogonal
vectors, and the example here becomes intractable. This is an artificial
example as the forecast error variance $\bm{W}_{h}$ can hardly be an identity
matrix in practice. It is likely that the forecast error of a linear
combination of series will be correlated to the forecast error of forecasting
these series directly. Nonetheless, the aim of the example is to demonstrate
how the forecast error variance can be reduced as a result of the component
forecasts bringing new information about the original series. The forecast
error variance reduction becomes larger as we increase the number of
components $p$. This is not a coincidence but a desirable property of FLAP, as
shown in the next section.

### 2.3 Monotonicity

In the results that follow, we break the base forecast error covariance matrix
into smaller blocks.

$\bm{W}_{h}=\begin{bmatrix}\bm{W}_{y,h}&\bm{W}_{yc,h}\\\
\bm{W}_{yc,h}^{\prime}&\bm{W}_{c,h}\end{bmatrix},$

where $\bm{W}_{y,h}$ is the forecast error covariance matrix of
$\hat{\bm{y}}_{t+h}$, $\bm{W}_{c,h}$ is the forecast error covariance matrix
of $\hat{\bm{c}}_{t+h}$, and $\bm{W}_{yc,h}$ contains error covariances
between elements of $\hat{\bm{y}}_{t+h}$ and $\hat{\bm{c}}_{t+h}$.

######  Theorem 2.2 (Monotonicity).

The forecast error variance reductions for each series, i.e., the diagonal
elements in the matrix

$\operatorname{Var}(\bm{y}_{t+h}-\hat{\bm{y}}_{t+h})-\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h})=\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime}$

are non-decreasing as $p$ increases. In particular, the sum of forecast error
variance reductions

${\operatorname{tr}(\operatorname{Var}(\bm{y}_{t+h}-\hat{\bm{y}}_{t+h})-\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h}))=\operatorname{tr}(\bm{J}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\bm{W}_{h}\bm{J}^{\prime})}$
(4)

is non-decreasing as $p$ increases.

Proof in Appendix, page 2.3.

#### Proof of Theorem 2.2

Suppose now that we want to include $q$ more components
$\bm{c}_{t}^{*}=\bm{\Phi}^{*}\bm{y}_{t}$ in the projection. We define
$\bm{z}_{t}^{*}=\begin{bmatrix}\bm{z}_{t}\\\ \bm{c}_{t}^{*}\end{bmatrix}$, the
constraint matrix

${\bm{C}^{*}=\begin{bmatrix}\lx@intercol\hfil\bm{C}\hfil\lx@intercol&\underset{p\times
q}{\bm{0}}\\\ -\underset{q\times m}{\bm{\Phi}^{*}}&\underset{q\times
p}{\bm{0}}&\bm{I}_{q}\end{bmatrix}=\begin{bmatrix}-\underset{p\times
m}{\bm{\Phi}}&\bm{I}_{p}&\underset{p\times q}{\bm{0}}\\\ -\underset{q\times
m}{\bm{\Phi}^{*}}&\underset{q\times
p}{\bm{0}}&\bm{I}_{q},\end{bmatrix}=\begin{bmatrix}\overline{\bm{C}}\\\
\underline{\bm{C}}\end{bmatrix}}$ (5)

where $\overline{\bm{C}}$ contains the first $p$ rows of $\bm{C}^{*}$ and
$\underline{\bm{C}}$ contains the remaining $q$ rows of $\bm{C}^{*}$, the
forecast error variance matrix

$\operatorname{Var}(\hat{\bm{z}}_{t+h}^{*}-\bm{z}_{t+h}^{*})=\bm{W}_{h}^{*}=\begin{bmatrix}\bm{W}_{h}&\bm{W}_{yc,h}^{*}\\\
\bm{W}_{cy,h}^{*}&\bm{W}_{c,h}^{*}\end{bmatrix}.$

where $\hat{\bm{z}}_{t+h}^{*}$ is the $h$-step-ahead base forecasts of
$\bm{z}_{t}^{*}$:

$\hat{\bm{z}}_{t+h}^{*}=\begin{bmatrix}\hat{\bm{z}}_{t+h}\\\
\hat{\bm{c}}_{t+h}^{*}\end{bmatrix},$

and the corresponding

$\bm{M}^{*}=\bm{I}-\bm{W}_{h}^{*}\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*}.$

Proving Theorem 2.2 requires proving the following two items.

  1. 1.

Including additional components in the mapping without including corresponding
component constraints is equivalent to not including these additional
components at all.

  2. 2.

For a fixed set of components to be included in the mapping, adding
constraints will reduce forecast error variance.

We start by proving the first statement. Consider the case where we include
the additional series $\bm{c}_{t}^{*}$ without using the additional constraint
$\bm{\Phi}^{*}$. Defining $\bm{M}^{+}$ only with $\overline{\bm{C}}$:

${\bm{M}^{+}=\bm{I}_{m+p+q}-\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}},}$
(6)

we have $\tilde{\bm{z}}_{t+h}^{+}=\bm{M}^{+}\hat{\bm{z}}_{t+h}^{*}$. Further,
we obtain

$\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}=\begin{bmatrix}\bm{W}_{h}&\bm{W}_{yc,h}^{*}\\\
\bm{W}_{cy,h}^{*}&\bm{W}_{c,h}^{*}\end{bmatrix}\begin{bmatrix}\bm{C}^{\prime}\\\
\underset{q\times p}{\bm{0}}\\\ \end{bmatrix}\\\
=\begin{bmatrix}\bm{W}_{h}\bm{C}^{\prime}\\\
\bm{W}_{cy,h}^{*}\bm{C}^{\prime}\end{bmatrix}$

and

$\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}=\begin{bmatrix}\bm{C}&\underset{p\times
q}{\bm{0}}\\\ \end{bmatrix}\begin{bmatrix}\bm{W}_{h}\bm{C}^{\prime}\\\
\bm{W}_{cy,h}^{*}\bm{C}^{\prime}\end{bmatrix}\\\
=\bm{C}\bm{W}_{h}\bm{C}^{\prime},$

which gives

$\displaystyle\bm{M}^{+}$
$\displaystyle=\bm{I}_{m+p+q}-\begin{bmatrix}\bm{W}_{h}\bm{C}^{\prime}\\\
\bm{W}_{cy,h}^{*}\bm{C}^{\prime}\end{bmatrix}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\begin{bmatrix}\bm{C}&\underset{p\times
q}{\bm{0}}\\\ \end{bmatrix}$
$\displaystyle=\bm{I}_{m+p+q}-\begin{bmatrix}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}&\bm{0}\\\
\bm{W}_{cy,h}^{*}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}&\bm{0}\end{bmatrix},$

and

$\displaystyle\tilde{\bm{z}}_{t+h}^{+}$
$\displaystyle=\bm{M}^{+}\hat{\bm{z}}_{t+h}^{*}$
$\displaystyle=\left(\bm{I}_{m+p+q}-\begin{bmatrix}\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}&\bm{0}\\\
\bm{W}_{cy,h}^{*}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}&\bm{0}\end{bmatrix}\right)\begin{bmatrix}\hat{\bm{z}}_{t+h}\\\
\hat{\bm{c}}_{t+h}^{*}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}(\bm{I}_{n+p}-\bm{W}_{h}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C})\hat{\bm{z}}_{t+h}\\\
\hat{\bm{c}}_{t+h}^{*}-\bm{W}_{cy,h}^{*}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\hat{\bm{z}}_{t+h}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bm{M}\hat{\bm{z}}_{t+h}\\\
\hat{\bm{c}}_{t+h}^{*}-\bm{W}_{cy,h}^{*}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\hat{\bm{z}}_{t+h}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\tilde{\bm{z}}_{t+h}\\\
\hat{\bm{c}}_{t+h}^{*}-\bm{W}_{cy,h}^{*}\bm{C}^{\prime}(\bm{C}\bm{W}_{h}\bm{C}^{\prime})^{-1}\bm{C}\hat{\bm{z}}_{t+h}\end{bmatrix}.$

If we only consider the forecast performance relevant to $\bm{y}_{t+h}$, and
define
$\bm{J}^{*}=\bm{J}_{m,p+q}=\begin{bmatrix}\bm{I}_{m}&\bm{0}_{m\times(p+q)}\end{bmatrix}$,
we have

$\tilde{\bm{y}}_{t+h}^{+}=\bm{J}^{*}\tilde{\bm{z}}_{t+h}^{+}=\bm{J}\tilde{\bm{z}}_{t+h}=\tilde{\bm{y}}_{t+h}.$

This means adding additional components without imposing the corresponding
constraints will yield the same projected forecasts as if these additional
components are not added, which implies that the forecast error variance stays
the same:

${\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{+}-\bm{y}_{t+h})=\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h})=\bm{J}\bm{M}\bm{W}_{h}\bm{J}^{\prime}.}$
(7)

This finishes the proof of the first statement. Now we move on to proving the
second statement. We have the forecast error variance matrices

$\displaystyle\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{+}-\bm{z}_{t+h}^{*})$
$\displaystyle=\bm{M}^{+}\bm{W}_{h}^{*}=(\bm{I}_{m+p+q}-\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}})\bm{W}_{h}^{*}$
$\displaystyle\text{and}\qquad\qquad\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{*}-\bm{z}_{t+h}^{*})$
$\displaystyle=\bm{M}^{*}\bm{W}_{h}^{*}=(\bm{I}_{m+p+q}-\bm{W}_{h}^{*}\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*})\bm{W}_{h}^{*}.$

Taking the difference, we have

$\displaystyle\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{+}-\bm{z}_{t+h}^{*})-\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{*}-\bm{z}_{t+h}^{*})$
$\displaystyle=(\bm{W}_{h}^{*}\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*}-\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}})\bm{W}_{h}^{*}$
$\displaystyle=\bm{W}_{h}^{*}(\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*}-\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}})\bm{W}_{h}^{*}.$

Using block matrix inversion, we have

$\displaystyle\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*}$
$\displaystyle=\begin{bmatrix}\overline{\bm{C}}^{\prime}&\underline{\bm{C}}^{\prime}\end{bmatrix}\begin{bmatrix}\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}&\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}\\\
\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}&\underline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}\end{bmatrix}^{-1}\begin{bmatrix}\overline{\bm{C}}\\\
\underline{\bm{C}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\overline{\bm{C}}^{\prime}&\underline{\bm{C}}^{\prime}\end{bmatrix}\begin{bmatrix}a&b\\\
c&d\end{bmatrix}\begin{bmatrix}\overline{\bm{C}}\\\
\underline{\bm{C}}\end{bmatrix}$
$\displaystyle=\overline{\bm{C}}^{\prime}a\overline{\bm{C}}+\overline{\bm{C}}^{\prime}b\underline{\bm{C}}+\underline{\bm{C}}^{\prime}c\overline{\bm{C}}+\underline{\bm{C}}^{\prime}d\underline{\bm{C}},$

where

$\displaystyle a$
$\displaystyle=(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}+(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}$
$\displaystyle\mbox{}\hskip
28.45274pt(\underline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}-\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}$
$\displaystyle=(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}+(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1},$
$\displaystyle b$
$\displaystyle=-(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}-\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}$
$\displaystyle=-(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1},$
$\displaystyle c$
$\displaystyle=-(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1},$
$\displaystyle d$
$\displaystyle=(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}.$

Thus,

$\displaystyle\bm{C}^{*\prime}(\bm{C}^{*}\bm{W}_{h}^{*}\bm{C}^{*\prime})^{-1}\bm{C}^{*}$
$\displaystyle=\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}$
$\displaystyle\mbox{}\hskip
28.45274pt+\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}$
$\displaystyle\mbox{}\hskip
28.45274pt-\overline{\bm{C}}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}$
$\displaystyle\mbox{}\hskip
28.45274pt-\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}^{\prime}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}$
$\displaystyle\mbox{}\hskip
28.45274pt+\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}$
$\displaystyle=\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}$
$\displaystyle\mbox{}\hskip
28.45274pt-\overline{\bm{C}}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+}$
$\displaystyle\mbox{}\hskip
28.45274pt+\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+}$
$\displaystyle=\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}+\bm{M}^{+\prime}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+}.$

Therefore,

$\displaystyle\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{+}-\bm{z}_{t+h}^{*})-\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{*}-\bm{z}_{t+h}^{*})$
$\displaystyle=\bm{W}_{h}^{*}(\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}}+\bm{M}^{+\prime}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+}-\overline{\bm{C}}^{\prime}(\overline{\bm{C}}\bm{W}_{h}^{*}\overline{\bm{C}}^{\prime})^{-1}\overline{\bm{C}})\bm{W}_{h}^{*}$
$\displaystyle=\bm{W}_{h}^{*}(\bm{M}^{+\prime}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+})\bm{W}_{h}^{*}$

is positive semi-definite. This concludes the proof of the second statement.
Combining the results above, we have

$\displaystyle\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{*}-\bm{y}_{t+h})$
$\displaystyle=\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{+}-\bm{y}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{*}-\bm{y}_{t+h})$
(8)
$\displaystyle=\bm{J}^{*}\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{+}-\bm{z}_{t+h}^{*})\bm{J}^{*\prime}-\bm{J}^{*}\operatorname{Var}(\tilde{\bm{z}}_{t+h}^{*}-\bm{z}_{t+h}^{*})\bm{J}^{*\prime}$
$\displaystyle=\bm{J}^{*}\bm{W}_{h}^{*}\bm{M}^{+\prime}\underline{\bm{C}}^{\prime}(\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\underline{\bm{C}}^{\prime})^{-1}\underline{\bm{C}}\bm{M}^{+}\bm{W}_{h}^{*}\bm{J}^{*\prime}$

being positive semi-definite. Finally, we have

$(\operatorname{Var}(\hat{\bm{y}}_{t+h}-\bm{y}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{*}-\bm{y}_{t+h}))-(\operatorname{Var}(\hat{\bm{y}}_{t+h}-\bm{y}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h}))\\\
=\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h})-\operatorname{Var}(\tilde{\bm{y}}_{t+h}^{*}-\bm{y}_{t+h})$

being a positive semi-definite matrix where the diagonal terms are non-
negative, whose trace, therefore, is non-negative. This means using a larger
number of components in the mapping achieves lower or equal forecast error
variances. In other words, the reduction in forecast error variance of each
series is non-decreasing as more components are added, giving Theorem 2.2.

Theorem 2.2 is the key result that demonstrates the usefulness of FLAP. It
means that we can keep increasing the number of components to reduce forecast
error variance, even when the number of components exceeds the number of
original series. It requires $\bm{C}$ to be
$\big{[}-\bm{\Phi}~{}~{}~{}\bm{I}_{p}\big{]}$ or
$\big{[}-\bm{\Phi}~{}~{}~{}\bm{L}\big{]}$ where $\bm{L}$ is a lower triangular
matrix. This implies that the components can also be constructed from existing
components, not only from the original series. This has little significance
since a linear combination of components of the original series, is just a
linear combination of the original series. This of course assumes that the
forecast error covariances are known, something we explore in Section 2.5 and
Section 3.

Extending the proof of Theorem 2.2, we can state the condition required for
the reduction sum of forecast error variances to be strictly positive. Denote
$\bm{\phi}_{i}$ as the row vector containing the weights associated with the
$i$th component, so that with $p$ components, the weights matrix is
$\bm{\Phi}=\big{[}\bm{\phi}_{1}^{\prime}~{}~{}~{}\bm{\phi}_{2}^{\prime}~{}~{}~{}\dots~{}~{}~{}\bm{\phi}_{p}^{\prime}\big{]}^{\prime}$.
Let $\bm{W}^{(i-1)}_{\tilde{\bm{y}},h}$ denote the error covariance matrix of
the projected forecasts of the original series based on the first $i-1$
components, $\bm{w}_{c_{1}\hat{\bm{y}},h}$ denote a vector of covariances of
the first component and the base forecasts of the original series, and
$\bm{w}_{c_{i}\tilde{y},h}^{(i-1)}$ denote a vector of covariances of the
projected $i$th component and the projected forecasts of the original series,
based on the first $i-1$ components.

######  Theorem 2.3 (Positive Forecast Error Variance Reduction Condition).

For the first component to achieve a guaranteed reduction of forecast error
variance (for the matrix in Theorem 2.1 to have strictly positive trace),

${\bm{\phi}_{1}\bm{W}_{y,h}\neq\bm{w}_{c_{1}y,h}.}$ (9)

For the $i$th component to have a positive reduction on forecast error
variance of the original series,

${\bm{\phi}_{i}{\bm{W}}^{(i-1)}_{\tilde{y},h}\neq\bm{w}_{c_{i}\tilde{y},h}^{(i-1)}.}$
(10)

Proof in Appendix, page 2.3.

#### Proof of Theorem 2.3

Denote
$\bm{\psi}_{i}=\begin{bmatrix}-\bm{\phi}_{i}&\bm{0}_{1\times(i-1)}&1\end{bmatrix}$
and $\bm{W}_{h}^{(i)}$ to be the base forecast error variance of the original
series and the first $i$ components. Starting with the first component,
Equation 4 becomes

${\operatorname{tr}(\bm{J}_{m,1}\bm{W}_{h}^{(1)}\bm{\psi}_{1}^{\prime}(\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{\psi}_{1}^{\prime})^{-1}\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{J}_{m,1}^{\prime})=(\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{\psi}_{1}^{\prime})^{-1}\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{J}_{m,1}^{\prime}\bm{J}_{m,1}\bm{W}_{h}^{(1)}\bm{\psi}_{1}^{\prime},}$
(11)
$\displaystyle\text{where}\qquad\qquad\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{J}_{m,1}^{\prime}=$
$\displaystyle\begin{bmatrix}-\bm{\phi}_{1}&1\end{bmatrix}\begin{bmatrix}\bm{W}_{z,h}&\bm{w}^{\prime}_{c_{1}z,h}\\\
\bm{w}_{c_{1}z,h}&\bm{W}_{c_{1},h}\end{bmatrix}\begin{bmatrix}\bm{I}_{m}\\\
0\end{bmatrix}$ $\displaystyle=$
$\displaystyle-\bm{\phi}_{1}\bm{W}_{z,h}+\bm{w}_{c_{1}z,h}.$

Equation 11 is obviously non-negative. For it to be larger than $0$, we need
$\bm{\psi}_{1}\bm{W}_{h}^{(1)}\bm{J}_{m,1}^{\prime}\neq 0$, which gives
$\bm{\phi}_{1}\bm{W}_{z,h}\neq\bm{w}_{c_{1}z,h}$.

When it comes to adding the $i$th component on top of the first $i-1$
components, we define

$\overline{\bm{C}}_{i}=\begin{bmatrix}\bm{\psi}_{1}&\bm{0}_{1\times i}\\\
\bm{\psi}_{2}&\bm{0}_{1\times(i-1)}\\\ \vdots&\vdots\\\ \bm{\psi}_{i}&0\\\
\end{bmatrix}$

and

$\bm{M}^{+}_{i}=\bm{I}_{m+i}-\bm{W}_{h}^{(i)}\overline{\bm{C}}_{i-1}^{\prime}(\overline{\bm{C}}_{i-1}\bm{W}_{h}^{(i)}\overline{\bm{C}}_{i-1}^{\prime})^{-1}\overline{\bm{C}}_{i-1}$

analogously to Equation 5 and Equation 6. Following Equation 8, the additional
reduction of forecast error variance when adding the $i$th component becomes

$\bm{J}_{m,i}\bm{W}_{h}^{(i)}\bm{M}^{+\prime}_{i}\bm{\psi}_{i}^{\prime}(\bm{\psi}_{i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{\psi}_{i}^{\prime})^{-1}\bm{\psi}_{i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{J}_{m,i}^{\prime}=(\bm{\psi}_{i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{\psi}_{i}^{\prime})^{-1}\bm{\psi}_{i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{J}_{m,i}^{\prime}\bm{J}_{m,i}\bm{W}_{h}^{(i)}\bm{M}^{+\prime}_{i}\bm{\psi}_{i}^{\prime}.$

Similar to before, we would want
$\bm{\psi}_{i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{J}_{m,i}^{\prime}\neq\bm{0}$.
Note that $\bm{\psi}_{i}$ concerns the first $m$ rows and the last row of
$\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}$, and $\bm{J}_{m,i}^{\prime}$ concerns the
first $m$ columns. Combined with the implication from Equation 7 that the
$m\times m$ leading principal submatrix in equation
$\bm{J}_{m,i}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{J}_{m,i}^{\prime}=\bm{J}_{m,i-1}\bm{M}_{i-1}\bm{W}_{h}^{(i-1)}\bm{J}_{m,i-1}^{\prime}$
is the same, we suppress the straightforward yet tiresome details, and obtain

$\bm{\phi}_{i}\bm{W}_{\tilde{z},h}^{(i-1)}\neq\big{[}\bm{0}_{1\times
m+i-1}~{}~{}~{}1\big{]}\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}\bm{J}_{m,i}^{\prime},$

where
$\bm{W}_{\tilde{z},h}^{(i-1)}=\bm{J}_{m,i-1}\bm{M}_{i-1}\bm{W}_{h}^{(i-1)}\bm{J}_{m,i-1}^{\prime}$
is the projected forecast error variance of the original series using the
first $i-1$ components, and the right hand side of the inequality is simply a
one-row matrix consisting of the first $m$ elements in the last row of
$\bm{M}^{+}_{i}\bm{W}_{h}^{(i)}$, which can be denoted as
$\bm{w}_{\tilde{c}_{i}\tilde{z},h}^{(i-1)}$ and interpreted as the covariance
between the projected forecast of the original series using the first $i-1$
components, and the projected forecast of the $i$th component using the first
$i-1$ components.

The condition of Equation 9 demonstrates that for a new component to be
beneficial, the information introduced by this new component, reflected in the
error covariance, cannot be a linear combination of already existing
information.

Theorem 2.3 can potentially provide insights into the selection of component
weights and forecast models to satisfy the conditions. Heuristically, we want
the components to be introducing new information, and this is more likely when
they are uncorrelated, which motivates out choice of principal components
which are uncorrelated by constructions. Optimal choices of weights we leave
to future research, but note that practically speaking, the conditions in
Theorem 2.3 are either almost always satisfied if the weights are simulated
randomly on a continuous scale, or the loss associated with the rare occasions
where the conditions are not satisfied is negligible compared to the
estimation error imposed by the limited sample size as the number of
components increases, as discussed in Section 3 and Section 4.

### 2.4 Optimality of the projection

Equation 1 can be seen as a solution to the optimisation problem,

$\underset{\check{\bm{z}}_{T+h}}{\arg\min}(\hat{\bm{z}}_{T+h}-\check{\bm{z}}_{T+h})^{\prime}\bm{W}_{h}^{-1}(\hat{\bm{z}}_{T+h}-\check{\bm{z}}_{T+h})\qquad\text{s.t.
}\bm{C}\check{\bm{z}}_{T+h}=0.$

If we consider the transformed space where all the vectors are first
transformed via pre-multiplying by $\bm{W}_{h}^{-1/2}$, where
$\bm{W}_{h}^{-1}=(\bm{W}_{h}^{-1/2})^{\prime}\bm{W}_{h}^{-1/2}$, then this
optimisation problem can be interpreted as finding the set of forecasts that
are closest to the base forecasts on the transformed space, while satisfying
the linear constraints imposed by the components.

An alternative way of characterising the same constraints is

$\bm{\Phi}\check{\bm{y}}_{T+h}=\check{\bm{c}}_{T+h},$

where $\check{\bm{c}}_{T+h}$ is the vector of the last $p$ elements of
$\check{\bm{z}}_{T+h}$, corresponding to the forecast of the components as
part of the solution. This equivalence is discussed in [WicEtAl2019], where
the authors find the solution by minimising the sum of forecast error variance
of all series (See [AndNar2024-F] for an alternative proof). The result is

${\tilde{\bm{z}}_{t+h}=\bm{S}\bm{G}\hat{\bm{z}}_{t+h},}$ (12)

where $\bm{S}=\begin{bmatrix}\bm{I}_{m}\\\ \bm{\Phi}\end{bmatrix}$ contains
the constraints, so that $\bm{z}_{t}=\bm{S}\bm{y}_{t}$, and

${\bm{G}=(\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S})^{-1}\bm{S}^{\prime}\bm{W}_{h}^{-1}.}$
(13)

In Equation 12, $\bm{G}\hat{\bm{z}}_{t+h}$ can be viewed as mapping all of the
series to a selected few. In the forecast reconciliation context, this is a
mapping of all series to the “bottom level”. In our multivariate forecasting
context, this is a mapping of all series including the components, to the
space of the original series. This leads to the solution

${\tilde{\bm{y}}_{t+h}=\bm{G}\hat{\bm{z}}_{t+h},}$ (14)

as equivalent to Equation 3. Recognising that Equation 12 is equivalent to
Equation 1, it is the solution that minimises the sum of forecast error
variances of the original series and all the components. We go further in
Theorem 2.4, and show that Equation 12 is also the solution to minimise each
individual forecast error variance of the original series, and their sum. This
can be viewed as a special case of Theorem 3.3 in [PanEtAl2021] applied in a
forecast projection context, by taking
$\bm{W}=\bm{S}(\bm{S}^{\prime}\bm{S})^{-1}(\bm{S}^{\prime}\bm{S})^{-1}\bm{S}^{\prime}$
in [PanEtAl2021, Thm. 3.3], as illustrated by [AndNar2024-F]. The earliest
work we can find that noted this interpretation in a non-forecasting context
is [Lue1969, p. 85]. We establish a few basic results leading to the
optimality of this solution first, also to check that Lemma 2.1, Corollary 2.1
and Lemma 2.2 hold under this alternative representation.

######  Lemma 2.3.

The matrix $\bm{S}\bm{G}$ is a projection onto the space where the constraint
$\bm{C}\bm{z}_{t}=\bm{0}$ is satisfied, provided that $\bm{G}\bm{S}=\bm{I}$.

Proof in Appendix, page 2.4.

#### Proof of Lemma 2.3

If $\bm{G}\bm{S}=\bm{I}$, $\bm{S}\bm{G}$ is a projection matrix:
$\bm{S}\bm{G}\bm{S}\bm{G}=\bm{S}\bm{G}$.

For any $\bm{z}$ such that $\bm{S}\bm{G}\bm{z}=\bm{y}$ for some $\bm{y}$, we
have $\bm{C}\bm{y}=\bm{C}\bm{S}\bm{G}\bm{z}=\bm{0}$ because
$\bm{C}\bm{S}=\big{[}-\bm{\Phi}~{}~{}~{}\bm{I}\big{]}\big{[}\bm{I}~{}~{}~{}\bm{\Phi}^{\prime}\big{]}^{\prime}=\bm{0}$.
Similarly to $\bm{M}$, $\bm{S}\bm{G}$ projects a vector to the same space
where $\bm{C}$ is satisfied.

######  Corollary 2.2.

Provided that $\bm{G}\bm{S}=\bm{I}$, the following results hold.

  1. 1.

The projected forecast in Equation 14 satisfies the constraint
$\bm{C}\tilde{\bm{z}}_{t+h}=\bm{C}\bm{S}\tilde{\bm{y}}_{t+h}=\bm{0}$.

  2. 2.

For $\bm{z}_{t+h}$ that already satisfies the constraint, the mapping does not
change its value, i.e., $\bm{G}\bm{z}_{t+h}=\bm{y}_{t+h}$.

  3. 3.

If the base forecasts are unbiased such that
$\operatorname{E}(\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t})$,
then the projected forecasts in Equation 14 are also unbiased:
$\operatorname{E}(\tilde{\bm{y}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{y}_{t+h}|\mathcal{I}_{t})$.

Proof in Appendix, page 2.4.

#### Proof of Corollary 2.2

Item 1 is an direct application of Lemma 2.3. From Lemma 2.3 and Lemma 2.4 in
[Rao1974-JRSSSBSM], we have

$\bm{S}\bm{G}\bm{z}_{t+h}=\bm{z}_{t+h}=\bm{S}\bm{y}_{t+h}.$

Left multiplying by $\bm{G}$ on both sides, we have
$\bm{G}\bm{z}_{t+h}=\bm{y}_{t+h}$ and item 2 is proven. To prove Item 3, we
have

$\operatorname{E}(\tilde{\bm{y}}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{G}\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\bm{G}\operatorname{E}(\hat{\bm{z}}_{t+h}|\mathcal{I}_{t})=\bm{G}\operatorname{E}(\bm{z}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{G}\bm{z}_{t+h}|\mathcal{I}_{t})=\operatorname{E}(\bm{y}_{t+h}|\mathcal{I}_{t}).$

######  Lemma 2.4.

The covariance matrix of the projected forecasts from Equation 14 is given by

$\operatorname{Var}(\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h})=\bm{G}\bm{W}_{h}\bm{G}^{\prime}.$

Proof in Appendix, page 2.4.

#### Proof of Lemma 2.4

Let the base and projected forecast errors be given as

$\displaystyle\hat{\bm{e}}_{y,t+h}$
$\displaystyle=\bm{y}_{t+h}-\hat{\bm{y}}_{t+h},$
$\displaystyle\hat{\bm{e}}_{z,t+h}$
$\displaystyle=\bm{z}_{t+h}-\hat{\bm{z}}_{t+h},$
$\displaystyle\tilde{\bm{e}}_{y,t+h}$
$\displaystyle=\bm{y}_{t+h}-\tilde{\bm{y}}_{t+h},$
$\displaystyle\text{and}\qquad\tilde{\bm{e}}_{z,t+h}$
$\displaystyle=\bm{z}_{t+h}-\tilde{\bm{z}}_{t+h}=\bm{S}\bm{y}_{t+h}-\bm{S}\tilde{\bm{y}}_{t+h}=\bm{S}\tilde{\bm{e}}_{y,t+h}.$
$\displaystyle\text{Then we have}\qquad\tilde{\bm{e}}_{z,t+h}$
$\displaystyle=\hat{\bm{e}}_{z,t+h}+\hat{\bm{z}}_{t+h}-\tilde{\bm{z}}_{t+h}$
$\displaystyle=\hat{\bm{e}}_{z,t+h}+\hat{\bm{z}}_{t+h}-\bm{S}\bm{G}\hat{\bm{z}}_{t+h}$
$\displaystyle=\hat{\bm{e}}_{z,t+h}+(\bm{I}-\bm{S}\bm{G})(\bm{z}_{t+h}-\hat{\bm{e}}_{z,t+h})$
and
$\displaystyle=\bm{S}\bm{G}\hat{\bm{e}}_{z,t+h}+(\bm{I}-\bm{S}\bm{G})\bm{S}\bm{y}_{t+h}$
$\displaystyle\bm{S}\tilde{\bm{e}}_{y,t+h}$
$\displaystyle=\bm{S}\bm{G}\hat{\bm{e}}_{z,t+h},$

where the last line comes from $\bm{G}\bm{S}=\bm{I}$. Left multiplying by
$\bm{G}$ on both sides, we have

$\bm{G}\bm{S}\tilde{\bm{e}}_{y,t+h}=\bm{G}\bm{S}\bm{G}\hat{\bm{e}}_{z,t+h}\qquad\text{and}\qquad\tilde{\bm{e}}_{y,t+h}=\bm{G}\hat{\bm{e}}_{z,t+h},$

and therefore

$\operatorname{Var}(\tilde{\bm{y}}_{t+h}-\bm{y}_{t+h})=\operatorname{Var}(\tilde{\bm{e}}_{y,t+h})=\operatorname{Var}(\bm{G}\hat{\bm{e}}_{z,t+h})=\bm{G}\operatorname{Var}(\hat{\bm{e}}_{z,t+h})\bm{G}^{\prime}=\bm{G}\bm{W}_{h}\bm{G}^{\prime}.$

We are now ready to present the following theorem.

######  Theorem 2.4 (Minimum Variance Unbiased Projected Forecast).

The solution to

${\underset{\bm{G}}{\arg\min}\
\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})\qquad\text{s.t.
}\bm{G}\bm{S}=\bm{I}}$ (15)

is Equation 13. This problem can be effectively split into independent
subproblems such that
$\bm{G}=\big{[}\bm{g}_{1}~{}~{}\bm{g}_{2}~{}~{}\dots~{}~{}\bm{g}_{m}\big{]}^{\prime}$,
where $\bm{g}_{i}$ is the solution to the subproblem of the $i$th series

${\underset{\bm{g}_{i}}{\arg\min}\
\bm{g}_{i}^{\prime}\bm{W}_{h}\bm{g}_{i}\qquad\text{s.t.
}\bm{g}_{i}^{\prime}\bm{s}_{j}=\delta_{ij},\quad j=1,2,.\ldots,m,}$ (16)

where $\bm{s}_{j}$ is the $j$th column of $\bm{S}$, and $\delta_{ij}$ is the
Kronecker delta function taking value 1 if $i=j$ and 0 otherwise.

Proof in Appendix, page 2.4.

#### Proof of Theorem 2.4

This can be proved in a few different ways. We adopt the approach of
[AndNar2024-F] to obtain the solution to Equation 15, but the procedure from
[Lue1969, p. 85] can also be used, where the problem is divided to Equation 16
and reconstructed to find the solution to Equation 15.

There exists a Lagrange multiplier $\bm{\Lambda}$ such that

$L(\bm{G})=\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})+\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{I}-\bm{G}\bm{S}))$

is stationary at an extremum $\bm{G}$ [Lue1969, p. 243, Theorem 1]. We find
the numerator of the Gateaux differential [Lue1969, p. 171]

$\lim_{\alpha\to 0}\frac{L(\bm{G}+\alpha\bm{H})-L(\bm{G})}{\alpha}$

to be

$\displaystyle L(\bm{G}+\alpha\bm{H})-L(\bm{G})=$
$\displaystyle\operatorname{tr}((\bm{G}+\alpha\bm{H})\bm{W}_{h}(\bm{G}+\alpha\bm{H})^{\prime})+\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{I}-(\bm{G}+\alpha\bm{H})\bm{S}))-\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})-\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{I}-\bm{G}\bm{S}))$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})+\operatorname{tr}(\alpha^{2}\bm{H}\bm{W}_{h}\bm{H}^{\prime})+\operatorname{tr}(\alpha\bm{G}\bm{W}_{h}\bm{H}^{\prime})+\operatorname{tr}(\alpha\bm{H}\bm{W}_{h}\bm{G}^{\prime})+\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{I}-\bm{G}\bm{S}))-\operatorname{tr}(\alpha\bm{\Lambda}^{\prime}\bm{H}\bm{S})$
$\displaystyle-\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})-\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{I}-\bm{G}\bm{S}))$
$\displaystyle=$
$\displaystyle\alpha^{2}\operatorname{tr}(\bm{H}\bm{W}_{h}\bm{H}^{\prime})+\alpha\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{H}^{\prime})+\alpha\operatorname{tr}(\bm{H}\bm{W}_{h}\bm{G}^{\prime})-\alpha\operatorname{tr}(\bm{\Lambda}^{\prime}\bm{H}\bm{S}).$

Thus, the Gateaux differential becomes

$\displaystyle\lim_{\alpha\to
0}\frac{L(\bm{G}+\alpha\bm{H})-L(\bm{G})}{\alpha}=$
$\displaystyle\lim_{\alpha\to
0}\alpha\operatorname{tr}(\bm{H}\bm{W}_{h}\bm{H}^{\prime})+\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{H}^{\prime})+\operatorname{tr}(\bm{H}\bm{W}_{h}\bm{G}^{\prime})-\operatorname{tr}(\bm{\Lambda}^{\prime}\bm{H}\bm{S})$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{H}^{\prime})+\operatorname{tr}(\bm{H}\bm{W}_{h}\bm{G}^{\prime})-\operatorname{tr}(\bm{\Lambda}^{\prime}(\bm{H}\bm{S}))$
$\displaystyle=$
$\displaystyle\operatorname{tr}(2\bm{H}\bm{W}_{h}\bm{G}^{\prime}-\bm{\Lambda}^{\prime}\bm{H}\bm{S})$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\bm{H}(2\bm{W}_{h}\bm{G}^{\prime}-\bm{S}\bm{\Lambda}^{\prime})),$

which we set to zero with a value of $\bm{G}^{*}$

$\displaystyle\operatorname{tr}(\bm{H}(2\bm{W}_{h}\bm{G}^{*\prime}-\bm{S}\bm{\Lambda}^{\prime}))$
$\displaystyle=0$ $\displaystyle 2\bm{W}_{h}\bm{G}^{*}$
$\displaystyle=\bm{S}\bm{\Lambda}^{\prime}$ $\displaystyle\bm{G}^{*\prime}$
$\displaystyle=\frac{1}{2}\bm{W}_{h}^{-1}\bm{S}\bm{\Lambda}^{\prime}.$

Multiplying $\bm{S}^{\prime}$ to the left of both sides. we have

$\bm{S}^{\prime}\bm{G}^{*\prime}=\bm{I}=\frac{1}{2}\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S}\bm{\Lambda}^{\prime}\qquad\text{and}\qquad\bm{\Lambda}^{\prime}=2(\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S})^{-1}$

because $\bm{G}^{*}\bm{S}=\bm{I}$. Putting it back in, we have

$\bm{G}^{*\prime}=\bm{W}_{h}^{-1}\bm{S}(\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S})^{-1}\qquad\text{and}\qquad\bm{G}^{*}=(\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S})^{-1}\bm{S}^{\prime}\bm{W}_{h}^{-1}.$

To see how Equation 15 can be split into separate problems, recognise

$\operatorname{tr}(\bm{G}\bm{W}_{h}\bm{G}^{\prime})=\sum^{m}_{i=1}\bm{g}_{i}^{\prime}\bm{W}_{h}\bm{g}_{i},$

where $\bm{W}_{h}$ is a positive definite variance covariance matrix, which
makes each $\bm{g}_{i}^{\prime}\bm{W}_{h}\bm{g}_{i}$ positive. Additionally,
the element in the $i$th row and $j$th column of $\bm{G}\bm{S}$ is
$\bm{g}_{i}^{\prime}\bm{s}_{j}$, and the element in the $i$th row and $j$th
column of an identity matrix is $\delta_{ij}$. Therefore, the problem in
Equation 15 is $m$ separate problems.

In other words, the forecast projection method gives optimal projected
forecast for a given set of components, in the sense that the unbiased
forecast of each series has minimum variance.

### 2.5 Estimation of $\bm{W}_{h}$

In practice, the base forecast error variance $\bm{W}_{h}$ is unknown and
needs to be estimated. Denote
$\hat{\bm{e}}_{t,h}=\bm{z}_{t}-\hat{\bm{z}}_{t|t-h}$ as the $h$-step-ahead
base forecast in-sample residual. The conventional forecast error variance
matrix estimator

$\widehat{\bm{W}}_{h}=\frac{1}{T-h-1}\sum^{T}_{t=h+1}\hat{\bm{e}}_{t,h}\hat{\bm{e}}_{t,h}^{\prime},$

albeit unbiased, is not considered a good approximation to the true forecast
error variance in a finite sample when $(m+p)\approx T-h$. It is even singular
when $(m+p)>T-h$, which makes the quantities discussed in the previous
sections impossible to calculate. For this reason, while other shrinkage
estimators can be considered, we adopt the covariance shrinkage method of
[SchStr2005] and the variance shrinkage method of [OpgStr2007]. Then, the
estimated forecast error variance matrix is guaranteed to be positive definite
with no numerical problems when computing their inverse. This estimator is
denoted as $\widehat{\bm{W}}_{h}^{shr}=(\hat{w}_{ij,h}^{shr})_{1\leq i,j\leq
m+p}$ with the element in row $i$ and column $j$

$\hat{w}_{ij,h}^{shr}=\hat{r}_{ij,h}^{shr}\sqrt{\hat{v}_{i,h}\hat{v}_{j,h}},$

where $\hat{r}_{ij,h}^{shr}=(1-\hat{\lambda}_{cor})\hat{r}_{ij,h}$ and
$\hat{v}_{i,h}=\hat{\lambda}_{var}\hat{w}_{h,median}+(1-\hat{\lambda}_{var})\hat{w}_{i,h}$,
with $\hat{\lambda}_{cor}$ being the shrinkage intensity parameter for the
correlation

$\hat{\lambda}_{cor}=\min\left(1,\frac{\sum_{i\neq
j}\widehat{\operatorname{var}}(\hat{r}_{ij,h})}{\sum_{i\neq
j}\hat{r}_{ij,h}^{2}}\right),$

and $\hat{\lambda}_{var}$ being the shrinkage intensity parameter for the
variance

$\hat{\lambda}_{var}=\min\left(1,\frac{\sum_{i=1}^{m+p}\widehat{\operatorname{var}}(\hat{w}_{i,h})}{\sum_{i=1}^{m+p}(\hat{w}_{i,h}-\hat{w}_{h,median})^{2}}\right),$

$\hat{r}_{ij,h}$ the sample correlation of the $h$-step-ahead forecast error
between the $i$th and the $j$th series (component) in $\bm{z}_{t}$,
$\hat{w}_{i,h}$ the $h$-step-ahead sample base forecast error variance
associated with the $i$th series (the $i$th diagonal element of
$\widehat{\bm{W}}_{h}$), and $\hat{w}_{h,median}$ the median of the $h$-step-
ahead sample forecast error variance of the series and components (the median
of the diagonal elements of $\widehat{\bm{W}}_{h}$). The estimation of
$\widehat{\bm{W}}_{h}^{shr}$ in the following sections are implemented using
the package corpcor [corpcor] in R [R].

Estimating $\widehat{\bm{W}}_{h}^{shr}$ for each forecast horizon $h$ is
desirable but computationally intensive. It involves the calculation of multi-
step-ahead in-sample residuals of the forecast models, which is especially
challenging for iterative forecasts. Because of this, in practice it is not
unreasonable to assume the $h$-step forecast error variance is proportional to
the $1$-step forecast error variance by a constant $\eta_{h}$, as do
[WicEtAl2019]:

$\widehat{\bm{W}}_{h}^{shr}=\eta_{h}\widehat{\bm{W}}_{1}^{shr}.$

Under this assumption, when $\widehat{\bm{W}}_{h}^{shr}$ is used in Equation
2, the proportionality constant $\eta_{h}$ cancels out regardless of the value
of $h$. We can effectively use only the one-step forecast error variance in
forecast projection, if we only need point forecasts. We calculate
$\widehat{\bm{W}}_{h}^{shr}$ for each $h$ for the simulation example in
Section 3, but assume this proportionality for the application in Section 4.

## 3 Simulation

In this section, we illustrate the performance of FLAP in a simulation
setting. We generate time series of length $T=400$ from a $m=70$ variable
VAR($3$) data generating process (DGP). The coefficients for the VAR DGP are
estimated from the first $70$ series in the Australian tourism data set used
in Section 4.1. The innovations are simulated from a multivariate normal
distribution with an identity covariance matrix. The estimation and simulation
is performed using the tsDyn package [tsDyn].

For each simulated sample we generate $h=1$ to $12$-step-ahead base forecasts
from benchmark models. We implement FLAP with a varying number of components
and component construction methods. This process is repeated $220$ times. The
improvements of FLAP forecasts over base forecasts is assessed, and the
statistical significance of these improvements is evaluated.

### 3.1 Generating base and FLAP forecasts

We generate base forecasts from two benchmarks. The first benchmark is the
univariate ARIMA model. For each series, we generate base forecasts from an
ARIMA model using the auto.arima() function from the forecast package
[forecast] with the default settings.

The second benchmark is the dynamic factor model (DFM). Following
[StoWat2002a], base forecasts

$\hat{y}_{T+h}=\hat{\alpha}_{h}+\sum^{n}_{j=1}\hat{\bm{\beta}}_{hj}^{\prime}\hat{\bm{F}}_{T-j+1}+\sum^{s}_{j=1}\hat{\gamma}_{hj}y_{T-j+1},$

where $\hat{\bm{F}}_{t}$ is the vector of $k$ estimated factors, and
$\hat{y}_{t}$ is the target series to forecast. The factors are estimated
using PCA on demeaned and scaled data. The optimal model is selected for each
series based on the Bayesian information criterion (BIC) from models fitted
using different combinations of meta-parameters in their corresponding range:
$1\leq k\leq 6$, $1\leq n\leq 3$ and $1\leq s\leq 3$. We note that the DFM
produces direct forecasts in the sense that a different model is fitted for
each forecast horizon $h$, in contrast to indirect or iterative forecasts
generated by the ARIMA models.

We use several sets of weights to construct linear components for the FLAP
method. The types of components are listed in Table LABEL:tbl-comp. Principal
components from PCA are established using the prcomp() function in the stats
package [R]. We generate random components using weights simulated from a
standard normal (Norm) distribution and a uniform (Unif) distribution with
range $(-1,1)$. We normalise the weights of all randomly generated components
into unit vectors to ensure numerically stable computation, with
$\bm{\phi}_{i}/\sqrt{\sum_{j}(\phi_{ij}^{2})}$ where $\phi_{ij}$ is the $j$th
value in the weight vector of the $i$th component. The total number of
components $p$ is selected to be either $m=70$ or $300$.

Table 1: Component construction methods for FLAP |
---|---
Component |  Description
PCA |  Principal components from PCA.
Norm |  The weights of all components are simulated from a standard normal distribution.
Unif |  The weights of all components are simulated from a uniform distribution
PCA+Norm |  $m$ principal components from PCA are complemented with random components whose weights are simulated from a standard normal distribution.
PCA+Unif |  $m$ principal components from PCA are complemented with random components whose weights are simulated from a uniform distribution.
Ortho |  A random orthonormal matrix is generated using package pracma [pracma] as the weight matrix.
Ortho+Norm |  A random orthonormal matrix is generated using package pracma, forming the weights of the first $m$ components. Weights of the additional components are simulated from a standard normal distribution.

### 3.2 Forecast evaluation

We employ the Friedman test [Fri1937, Fri1939-JASA] along with post-hoc
Nemenyi tests [Nem1963, HolEtAl2013] using the tsutils package [tsutils] to
compare forecast performance between different methods. The analysis involves
the use of Multiple Comparisons with the Best (MCB) plot introduced by
[KonEtAl2005] to visualise the comparison. The MSE of each series over
different samples is calculated, and the MSEs of all the series are treated as
observations in the Nemenyi test. The average ranks are plotted in Figure 1
for forecast horizons $1$, $6$ and $12$. The methods using FLAP are labelled
“Model – Component Weights – Number of Components”. The benchmarks are
labelled “Model – Benchmark”. These points are marked with triangles. The
shaded region is the confidence interval of the best-performing method.
Methods outside the shaded region are significantly worse than the best model.

In Figure 2, we plot the out-of-sample MSE values as a function of the number
of components $p$ included in the FLAP method. Here we also include the
performance of the DGP VAR model (VAR – DGP), and an estimated VAR model with
the correct specification (VAR – Est.), as well as their corresponding FLAP
generated forecasts. We have not presented the FLAP forecasts using components
generated from a uniform distribution or random orthonormal matrices, as these
are visually identical to the methods with random weights generated from a
standard normal distribution. The vertical black line indicates when $p=m=70$,
the number of original series.

Figure 1: Average ranks of $1$-, $6$\- and $12$-step-ahead MSE of different
model and component specifications in the simulation. The methods using FLAP
are labelled as “Model – Component Weights – Number of Components”. The two
benchmark models generating base forecasts are labelled “Model – Benchmark”.
Their MSEs are marked with triangles. The shaded region is the confidence
interval of the best performing method. Methods outside the shaded region are
significantly worse than the best method. Figure 2: Out-of-sample MSE for base
and FLAP forecasts as the number of components $p$ increases, for forecast
horizons $1$, $6$ and $12$. “VAR – DGP” indicates the performance of the true
data generating VAR model. “VAR – Est.” indicates the performance of the VAR
model with the same structure as the true model but with estimated parameters.
The solid horizontal lines show the MSE for the base forecasts (with no
projection) while the dashed lines show the MSEs of the FLAP forecasts. The
vertical black line indicates the location of $p=m=70$, the number of series.

#### Evaluating base forecasts

Having simulated data from a VAR DGP, we expect the DFM to capture the
correlations between series, something not possible with univariate ARIMA
models. Hence, we expect the DFM benchmark to perform better than the ARIMA
benchmark. This is indeed the case. Figure 1 shows that the base DFM forecasts
are significantly better than base ARIMA forecasts, with the exception for
$h=1$, where the difference in rank is statistically insignificant but only
marginally.

This is also observed in Figure 2. The solid horizontal line representing the
MSE of the base DFM forecasts is always much lower than the horizontal solid
line of the base ARIMA forecasts.

#### Evaluating FLAP forecasts

The most important observation from Figure 1 are the statistically significant
improvements of the FLAP forecasts over their corresponding benchmark base
forecasts. Note that we are referring here to forecasts corresponding to the
same model, i.e., FLAP ARIMA forecasts, compared to base ARIMA forecasts and
FLAP DFM forecasts, compared to base DFM forecasts. The average ranks of the
FLAP forecasts are better than the corresponding base forecasts for all
forecast horizons, and the differences are all statistically significant. The
only exception is for ARIMA base forecasts with FLAP using only the $p=m=70$
principal components, where the difference in rank is statistically
insignificant.

Hence, the number of components seems to be important. The best-performing
models all comprise the maximum $p=300$ components. For $h=1$-step-ahead
forecasts, the differences between the methods with $p=300$ components are not
significantly significant, regardless of the model that generated the base
forecasts or the method used to construct the components.

Indeed, Figure 2 shows that the MSEs for the FLAP ARIMA and DFM forecasts,
represented by the dashed lines, decrease monotonically relative to the base
forecasts, as the number of components increases. This supports Theorem 2.2,
demonstrating that increasing the number of components in FLAP reduces the
forecast error variance. Of course this is only feasible in this ideal setting
where each time series has $400$ observations and we use no more than $300$
components in FLAP. The relatively large number of observations coupled with
the relatively simple VAR DGP, has potentially eased the challenge of
estimation. This continuous reduction in forecast error variance is not always
achievable with real data, as demonstrated in the empirical applications in
Section 4, especially with FRED-MD dataset in Section 4.2.

Worth noting from Figure 2 is the performance of the FLAP ARIMA forecasts
relative to DFM forecasts as the number of components $p$ increases. The MSE
of the FLAP ARIMA forecasts becomes lower than the DFM base forecasts for
$p<70$ for $h=1$, for $70<p<100$ for $h=6$ and for a larger $p$ for $h=12$.
Hence, FLAP is able to enhance univariate ARIMA forecasts with shared
information between series by capturing and projecting these cross-
correlations with the components.

Interestingly, for $h=1$, the MSEs of FLAP ARIMA and FLAP DFM forecasts
converge to the same value as $p$ reaches $300$, no matter how the components
are constructed. For the longer forecast horizons, the MSEs again converge;
they do not reach the same value for $p=300$, although they are closer for
$h=6$ than $h=12$. The reason for this is not the contribution of FLAP, but
the starting point of the two sets of base forecasts. As the forecast horizon
increases from $h=1$ to $6$ and $12$ the ARIMA base forecasts MSE increases.
In contrast the MSEs for DFM base forecasts remains at a relatively constant
level. Thus, it seems that the direct forecasts from the DFM model,
particularly fitting a different model for each forecast horizon $h$, is
advantageous for the longer forecast horizons, in this setting of forecasting
data generated from a simple VAR DGP. This is in contrast to the iterative
nature of ARIMA generated forecasts.

Note that the same forecasts of the components, generated from univariate
ARIMA models, are used for both the FLAP ARIMA and the FLAP DFM methods. It is
the increase in the number of the components, that leads to the dominant
performance of the FLAP forecasts. Hence, there is valuable information in the
series that is not captured by either the ARIMA model or the DFM, but is
captured by the components. Once again, this emphasises the importance of the
components and FLAP. In this extreme case, a simple model after FLAP, performs
as well as a more complicated model, while the forecast model itself is not as
important as FLAP.

### 3.3 Forming components for FLAP

The construction of components is obviously a primary element in the proposed
FLAP method. However, the simulation results show that the method of
generating the components may not be as important as one may have expected.

The most important feature in Figure 1, as discussed in the previous section,
is the number of components used in this simulation setting. The results
clearly show that the FLAP forecasts using the maximum considered $p=300$
components outrank the FLAP forecasts using $p=m=70$ components, no matter how
the components are generated. The only exception seems to be for $h=12$ and
for the DFM base forecasts when using all $p=m=70$ PCA components or
complementing these with randomly generated components from normal weights.
Our conjecture, as also stated above, is that the key difference here lies in
the quality of the direct DFM base forecasts for the longer forecast horizons,
rather than in the component generating mechanism itself. Furthermore, between
the methods using the maximum considered $p=300$, there seems to be no
evidence that the component generating mechanism makes any difference to the
ranks, and especially there seems to be no difference between the distribution
used to generate the random components.

Hence, in Figure 2, we present only the results for principal components and
randomly generated components from a standard normal distribution. The results
from using components generated from a uniform distribution or random
orthonormal matrices are omitted, as these are visually identical to the
latter.

Focusing on the ARIMA base forecasts it seems that for a relatively small $p$,
the MSE decreases at a faster rate when principal components are used compared
to random components. However, FLAP with random components seems to achieve a
lower MSE than PCA as $p$ increases. The gap between the two diminishes for
large $p$. We conjecture that the initial favourable performance of PCA stems
from the construction of the orthogonal components aimed at capturing the
maximum variance in the data. These components are then ranked from largest to
smallest variance, which likely contributes to the method’s effectiveness for
smaller $p$ compared to random components. We comment more on PCA when
analysing the results in Section 4.

The results show that the performance of FLAP using orthonormal weights is
almost identical to FLAP using weights simulated from a normal distribution.
This leads us to infer that it is not the orthogonality of the principal
components that leads to the strong performance of FLAP using PCA. It may also
suggest using simple random weights if one is prepared to include a relatively
large number of components in FLAP; and using PCA when the intention is to use
a small number of components. However, as we will see in Section 4, this is
not the case with real data. In this case, PCA seems to be the preferred
approach even when the number of components is large.

### 3.4 Sources of uncertainty

The MSEs of the true VAR data generating process, VAR-DGP, and the correctly
specified but estimated VAR, VAR-Est, are plotted at the bottom of each panel
in Figure 2. As expected, this confirms that they are the best performing
models. As the forecast horizon increases so does the MSE difference between
these, with estimation error for the VAR-Est accumulating as forecasts are
iteratively generated for the longer forecast horizons.

For both these models, FLAP cannot improve on the base forecasts. More
importantly, in the case of VAR-Est, FLAP does not seem to be able to reduce
the estimation error or the parameter uncertainty. On the other hand, FLAP
shows significant improvements over base forecasts generated from misspecified
models, such as the two benchmarks. This implies that the uncertainty which
FLAP can account for and reduce, is mainly due to model misspecification. It
seems that FLAP operates similarly to bagging, as bagging also reduces
variance by controlling model uncertainty [PetEtAl2018]. The uncertainty in
how the components are constructed, unique in FLAP problems, is discussed in
Section 3.3 and awaits future research.

## 4 Empirical applications

In this section, we apply the FLAP method in an empirical setting using two
diverse datasets. The analysis confirms many of the results and conclusions
drawn from the simulation setting, but it also identifies a few key
differences.

### 4.1 Australian domestic tourism

The Australian Tourism Data Set compiled from the National Visitor Survey by
Tourism Research Australia contains the total number of nights spent by
Australians away from home. We refer to these as visitor nights. The monthly
visitor nights are recorded for $m=77$ regions around Australia, covering the
period January 1998 to December 2019. To evaluate the forecast performance, we
conduct am expanding window time series cross-validation [fpp3]. More
specifically, the first $T=84$ observations are used as the first training set
and the following $12$ months are used as the test set for evaluation. We
repeat the evaluation for the rest of the data, by expanding each training
sample by one observation at a time. This generates $169$ forecasts for
evaluation, for each of the forecast horizons, $h=1$ to $12$. Base forecasts
for each series and for the components are generated using univariate ETS
models selected and fitted using the ets() function in the forecast package
[forecast, HynKha2008]. See [fpp2, Chapter 7] for more details.

Figure 3: Average ranks of $1$-, $6$\- and $12$-step-ahead cross-validation
MSE of different model and component specifications on the visitor nights
data. The methods using forecast projection are named as “Model – Component
Weights – Number of Components”. The base models are named as “Model –
Benchmark” and these points are marked with triangles. The shaded region is
the confidence interval of the best performing model. Methods outside the
shaded region are significantly worse than the best model. Figure 4: Out-of-
sample MSE of base and FLAP forecasts as the number of components $p$
increases, for forecast horizons $1$, $6$ and $12$, using the visitor nights
data. The solid horizontal lines show the MSE for the base forecasts while the
dashed lines show the MSEs of the FLAP forecasts. The vertical black line
indicates the location of $p=m$, the number of series.

MCB and MSE plots for $h=1$, $6$ and $12$ (as described in Section 3.2) are
presented in Figure 3 and Figure 4, respectively. The FLAP forecasts in
general outperform the base forecasts, which is consistent with the findings
in Section 3. In Figure 3, the base forecasts are always ranked last.
Interestingly, FLAP forecasts improve over base forecasts, even when FLAP is
using either one or at most two principal components. These improvements are
statistically significant across all forecast horizons.

We highlight two findings that diverge from those observed in Section 3.

First, in contrast to the results from Section 3, the method implemented to
construct components matters. Using principal components in FLAP improves base
forecasts more compared to using random components with normal weights. Figure
3 shows that the best ranked forecasts are the FLAP forecasts using either the
$p=m=77$ principal components, or complementing these with random components
so that $p=200$. These improve on FLAP forecasts using only $p=200$ random
components and the improvements are statistically significant.

This is also clearly reflected in Figure 4, where the reduction in terms of
MSE from using PCA is always lower, even after the $p=m$ principal components
are exhausted and random components with weights generated from a normal
distribution are added.

PCA aims to find components along the direction where the data varies the
most. The first principal components accounts for the maximum variance in the
data. Subsequent principal components are orthogonal to the previous ones and
capture the maximum remaining variance. Based on this variance maximisation
and ranking of PCA, we propose two potential explanations for its superior
performance.

  1. 1.

Optimality: By maximising and ranking the variance of PCs from largest to
smallest, we ensure that the projection utilises components containing
significantly more information (as measured by variance) compared to randomly
weighted components.

  2. 2.

Diversity: Actively seeking principal components with the highest variance,
results in the incorporation of a more diverse set of components into the
projection.

Second, consistent with the results in Section 3, the MSEs of FLAP forecasts
shown in Figure 4 decrease as the number of components increases. However,
this decreases ceases, especially for $h=1$, approximately for $p>150$. This
is also reflected in the results presented in Figure 3, where the difference
in ranking of the two best methods is not statistically significant, even
though they use very different numbers of components ($p=m=77$ and $p=200$).

Choosing the number of components is a trade-off between the increasing
estimation error as the dimension of forecast error variance $\bm{W}_{h}$
increases, and the additional benefit brought by the information embedded in
the new components. Of course this all depends on the complexity of the DGP.
For the visitor nights data, the benefit of components over and above the
estimation error, diminishes after the number of components reaches $p=m=77$.

### 4.2 FRED-MD

The FRED-MD [McCNg2016] is a popular monthly data set of macroeconomic
variables, and shares similar properties with the [StoWat2002] data. We
downloaded and transformed the data using the fbi package [fbi]. The period we
use for this exercise is from January 1959 to September 2023, containing $777$
observations. Following [McCNg2016], we replace observations that deviate from
the sample median by more than 10 interquartile ranges (which are recognised
as outliers), with missing values. We then drop any series with more than $5$%
observations missing. This results in $m=122$ series. We fill in the missing
values using the expectation-maximization (EM) algorithm described in
[StoWat2002a] with $8$ factors. The number $8$ is identified by [McCNg2016],
albeit with a different time span.

In order to relate to the theoretical forecast error variance reduction, we
use MSE as the error measure, instead of other scaled or percentage error
measures. To reliably calculate MSE over series with different scales, we
demean the series and scale them to have variance $1$. The MSEs are calculated
on this standardised scale without back-transformation.

We evaluate the performance of forecasts using time series cross-validation.
We start with $300$ observations in the first training set and the following
$12$ observations as the test set. We repeat the exercise by expanding the
training set increasing by $1$ observation at a time. This provides us with
$466$ forecasts for each of the forecast horizons from $h=1$ to $12$ for
evaluation.

We generate base forecasts from ARIMA models using auto.arima() function from
the forecast package [forecast] with the default settings, and DFMs. The
ranges of the meta-parameters in the DFMs are $1\leq k\leq 8$ (since $8$
factors are identified and used to fill in the missing values), $1\leq n\leq
3$ and $0\leq s\leq 6$. For more details see Section 3.1. The MCB and MSE
plots for $h=1$, $6$ and $12$ are presented in Figure 5 and Figure 6,
respectively.

Figure 5: Average ranks of $1$-, $6$\- and $12$-step-ahead cross-validation
MSE of different model and component specifications on the FRED-MD data. The
methods using forecast projection are named as “Model – Component Weights –
Number of Components”. The base models are named as “Model – Benchmark” and
these points are marked with triangles. The shaded region is the confidence
interval of the best performing model. Methods outside the shaded region are
significantly worse than the best model. Figure 6: Out-of-sample MSE for base
and FLAP forecasts as the number of components $p$ increases, for forecast
horizons $1$, $6$ and $12$, using the FRED-MD data. The solid horizontal lines
show the MSE for the base forecasts while the dashed lines show the MSEs of
the FLAP forecasts. The vertical black line indicates the location of $p=m$
the number of series.

The results presented in Figure 5 show that the ARIMA base forecasts rank
better that the DFM base forecasts. This difference is statistically
significant for $h=6$. Likewise, ARIMA base forecasts return a lower MSE
compared to DFM as shown in Figure 6 across all forecast horizons.

The best ranked forecasts across all forecast horizons are once again FLAP
forecasts generated using principal components. The improvements over the
ARIMA base forecasts are significantly significant for $h=1$ and $6$.
Furthermore, FLAP forecasts using PCA rank better than FLAP forecasts using
components with random weights, and these differences are statistically
significant. This is also observed in terms of MSE in Figure 6, reaffirming
our findings based on the tourism data set discussed in Section 4.1. Figure 6
shows that both FLAP ARIMA and FLAP DFM forecasts, seem to be worse than the
base forecasts when using a small number of components $p$ for $h=6$ and $12$.
However, this improves as $p$ increases and FLAP forecasts outperform the base
forecasts gradually. The trade-off between the benefit of including more
components and the loss in degrees of freedom due to estimation in a large
dimension is once again observed. However in this case it seems to be more
extreme, compared to the tourism example, with MSEs starting to increase for
$p>m$.

As we have seen in Section 3 and Section 4.1, $m$ is not always a clear cut-
off point in determining the number of components to use. The performance of
FLAP is jointly determined by the number of series $m$, the sample size $T$,
the component construction method, and the DGP. The example of FRED-MD shows
the importance of PCA, as $m$ is the point that the component changes from PCA
to random normal weighted linear combinations, implying that PCA can exploit
the information in the data while random weights cannot. This is more obvious
for $h=1$ and $h=6$, as PCA works when $p<m$, but random normal weights do not
seem to work from the beginning.

## 5 Conclusion

The proposed forecast linear augmented projection (FLAP) method has been shown
to be a simple but effective way to reduce forecast error variance of any
multivariate forecasting problem. It simply involves augmenting the data with
linear combinations, forecasting these and then projecting the augmented
vector of forecasts. We have shown theoretically that FLAP will continue to
reduce forecast error covariance as more components are added, assuming that
the forecast error covariance matrix is known. In practice, a plug-in estimate
of this covariance matrix can be used, and in both simulated and empirical
data we demonstrate that a simple shrinkage estimator does indeed lead to
improvements in forecast accuracy. Regarding the construction of components,
we find that PCA in practice achieves significant improvements in forecast
accuracy. Another appealing property of FLAP is that the projection step can
even compensate for a poor choice of base forecasting model. This is
particularly attractive since it makes FLAP robust against model
misspecification in the base forecasting step.

One outstanding issue is to find alternatives to PCA to select component
weights. For example, [Goerg2013-dg] proposed “forecastable components” that
are optimal in the sense of minimising the forecast error variance of the
components, while [Matteson2011-jf] proposed “dynamic orthogonal components”
that reduce a multivariate time series to a set of uncorrelated univariate
time series. It would be interesting to explore whether these components (or
other similar suggestions) can be used effectively in FLAP. Another route to
improving FLAP may be found by optimizing Equation 15 over $\bm{\Phi}$ and
$\bm{G}$ rather than just $\bm{G}$.

Component construction should be studied together with the selection of the
forecast model since both the weight matrix $\bm{\Phi}$ and the base forecast
error variance $\bm{W}_{h}$ can affect the projection simultaneously in
Equation 2. This is likely to be related to forecast combination methods,
focusing on the properties of the base forecasts, and the diversity and
robustness of the forecast model and components. Examples of studies on this
issue in the forecast combination literature include [BatDua1995],
[KanEtAl2022] and [LicWin2020].

Finally, while FLAP is motivated by the forecast reconciliation literature,
our focus here is very much on multivariate time series with no constraints.
However, it would be possible to use both forecast reconciliation and forecast
projection together. This may be particularly useful when there are
relationships between series that are not captured in the known hierarchical
structure.

## Acknowledgements

We thank Daniele Girolimetto for contributing to the initial proof of Theorem
2.2 in his unpublished work.

## References

  * [1] Sakai Ando and Futoshi Narita  “An Alternative Proof of Minimum Trace Reconciliation”  In _Forecasting_ 6 Multidisciplinary Digital Publishing Institute, 2024, pp. 456–461  DOI: 10.3390/forecast6020025
  * [2] V Assimakopoulos and K Nikolopoulos  “The theta model: a decomposition approach to forecasting”  In _International Journal of Forecasting_ 16 Elsevier, 2000, pp. 521–530  DOI: 10.1016/S0169-2070(00)00066-2
  * [3] George Athanasopoulos, Rob J Hyndman, Nikolaos Kourentzes and Anastasios Panagiotelis  “Forecast reconciliation: A review”  In _International Journal of Forecasting_ Elsevier BV, 2023  DOI: 10.1016/j.ijforecast.2023.10.010
  * [4] George Athanasopoulos, Rob J Hyndman, Nikolaos Kourentzes and Fotios Petropoulos  “Forecasting with temporal hierarchies”  In _European Journal of Operational Research_ 262 Elsevier BV, 2017, pp. 60–74  DOI: 10.1016/j.ejor.2017.02.046
  * [5] Roy Batchelor and Pami Dua  “Forecaster diversity and the benefits of combining forecasts”  In _Management Science_ 41 Institute for Operations Researchthe Management Sciences (INFORMS), 1995, pp. 68–75  DOI: 10.1287/mnsc.41.1.68
  * [6] Christoph Bergmeir, Rob J Hyndman and José M Benítez  “Bagging exponential smoothing methods using STL decomposition and Box–Cox transformation”  In _International Journal of Forecasting_ 32, 2016, pp. 303–312  DOI: 10.1016/j.ijforecast.2015.07.002
  * [7] Ben S Bernanke, Jean Boivin and Piotr Eliasz  “Measuring the effects of monetary policy: a factor-augmented vector autoregressive (FAVAR) approach”  In _The Quarterly journal of economics_ 120.1 MIT Press, 2005, pp. 387–422
  * [8] Hans W. Borchers  “pracma: Practical Numerical Math Functions” R package version 2.4.4, 2023  URL: https://CRAN.R-project.org/package=pracma
  * [9] Leo Breiman  “Bagging predictors”  In _Machine learning_ 24 Springer Nature, 1996, pp. 123–140  DOI: 10.1007/bf00058655
  * [10] Andrea Carriero, Ana Beatriz Galvao and George Kapetanios  “A comprehensive evaluation of macroeconomic forecasting methods”  In _International Journal of Forecasting_ 35.4 Elsevier, 2019, pp. 1226–1239
  * [11] Yankang Chen, Serena Ng and Jushan Bai  “fbi: Factor-Based Imputation and FRED-MD/QD Data Set” R package version 0.7.0, 2023  URL: https://github.com/cykbennie/fbi
  * [12] Jacopo De Stefani et al.  “Batch and incremental dynamic factor machine learning for multivariate and multi-step-ahead forecasting”  In _International Journal of Data Science and Analytics_ 7 Springer ScienceBusiness Media LLC, 2019, pp. 311–329  DOI: 10.1007/s41060-018-0150-x
  * [13] Tommaso Di Fonzo and Daniele Girolimetto  “Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives”  In _International Journal of Forecasting_ 39 Elsevier BV, 2023, pp. 39–57  DOI: 10.1016/j.ijforecast.2021.08.004
  * [14] Stephen Michael Disney and Fotios Petropoulos  “Forecast combinations using multiple starting points”, 2015
  * [15] Antonio Fabio Di Narzo, Jose Luis Aznarte and Matthieu Stigler  “tsDyn: Time series analysis based on dynamical systems theory” R package version 0.7, 2009  URL: https://CRAN.R-project.org/package=tsDyn
  * [16] Milton Friedman  “The Use of Ranks to Avoid the Assumption of Normality Implicit in the Analysis of Variance”  In _Journal of the American Statistical Association_ 32 Taylor & Francis, 1937, pp. 675–701  DOI: 10.1080/01621459.1937.10503522
  * [17] Milton Friedman  “A Correction”  In _Journal of the American Statistical Association_ 34 Informa UK Limited, 1939, pp. 109–109  DOI: 10.1080/01621459.1939.10502372
  * [18] Georg Goerg  “Forecastable Component Analysis”  In _Proceedings of The 30th International Conference on Machine Learning_ , 2013, pp. 64–72  URL: http://jmlr.org/proceedings/papers/v28/goerg13.pdf
  * [19] Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The elements of statistical learning: Data mining, inference, and prediction”, Springer series in statistics  New York, NY: Springer, 2003  DOI: 10.1007/978-0-387-21606-5
  * [20] Myles Hollander, Douglas A Wolfe and Eric Chicken  “Nonparametric Statistical Methods”  John Wiley & Sons, 2013
  * [21] Ross Hollyman, Fotios Petropoulos and Michael E Tipping  “Understanding forecast reconciliation”  In _European Journal of Operational Research_ 294 Elsevier BV, 2021, pp. 149–160  DOI: 10.1016/j.ejor.2021.01.017
  * [22] Rob Hyndman et al.  “forecast: Forecasting functions for time series and linear models” R package version 8.21.1, 2023  URL: https://pkg.robjhyndman.com/forecast/
  * [23] Rob J Hyndman and George Athanasopoulos  “Forecasting: principles and practice”  Melbourne, Australia: OTexts, 2018  URL: http://OTexts.org/fpp2
  * [24] Rob J Hyndman and George Athanasopoulos  “Forecasting: principles and practice”  Melbourne, Australia: OTexts, 2021  URL: http://OTexts.org/fpp3
  * [25] Rob J Hyndman and Yeasmin Khandakar  “Automatic time series forecasting: The forecast package for R”  In _Journal of Statistical Software_ 27 Foundation for Open Access Statistic, 2008, pp. 1–22  DOI: 10.18637/jss.v027.i03
  * [26] I T Jolliffe  “Principal Component Analysis”  Springer, New York, NY, 2002  DOI: 10.1007/b98835
  * [27] Yanfei Kang, Wei Cao, Fotios Petropoulos and Feng Li  “Forecast with forecasts: Diversity matters”  In _European Journal of Operational Research_ 301 Elsevier BV, 2022, pp. 180–190  DOI: 10.1016/j.ejor.2021.10.024
  * [28] Alex J Koning, Philip Hans Franses, Michèle Hibon and H O Stekler  “The M3 competition: Statistical tests of the results”  In _International Journal of Forecasting_ 21 Elsevier, 2005, pp. 397–409  DOI: 10.1016/j.ijforecast.2004.10.003
  * [29] Nikolaos Kourentzes  “tsutils: Time Series Exploration, Modelling and Forecasting” R package version 0.9.4, 2023  URL: https://CRAN.R-project.org/package=tsutils
  * [30] Nikolaos Kourentzes, Fotios Petropoulos and Juan R Trapero  “Improving forecasting by estimating time series structural components across multiple frequencies”  In _International Journal of Forecasting_ 30 Elsevier BV, 2014, pp. 291–302  DOI: 10.1016/j.ijforecast.2013.09.006
  * [31] Xixi Li, Fotios Petropoulos and Yanfei Kang  “Improving forecasting by subsampling seasonal time series”  In _International Journal of Production Research_ 61 Informa UK Limited, 2022, pp. 1–17  DOI: 10.1080/00207543.2021.2022800
  * [32] Kenneth C Lichtendahl and Robert L Winkler  “Why do some combinations perform better than others?”  In _International Journal of Forecasting_ 36 Elsevier BV, 2020, pp. 142–149  DOI: 10.1016/j.ijforecast.2019.03.027
  * [33] David G Luenberger  “Optimization by vector space methods”  Nashville, TN: John Wiley & Sons, 1969
  * [34] David S Matteson and Ruey S Tsay  “Dynamic Orthogonal Components for Multivariate Time Series”  In _Journal of the American Statistical Association_ 106.496, 2011, pp. 1450–1463  DOI: 10.1198/jasa.2011.tm10616
  * [35] Michael W McCracken and Serena Ng  “FRED-MD: A Monthly Database for Macroeconomic Research”  In _Journal of Business & Economic Statistics_ 34 Taylor & Francis, 2016, pp. 574–589  DOI: 10.1080/07350015.2015.1086655
  * [36] Peter Bjorn Nemenyi  “Distribution-free multiple comparisons”  ProQuest DissertationsTheses, 1963
  * [37] Rainer Opgen-Rhein and Korbinian Strimmer  “Accurate ranking of differentially expressed genes by a distribution-free shrinkage approach”  In _Statistical Applications in Genetics and Molecular Biology_ 6 Walter de Gruyter GmbH, 2007  DOI: 10.2202/1544-6115.1252
  * [38] Anastasios Panagiotelis, George Athanasopoulos, Puwasala Gamakumara and Rob J Hyndman  “Forecast reconciliation: A geometric view with new insights on bias correction”  In _International Journal of Forecasting_ 37 Elsevier BV, 2021, pp. 343–359  DOI: 10.1016/j.ijforecast.2020.06.004
  * [39] Fotios Petropoulos, Rob J Hyndman and Christoph Bergmeir  “Exploring the sources of uncertainty: Why does bagging for time series forecasting work?”  In _European Journal of Operational Research_ 268 Elsevier BV, 2018, pp. 545–554  DOI: 10.1016/j.ejor.2018.01.045
  * [40] Fotios Petropoulos and Evangelos Spiliotis  “The wisdom of the data: Getting the most out of univariate time series forecasting”  In _Forecasting_ 3 MDPI AG, 2021, pp. 478–497  DOI: 10.3390/forecast3030029
  * [41] R Core Team  “R: A Language and Environment for Statistical Computing”, 2023  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [42] C Radhakrishna Rao  “Projectors, generalized inverses and the Blue’S”  In _Journal of the Royal Statistical Society. Series B, Statistical methodology_ 36 Oxford University Press (OUP), 1974, pp. 442–448  DOI: 10.1111/j.2517-6161.1974.tb01019.x
  * [43] Juliane Schafer et al.  “corpcor: Efficient Estimation of Covariance and (Partial) Correlation” R package version 1.6.10, 2021  URL: https://CRAN.R-project.org/package=corpcor
  * [44] Juliane Schäfer and Korbinian Strimmer  “A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics”  In _Statistical Applications in Genetics and Molecular Biology_ 4 Walter de Gruyter GmbH, 2005  DOI: 10.2202/1544-6115.1175
  * [45] James H Stock and Mark W Watson  “Forecasting Using Principal Components From a Large Number of Predictors”  In _Journal of the American Statistical Association_ 97 Taylor & Francis, 2002, pp. 1167–1179  DOI: 10.1198/016214502388618960
  * [46] James H Stock and Mark W Watson  “Macroeconomic Forecasting Using Diffusion Indexes”  In _Journal of Business & Economic Statistics_ 20 Taylor & Francis, 2002, pp. 147–162  DOI: 10.1198/073500102317351921
  * [47] James H Stock and Mark W Watson  “Generalized shrinkage methods for forecasting using many predictors”  In _Journal of Business & Economic Statistics_ 30 Informa UK Limited, 2012, pp. 481–493  DOI: 10.1080/07350015.2012.715956
  * [48] Ruey S Tsay  “Multivariate time series analysis: with R and financial applications”  John Wiley & Sons, 2013
  * [49] Xiaoqian Wang, Rob J Hyndman, Feng Li and Yanfei Kang  “Forecast combinations: An over 50-year review”  In _International Journal of Forecasting_ 39 Elsevier BV, 2023, pp. 1518–1547  DOI: 10.1016/j.ijforecast.2022.11.005
  * [50] Shanika L Wickramasuriya, George Athanasopoulos and Rob J Hyndman  “Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization”  In _Journal of the American Statistical Association_ 114 Informa UK Limited, 2019, pp. 804–819  DOI: 10.1080/01621459.2018.1448825
  * [51] Yangzhuoran Fin Yang  “flap: Forecast Linear Augmented Projection” R package version 0.1.0, 2024  URL: https://CRAN.R-project.org/package=flap