The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
114,266 characters
Forecasting with panel data: Estimation uncertainty versus parameter heterogeneity
\title{Forecasting with panel data: Estimation uncertainty versus parameter
heterogeneity\thanks{
We thank the three anonymous reviewers and the Co-editor, Stephane Bonhomme, for their helpful and
constructive comments. We also thank Laura Liu, Mahrad Sharifvaghefi, Ron
Smith, Cynthia Yang, Liying Yang, and seminar participants at the University of Pittsburgh, the Board of the Federal Reserve, ESEM, VTSS, IAAE, and at 2024
NBER-NSF Time Series Conference held at University of Pennsylvania for
helpful comments.}}
\author{M. Hashem Pesaran\thanks{
University of Cambridge, UK, and University of Southern California, USA
Email: [email removed]} \and Andreas Pick\thanks{
Erasmus University Rotterdam, Erasmus School of Economics, Burgemeester
Oudlaan 50, 3000DR Rotterdam, and Tinbergen Institute. Email:
[email removed]} \and Allan Timmermann\thanks{
UC San Diego, Rady School of Management, 9500 Gilman Drive, La Jolla CA
92093-0553. Email: [email removed].}}
\maketitle
\begin{abstract}
We provide a comprehensive examination of the predictive accuracy of panel
forecasting methods based on individual, pooling, fixed effects, and empirical
Bayes estimation, and propose optimal weights for forecast combination
schemes. We consider linear panel data models, allowing for weakly exogenous
regressors and correlated heterogeneity. We quantify the gains from exploiting
panel data and demonstrate how forecasting performance depends on the degree
of parameter heterogeneity, whether such heterogeneity is correlated with
the regressors, the goodness-of-fit of the model, and the dimensions of
the data. Monte Carlo simulations and empirical applications to house prices
and CPI inflation show that empirical Bayes and forecast combination methods
perform best overall and rarely produce the least accurate forecasts for
individual series. \enlargethispage{3em} \newline
\emph{JEL codes: C33, C53}\newline
\emph{Keywords: Forecasting, Panel data, Heterogeneity, Pooled estimation,
Empirical Bayes; Forecast combination.}
\end{abstract}
\thispagestyle{empty}
\newpage
\setcounter{page}{1}
\section{Introduction}
Panel data sets on economic and financial variables are widely available
at individual, firm, industry, regional, and country granularities and
have been extensively used for estimation and inference. Yet, panel estimation
methods have had a comparatively lower impact on common practices in economic
forecasting, which remain dominated by unit-specific forecasting models
or low-dimensional multivariate models such as vector autoregressions \citep{Hsi2022}.
The relative shortage of panel applications in the economic forecasting
literature is, in part, a result of the absence of a deeper understanding
of the determinants of forecasting performance for different panel estimation
methods and the absence of guidelines on which methods work well in different
settings.
In this paper, we examine existing approaches and develop novel forecast
combination methods for panel data with possibly correlated heterogeneous
parameters. We conduct a systematic comparison of their predictive accuracy
in settings with different cross-sectional ($N$) and time ($T$) dimensions
and varying degrees of parameter heterogeneity, whether correlated or not.
Our analysis provides a deeper understanding of the determinants of the
performance of these methods across a variety of settings chosen for their
relevance to economic forecasting problems. This includes the important
choice of whether to use pooled versus individual estimates, or perhaps
a combination of the two approaches, with a focus on forecasting rather
than parameter estimation and inference.
We begin by exploring analytically the bias-variance trade-off between
individual, fixed effects (FE), and pooled estimation for forecasting.
Our analysis is conducted in a general setting that allows for weakly exogenous
regressors and correlated heterogeneity, consistent with the type of dynamic
panel models commonly used in empirical applications. We show how such
effects contribute to the mean squared forecast error (MSFE) of forecasts
based on individual, FE, and pooled estimates.
We next examine forecast combination methods. Estimation errors are well
known to lead to imprecisely estimated combination weights for data with
a small time-series dimension. Our main combination scheme assumes homogeneous
weights across individual variables. This allows us to use cross-sectional
information to reduce the effect of estimation error on the combination
weights compared to the conventional combination scheme that lets the weights
be individual-specific, which we also consider. To handle cases where the
pooling estimator imposes too much homogeneity, we also consider combinations
based on forecasts from the individual-specific and fixed effect estimators.
Our theoretical analysis of the individual and pooled estimation schemes
focuses on the case with finite $T$ and $N\rightarrow \infty $ and does
not require that $\sqrt{N}/T\rightarrow 0$ as $N$ and
$T\rightarrow \infty $, jointly, which is often assumed in the literature.
The estimation of the combination weights, however, requires $T$
$\rightarrow \infty $, but at a much slower rate compared to $N$.
Finally, we consider forecasts based on the empirical Bayesian (EB) approach
of \cite{Hsietal1999}. These are related to forecast combination and we
show for the empirical Bayes estimator that it can be thought of as a weighted
average of an estimator that allows for full heterogeneity and a pooled
mean group estimator. The empirical Bayes scheme assigns greater weight
to the pooled estimator, the lower the estimated degree of parameter heterogeneity
and so adapts to the degree of parameter heterogeneity characterizing a
given data set.\footnote{In the Supplemental Appendix, we
also report results based on the hierarchical Bayesian approach of \citet{LinSmi1972},
\cite{LeeGri1979}, and \cite{Madetal1997}.}
We evaluate the predictive accuracy of these alternative panel forecasting
methods through Monte Carlo simulations. The simulations explore the importance
to forecasting performance of the degree of parameter heterogeneity, how
correlated it is with the regressors, whether it affects intercepts or slopes, the value of the regressors in the forecast period, and dimensions of $N$ and $T$. In
the scenario with homogeneous parameters, forecasts based on pooled estimates
are most accurate. Forecasts based on fixed or random effect estimates
perform well, relative to other methods, when parameter heterogeneity is
confined to the intercepts and does not affect slopes. Outside these cases,
empirical Bayes and forecast combinations produce the most accurate forecasts
and are better able to handle parameter heterogeneity, whether correlated
or not, while being more robust in cases with a small $T$ than the individual-specific
approach.
Next, we consider two empirical applications selected to represent varying
degrees of heterogeneity and predictive power of the underlying forecasting
models. Our first application considers predictability of house prices
across 362 US metropolitan statistical areas (MSAs). In this application,
individual-specific forecasts perform quite poorly, producing the highest
MSFE values among all methods for more than 50\% of the MSAs. Forecasts
based on pooled estimates perform notably better and, across all forecasts,
reduce the average MSFE value by 8\% relative to the forecasts based on
individual estimates, though this gets reversed if the regressor set is
close to the sample average. Empirical Bayes and forecast combinations
work even better in this application, beating forecasts based on individual
estimates for over 90\% of MSAs while almost never generating the least
accurate forecasts for individual series.
Our second application considers forecasts for a panel containing 187 subcategories
of CPI inflation. In this application, forecasts based on individual estimates
generate the highest MSFE values for 44\% of the series. Forecasts based
on pooled estimates produce the highest MSFE values for 40\% of the individual
series but, conversely, generate the lowest MSFE-values for 20\% of the
series. Combination forecasts produce lower MSFE values for 73--97\% of
the individual CPI series than the individual forecasts while almost never
generating the largest MSFE value. Even better inflation forecasts are
produced by the empirical Bayes method, which is more accurate (in the
MSFE sense) than the individual forecasts for 98\% of the series and generates
the lowest MSFE values for 39\% of the individual variables while never
generating the largest MSFE value.
Overall, forecasts that use only the information on a given unit tend to
have loss distributions with wide dispersion across units. Their associated
forecasts are therefore sometimes the best but far more often the worst,
and their distribution of MSFE performance is often shifted to the right,
implying larger losses on average than for other methods. Forecasts based
on pooling, random effects, or fixed effects estimation tend to perform
better, on average, than the individual forecasts whose accuracy they beat
for the majority of series. However, relative to the individual-specific
forecasts, these approaches also tend to have a right-skewed MSFE distribution,
suggesting a high risk of poor forecasting performance for individual series
whose model parameters are very different from the average. Combinations
and empirical Bayes forecasts have narrower MSFE distributions across units,
often shifted to the left as they are centered around a smaller average
loss, and rarely produce the largest squared forecast error among all methods
we consider.
\paragraph*{Related literature} The review articles by \citeauthor{Bal2008} (\citeyear{Bal2008,Bal2013})
consider the forecasting performance of the best linear unbiased predictor
(BLUP) of \citet{Gol1962} in models with either fixed effects or random
effects. The BLUP estimator gives rise to a generalized least squares (GLS)
predictor, which Baltagi compares to models that allow for autoregressive
moving average (ARMA) dynamics in innovations as well as models with spatial
dependencies in the errors. \citet{TraUrg2009} use Monte Carlo simulations
to assess the forecasting performance of pooled, individual, and shrinkage
estimators and find that parameter heterogeneity is a key determinant of
the accuracy of different forecasts. \citet{BruSil2006}
consider a similar group of methods to forecast migration data and find
that fixed effects and shrinkage estimators perform best; see \citet{PicTim2024} for a review of the literature.
\cite{Wanetal2019} also propose forecast combination methods. However,
their analysis does not allow for correlation of regressors and parameters
or dynamics in the model. Additionally, their combination weights are determined
from in-sample test statistics rather than the expected out-of-sample performance
that we propose. In this sense, our approach is closer to the forecast-based
test for a structural break of \cite{Pesetal2013} and \citet{BooPic2020},
where the target is also significant improvements in forecast accuracy
rather than a significant change in parameters.
\citet{Liuetal2020} study forecasting for dynamic panel data
models with a short time-series dimension. Though $T$ exceeds the number
of parameters that have to be estimated for each series, such estimates
are typically very noisy and not consistent under large $N$, fixed
$T$ asymptotics. To handle estimation noise, they adopt a nonparametric
Bayesian approach that shrinks the heterogeneous parameters towards local
patterns in the distribution. This is closely related to the idea of using
forecast combinations to reduce the effect on the forecasts of noisy estimates
of individual-specific parameters.
\paragraph*{Outline} The rest of the paper is organized as follows. Section~\ref{SETUP} introduces the model setup and our assumptions, while Section~\ref{sec:theory_msfe} derives
analytical results on the predictive accuracy of individual, pooled, and
FE forecasting schemes. Section~\ref{sec:combinations} introduces our forecast combination schemes.
Section~\ref{sec:MC} describes the empirical Bayes estimator. Section~\ref{sec:applications} presents Monte
Carlo experiments, Section~7 reports our empirical applications,
and Section~8 concludes. Technical details are provided in appendices at the end of
the paper and in the Supplemental Appendix.
\section{Setup and assumptions\label{SETUP}}
We begin by describing the panel regression setup and assumptions used
in our analysis.
\subsection{Panel regression model\label{PanelReg}}
Our analysis considers the following linear panel regression model:
\begin{equation}
y_{it}=\alpha _{i}+\boldsymbol{\beta
}_{i}^{\prime }\boldsymbol{x}
_{it}+
\varepsilon _{it}=\boldsymbol{\theta }_{i}^{\prime }
\boldsymbol{w}
_{it}+\varepsilon _{it},\quad
\varepsilon _{it}\sim \bigl(0,\sigma _{i}^{2}
\bigr), \label{eq:model1}
\end{equation}
where $i=1,2,\ldots ,N$ refers to the individual units and
$t=1,2,\ldots ,T$ refers to the time period, $y_{it}$ is the outcome of
unit $i$ at time $t$, $\boldsymbol{x}_{it}$ is a $k\times 1$ vector of
regressors---or predictors---used
to forecast $y_{it}$ (including, possibly, latent factors),
$\boldsymbol{\beta }_{i}$ is the associated vector of regression coefficients,
and $\varepsilon _{it}$ is the disturbance of unit $i$ in period $t$. The
second equality in (\ref{eq:model1}) introduces the notation
$\boldsymbol{\theta }_{i}=(\alpha _{i},\boldsymbol{\beta }
_{i}^{\prime })^{\prime }$ and
$\boldsymbol{w}_{it}=(1,\boldsymbol{x}
_{it}^{\prime })^{\prime }$, which have dimensions $K\times 1$, with
$K=k+1$. For simplicity, we use the time subscript $t$ for
$\boldsymbol{x}_{it}$ and $
\boldsymbol{w}_{it}$, but it is important to emphasize that this refers
to the predicted time for the outcome variable, $y_{it}$. For a forecast
horizon of $h$ periods, all variables in $\boldsymbol{x}_{it}$ must therefore
be known at time $t-h$. Our notation avoids explicitly referring to
$h$ everywhere, but it should be recalled throughout the analysis that
$\boldsymbol{x}_{it}$ includes suitably lagged predictors. We will focus
on the case of $h=1$ but extensions to larger $h$ are straightforward.
\paragraph*{Notation} Stacking the time series of outcomes, regressors, and
disturbances, define
$\boldsymbol{y}_{i}=(y_{i1},y_{i2},\ldots ,y_{iT})^{\prime }$,
$\boldsymbol{X}_{i}=(\boldsymbol{x}_{i1},
\boldsymbol{x}_{i2},\ldots ,\boldsymbol{x}_{iT})^{
\prime }$,
$\boldsymbol{W}_{i}= ( \boldsymbol{\tau }_{T},\boldsymbol{X}
_{i} ) $, where $\boldsymbol{\tau }_{T}$ is a $T\times 1$ vector
of ones, and
$\boldsymbol{\varepsilon }_{i}=(\varepsilon _{i1},\varepsilon _{i2},
\ldots ,\varepsilon _{iT})^{\prime }$. Further, let
$\boldsymbol{y}=(
\boldsymbol{y}_{1}^{\prime },\boldsymbol{y}_{2}^{\prime },\ldots ,
\boldsymbol{y}_{N}^{\prime })^{\prime }$,
$\boldsymbol{X}=(\boldsymbol{X}
_{1}^{\prime },\boldsymbol{X}_{2}^{\prime },\ldots ,\boldsymbol{X}
_{N}^{\prime })^{\prime }$,
$\boldsymbol{W}=(\boldsymbol{W}
_{1}^{\prime },\boldsymbol{W}_{2}^{\prime },\ldots ,\boldsymbol{W}
_{N}^{\prime })^{\prime }$, and
$\boldsymbol{\varepsilon }=(\boldsymbol{
\varepsilon }_{1}^{\prime },\boldsymbol{\varepsilon }_{2}^{\prime },
\ldots ,\boldsymbol{\varepsilon }_{N}^{\prime })^{\prime }$. Generic positive finite
constants are denoted by $C$ when large and $c$ when small. They can take
different values at different instances.
$\lambda _{\max } ( \boldsymbol{
A} ) $ and $\lambda _{\min } ( \boldsymbol{A} ) $ denote
the maximum and minimum eigenvalues of matrix $\boldsymbol{A}$.
$\boldsymbol{A}
\succ \boldsymbol{0}$ and $\boldsymbol{A}\succeq \boldsymbol{0}$ denote
that $\boldsymbol{A}$ is a positive definite and a nonnegative definite
matrix, respectively.
$ \Vert \boldsymbol{A} \Vert =\lambda _{\max }^{1/2}(
\boldsymbol{A}^{\prime }\boldsymbol{A)}$ and
$ \Vert \boldsymbol{A}
\Vert _{1}$ denote the spectral and column sum norms of matrix
$
\boldsymbol{A}$, respectively,
$ \Vert \boldsymbol{x} \Vert _{p}=
[ \mathrm{E} ( \Vert \boldsymbol{x} \Vert ^{p}
) ] ^{1/p}$. If
$ \{ f_{n} \} _{n=1}^{\infty }$ is any real sequence and
$ \{ g_{n} \} _{n=1}^{\infty }$ is a sequence of positive real
numbers, then $f_{n}=O(g_{n})$, if there exists a $C$ such that
$ \vert f_{n} \vert /g_{n}\leq C$ for all $n$ and
$f_{n}=o(g_{n})$ if $f_{n}/g_{n}\rightarrow 0$ as
$n\rightarrow \infty $. Similarly, $f_{n}=O_{p}(g_{n})$ if
$f_{n}/g_{n}$ is stochastically bounded and $f_{n}=o_{p}(g_{n})$ if
$f_{n}/g_{n}\overset{p}{\rightarrow}0$. The operator
$\overset{p}{\rightarrow}$ denotes convergence in probability, and
$\overset{d}{\rightarrow}$ denotes convergence in distribution.
\subsection{Assumptions}
Our theoretical analysis builds on a set of standard assumptions about
the underlying data generating process.
\begin{assumption}
\label{ass:1}
$\varepsilon _{it}$ is serially independent with mean zero, a fixed variance
$\sigma _{i}^{2}$ $(0<c<\sigma _{i}^{2}<C<\infty)$, and with
$\sup_{i,t}\mathrm{E} \vert \varepsilon _{it} \vert ^{4}<C<
\infty $.
\end{assumption}
\begin{assumption}
\label{ass:weak_exogeneity_a}
$ \{ \varepsilon _{it} \} $ for $i$ $
=1,2,\ldots ,N$ are martingale difference processes with respect to the
filtration,
$\mathcal{I}_{it}= ( \boldsymbol{w}_{it},\boldsymbol{w}
_{i,t-1},\ldots ) $, so that
\begin{equation*}
\mathrm{E}\left( \varepsilon _{it}\vert \boldsymbol{w}_{is}
\right. ) =0,\quad \text{for }t\geq s,\text{ for }t=1,2,\dots ,T,T+1.
\end{equation*}
\end{assumption}
\begin{assumption}
\label{ass:3}
(a) $ \{ \boldsymbol{w}_{it} \} $ for $i=1,2,\ldots ,N $ are
covariance stationary with
$\mathrm{E}(\boldsymbol{w}_{it}
\boldsymbol{w}_{it}^{\prime })=\boldsymbol{Q}_{i}$,
$\sup_{i,t= \{ 1,2,\ldots ,T \} }\mathrm{E} \Vert
\boldsymbol{w}_{it} \Vert ^{4}<C$,
$\sup_{i,T} \Vert \boldsymbol{w}_{i,T+1} \Vert <C$, and
\begin{equation}
\sup_{i}\mathrm{\lambda }_{\max } (
\boldsymbol{Q}_{i } ) <C<\infty,\quad \text{and}\quad \sup
_{i}\mathrm{\lambda }_{\max } \bigl(
\boldsymbol{Q}
_{i }^{-1} \bigr) <C<\infty
. \label{CQi}
\end{equation}
(b) The sample covariance matrices
$\boldsymbol{Q}_{iT }=T^{-1}\boldsymbol{W
}_{i}^{\prime }\boldsymbol{W}_{i}=T^{-1}\sum_{t=1}^{T}\boldsymbol{w}_{it}
\boldsymbol{w}_{it}^{\prime }$, for $i=1,2,\ldots ,N$, satisfy the conditions
$\sup_{i}\mathrm{\lambda }_{\max } ( \boldsymbol{Q}_{iT
} ) <C<\infty $, and
$\sup_{i}\mathrm{\lambda }_{\max } ( \boldsymbol{Q}_{iT }^{-1}
) <C<\infty $.
\end{assumption}
\begin{assumption}
\label{ass:weak_exogeneity_b}
There exists a fixed $T_{0}$ such that for all $T>T_{0}$,
\begin{align}
\sup_{i}\mathrm{E} \bigl\Vert T^{-1/2}
\boldsymbol{W}_{i}^{\prime } \boldsymbol{
\varepsilon
}_{i} \bigr\Vert ^{4}&<C<\infty , \label{CWe}
\\
\sup_{i}\mathrm{E} \bigl[ \lambda _{\max }^{4}
( \boldsymbol{Q}_{iT
} ) \bigr] &<C<\infty, \quad \text{and}\quad \sup
_{i}\mathrm{E} \bigl[ \lambda _{\max }^{4}
\bigl( \boldsymbol{Q}_{iT }^{-1} \bigr) \bigr] <C<\infty .
\label{CQiT}
\end{align}
\end{assumption}
Under Assumption~\ref{ass:1}, the optimal forecast of $y_{i,T+1}$, in a
mean squared error sense, is given by
$\mathrm{E} ( y_{i,T+1} \vert \boldsymbol{w}_{i,T+1},
\boldsymbol{W}_{i} ) =\boldsymbol{\theta }_{i}^{
\prime }\boldsymbol{w}_{i.T+1}$. Note that $\boldsymbol{w}_{i.T+1}$ is
known at time $T$, and is bounded under Assumption~\ref{ass:3}. Assumption~\ref{ass:weak_exogeneity_a} allows the regressors to be weakly exogenous
with respect to $\boldsymbol{\varepsilon }_{i}$ and, therefore, permits
the inclusion of lagged dependent variables such as $y_{i,T}$ in
$
\boldsymbol{w}_{i,T+1}$. Part (a) of Assumption~\ref{ass:3} is standard
in the forecasting literature and requires the regressors to be stationary.
Part (b) is an identification assumption that allows estimation of individual
slope coefficients, $\boldsymbol{\theta }_{i}$, by least squares. Assumption~\ref{ass:weak_exogeneity_b} is required when we compare average MSFEs based
on individual and pooled estimators. It provides sufficient conditions
under which (see Lemma~\ref{Lemma_2_VTEX1})
\begin{equation}
\mathrm{E} \bigl\Vert \sqrt{T} ( \hat{\boldsymbol{\theta }_{i}}-
\boldsymbol{\theta }_{i} ) \bigr\Vert ^{2}=\mathrm{E}
\bigl\Vert \boldsymbol{Q}_{iT }^{-1} \bigl(
T^{-1/2}\boldsymbol{W}_{i}^{
\prime }
\boldsymbol{\varepsilon }_{i} \bigr) \bigr\Vert ^{2}<C<
\infty , \label{Momentconthetai}
\end{equation}
where
$\hat{\boldsymbol{\theta }_{i}}= ( \boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}_{i} ) ^{-1}\boldsymbol{W}_{i}^{\prime }
\boldsymbol{y}
_{i} $ is the least squares estimator of $\boldsymbol{\theta }_{i}$. The
moment conditions in Assumption~\ref{ass:weak_exogeneity_b} can be relaxed
when $\boldsymbol{w}_{it}$ is strictly exogenous.
Under weakly exogenous regressors the least squares estimator has a small
$T$ bias, and
$\mathrm{E} ( \hat{\boldsymbol{\theta }}_{i}-
\boldsymbol{
\theta }_{i} ) =O ( T^{-1} ) $. Under strictly exogenous
regressors, in contrast,
$\mathrm{E} ( \hat{\boldsymbol{\theta }_{i}}-
\boldsymbol{\theta }_{i} ) =\mathbf{0}$. We also note that, under
Assumptions \ref{ass:3} and \ref{ass:weak_exogeneity_b},
$ \Vert \boldsymbol{Q}_{iT }-\boldsymbol{Q}_{i} \Vert =O_{p}(T^{-1/2})$,
and
$
\Vert \boldsymbol{Q}_{iT }^{-1}-\boldsymbol{Q}_{i}^{-1} \Vert =O_{p}(T^{-1/2})$. These results, which hold for each $i$, are used
in the implementation of our combination forecasts below. For proof of
consistency of the weights in the combined forecasts discussed in Section~\ref{sec:combinations}
below, we need the stronger conditions
\begin{equation}
\sup_{i} \Vert \boldsymbol{Q}_{iT }-
\boldsymbol{Q}_{i} \Vert =O_{p} \biggl( \frac{\ln (N)}{
\sqrt{T}} \biggr) ,\quad \text{and}\quad \sup_{i} \bigl\Vert
\boldsymbol{Q}_{iT }^{-1}-\boldsymbol{Q}_{i}^{-1}
\bigr\Vert =O_{p} \biggl( \frac{\ln (N)}{\sqrt{T}} \biggr) ,
\label{supQinv}
\end{equation}
still allowing $N$ to rise much faster than $T$.\footnote{As noted by
\citeauthor{Fanetal2015} (\citeyear{Fanetal2015}, Section~\ref{ForecastsIP}), this stronger condition is
typically satisfied for strictly stationary data that satisfy strong mixing
conditions.}
Finally, let
$\boldsymbol{g}_{it}=\boldsymbol{w}_{it}\varepsilon _{it}$, and note that
$T^{-1/2}\boldsymbol{W}_{i}^{\prime }\boldsymbol{\varepsilon }
_{i}=T^{-1/2}\sum_{t=1}^{T}\boldsymbol{g}_{it}$. Also, under Assumption~\ref{ass:weak_exogeneity_a}
$\boldsymbol{g}_{it}$ is a martingale difference process with respect to
$\mathcal{I}_{it}= ( \boldsymbol{w}_{it},
\boldsymbol{w}_{i,t-1},\ldots ) $, and we have
$\mathrm{E} ( \boldsymbol{g}_{it} ) =\mathbf{0}$,
\begin{equation*}
\operatorname{Var} \Biggl( T^{-1/2}\sum_{t=1}^{T}
\boldsymbol{g}_{it} \Biggr) =T^{-1} \sum
_{t=1}^{T}\mathrm{E} \bigl( \boldsymbol{g}_{it}
\boldsymbol{g}
_{it}^{\prime } \bigr)
=T^{-1}\mathrm{E} \bigl( \boldsymbol{W}_{i}^{
\prime }
\boldsymbol{\varepsilon }_{i}\boldsymbol{\varepsilon
}_{i}^{\prime }
\boldsymbol{W}_{i}
\bigr) =T^{-1}\sum_{t=1}^{T}
\sigma _{i}^{2} \mathrm{E}
\bigl(
\boldsymbol{w}_{it}\boldsymbol{w}_{it}^{\prime }
\bigr) .
\end{equation*}
Further, under Assumption~\ref{ass:3},
$\mathrm{E} ( \boldsymbol{w}_{it}
\boldsymbol{w}_{it}^{\prime } ) =\boldsymbol{Q}_{i}$, and it follows
that
\begin{equation}
\mathrm{E} \bigl( T^{-1}\boldsymbol{W}_{i}^{\prime }
\boldsymbol{\varepsilon }
_{i}\boldsymbol{\varepsilon
}_{i}^{\prime }\boldsymbol{W}_{i} \bigr) = \sigma
_{i}^{2}\boldsymbol{Q}_{i}.
\label{Gi}
\end{equation}
We next introduce assumptions that are required primarily for establishing
the properties of pooled and fixed effects predictors.
\begin{assumption}
\label{ass:4}
(a)
$\boldsymbol{\theta}_{i}=\boldsymbol{\theta}+\boldsymbol{
\eta }_{i}$ with $ \Vert \boldsymbol{\theta } \Vert <C$,
$\mathrm{E
} \Vert \boldsymbol{\eta }_{i} \Vert <C$,
$\mathrm{E} ( \boldsymbol{\eta }_{i} ) =0$,
$\mathrm{E} ( \boldsymbol{\eta }_{i}
\boldsymbol{\eta }_{i}^{\prime } ) =$
$\boldsymbol{\Omega }_{\eta }$, and
$ \Vert \boldsymbol{\Omega }_{\eta } \Vert <C$. (b) Let
$
\boldsymbol{q}_{it}=\boldsymbol{w}_{it}\boldsymbol{w}_{it}^{\prime }
\boldsymbol{\eta }_{i}$, then
$\mathrm{E} ( \boldsymbol{q}_{it} ) =
\boldsymbol{q}_{i}$ (fixed),
$\sup_{i} \Vert \boldsymbol{q}
_{i} \Vert <C$,
$\sup_{i,t}\mathrm{E} \Vert \boldsymbol{q}
_{it} \Vert ^{2}<C$, and
$\sup_{i}\mathrm{E } \Vert \boldsymbol{w}
_{i,T+1}^{\prime }\boldsymbol{\eta }_{i} \Vert ^{2}<C$.
\end{assumption}
\begin{assumption}
\label{ass:5}
$\boldsymbol{\eta }_{i}$ is distributed independently of
$
\boldsymbol{\varepsilon }_{i}$, for all $i$.
\end{assumption}
\begin{assumption}
\label{ass:2.2}
$ \bar{\boldsymbol{\xi }}_{NT}=N^{-1}\sum_{i=1}^{N}
\boldsymbol{\xi }_{iT}=O_{p} ( N^{-1/2}T^{-1/2} ) $, where
$
\boldsymbol{\xi }_{iT }=T^{-1}\boldsymbol{W}_{i}^{\prime }
\boldsymbol{
\varepsilon }_{i}=T^{-1}\sum_{t=1}^{T}\boldsymbol{w}_{it}
\varepsilon _{it}=O_{p}(T^{-1/2})$.
\end{assumption}
\begin{assumption}
\label{ass:3b}
There exists a fixed $T_{0}$ such that for all $T>T_{0}$ and
$
N=1,2,\ldots $, the pooled covariance matrices
$\boldsymbol{\bar{Q}}_{NT}$ and $\boldsymbol{\bar{Q}}_{N}$, defined in
terms of
$\boldsymbol{Q}_{iT
}=T^{-1}\boldsymbol{W}_{i}^{\prime }\boldsymbol{W}_{i}$ and
$\boldsymbol{Q}
_{i}=\mathrm{E} ( \boldsymbol{Q}_{iT } ) $,
\begin{equation}
\boldsymbol{\bar{Q}}_{NT}=N^{-1}\sum
_{i=1}^{N}\boldsymbol{Q}_{iT}
, \quad \text{and}\quad \boldsymbol{\bar{Q}}_{N}=\mathrm{E} ( \boldsymbol{
\bar{Q}}_{NT} ) =N^{-1}\sum_{i=1}^{N}
\boldsymbol{Q}_{i}, \label{Qbar}
\end{equation}
are positive definite,
$ \Vert \boldsymbol{\bar{Q}}_{N}^{-1} \Vert <C$, and
\begin{equation*}
\sup_{N,T}\mathrm{E} \bigl[ \lambda _{\max }^{2}
( \boldsymbol{\bar{Q}}
_{NT} ) \bigr] <C<\infty
,\quad \text{and}\quad \sup_{N,T}\mathrm{E} \bigl[ \lambda
_{\max }^{2} \bigl( \boldsymbol{\bar{Q}}_{NT}^{-1}
\bigr) \bigr] <C<\infty .
\end{equation*}
\end{assumption}
\begin{assumption}
\label{ass:6}
$(\boldsymbol{\varepsilon }_{i},\boldsymbol{W}_{i},
\boldsymbol{
\eta }_{i})$ are distributed independently over $i$.
\end{assumption}
For pooled estimation of $\boldsymbol{\theta }$, the conditions on
$
\boldsymbol{Q}_{iT}$ can be relaxed and it is sufficient that
$\bar{
\boldsymbol{Q}}_{NT}$ is positive definite, and
$\sup_{N,T}\mathrm{E}
\Vert \boldsymbol{Q}_{NT}^{-1} \Vert ^{2}<C$. Assumptions
\ref{ass:4} and \ref{ass:5} identify the population mean of
$\boldsymbol{\theta }
_{i}$ denoted by $\boldsymbol{\theta }$, but allow for correlated heterogeneity.\footnote{We simplify the notation and use $\boldsymbol{\theta }$, rather than
$
\boldsymbol{\theta }_{0}$, to denote the population mean, which is technically
more appropriate.} The degree of parameter heterogeneity is measured by
the norm of $\boldsymbol{\Omega }_{\eta }$, and the extent to which heterogeneity
is correlated is measured by the norm of $\boldsymbol{q}
_{i}$.\footnote{Under Assumption~\ref{ass:weak_exogeneity_a},
$\mathrm{E} ( \boldsymbol{
\xi }_{iT} ) =T^{-1}\sum_{t=1}^{T}\mathrm{E} (
\boldsymbol{w}
_{it}\varepsilon _{it} ) =\boldsymbol{0}$, and
$\mathrm{E} ( \boldsymbol{\xi }_{NT} ) =\boldsymbol{0}$. Note
that $\varepsilon _{it}$ and $\boldsymbol{w}_{it}$ are uncorrelated but
not independently distributed. Under Assumption~\ref{ass:3},
$ \Vert \boldsymbol{\bar{Q}}
_{NT} \Vert \leq \sup_{i} \Vert \boldsymbol{Q}_{iT}
\Vert <C$, and
$ \Vert \boldsymbol{\bar{Q}}_{N} \Vert \leq \sup_{i}
\Vert \boldsymbol{Q}_{i} \Vert <C$.}
Assumptions \ref{ass:4}--\ref{ass:3b} are not required for forecasts based
on the individual estimates and the associated MSFE. Assumption~\ref{ass:6} of cross-sectional independence for $\varepsilon _{it}$ (or
$\boldsymbol{w}_{it}$) is not needed to establish results on the MSFE of individual
forecasts. However, we do require some degree of uncorrelatedness over
$i$ when the objective is to compute the MSFE averaged across all
$N$ units under consideration or over a subgroup of the units. In particular,
to ensure that the cross-sectional average MSFE tends to a nonrandom limit,
the units under consideration must satisfy the law of large numbers. To
this end, we need the units to be cross-sectionally weakly correlated,
possibly conditional on known (or estimated) common factors. The situation
is different when we consider pooled or Bayesian forecasts. Optimality
of these forecasts \textit{does} depend on the assumption of cross-sectional
independence, or at least some form of weak cross-sectional dependence.
A comprehensive analysis of the implications of cross-sectional dependence
for forecast combinations and comparisons of predictive accuracy are beyond
the scope of the present paper, however.\footnote{Cross-sectional dependence in forecast errors can be exploited by using
interactive time effects (latent factors) or spatial (network) effects;
see, for example, \cite{Chuetal2016}.}
We measure the degree of correlated heterogeneity for unit $i$ at time
$t$ by
$\boldsymbol{q}_{i}=\mathrm{E} ( \boldsymbol{w}_{it}
\boldsymbol{w}
_{it}^{\prime }\boldsymbol{\eta }_{i} ) $ and, on average, by
\begin{equation}
\boldsymbol{\bar{q}}_{NT}=N^{-1}T^{-1}\sum
_{i=1}^{N}\boldsymbol{W}
_{i}^{\prime }
\boldsymbol{W}_{i}\boldsymbol{\eta }_{i}=N^{-1}T^{-1}
\sum_{i=1}^{N}\sum
_{t=1}^{T}\boldsymbol{w}_{it}
\boldsymbol{w}_{it}^{
\prime }
\boldsymbol{\eta
}_{i}. \label{qbarNT}
\end{equation}
Taking expectations,
\begin{equation}
\mathrm{E} ( \boldsymbol{\bar{q}}_{NT} ) = \boldsymbol{
\bar{q}}
_{N}=N^{-1}\sum
_{i=1}^{N}\boldsymbol{q}_{i}.
\label{qbar}
\end{equation}
Assumptions \ref{ass:4} and \ref{ass:5} accommodate correlated heterogeneity
and allow for nonzero values of
$\mathrm{E} ( \boldsymbol{W}
_{i}^{\prime }\boldsymbol{W}_{i}\boldsymbol{\eta }_{i} ) $. In the
context of fixed effects models, the intercepts $\alpha _{i}$ in (\ref{eq:model1})
are allowed to have nonzero correlation with the regressors, but optimality
of forecasts based on pooled estimates of $\boldsymbol{\beta } $ requires
Assumption~\ref{ass:5} and the condition
$\lim_{n\rightarrow \infty }n^{-1}\*\sum_{i=1}^{n}\mathrm{E} (
\boldsymbol{X}_{i}^{\prime }
\boldsymbol{M}_{T}\boldsymbol{X}_{i}\boldsymbol{\eta }_{i\beta }
) =
\boldsymbol{0}$, where
$\boldsymbol{\eta }_{i\beta }=\boldsymbol{\beta }_{i}-
\boldsymbol{\beta }$,
$\boldsymbol{M}_{T}=\boldsymbol{I}_{T}-\boldsymbol{
\tau }_{T} ( \boldsymbol{\tau }_{T}^{\prime }
\boldsymbol{\tau }
_{T} ) ^{-1}\boldsymbol{\tau }_{T}^{\prime }$,
$\boldsymbol{\tau }_{T}$ is a $T\times 1$ vector of ones, and
$\boldsymbol{I}_{T}$ is a $T\times T$ identity matrix.\footnote{See \citet{Pesetal2024b}. Note that
$\mathrm{E} ( \boldsymbol{X}
_{i}^{\prime }\boldsymbol{M}_{T}\boldsymbol{X}_{i}\boldsymbol{\eta }_{i
\beta } ) =\boldsymbol{0}$ is sufficient but not necessary for the
validity of fixed effects estimation. This condition is not met if
$\boldsymbol{x}
_{it}$ includes lagged values of $y_{it}$, even if
$T\rightarrow \infty $.}
\section{Theoretical results on forecasting performance}
\label{sec:theory_msfe}
We next use the setup and assumptions from Section~\ref{SETUP} to establish
theoretical results on the forecasting performance of different modeling
approaches. Section~\ref{ForecastsIP} discusses forecasts based on individual
and pooled estimation, and building on this, Section~\ref{sec:FE_theory}
covers fixed effects forecasts.
Note that our theoretical framework can be equally applied to forecasts
across groups instead of individuals, when there are \textit{a priori} known groups
such as industries or states within a given country. Pooled regressions
can be applied to any given, \textit{a priori} known group, so long as the number
of units within the group is sufficiently large and the cross-sectional
dependence of units within the group is sufficiently weak. Failure of the
latter assumption implies that there are missing pervasive (strong) common
factors that must also be taken into account but such an extension lies beyond the scope
of the present paper.
\subsection{Forecasts based on individual and pooled estimation\label
{ForecastsIP}}
We are interested in forecasting $y_{i,T+1}$ conditional on the information
known at time $T$, which we denote by $\boldsymbol{w}_{i,T+1}$ to clarify
the correspondence to $y_{i,T+1}$. Without loss of generality, given the
conditional nature of the forecasting exercise, we assume that $
\sup_{i,T}\left\Vert \boldsymbol{w}_{i,T+1}\right\Vert <C$.\footnote{
See part (a) of Assumption \ref{ass:3}.} Forecasts based on individual
estimators take the form
\begin{equation}
\hat{y}_{i,T+1}=\hat{\boldsymbol{\theta }}_{i}^{\prime }\boldsymbol{w}
_{i,T+1},\text{ }i=1,2,\ldots ,N, \label{eq:indFore}
\end{equation}
where $\hat{\boldsymbol{\theta }}_{i}=(\boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}_{i})^{-1}\boldsymbol{W}_{i}^{\prime }\boldsymbol{y}_{i},$ is
the least squares estimator of $\boldsymbol{\theta }_{i}$. Similarly,
forecasts based on the pooled estimator are given by
\begin{equation}
\tilde{y}_{i,T+1}=\tilde{\boldsymbol{\theta }}^{\prime }\boldsymbol{w}
_{i,T+1},\text{ }i=1,2,\ldots ,N, \label{eq:poolFore}
\end{equation}
where $\tilde{\boldsymbol{\theta }}=(\boldsymbol{W}^{\prime }\boldsymbol{W}
)^{-1}\boldsymbol{W}^{\prime }\boldsymbol{y}$. Using (\ref{Qbar}), (\ref
{qbarNT}) and the definition of $\boldsymbol{\bar{\xi}}_{NT}$ in Assumption
\ref{ass:2.2},
\begin{equation}
\tilde{\boldsymbol{\theta }}-\boldsymbol{\theta }_{i}=-\boldsymbol{\eta }
_{i}+\boldsymbol{\bar{Q}}_{NT}^{-1}\boldsymbol{\bar{q}}_{NT}+\boldsymbol{
\bar{Q}}_{NT}^{-1}\boldsymbol{\bar{\xi}}_{NT}. \label{bpooled}
\end{equation}
Forecast errors from these schemes take the form
\begin{eqnarray}
\hat{e}_{i,T+1} &=&y_{iT+1}-\hat{y}_{i,T+1}=\varepsilon _{i,T+1}-(\hat{
\boldsymbol{\theta }_{i}}-\boldsymbol{\theta }_{i})^{\prime }\boldsymbol{w}
_{i,T+1}, \label{ehat} \\
\tilde{e}_{i,T+1} &=&y_{iT+1}-\tilde{y}_{i,T+1}=\varepsilon _{i,T+1}-(\tilde{
\boldsymbol{\theta }}-\boldsymbol{\theta }_{i})^{\prime }\boldsymbol{w}
_{i,T+1}. \label{etilda}
\end{eqnarray}
\subsubsection*{Forecasts based on individual estimation}
Noting that
$(\hat{\boldsymbol{\theta }}_{i}-\boldsymbol{\theta }
_{i})^{\prime }\boldsymbol{w}_{i,T+1}=\boldsymbol{\varepsilon }_{i}^{
\prime }
\boldsymbol{W}_{i}(\boldsymbol{W}_{i}^{\prime }\*\boldsymbol{W}_{i})^{-1}
\boldsymbol{w}_{i,T+1}$, it is easily seen that the forecasts based on
the individual estimates generate the following average MSFE:
\begin{equation}
N^{-1}\sum_{i=1}^{N}
\hat{e}_{i,T+1}^{2}=N^{-1}\sum
_{i=1}^{N} \varepsilon _{i,T+1}^{2}+T^{-1}S_{NT}-2R_{NT},
\label{MSFE_i}
\end{equation}
where $S_{NT}=N^{-1}\sum_{i=1}^{N}s_{iT}$,
$R_{NT}=N^{-1}\sum_{i=1}^{N}r_{iT
}$, with elements
\begin{equation}
r_{iT }= \bigl( \boldsymbol{\varepsilon }_{i}^{\prime }
\boldsymbol{W}_{i}\bigl(
\boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}_{i}\bigr)^{-1}\boldsymbol{w}
_{i,T+1}
\bigr) \varepsilon _{i,T+1}, \label{r_iT}
\end{equation}
and
\begin{equation}
s_{iT}=\boldsymbol{w}_{i,T+1}^{\prime }
\boldsymbol{Q}_{iT}^{-1} \bigl( T^{-1}
\boldsymbol{W}_{i}^{\prime }\boldsymbol{\varepsilon
}_{i} \boldsymbol{
\varepsilon }_{i}^{\prime }
\boldsymbol{W}_{i} \bigr) \boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}_{i,T+1}. \label{s_iT}
\end{equation}
Under Assumptions \ref{ass:1} and \ref{ass:3},
$\mathrm{E} ( r_{iT
} ) =0$ and
$\sup_{i,T}\mathrm{E} \vert r_{iT} \vert <C$, and under cross-sectional
independence (Assumption~\ref{ass:6}) we have
$
R_{NT}=O_{p}(N^{-1/2})$. Similarly,
$\sup_{i,T}\mathrm{E} \vert s_{iT} \vert <C$,
\begin{equation*}
\mathrm{E} ( s_{iT } ) =\mathrm{E} \biggl[ \boldsymbol{w}
_{i,T+1}^{\prime }
\boldsymbol{Q}_{iT}^{-1} \biggl( \frac{
\boldsymbol{W}
_{i}^{\prime }\boldsymbol{\varepsilon
}_{i}\boldsymbol{\varepsilon }
_{i}^{\prime }
\boldsymbol{W}_{i}}{T} \biggr) \boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}_{i,T+1} \biggr] ,
\end{equation*}
$S_{NT}=\mathrm{E} ( S_{NT} ) +O_{p}(N^{-1/2})$, and we obtain
the results summarized in the following proposition for the average MSFE
of the forecasts based on the individual estimates (for a detailed proof,
see Section~\ref{Proof_individual} of the Appendix):
\begin{proposition}
\label{prop:individual}
\begin{enumerate}
\item[(a)] Suppose that Assumptions \ref{ass:1}--\ref{ass:weak_exogeneity_b} and
\ref{ass:6} hold. Then, for a fixed $T_{0}$ such that $T>T_{0}$, the average
MSFE resulting from individual-specific estimation of the parameters, given
by (\ref{MSFE_i}), has the following representation:
\begin{equation}
N^{-1}\sum_{i=1}^{N}
\hat{e}_{i,T+1}^{2}=N^{-1}\sum
_{i=1}^{N} \varepsilon _{i,T+1}^{2}+T^{-1}h_{NT}+O_{p}
\bigl(N^{-1/2}\bigr), \label{MSFEI}
\end{equation}
where
\begin{equation}
h_{NT}=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \biggl[ \boldsymbol{w}_{i,T+1}^{
\prime }
\boldsymbol{Q}_{iT}^{-1} \biggl( \frac{
\boldsymbol{W}_{i}^{\prime }\boldsymbol{
\varepsilon
}_{i}\boldsymbol{\varepsilon }_{i}^{\prime }
\boldsymbol{W}_{i}}{T
} \biggr) \boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}_{i,T+1} \biggr] , \label{h_NT}
\end{equation}
$\boldsymbol{Q}_{iT}=T^{-1}\boldsymbol{W}_{i}^{\prime }\boldsymbol{W}_{i}$,
$
h_{NT}>0$, and $h_{NT}=O(1)$.\vadjust{\goodbreak}
\item[(b)] If $\boldsymbol{W}_{i}$ is strictly exogenous, $h_{NT}$ simplifies
to
$ h_{NT}=N^{-1}\sum_{i=1}^{N}\sigma _{i}^{2}\mathrm{E} (
\boldsymbol{w}
_{i,T+1}^{\prime }\boldsymbol{Q}_{iT}^{-1}\*\boldsymbol{w}_{i,T+1}
) $.
\end{enumerate}
\end{proposition}
The $h_{NT}$ term captures the cost associated with the error in estimation
of $\hat{\boldsymbol{\theta }}_{i}$. For typical panel data sets,
$T$ is not large and parameter estimation uncertainty captured by the
$O ( T^{-1} ) $ term $T^{-1}h_{NT}$ in (\ref{MSFEI}) can therefore
be important. Parameter heterogeneity, in contrast, does not affect the
accuracy of the forecasts in (\ref{MSFEI}). The magnitude of $h_{NT}$ plays
an important role in the comparisons of forecasts based on individual and
pooled estimates and depends on how far the predictors are from their mean.
For example, when $\boldsymbol{w}_{it}=(1,x_{it})^{\prime }$ and
$x_{it}$ is strictly exogenous,
\begin{equation*}
h_{NT}=\bar{\sigma}_{N}^{2}+N^{-1}
\sum_{i=1}^{N}\sigma
_{i}^{2} \text{E} \biggl[ \frac{ (
x_{i,T+1}-\bar{x}_{iT} ) ^{2}}{s_{iT}^{2}}
\biggr] ,
\end{equation*}
where $\bar{\sigma}_{N}^{2}=N^{-1}\sum_{i=1}^{N}\sigma _{i}^{2}$,
$
s_{iT}^{2}=T^{-1}\sum_{t=1}^{T}(x_{it}-\bar{x}_{iT})^{2}$, and
$\bar{x}
_{iT}=T^{-1}\sum_{t=1}^{T}x_{it}$. Hence, $h_{NT}$ is minimized when
$
x_{i,T+1}=\bar{x}_{iT}$, for all $i$. When
$x_{i,T+1}\neq \bar{x}_{iT}$ for most $i$, $T$ must be sufficiently large
such that $\sup_{i} \mathrm{E} [ ( x_{i,T+1}-\bar{x}_{iT} ) ^{2}/s_{iT}^{2}
] <C$.
\subsubsection*{Forecasts based on pooled estimation}
While the forecast accuracy results for the individual regressions do not
depend on the degree of parameter heterogeneity, whether correlated or
not, the degree of correlated heterogeneity does matter for consistency
of the pooled estimator. Using (\ref{bpooled}) in (\ref{etilda}), we can
express the squared forecast error when pooled estimates are used as follows:
\begin{equation*}
\tilde{e}_{i,T+1}^{2}=\varepsilon _{i,T+1}^{2}+
\boldsymbol{w}
_{i,T+1}^{\prime }\boldsymbol{d}_{i,NT}
\boldsymbol{d}_{i,NT}^{\prime }
\boldsymbol{w}_{i,T+1}-2
\boldsymbol{d}_{i,NT}^{\prime } \boldsymbol{w}
_{i,T+1}
\varepsilon _{i,T+1},
\end{equation*}
where
$\boldsymbol{d}_{i,NT}=-\boldsymbol{\eta }_{i}+\boldsymbol{\bar{Q}}
_{NT}^{-1}\boldsymbol{\bar{q}}_{NT}+\boldsymbol{\bar{Q}}_{NT}^{-1}
\boldsymbol{\bar{\xi}}_{NT}$, $\boldsymbol{\bar{Q}}_{NT}$, and
$\boldsymbol{
\bar{q}}_{NT}$ are defined by (\ref{Qbar}) and (\ref{qbarNT}), and
$
\boldsymbol{\bar{\xi}}_{NT}$ is defined under Assumption~\ref{ass:2.2}. After some algebra, and averaging over $i$, we have
\begin{equation}
N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2} = N^{-1}\sum
_{i=1}^{N} \varepsilon _{i,T+1}^{2}+N^{-1}
\sum_{i=1}^{N}\boldsymbol{w}_{i,T+1}^{
\prime }
\boldsymbol{
\eta }_{i}\boldsymbol{\eta
}_{i}^{\prime }\boldsymbol{w}_{i,T+1}
+\tilde{S}_{N,T+1}+2\tilde{R}_{N,T+1},\label{MSFEPool}
\end{equation}
where $\tilde{S}_{N,T+1}$, and $\tilde{R}_{N,T+1}$ are defined by equations
(
\ref{stilda}) and (\ref{rtilda}) in Section~\ref{Proof_pooled} of the Appendix.
It can be shown that $\tilde{R}_{N,T+1}=O_{p}(N^{-1/2})$, and
$
\tilde{S}_{N,T+1}=-\boldsymbol{\bar{q}}_{N}^{\prime }
\boldsymbol{\bar{Q}}
_{N}^{-1}\boldsymbol{\bar{q}}_{N}+O_{p} ( N^{-1/2} ) $.
The limiting properties of the average MSFE based on pooled estimates are summarized
in the following proposition.
\begin{proposition}
\label{prop:pooled}
\begin{enumerate}
\item[(a)] Under Assumptions \ref{ass:1}--\ref{ass:6}, the MSFE for the forecasts
based on pooled estimation of the parameters, given by (\ref{MSFEPool}),
is
\begin{equation}
N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2}=N^{-1}\sum
_{i=1}^{N} \varepsilon _{i,T+1}^{2}+
\Delta _{NT}+O_{p}\bigl(N^{-1/2}\bigr),
\label{MSFEP1}
\end{equation}
where
\begin{equation}
\Delta _{NT}=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl( \boldsymbol{w}
_{i,T+1}^{\prime }
\boldsymbol{\eta }_{i}\boldsymbol{\eta }_{i}^{
\prime }
\boldsymbol{w}_{i,T+1} \bigr) -\boldsymbol{\bar{q}}_{N}^{\prime }
\boldsymbol{
\bar{Q}}_{N}^{-1}\boldsymbol{
\bar{q}}_{N}. \label{Delta_NT}
\end{equation}
\item[(b)] Parameter heterogeneity (whether correlated or uncorrelated) increases
the MSFE of the forecasts based on the pooled estimator, namely
$\Delta _{NT}>0$.
\end{enumerate}\end{proposition}
Note that the impact on the MSFE from neglected heterogeneity,
$\Delta _{NT}$
, does not vanish even if both $N$ and $T\rightarrow \infty $, which is
similar to the finding by \citet{PesSmi1995} for heterogeneous dynamic
panels since heterogeneity is always correlated in dynamic panels.\footnote{This latter property is illustrated by a simple panel AR(1) model with
heterogeneous AR coefficients in Section~\ref{PanelAR} of the Appendix.
See also \citet{Pesetal2024a} where estimation of such models with
short $T$ panels is considered.}
\subsubsection*{A comparison of forecasts based on individual and pooled
estimates}
Next, we consider the difference in the average MSFE performance of the
forecasts based on the pooled versus individual parameter estimates. Proposition~\ref{prop:individual}
shows that the MSFE from the forecasts based on the individual estimates
will be affected by an estimation error term of the form
\begin{equation*}
h_{NT}=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \biggl[ \boldsymbol{w}_{i,T+1}^{
\prime }
\boldsymbol{Q}_{iT}^{-1} \biggl( \frac{
\boldsymbol{W}_{i}^{\prime }\boldsymbol{
\varepsilon
}_{i}\boldsymbol{\varepsilon }_{i}^{\prime }
\boldsymbol{W}_{i}}{T
} \biggr) \boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}_{i,T+1} \biggr] >0.
\end{equation*}
While the forecasts from the pooled estimates are more robust to estimation
errors, they are in turn affected by correlated and uncorrelated heterogeneity
as captured by the term
\begin{equation*}
\Delta _{NT}=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl( \boldsymbol{w}
_{i,T+1}^{\prime }
\boldsymbol{\eta }_{i}\boldsymbol{\eta }_{i}^{
\prime }
\boldsymbol{w}_{i,T+1} \bigr) -\boldsymbol{\bar{q}}_{N}^{\prime }
\boldsymbol{
\bar{Q}}_{N}^{-1}\boldsymbol{
\bar{q}}_{N}.
\end{equation*}
We compare the difference in the average MSFE of the forecasts from the
pooled versus individual estimates as a ratio measured relative to the
MSFE of the forecasts from the individual estimates:
\begin{equation*}
\frac{N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2}-N^{-1}\sum
_{i=1}^{N} \hat{e}_{i,T+1}^{2}}{N^{-1}
\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}}=
\frac{\Delta _{NT} -T^{-1}h_{NT}+O_{p}
\bigl(N^{-1/2}\bigr)}{N^{-1}\sum
_{i=1}^{N}\varepsilon _{i,T+1}^{2}
+T^{-1}h_{NT}+O_{p}\bigl(N^{-1/2}
\bigr)}.
\end{equation*}
Hence, there exists a $T_{0}$ such that, for a fixed $T>T_{0}$, and as
$N\rightarrow \infty $,
\begin{equation}
\frac{N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2}-N^{-1}\sum
_{i=1}^{N}\hat{e}
_{i,T+1}^{2}}{N^{-1}
\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}}
\overset{p}{
\rightarrow }\frac{\Delta -T^{-1}h_{T}}{
\bar{\sigma}^{2}+T^{-1}h_{T}},
\label{finiteT}
\end{equation}
where $h_{T}=\lim_{N\rightarrow \infty }h_{NT}\geq 0$,
$\Delta =\lim_{N\rightarrow \infty }\Delta _{N}\geq 0$, and
$\bar{\sigma}
^{2}=\lim_{N\rightarrow \infty }N^{-1}\sum_{i=1}^{N}\sigma _{i}^{2}>0$
\textbf{.} It follows that when $T$ is fixed and $N$ is large, the ranking
of the two forecasting schemes will depend on the sign and magnitude of
$\Delta -T^{-1}h_{T}$.\footnote{In comparing $\Delta _{T}$ with $T^{-1}h_{T}$, it is also important to
bear in mind that $h_{T}$ is well-defined if moments of
$\boldsymbol{\hat{\theta}}
_{i}$ (at least up to second order) exist (see the moment condition (\ref{Momentconthetai})).
This in turn requires that $T>T_{0}$ for some finite $
T_{0}$. The value of $T_{0}$ depends on the nature of the
$(\boldsymbol{w}
_{it},\varepsilon _{it}$) process and its distributional properties.}
For large values of $T$, however, the individual forecasts generate the
lowest MSFE values. Specifically, for a fixed $N$ and as
$T\rightarrow \infty $,
\begin{equation*}
\frac{N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2}-N^{-1}\sum
_{i=1}^{N}\hat{e}
_{i,T+1}^{2}}{N^{-1}
\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}}
\overset{p}{
\rightarrow }\frac{\Delta _{N}}{\bar{
\sigma}^{2}}+O_{p}\bigl(N^{-1/2}\bigr).
\end{equation*}
Similarly, when both $N$ and $T\rightarrow \infty $ (in any order)
\begin{equation*}
\frac{N^{-1}\sum_{i=1}^{N}
\tilde{e}_{i,T+1}^{2}-N^{-1}\sum
_{i=1}^{N}\hat{e}
_{i,T+1}^{2}}{N^{-1}
\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}}
\overset{p}{
\rightarrow }\Delta /\bar{\sigma}^{2}\geq 0,
\end{equation*}
where $\Delta =\lim_{T\rightarrow \infty }(\Delta _{T})$. Therefore, vanishing
estimation uncertainty implied by large $T$ means that, on average, individual
forecasts are at least as precise as pooled forecasts irrespective of
$N$.
\subsection{Forecasts based on fixed effects estimation\label{sec:FE_theory}}
The comparison of forecasts based on individual or pooled estimates can
be extended to intermediate cases where a subset of the parameters are
allowed to vary across units. A prominent example is the FE forecast
\begin{equation}
\hat{y}_{i,T+1}^{\text{FE}}=\hat{\alpha}_{i,\text{FE}}+
\boldsymbol{\hat{\beta
}}_{\text{FE}}^{\prime }
\boldsymbol{x}_{i,T+1}, \label{ForFE}
\end{equation}
where
$\hat{\alpha}_{i,\text{FE}}=\boldsymbol{\tau }_{T}^{\prime }(
\boldsymbol{y}_{i}-\boldsymbol{\hat{\beta}}_{\text{FE}}^{\prime }
\boldsymbol{
X}_{i})/T$ and
$\boldsymbol{\hat{\beta}}_{\text{FE}}= ( \sum_{i=1}^{N}
\boldsymbol{X}_{i}^{\prime }\boldsymbol{M}_{T}\boldsymbol{X}_{i}
) ^{-1}\sum_{i=1}^{N}\boldsymbol{X}_{i}^{\prime }
\boldsymbol{M}_{T}\boldsymbol{
y}_{i}$. The associated FE forecast error is given by
\begin{equation}
\hat{e}_{i,T+1}^{\text{FE}}=\bar{\bar{\varepsilon}}_{i,T+1}-(
\hat{
\boldsymbol{\beta }}_{\text{FE}}-\boldsymbol{\beta
}_{i})^{\prime } \bar{\bar{
\boldsymbol{x}}}_{i,T+1},
\label{eiFE}
\end{equation}
where
$\bar{\bar{\varepsilon}}_{i,T+1}=\varepsilon _{i,T+1}-
\bar{\varepsilon}
_{iT}$,
$\bar{\bar{\boldsymbol{x}}}_{i,T+1}=\boldsymbol{x}_{i,T+1}-
\boldsymbol{\bar{x}}_{iT}$, $\bar{\varepsilon}_{iT}$
$=T^{-1}\sum_{t=1}^{T}
\varepsilon _{it}$, and
$\bar{\boldsymbol{x}}_{iT}=T^{-1}\*\sum_{t=1}^{T}
\boldsymbol{x}_{it}$. Section S.5 in the Online Supplement provides details
of the derivation of the MSFE under fixed effects estimation:
\begin{equation}
N^{-1}\sum_{i=1}^{N}
\bigl( \hat{e}_{i,T+1}^{\text{FE}} \bigr) ^{2}=N^{-1}
\sum_{i=1}^{N}\bar{\bar{
\varepsilon}}_{i,T+1}^{2}+\Delta _{NT}^{
\text{FE}}-2c_{NT}^{\text{FE}}+O_{p}
\bigl(N^{-1/2}\bigr), \label{FEmsfe}
\end{equation}
where
\begin{equation}
\Delta _{NT}^{\text{FE}}=N^{-1}\sum
_{i=1}^{N}\mathrm{E}\bigl( \bar{\bar{
\boldsymbol{x}}}_{i,T+1}^{\prime }\boldsymbol{\eta
}_{i,\beta } \boldsymbol{
\eta }_{i,\beta }^{\prime }
\bar{\bar{\boldsymbol{x}}}_{i,T+1}\bigr)- \bar{
\boldsymbol{q}}_{N,\beta }^{\prime }\bar{\boldsymbol{Q}}_{N,\beta }^{-1}
\bar{
\boldsymbol{q}}_{N,\beta }, \label{eq:Delta_FE}
\end{equation}
$\boldsymbol{\eta }_{i,\beta }=\boldsymbol{\beta }_{i}-
\boldsymbol{\beta }$,
$\bar{\boldsymbol{\xi }}_{NT,\beta }=N^{-1}\sum_{i=1}^{N}T^{-1}
\boldsymbol{X}
_{i}^{\prime }\boldsymbol{M}_{T}\boldsymbol{\varepsilon }_{i}$, $\bar{
\boldsymbol{Q}}_{NT,\beta }=N^{-1}\sum_{i=1}^{N}T^{-1}
\boldsymbol{X}
_{i}^{\prime }\boldsymbol{M}_{T}\boldsymbol{X}_{i}$,
$\bar{\boldsymbol{q}}
_{NT,\beta }=N^{-1}\sum_{i=1}^{N} ( T^{-1}\boldsymbol{X}_{i}^{
\prime }
\boldsymbol{M}_{T}\boldsymbol{X}_{i} ) \boldsymbol{\eta }_{i,
\beta }$ and
\begin{equation}
c_{NT}^{\text{FE}}=-N^{-1}\sum
_{i=1}^{N}\mathrm{E} \bigl( \boldsymbol{\eta
}
_{i,\beta }^{\prime }\bar{\bar{\boldsymbol{x}}}_{i,T+1}
\bar{\varepsilon}
_{iT} \bigr) +\bar{\boldsymbol{q}}_{N,\beta }^{\prime }
\bar{\boldsymbol{Q}}
_{N,\beta }^{-1} \Biggl[
N^{-1}\sum_{i=1}^{N}
\mathrm{E} ( \bar{\boldsymbol{
x}}_{iT}\bar{
\varepsilon}_{iT} ) \Biggr] . \label{cNT-FE}
\end{equation}
$c_{NT}^{\text{FE}}$ tends to zero for $T$ sufficiently large or if
$
\boldsymbol{x}_{it}$ is strictly exogenous.
Similar to the case of the individual and pooled forecasts, for $T$ finite
and $N$ large, the ranking of the individual and FE forecasts will depend
on the relative magnitudes of estimation error and parameter heterogeneity.
Precise expressions can be found in the Supplemental Appendix. For
$T\rightarrow \infty $ the individual forecasts will be more precise than
the FE forecasts.
\section{Forecast combinations}\label{sec:combinations}
We next consider approaches that combine the forecasts from Section~\ref{sec:theory_msfe}
to minimize the MSFE.
\subsection{Combinations of individual and pooled forecasts}\label
{com_Ind_pool}
Given the MSFE trade-off associated with the forecasts in (\ref{eq:indFore})
and (\ref{eq:poolFore}), combining the forecasts based on the individual
and pooled estimates, $\hat{y}_{i,T+1}$ and $\tilde{y}_{i,T+1}$, may be
desirable. As noted in the literature (e.g., \citet{Tim2006}), forecast
combinations tend to perform particularly well, relative to the underlying
forecasts, if the forecast errors are weakly correlated and have MSFE values
of a similar magnitude. Correlations between forecast errors based on the
individual and pooled estimation schemes tend to be lower for (i) greater
differences in the estimates of $\boldsymbol{\theta }_{i}$ resulting from
larger estimation errors (small $T$); (ii) greater heterogeneity (large
$
\Vert \boldsymbol{\Omega }_{\eta } \Vert $), and (iii) greater
bias of the pooled estimator due to correlated heterogeneity.
If the level of parameter heterogeneity is either very large or very small,
one of the individual or pooled estimation approaches will be dominant,
reducing potential gains from forecast combination. Similarly, if
$T$ is very small but $N$ is large and there is little parameter heterogeneity,
we would expect pooled estimation to dominate individual estimation by
a sufficiently large margin that forecast combination offers small, if
any, gains. Conversely, if $T$ is very large, forecasts using individual
estimates will dominate forecasts based on pooled estimates by a sufficient
margin that renders forecast combination less attractive. Building on these
observations, we combine the two forecasts $\hat{y}_{i,T+1}$ and
$\tilde{y}
_{i,T+1}$ using common weights, $\omega $, to obtain\footnote{We focus here on a simple constant-coefficient linear combination scheme.
\cite{Lahetal2017} discuss a broader range of combination methods and
\citet{Ell2017} provides an analysis of the effect on the combination weights
and forecasting performance from having a large common component in the
forecast errors.}
\begin{equation}
y_{i,T+1}^{\ast }(\omega )=\omega \hat{y}_{i,T+1}+(1-
\omega ) \tilde{y}
_{i,T+1}, \label{combination}
\end{equation}
with associated forecast error
$e_{i,T+1}^{\ast }(\omega )=\omega \hat{e}
_{i,T+1}+(1-\omega )\tilde{e}_{i,T+1}$. The average MSFE of the combined
forecast is given by
\begin{eqnarray*}
N^{-1}\sum_{i=1}^{N}e_{i,T+1}^{\ast 2}(
\omega ) &=&\omega ^{2} \Biggl( N^{-1}\sum
_{i=1}^{N}\hat{e}_{i,T+1}^{2}
\Biggr) +(1-\omega )^{2} \Biggl( N^{-1}\sum
_{i=1}^{N}\tilde{e}_{i,T+1}^{2}
\Biggr)
\\
&&{}+2\omega (1-\omega ) \Biggl( N^{-1}\sum
_{i=1}^{N}\hat{e}_{i,T+1}
\tilde{e}
_{i,T+1} \Biggr) .
\end{eqnarray*}
The value of $\omega $ that minimizes the average MSFE is therefore given
by
\begin{equation}
\omega _{NT}^{\ast }= \frac{N^{-1}\sum
_{i=1}^{N}\tilde{e}_{i,T+1}^{2}-
\Biggl( N^{-1}\sum_{i=1}^{N}
\hat{e}_{i,T+1}\tilde{e}_{i,T+1} \Biggr) }{ \Biggl(
N^{-1}\sum_{i=1}^{N}
\hat{e}_{i,T+1}^{2} \Biggr) + \Biggl( N^{-1}\sum
_{i=1}^{N}
\tilde{e}_{i,T+1}^{2} \Biggr) -2 \Biggl( N^{-1}
\sum_{i=1}^{N}\hat{e}_{i,T+1}
\tilde{e}_{i,T+1} \Biggr) }. \label{wstar}
\end{equation}
Expressions for $N^{-1}\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}$ and
$
N^{-1}\sum_{i=1}^{N}\tilde{e}_{i,T+1}^{2}$ are given by (\ref{MSFEI})
and (
\ref{MSFEP1}), respectively. We obtain a similar expression for
$
N^{-1}\sum_{i=1}^{N}\hat{e}_{i,T+1}\tilde{e}_{i,T+1}$, with
$
N^{-1}\sum_{i=1}^{N}\varepsilon _{i,T+1}^{2}$ canceling out from
$\omega _{NT}^{\ast }$. The result is summarized in the following proposition
(proven in Appendix Section~\ref{ProofCombined}).
\begin{proposition}
\label{prop:combined}
\begin{enumerate}
\item[(a)] Under Assumptions \ref{ass:1}--\ref{ass:6}, and for a given value of
$\boldsymbol{w}_{i,T+1}$, the optimal combination weight that minimizes
the MSFE of the forecast combination in (\ref{combination}) is given by
\begin{equation}
\omega _{NT}^{\ast }= \frac{\Delta _{NT}-T^{-1}{
\psi }_{NT}}{\Delta _{NT}+T^{-1}h_{NT}-2T^{-1}{
\psi }_{NT}}+O_{p}\bigl(N^{-1/2}\bigr),
\label{w*NT}
\end{equation}
where
\begin{align}
h_{NT}&=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \biggl[ \boldsymbol{w}_{i,T+1}^{
\prime }
\boldsymbol{Q}_{iT}^{-1} \biggl( \frac{
\boldsymbol{W}_{i}^{\prime }\boldsymbol{
\varepsilon
}_{i}\boldsymbol{\varepsilon }_{i}^{\prime }
\boldsymbol{W}_{i}}{T
} \biggr) \boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}_{i,T+1} \biggr] >0, \label{hNT}
\\
\Delta _{NT}&=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl( \boldsymbol{w}
_{i,T+1}^{\prime }
\boldsymbol{\eta }_{i}\boldsymbol{\eta }_{i}^{
\prime }
\boldsymbol{w}_{i,T+1} \bigr) -\boldsymbol{\bar{q}}_{N}^{\prime }
\boldsymbol{
\bar{Q}}_{N}^{-1}\boldsymbol{
\bar{q}}_{N}>0, \label{DeltaNT}
\end{align}
and
\begin{align}
\psi _{NT} ={}&TN^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ \boldsymbol{\varepsilon }
_{i}^{\prime }
\boldsymbol{W}_{i}\bigl(\boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}
_{i}\bigr)^{-1}
\boldsymbol{w}_{i,T+1}\boldsymbol{w}_{i,T+1}^{\prime }
\bigr] \boldsymbol{\bar{Q}}_{N}^{-1}\boldsymbol{
\bar{q}}_{N} \notag
\\
& {}-TN^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ \boldsymbol{\varepsilon }
_{i}^{\prime }
\boldsymbol{W}_{i}\bigl(\boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}
_{i}\bigr)^{-1}
\boldsymbol{w}_{i,T+1}\boldsymbol{w}_{i,T+1}^{\prime }
\boldsymbol{
\eta }_{i} \bigr] .\label{epsi_NT}
\end{align}
\item[(b)] Under strict exogeneity, irrespective of whether heterogeneity is correlated,
we have $\psi _{NT}=0$,
$h_{NT}=N^{-1}
\sum_{i=1}^{N}\sigma ^{2}_{i}\mathrm{E} ( \boldsymbol{w}_{i,T+1}^{
\prime }\boldsymbol{Q}
_{iT}^{-1}\boldsymbol{w}_{i,T+1} ) $, and
\[
\Delta _{NT}=N^{-1}\sum_{i=1}^{N}\mathrm{E} \bigl( \boldsymbol{w}_{i,T+1}^{
\prime }
\boldsymbol{\Omega }_{\eta }\boldsymbol{w}_{i,T+1} \bigr) .
\]
\end{enumerate}\end{proposition}
For small to moderate values of $T$ and large $N $, we expect
$\omega _{NT}^{\ast }<1$, with a nonzero weight placed on forecasts
based on the pooled estimate.
\subsubsection*{Forecast combinations with individual weights}
\cite{Pesetal2022} show that, under strict exogeneity of the regressors
and uncorrelated heterogeneity, optimal weights can be obtained that are
specific to the individual unit. The combination of individual and pooled
forecast is then
\begin{equation*}
y_{i,T+1}^{\ast }=\omega _{i}
\hat{y}_{i,T+1}+(1-\omega _{i}) \tilde{y}
_{i,T+1},
\end{equation*}
where the optimal value of $\omega _{i}$ is given by
\begin{equation}
\omega _{i}^{\ast }= \frac{\boldsymbol{w}_{i,T+1}^{\prime }
\boldsymbol{\Omega }_{\eta }\boldsymbol{w}_{i,T+1}}{
\boldsymbol{w}_{i,T+1}^{\prime }\bigl(T^{-1}\sigma
_{i}^{2}\boldsymbol{Q}_{iT}^{-1}+
\boldsymbol{\Omega }_{\eta }\bigr)
\boldsymbol{w}_{i,T+1}}.
\label{eq:weights_per_i}
\end{equation}
The weights again depend on the variances and covariances of the underlying
forecast errors. Related to this, \cite{Giaetal2023} develop a random
effects approach for linear panels that similarly combines univariate and
pooled forecasts in a way that minimizes minimax-regret and MSFE.
\subsection{Combining individual and fixed effect forecasts}
\label{com_Ind_FE}
Combination weights can also be determined for the case where the pooled
forecast is replaced with the FE forecasts. In this case, the combined
forecast is given by
\begin{equation}
y_{i,T+1}^{\ast }(\omega _{\text{FE}})=\omega
_{\text{FE}}\hat{y}
_{i,T+1}+(1-\omega _{\text{FE}})
\hat{y}_{i,T+1,\text{FE}}, \label{eq:combinationFE}
\end{equation}
yielding the optimal pooled weight
\begin{equation}
\omega _{\text{FE},NT}^{\ast }= \frac{N^{-1}\sum
_{i=1}^{N} \bigl( \hat{e}
_{i,T+1}^{\text{FE}}
\bigr) ^{2}- \Biggl( N^{-1}\sum
_{i=1}^{N}\hat{e}_{i,T+1}^{
\text{FE}}
\hat{e}_{i,T+1} \Biggr) }{ \Biggl( N^{-1}\sum
_{i=1}^{N}\hat{e}
_{i,T+1}^{2}
\Biggr) +N^{-1}\sum_{i=1}^{N}
\bigl( \hat{e}_{i,T+1}^{\text{FE}
} \bigr) ^{2}-2
\Biggl( N^{-1}\sum_{i=1}^{N}
\hat{e}_{i,T+1}^{\text{FE}}\hat{e}
_{i,T+1}
\Biggr) }. \label{wNT-FE}
\end{equation}
The expressions for
$N^{-1}\sum_{i=1}^{N} ( \hat{e}_{i,T+1}^{\text{FE}
} ) ^{2}$ and $N^{-1}\sum_{i=1}^{N}\hat{e}_{i,T+1}^{2}$ are given
by (
\ref{FEmsfe}) and (\ref{MSFEI}), respectively, and the expression for
$
N^{-1}\sum_{i=1}^{N}\hat{e}_{i,T+1}^{\text{FE}}\hat{e}_{i,T+1}$ can be
similarly obtained. In this case, the shared term
$\sum_{i=1}^{N}(
\varepsilon _{i,T+1}-\bar{\varepsilon}_{iT})^{2}/N$ cancels out and we
have the result summarized in the following proposition with proofs provided
in Section~\ref{ProofCombinedFE} of the Appendix.
\begin{proposition}
\label{prop:combinedFE}
\begin{enumerate}
\item[(a)] Under Assumptions \ref{ass:1}--\ref{ass:6}, the optimal combination
weight that minimizes the MSFE of the forecast combination in (\ref{eq:combinationFE})
is given by
\begin{equation}
\omega _{\text{FE},NT}^{\ast }= \frac{\Delta _{NT}^{\text{FE}}-T^{-1}
\psi _{NT}^{\text{FE}}- \bigl( c_{NT}^{\text{FE}}-c_{NT,\beta }
\bigr) }{\Delta _{NT}^{\text{FE}}+T^{-1}h_{NT,\beta }-2T^{-1}
\psi _{NT}^{\text{FE}}}
+O_{p}
\bigl(N^{-1/2}\bigr), \label{wFE}
\end{equation}
where $\Delta _{NT}^{\text{FE}}$ and $c_{NT}^{\text{FE}}$ are defined in (
\ref{eq:Delta_FE}) and (\ref{cNT-FE}), respectively,
\begin{eqnarray*}
h_{NT,\beta }&=&N^{-1}\sum_{i=1}^{N}
\mathrm{E} \biggl[ \bar{\bar{\boldsymbol{x}}}
_{i,T+1}^{\prime }
\boldsymbol{Q}_{iT,\beta }^{-1} \biggl( \frac{
\boldsymbol{X}
_{i}^{\prime }\boldsymbol{M}_{T}
\boldsymbol{\varepsilon }_{i}\boldsymbol{
\varepsilon
}_{i}^{\prime }\boldsymbol{M}_{T}
\boldsymbol{X}_{i}}{T} \biggr) \boldsymbol{Q}_{iT,\beta }^{-1}
\bar{\bar{\boldsymbol{x}}}_{i,T+1} \biggr] , \label{eq:hNT_FEapp2}
\\
\psi _{NT}^{\text{FE}} &=&TN^{-1}\sum
_{i=1}^{N}\mathrm{E} \bigl[ (\hat{
\boldsymbol{\beta}}_{\text{FE}}-\boldsymbol{\beta}_{i}
)^{
\prime }\bar{\bar{
\boldsymbol{x}}}_{i,T+1}\bar{\bar{
\boldsymbol{x}}}_{i,T+1}^{\prime } \bigr] \bar{
\boldsymbol{Q}}_{N,\beta }^{-1}\bar{\boldsymbol{q}}_{N,
\beta }
\notag
\\
&&{}-TN^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ ( \hat{\boldsymbol{\beta}}_{
\text{FE}}-\boldsymbol{
\beta}_{i} )^{\prime } \bar{\bar{\boldsymbol{x}}}
_{i,T+1}
\bar{\bar{\boldsymbol{x}}}_{i,T+1}^{\prime } \boldsymbol{\eta
}
_{i,\beta } \bigr], \label{FEepsi}
\end{eqnarray*}
and
\begin{equation*}
c_{NT,\beta }=N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ \bar{\bar{\boldsymbol{x}}}
_{i,T+1}^{\prime }
\bigl(\boldsymbol{X}_{i}^{\prime }\boldsymbol{M}_{T}
\boldsymbol{
X}_{i}\bigr)^{-1}
\boldsymbol{X}_{i}^{\prime }\boldsymbol{M}_{T}
\boldsymbol{
\varepsilon }_{i}\bar{\varepsilon}_{iT}
\bigr] . \label{cNTbeta}
\end{equation*}
\item[(b)] Under uncorrelated heterogeneity, $\psi _{NT}^{\text{FE}}=0$, and
$
\Delta _{NT}^{\text{FE}}$ and $h_{NT,\beta}$ will be affected accordingly.
\end{enumerate}
\end{proposition}
\subsection{Estimation of combination weights}
\label{comW}
Estimates of the weights for the forecast combination in Proposition~\ref{prop:combined}
require estimates of $\Delta _{NT}$, $h_{NT}$, and $\psi _{NT}$. Under
Assumption~\ref{ass:6}, these terms can be estimated by their sample means
with unknown parameters replaced by their estimates. We summarize the estimators
here with details in Appendix~\ref{app:estimationOfWeights}
. Using (\ref{hNT}) and (\ref{DeltaNT}), the estimators of
$\Delta _{NT}$ and $h_{NT}$ are given by
\begin{equation}
\hat{\Delta}_{NT}=N^{-1}\sum
_{i=1}^{N}\boldsymbol{w}_{i,T+1}^{
\prime }
\tilde{
\boldsymbol{\eta }}_{i}\tilde{\boldsymbol{\eta
}}_{i}^{\prime } \boldsymbol{w}
_{i,T+1},
\label{eq:Deltahat}
\end{equation}
where
$\tilde{\boldsymbol{\eta }}_{i}=\tilde{\boldsymbol{\theta }}-
\boldsymbol{\hat{\theta}}_{i}$, and
\begin{equation}
\hat{h}_{NT}=N^{-1}\sum_{i=1}^{N}
\boldsymbol{w}_{i,T+1}^{\prime } \boldsymbol{
Q}_{iT}^{-1}
\hat{\boldsymbol{H}}_{iT}\boldsymbol{Q}_{iT}^{-1}
\boldsymbol{w}
_{i,T+1}, \label{eq:hhat}
\end{equation}
where
$\hat{\boldsymbol{H}}_{iT}=\hat{\sigma}_{i}^{2}T^{-1}\sum_{t=1}^{T}
\boldsymbol{w}_{it}\boldsymbol{w}_{it}^{\prime }$,
$\hat{\sigma}_{i}^{2}=$
$
\sum_{t=1}^{T}\hat{\varepsilon}_{it}^{2}/(T-K)$, and
$\hat{\varepsilon}
_{it}=y_{it}-\hat{\boldsymbol{\theta }}_{i}^{\prime }\boldsymbol{w}_{it}$.
We show in Appendix~\ref{app:estimationOfWeights} that
\begin{eqnarray*}
\hat{\Delta}_{NT}-\Delta _{NT} &=&O_{p}
\bigl( N^{-1/2} \bigr) +O_{p}\bigl(T^{-1}\bigr),
\\
\hat{h}_{NT}-h_{NT} &=&O_{p}
\bigl(N^{-1/2}\bigr)+O_{p} \biggl( \frac{\ln (N)}{
\sqrt{T}}
\biggr) .
\end{eqnarray*}
In the case of strictly exogenous regressors, $\hat{\Delta}_{NT}$ is a
consistent estimator of $\Delta _{NT}$ for fixed $T$ as
$N\rightarrow \infty $.
Consider now $\psi _{NT}$, given by (\ref{epsi_NT}), and recall that
$\psi _{NT}=0$ under uncorrelated heterogeneity. To estimate
$\psi _{NT}$ under correlated heterogeneity, we first note that the approach
of replacing expectations by sample moments and then estimating
$\boldsymbol{\varepsilon}_{i}$ from the OLS residuals,
$\boldsymbol{\hat{\varepsilon}}_{i}= (\mathbf{y}_{i}-\mathbf{W}_{i}
\hat{\boldsymbol{\theta }}_{i} ) $ will not work in the case of
$\psi _{NT}$, since
$\mathbf{W}_{i}^{\prime} \hat{\boldsymbol{\varepsilon}}_{i} = \mathbf{W}_{i}^{
\prime } (\mathbf{y}_{i} -\mathbf{W}_{i}
\hat{\boldsymbol{\theta }}_{i} ) =\mathbf{0}$ for all $i$. If used
in (\ref{epsi_NT}), this results in $\hat{\psi}_{NT}=0$, which is not a
consistent estimator of $\psi _{NT}$ under correlated heterogeneity. To
overcome this problem, we replace
$\boldsymbol{\varepsilon }_{i}^{\prime} \boldsymbol{W}_{i}(
\boldsymbol{W}_{i}^{\prime }\boldsymbol{W}_{i})^{-1}$ by
$ ( \hat{\boldsymbol{\theta }}_{i}-\boldsymbol{\theta }_{i}
)^{\prime}$ and note that $\psi _{NT}$ can be written equivalently
as (noting that $\boldsymbol{w}_{i,T+1} $ for $i=1,2,\ldots ,N$ are given)
\begin{align}
\psi _{NT} ={}&N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ T ( \hat{\boldsymbol{
\theta }}_{i}-
\boldsymbol{\theta }_{i} ) ^{\prime } \bigr]
\mathrm{E}
\bigl( \boldsymbol{w}_{i,T+1}\boldsymbol{w}_{i,T+1}^{\prime }
\bigr) \boldsymbol{\bar{Q}}_{N}^{-1}\boldsymbol{
\bar{q}}_{N} \notag
\\
&{} -N^{-1}\sum_{i=1}^{N}
\mathrm{E} \bigl[ T ( \hat{\boldsymbol{\theta }}
_{i}-
\boldsymbol{\theta }_{i} ) ^{\prime } \bigr] \mathrm{E} \bigl(
\boldsymbol{w}_{i,T+1}\boldsymbol{w}_{i,T+1}^{\prime }
\boldsymbol{\eta }
_{i} \bigr) .\label{epsi_NT2}
\end{align}
We now employ a half-jackknife estimator of
$\boldsymbol{\theta }_{i}$ (\citet{DhaJoc2015}, \cite{Chuetal2018})
to estimate
$T\mathrm{E}
( \hat{\boldsymbol{\theta }}_{i}-\boldsymbol{\theta }_{i}
) $, which is the small sample bias of
$\hat{\boldsymbol{\theta }}_{i}$ in the case of weakly exogenous regressors.
The half-jackknife estimator of $
\boldsymbol{\theta }_{i}$ is defined by
$\hat{\boldsymbol{\theta }}_{i,\mathit{JK}}=2
\hat{\boldsymbol{\theta }}_{i}-\frac{1}{2} (
\hat{\boldsymbol{\theta }}
_{ia}+\hat{\boldsymbol{\theta }}_{ib} ) $, where
$\hat{\boldsymbol{
\theta }}_{ia}$ and $\hat{\boldsymbol{\theta }}_{ib}$ are least squares
estimators of $\boldsymbol{\theta }_{i}$ based on two equal halves of the
sample of size $T_{h}=T/2$ (omitting an observation in the case of uneven
$T$
), namely
$\hat{\boldsymbol{\theta }}_{ia}= ( \sum_{t=1}^{T_{h}}
\boldsymbol{w}_{it}\boldsymbol{w}_{it}^{\prime } ) ^{-1}\sum_{t=1}^{T_{h}}
\boldsymbol{w}_{it}y_{it}$ and
$\hat{\boldsymbol{
\theta }}_{ib}= ( \sum_{t=T_{h}+1}^{T}\boldsymbol{w}_{it}
\boldsymbol{w}
_{it}^{\prime } ) ^{-1}\sum_{t=T_{h}+1}^{T}\boldsymbol{w}_{it}y_{it}$.
Then
$\mathrm{E} [ T ( \hat{\boldsymbol{\theta }}_{i}-
\boldsymbol{
\theta }_{i} ) ] $ can be estimated by
$T [ \frac{1}{2} ( \hat{\boldsymbol{\theta }}_{ia}+
\hat{\boldsymbol{\theta }}_{ib} ) -\hat{
\boldsymbol{\theta }}_{i} ] $ and $\psi _{NT}$ by
\begin{eqnarray}
{\hat{\psi}_{NT}} &=& \Biggl[ TN^{-1}\sum
_{i=1}^{N} \biggl[ \frac{1}{2} ( \hat{
\boldsymbol{\theta }}_{ia}+ \hat{\boldsymbol{\theta
}}_{ib} ) -\hat{
\boldsymbol{\theta }}_{i}
\biggr] ^{\prime }\boldsymbol{w}_{i,T+1} \boldsymbol{
w}_{i,T+1}^{\prime }
\Biggr] \boldsymbol{\bar{Q}}_{NT}^{-1} \boldsymbol{
\bar{q}
}_{NT} ( \hat{\boldsymbol{\eta }} )\notag
\\
&&{}-TN^{-1}\sum_{i=1}^{N}
\biggl[ \frac{1}{2} ( \hat{\boldsymbol{\theta }}
_{ia}+
\hat{\boldsymbol{\theta }}_{ib} ) - \hat{\boldsymbol{\theta
}}_{i}
\biggr] ^{\prime }\boldsymbol{w}_{i,T+1}
\boldsymbol{w}_{i,T+1}^{
\prime }\hat{
\boldsymbol{\eta
}}_{i},\label{eq:psihat}
\end{eqnarray}
where
$\hat{\boldsymbol{\eta }}_{i}=\hat{\boldsymbol{\theta }}
_{i}-N^{-1}\sum_{i=1}^{N}\hat{\boldsymbol{\theta }}_{i}$, and
$\boldsymbol{
\bar{q}}_{NT} ( \hat{\boldsymbol{\eta }} ) =(NT)^{-1}
\sum_{i=1}^{N}\boldsymbol{W}_{i}^{\prime }\boldsymbol{W}_{i}
\hat{
\boldsymbol{\eta }}_{i}$. Consistency of ${\hat{\psi}_{NT}}$ as an estimator
of $\psi _{NT}$ is established as $N,T\rightarrow \infty $, since
$T\mathrm{E
} ( \hat{\boldsymbol{\theta }}_{i,\mathit{JK}}-
\boldsymbol{\theta }_{i} ) =O ( T^{-1} ) $.\footnote{Note that ${\hat{\psi}_{NT}-\psi _{NT}}$ depends on
$\mathrm{E} [ ( \hat{\boldsymbol{\theta }}_{i}-
\boldsymbol{\theta }
_{i} ) - ( \frac{1}{2} ( \hat{\boldsymbol{\theta }}_{ia}+
\hat{
\boldsymbol{\theta }}_{ib} ) -\hat{\boldsymbol{\theta }}_{i}
) ] =\mathrm{E} ( \hat{\boldsymbol{\theta }}_{i,
\mathit{JK}}-\boldsymbol{
\theta }_{i} )$.} Thus, to use the half-jackknife method for models
with weakly exogenous regressors we need $T$ large, although it is not
required that $\sqrt{T}/N$ tends to zero, as it is in the case of large
$N$ and $T$ asymptotics.
The components of the weights in Proposition~\ref{prop:combinedFE} that
combine individual and fixed effects forecasts can be estimated in a similar
fashion, with details provided in Appendix~\ref{app:estimationOfWeights}.
\subsection{Empirical Bayes forecasts}
Bayesian panel forecasts are becoming increasingly common in empirical
applications and constitute an alternative approach to the frequentist
forecasts discussed so far. Due to their resemblance to our forecast combination
schemes and their recent popularity (e.g., \cite{Armetal2022} and
\citet{Efr2016}), we focus on empirical Bayes (EB) methods. The EB forecast
uses the estimator of \cite{Hsietal1999} and takes the form
$ \hat{y}_{i,T+1}^{\mathit{EB}}=\hat{\boldsymbol{\theta }}_{i,
\mathit{EB}
}^{\prime }\boldsymbol{w}_{i,T+1}$, where
\begin{equation}
\hat{\boldsymbol{\theta }}_{i,\mathit{EB}}=\bigl(\hat{\sigma}_{i}^{-2}
\boldsymbol{
W}_{i}^{\prime }\boldsymbol{W}_{i}+
\boldsymbol{\hat{\Omega}}_{\eta }^{-1}\bigr)^{-1}
\bigl( \hat{\sigma}_{i}^{-2}\boldsymbol{W}_{i}^{\prime }
\boldsymbol{y}
_{i}+\boldsymbol{\hat{\Omega}}_{\eta }^{-1}
\bar{\hat{\boldsymbol{
\theta }}}\bigr), \label{EB}
\end{equation}
$
\bar{\hat{\boldsymbol{\theta }}}=N^{-1}\sum_{i=1}^{N}
\hat{\boldsymbol{\theta }}_{i}$, $ \hat{\sigma}_{i}^{2}=(T-K)^{-1}
\hat{\boldsymbol{\varepsilon }}
_{i}^{\prime }\hat{\boldsymbol{\varepsilon }}_{i}$, and
$\boldsymbol{\hat{\Omega}}_{\eta }=\frac{1}{N}\sum_{i=1}^{N}(
\hat{
\boldsymbol{\theta }}_{i}-\bar{\hat{\boldsymbol{\theta }}})(
\hat{\boldsymbol{
\theta }}_{i}-\bar{\hat{\boldsymbol{\theta }}})^{\prime }$, where
$\hat{
\boldsymbol{\varepsilon }}_{i}=\boldsymbol{y}_{i}-\boldsymbol{W}_{i}
\hat{
\boldsymbol{\theta }}_{i}$, and
$\hat{\boldsymbol{\theta }}_{i}= ( \boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}_{i} ) ^{-1}\boldsymbol{W}
_{i}^{\prime }\boldsymbol{y}_{i}$.\footnote{It is necessary that $N>T$ for $\boldsymbol{\hat{\Omega}}_{\eta }$ to be
positive definite.} $\hat{\boldsymbol{\theta }}_{i,\mathit{EB}}$ can also
be written as a weighted average of
$\hat{\boldsymbol{\theta }}_{i}$, which allows for full heterogeneity,
and the mean group estimator, $\bar{\hat{
\boldsymbol{\theta }}}$, namely
$\hat{\boldsymbol{\theta }}_{i,\mathit{EB}}=
\boldsymbol{\mathcal{W}}_{iT}\hat{\boldsymbol{\theta }}_{i}+ (
\boldsymbol{I}_{k}-\boldsymbol{\mathcal{W}}_{iT} )
\bar{\hat{
\boldsymbol{\theta }}}$, with the weight matrix
$\boldsymbol{\mathcal{W}}
_{iT}$ given by
\begin{equation}
\boldsymbol{\mathcal{W}}_{iT}= \bigl( \boldsymbol{I}_{k}+T^{-1}
\hat{\sigma}
_{i}^{2}\boldsymbol{Q}_{iT }^{-1}
\boldsymbol{\hat{\Omega}}_{\eta }^{-1} \bigr)
^{-1}, \label{eq:empBayweights}
\end{equation}
recalling that
$\boldsymbol{Q}_{iT }=T^{-1}\boldsymbol{W}_{i}^{\prime }
\boldsymbol{W}_{i}$ is invertible under Assumption~\ref{ass:3}. The weights
on the heterogeneous estimates are larger, the greater the degree of heterogeneity,
as measured by the norm of $\boldsymbol{\hat{\Omega}}_{\eta }$
, with
$\hat{\boldsymbol{\theta }}_{i,\mathit{EB}}\rightarrow
\hat{
\boldsymbol{\theta }}_{i}$ as
$ \Vert \boldsymbol{\hat{\Omega}}_{\eta } \Vert
\rightarrow \infty $. Also, since
$\hat{\sigma}_{i}^{2}
\boldsymbol{Q}_{iT }^{-1}\boldsymbol{\hat{\Omega}}_{\eta }^{-1}$ is bounded
in $T$, $\hat{\boldsymbol{\theta }}_{i,\mathit{EB}}$ converges
\textit{
numerically} to $\hat{\boldsymbol{\theta }}_{i}$, as
$T\rightarrow \infty $. Hence, one would expect the EB estimator to perform
well even when $T$ is relatively small and the degree of heterogeneity
is not too large. For large $T$, EB and individual forecasts coincide and
both methods will work well.
The EB weights do \textit{not} depend on $\boldsymbol{w}_{i,T+1}$ and are
derived assuming uncorrelated heterogeneity and strictly exogenous regressors.
They have the desirable feature of placing more weights on individual estimates
if they are precisely estimated relative to the degree of parameter heterogeneity
measured by $\boldsymbol{\hat{\Omega}}_{\eta }$. The individual optimum
weights in (\ref{eq:weights_per_i}) fall somewhere between the common optimal
weights and the EB weights.\footnote{While the EB estimator in (\ref{EB}) is fully parametric, other studies
pursue a nonparametric approach to the distribution of
$\hat{\boldsymbol{
\theta }}_{i}$; see, for example, \citet{BroGre2009} and \citet{GuKoe2017}, and more recently, \citet{Liu2023}
and \cite{Liuetal2023}.} Like
the EB weights, consistent estimation of individual weights require strict
exogeneity and uncorrelated heterogeneity.
The EB weights are comparable to the unit-specific weights given by (\ref{eq:weights_per_i})
and the two sets of weights coincide only when
$\mathbf{w
}_{i,T+1}$ is an scalar. To see this, note that the estimates of the unit
specific weights can be written as
\begin{equation}
\hat{\omega}_{iT}^{\ast }= \biggl[ 1+\hat{
\sigma}_{i}^{2}T^{-1} \frac{
\mathbf{w}
_{i,T+1}^{\prime }\mathbf{Q}_{iT}^{-1}
\mathbf{w}_{i,T+1}}{\mathbf{w}
_{i,T+1}^{\prime }
\hat{\boldsymbol{\Omega}}_{\eta }\mathbf{w}_{i,T+1}}
\biggr] ^{-1}, \label{eq:indOptimalWeights}
\end{equation}
and reduces to the EB weights only when $K=1$, and
$\hat{\omega}_{iT}^{\ast }$ no longer depend on $\mathbf{w}_{i,T+1}$. But
in general the estimates of the unit-specific weights differ from the EB
weights.
An alternative to the EB forecast is a hierarchical Bayesian approach as
proposed by \citet{LinSmi1972} and further explored by \citet{Geletal1996}.
The full Bayesian treatment would
require choices of the priors of each component, including the parameter
covariance matrix. In the Supplemental Appendix, we provide Monte Carlo
results that shows that the resulting forecast performance is highly sensitive
to the choice of priors.
\section{Monte Carlo experiments}
\label{sec:MC}
We examine the finite-sample performance of the panel forecasting schemes in
the context of a dynamic heterogeneous panel data model using Monte Carlo
experiments.\footnote{
Further analytical results for a simple panel AR(1) model are provided in
Section~\ref{PanelAR} of the Appendix.} We allow for dynamics, parameter
heterogeneity, and correlations between the regressors and coefficients. The
forecasting methods are: (1) individual estimation which serves as the
benchmark against which other methods are compared, (2) pooled estimation,
(3) random effects, (4) fixed effects, (5) combination of individual and
pooled forecasts using the weights in~(\ref{w*NT}), (6) combination of
individual and FE forecasts using the weights in~(\ref{wFE}), (7) individual
forecast combination weights, and (8) EB forecasts.\footnote{
Additional results for equal weighted combinations and oracle weights are in
Section~S.4 of the Supplemental Appendix.}
Results do not vary greatly along the $N$ dimension, so we focus on the case
with $N=100$ and provide results for $N=1000$ in the Supplemental Appendix. The $
T$ dimension of the panel is more important, so we consider three different
values, $T=\{20,50,100\}$. The values of the parameters used in the
simulations are reported in Table~S.1 in Appendix~S.2.
\subsection{Data generating process}
\label{sec:MCDesignARX}
Our DGP augments a panel AR(1) model with an additional regressor,
\begin{equation}
y_{it}=\alpha _{i}+\beta _{i}y_{i,t-1}+
\gamma _{i}x_{it}+\varepsilon _{it},
\label{eq:MC1}
\end{equation}
where $\varepsilon _{it}=\sigma _{i}(z_{it}^{2}-1)/\sqrt{2}$ with
$
z_{it}\sim \mathit{iid}\mathrm{N}(0,1)$,
$\sigma _{i}^{2}\sim \mathrm{iid}
( 1+\chi _{1}^{2} ) /2$, and $x_{it}$ is generated as
\begin{equation}
x_{it}=\mu _{xi}+\xi _{it},
\label{xit}
\end{equation}
where
$
\xi _{it}=\rho _{xi}\xi _{i,t-1}+\sigma _{xi} ( 1-\rho _{xi}^{2}
) ^{1/2}\nu _{it}$, $ \nu _{it}\sim \mathrm{iidN} ( 0,1
) $, $\mu _{xi}=(z_{i}^{2}-1)/\sqrt{2}$,
$z_{i}\sim \mathrm{iidN} ( 0,1 ) $, and
$\sigma _{xi}^{2}\sim \mathrm{iid} ( 1+\chi _{1}^{2} ) /2$,
for individual units $i=1,2,\ldots ,N$, and observation periods
$t=1,2,\ldots ,T$. The autocorrelation coefficient of $x_{it}$ is
$
\rho _{xi}\sim \mathrm{iid} \text{Uniform}(0,0.95)$, allowing for a high
degree of dynamic heterogeneity in the regressors.
The coefficients of the lagged dependent variables, $y_{i,t-1}$, are generated
as $ \beta _{i}=\beta _{0}+\eta _{i\beta }$, with
$ \eta _{i\beta }\sim \mathrm{iid}\operatorname{Uniform} (-a_{\beta }/2,a_{\beta }/2)$
and $0\leq a_{\beta }<2(1- \vert \beta _{0} \vert )$.
To allow for correlated heterogeneity, we set
\begin{equation}
\alpha _{i}=\alpha _{0i}+\phi \mu _{xi}+
\sigma _{\eta }\eta _{i},\quad \text{and}\quad
\gamma
_{i}=\gamma _{0i}+\pi \mu _{xi}+\sigma
_{\zeta }\zeta _{i}, \label{aci}
\end{equation}
where $\eta _{i},\zeta _{i}\sim \mathrm{iidN}(0,1)$ and
$\alpha _{0}=\mathrm{
E} ( \alpha _{i} ) =\alpha _{0i}+\phi \mathrm{E} (
\mu _{xi} ) =\alpha _{0i}$. We examine three settings:
\begin{itemize}
\item $\alpha _{0i}=2/3$ if $i\leq N/2$, $\alpha _{0i}=4/3$ if
$i>N/2$, $
\sigma _{\alpha }^{2}=0.5$, $\gamma _{0i}=0.1$, and
$\sigma _{\gamma }^{2}=a_{\beta }=0$
\item $\alpha _{0i}=2/3$ if $i\leq N/2$, $\alpha _{0i}=4/3$ if
$i>N/2$, $
\sigma _{\alpha }^{2}=0.5$, $\gamma _{0i}=0.2/3$ if $i\leq N/2$,
$\gamma _{0i}=0.4/3$ if $i>N/2$, $\sigma _{\gamma }^{2}=0.1$, and
$a_{\beta }=0.5$
\item $\alpha _{0i}=2/3$ if $i\leq N/2$, $\alpha _{0i}=4/3$ if
$i>N/2$, $
\sigma _{\alpha }^{2}=1$, $\gamma _{0i}=0.2/3$ if $i\leq N/2$,
$\gamma _{0i}=0.4/3$ if $i>N/2$, $\sigma _{\gamma }^{2}=0.2$, and
$a_{\beta }=1$
\end{itemize}
\noindent Note that nonzero correlations need not bias the pooled estimates. What
matters for pooled estimates is the correlation between
$y_{i,t-1}^{2}, x_{it}^{2}$ and the individual coefficients.
Using (\ref{xit}) and (\ref{aci}), we have
\begin{align*}
\mathrm{E} \bigl[ x_{it} ( \gamma _{i}-\gamma
_{0} ) \bigr] &=
\mathrm{E} \bigl[ (\mu _{xi}+
\xi _{it}) ( \pi \mu _{xi}+\sigma _{
\zeta }\zeta
_{i} ) \bigr] =\pi \mathrm{E} \bigl( \mu _{xi}^{2}
\bigr) \neq 0,
\\
\mathrm{E} \bigl[ x_{it}^{2} ( \gamma _{i}-
\gamma _{0} ) \bigr] &=
\mathrm{E} \bigl[ (\mu
_{xi}+\xi _{it})^{2} ( \pi \mu
_{xi}+ \sigma _{\zeta }\zeta _{i} ) \bigr] =
\pi \mathrm{E} \bigl( \mu _{xi}^{3} \bigr) .
\end{align*}Therefore,
$\mathrm{E} [ x_{i,t-1}^{2} ( \gamma _{i}-\gamma _{0}
) ] =0$ if $\mu _{xi}$ are draws from a symmetric distribution
around $0$. To rule out this possibility, we draw $\mu _{xi}$ from a chi-square
distribution. To control the degree of correlated heterogeneity, we first
note that (taking expectations with respect to both $
i$ and $t$)
\begin{align*}
\mathrm{E} ( \gamma _{i} ) & =\gamma _{0},\qquad
\operatorname{Var}(\gamma _{i})=\pi ^{2}+\sigma
_{\zeta }^{2},
\\
\mathrm{E} ( x_{it} ) & =\mathrm{E} ( \mu _{xi}+\xi
_{it} ) =0, \qquad \operatorname{Var} ( x_{it} ) =\mathrm{E} (
x_{it}-\mu _{xi} ) ^{2}=\sigma
_{xi}^{2},
\end{align*}
and
$\mathrm{E} [ \operatorname{Var} ( x_{it} ) ] =
\mathrm{E}
( 1+\chi _{1}^{2} ) /2=1$. Also, since $\nu _{it}$ is distributed
independently of $\eta _{j}$ and $\zeta _{j}$ for all $t$, $i$, and
$j$, $
\operatorname{Cov} ( \gamma _{i},x_{it} ) =\pi $ and
$\operatorname{Corr} ( \gamma _{i},x_{it} ) =\pi ( \sigma _{
\zeta }^{2}+\pi ^{2} ) ^{-1/2}$. While heterogeneity is generally
correlated in AR panel models (\citet{PesSmi1995}), this setup allows
us to study further the role of correlated heterogeneity by varying the
correlation between the coefficient $
\gamma _{i}$ and $x_{it}$ as measured by $\rho _{\gamma x}$. To achieve
a given level of $\operatorname{Corr} ( \gamma _{i},x_{it} ) =\rho _{
\gamma x}$, we set
\begin{equation}
\pi = \frac{\rho _{\gamma x}\sigma _{\zeta }}{ \bigl( 1-\rho
_{\gamma x}^{2} \bigr) ^{1/2}}. \label{eq:pi}
\end{equation}
Similarly, to achieve $\operatorname{Corr} ( \alpha _{i},x_{i,t-1} ) =
\rho _{\alpha x}$, we set
\begin{equation}
\phi = \frac{\rho _{\alpha x}\sigma _{\eta }}{ \bigl( 1-\rho
_{\alpha x}^{2} \bigr) ^{1/2}}. \label{eq:phi}
\end{equation}
Defining
$\sigma _{\gamma }^{2}=\operatorname{Var}(\gamma _{i})=\pi ^{2}+\sigma _{
\zeta }^{2}$, we can use (\ref{eq:pi}) to see that
$\pi =\rho _{\gamma x}\sigma _{\gamma }$. An equivalent result emerges
for $\phi $ where, for
$
\sigma _{\alpha }^{2}=\operatorname{Var}(\alpha _{i})$, we have
$\phi =\rho _{\alpha x}\sigma _{\alpha }$. We thus use the parameters
$\sigma _{\alpha }^{2}$, $\sigma _{\gamma }^{2}$, and $a_{\beta }$ to vary
the degree of parameter heterogeneity in $\alpha _{i}$,
$\gamma _{i}$, and $\beta _{i}$, respectively.
We set $\xi _{i0}=0$ and initialize $y_{i0}$ as
$y_{i0}\sim \mathrm{iidN}
( \mu _{iy0},\sigma _{iy0}^{2} ) $ with
$
\mu _{iy0}=\frac{\alpha _{i}+\gamma _{i}\mu _{xi}}{1-\beta _{i}^{2}}$, $
\sigma _{iy0}^{2}=
\frac{\gamma _{i}^{2}\sigma _{xi}^{2}+\sigma _{i}^{2}}{
1-\beta _{i}^{2}}$, We also experimented with initialization schemes that
started the DGP on values away from the long run equilibrium, which did
not change the results qualitatively.
Since the forecast combinations use
$\boldsymbol{w}
_{i,T+1}=(1,y_{iT},x_{i,T})^{\prime }$ as an input, in the simulations
we set $\boldsymbol{w}_{i,T+1}$ as
$\boldsymbol{w}_{i,T+1}= ( 1,E ( y_{it} ) +\kappa _{i}
\sqrt{\operatorname{Var}(y_{it})},\mu _{xi}+\kappa _{i}\sigma _{xi} ) ^{
\prime }$, where $E ( y_{it} ) $ and
$
\operatorname{Var}(y_{it})$ are derived by assuming $y_{it}$ is stationary and
conditional on the model's parameters.\footnote{It is easily established that
$\mathrm{E} ( y_{it} ) =
\frac{\alpha _{i}+\beta _{i}\mu _{xi}}{1-\beta _{i}}$ and
$\operatorname{Var}(y_{it})=\frac{
\sigma _{i}^{2}}{1-\beta _{i}^{2}}+ (
\frac{\gamma _{i}^{2}\sigma _{xi}^{2}}{1-\beta _{i}^{2}} )
( 1+\frac{2\beta _{i}\rho _{xi}}{
1-\beta _{i}\rho _{xi}} ) $.}
The panel forecasts are evaluated using the ratio of the average MSFE of
method $j$ (pooled, fixed effects, random effects, empirical Bayes, and
the forecast combinations) measured relative to that of the reference individual
forecasts
\begin{equation*}
\text{rMSFE}_{j}= \frac{\frac{1}{NR}\sum
_{i=1}^{N}\sum_{r=1}^{R}(y_{i,T+1,r}-
\hat{y}_{i,T+1,j,r})^{2}}{\frac{1}{NR}\sum
_{i=1}^{N}
\sum
_{r=1}^{R}(y_{i,T+1,r}-
\hat{y}_{i,T+1,b,r})^{2}},
\end{equation*}
where $b$ denotes the benchmark forecast, which is the individual forecast.
Replications are denoted by $r=1,2,\ldots ,R$, where $R=10{,}000$.
\subsection{Simulation results}
Monte Carlo simulation results are reported in Table~\ref{tbl:MC_ARX_N100}.
Our theoretical analysis shows that the term $h_{NT}$ that adversely affects
forecasts from the individual estimates depends on the value of
$ \Vert \boldsymbol{w}_{i,T+1}-\mathrm{E}[\boldsymbol{w}_ {i,T+1}]
\Vert $, with small values of these deviations leading to better
forecasting performance for the individual estimates. To examine this effect,
we present two sets of conditional forecasting performance results, namely
for $\kappa _{i}=0$, that is, when $\boldsymbol{w}_{i,T+1}$ is set to its
mean
$\mathrm{E}
(\boldsymbol{w}_{it})= ( 1,\mathrm{E} ( y_ {it} ) ,
\mu _{xi}+\kappa _{i}\sigma _{xi} ) $ in the top panel and when
$\boldsymbol{
w}_{i,T+1}$ deviates from its mean by generating forecasts conditional
on
$
\boldsymbol{w}_{i,T+1}= ( 1,\mathrm{E} ( y_ {it} ) +
\kappa _{i}
\sqrt{\operatorname{Var}(y_{it})},\mu _{xi}+\kappa _ {i}\sigma _{xi} )
^{\prime }$ in the bottom panel. We set $\kappa _ {i}=1$ for
$i\leq N/2$, and $\kappa _{i}=-1$, for $i>N/2$.
We vary the parameter that controls the degree of correlated heterogeneity ($
\rho _{\gamma x}$) across three blocks of results and examine different
combinations of the two hyperparameters that determine the degree of heterogeneity,
$
a_{\beta }$ and $\sigma _{\alpha }^{2}$. Finally, we vary the time-series
dimension ($T$) along the columns.
\begin{sidewaystable}\thisfloatpagestyle{empty}
\caption{Monte Carlo results}
\label{tbl:MC_ARX_N100}
\hspace{-3em}
\resizebox{23cm}{!}{
\setlength\tabcolsep{3pt}
\begin{tabular}{lllllllllllllllllllllllllllllll}
\hline\hline
$a_\beta$&$\sigma^2_\alpha$ &&
\multicolumn{3}{c}{Pooled}&&\multicolumn{3}{c}{RE}&&\multicolumn{3}{c}{FE} &&
\multicolumn{3}{c}{Empirical Bayes}&&
\multicolumn{3}{c}{Comb.\ (pool)} &&
\multicolumn{3}{c}{Comb.\ (FE)}&&
\multicolumn{3}{c}{Comb.\ $\omega_i^*$}\\
\cline{1-2}\cline{4-6}\cline{8-10}\cline{12-14}\cline{16-18}\cline{20-22}\cline{24-26}\cline{28-30}
&$T$&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}&&\multicolumn{1}{c}{20}&\multicolumn{1}{c}{50}&\multicolumn{1}{c}{100}\\
\hline
\vspace*{-0.975 em}\\
\multicolumn{30}{c}{Conditional on $\kappa_i=0$}\\
\hline
\multicolumn{30}{c}{$\rho_{\gamma x}=0$}\\
0.0&0.5&&0.864&0.985&1.010&&0.911&0.985&0.996&&0.923&0.987&0.997&&0.935&0.989&0.997&&0.913&0.982&0.995&&0.955&0.992&0.998&&0.947&0.994&0.999\\
0.5&0.5&&0.860&0.981&1.006&&0.953&1.004&1.005&&0.978&1.009&1.007&&0.944&0.992&0.998&&0.911&0.981&0.995&&0.979&0.999&1.000&&0.939&0.994&0.999\\
1.0&1.0&&0.819&0.964&0.994&&1.089&1.096&1.061&&1.167&1.119&1.068&&0.937&0.990&0.998&&0.895&0.975&0.992&&1.022&1.006&1.000&&0.900&0.986&0.998
\vspace*{0.1 em}\\
\multicolumn{30}{c}{$\rho_{\gamma x}=0.5$}\\
0.0&0.5&&0.862&0.982&1.007&&0.910&0.985&0.996&&0.923&0.987&0.997&&0.937&0.989&0.997&&0.912&0.981&0.995&&0.955&0.992&0.998&&0.950&0.994&0.999\\
0.5&0.5&&0.856&0.977&1.001&&0.950&1.003&1.005&&0.977&1.009&1.007&&0.946&0.993&0.998&&0.910&0.980&0.994&&0.979&0.999&1.000&&0.942&0.994&0.999\\
1.0&1.0&&0.820&0.964&0.994&&1.098&1.103&1.065&&1.174&1.125&1.073&&0.941&0.991&0.998&&0.895&0.975&0.992&&1.024&1.006&1.000&&0.900&0.986&0.998
\vspace*{0.25 em}\\
\multicolumn{30}{c}{Conditional on $\kappa_i=\pm 1$}\\
\hline
\multicolumn{30}{c}{$\rho_{\gamma x}=0$}\\
0.0&0.5&&0.744&0.997&1.065&&0.733&0.924&0.971&&0.752&0.928&0.972&&0.804&0.945&0.979&&0.806&0.948&0.985&&0.842&0.951&0.980&&0.800&0.957&0.989\\
0.5&0.5&&0.867&1.163&1.243&&0.910&1.113&1.162&&0.949&1.122&1.164&&0.849&0.969&0.992&&0.832&0.964&0.991&&0.913&0.985&0.996&&0.808&0.965&0.992\\
1.0&1.0&&1.049&1.455&1.572&&1.274&1.504&1.529&&1.388&1.540&1.540&&0.881&0.977&0.995&&0.855&0.975&0.995&&0.983&0.999&1.000&&0.786&0.963&0.993\vspace*{0.1 em}\\
\multicolumn{30}{c}{$\rho_{\gamma x}=0.5$}\\
0.0&0.5&&0.764&1.023&1.093&&0.731&0.923&0.971&&0.752&0.928&0.972&&0.803&0.945&0.979&&0.809&0.950&0.986&&0.842&0.951&0.980&&0.801&0.958&0.989\\
0.5&0.5&&0.894&1.196&1.278&&0.912&1.120&1.168&&0.955&1.130&1.171&&0.846&0.968&0.991&&0.836&0.966&0.992&&0.914&0.986&0.997&&0.809&0.965&0.992\\
1.0&1.0&&1.068&1.479&1.597&&1.287&1.515&1.535&&1.401&1.552&1.547&&0.878&0.974&0.993&&0.856&0.975&0.995&&0.986&0.999&1.000&&0.785&0.960&0.992\\
\hline\hline
\multicolumn{30}{p{25.5cm}}{\footnotesize{Notes:
The table reports the ratio of average MSFE for a given forecasting method over the average MSFE of the forecasts based on
individual estimates.
The forecasts are:
`Pooled' based on pooled estimation, `RE' based on the random effects estimation,
`FE' based on the fixed effects estimation, `Empirical Bayes' based on the empirical Bayes estimation,
`Comb.~(pool)' refers to the combination
of forecasts based on individual and pooled estimation,
`Comb.~(FE)' the combination of forecasts based on
individual and fixed effects estimation,
and `Comb.\ $\omega_i^*$' the combination forecasts using the individual weights of Pesaran et al.~(2022).
The parameters
$a_\beta$ and $\sigma_\alpha^2$ determine the heterogeneity of the slope coefficient and the intercept.
The results in the upper panel are for $\kappa_i = 0$ where $\bs w_{i,T+1}$ equals the expected values of the regressors.
The results in the lower panel are for $\kappa_i=\pm 1$ where $\bs w_{i,T+1}$ equals the expected values of the regressors plus or minus one standard deviation.
Results are for \textit{PR}$^2$ of approximately 0.6 and $N=100$.
The DGP is set out in Section~\ref{sec:MCDesignARX}.}}
\end{tabular}
}
\end{sidewaystable}
With little heterogeneity and a small time-series dimension, $T=20$, consistent
with Propositions \ref{prop:individual} and \ref{prop:pooled}, pooling
yields an MSFE up to 25\% lower than the individual forecast with the gain
being largest when the predictor is far from its mean ($\kappa _{i}=
\pm 1$). However, the advantage of the pooled forecasts over the individual
forecasts vanishes quickly for the two larger values of $T$ and turns to
distinctly worse performance under larger parameter heterogeneity---particularly
when the predictors are away from their means.
The RE estimator produces the most accurate forecasts when parameter heterogeneity
is limited to the intercept ($a_{\beta }=0$, $\sigma _{\alpha }^{2}=0.5$)
and the predictor is far from its mean. When slope coefficients are heterogeneous,
this method yields quite poor forecasting performance that deteriorates
with $T$. Similar findings hold for the forecasts based on the FE method.
Forecast accuracy for both RE and FE methods tend to worsen (relative to the benchmark
forecasts) under correlated heterogeneity.
Regardless of the level of heterogeneity in parameters (whether correlated
or not), the empirical Bayes forecasts perform very well particularly for
the smallest sample size ($T=20$). Unlike forecasts based on the pooled,
RE or FE estimators, the empirical Bayes forecasts have the attractive
feature that they never perform worse, on average, than the benchmark.
These forecasts perform particularly well when the predictor is away from
its mean value.
Among the three forecast combinations, the cross-sectional averaging scheme
that combines the pooled and individual forecasts generally performs better
than the fixed effect combination scheme and also, in some cases, improves
on the EB forecasts. When $\boldsymbol{w}_{i,T+1}$ is far away from its
mean, $T$ is small, and parameter heterogeneity is high, the combination
scheme with individual weights performs particularly well, including relative
to the EB forecast.
\section{Empirical applications}
\label{sec:applications}
We next apply our set of panel forecasting methods to two empirical applications
on house price inflation in U.S. metropolitan areas and inflation in CPI
subindices. These applications represent quite different levels of in-sample
fit: For the CPI data, the pooled $\mathrm{R}^{2}$ ($
\mathrm{PR}^{2}$) of our models is around 0.2 while for house prices it
exceeds 0.8.
\subsection{Measures of forecasting performance}
Our empirical applications compute the out-of-sample MSFE as
$\mathrm{MSFE}
_{ij}=(T-T_{1})^{-1}\*\sum_{t=T_{1}}^{T-1}(y_{i,t+1}-\hat{y}_{i,j,t+1})^{2}$,
where $\hat{y}_{i,j,t+1}$ is the forecast of $y_{i,t+1}$ using method
$j$ and information known at time $t$. Each forecast in the test sample,
$\hat{y}
_{i,j,t+1}$, is generated using a rolling estimation window of observations
$
t-w+1,t-w,\ldots ,t$, where $w$ is the length of the rolling window, which
we set to $w=60$ in both applications. As in the simulations, we report
the ratio of the average MSFE of method $j$ relative to the average MSFE
for the benchmark forecasts ($b$) from the individual-specific model
$\mathrm{rMSFE}
_{j}= ( N^{-1}\sum_{i=1}^{N}\mathrm{MSFE}_{ij} ) / ( N^{-1}
\sum_{i=1}^{N}\mathrm{MSFE}_{{ib}} ) $. We also report the proportion
of units in the cross-section for which each method produces a smaller
MSFE than the benchmark along with the proportion of units in the cross-section
for which each method has the smallest or largest MSFE value.
Similar to the simulation study, we distinguish between forecasts where
the regressors are close to their means and when they are one standard deviation away from their means. Unlike in the Monte Carlo experiments, the parameters are unknown
in the two applications, and we therefore select forecasts based on
$d_{i,T+1}=
\hat{\boldsymbol{\theta }}_{i}^{\prime }\boldsymbol{w}_{i,T+1}$. Regressors
are said to be in the neighborhood of the mean of $d_{it}$ when
$|d_{i,T+1}-
\bar{d}_{i}-\kappa _{i}s_{d}|<c\sigma _{d}$, where $\bar{d}_{i}$ is the
mean and $s_{d}$ the standard deviation of $d_{it}$ in the estimation sample,
$
t=1,2,\ldots ,T$, and $c=0.1$. $\kappa _{i}=0$ then gives the results where
the predictors are close to their mean and $\kappa _{i}=\pm 1$ shows the
results when the predictors are one standard deviation away from the mean.
Additionally, we report results for all forecasts.
We examine the significance of any differences in forecast accuracy using
the \citet{DieMar1995} (DM) test of predictive accuracy both for
the panel as a whole and for the individual series. First, we use the panel
version of the DM test proposed by \cite{Pesetal2013}, which tests the
null that the MSFE generated by the individual forecasts, averaged both
across time and units, is equal in expectation to the equivalent MSFE generated
by the panel models.\footnote{The panel DM test first computes the difference between the cross-sectional
average squared forecast error at a given point in time for the benchmark
versus competing model. It then uses the time series of these average squared
forecast errors to compute Newey--West HAC standard errors that account
for serial dependencies.} Second, we apply the DM test to the $N$ forecasts
for individual units in the sample and report the number of significant
values in either direction and the number of insignificant test statistics.
The tests are set up so that negative values indicate that the panel forecasts
are more accurate than the individual forecasts, while positive values
of the DM tests indicate that the individual forecasts are more accurate.
For simplicity, we report results for all forecasts.
\subsection{U.S.\ house prices}
Our first application uses quarterly data on real house price inflation
in 377 U.S. Metropolitan Statistical Areas (MSAs) from the first quarter
of 1975 to the first quarter of 2023, which we obtain from the Freddie
Mac website.\footnote{For each MSA, house prices are calculated by deflating the Freddie Mac
house price index by the CPI.} Our forecasts target the one-quarter-ahead
MSA-level rate of house price log changes. After accounting for the necessary
presample and the estimation window, the first forecast is for 1991Q2 and
the last for 2023Q1, a total of 128 forecasts per MSA.
Our prediction model for the house price inflation rate in quarter
$t$ for MSA $i$, $y_{it}$, takes the form
\begin{equation}
y_{it}=\alpha _{i}+\beta _{i}y_{i,t-1}+
\beta _{i}^{\ast }y_{i,t-1}^{
\ast }+\gamma
_{Ri}\bar{y}_{i,t-1}^{(R)}+\gamma
_{Ci}\bar{y}_{t-1}^{(C)}+
\varepsilon _{it}, \label{eq:spatial_model}
\end{equation}
where $i=1,2,\ldots ,N$ denotes individual MSAs and
$t=1,2,\ldots ,T$ refers to the time period,
$y_{it}^{\ast }=\sum_{k=1,k\neq i}^{N}\omega _{ik}^{s}y_{kt}$ is the spatial
effect for a set of spatial weights $\omega _{ik}^{s}$,
$\bar{y}_{it}^{(R)}$ is the average house price inflation in the region
of unit $i$, and $\bar{y}_{t}^{(C)}$ is the countrywide average house price
inflation. The weights, $\omega _{ik}$, measure the spatial effect of house
prices in MSA $k$ on house prices in MSA $i$ and are based on geographic
distance, that is, $\omega _{ik}^{s}=v_{ik}/\sum_{k=1}^{N}v_{ik} $ and
$v_{ik}=1$ if MSAs $ ( i,k ) $ are at most 100 miles apart and
is zero otherwise. We obtain the weights from the data set of \citet{Yan2021}
and exclude MSAs without neighbors within 100 miles, which leaves 362 MSAs
in our sample.
The top panel in Table~\ref{tbl:applications_msfe} reports the results.
The column labeled ``all'' shows results
averaged across the full test sample, while columns labeled
$\kappa _{i}=0$ and $\kappa _{i}=\pm 1$ show results for subsamples in
which the predictor vector is close to the mean and one standard deviation
away from the mean, respectively. In the first three columns, the first
row shows the cross-sectional average MSFE value for the forecasts based
on individual estimates. Subsequent rows report ratios of the mean of the
individual MSFE for the respective methods relative to the benchmark forecasts.
Values below unity show that the ratio of average MSFE performance (across
MSAs) is better for the method listed in the row than for the benchmark
while values above unity indicate the opposite. The next three columns
headed ``freq. beating benchmark'' report the proportion of MSAs for which the respective
methods have a smaller MSFE than the benchmark, while the columns headed
``freq. smallest MSFE'' and ``freq. largest MSFE'' show the proportion
of MSAs for which the respective methods have the smallest or largest MSFE
among all forecasting methods.
\begin{sidewaystable}\thisfloatpagestyle{empty}
\caption{Results for the applications}
\label{tbl:applications_msfe}
\centering
{\footnotesize
\hspace*{-1.5cm}
\begin{tabular}{llllllllllllllll}
\hline\hline
& \multicolumn{3}{l}{Ratio of} && \multicolumn{3}{l}{Freq.\ beating} &&\multicolumn{3}{l}{Freq.\ smallest}&&\multicolumn{3}{l}{Freq.\ largest}\\
& \multicolumn{3}{l}{ave.~MSFE} &&\multicolumn{3}{l}{benchmark} &&\multicolumn{3}{l}{MSFE}&&\multicolumn{3}{l}{MSFE}\\
\cline{2-4}\cline{6-8}\cline{10-12}\cline{14-16}
Observations & all & $\kappa_i=0$ & $\kappa_i=\pm 1$ && all & $\kappa_i=0$ & $\kappa_i=\pm 1$&& all & $\kappa_i=0$ & $\kappa_i=\pm 1$ && all& $\kappa_i=0$ & $\kappa_i=\pm 1$ \\
\hline
\vspace*{-.5 em}\\
\multicolumn{16}{l}{House price inflation forecasts}\\
\hline
Individual & 2.822&2.520 & 3.542 && -- & -- & -- && 0.008&0.273 & 0.146 && 0.569 &0.282 & 0.420\vspace*{.2 em}\\
Pooled & 0.920&1.162 & 0.947 && 0.613&0.381 & 0.536 && 0.116&0.171 & 0.249 && 0.157 &0.188 & 0.160\\
RE & 0.924&1.166 & 0.960 && 0.619&0.376 & 0.528 && 0.108&0.022 & 0.047 && 0.003 &0.019 & 0.008\\
FE & 0.936&1.186 & 0.980 && 0.591&0.381 & 0.517 && 0.055&0.108 & 0.105 && 0.251 &0.376 & 0.296\\
Emp.Bayes & 0.901&0.955 & 0.881 && 0.942&0.519 & 0.652 && 0.185&0.157 & 0.127 && 0.003 &0.064 & 0.052\\
Comb.\ (pool) & 0.920&0.961 & 0.932 && 0.939&0.522 & 0.688 && 0.157&0.077 & 0.099 && 0.000 &0.011 & 0.017\\
Comb.\ (FE) & 0.937&0.977 & 0.940 && 0.917&0.494 & 0.677 && 0.019&0.072 & 0.069 && 0.011 &0.039 & 0.044\\
Comb.\ ($\omega^*_i$) & 0.921&0.957 & 0.909 && 0.936&0.541 & 0.713 && 0.044&0.028 & 0.052 && 0.006 &0.003 & 0.006\vspace*{.5 em}\\
\multicolumn{16}{l}{CPI inflation forecasts}\\
\hline
Individual &15.501&10.451 &11.295&& -- & -- & -- && 0.005&0.134 & 0.070 && 0.439&0.214 & 0.316\vspace*{.2 em}\\
Pooled & 0.878& 1.013& 0.971 && 0.444&0.374 & 0.417 && 0.203&0.118 & 0.112 && 0.396&0.406 & 0.380\\
RE & 0.880& 1.001& 0.957 && 0.508&0.390 & 0.401 && 0.016&0.070 & 0.064 && 0.000&0.037 & 0.032\\
FE & 0.883& 0.992& 0.959 && 0.508&0.401 & 0.401 && 0.000&0.086 & 0.102 && 0.166&0.230 & 0.225\\
Emp.Bayes & 0.892& 0.991& 0.926 && 0.984&0.652 & 0.818 && 0.390&0.278 & 0.278 && 0.000&0.070 & 0.011\\
Comb.\ (pool) & 0.930& 0.987& 0.953 && 0.733&0.481 & 0.572 && 0.128&0.091 & 0.123 && 0.000&0.027 & 0.016\\
Comb.\ (FE) & 0.935& 0.980& 0.967 && 0.791&0.524 & 0.583 && 0.053&0.064 & 0.070 && 0.000&0.016 & 0.021\\
Comb.\ ($\omega^*_i$) & 0.897& 0.972& 0.931 && 0.973&0.695 & 0.813 && 0.203&0.160 & 0.182 && 0.000&0.000 & 0.000\\
\hline\hline
\multicolumn{16}{p{19.3cm}}{\footnotesize{
results for the house price application and the bottom panel reports the
results for the CPI subindices application. The first three columns report
the ratio of average MSFE of the respective method in the row relative to
that of the individual forecast. The exception is the individual forecast,
which reports the average MSFE (times $10^5$ in the case of CPI). The second
three columns report the proportion of cross-section units for which the
respective method in the rows have a lower MSFE than the individual
forecast. The third three columns report the proportion of cross-section
units for which the respective method in the row has the lowest MSFE. The
last three columns report the proportion of units for which the respective
method in the row has the highest MSFE. In each block, the first column
averages over all forecasts, the second over the forecasts for which $d_
{i,T+1}=\hat{\boldsymbol{\theta }}_{i}^{\prime }\boldsymbol{w}_{i,T+1}$ is
close to its mean in the estimation sample. The third column averages over
the forecast for which $d_{i,T+1}$ is close to plus or minus one standard
deviation from its mean in the estimation sample. The methods in the rows
are listed in the footnote of Table~\ref{tbl:MC_ARX_N100}.}}
\end{tabular}
}
\end{sidewaystable}
Across the full sample, the average MSFE ratio below one for the pooled,
RE, and FE forecasts. However, these methods do notably worse than the
forecasts based on individual estimates when the predictors are close to
their mean ($
\kappa _{i}=0$). Empirical Bayes forecast produce the best overall MSFE
performance, reducing the MSFE of the benchmark by 10\%, followed by reductions
of 6--8\% among the three forecast combination schemes. The EB forecasts
perform particularly well when the predictors are far away from their mean.
While the proportional reductions in MSFE ratios may not seem very large,
they translate into very high frequencies of beating the benchmark. The
EB forecasts produce lower MSFE values than the benchmark for 94\% of the
housing price series followed by 92--94\% for the forecast combinations
but only 59--62\% for the pooled, RE, and FE forecasts.
Turning to evidence of individual forecasts being ``best'' or ``worst,'' for the full test sample
the benchmark forecasts only produce the smallest MSFE for 1\% of the variables
versus 19\% for the EB and 16\% for pooled forecast combination schemes.
Using this metric, again the benchmark forecasts perform much better when
the predictors are close to their sample mean for which they are most accurate
for 27\% of the MSAs versus 16\% and 8\% for the EB and pooled combinations,
respectively. Conversely, forecasts based on individual estimates are worst
overall for 57\% of the variables versus 1\% or less for the EB and forecast
combination schemes.
These results show that the EB and combination approaches offer the attractive
feature of not only improving on the MSFE values of the baseline
``on average'' but, equally importantly,
rarely producing markedly worse forecasts than the baseline and often generating
substantially better results. Interestingly, the risk of producing the
highest MSFE value is notably lower for the pooled combination and individual
weighted combination than for the EB forecasts when the predictors are
close to their mean.\footnote{Equal-weighted combinations also performs quite well in both of the empirical
applications, which is a known feature in the forecast combination literature.}
Figure~\ref{fig:densities} summarizes our findings visually through density
plots fitted to the cross-sectional distribution of MSFE ratios for our
forecasting methods.\footnote{To reduce the number of lines, we do not plot the densities for the FE
and RE approaches, which are very similar to those from pooling.} MSFE
ratios have a widely dispersed, right-skewed distribution for the pooled
forecasts compared to the Bayesian and combination approaches whose distributions
are far more peaked and centered just below unity. This feature is highly
undesirable as it raises the likelihood of very poor forecasts for an individual
housing price series compared with that of the Bayesian and combination
approaches.\footnote{The impressive performance of the EB approach for the tail groups is consistent
with \citet{Efr2011}.}
\begin{figure}[tbp]
\caption{Distributions of ratios of MSFEs}
\label{fig:densities}\centering
\hspace*{1.5cm}{\footnotesize {House price forecasts \hspace{3cm} CPI
subindices forecasts}\newline
\hspace*{.3cm}
\includegraphics[scale=.4,trim={2.1cm 7.5cm 2.1cm
7.5cm},clip]{./figures/house_density_all_pdf_26May2025.pdf}
\includegraphics[scale=.4,trim={2.1cm 7.5cm 2.1cm
7.5cm},clip]{./figures/CPI_pdf_all_26May2025.pdf}\newline
\hspace*{.3cm}
\includegraphics[scale=.4,trim={2.1cm
7.5cm 2.1cm 7.5cm},clip]{./figures/house_density_mean_pdf_26May2025.pdf}
\includegraphics[scale=.4,trim={2.1cm 7.5cm 2.1cm
7.5cm},clip]{./figures/CPI_pdf_mean_26May2025.pdf}\newline
\hspace*{.3cm}
\includegraphics[scale=.4,trim={2.1cm 7.2cm 2.1cm
7.5cm},clip]{./figures/house_density_std_pdf_26May2025.pdf}
\includegraphics[scale=.4,trim={2.1cm 7.2cm 2.1cm
7.5cm},clip]{./figures/CPI_pdf_std_26May2025.pdf}\newline
}
{\footnotesize \vspace{1em}
\parbox{13cm}{\footnotesize{Notes:
The graphs show density plots of the ratios of MSFEs for the house price
application in the left column and those for the CPI subindices application
in the right column. In the first row are the density plots for the MSFEs
from all forecasts, in the second row for the forecasts for which
$d_{i,T+1} = \hat{\boldsymbol{\theta}}_{i}^{\prime}\boldsymbol{w}_{i,T+1}$
is close to its mean in the estimation sample, and in
the third row for the forecast for which $d_{i,T+1}$ is close to plus or minus one
standard deviation from its mean in the estimation sample. The
density estimates use a normal kernel with a bandwidth 0.04. The forecasting
methods are listed in the footnote of Tables~\ref{tbl:MC_ARX_N100}.
}} }
\end{figure}
The first and second rows of Table~\ref{tbl:DM} reports panel DM test statistics
and the number of cross-sectional units with a DM test below
$
-1.96$ (panel forecasts are significantly more accurate) or above 1.96
(individual-specific forecasts are significantly more accurate), respectively,
for each application.
\begin{table}
\caption{Diebold-Mariano test statistics for equal predictive accuracy}
\label{tbl:DM}\centering
\resizebox{17cm}{!}{
\begin{tabular}{lrrrrrrr}
\hline\hline
& Pooled & RE & FE & Emp.Bay.& Comb(pool)& Comb(FE) & Comb($\omega^*_i$) \\
\hline
\multicolumn{8}{l}{House Prices: all forecasts}\\
\hline
Panel DM & $-$9.45 & $-$9.11 & $-$7.57& $-$24.63 & $-$22.93 & $-$21.19 &$-$27.59\\
$\text{DM}<-1.96$/$\text{DM}>1.96$ & 60/6 & 62/6 & 57/8 & 209/2& 189/0 & 169/0 & 240/0
\vspace*{0.25 em}\\
\multicolumn{8}{l}{CPI: all forecasts}\\
\hline
Panel DM & $-$7.95& $-$7.78& $-$7.56& $-$11.25&$-$11.67& $-$10.59 & $-$11.48\\
$\text{DM}<-1.96$/$\text{DM}>1.96$ & 35/60 & 33/45 & 32/42 & 134/0 & 56/23 & 50/10 & 137/0\\
\hline\hline
\multicolumn{8}{p{17.3cm}}{\footnotesize{Notes:
The row ``Panel DM'' reports the results of the panel version of the
Diebold-Mariano test of Pesaran et al.~(2013). The second row report unit by
unit Diebold-Mariano test results: ``$\text{DM}<-1.96$'' reports the number of units
with a DM test statistic smaller than $-1.96$ and ``$\text{DM}>1.96$'' shows the
number of units whose test statistic exceeds 1.96. The remaining units have
insignificant DM test statistics. In total the house prices panel consist of
362 units and the CPI panel of 187 units. Each test is for the null
hypothesis that the forecasting method in the columns has equal forecast
accuracy as the forecasts based on individual estimates. The forecasting
methods are listed in the footnote of Table~\ref{tbl:MC_ARX_N100}.
}}
\end{tabular}}
\end{table}
The panel DM tests show that the EB and combination forecasts are significantly
more accurate than the individual forecasts ``on average''
as well as for a large portion of the individual series (between 169 and
240 MSAs), while the opposite only happens for two individual MSAs in the
case of the EB forecasts. Pooled, RE and FE panel forecasts are also significantly
more accurate than the individual forecasts on average as well as for between
57 and 62 of the individual MSAs and significantly less accurate for very
few MSAs.
\subsection{CPI inflation of sub-indices}
Our second application covers inflation rates for up to 187 subindices
of the U.S. consumer price index (CPI) obtained from the FRED database.
The data is measured at the monthly frequency and spans the period from
January 1967 to December 2022. Again, we use rolling estimation windows
with 60 observations and require each estimation sample to be balanced,
excluding individual series without a complete set of observations in a
given window. After accounting for the necessary pre-samples, we generate
up to 599 forecasts for each series, with the first forecast computed for
February 1973.
We consider an autoregressive forecasting specification with lags 1, 2,
and 12 augmented with lagged values of the first principal component of
the data, the default yield and term spread.
The bottom panel of Table~\ref{tbl:applications_msfe} shows that, for the
full test sample, all forecasting methods produce lower MSFE values than
the benchmark. The pooled, RE, FE, and EB forecasts reduce the average
MSFE of the benchmark by around 12\%, while the forecast combination methods
reduce it by 7--10\%. Interestingly, when the predictors are close to their
sample mean, the lowest MSFE ratios are produced by the three forecast
combination methods, while conversely the EB scheme performs best when
the predictors are further removed from their mean.
The EB forecasting scheme performs particularly well overall, beating the
benchmark model's accuracy for 98\% of the variables followed by 97\% for
the individual weights, 73--79\% for the pooled and FE combinations and
around 50\% for the RE and FE schemes. As in the first application, these
percentages are notably lower for predictors close to their mean and higher
further away.
The EB forecasts also produce by far the highest frequency with the smallest
MSFE values overall (39\%) followed by 20\% for the pooled and individual
forecast combination scheme. This is matched by very low probabilities
of producing the worst forecast, which never occurs in our sample for the
EB method or any of the three forecast combination schemes but is far more
likely to occur for the benchmark (43.9\%) and pooled forecasts (39.6\%).
Our evidence is summarized by the probability density plots for the MSFE ratios
in the right panels of Figure~\ref{fig:densities}. The figure clearly highlights
the pronounced dispersion and thick right tails of the MSFE-ratio distribution
for the pooled forecasts. The distributions of MSFE ratios of the EB and
combination approaches are far more concentrated and less asymmetrical.
For values of the predictors farther away from the mean, the tails of the
densities are somewhat thicker, with the EB approach standing out as having
the thinnest right tail, and hence, the lowest probability of generating
forecasts less accurate than those from the individual-specific benchmark.
Turning to the DM test results for the CPI inflation data in Table~\ref{tbl:DM}, all panel models generate significantly negative DM panel
test statistics and so their associated forecasts are significantly more
accurate, on average, than the individual forecasts. The pooled, RE, and
FE models perform somewhat worse in this application, as the number of
individual CPI series for which their forecasts are significantly more
accurate than the individual-specific forecasts is smaller than those for
which the opposite holds. Conversely, the EB and combination forecasts
continue to be significantly more accurate than the benchmark forecasts
for between 50 and 137 of the individual CPI series and are only significantly
less accurate for between zero and 23 series. The EB and individual combination
approaches perform particularly well in this application.
\section{Conclusion}
We provide a comprehensive examination of the out-of-sample predictive
accuracy of a large set of novel and existing panel forecasting methods,
including individual estimation, pooled estimation, random effects, fixed
effects, empirical Bayes, and forecast combinations.
Our main findings can be summarized in three points. First, we find that
many panel forecasting approaches perform systematically better than forecasts
based on individual estimates. For panels with a small or medium-sized
time-series dimension $T$---a setting relevant to many empirical applications
in economics---our Monte Carlo simulations and empirical applications demonstrate
sizeable gains both on average and for the majority of individual units
from exploiting panel information.
Second, our analytical results and Monte Carlo simulations show that one
should not expect a single forecasting approach to be uniformly dominant
across applications that differ in terms of the cross-sectional and time-series
dimensions, strength of predictive power, and degree of heterogeneity in
intercept and slope coefficients along with how correlated this heterogeneity
is.
Forecasts based on pooled estimates are most accurate only in situations
with little or no parameter heterogeneity and a small $T$ dimension, while
forecasts based on FE and RE estimates perform relatively well mainly when
heterogeneity is confined to model intercepts and $T$ is small. Neither
of these approaches perform well in settings with high levels of heterogeneity
where individual-specific forecasts tend to perform better, particularly
if $
T$ is relatively large. By overweighting forecasts that perform well and
underweighting forecasts that perform poorly, forecast combination and
empirical Bayes methods manage to produce the most accurate forecasts across
a broad range of settings.
Third, the panel forecasting methods differ in terms of their ability to
reduce the probability of generating very poor forecasts for individual
units in a cross-section. While the individual, pooled, random and fixed
effect estimation methods perform poorly in some of the simulations and
empirical applications, the forecast combination and empirical Bayes methods
rarely generate the least accurate forecasts for individual units and retain
some probability of being the best forecasting method. These panel forecasting
approaches therefore come out on top of our analysis.
In a nutshell, our simulations and empirical applications suggest that
forecast combinations and Bayesian panel methods offer insurance against
poor performance. Compared to the alternative forecasting methods we consider,
this better ``risk-return'' trade-off makes the combination and Bayes methods attractive in forecast
applications with panel data.
\newpage