Forecasting with panel data: Estimation uncertainty versus parameter heterogeneity

Forecasting with panel data: Estimation uncertainty versus parameter heterogeneity

\thispagestyle{empty}

\setcounter{page}{1}

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 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 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 LinSmi1972, LeeGri1979, and 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 Bal2008 (Bal2008,Bal2013) consider the forecasting performance of the best linear unbiased predictor (BLUP) of 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. 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. BruSil2006 consider a similar group of methods to forecast migration data and find that fixed effects and shrinkage estimators perform best; see PicTim2024 for a review of the literature.

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 Pesetal2013 and BooPic2020, where the target is also significant improvements in forecast accuracy rather than a significant change in parameters.

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) introduces the model setup and our assumptions, while Section (ref) derives analytical results on the predictive accuracy of individual, pooled, and FE forecasting schemes. Section (ref) introduces our forecast combination schemes. Section (ref) describes the empirical Bayes estimator. Section (ref) 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.

Setup and assumptions

We begin by describing the panel regression setup and assumptions used in our analysis.

Panel regression model

Our analysis considers the following linear panel regression model:

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)) 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.

Assumptions

Our theoretical analysis builds on a set of standard assumptions about the underlying data generating process.

Under Assumption (ref), 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). Assumption (ref) 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) 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) is required when we compare average MSFEs based on individual and pooled estimators. It provides sufficient conditions under which (see Lemma (ref))

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) 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) and (ref), $ \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) below, we need the stronger conditions

still allowing $N$ to rise much faster than $T$.\footnote{As noted by Fanetal2015 (Fanetal2015, Section (ref)), 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) $\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}$,

Further, under Assumption (ref), $\mathrm{E} ( \boldsymbol{w}_{it} \boldsymbol{w}_{it}^{\prime } ) =\boldsymbol{Q}_{i}$, and it follows that

We next introduce assumptions that are required primarily for establishing the properties of pooled and fixed effects predictors.

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) and (ref) 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), $\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), $ \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)--(ref) are not required for forecasts based on the individual estimates and the associated MSFE. Assumption (ref) 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 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, 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

Taking expectations,

Assumptions (ref) and (ref) 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)) are allowed to have nonzero correlation with the regressors, but optimality of forecasts based on pooled estimates of $\boldsymbol{\beta } $ requires Assumption (ref) 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 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 $.}

Theoretical results on forecasting performance

We next use the setup and assumptions from Section (ref) to establish theoretical results on the forecasting performance of different modeling approaches. Section (ref) discusses forecasts based on individual and pooled estimation, and building on this, Section (ref) covers fixed effects forecasts.

Note that our theoretical framework can be equally applied to forecasts across groups instead of individuals, when there are a priori known groups such as industries or states within a given country. Pooled regressions can be applied to any given, 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.

Forecasts based on individual and pooled estimation

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).} Forecasts based on individual estimators take the form

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

where $\tilde{\boldsymbol{\theta }}=(\boldsymbol{W}^{\prime }\boldsymbol{W} )^{-1}\boldsymbol{W}^{\prime }\boldsymbol{y}$. Using ((ref)), ((ref)) and the definition of $\boldsymbol{\bar{\xi}}_{NT}$ in Assumption (ref),

Forecast errors from these schemes take the form

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:

where $S_{NT}=N^{-1}\sum_{i=1}^{N}s_{iT}$, $R_{NT}=N^{-1}\sum_{i=1}^{N}r_{iT }$, with elements

and

Under Assumptions (ref) and (ref), $\mathrm{E} ( r_{iT } ) =0$ and $\sup_{i,T}\mathrm{E} \vert r_{iT} \vert <C$, and under cross-sectional independence (Assumption (ref)) we have $ R_{NT}=O_{p}(N^{-1/2})$. Similarly, $\sup_{i,T}\mathrm{E} \vert s_{iT} \vert <C$,

$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) of the Appendix):

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)) can therefore be important. Parameter heterogeneity, in contrast, does not affect the accuracy of the forecasts in ((ref)). 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,

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$.

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)) in ((ref)), we can express the squared forecast error when pooled estimates are used as follows:

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)) and ((ref)), and $ \boldsymbol{\bar{\xi}}_{NT}$ is defined under Assumption (ref). After some algebra, and averaging over $i$, we have

where $\tilde{S}_{N,T+1}$, and $\tilde{R}_{N,T+1}$ are defined by equations ( (ref)) and ((ref)) in Section (ref) 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.

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 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) of the Appendix. See also Pesetal2024a where estimation of such models with short $T$ panels is considered.}

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) shows that the MSFE from the forecasts based on the individual estimates will be affected by an estimation error term of the form

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

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:

Hence, there exists a $T_{0}$ such that, for a fixed $T>T_{0}$, and as $N\rightarrow \infty $,

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$ . 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))). 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 $,

Similarly, when both $N$ and $T\rightarrow \infty $ (in any order)

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$.

Forecasts based on fixed effects estimation

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

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

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:

where

$\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

$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.

Forecast combinations

We next consider approaches that combine the forecasts from Section (ref) to minimize the MSFE.

Combinations of individual and pooled forecasts

Given the MSFE trade-off associated with the forecasts in ((ref)) and ((ref)), 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., 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. Lahetal2017 discuss a broader range of combination methods and Ell2017 provides an analysis of the effect on the combination weights and forecasting performance from having a large common component in the forecast errors.}

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

The value of $\omega $ that minimizes the average MSFE is therefore given by

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)) and ( (ref)), 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)).

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.

Forecast combinations with individual weights

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

where the optimal value of $\omega _{i}$ is given by

The weights again depend on the variances and covariances of the underlying forecast errors. Related to this, 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.

Combining individual and fixed effect forecasts

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

yielding the optimal pooled weight

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)) and ((ref)), 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) of the Appendix.

Estimation of combination weights

Estimates of the weights for the forecast combination in Proposition (ref) require estimates of $\Delta _{NT}$, $h_{NT}$, and $\psi _{NT}$. Under Assumption (ref), 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) . Using ((ref)) and ((ref)), the estimators of $\Delta _{NT}$ and $h_{NT}$ are given by

where $\tilde{\boldsymbol{\eta }}_{i}=\tilde{\boldsymbol{\theta }}- \boldsymbol{\hat{\theta}}_{i}$, and

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) that

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)), 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)), 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)

We now employ a half-jackknife estimator of $\boldsymbol{\theta }_{i}$ (DhaJoc2015, 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

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) that combine individual and fixed effects forecasts can be estimated in a similar fashion, with details provided in Appendix (ref).

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., Armetal2022 and Efr2016), we focus on empirical Bayes (EB) methods. The EB forecast uses the estimator of 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

$ \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

recalling that $\boldsymbol{Q}_{iT }=T^{-1}\boldsymbol{W}_{i}^{\prime } \boldsymbol{W}_{i}$ is invertible under Assumption (ref). 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 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 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)) fall somewhere between the common optimal weights and the EB weights.\footnote{While the EB estimator in ((ref)) is fully parametric, other studies pursue a nonparametric approach to the distribution of $\hat{\boldsymbol{ \theta }}_{i}$; see, for example, BroGre2009 and GuKoe2017, and more recently, Liu2023 and 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)) 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

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 LinSmi1972 and further explored by 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.

Monte Carlo experiments

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) 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)), (6) combination of individual and FE forecasts using the weights in ((ref)), (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.

Data generating process

Our DGP augments a panel AR(1) model with an additional regressor,

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

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

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:

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)) and ((ref)), we have

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$)

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 (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

Similarly, to achieve $\operatorname{Corr} ( \alpha _{i},x_{i,t-1} ) = \rho _{\alpha x}$, we set

Defining $\sigma _{\gamma }^{2}=\operatorname{Var}(\gamma _{i})=\pi ^{2}+\sigma _{ \zeta }^{2}$, we can use ((ref)) 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

where $b$ denotes the benchmark forecast, which is the individual forecast. Replications are denoted by $r=1,2,\ldots ,R$, where $R=10{,}000$.

Simulation results

Monte Carlo simulation results are reported in Table (ref). 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.

With little heterogeneity and a small time-series dimension, $T=20$, consistent with Propositions (ref) and (ref), 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.

Empirical 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.

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 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 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.

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

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 Yan2021 and exclude MSAs without neighbors within 100 miles, which leaves 362 MSAs in our sample.

The top panel in Table (ref) 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.

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) 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 Efr2011.}

The first and second rows of Table (ref) 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.

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.

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) 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). 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), 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.

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.