The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
55,970 characters
Inference for Forecasting Accuracy: Pooled versus Individual Estimators in High-dimensional Panel Data
\title{Inference for Forecasting Accuracy: Pooled versus Individual Estimators in High-dimensional Panel Data}
\author{
Tim Kutta\thanks{Department of Mathematics, Aarhus University, \texttt{[email removed]}; corresponding author} \and
Martin Schumann \thanks{School of Business and Economics, Maastricht University, \texttt{[email removed]}
} \and
Holger Dette\thanks{Fakultät für Mathematik, Ruhr-Universität Bochum, \texttt{[email removed]}
}
}
\maketitle
\begin{abstract}
Panels with large time $(T)$ and cross-sectional $(N)$ dimensions are a key data structure in social sciences and other fields. A central question in panel data analysis is whether to pool data across individuals or to estimate separate models. Pooled estimators typically have lower variance but may suffer from bias, creating a fundamental trade-off for optimal estimation.
We develop a new inference method to compare the forecasting performance of pooled and individual estimators. Specifically, we propose a confidence interval for the difference between their forecasting errors and establish its asymptotic validity. Our theory allows for complex temporal and cross-sectional dependence in the model errors and covers scenarios where $N$ can be much larger than $T$—including the independent case under the classical condition $N/T^2 \to 0$. The finite-sample properties of the proposed method are examined in an extensive simulation study.
\end{abstract}
\section{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, \citealp{hsiao2008random}, \citealp{baltagi2008pool}, \citealp{browning2007heterogeneity} among many others).
To make data-driven decisions about whether to pool or not, one strand of research—following the seminal work of \cite{Swamy}—has developed tests for cross-sectional slope homogeneity. Examples include \cite{phillips-sul}, \cite{Pesaran}, and \cite{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 \cite{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 \citealp{su2016identifying}; \citealp{wang2018homogeneity}). A recent extension by \cite{wang2024homogeneity} considers models with a large number of covariates. A quantitative measure for slope homogeneity in large-$N$, large-$T$ settings was proposed by \cite{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 \cite{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 \cite{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 \cite{pesaran:pick:timmermann:2022}. The most important conceptual distinction is that we allow for fixed (non-random) slopes, whereas \cite{pesaran:pick:timmermann:2022} derive their inference methods within the random coefficient model of \cite{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, \cite{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{rem:compa} below.
We summarize the main contributions of this paper as follows:
\begin{itemize}
\item[(1)] We develop the first asymptotic confidence interval for the difference between pooled and individual prediction errors in large data panels.
\item[(2)] Our methodology accommodates general error distributions and simultaneous temporal and cross-sectional dependence.
\item[(3)] We establish validity for large samples with $T, N \to \infty$, allowing for $N \gg T$. New asymptotics are derived under a “moderately heterogeneous” model, where the decision whether to pool is most challenging and statistical inference most valuable.
\end{itemize}
The remainder of this paper is organized as follows: Statistical methodology is developed and theoretically justified in Section \ref{section-2}. Finite sample properties of our approach are studied in Section \ref{section-simulations}. Proofs and additional simulations are located in the Appendix.
\section{Inference for the prediction error in panel data}\label{section-2}
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 \textit{moderate
heterogeneity} (Assumption \ref{ass_3}), 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{theo_main}, which implies confidence intervals for the forecasting error \eqref{e:defCI}. We begin this section, by introducing some mathematical concepts which are important for our subsequent mathematical analysis.
\subsection{Mathematical preliminaries} \label{sec_prelim}
\noindent \textbf{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 \cite[see][]{doukhan:1994}. For an overview of mixing concepts, we refer the reader to \cite{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
\begin{align*}
\alpha(r)=\sup &\Big\{ |\mathbb{P}(A \cap B)- \mathbb{P}(A) \mathbb{P}(B)|: A \in \mathcal{F}_\mathcal{I}, ~
B \in \mathcal{F}_\mathcal{J},~
dist(\mathcal{I}, \mathcal{J})\ge r,~ \mathcal{I}, \mathcal{J} \subset \mathcal{M}\Big\}.
\end{align*}
The panel is called $\alpha$-mixing if $\alpha(r) \to 0$ as $r \to \infty$.
\noindent In the case of this work, we will consider the dimension $d=2$, even though extensions to higher dimensional panels are possible.
\smallskip
\textbf{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
\begin{equation} \label{e:stochLan}
\begin{split}
X_n =\, &\mathcal{O}_P^{|Y_n}(a_n) \,\, \Leftrightarrow \,\, \forall \delta>0: \lim_{C \to \infty}\limsup _n \mathbb{P}\Big (\mathbb{P} \big (|X_n|/a_n > C~|~Y_n \big )>\delta\Big ) =0,\\
X_n =\, & o_P^{|Y_n}(a_n) \,\,\,\, \Leftrightarrow \,\, \forall C, \delta>0:\,\,\lim _n \mathbb{P}\Big (\mathbb{P} \big (|X_n|/a_n > C~|~ Y_n \big )>\delta\Big ) =0.
\end{split}
\end{equation}
\subsection{Statistical methodology}
\noindent \textbf{Panel model} We consider the linear regression panel
\begin{equation} \label{model_1}
y_{i,t} =x_{i,t}' \beta_i+\varepsilon_{i,t} \quad \quad i=1,...,N, \,\, t=1,...,T,
\end{equation}
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 \textit{time component} of the panel and $i$ the \textit{individual} or \textit{cross-sectional component}. We collect all equations concerning the $i$th individual in the model
\begin{equation} \label{model_2}
y_{i} =X_{i} \beta_i+\varepsilon_{i} \quad \quad i=1,...,N,
\end{equation}
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 \eqref{model_2} for simplicity, and would practically be removed. We illustrate this case in our simulation study in Section \ref{section-simulations}.
\smallskip
\textbf{Estimators} For model \eqref{model_1}, 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.
\smallskip
\textbf{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 \cite{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
\begin{align}\label{E_ind}
E^{ind}:= & \frac{1}{N}\sum_{i=1}^N E_i^{ind} := \frac{1}{N}\sum_{i=1}^N{\mathbb E} \big[ ( x_{i,T+1}'\hat \beta_i -y_{i,T+1})^2\big| \mathbf{X}\big],\\
E^{pool}:=& \frac{1}{N}\sum_{i=1}^N E_i^{pool} := \frac{1}{N}\sum_{i=1}^N{\mathbb E} \big[ (x_{i,T+1}'\hat{\beta}^{pool} -y_{i,T+1})^2\big| \mathbf{X}\big]. \label{E_pool}
\end{align}
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$.
\smallskip
\textbf{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.
\begin{assumption} \label{ass_1}
\begin{itemize} $ $
\item[i)] (Moments) For some $M \in \mathbb{N}$ with $M\ge 16$ and a constant $C>0$, the errors $\varepsilon_{i,t}$ and regressors $x_{i,t}$ satisfy
\[
\max_{i,t}\mathbb{E}|\varepsilon_{i,t}|^M \le C<\infty, \qquad \max_{i,t}\mathbb{E}\|x_{i,t}\|_2^M \le C<\infty.
\]
\item[ii)] (Error covariance) The errors $\varepsilon_{i,t}$ are centered. The vector of errors $\varepsilon = (\varepsilon_1',....,\varepsilon_N')' \in \mathbb{R}^{NT}$ has
a covariance matrix $\Sigma = \Sigma_N \otimes \Sigma_T$, where "$\otimes$" denotes the Kronecker product, $\Sigma_T \in \mathbb{R}^{T \times T}$ is a covariance matrix of a stationary process and $\Sigma_N \in \mathbb{R}^{N \times N}$ is another covariance matrix.
\item[iii)] (Exogeneity) The vector of errors $\varepsilon$ and the matrix of regressors $\mathbf{X}$ are independent of each other.
\item[iv)] (Prediction) The regressors $x_{i,T+1}$ are fixed, non-random vectors. The vector $(\varepsilon_{1,T+1},...,\varepsilon_{N,T+1})$ is independent of $(\varepsilon, \mathbf{X})$, centered and has covariance matrix $(\Sigma_T)_{1,1}\cdot \Sigma_N$.
\end{itemize}
\end{assumption}
\noindent 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 \cite{Zellner}, is incorporated by $\Sigma_N \otimes I_T$ (this is also closely related to \cite{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{ass_1}, we can give a closed form for the prediction errors defined in \eqref{E_ind} and \eqref{E_pool}.
\begin{lem} \label{lem_1} Suppose that Assumption \ref{ass_1} holds, then
\begin{align*}
E_i^{ind} = & (\Sigma_N)_{i,i} Tr \Big[ \Sigma_T \Big(X_i [X_i' X_i]^{-1}x_{i,T+1}x_{i,T+1}'[X_i' X_i]^{-1}X_i' \Big)\Big]+(\Sigma_N)_{i,i}(\Sigma_T)_{1,1}. \\
E_i^{pool} = &\bigg(x_{i,T+1}'\Big( \sum_{j=1}^N X_j' X_j \Big)^{-1} \sum_{j =1}^N X_j' X_j (\beta_j-\beta_i)\bigg)^2\\
& +Tr \bigg[ \Big( \sum_{j=1}^N X_j' X_j \Big)^{-1} x_{i,T+1} x_{i,T+1}'\Big( \sum_{j=1}^N X_j' X_j \Big)^{-1}\sum_{j,k}(\Sigma_N)_{j,k} X_j' \Sigma_T X_k\bigg]+(\Sigma_N)_{i,i} (\Sigma_T)_{1,1}.
\end{align*}
and in particular $E^{ind}-E^{pool} = E_1 -E_2 -E_3$, where
\begin{align}
E_1&:=
Tr \Big[ \Sigma_T \sum_{i=1}^N\Big(\frac{1}{N}(\Sigma_N)_{i,i} X_i [X_i' X_i]^{-1}x_{i,T+1}x_{i,T+1}'[X_i' X_i]^{-1}X_i' \Big)\Big] \label{e:E_1}\\
E_2&:=\frac{1}{N}\sum_{i=1}^N\bigg(x_{i,T+1}'\Big( \sum_{j=1}^N X_j' X_j \Big)^{-1} \sum_{j =1}^N X_j' X_j (\beta_j-\beta_i)\bigg)^2 \label{e:E_2} \\
E_3& :=Tr \bigg[ \Big( \sum_{j=1}^N X_j' X_j \Big)^{-1}\frac{1}{N}\sum_{i=1}^N \big\{x_{i,T+1} x_{i,T+1}'\big\}\Big( \sum_{j=1}^N X_j' X_j \Big)^{-1}\sum_{j,k}(\Sigma_N)_{j,k} X_j' \Sigma_T X_k\bigg]\label{e:E_3}
.
\end{align}
\end{lem}
\noindent 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{lem_1} now implies that to understand $E^{ind}-E^{pool}$ we may study the terms $E_1,E_2,E_3$. \\
\smallskip
\textbf{Analysis of $E_i$} We first have to impose some additional assumptions.
\begin{assumption}$ $ \label{ass_2}
\begin{itemize}
\item[i)] (Regressor convergence) For each $i=1,...,N$ there exists a positive definite matrix $Q_i$, such that
\begin{align*}
\max_i \|X_i'X_i/T-Q_i\|_2 \overset{\mathbb{P}}{\to} 0, \quad 0< c_1 \le \min_i \lambda_{min}(Q_i), \quad \max_i \|Q_i\|_2 \le c_2.
\end{align*}
Here, $\lambda_{min}$ denotes the smallest eigenvalue of a matrix and $c_1, c_2$ some fixed positive constants.
\item[ii)] (Boundedness $x_{i,T+1}$) There exists a constant $c_3>0$, such that $\max_i \|x_{i,T+1}\|_2 \le c_3$.
\item[iii)] (Errors) The errors $\{\varepsilon_{i,t}: 1 \le i \le N, 1 \le t \le T\}$ form a strongly mixing field, with mixing coefficients satisfying for all $r$ sufficiently large, and some constant $\psi>0$ that
$\alpha(r) \le \psi^{-r}$.
\end{itemize}
\end{assumption}
\noindent Assumptions of the sort imposed by Condition $i)$ are typical in the study of large panel data (see, e.g., Assumption 2 in \cite{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 \cite[see][]{Mo88}. We also point out that in the Appendix we prove our results for even weaker, polynomial mixing conditions (see condition \eqref{e:mix:app}). Such conditions come at the cost of a more restrictive relation between $N$ and $T$ (see condition \eqref{e:eta:app}) and are therefore not further discussed here.
\smallskip
\textbf{Analysis of the variance terms} As a first step to analyze the error terms $E_1, E_2, E_3$ from Lemma \ref{lem_1}, 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.
\begin{lem} \label{lem_2}
Suppose that Assumptions \ref{ass_1} and \ref{ass_2} hold, then
\[
E_1 = \mathcal{O}^{|\mathbf{X}}_P(T^{-1}), \qquad \textnormal{and} \qquad E_3 = \mathcal{O}_P^{|\mathbf{X}}((NT)^{-1}).
\]
\end{lem}
\noindent 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.
\begin{assumption} \label{ass_3}
\begin{itemize} $ $
\item[i)] (Bounded slopes) The slopes $\beta_1, \beta_2,...$ are uniformly bounded $\max_i\|\beta_i\|\le c_4$ for some $c_4>0$.
\item[ii)] (Moderate heterogeneity) $\max_{i,j} \|\beta_i-\beta_j\| \le c_5/\sqrt{T}$ for some $c_5>0$.
\item[iii)] (Intersect-time relation) There exists a constant $\eta \in (0,2)$ such that $N/T^\eta\to 0$.
\end{itemize}
\end{assumption}
\noindent 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{lem_2} 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{ass_2} part i), we have
\begin{align*}
E_2= & \frac{1}{NT}\sum_{i=1}^N\bigg(x_{i,T+1}'\Big[ \frac{1}{NT}\sum_{j=1}^N X_j' X_j \Big]^{-1} \frac{1}{N}\sum_{j =1}^N \frac{X_j' X_j }{T}[\sqrt{T}(\beta_j-\beta_i)] \bigg)^2 \\
\approx & \frac{1}{T}\bigg\{\frac{1}{N}\sum_{i=1}^N\bigg(x_{i,T+1}'\Big[ \frac{1}{N}\sum_{j=1}^N Q_j\Big]^{-1} \frac{1}{N}\sum_{j =1}^N Q_j[\sqrt{T}(\beta_j-\beta_i)] \bigg)^2\bigg\}.
\end{align*}
Now, consider the right side: Under Assumption \ref{ass_3} 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{lem_2}).
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{ass_3} 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 \cite[see the discussion in][] {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{ass_3} 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.
\smallskip
\textbf{Statistical inference} We now develop an inference method for the difference of prediction errors. As a first step, we propose the estimator
\begin{align}\label{e:Ehat}
\hat E := & \frac{1}{N}\sum_{i=1}^N\bigg(x_{i,T+1}'\Big( \sum_{j=1}^N X_j' X_j \Big)^{-1} \sum_{j =1}^N X_j' X_j (\hat \beta_j-\hat \beta_i)\bigg)^2.
\end{align}
On the first glance, $\hat E$ may seem like a simple plug-in estimator for $E_2$ (see eq. \eqref{e:E_2}). 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$.
\begin{lem} \label{lem_3}
Suppose that Assumptions \ref{ass_1} and \ref{ass_2} hold. Then
\[
|\hat E-(E_1+E_2)|= \mathcal{O}^{|\mathbf{X}}_P(T^{-1}N^{-1/2}).\qquad
\]
\end{lem}
\noindent This lemma is a consequence of the much more precise analysis in the proof of Theorem \ref{theo_main}, 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
\begin{equation} \label{def_hat_Sigma_i}
\hat \Sigma^{(i,k)}(b) := (\hat \xi^{(i,k)}(|s-t|) \mathbb{I}\{|s-t| < b\})_{1\le s,t\le T },
\end{equation}
where the entry
\begin{eqnarray} \label{def_hat_xi}
\hat \xi^{(i,k)}(h) &:=& \frac{(y_i-X_i\hat \beta_i )_{1:T-h}' (y_k-X_k\hat \beta_k)_{h+1:T}}{T -h-K} \\
\nonumber &=& \frac{1}{T -h-K} \sum_{t=1}^{T -h} (y_i-X_i\hat \beta_i )_t (y_k-X_k\hat \beta_k )_{t+h}
\end{eqnarray}
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
\begin{equation} \label{e:hE1}
\hat E_1:= \frac{1}{N}\sum_{i=1}^NTr\Big[
\hat \Sigma^{(i,i)}( b) X_i [X_i' X_i]^{-1}x_{i,T+1}x_{i,T+1}'[X_i' X_i]^{-1}X_i' \Big ] =: \frac{1}{N}\sum_{i=1}^N \hat E_1^{(i)}.
\end{equation}
Notice that $ \hat E_1$ is the direct, empirical analogue to $E_1$ defined in \eqref{e:E_1}.
\begin{lem}\label{lem_4}
Under Assumptions \ref{ass_1}, \ref{ass_2} and \ref{ass_3} it holds that
\[
|\hat E_1-E_1|=\mathcal{O}^{|\mathbf{X}}_P\Big(\frac{b}{T^2}+\frac{ b^{7/4}}{\sqrt{N}T^{3/2}}\Big).
\]
If we choose $b=T^\rho$ with
\[
0 <\rho < \min\Big(1-\frac{\eta}{2},\frac{2}{7} \Big)
\]
the right side is of order $o^{|\mathbf{X}}_P(1/(\sqrt{N}T))$.
\end{lem}
\noindent 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
\begin{align*}
\Lambda := & \Big( \sum_{j=1}^N \frac{X_j' X_j}{NT} \Big)^{-1} \Big\{ \sum_{i=1}^N \frac{x_{i,T+1}x_{i,T+1}'}{N}\Big( \sum_{j=1}^N \frac{X_j' X_j}{NT} \Big)^{-1} \sum_{j =1}^N \frac{X_j' X_j (\beta_j-\beta_i)}{N\sqrt{T}}\Big\},\\
\Lambda_k := & \Big(\frac{X_k'X_k}{T}\Big)^{-1} x_{k,T+1}x_{k,T+1}'\Big( \sum_{j=1}^N \frac{X_j' X_j}{NT} \Big)^{-1} \sum_{j =1}^N \frac{X_j' X_j (\beta_j-\beta_k)}{N \sqrt{T}}.
\end{align*}
Therewith, we define the conditional asymptotic variance
\begin{align} \label{e:deftau}
\tau_N^2 = & \frac{1}{N}\sum_{i,k=1}^N \bigg\{2(\Sigma_N)_{i,k}^2 (x_{i,T+1}'(X_i'X_i/T)^{-1}[X_i'\Sigma_T X_k'/T](X_k'X_k/T)^{-1}x_{k,T+1})^2 \\
&\qquad \qquad \quad+4 ( \Lambda+ \Lambda_i)' \frac{X_i' (\Sigma_N)_{i,k} \Sigma_T X_k}{T} ( \Lambda+ \Lambda_k)\bigg\} \nonumber
\end{align}
As we show in the proof of Theorem \ref{theo_main}, $\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.
\begin{assumption} \label{ass_var}$ $
\begin{itemize}
\item[] (Non-degenerate variance) The random variable $\tau_N$ satisfies the two conditions
\[
\tau_N = \mathcal{O}_P(1), \qquad \tau_N^{-1} = \mathcal{O}_P(1).
\]
\end{itemize}
\end{assumption}
\noindent Let us define the standardized estimator
\begin{align} \label{e:ebreve}
\breve E:= \frac{\sqrt{N}T}{\tau_N}\big\{(\hat E-2\hat E_1)-(E^{pool}-E^{ind})\big\} \sim F_{\breve E},
\end{align}
where $F_{\breve E}$ refers to the cumulative distribution function of $\breve E$, conditional on $\mathbf{X}$ and $\tau_N$ is defined in \eqref{e:deftau}.
\begin{theo} \label{theo_main}
Suppose that Assumptions \ref{ass_1}, \ref{ass_2}, \ref{ass_3} and \ref{ass_var} hold. Denoting by $\Phi $ the cumulative distribution function of a standard normal distribution, we have
\[
\mathbb{E} \Big [ \sup_{z \in \mathbb{R}} |F_{\breve E}(z)-\Phi(z)| \big ] =o(1).
\]
\end{theo}
\noindent The proof of Theorem \ref{theo_main} 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.
\smallskip
\textbf{A conditional confidence interval} Theorem \ref{theo_main} 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
\begin{equation}\label{e:defCI}
\mathcal{C}_{1-\alpha}:= \Big[ \frac{\tau_N\Phi^{-1}(\alpha)}{\sqrt{N}T} +\hat E-2\hat E_1, \frac{\tau_N\Phi^{-1}(1-\alpha)}{\sqrt{N}T}+\hat E-2\hat E_1\Big].
\end{equation}
Notice that $\mathcal{C}_{1-\alpha}$ is an approximate $1-\alpha$ confidence interval in the sense that
\begin{align} \label{e:conf:guar}
\mathbb{P}\big((E^{pool}-E^{ind})\in \mathcal{C}_{1-\alpha}|\mathbf{X}\big) \overset{\mathbb{P}}{\to} 1-\alpha.
\end{align}
\begin{rem}{\rm (Variance estimation)
To conduct statistical inference it is necessary to estimate the long-run variance $\tau_N^2$. In practice, an estimator might use additional prior knowledge available to the user, such as independence of the individuals across $i$ or the like. Yet, it is possible to construct fully non-parametric estimators for $\tau_N^2$. For this purpose, define $\hat \Lambda, \hat \Lambda_k$ as the sample versions of $\Lambda, \Lambda_k$, where the slopes $\beta_i$ have been replaced by the OLS estimators $\hat \beta_i$. Moreover, recall the covariance estimator $\hat \Sigma^{(i,k)}$. Finally, let $\mathcal{K}$ be a continuous, symmetric, non-negative function, such that $\mathcal{K}(x)$ is non-increasing for $x\ge 0$. We assume that $\mathcal{K}(0)=1$ and $\mathcal{K}(1)=0$ and let $b'>0$ be a bandwidth choice. Then, we define
\begin{align}\label{e:def-tauhat}
\hat \tau_N^2:= & \frac{1}{N }\sum_{i,k=1}^N \mathcal{K}\bigg(\frac{|i-k|}{b'}\bigg)\bigg\{2 \Big( (x_{i,T+1}'(X_i'X_i/T)^{-1}\big(X_{i}'\hat \Sigma^{(i,k)}(b) X_k)_{s,t}/T\big)(X_k'X_k/T)^{-1}x_{k^,T+1}) \Big)^2 \notag\\
&\qquad \qquad \quad+4 (\hat \Lambda+ \hat \Lambda_i)' \frac{X_i' \hat \Sigma^{(i,k)}(b) X_k}{T} (\hat \Lambda+ \hat \Lambda_k)\bigg\}.
\end{align}
Two bandwidth parameters have to be chosen in this estimator, namely $b, b'$, which quantify dependence across time and individuals respectively. Stronger dependence requires larger values of the bandwidth choices. In the case of temporal independence, one can simply choose $b=1$ and in the case of spatial independence $b'=1$.}
\end{rem}
\begin{rem}\label{rem:compa}{\rm (Relation to existing work) We discuss certain aspects of the contribution, compared to \cite{pesaran:pick:timmermann:2022}. We have already mentioned that work in our Introduction, where we have pointed out its pivotal role of defining the forecasting error as a pooling criterion. We now focus on more technical aspects of that work.
\begin{itemize}
\item[i)] \cite{pesaran:pick:timmermann:2022} developed their procedures in a different model from ours: slopes are random and a procedure is developed for fixed $T$ and as $N \to \infty$. The proposed approach is technically speaking an oracle procedure and $T \to \infty$ is needed to estimate variances for practical inference. In our paper, by contrast, we have modeled slopes as fixed. The reason is that a random model is not helpful when assessing the forecasting error for large $N,T$ scenarios. The forecasting error can be estimated at a precision of $1/\sqrt{NT}$, while the average of $N$ slopes fluctuates at a larger order of magnitude $1/\sqrt{N}$. So, with random slopes and large $N,T$, one either has to condition on the slopes (rendering them again fixed) or one has to see an inference procedure dominated not by the model errors, but by the "random variation of the slopes" which may be undesirable. For this reason, in our large $N,T$ panel setup, we believe that fixed slopes are the preferable model. We also notice that fixing slopes avoids imposing assumptions on the joint slope distribution, which may be restrictive in practice.
\item[ii)] The inference tools developed in this paper and that of \cite{pesaran:pick:timmermann:2022} are different - we develop a confidence interval for the difference of forecasting errors, while they develop a hypothesis test. It turns out that these results are not equivalent, and we have made a confidence interval specifically to avoid challenges associated with conditional hypothesis testing. To be precise, \cite{pesaran:pick:timmermann:2022} develop a test for the point hypothesis of equivalent forecasting errors (eq. (19) of that work). This hypothesis is strictly speaking random and under standard assumptions only holds with probability $0$. Even when interpreted in a one-sided sense, it seems difficult to formulate standard asymptotic criteria for statistical tests, such as an asymptotic error rate or consistency. In contrast, we develop an asymptotic confidence interval that, in a clearer way, holds the asymptotic level $1-\alpha$ (see eq. \eqref{e:conf:guar} above).
\item[iii)] On a more basic level, theory in \cite{pesaran:pick:timmermann:2022} is developed for iid Gaussian model errors, while our approach can accommodate temporal and spatial dependence, and only imposes some weak moment assumptions on the error distributions.
\end{itemize}
}
\end{rem}
\section{Finite sample properties}\label{section-simulations}
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 \eqref{e:defCI} is replaced by the square root of $\hat{\tau}^2_N$ defined in \eqref{e:def-tauhat}. 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{sim:homogeneous}). 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{sim:heterogeneous}). This is followed by results for homogeneous and heterogeneous slopes under serial correlation in the model errors (Tables \ref{sim:homogeneous-ar1phi-03} and \ref{sim:heterogeneous-ar1-phi03}), 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{lem_4}, 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{sim:homogenoeus-fixedeffect} and \ref{sim:heterogeneous-fixed}), 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{sim:heteroskedastic-homogeneous} - \ref{sim:heteroskedastic-heterogeneous-parametric}). 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{sim:heterogeneous-fixedeffect-random}). We find that our CI performs well, which is not surprising, as our theory does not impose any distributional assumptions on the slope coefficients.
\subsection{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{sim:homogeneous} 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 \eqref{e:deftau} whereas $\hat \Lambda$ and $\hat \Lambda_i$ are non-zero in \eqref{e:def-tauhat}. 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{sim:heterogeneous}, 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{sim:homogeneous-ar1phi-03} and \ref{sim:heterogeneous-ar1-phi03}), 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{sim:homogeneous-ar1phi-03} and \ref{sim:heterogeneous-ar1-phi03} 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{sim:homogeneous-ar1phi-05} and \ref{sim:heterogeneous-ar1-phi05}) 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{sim:homogeneous-parametric-ar1phi-05} and \ref{sim:heterogeneous-parametric-ar1-phi05}).
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{sim:homogenoeus-fixedeffect} and \ref{sim:heterogeneous-fixed} resemble the ones in \ref{sim:homogeneous-ar1phi-03} and \ref{sim:heterogeneous-ar1-phi03}. 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{sim:heterogeneous-fixedeffect-random}, 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$.
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.9922 &1.6420 &0.9498 &0.9006 && 0.9890 &1.2841 &0.9530 &0.7388\\
15& 0.9972 &0.8861 &0.9572 &0.4983 && 0.9968 &0.3134 &0.9510 &0.1755\\
20& 0.9980 &0.4761 &0.9566 &0.2701 && 0.9982 &0.1905 &0.9502 &0.1075\\
25& 0.9988 &0.3169 &0.9580 &0.1802 && 0.9974 &0.1454 &0.9518 &0.0826\\
30& 0.9982 &0.2660 &0.9514 &0.1517 && 0.9992 &0.1084 &0.9488 &0.0617\\
40& 0.9990 &0.1820 &0.9518 &0.1038 && 0.9986 &0.0741 &0.9526 &0.0424\\
60& 0.9990 &0.1046 &0.9538 &0.0600 && 0.9986 &0.0482 &0.9480 &0.0276\\
80& 0.9990 &0.0794 &0.9474 &0.0456 && 0.9988 &0.0339 &0.9512 &0.0195\\
\bottomrule
\end{tabular}
\caption{Coverage rates and length of feasible and infeasible confidence interval with homogeneous slope coefficients.}
\label{sim:homogeneous}
\end{table}
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.9652 &3.0652 &0.9536 &2.7388 && 0.9732 &1.7765 &0.9524 &1.4231\\
15& 0.9608 &2.1541 &0.9542 &2.0317 && 0.9604 &0.8114 &0.9512 &0.7689\\
20& 0.9660 &1.5388 &0.9598 &1.4877 && 0.9552 &0.6446 &0.9466 &0.6250\\
25& 0.9564 &1.3201 &0.9554 &1.2948 && 0.9556 &0.5575 &0.9494 &0.5444\\
30& 0.9598 &1.1893 &0.9568 &1.1693 && 0.9528 &0.4808 &0.9494 &0.4725\\
40& 0.9596 &0.9772 &0.9568 &0.9668 && 0.9556 &0.3978 &0.9536 &0.3931\\
60& 0.9602 &0.7370 &0.9590 &0.7324 && 0.9514 &0.3186 &0.9490 &0.3161\\
80& 0.9538 &0.6435 &0.9528 &0.6403 && 0.9526 &0.2699 &0.9486 &0.2684\\
\bottomrule
\end{tabular}
\caption{Coverage rates and length of feasible and infeasible confidence interval with heterogeneous slope coefficients.}
\label{sim:heterogeneous}
\end{table}
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.9760 &2.8087 &0.9506 &1.6053 && 0.9298 &1.5023 &0.9510 &0.8435\\
15& 0.9792 &1.1745 &0.9562 &0.6975 && 0.8740 &0.5415 &0.9512 &0.3168\\
20& 0.9910 &0.7589 &0.9600 &0.4471 && 0.9102 &0.3393 &0.9554 &0.1998\\
25& 0.9942 &0.5597 &0.9622 &0.3298 && 0.9424 &0.2485 &0.9514 &0.1460\\
30& 0.9958 &0.4463 &0.9576 &0.2627 && 0.9594 &0.1961 &0.9522 &0.1152\\
40& 0.9980 &0.3149 &0.9574 &0.1852 && 0.9782 &0.1390 &0.9502 &0.0816\\
60& 0.9994 &0.1988 &0.9574 &0.1167 && 0.9918 &0.0880 &0.9536 &0.0515\\
80& 0.9992 &0.1444 &0.9538 &0.0843 && 0.9962 &0.0647 &0.9526 &0.0378\\
\bottomrule
\end{tabular}
\caption{ Coverage rates and length of feasible and infeasible confidence interval with homogeneous slope coefficients and AR(1) model errors with $\phi=0.3$.}
\label{sim:homogeneous-ar1phi-03}
\end{table}
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.9592 &4.1472 &0.9494 &3.5791 && 0.9366 &2.0391 &0.9480 &1.6667\\
15& 0.9494 &2.4120 &0.9516 &2.3537 && 0.9184 &1.0872 &0.9482 &1.0488\\
20& 0.9512 &1.9234 &0.9536 &1.9169 && 0.9332 &0.8457 &0.9500 &0.8415\\
25& 0.9494 &1.6337 &0.9504 &1.6393 && 0.9334 &0.7218 &0.9494 &0.7243\\
30& 0.9534 &1.4501 &0.9560 &1.4603 && 0.9392 &0.6382 &0.9492 &0.6427\\
40& 0.9498 &1.2086 &0.9548 &1.2217 && 0.9422 &0.5380 &0.9508 &0.5434\\
60& 0.9486 &0.9627 &0.9522 &0.9749 && 0.9454 &0.4291 &0.9482 &0.4343\\
80& 0.9554 &0.8286 &0.9586 &0.8364 && 0.9466 &0.3674 &0.9516 &0.3709\\
\bottomrule
\end{tabular}
\caption{Coverage rates and length of feasible and infeasible confidence interval with heterogeneous slope coefficients and AR(1) model errors with $\phi=0.3$.}
\label{sim:heterogeneous-ar1-phi03}
\end{table}
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.8386 &5.0807 &0.9562 &3.3365 && 0.3550 &2.4451 &0.9464 &1.5662\\
15& 0.9790 &1.4493 &0.9620 &0.8762 && 0.7832 &0.6764 &0.9480 &0.4065\\
20& 0.9910 &0.9164 &0.9522 &0.5432 && 0.9422 &0.4005 &0.9576 &0.2373\\
25& 0.9944 &0.6681 &0.9534 &0.3937 && 0.9742 &0.2821 &0.9508 &0.1658\\
30& 0.9988 &0.4909 &0.9570 &0.2872 && 0.9872 &0.2171 &0.9528 &0.1271\\
40& 0.9990 &0.3455 &0.9568 &0.2015 && 0.9940 &0.1513 &0.9496 &0.0881\\
60& 0.9990 &0.2076 &0.9596 &0.1209 && 0.9976 &0.0941 &0.9464 &0.0546\\
80& 0.9996 &0.1485 &0.9482 &0.0862 && 0.9992 &0.0678 &0.9564 &0.0393\\
\bottomrule
\end{tabular}
\caption{Coverage rates and length of feasible and infeasible confidence interval with homogeneous slopes and individual fixed effects. }
\label{sim:homogenoeus-fixedeffect}
\end{table}
\begin{table}[H]\centering
\begin{tabular}{@{}lrrrrcrrrr@{}}\toprule
& \multicolumn{4}{c}{$N=100$} & \phantom{abc}& \multicolumn{4}{c}{$N=500$} \\
\cmidrule{2-5} \cmidrule{7-10}
$T$& $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$ && $\hat{\mathcal{C}}_{0.95}$ & $L(\hat{\mathcal{C}}_{0.95})$ & $\mathcal{C}_{0.95}^*$ &$L(\mathcal{C}_{0.95}^*)$\\
\midrule
10& 0.8894 &8.5260 &0.9544 &7.7738 && 0.6284 &4.8723 &0.9494 &3.9320\\
15& 0.9522 &4.3905 &0.9540 &4.1887 && 0.9002 &1.9939 &0.9504 &1.8926\\
20& 0.9510 &3.1799 &0.9496 &3.0716 && 0.9446 &1.4991 &0.9512 &1.4498\\
25& 0.9560 &2.8129 &0.9520 &2.7550 && 0.9514 &1.2308 &0.9538 &1.2063\\
30& 0.9556 &2.3995 &0.9532 &2.3718 && 0.9500 &1.0865 &0.9468 &1.0696\\
40& 0.9562 &1.9759 &0.9558 &1.9589 && 0.9482 &0.8913 &0.9474 &0.8818\\
60& 0.9518 &1.5410 &0.9524 &1.5321 && 0.9504 &0.7028 &0.9488 &0.6984\\
80& 0.9492 &1.3044 &0.9490 &1.3010 && 0.9512 &0.5944 &0.9518 &0.5918\\
\bottomrule
\end{tabular}
\caption{Coverage rates and length of feasible and infeasible confidence interval with slope heterogeneity and individual fixed effects.}
\label{sim:heterogeneous-fixed}
\end{table}
\section{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.
\clearpage
\bibliographystyle{chicago}
\bibliography{ref_separability-supnorm}
\newpage