Inference for Forecasting Accuracy: Pooled versus Individual Estimators in High-dimensional Panel Data

Inference for Forecasting Accuracy: Pooled versus Individual Estimators in High-dimensional Panel Data

Introduction

Panel data arise in settings where information is observed across both individuals and time. In many applications, panels feature a large number of individuals $N$ and a non-negligible time dimension $T$, enabling researchers to study both cross-sectional heterogeneity and temporal dynamics. Most practically used panel models account for individual heterogeneity through unit-specific intercepts, but the vast majority impose (complete) homogeneity on the slope coefficients of interest. This modeling choice effectively implies the use of pooled estimators that aggregate data across individuals. While pooling can improve efficiency by reducing variance, it may lead to misleading conclusions when slope coefficients differ across units. Such heterogeneity is frequently observed in empirical work (see, for discussions, hsiao2008random, baltagi2008pool, browning2007heterogeneity among many others).

To make data-driven decisions about whether to pool or not, one strand of research—following the seminal work of Swamy—has developed tests for cross-sectional slope homogeneity. Examples include phillips-sul, Pesaran, and ando2015. These tests are typically derived under the assumption that model errors are independent across both individuals and time, with a few exceptions such as Blomquist, who consider scenarios with small $N$ and $T \to \infty$. While cross-sectional slope homogeneity is certainly sufficient to justify pooling, it is not a necessary condition. In large-$N$ panels, complete slope homogeneity is often unrealistic, and consequently a number of less restrictive criteria have been explored in the literature. One important example involves latent group structures, where slopes are homogeneous within but not across groups (see su2016identifying; wang2018homogeneity). A recent extension by wang2024homogeneity considers models with a large number of covariates. A quantitative measure for slope homogeneity in large-$N$, large-$T$ settings was proposed by kutta:dette:2024, where small values of this measure suggest the use of pooled estimators. Another, more direct approach to assessing the cost of pooling was introduced by Campello, who developed tests for the hypothesis of no slope heterogeneity bias. While this hypothesis can still be practically restrictive for large $N$, this contribution is important by widening the discussion about what counts as a reasonable pooling criterion.

In this paper, we investigate a different pooling criterion that is directly motivated by prediction—a primary concern for applied users. Rather than testing for precise slope homogeneity (possibly after grouping), we ask which estimators—individual or pooled—yield lower prediction errors. This approach is motivated by the insight that in most practical applications, some degree of slope heterogeneity is unavoidable, even within groups. The key question is therefore whether heterogeneity can be safely ignored or whether it meaningfully affects forecasting accuracy. Our pooling criterion is based on the mean squared forecasting error (MSFE), first suggested for this purpose by pesaran:pick:timmermann:2022. While we adopt minimal MSFE as a decision rule, the statistical model and asymptotic framework developed here differ substantially from those in pesaran:pick:timmermann:2022. The most important conceptual distinction is that we allow for fixed (non-random) slopes, whereas pesaran:pick:timmermann:2022 derive their inference methods within the random coefficient model of Swamy. The latter framework is particularly suitable when the goal is consistent estimation of a mean slope coefficient. However, as we argue in this paper, assuming that individual slopes—or their deviations from the mean—follow a specific distribution is both unnecessary and potentially restrictive. Moreover, pesaran:pick:timmermann:2022 develop, technically speaking, an oracle procedure for fixed $T$ and $N \to \infty$, though in practice a large $T$ seems necessary for variance estimation. Our approach, in contrast, is fully feasible and provably consistent for large time and cross-sectional dimension. We defer a detailed comparison to Remark (ref) below. We summarize the main contributions of this paper as follows:

The remainder of this paper is organized as follows: Statistical methodology is developed and theoretically justified in Section (ref). Finite sample properties of our approach are studied in Section (ref). Proofs and additional simulations are located in the Appendix.

Inference for the prediction error in panel data

In this section, we develop new methodology to quantify the difference of prediction errors between pooled and individual slope estimators. Methods are theoretically justified in an asymptotic regime, where both the length of the time series $T$ and the size of the cross-section $N$ are large. To derive meaningful probabilistic limits, theory is developed for regimes of moderate heterogeneity (Assumption (ref)), where the decision for the superior forecasting method is most challenging. The main statistical result is the Gaussian approximation for the distribution of the estimated forecasting error in Theorem (ref), which implies confidence intervals for the forecasting error (ref). We begin this section, by introducing some mathematical concepts which are important for our subsequent mathematical analysis.

Mathematical preliminaries

Strongly mixing panels In this work, we study panels of dependent random variables. Dependence is captured by $\alpha$-mixing, with weak dependence expressed by fast decaying mixing coefficients. In the following, we provide a general outline of $\alpha$-mixing for multivariate panels of random variables doukhan:1994. For an overview of mixing concepts, we refer the reader to bradley:2007.\\ Let $\mathcal{M}\subset {\mathbb Z}^d$ for some $d \in {\mathbb N}$ be endowed with the maximum norm $\|\cdot \|_\infty$ and $(\Xi_{z})_{z\in \mathcal{M} }$ be a panel of random variables indexed in $\mathcal{M}$. For index sets $\mathcal{I}$ and $\mathcal{J}$ we define the distance w.r.t. the maximum norm as $$ \mbox{dist}(\mathcal{I}, \mathcal{J}) = \min(\|i-j\|_\infty: i \in \mathcal{I}, j \in \mathcal{J}) $$ and denote by $\mathcal{F}_\mathcal{I} =\sigma\big(\Xi_{z}: z \in \mathcal{I}\big)$ the sigma field generated by the random variables $\{ \Xi_{z}: z \in \mathcal{I}\} $. For $r \in {\mathbb N}_0$ the $r$th $\alpha$-mixing coefficient is defined by

The panel is called $\alpha$-mixing if $\alpha(r) \to 0$ as $r \to \infty$. In the case of this work, we will consider the dimension $d=2$, even though extensions to higher dimensional panels are possible.

Conditional Landau symbols One characteristic of this work is the development of conditional inference methods for panel data. We therefore have to clarify notions of conditional convergence of random variables. Suppose that $(X_n)_n, (Y_n)_n$ are sequences of random variables and $(a_n)_n$ is a sequence of positive, real numbers. Then we say, that

Statistical methodology

Panel model We consider the linear regression panel

where $x_{i,t}$ is a $K$-dimensional vector of regressors, $\beta_i$ is a $K$-dimensional vector of slope coefficients and $\varepsilon_{i,t}$ a centered, real model error with unknown distribution. We call $t$ the time component of the panel and $i$ the individual or cross-sectional component. We collect all equations concerning the $i$th individual in the model

where $X_{i}=(x_{i,1},...,x_{i,T})' $ is a regression matrix of dimension $T \times K$ and $y_{i}=(y_{i,1},...,y_{i,T})'$ and $\varepsilon_{i}=(\varepsilon_{i,1},...,\varepsilon_{i,T})'$ are $T-$dimensional vectors. Sometimes we refer to all regressors collectively and therefore define the compounded matrix $\mathbf{X}=(X_1,...,X_N)$. Often panel models comprise constant, individual specific intercepts, which are omitted in model (ref) for simplicity, and would practically be removed. We illustrate this case in our simulation study in Section (ref).

Estimators For model (ref), we define the individual ordinary least squares (OLS) slope estimator for the $i$th individual as \[ \hat \beta_i := [X_i' X_i]^{-1} X_i' y_i. \] Next, drawing on data from all $N$ individuals, we define the pooled version \[ \hat \beta^{pool}:= \Big( \sum_{i=1}^N X_i' X_i \Big)^{-1} \sum_{i=1}^N X_i' y_i. \] The estimator $\hat \beta^{pool}$ is commonly used under the assumption of slope homogeneity $\beta_1=...=\beta_N$, where it is substantially more efficient than the individual estimators $\hat \beta_1,...,\hat \beta_N$ due to its smaller variance. However, the effectiveness of the pooled estimator relies on the individual slopes being, if not the same, at least very similar - otherwise, it can be severely biased. A standard way to assess slope homogeneity is the use of slope homogeneity tests as discussed in the Introduction. While these tests are very powerful for larger $N$, their very power often makes them oversensitive to even minor inhomogeneities of slopes, discouraging pooling even when practically beneficial. In the below discussion, we therefore try to shift the subject, away from somewhat stylized assumptions on the true slopes, towards the comparative merits of the estimators in terms of forecasting.

Individual and pooled MSE The performance of the individual estimators $\hat \beta_1,...,\hat \beta_N$ and the pooled version $\hat \beta^{pool}$ can be assessed by comparing their prediction accuracy for a new vector of predictor-response pairs $(x_{i,T+1}, y_{i,T+1})_{i=1,...,N}$. Thus, following pesaran:pick:timmermann:2022, we invoke the mean squared prediction error (MSPE) in $(x_{i,T+1},y_{i,T+1})_{i=1,...,N}$ conditionally on the known matrix of regressors $\mathbf{X}=(X_1,...,X_N)$. More precisely, we define respectively the individual and pooled prediction errors as

Notice that in this formulation, the vector of predictors $(x_{i,T+1})_{i=1,...,N}$ is non-random - the user determines in which predictors they would like to make the comparison. Selecting $(x_{i,T+1})_{i=1,...,N}$ means specifying a scenario where forecasts are of interest and we seek to determine which forecasting method is most suitable for it. The responses at time $T+1$ are defined as $y_{i,T+1}:=x_{i,T+1}'\beta_i+\varepsilon_{i,T+1}$ with model errors $\varepsilon_{i,T+1}$ and the expectations in $E^{ind}$ and $E^{pool}$ are taken over all model errors $\varepsilon_{i,t}$, $i=1,...,N$ and $t=1,...,T+1$.

A closed form for the prediction error For our formal analysis of the prediction error, we impose some mathematical assumptions. In the following, let $\|\cdot\|_2$ denote the Euclidean norm (Frobenius norm) for vectors and matrices.

Conditions $i)$ and $ii)$ permit the existence of complex error structures. Error distributions are non-parametric and only some polynomial moments are required. Moreover, the errors can be dependent across space and time simultaneously. More precisely, the covariance matrix of the errors $(\varepsilon_1^\prime, \ldots ,\varepsilon_N^\prime)^\prime$ is the Kronecker product $\Sigma_N \otimes \Sigma_T$, where the first factor captures cross-sectional and the second factor temporal dependence. This structure is more general than those typically assumed in the literature and which are special cases of this setting. For example, the traditional assumptions of independent, homoscedastic errors is captured by $\sigma^2 \cdot I_N \otimes I_T$, where $\sigma^2>0$ is the variance and $I_N, I_T$ are the identity matrices of dimension $N$ and $T$, respectively. Independent errors with individual-specific variances are similarly captured by $D \otimes I_T$, where $D=diag(\sigma_1^2,...,\sigma_N^2)$. The case of cross-sectional dependence only, as used e.g. by Zellner, is incorporated by $\Sigma_N \otimes I_T$ (this is also closely related to ando2015). Allowing both factors to differ from the identity matrix (simultaneous temporal and cross-sectional dependence) obviously encompasses much richer models than have been treated before, particularly in high dimensional panels. Separable covariance structures are a standard tool in the analysis of spatio-temporal data and often appropriate for panels, where the individual component can, in a broad sense, be interpreted as a location. Condition $iii)$ is standard in the study of linear models and stronger than the common exogeneity assumptions used in the literature on panel data. We require it, to rigorously formulate weak convergence results conditionally on the regressors $\mathbf{X}$, even though we expect similar results to be true if errors and regressors are weakly dependent. Finally, Condition $iv)$ implies that the errors in the (hypothetical) time period $T+1$ are independent of all previous errors and maintain the same covariance structure. \\ Using Assumption (ref), we can give a closed form for the prediction errors defined in (ref) and (ref).

Our main object of interest is the difference $E^{ind}-E^{pool}$, because, e.g., $E^{ind}-E^{pool}<0$ implies that the pooled estimator is outperformed by the individual estimators. Lemma (ref) now implies that to understand $E^{ind}-E^{pool}$ we may study the terms $E_1,E_2,E_3$. \\

Analysis of $E_i$ We first have to impose some additional assumptions.

Assumptions of the sort imposed by Condition $i)$ are typical in the study of large panel data (see, e.g., Assumption 2 in Pesaran). The regressor matrices converge to respective limits $Q_i$ that are positive definite (smallest eigenvalue is positive) and have uniformly bounded norm. The second condition ensures that the forecasting error is calculated for $x_{i,T+1}$ that are not too large, which can be seen as a condition to avoid extrapolating into areas that are far apart from the original data. Finally, we assume that the model errors are weakly dependent along the space and time dimensions. This structure allows for general dependence patterns that are more realistic than traditional models with independent errors. The exponential decay condition on mixing coefficients is satisfied by most of the typical time series models such as ARMA processes Mo88. We also point out that in the Appendix we prove our results for even weaker, polynomial mixing conditions (see condition (ref)). Such conditions come at the cost of a more restrictive relation between $N$ and $T$ (see condition (ref)) and are therefore not further discussed here.

Analysis of the variance terms As a first step to analyze the error terms $E_1, E_2, E_3$ from Lemma (ref), we investigate their (asymptotic) order of magnitudes. We begin by studying the two terms $E_1$ and $E_3$ that are independent of the regression slopes and represent the variance of the individual estimators and the pooled estimator, respectively. Convergence is formulated conditionally on $\mathbf{X}$ to harmonize this with later results, but obviously, $E_1, E_3$ are $\mathbf{X}$-measurable and hence the derived rates hold conditionally and unconditionally.

The derived rates of convergence are highly intuitive. For the individual slope estimators that are based (each) on $T$ observations, the variance of predictions are of size $\approx C/T$. Similarly, the predictions based on the pooled estimator have variance of size $\approx C/(NT)$, which is much smaller. While the variance of the pooled predictions are much smaller, the pooled estimator can produce biased results if the slopes are different. Here, the bias is measured by $E_2$. Under the classical hypothesis of slope homogeneity ($\beta_1= \ldots =\beta_N$) it directly follows that $E_2=0$. In other, very heterogeneous scenarios, $E_2=\mathcal{O}_P(1)$ is possible such that $E_2$ dominates both $E_1$ and $E_3$. These two extremes illustrate situations, where a pooled estimator is either evidently superior (because of high homogeneity) or evidently inferior (because of high heterogeneity) compared to individual estimators. In applications, the situation is often less clear-cut. For this reason, we investigate in the following the practically more relevant intermediate case, where $E_2$ is of the same order of magnitude as $E_1$. In this scenario, statistical inference can be invoked to make an informed decision on which prediction method is better. We hence impose further assumptions.

The first condition is standard in the literature on slope homogeneity and guarantees that no single slope dominates the final test statistics. The second condition helps us focus on the main case of interest, where the two errors, individual and pooled, are of the same order of magnitude. Lemma (ref) entails that $E_3$ is negligible compared to $E_1$ and hence, for $E_1-E_2-E_3$ to be close to $0$, $E_2$ has to be of the same order as $E_1$. We provide a short calculation to illustrate that under moderate slope heterogeneity $E_2=\mathcal{O}_P(T^{-1})$ (just as $E_1$). Using Assumption (ref) part i), we have

Now, consider the right side: Under Assumption (ref) part ii), the object on the inside of the round bracket is of order $\mathcal{O}(1)$ and hence the entire term inside the curved bracket is also of size $\mathcal{O}(1)$, suggesting that $E_2 = \mathcal{O}_P(T^{-1})$ (see Lemma (ref)). It is clear that the formulation of Condition $ii)$ can be relaxed such that not all slopes have to be close to one another, but only the majority. For instance, one may impose that there exists an index set $\mathcal{N}\subset \{1,...,N\}$ of exceptions such that only the weaker assumption \[ \max_{i,j \in \{1,...,N\}\setminus \mathcal{N}}\|\beta_i-\beta_j\|\le c_5/\sqrt{T} \] holds. As long as $\mathcal{N}$ is small enough, say if $|\mathcal{N}|=o(\sqrt{N})$, the results in this section remain valid. The final Condition iii) in Assumption (ref) moderates the size of the temporal dimension of the panel, compared to its cross-sectional component. The assumption $N/T^2 \to 0$ has been used in classical slope homogeneity tests Pesaran and the fact that in our case $N/T^\eta \to 0$ for $\eta$ arbitrarily close to $2$ is required is a small additional price that we pay for allowing temporal and cross-sectional dependence. Finally, notice that Assumption (ref) does not impose any distribution on the slope coefficients as is typically the case in random coefficient models. While distributional assumptions are convenient when the focus is on consistent estimation of the average slope coefficient, these assumptions can prove to be restrictive when the main objective is prediction.

Statistical inference We now develop an inference method for the difference of prediction errors. As a first step, we propose the estimator

On the first glance, $\hat E$ may seem like a simple plug-in estimator for $E_2$ (see eq. (ref)). Yet, a careful analysis reveals that $\hat E$ actually approximates the error sum $E_1+E_2$. It can therefore serve as a building block to approximate our true object of interest $E_1-E_2$.

This lemma is a consequence of the much more precise analysis in the proof of Theorem (ref), where we study the weak convergence behavior of an appropriately standardized version of $\hat E$. We state the lemma here, to make the next step of our procedure more understandable. We recall that for our inference method, we do not need an estimator for $(E_1+E_2)$, but rather for $(E_1-E_2)$. Therefore, we supplement $\hat E$ with additional estimator $\hat E_1$ of $E_1$. We will then have \[ \hat E-2\hat E_1\approx E_2-E_1. \] Let us define the $T \times T$ matrix

where the entry

is an estimator of the autocovariance of lag $h$. Here $b$ is a regularization parameter (all $h$-diagonals with $h>b$ are set equal to $0$). Finally, we define the estimate

Notice that $ \hat E_1$ is the direct, empirical analogue to $E_1$ defined in (ref).

Notice that the above Lemmata now imply that \[ \sqrt{N}T\big\{(\hat E-2\hat E_1)-(E^{ind}-E^{pool})\big\} = \sqrt{N}T \big\{\hat E-(E_1-E_2)\big\}+o^{|\mathbf{X}}_P(1). \] We can thus develop statistical inference for the difference $E^{ind}-E^{pool}$, by studying the weak convergence of the statistic $\sqrt{N}T \{\hat E-(E_1-E_2)\}$. For this purpose, we define the following two terms

Therewith, we define the conditional asymptotic variance

As we show in the proof of Theorem (ref), $\tau_N^2$ is asymptotically close to the variance of $\hat E$ and can thus be used for standardization. As common in the study of dependent time series, we have to ensure that the variance does not asymptotically degenerate.

Let us define the standardized estimator

where $F_{\breve E}$ refers to the cumulative distribution function of $\breve E$, conditional on $\mathbf{X}$ and $\tau_N$ is defined in (ref).

The proof of Theorem (ref) is technically challenging in two respects. First, the derivation of conditional weak convergence requires a non-standard application of a Berry-Esseen theorem for random variables conditional on the regressors. In the scenario of moderate heterogeneity an unusual linearization of $\breve E$ occurs that features both linear and squared error terms. These different terms influence the variance $\tau_N^2$, where two terms in the curly brackets occur, one for the squared errors (first term) and one for the non-squared errors (last term). Second, due to the complex dependence structure, it is difficult to prove that standardizing by $\tau_N^2$ is asymptotically equal to standardizing by the true variance of $\hat E$. In the special case of Gaussian data, the proof turns out to be much simpler because the squared and linear errors are uncorrelated; if $Z_1, Z_2$ are jointly normally distributed, $\mathbb{E}[Z_1Z_2^2]=0$ regardless of the covariance structure. Yet, for general error distribution, the analysis is substantially more intricate.

A conditional confidence interval Theorem (ref) can be used to construct confidence intervals for the difference of prediction errors $E^{pool}-E^{ind}$. Suppose that the conditional variance $\tau_N^2$ is known. Then, an asymptotic $(1-\alpha)$ confidence interval for $E^{pool}-E^{ind}$ is given by

Notice that $\mathcal{C}_{1-\alpha}$ is an approximate $1-\alpha$ confidence interval in the sense that

Finite sample properties

We investigate the finite sample properties of our new confidence interval (CI) by means of a Monte Carlo study based on 5,000 iterations, with sample sizes $N\in\{100,500\}$ and $T\in\{10,15, 20,25,30, 40,60,80\}$. In each setup, we simulate $x_{i,t}$ as a $5$-dimensional vector, where each element is independently drawn from $\mathcal{N}(1,1)$. We report the empirical coverage rates of the feasible confidence interval $\hat{\mathcal{C}}_{1-\alpha}$, where $\tau_N$ in (ref) is replaced by the square root of $\hat{\tau}^2_N$ defined in (ref). As an infeasible benchmark, denoted as $\mathcal{C}_{1-\alpha}^*$, we compute $\mathcal{C}_{1-\alpha}$ under the assumption that the true values for $\tau_N$, $\Sigma_N$ and $\Sigma_T$ are known. Therefore, randomness in $\mathcal{C}_{1-\alpha}^*$ stems only from the estimation of the slope parameters. For both the feasible and the infeasible confidence interval, we report their average lengths $L(\hat{\mathcal{C}}_{1-\alpha})$ and $L(\mathcal{C}_{1-\alpha}^*)$, respectively. Across all simulations, $\alpha=0.05$.

We begin by simulating homogeneous slope parameters with spherical errors by setting $\beta_i=1$ for all individuals, $\Sigma_N=\mathrm{I}_{N\times N}$ and $\Sigma_T=\mathrm{I}_{T\times T}$ (Table (ref)). Next, we analyze the properties of our CI under heterogeneous slopes, i.e., $\beta_i=1$ for the first half of the sample and $\beta_i=2$ for the remainder, while maintaining spherical errors (Table (ref)). This is followed by results for homogeneous and heterogeneous slopes under serial correlation in the model errors (Tables (ref) and (ref)), where $\varepsilon_{i,t}=\phi\varepsilon_{i,t-1} + u_{i,t}$ with $\phi=0.3$ and $u_{it}\sim\mathcal{N}(0,1)$. Consequently, $\Sigma_T$ is a Toeplitz matrix with $(\Sigma_{T})_{s,t}=\frac{1}{1-\phi^2}\phi^{|t-s|}$. Motivated by the upper bound on $\rho$ in Lemma (ref), the bandwidth is set to $b=T^{2/7}$ in the simulations with dynamic errors. Finally, we report results for a data generating process containing unobserved individual fixed effects (Tables (ref) and (ref)), i.e., $Y_{i,t}=x_{i,t}'\beta_i+\alpha_i+\varepsilon_{i,t}$, where $\alpha_i\sim\mathcal{N}(\bar{x}_i,1)$ and $\bar{x}_i=T^{-1}\sum_{t=1}^T x_{i,t}$. In this case, our methodology can be applied by replacing $x_{i,t}$ with the demeaned version $\ddot{x}_{i,t}=x_{i,t}-\bar{x}$ and by adapting the estimator of the error covariance matrix to the demeaning operation by multiplying $M_T=\mathbf{I}_{T}- \mathbf{1}_T\mathbf{1}_T'/T$, where $\mathbf{1}_T$ denotes a vector of ones, from both sides to $\hat{\Sigma}^{(i,i)}(b)$.

In the Appendix, we report simulation results for heteroskedastic model errors, which are not covered by our theory (Tables (ref) - (ref)). Here, we set $\mathrm{var}(\varepsilon_{i,t})= |(x_{i,t})_1|$, i.e., the variance of $\varepsilon_{i,t}$ depends on the absolute value of the first component of $x_{i,t}$. We then illustrate that our CI can be made robust to heteroskedasticity by adapting the estimator of the covariance matrix of the model errors to this feature of the data generating process. Furthermore, we apply a similar parametric estimation strategy to models with dynamic model errors and show that this approach can lead to a substantial improvement in the empirical coverage rate of your CI. While we do not provide a formal proof, we further illustrate that our CI can be made robust to both serial correlation and heteroskedasticity in the model errors by using a heteroskedasticity and autocorrelation consistent (HAC) type estimator for the covariance matrix of the model errors. Finally, we present results for the heterogeneous model with individual fixed effects, where the slope coefficients $\beta_i$ are $\mathrm{i.i.d.}$ draws from a standard normal distribution (see Table (ref)). We find that our CI performs well, which is not surprising, as our theory does not impose any distributional assumptions on the slope coefficients.

Results

Across all simulations, the empirical coverage rate of the infeasible confidence interval $\mathcal{C}_{0.95}^*$ practically coincides with the desired nominal level so that $L(\mathcal{C}_{0.95}^*)$ provides a sensible benchmark for the length of the feasible confidence interval. Table (ref) illustrates that our feasible confidence interval is conservative in the homogeneous slopes model when the model errors are uncorrelated. Unsurprisingly, it is thus substantially wider than the infeasible CI for any sample size. We find that this is due to the overestimation of the conditional asymptotic variance $\tau_N^2$ under slope homogeneity, as $\Lambda=\Lambda_k=0$ in (ref) whereas $\hat \Lambda$ and $\hat \Lambda_i$ are non-zero in (ref). Arguably, when choosing between the pooled and the heterogeneous model for prediction, being conservative under slope homogeneity is somewhat desirable, as it implies that the right endpoint of the CI is rarely below zero, correctly indicating that the pooled estimator is preferable in terms of the prediction error. When slope coefficients are heterogeneous, the feasible CI closely approximates the infeasible CI, as shown in Table (ref), even when $T$ is only moderately large. A similar behavior can be observed when the model errors follow an AR(1) process (see Tables (ref) and (ref)), albeit with the additional requirement that $N$ should not be too large relative to $T$, as expected in the presence of time dependence. For instance, Tables (ref) and (ref) show that $\hat{\mathcal{C}}_{0.95}$ can undercover the true difference of prediction errors when $N$ is very large relative to $T$, since estimation error in $\hat{\Sigma}^{(i,i)}(b)$ affects the accuracy of the point estimate $\hat{E}-2\hat{E}_1$. However, $T$ only needs to be moderately large relative to $N$ (e.g., $N=500$ and $T\approx 30$), to make this effect negligible. As we illustrate in the Appendix, distortions in the empirical coverage rate can be more severe with higher levels of serial correlation (see Tables (ref) and (ref)) so that longer panels might be necessary for a satisfactory coverage rate. However, we also provide numerical evidence that the empirical coverage rate can be further improved when a parametric estimator is used to estimate the covariance matrix of the model errors instead of the nonparametric estimator $\hat{\Sigma}^{(i,i)}(b)$ (see Tables (ref) and (ref)). When the model errors are independent but the model contains an unobserved individual fixed effect, the data must be demeaned across time, leading to serial correlation in the demeaned model errors. Consequently, the results in Tables (ref) and (ref) resemble the ones in (ref) and (ref). Again, the feasible CI undercovers the true difference in prediction errors only when $N$ is very large relative to $T$. However, already when $T$ is as small as 20, the empirical level is close to the nominal level, and the length of the feasible CI closely approximates the length of the infeasible benchmark. As illustrated by Table (ref), our approach is also robust to individual fixed effects when the slope heterogeneity is not fixed but rather follows a random coefficients specification where $\beta_i$ is drawn from a standard normal distribution. As before, the empirical coverage rate is very accurate when $T$ is moderately large, i.e., $T\approx 20$.

Conclusion

Researchers frequently face a bias-variance trade-off when choosing between pooled and individual-specific estimators for the slope coefficients in panel data analysis. A sensible strategy is to use the estimator that yields a smaller mean squared prediction error. Here, we have derived a closed form expression of the difference in mean squared prediction errors between the pooled and the individual-specific OLS estimators for panel data models with potentially heterogeneous slopes. We have then constructed a novel confidence interval for the said difference in mean squared prediction errors. Our asymptotic analysis shows that our confidence interval has asymptotically the correct coverage as $N,T\to\infty$, while allowing for the cross-sectional dimension to grow at a faster rate than the time-dimension. By means of an extensive simulation study, we have demonstrated that the empirical coverage rate is close to the nominal coverage rate in sufficiently long panels, even when the cross-sectional dimension is much larger than the time series dimension. Finally, we have illustrated that our confidence interval can be flexibly adapted to features of the data generating process to further enhance its small-sample performance.

{\bf Acknowledgments. } Holger Dette has been partially supported by the Deutsche Forschungsgemeinschaft (DFG), project number 45723897, and by TRR 391 {\it Spatio-temporal Statistics for the Transition of Energy and Transport}, project number 520388526 (DFG). Tim Kutta's work has been partially funded by AUFF grants 47331 and 47222.