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.
70,712 characters
Improved Bootstrap Inference for Dynamic Panel Data models with Interactive Effects
{
\title{\bf Improved Bootstrap Inference for Dynamic Panel Data models with Interactive Effects}
\author[1]{Artūras Juodis }
\author[2]{Ovidijus Stauskas\thanks{Corresponding author. Financial support from the Dutch Research Council (NWO) under research grant VI.Vidi.231E.030 (Juodis) is gratefully acknowledged by all authors.}}
\author[1]{Sander Tromp}
\affil[1]{Amsterdam School of Economics, University of Amsterdam and Tinbergen Institute}
\affil[2]{BI Norwegian Business School, Department of Economics}
\maketitle
}
\begin{abstract}
We study recursive-design wild bootstrap inference for dynamic panel data models with unobserved common factors estimated by Common Correlated Effects. In the large $N,T$ setting, the bootstrap reproduces the biased limiting distribution in pure autoregressive models, but fails to capture all bias and factor-estimation variance components in models with additional regressors, particularly under weak exogeneity. We trace this failure to holding regressors fixed across bootstrap replications. We propose to combine bootstrap procedure with available bias-correction methods to conduct adjusted inference. Monte Carlo evidence shows substantial improvements over conventional strategies of using bias-correction paired with cross-sectional bootstrap methods.
\end{abstract}
\textbf{Keywords}: Common correlated effects, wild bootstrap, dynamic panel data models.
\section{Introduction}
Dynamic panel data models with unobserved common shocks are now standard in empirical macroeconomics and finance, where $N$ cross-sectional units are exposed to global factors and exhibit substantial levels of persistence. In such environments, the common correlated effects (CCE) estimator of \citet{Pesaran2006} has become a central tool for dealing with cross-sectional dependence by augmenting the model with cross-sectional averages of the observed variables. Subsequent work has extended and refined this approach in multiple directions; see, e.g., \citet{JUODIS2026106120} for a recent review.
Despite these advances, valid inference in dynamic CCE panels remains challenging, as emphasized in Lesson 2 of \citet{JUODIS2026106120}. When lagged dependent variables and individual effects are present, the CCE pooled (CCEP) estimator inherits the well-known ``Nickell bias'', whose magnitude can be amplified by the presence of unobserved factors; see \citet{de2021bias} and \citet{juodis2021robustness}. In addition to the ``Nickell bias'', which is inversely proportional to the time-series dimension $T$, the CCEP estimator also suffers from a cross-sectional-averages-induced bias, that is inversely related to the cross-sectional dimension $N$; see e.g., \citet{Westerlund2015}, \citet{juodis2021robustness}, and \citet{de2024cross}.
Due to the presence of the ``Nickell bias'', the CCE method becomes inconsistent for a small number of time-series observations, in the so-called fixed-$T$ regime. In many situations, analytical bias correction, similar to that of \citet{moon2017dynamic}, or the half-panel jackknife (HPJ) of \citet{dhaene2015split}, paired with the cross-sectional (pairs) bootstrap, is used. These methods are based on an asymptotic approximation in which both $N$ and $T$ are large, subject to some rate restrictions.
In response to the unsatisfactory finite-sample properties of such large-$N$, large-$T$ bias-correction approaches, \citet{de2021bias} proposed a modified CCEP-based methodology in which only $N$ is large, while $T$ can be fixed. This approach builds upon the insights of, among others, \citet{Bun2005} and \citet{Dhaene}, that, for the panel AR(1) model with strictly exogenous regressors, a fixed-$T$ bias-corrected estimator can be constructed using the formula for the corresponding ``Nickell bias''. While the approach works well within the class of models studied by \citet{de2021bias}, it is generally cumbersome to extend to settings with higher-order dynamics, unrestricted forms of heteroskedasticity, and unrestricted weakly exogenous regressors.
In this paper, we use the finite-sample evidence in \citet{de2021bias} as motivation, but instead leverage the recent panel-data literature that highlights the benefits of using bootstrap-based methods; see, e.g., \citet{Higgins2024Bootstrap,HigginsJochmans2025}. Given the dynamic nature of the autoregressive model, this paper focuses on the properties of the recursive-design wild bootstrap (RD-WB) of \citet{gonccalves2015bootstrap} in the context of dynamic panel data models estimated using the CCEP estimator. We consider both a pure panel AR(1) specification and a more general ARX(1) model with additional regressors. To the best of our knowledge, this is the first paper in the literature that focuses on these settings.
Within a double-asymptotic framework with large $N$ and $T$, we investigate conditions under which the procedure of combining dynamic CCEP with RD-WB yields valid inference for CCEP estimators. Our first main result shows that, in the AR(1) case, the RD-WB fully replicates the asymptotic distribution of the CCEP estimator: the bootstrap distribution centers at the same biased value as the original estimator, mirroring the behavior documented by \citet{gonccalves2015bootstrap} for the fixed effects estimator. In contrast, for ARX(1) models, the RD-WB neither fully replicates the bias nor the variances originating from the factor estimates; see \citet{Juodis2022CCER}. These problems are generally amplified when some of the regressors are only weakly exogenous, as the ``Nickell bias'' is then also not fully replicated.
The failure of the bootstrap in this general design can be fully attributed to the fact that additional regressors are kept fixed throughout the bootstrap replications, similar to the fixed-design bootstrap studied by \citet{gonccalves2015bootstrap}. Hence, the bootstrap implementation shuts down any channels through which these regressors can affect the first-order asymptotic distribution of the CCEP estimator. In the Supplementary Online Appendix, we show how discrepancies between the asymptotic and bootstrap distributions diminish if some structure is imposed on the additional regressors, opening up channels for bootstrap variation originating from these regressors.
Given the shortcomings of the RD-WB method in the general setting, we show how the bootstrap method can be combined with the available statistical CCEP toolkit. The proposed solution involves two main ingredients: (i) appropriate studentization that accounts for variance discrepancies between the true and bootstrap distributions; and (ii) the use of bias-corrected estimators, either analytical or jackknife-based, in the bootstrap recursions. Monte Carlo results, presented in the Supplementary Online Appendix, show that this combined approach substantially improves upon the naive approach based on bias correction with the cross-sectional (pairs) bootstrap. Improvements are especially stark for designs with small to moderate values of $T$.
The remainder of the paper is organized as follows. Section \ref{section::model} introduces the model and the pooled CCE estimator. Section \ref{section::bootstrap} summarizes the main theoretical results for both the AR(1) and ARX(1) models. Section \ref{section::takeaways} translates these theoretical results into takeaways for empirical researchers. Section \ref{section::empirical} provides an empirical illustration. Section \ref{section::conclusions} concludes.
\section{Dynamic linear regression with common factors}
\setcounter{equation}{0}
\label{section::model}
\subsection{The model}
We consider the following first-order dynamic panel data model with regressors and unobserved interactive effects (common factors)
\begin{equation}
\label{eq::model_main}
y_{i,t}=\alpha y_{i,t-1}+\+\beta'\*\mathbf{x}_{i,t}+\+\gamma_{y,i}'\+\mathbf{f}_{t}+\varepsilon_{i,t},
\end{equation}
for $i=1,\ldots,N$ cross-sectional units and $t=1,\ldots,T$ time periods. Here $y_{i,t}$ is the observed target variable, $\+\mathbf{x}_{i,t}$ is a $[k\times 1]$ vector of observed regressors (covariates), $\+\gamma_{y,i}'\+\mathbf{f}_{t}$ is the unobserved common factor, while $\varepsilon_{i,t}$ is the idiosyncratic error term. For simplicity, we assume that initial observations $(y_{i,0},i=1,\ldots,N)$ are observed. As standard in the literature, $\+\gamma_{y,i}$ is an $[R\times 1]$ vector of factor loadings, while $\+\mathbf{f}_{t}$ is the corresponding vector of common factors. Neither the factor loadings, nor the corresponding factors are observed, constituting the interactive fixed effects model.\footnote{For the sake of notational simplicity, we omit the individual effects $\eta_{i}$ from the model in Eq. \eqref{eq::model_main}. All theoretical results in this paper continue to hold in the presence of individual-specific intercepts and/or individual-specific deterministic trends.}
In what follows, the main object of interest is the full vector of common coefficient $\+\delta:=[\alpha,\+\beta']'$. Not accounting for unobserved factors results in inconsistent estimates for $\+\delta$ when the omitted factors are correlated with the included regressors, irrespective of the length of the time-series.
In the literature multiple estimation methods for $\+\delta$ have been proposed. The suggested approaches mainly differ in whether $T$ is assumed to be small (the so-called fixed-$T$ regime), or that $T$ is allowed to be large and of similar magnitude to $N$ (the cross-sectional dimension). Contributions of \citet{Holtz-Eakin1988}, \citet{Ahn2013}, \citet{Robertson2015}, and \citet{JuodisSarafidis2019} (among others) fall into the fixed-$T$ category. The methods use the Generalized Method of Moments (GMM) methodology similar to the seminal work of \citet{Arellano1991}; see, e.g., \citet{JuodisSarafidis2016ER} for a detailed review.
The large $T$ methods, on the other hand, use some least-squares objective function to estimate $\+\delta$ jointly with $\+\gamma_{y,i}$ and $\+\mathbf{f}_{t}$. The Interactive Fixed Effects (or the Principal components) estimator of \citet{bai2009panel}, and \citet{moon2015linear}, and the Common Correlated Effects (CCE) estimator of \citet{Pesaran2006} are the two prominent estimators of this class. Due to the dynamic nature of the model, these estimators generally suffer from the weak-exogeneity (or the \citealp{Nickell1981}) bias. As a result, some form of bias-correction is required to guarantee asymptotic validity of the corresponding inference procedures; see, e.g., \citet{moon2017dynamic} and \citet{juodis2021robustness}.
In this paper, we contribute to the large $T$ literature and focus on bootstrap-based inference for the CCE estimator of \citet{Pesaran2006}.
\subsection{The Common Correlated Effects (CCE) estimator}
In this section, we introduce the pooled CCE (CCEP) estimator of \citet{Pesaran2006} and discuss the theoretical results available to this estimator in the context of the model in Eq. \eqref{eq::model_main}.
Let $\+\mathbf{w}_{i,t}:=[y_{i,t-1},\+\mathbf{x}_{i,t}']'$, such that the original model can be expressed as
\begin{equation}
y_{i,t}=\+\delta'\*\mathbf{w}_{i,t}+\+\gamma_{y,i}'\+\mathbf{f}_{t}+\varepsilon_{i,t},
\end{equation}
or using the stacked matrix notation as
\begin{equation}
\+\mathbf{y}_{i}=\*\mathbf{W}_{i}\+\delta+\+\mathbf{F} \+\gamma_{y,i}+\+\varepsilon_{i}.
\end{equation}
Here $\+\mathbf{y}_{i}$ is a $[T\times 1]$ vector of time-series stacked observations for the unit $i$, and similarly for all other quantities.
If all elements of $\+\mathbf{F}$ were known, the natural estimator for $\+\delta$ was the fixed effects (within group) type estimator that is defined as the joint minimizer of the least squares objective function together with the factor loadings $(\+\gamma_{y,1},\ldots,\+\gamma_{y,N})$. By the Frisch-Waugh-Lowell theorem the estimator $\widehat{\+\delta}_{LS}$ is given by
\begin{equation}
\label{eq::LS_infeasbile}
\widehat{\+\delta}_{LS}=\left(\frac{1}{NT}\sum_{i=1}^{N}\+\mathbf{W}_{i}'\+\mathbf{M}_{\+\mathbf{F}}\+\mathbf{W}_{i}\right)^{-1}\left(\frac{1}{NT}\sum_{i=1}^{N}\+\mathbf{W}_{i}'\+\mathbf{M}_{\+\mathbf{F}}\+\mathbf{y}_{i}\right),
\end{equation}
where $\+\mathbf{M}_{\+\mathbf{F}}:= \+\mathbf{I}_{T} - \+\mathbf{F} \left(\+\mathbf{F}'\+\mathbf{F}\right)^{-1} \+\mathbf{F}'$. Unfortunately, not all columns of $\+\mathbf{F}$ are known; therefore, they need to be estimated from the data.
For this, \citet{Pesaran2006} suggested estimating $\+\mathbf{F}$ by assuming that
\begin{equation}
\label{eq::model_w}
\+\mathbf{W}_{i}=\+\mathbf{F} \+\Gamma_{w,i}+\+\mathbf{U}_{i},
\end{equation}
where $\+\mathbf{U}_{i}$ is the idiosyncratic error vector independent (or at least uncorrelated) with the factor component $\+\mathbf{F} \+\gamma_{y,i}$. The corresponding estimator $\widehat{\+\mathbf{F}}$ is formed by taking cross-sectional averages of all observed quantities (that are linear in common factors), i.e.,
\begin{equation*}
\widehat{\+\mathbf{F}}=[\overline{\+\mathbf{y}},\overline{\+\mathbf{W}}].
\end{equation*}
The resulting Common Correlated Effects Pooled (CCEP) estimator is just a feasible (plug-in) version of Eq. \eqref{eq::LS_infeasbile}
\begin{equation}
\label{eq::LS_feasbile}
\widehat{\+\delta}_{CCEP}=\left(\frac{1}{NT}\sum_{i=1}^{N}\+\mathbf{W}_{i}'\+\mathbf{M}_{\widehat{\+\mathbf{F}}}\+\mathbf{W}_{i}\right)^{-1}\left(\frac{1}{NT}\sum_{i=1}^{N}\+\mathbf{W}_{i}'\+\mathbf{M}_{\widehat{\+\mathbf{F}}}\+\mathbf{y}_{i}\right).
\end{equation}
As long as $\widehat{\+\mathbf{F}}$ is a consistent estimator of $\+\mathbf{F}$ (the requirement usually formalized in the form of the rank condition on cross-sectional averages of all factor loadings), the CCEP estimator can be shown to be consistent as $N,T\to\infty$; see, e.g., \citet{Pesaran2006} and \citet{juodis2021robustness}.
\subsection{CCE framework for dynamic models}
While appealing, the dynamic nature of the model in Eq. \eqref{eq::model_main}, is generally incompatible for all regressors $\+\mathbf{W}_{i}$ in Eq. \eqref{eq::model_w}. This was first recognized in \citet{Chudik2015}, \citet{de2021bias}, and later formalized in \citet{juodis2021robustness} and \citet{Juodis2022CCER}. In particular, in dynamic models the following three features might appear: i) only rotation of $\widehat{\+\mathbf{F}}$ identifies $\+\mathbf{F}$; ii) factors different to $\+\mathbf{F}$ (nuisance/distinct factors) can enter $\+\mathbf{W}_{i}$; iii) some of these \emph{distinct} factors are not consistently estimable by cross-sectional averages.
Below, we illustrate these three features using two stylized Data Generating Processes (DGPs). First, consider the panel AR(1) model
\begin{equation}
\label{eq::model_main_ar}
y_{i,t}=\alpha y_{i,t-1}+\+\gamma_{y,i}'\+\mathbf{f}_{t}+\varepsilon_{i,t},
\end{equation}
where $\varepsilon_{i,t} \sim i.i.d.(0,\sigma_{\varepsilon}^{2})$. In this model, the unobserved factor component $\+\gamma_{y,i}'\+\mathbf{f}_{t}$ can have at most one unobserved factor (thus, it collapses to a scalar $\gamma_{y,i} f_{t}$) for factors to be estimable by cross-sectional averages. For the true DGP we have
\begin{equation*}
(\overline{y}_{t}-\alpha_{0}\overline{y}_{t-1})\overline{\gamma}_{y}^{-1}=f_{t}+O_P(N^{-1/2}).
\end{equation*}
Hence, the linear combination of the two cross-sectional averages consistently estimates the true factor $f_{t}$ as long as $\overline{\gamma}_{y}$ is non-zero in the limit. This is the property $i)$ mentioned above.
Next, consider the cross-sectional average of the only regressor itself - $y_{i,t-1}$. It can be expanded
\begin{equation}
\label{eq::average_lagy_expansion_AR1}
\overline{y}_{t-1}=\overline{\gamma}_{y}g_{t-1} + O_P(N^{-1/2}).
\end{equation}
Here, assuming infinite initialization of the process $y_{i,t}$, $g_{t-1}=\sum_{j=0}^\infty (\alpha_0)^jf_{t-j-1}$, and so we can expand $g_{t}=\alpha_{0}g_{t-1}+f_{t}$. Evidently, $\overline{y}_{t-1}$ will not identify $f_{t}$ (the factor that enters the model in Eq. \eqref{eq::model_main_ar}). Instead, the new factor - $g_{t-1}$ - is identified. This factor is a nuisance factor, as long as identification of $f_{t}$ is considered. This is the property $ii)$ mentioned above.
The third property is best illustrated by adding an additional regressor $x_{i,t}$ into the above model
\begin{align}
\label{example:2eq1}
y_{i,t}&=\alpha y_{i,t-1}+\beta x_{i,t}+\+\gamma_{y,i}'\+\mathbf{f}_{t}+\varepsilon_{i,t},\\
x_{i,t}&=\+\gamma_{x,i}'\+\mathbf{f}_{t}+u_{i,t}.
\label{example:2eq2}
\end{align}
Here we assume that $\+\mathbf{f}_{t}$ now has two unobserved factors, while $u_{i,t}$ is possibly serially correlated sequence independent of $\varepsilon_{i,t}$ for $t=1,\ldots,T$. Using similar algebraical manipulations as above, it can be easily seen that both factors $\+\mathbf{f}_{t}$ are identified from cross-sectional averages as long as the rank of $[\overline{\+\gamma}_{y},\overline{\+\gamma}_{x}]$ is full in the limit as $N\to\infty$. Next, consider a cross-sectional of $y_{i,t-1}$. It admits expansion of the form
\begin{equation}
\label{eq::average_lagy_expansion_ARX}
\overline{y}_{t-1}=\overline{\+\gamma}_{y^{+}}'\+\mathbf{g}_{t-1} + O_P(N^{-1/2}),
\end{equation}
where (vector-valued) $\+\mathbf{g}_{t-1}$ is defined recursively as above and $\overline{\+\gamma}_{y^{+}}:=\overline{\+\gamma}_{y}+\beta_{0}\overline{\+\gamma}_{x}$. While Eq. \eqref{eq::average_lagy_expansion_ARX} looks similar to Eq. \eqref{eq::average_lagy_expansion_AR1}, the implications for the properties of the CCEP estimator are different: in the former case $\overline{y}_{t-1}$ is a scalar while $\+\mathbf{g}_{t-1}$ is a vector. As a result, $\overline{y}_{t-1}$ will not consistently estimate both factors, but only their linear combination - $\widetilde{g}_{t-1}$. This is the property $iii)$ mentioned above.
In particular, upon adding and subtracting $\widetilde{g}_{t-1}$
\begin{equation*}
\overline{y}_{t-1}=\underbrace{\widetilde{g}_{t-1}}_{I}+\underbrace{(\overline{\+\gamma}_{y^{+}}'\+\mathbf{g}_{t-1}-\tilde{g}_{t-1})}_{II} + O_P(N^{-1/2}).
\end{equation*}
Here $II$ is a mean zero term (at least asymptotically) and is generally of the same order as the remainder. Unlike the remainder term, $II$ is not cross-sectionally independent and can be correlated with both $\+\gamma_{y,i}'\+\mathbf{f}_{t}$ and $\+\gamma_{x,i}'\+\mathbf{f}_{t}$. As shown by \citet{de2021bias}, \citet{juodis2021robustness}, and \citet{Juodis2022CCER}, the presence of such additional (distinct/nuisance) factors has non-negligible effect on the asymptotic distribution of the CCEP estimator.
In particular, for the class of models considered in this paper, the CCEP estimator admits the asymptotic expansion of the from
\begin{equation}
\label{eq::asymptotic_distribution}
\sqrt{NT}(\widehat{\+\delta}_{CCEP}-\+\delta_{0})=\+\xi_{0}+\+\xi_{\perp}+\sqrt{\frac{N}{T}} \* b_{1}+\sqrt{\frac{T}{N}}\* b_{2}+\*r_{N,T},
\end{equation}
see also \citet{Juodis2022CCER}.
Here $\+\xi_{0}$ is asymptotically normal term driven by innovations $(\varepsilon_{i,t}, t=1,\ldots,T)$; $\+\xi_{\perp}$ is asymptotically normal term (uncorrelated with $\+\xi_{0}$), present when $II$ is non-negligible; $\+ b_{1}$ is the ``Nickell bias''; $\+ b_{2}$ is the factor-approximation bias (first derived in \citealp{Westerlund2015}); $\*r_{N,T}$ is an asymptotically negligible remainder term. Unless stated otherwise, all terms in Eq. \eqref{eq::asymptotic_distribution} are $O_P(1)$.
Hence, any asymptotically valid inference procedure should account for the two variance terms and the two bias terms present in Eq. \eqref{eq::asymptotic_distribution}. In the next section, we show how the Recursive Design Wild Bootstrap (RD-WB) of \citet{gonccalves2015bootstrap} can be adapted to our setting, and under which conditions it replicates the decomposition in Eq. \eqref{eq::asymptotic_distribution}.
\section{Recursive design wild bootstrap}
\setcounter{equation}{0}
\label{section::bootstrap}
\subsection{Implementation}
In the following, we summarize the main ingredients of the RD-WB procedure from \citet{gonccalves2015bootstrap} in the context of the CCEP estimator. Note that while innovations $\varepsilon_{i,t}$ are usually assumed to be serially uncorrelated, i.e., the model is dynamically complete, the dynamic structure of $\+\mathbf{x}_{i,t}$ is normally kept unrestricted. As a result, in what follows, we keep $\+\mathbf{x}_{i,t}$ fixed throughout all bootstrap replications.
Motivated by our empirical application, we assume that the factor estimates are given by
\begin{equation*}
\widehat{\+\mathbf{F}}=[\overline{\+\mathbf{y}},\overline{\+\mathbf{y}}_{-1},\overline{\+\mathbf{X}},\overline{\+\mathbf{X}}_{-1}].
\end{equation*}
The inclusion of $\overline{\+\mathbf{X}}_{-1}$ symmetrizes treatment of $y_{i,t}$ and $\+\mathbf{x}_{i,t}$, and aligns our suggested implementation to those in \citet{Chudik2015} and \citet{de2021bias}.
Let $\widehat{\+\delta}_{CCEP}$ be given as in Eq. \eqref{eq::LS_feasbile} then we set
\begin{align}
\label{eq::cce_fitted1}
\widehat{\+\gamma}_{y,i}&:= \left(\sum_{t=1}^{T}\widehat{\+\mathbf{f}}_{t}\widehat{\+\mathbf{f}}_{t}'\right)^{-1}\sum_{t=1}^{T}\widehat{\+\mathbf{f}}_{t}(y_{i,t}-\+\mathbf{w}_{i,t}'\widehat{\+\delta}_{CCEP}),\\
\widehat{\varepsilon}_{i,t}&:= y_{i,t}-\+\mathbf{w}_{i,t}'\widehat{\+\delta}_{CCEP}-\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t},
\label{eq::cce_fitted2}
\end{align}
where $\widehat{\+\mathbf{f}}_{t}$ is the $t$th column of $\widehat{\+\mathbf{F}}$.
\begin{algorithm}
\label{algo::naive}
\newcounter{bean}
\setcounter{bean}{0}
\begin{center}
\textnormal{
\begin{list}
{\textsc{Step} \arabic{bean}.}{\usecounter{bean}}
\item Obtain CCEP loadings and residuals as in Eqs. \eqref{eq::cce_fitted1}-\eqref{eq::cce_fitted2}.
\item Set $y_{i,0}^{*}=y_{i,0}$.
\item For $b=1,\ldots,B$ generate bootstrap dataset $y_{i,t}^{*}$ recursively for $t=1,\ldots,T$,
\begin{equation}
\label{eq::bootstrap_DGP}
y_{i,t}^{*}:=\widehat{\alpha}_{CCEP} y_{i,t-1}^{*}+\widehat{\+\beta}_{CCEP}'\+\mathbf{x}_{i,t}+\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t}+\omega_{i,t}^{*}\widehat{\varepsilon}_{i,t},
\end{equation}
for all $i=1,\ldots,N$, where $\omega_{i,t}^{*} \sim i.i.d.(0,1)$, $E^*[(\omega_{i,t}^{*})^4]$ finite.
\item Given $((y_{i,t}^{*},\+\mathbf{x}_{i,t}')',i=1,\ldots,N,t=1,\ldots, T)$ obtain $\widehat{\+\delta}_{CCEP}^{*}$ from Eq. \eqref{eq::LS_feasbile} using $\widehat{\+\mathbf{F}}^{*}=[\overline{\+\mathbf{y}}^{*},\overline{\+\mathbf{y}}^{*}_{-1},\overline{\+\mathbf{X}},\overline{\+\mathbf{X}}_{-1}]$
\item Repeat for $b=1,\ldots,B$.
\end{list}
}
\end{center}
\end{algorithm}
The choice in Step 2 can be different in practice due to large $T$. However, it simplifies the demonstration on how the bootstrap averages consistently estimate latent factors in the bootstrap realm (see Section \ref{ssection::discussion} for a precise argument). In practice, we implement Algorithm \ref{algo::naive} using $\omega_{i,t}^{*} \sim Rademacher(-1;1)$.
In what follows, we will refer to the bootstrap procedure described in Algorithm \ref{algo::naive} as ``naive'' bootstrap. This procedure does not impose any DGP-induced restrictions on cross-sectional averages of the data, see Section \ref{section:extensions} for further discussion.
\begin{remark}
\textnormal{As an alternative to re-estimating $\widehat{\+\mathbf{f}}_{t}^{*}$, we can alternatively consider the bootstrap scheme where $\widehat{\+\mathbf{f}}_{t}^{*}:=\widehat{\+\mathbf{f}}_{t}$; see e.g.,\citet{Westerlund2019}. Unfortunately, this implementation fails to replicate the variance term $\+\xi_{\perp}$, as well as the bias term $\* b_{2}$.}
\end{remark}
\subsection{Assumptions}
We assume that $\+\mathbf{x}_{i,t}$ are generated as follows
\begin{equation}
\label{eq::x}
\+\mathbf{x}_{i,t}=\+\theta_0 y_{i,t-1}+\+\Gamma_{x,i}'\*f_t+\+\nu_{i,t}.
\end{equation}
Such that the reduced form for $\+\mathbf{z}_{i,t}:=[y_{i,t},\+\mathbf{x}_{i,t}']'$ is of the form
\begin{align}
\label{z_VAR}
\+\mathbf{z}_{i,t}=\*A_0^{\prime}\+\mathbf{z}_{i,t-1}+\+\Gamma_i'\+\mathbf{f}_t+\+\mathbf{e}_{i,t},
\end{align}
where
\begin{align*}
&\*A_0:=\begin{bmatrix} \alpha_0+\+\beta_0^{\prime}\+\theta_0 & \+\theta_0^{\prime}\\
\*0_{k\times 1} & \*0_{k\times k}\end{bmatrix},\\
&\+\Gamma_i:=[\+\Gamma_{x,i}\+\beta_0+\+\gamma_{y,i}, \+\Gamma_{x,i}],\\
&\*e_{i,t}:=[\+\nu_{i,t}'\+\beta_0+\varepsilon_{i,t}, \+\nu_{i,t}']'.
\end{align*}
Below we summarize the set of assumptions used in the remainder of the paper. To fix the notation we set $K:=k+1$, and $0<\Delta<\infty$ is some arbitrary finite constant.
\begin{assumption}
\label{ass::1}(a) The error term $\varepsilon_{i,t}$ is i.i.d. over $i,t$ with $E[\varepsilon_{i,t}]=0$, $E[\varepsilon_{i,t}]=\sigma^2$ and $E[|\varepsilon_{i,t}|^{8}]<\Delta$; (b) The error terms $(\+\nu_t,t=1,\ldots,T)$ are covariance stationary for all $i=1,\ldots,N$; (c) $E[\+\nu_{i,t}]
=\*0_{k\times 1}$, $E[\+\nu_{i,t}\+\nu_{i,t}']=\+\Sigma_{\+\nu}(0)$, $E\left[\left\|\+\nu_{i,t} \right\|^{8}\right]<\Delta$; (d) The sequence $E\left[\+\nu_{i,h}\+\nu'_{i,0}\right]=\+\Sigma_{\+\nu}(h)$ is absolutely summable.
\end{assumption}
\begin{assumption}
\label{ass::2}(a) The factors $(\+\mathbf{f}_t,t=1,\ldots,T)$ are covariance stationary; (b) $E[\+\mathbf{f}_t\+\mathbf{f}_t']=\+\Sigma_{\+\mathbf{f}}(0)\in \mathbb{R}^{R\times R}$ is positive definite, and $E[\left\|\+\mathbf{f}_t \right\|^{8} ]<\Delta$; (c) The sequence $E[\+\mathbf{f}_{h}\+\mathbf{f}_{0}']=\+\Sigma_{\+\mathbf{f}}(h)$ is absolutely summable.
\end{assumption}
\begin{assumption}
\label{ass::3}(a) The loadings $\+\Gamma_{i}\in \mathbb{R}^{R\times K}$; (b) $\mathrm{rk}(\overline{\+\Gamma})=R=K$ as $N\to \infty$; (c) $\frac{1}{N}\sum_{i=1}^N\left\|\+\Gamma_i \right\|^{8}<\Delta$ for all $N$, including $N\to\infty$.
\end{assumption}
\begin{assumption}
\label{ass::4}Sequences $\+\mathbf{f}_t$, $\+\nu_{i,s}$ and $\varepsilon_{j,r}$ are independent for all $i,j,r,t$ and $s$.
\end{assumption}
\begin{assumption}
\label{ass::5}(a) The process $( \+\mathbf{z}_{i,t},t=1,\ldots,T)$ is initialized at an infinite past; (b) $|\alpha_0|<1$, and the spectral radius of $\*A_0$ is bounded by 1; (c) $\+\mathbf{z}_{i,0}$ is available for all $i=1,\ldots, N$.
\end{assumption}
\begin{assumption}
\label{ass::6}$NT^{-1}\to \kappa\in (0,\infty)$ as $N,T\to\infty$ jointly.
\end{assumption}
Assumptions \ref{ass::1}-\ref{ass::6} are fairly standard for the CCE literature; see e.g., \citet{Pesaran2006}, \citet{Westerlund2015}, and \citet{juodis2021robustness}. Assumption \ref{ass::3} is the rank condition of \citet{Pesaran2006}. Here by treating factor loadings as fixed non-stochastic quantities, we deviate from \citet{de2021bias} and \citet{Juodis2022CCER}, but expect that quantitatively similar results can be derived under the assumption of stochastic loadings (as long as the joint distribution of $\+\Gamma_i$ is left unrestricted). Assumption \ref{ass::6} is standard in the interactive fixed effects literature with weakly exogenous regressors; see e.g., \citet{moon2017dynamic}.
Finally, we assume that all stochastic quantities are covariance stationary and homoscedastic. This is mostly a technical assumption that substantially simplifies the derivations. Given that we consider wild bootstrap-based inference, all results are expected to extend beyond this restricted setting.
\begin{remark}
\textnormal{When analyzing the results for the AR(1) model, we will continue to specify the results assuming that all Assumptions \ref{ass::1}-\ref{ass::6} are satisfied. Then it is implicitly assumed that the corresponding parts associated with $\+\mathbf{x}_{i,t}$ are omitted.}
\end{remark}
\begin{remark}
\textnormal{The DGP in Eq. \eqref{eq::x} provides a convenient parameterization of $\+\mathbf{x}_{i,t}$ that accommodates three empirically relevant features: (i) a factor structure in the regressors; (ii) idiosyncratic serial correlation through $\+\nu_{i,t}$; and (iii) weak exogeneity with respect to $y_{i,t-1}$. We acknowledge that these features can be accommodated in other ways; see, for example, \citet{Juodis2022CCER} for an alternative DGP. However, our main conclusions do not depend critically on this particular specification.}
\end{remark}
\begin{remark}
\textnormal{In this paper, we assume that the rank condition is satisfied exactly, i.e. $\mathrm{rk}(\overline{\+\Gamma})=R=K$. \citet{Pesaran2006} allows for a less restricted version where $\mathrm{rk}(\overline{\+\Gamma})=R\leq K$. Allowing for such setup leads to additional non-trivial complications for asymptotic analysis under Assumption \ref{ass::6}; see e.g., \citet{Karabiyik2017} and \citet{de2024cross}. For this reason, we leave the analysis of this setup for future research.\footnote{Similarly to \citet{de2021bias}, the setting with $\mathrm{rk}(\overline{\+\Gamma})=R\leq K$ is expected to deliver quantitatively similar results as long as $\+\theta_{0}=\*0_{k}$ and $T/N\to 0$ and $N/T^{3}\to 0$. Assumption \ref{ass::6} needs to be modified appropriately to account for this possibility.} Alternatively, the number of factors can be selected using the tools suggested in the literature; see e.g., \citet{Juodis2022CCER}, \citet{10.1093/ectj/utad009}, and \citet{ditzen2026selection}.}
\end{remark}
\subsection{Results}
This section presents the paper’s two main results. We first discuss the properties of the recursive-design bootstrap in Algorithm \ref{algo::naive} for the special case of the panel AR(1) model; the results are stated in Theorem \ref{theorem::ar1}. The corresponding results for the more general ARX(1) model are presented in Theorem \ref{theorem::arx1}.
For what follows, let $\mathbb{P}^{*}(\cdot)$ be the bootstrap distribution function conditional on the realization of $\{(z_{i,t},i=1,\ldots,N; t=0,\ldots,T)\}$. In the AR(1) model the corresponding estimator is $\widehat{\alpha}_{CCEP}$ and $\alpha_{0}$ is the corresponding true value. For the general ARX(1) model, the corresponding quantities are given by $\widehat{\+\delta}_{CCEP}$ and $\+\delta_{0}$, respectively.
\begin{theorem}
\label{theorem::ar1}
If Assumptions \ref{ass::1}-\ref{ass::6} are satisfied, then in the AR(1) model
\begin{align*}
\sup_{x\in\mathbb{R}} \Big|\mathbb{P}^{*}\Big(\sqrt{NT}(\widehat{\alpha}_{CCEP}^*-\widehat{\alpha}_{CCEP}) \leq x \Big) - \mathbb{P}\Big(\sqrt{NT}(\widehat{\alpha}_{CCEP}-\alpha_0) &\leq x\Big)\Big| =o_P(1).
\end{align*}
\end{theorem}
The result in Theorem \ref{theorem::ar1} is positive: under Assumptions \ref{ass::1}-\ref{ass::6}, the RD-WB procedure fully replicates the first-order asymptotic distribution of the CCEP estimator. In particular, the bootstrap distribution has the same asymptotic variance and bias as the sampling distribution under the true DGP. This result extends Theorem 3.1 of \citet{gonccalves2015bootstrap} from the standard fixed effects setting to the CCEP setting.
The results for the CCEP estimator in the AR(1) model resemble those for the FE estimator because the CCEP estimator has a particularly simple asymptotic distribution in this case. In particular, following \citet{juodis2021robustness}, we have $\+\xi_{\perp}=\* b_{2}=0$ for decomposition in Eq. \eqref{eq::asymptotic_distribution}. Hence, factor-estimation error has no impact on the first-order asymptotic properties of the CCEP estimator. Moreover, the model is dynamically complete, as it contains no additional regressors whose dynamic properties are left unspecified.
Our next result summarizes the implications of departing from the ideal AR(1) model and considering the more general setting with regressors satisfying Eq. \eqref{eq::x}.
\begin{theorem}
\label{theorem::arx1}
If Assumptions \ref{ass::1}-\ref{ass::6} are satisfied, then in the ARX(1) model
\begin{align*}
\sup_{\*x\in\mathbb{R}^{K}}\Big|\mathbb{P}^{*}\Big(\+\Xi\sqrt{NT}(\widehat{\+\delta}_{CCEP}^*-\widehat{\+\delta}_{CCEP}-\Delta\*b^*)\leq \*x \Big) - \mathbb{P}\Big(\sqrt{NT}(\widehat{\+\delta}_{CCEP}-\+\delta_{0}) &\leq \*x\Big)\Big|=o_P(1),
\end{align*}
where
\begin{equation*}
\Delta\*b^{*}:= \frac{1}{T}\left(\*b_1^*-\*b_1\right) + \frac{1}{N}\left(\*b_2^*-\*b_2\right),
\end{equation*}
and $\+\Xi$ is some positive-definite matrix. All inequalities are interpreted coordinatewise.
\end{theorem}
The exact expressions of all terms provided in Theorem \ref{theorem::arx1} are provided in the corresponding proof in the Supplementary Online Appendix.
Overall, unlike in the AR(1) case studied in Theorem \ref{theorem::ar1}, the RD-WB procedure does not fully replicate the first-order asymptotic distribution of the CCEP estimator in the more general ARX(1) setting. In particular, the bootstrap and sampling distributions exhibit asymptotic discrepancies in both bias and variance components.
This negative conclusion is primarily driven by the fact that the regressors $\+\mathbf{x}_{i,t}$ are held fixed in the proposed bootstrap scheme. The next corollary summarizes the results for the case in which all regressors $\+\mathbf{x}_{i,t}$ are assumed to be strictly exogenous as in \citet{de2021bias}.
\begin{corollary}
If Assumptions \ref{ass::1}-\ref{ass::6} are satisfied, then in the ARX(1) model with $\+\theta_{0}=\*0_{k}$ then
\begin{align*}
\*b_1^*-\*b_1=\*0_{K}.
\end{align*}
\end{corollary}
Hence, when the regressors are assumed to be strictly exogenous, the RD-WB procedure fully replicates the corresponding ``Nickell bias'' of the CCEP estimator.
In Section \ref{section::takeaways}, we draw practical implications from Theorem \ref{theorem::arx1} for empirical researchers using CCEP estimators in dynamic panel data models. Among other things, we suggest that bias-correction methods proposed in the literature can be combined with the RD-WB procedure, similarly to the recommendation in \citet{gonccalves2015bootstrap}.
\begin{remark}
\textnormal{It is easy to see that if we were to extend the FE results in \citet{gonccalves2015bootstrap} to the setting with potentially weakly exogenous regressors, the results would be qualitatively similar to Theorem \ref{theorem::arx1} (except for $\*b_{2}^{*}-\*b_{2}=\*0_{K}$ and $\+\Xi=\mathbf{I}_{K}$ in that case).}
\end{remark}
In the next section, we intuitively explain the mechanisms behind the main conclusions of Theorem \ref{theorem::arx1}.
\subsection{Discussion}
\label{ssection::discussion}
The remaining negative aspects of Theorem \ref{theorem::arx1}
are not affected by exogeneity properties of regressors, and are solely driven by the fact that factor proxies $(\overline{\+\mathbf{x}}_{t},t=1,\ldots,T)$ remain fixed for all bootstrap replications. Hence, in the bootstrap world, these factors are no longer \emph{latent}, but, rather, observed.
We illustrate these features using the example presented in Eqs. \eqref{example:2eq1}-\eqref{example:2eq2}. For that model the bootstrap counterpart takes form
\begin{align*}
y_{i,t}^{*}&=\widehat{\alpha}_{CCEP} y_{i,t-1}^{*}+\widehat{\beta}_{CCEP} x_{i,t}+\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t}+\omega_{i,t}\widehat{\varepsilon}_{i,t},\\
x_{i,t}&=\+\gamma_{x,i}'\+\mathbf{f}_{t}+u_{i,t}.
\end{align*}
Note that, although the original model contains only $R=2$ latent factors, the bootstrap DGP contains four factors. Two of these-namely,
$\overline{\*x}$ and $\overline{\*x}_{-1}$ - are observed in the bootstrap world. Hence, only two factors remain latent.
As an intermediate step in the proof of Theorem \ref{theorem::arx1} we show that the loadings of the ($4-2=2$) ``excessive'' factors are asymptotically negligible, such that the bootstrap DGP is asymptotically equivalent to the DGP
\begin{align*}
y_{i,t}^{*}&=\widehat{\alpha}_{CCEP} y_{i,t-1}^{*}+\widehat{\beta}_{CCEP} x_{i,t}+\widetilde{\+\gamma}_{y,i}'\widetilde{\+\mathbf{f}}_{t}+\omega_{i,t}\widehat{\varepsilon}_{i,t},
\end{align*}
where $\widetilde{\+\mathbf{f}}_{t}$ is a $[\widetilde{R}^{*}\times 1]$ vector of rotated cross-sectional averages with $\widetilde{R}^{*}=2$. Here, the first element is given by $\widetilde{f}_{t}^{(1)}:=\overline{y}_{t}-\alpha_{0}\overline{y}_{t-1}-\beta_{0}\overline{x}_{t}$, while $\widetilde{f}_{t}^{(2)}:=\overline{x}_{t}$. Given that the factor proxies in the bootstrap world are given by $\widehat{\+\mathbf{f}}_{t}^{*}=[\overline{y}_{t}^{*},\overline{y}_{t-1}^{*},\overline{x}_{t},\overline{x}_{t-1}]$, $\widetilde{f}_{t}^{(1)}$ is the only latent factor that drives $y_{i,t}^{*}$.
This has implications for the bootstrap asymptotic distribution of the CCEP estimator. In particular, in the asymptotic distribution in Eq. \eqref{eq::asymptotic_distribution}, the two CCEP-specific components, $\+\xi_{\perp}$ and $\*b_{2}$, are determined by the factor loadings associated with the latent factors driving $y_{i,t}$. In the true DGP, there are generally $R=2$ such factors, whereas in the bootstrap DGP only $R^{*}=1$ factor remains latent. Consequently, the corresponding bootstrap terms $\+\xi_{\perp}^{*}$ and $\*b_{2}^{*}$ cannot be asymptotically equivalent to $\+\xi_{\perp}$ and $\*b_{2}$, respectively, because one factor is missing from the bootstrap DGP. The resulting bias discrepancy, $\*b_{2}^{*}-\*b_{2}$, is shown explicitly in the definition of $\Delta\*b^{*}$, whereas the variance discrepancy, $\+\xi_{\perp}^{*}-\+\xi_{\perp}$, determines the scaling factors $\+\Xi$.
The above discussion extends directly to the more general setting with an arbitrary number of regressors, $k$, and, subsequently, to settings with weakly exogenous regressors. These cases are fully considered in the Supplementary Online Appendix.
Finally, why in Algorithm \ref{algo::naive} we set $y_{i,0}^{*}=y_{i,0}$. This choice is not innocuous and significantly simplifies the asymptotic analysis. Using the example above, note that
\begin{equation*}
y_{i,t}^{*}=\widehat{\alpha}_{CCEP} y_{i,t-1}^{*}+\widehat{\beta}_{CCEP} x_{i,t}+\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t}+\omega_{i,t}\widehat{\varepsilon}_{i,t}.
\end{equation*}
Consider now what happens with the corresponding cross-sectional average $\overline{y}_{t}^{*}$:
\begin{equation}
\label{eq::boostrap_cs_average1}
\overline{y}_{t}^{*}=\widehat{\alpha}_{CCEP} \overline{y}_{t-1}^{*}+\widehat{\beta}_{CCEP} \overline{x}_{t}+\overline{\widehat{\+\gamma}}_{y}'\widehat{\+\mathbf{f}}_{t}+\overline{\omega\widehat{\varepsilon}}_{t}.
\end{equation}
Using the definition of $\widehat{\+\gamma}_{y,i}$ it is easy to see that $\overline{\widehat{\+\gamma}}_{y}'\widehat{\+\mathbf{f}}_{t}=\overline{y}_{t}-\widehat{\alpha}_{CCEP} \overline{y}_{t-1}-\widehat{\beta}_{CCEP} \overline{x}_{t}$.
Inserting this into (\ref{eq::boostrap_cs_average1}) together with $\overline{y}_{0}^{*}=\overline{y}_{0}$ gives
\begin{equation}
\label{eq::boostrap_cs_average2}
\overline{y}_{t}^{*}-\overline{y}_{t}=\widehat{\alpha}_{CCEP}( \overline{y}_{t-1}^{*}-\overline{y}_{t-1})+\overline{\omega\widehat{\varepsilon}}_{t}\quad \Rightarrow\quad \overline{y}_{t}^{*}=\overline{y}_{t}+\sum_{j=0}^{t-1}\widehat{\alpha}_{CCEP}^{j}\overline{\omega\widehat{\varepsilon}}_{t-j},
\end{equation}
for $t\geq 1$. It is evident that $\overline{y}_{t}^{*}$ is generally consistent for $\overline{y}_{t}$ (and likewise $\overline{y}_{t-1}^{*}$ is consistent for $\overline{y}_{t-1}$). This implies that the only latent factors in the bootstrap world - $\widetilde{f}_{t}^{(1)}$ - can be consistently estimated by a linear combination of $[\overline{y}_{t}^{*},\overline{y}_{t-1}^{*},\overline{x}_{t}]'$.
\section{Practical implications}
\setcounter{equation}{0}
\label{section::takeaways}
The negative result in Theorem \ref{theorem::arx1} raises the following question: \emph{Should the RD-WB procedure be used with the CCEP estimator at all?} We argue that the answer is yes, subject to appropriate modifications. We discuss these modifications below.
\subsection{Bias-correction}
\label{section:bias_correction}
First, we discuss what can be done in practice with the \emph{bias} wedge $\Delta \*b^{*}$ derived in Theorem \ref{theorem::arx1}. This term consists of two wedges - the ``Nickell bias'' wedge - $T^{-1}(\*b_{1}^{*}-\*b_{1})$, and the factor-approximation bias wedge - $N^{-1}(\*b_{2}^{*}-\*b_{2})$.
The factor-approximation bias of the CCEP estimator has generally received little attention of empirical researchers. As reviewed by \citet{JUODIS2026106120}, in many cases the bias itself is either assumed away (e.g. \citealp{de2021bias} assume $T/N\to 0$) or simply ignored. If one wishes to account for this bias, thus also account for the corresponding wedge between the distributions, it can be accounted for by using analytical bias-correction methods; see e.g., \citet{Westerlund2015} and \citet{Juodis2022CCER}. Suggested bias-correction methods generally do not fully remove this bias (hence, also the wedge), as part of the bias can be non-deterministic (see the corresponding discussion in \citealp{Juodis2022CCER}). On the other hand, the empirical consequence of the wedge $N^{-1}(\*b_{2}^{*}-\*b_{2})$ is expected to be limited, as the Monte Carlo results in the Supplementary Online Appendix indicate \footnote{All estimators considered in the Monte Carlo study do not explicitly account for the presence of the factor-estimation bias. This, however, has little impact both on reported biases as well as rejection rates.}
The wedge in the ``Nickell bias'' (when suspected to be present) generally should not be ignored in typical datasets where either $N>T$ or $N\approx T$ holds. As suggested in \citet{gonccalves2015bootstrap}, the RD-WB approach can be combined with any bias-corrected CCEP estimator that targets the ``Nickell bias'' of the CCEP estimator. The commonly used approaches are the Half Panel Jackknife (HPJ) bias-correction approach of \citet{dhaene2015split}, and the analytical bias-corrected estimator of \citet{HahnKursteiner2011BiasReduction} and \citet{moon2017dynamic}.\footnote{For the panel AR(1) model with additive fixed effects, \citet{gonccalves2015bootstrap} suggested using the analytical bias-corrected estimator of \citet{Hahn2002c}. Unfortunately, no similarly simple bias-corrected estimator is available for the factor-augmented setting considered here.}
Such bias-corrected CCEP estimators $\widehat{\+\delta}_{CCEP-bc}$ can be then used to re-estimate the factor loadings and the corresponding residuals
\begin{align}
\label{eq::cce_fitted1_bc}
\widehat{\+\gamma}_{y,i}&:= \left(\sum_{t=1}^{T}\widehat{\+\mathbf{f}}_{t}\widehat{\+\mathbf{f}}_{t}'\right)^{-1}\sum_{t=1}^{T}\widehat{\+\mathbf{f}}_{t}(y_{i,t}-\+\mathbf{w}_{i,t}'\widehat{\+\delta}_{CCEP-bc}),\\
\widehat{\varepsilon}_{i,t}&:= y_{i,t}-\+\mathbf{w}_{i,t}'\widehat{\+\delta}_{CCEP-bc}-\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t}.
\label{eq::cce_fitted2_bc}
\end{align}
The modified RD-WB algorithm is summarized below.
\begin{algorithm}
\label{algo::naive_bc}
\newcounter{bean2}
\setcounter{bean2}{0}
\begin{center}
\textnormal{
\begin{list}
{\textsc{Step} \arabic{bean}.}{\usecounter{bean}}
\item Obtain CCEP-bc loadings and residuals as in Eqs. \eqref{eq::cce_fitted1_bc}-\eqref{eq::cce_fitted2_bc}.
\item Set $y_{i,0}^{*}=y_{i,0}$.
\item For $b=1,\ldots,B$ generate bootstrap dataset $y_{i,t}^{*}$ recursively for $t=1,\ldots,T$,
\begin{equation}
\label{eq::bootstrap_DGP}
y_{i,t}^{*}:=\widehat{\alpha}_{CCEP-bc} y_{i,t-1}^{*}+\widehat{\+\beta}_{CCEP-bc}'\+\mathbf{x}_{i,t}+\widehat{\+\gamma}_{y,i}'\widehat{\+\mathbf{f}}_{t}+\omega_{i,t}^{*}\widehat{\varepsilon}_{i,t},
\end{equation}
for all $i=1,\ldots,N$, where $\omega_{i,t}^{*} \sim Rademacher(-1;1)$.
\item Given $((y_{i,t}^{*},\+\mathbf{x}_{i,t}')',i=1,\ldots,N,t=1,\ldots, T)$ obtain $\widehat{\+\delta}_{CCEP-bc}^{*}$ from Eq. \eqref{eq::LS_feasbile} using $\widehat{\+\mathbf{F}}^{*}=[\overline{\+\mathbf{y}}^{*},\overline{\+\mathbf{y}}^{*}_{-1},\overline{\+\mathbf{X}},\overline{\+\mathbf{X}}_{-1}]$
\item Repeat for $b=1,\ldots,B$.
\end{list}
}
\end{center}
\end{algorithm}
\subsection{Studentization}
\label{section:Studentization}
In general, the proposed bootstrap procedure cannot replicate the asymptotic variance of the CCEP estimator, since $\+\Xi\neq \mathbf{I}_{K}$. Below, we propose an empirical strategy that accounts for this feature. First, however, we show why the most natural, or naive, approach is not appropriate in this setting.
The most natural starting point is to use the usual sandwich variance estimator, with the clustered covariance matrix (CCM) of \citet{doi:10.1111/j.1468-0084.1987.mp49004006.x} as its middle component. This approach is intuitively appealing, as \citet{Cui03072023} show that the CCM consistently estimates the asymptotic variance of the IFE estimator of \citet{bai2009panel}.
Unfortunately, for the CCEP estimator, the CCM approach accounts only for variation arising from the $\+\xi_{0}$ term, but not for the variation associated with factor approximation, $\+\xi_{\bot}$; see the decomposition in Eq. \eqref{eq::asymptotic_distribution}. As discussed in Section \ref{ssection::discussion}, it is precisely the $\+\xi_{\bot}$ term that causes $\+\Xi\neq \mathbf{I}_{K}$.
This issue has been discussed by \citet{JuodisSarafidis2019}, \citet{Juodis2022CCER}, and \citet{Brown28052026}. In response, \citet{Juodis2022CCER} advocate the use of the pairs bootstrap,\footnote{In the context of their bias-corrected estimator, \citet{de2021bias} provide both a CCM-type estimator of the asymptotic variance and an estimator based on the pairs bootstrap. It can be shown that their proposed CCM-type estimator is generally inconsistent because it does not replicate the variance of $\boldsymbol{\+\xi}_\perp$ from (\ref{eq::asymptotic_distribution}). In contrast, the pairs bootstrap, also primarily advocated by \citet{de2021bias}, is consistent because it fully accounts for sampling uncertainty induced by factor estimation.} whereas \citet{JuodisSarafidis2019} and \citet{Brown28052026} propose modified CCM-type variance-matrix estimators. Unfortunately, these modifications do not translate directly to our setting, because the corresponding estimators are not consistent under fixed-$T$.
\footnote{\citet{Brown28052026} consider the CCEP setting with strictly exogenous regressors.}
As a solution, we suggest a jackknife-based variance estimator $\widehat{\+\delta}_{CCEP}$ (see \citealp{Tukey1958jackknife} and \citealp{Mackinnon2023leverage})
\begin{align}
\label{eq::jackknife_conventional}
\widehat{\text{Avar}}_{CCEP} = \frac{N-1}{N} \sum_{i=1}^N \left(\widehat{\+\delta}_{CCEP}^{(-i)} - \overline{\widehat{\+\delta}_{CCEP}}\right)\left(\widehat{\+\delta}_{CCEP}^{(-i)} - \overline{\widehat{\+\delta}_{CCEP}}\right)'.
\end{align}
Here $\widehat{\+\delta}_{CCEP}^{(-i)}$ is the CCEP estimate with $i-$th cross-sectional unit removed, while $\overline{\widehat{\+\delta}_{CCEP}}$ is the corresponding average of $N$ such estimates. For our purpose, it is critical, that within every $\widehat{\+\delta}_{CCEP}^{(-i)}$ also the factor estimates omit the $i-$th cross-sectional unit. The estimator in Eq. \eqref{eq::jackknife_conventional} is the ``conventional'' jackknife-based variance estimator using the terminology of \citet{Hansen2026Jackknife}. For the cross-sectional (clustered) setting \citet{Hansen2026Jackknife} suggested other variants of Eq. \eqref{eq::jackknife_conventional}. However, given relatively large $N$ available in a typical panel application, we do not expect any major differences between different jackknife estimators.
The jackknife approach in Eq. \eqref{eq::jackknife_conventional} while applicable for true data, cannot be directly used for bootstrap data $((y_{i,t}^{*},\+\mathbf{x}_{i,t}')',i=1,\ldots,N,t=1,\ldots, T)$. In particular, the delete-one analogue of Eq. \eqref{eq::boostrap_cs_average2} takes the form
\begin{equation}
\label{eq::boostrap_cs_average2_delete}
\overline{y}_{t}^{*(-i)}-\frac{N}{N-1}\overline{y}_{t}=\widehat{\alpha}_{CCEP}\left( \overline{y}_{t-1}^{*(-i)}-\frac{N}{N-1}\overline{y}_{t-1}\right)+\overline{\omega\widehat{\varepsilon}}_{t}^{(-i)}-\frac{1}{N-1}\left(\widehat{\beta}_{CCEP}x_{i,t}+\widehat{\+\gamma}_{y,i}'\widehat{\+
\mathbf{f}}_{t}\right).
\end{equation}
This extra term in Eq. \eqref{eq::boostrap_cs_average2_delete} prevents straightforward use of the jackknife methodology in this case. Intuitively, in the bootstrap realm there is no error from factor uncertainty that stems from $\+\mathbf{x}_{i,t}$, because we keep $\+\mathbf{x}_{i,t}$ fixed. The jackknife procedure reintroduces this error by omitting the $i-th$ cross-sectional unit, such that another wedge is created between the bootstrap realm variance and the jackknife estimate. Instead, as a solution to this problem, we use the idea recently highlighted in \citet{heller_jochmans_2026_iterated_bootstrap} and run a second round of $d=1,\ldots,D$ RD-WB bootstrap procedure within each $b=1,\ldots,B$ replications (though their application and motivation is different).
Iterated bootstrap of this form is applicable (even if do not attempt to formally prove validity), as long as regressors remain to be kept fixed in the second round of RD-WB iterations. The corresponding variance estimator for $\widehat{\+\delta}_{CCEP}^{*}$ is given by
\begin{align}
\label{eq::bootstrap_conventional}
\widehat{\text{Avar}}_{CCEP^{*}} = \frac{1}{D-1} \sum_{d=1}^D \left(\widehat{\+\delta}_{CCEP}^{**(d)} - \overline{\widehat{\+\delta}_{CCEP}^{**}}\right)\left(\widehat{\+\delta}_{CCEP}^{**(d)} - \overline{\widehat{\+\delta}_{CCEP}^{**}}\right)',
\end{align}
where $\widehat{\+\delta}_{CCEP}^{**(b)}$ is the CCEP estimate (given any bootstrap sample $b=1,\ldots,B$) in Algorithm \ref{algo::naive} for every second round iteration $d=1,\ldots,D$. In practice, we simply set $D=B$.\footnote{While theoretically Theorems \ref{theorem::ar1}-\ref{theorem::arx1} are not sufficient to prove consistency of the variance estimators of the form Eq. \eqref{eq::bootstrap_conventional}, and some form of trimming is needed to ensure consistency (see e.g., \citealp{Goncalves01092005}). We note that the simple proposal in Eq. \eqref{eq::bootstrap_conventional} works reasonably well in practice.}
\begin{remark}
\textnormal{Note that the studentization approach discussed above is only necessary in the context of models covered in Theorem \ref{theorem::arx1}. For the simple AR(1) model, while not necessary in practice, studentization can be nevertheless implemented. For that case (as well as other special cases we cover in Section \ref{section:extensions}) the simple CCM-based variance estimator suffices.}
\end{remark}
\subsection{Extensions}
\label{section:extensions}
Although the conclusions of Theorem \ref{theorem::arx1} are largely negative, the theorem also identifies conditions under which the positive conclusions of Theorem \ref{theorem::ar1} extend to more general settings. The two most straightforward extensions are panel autoregressive and panel vector autoregressive models of a finite order $p\geq 1$.
In the Supplementary Online Appendix, we sketch the argument establishing the validity of the RD-WB in the AR(p) setting. Unsurprisingly, the resulting conclusions fully mirror those of Theorem \ref{theorem::ar1}. For VAR(p) models, the failure to replicate the asymptotic distribution in the presence of additional regressors, as documented in Theorem \ref{theorem::arx1}, is driven solely by the fact that the regressors $\+\mathbf{x}_{i,t}$ and their corresponding cross-sectional averages, $\overline{\+\mathbf{x}}_{t}$, are held fixed throughout the bootstrap replications.
As a result, if we are willing to impose some structure on regressors and \emph{exploit} that structure explicitly in the construction of the bootstrap DGP, i.e. to make the full model for $\+\mathbf{z}_{i,t}$ dynamically complete, the original problem simplifies dramatically. For example, take the VAR(1) model in Eq. \eqref{z_VAR}
\begin{align*}
\+\mathbf{z}_{i,t}=\*A_0^{\prime}\+\mathbf{z}_{i,t-1}+\+\Gamma_i'\+\mathbf{f}_t+\+\mathbf{e}_{i,t}.
\end{align*}
The bootstrap counterpart takes the form
\begin{align}
\label{eq::var_z_star}
\+\mathbf{z}_{i,t}^{*}=\widehat{\+\mathbf{A}}_{CCEP}\+\mathbf{z}_{i,t-1}^{*}+\widehat{\+\Gamma}_i'\widehat{\+\mathbf{f}}_t+\omega_{i,t}\widehat{\mathbf{e}}_{i,t},
\end{align}
with $\+\mathbf{z}_{i,0}^{*}=\+\mathbf{z}_{i,0}$ as previously. Hence, unlike in the ARX(1) setting, all elements of $\*A_0$ must be estimated using CCEP-not only the equation for $y_{i,t}$- and the corresponding vector of CCEP-based residuals, $(\widehat{\mathbf{e}}_{i,t},i=1,\ldots,N;t=1,\ldots,T)$, must be constructed.\footnote{ DGP-consistent restrictions can be imposed when constructing $\widehat{\+\mathbf{A}}_{CCEP}$, such as known zero restrictions.} To replicate the asymptotic distribution of the CCEP estimator, it is crucial to use scalar weights $\omega_{i,t}$.\footnote{We provide the theoretical evidence on this case in Section 4 of Supplementary Online Appendix.}
Finally, although the RD-WB procedure summarized in Algorithm \ref{algo::naive} is expected to be consistent for dynamically complete models beyond the AR(1) model, it is not necessarily the most efficient CCEP-based implementation of the RD-WB. In Eq. \eqref{eq::var_z_star}, we use
\begin{equation*}
\widehat{\+\mathbf{f}}_t=[\overline{y}_{t},\overline{y}_{t-1}, \overline{\+\mathbf{x}}_{t}',\overline{\+\mathbf{x}}_{t-1}']',
\end{equation*}
so that $\widehat{\+\mathbf{f}}_t$ has $2K$ elements, whereas the original vector $\+\mathbf{f}_{t}$ has only $K=R$ elements under Assumptions \ref{ass::2}-\ref{ass::3}. Thus, after an appropriate rotation, the remaining $2K-K=K$ elements are asymptotically redundant. Although this rotation is generally unknown, as it depends on the full matrix $\*A_0$, it can be consistently estimated within the multi-equation VAR(1) model—and hence also in the single-equation AR(1) model.
Hence, the alternative version of the bootstrap algorithm uses
\begin{equation}
\label{eq::f_hat_restricted}
\ddot{\+\mathbf{f}}_t:=\overline{\+\mathbf{z}}_{t}-\widehat{\+\mathbf{A}}_{CCEP}'\overline{\+\mathbf{z}}_{t-1}.
\end{equation}
The factor loading matrix $\widehat{\+\Gamma}_i$ should be re-estimated accordingly using the factor estimates $\ddot{\+\mathbf{F}}$. For obvious reasons, the implementation of the bootstrap algorithm that uses the restricted factors Eq. \eqref{eq::f_hat_restricted} is referred to as the \emph{sophisticated}, while the original implementation as the \emph{naive} one.\footnote{Unlike \citet{Everaert2016}, we do not suggest a restricted (non-linear) version of the original CCEP estimator in the dynamically complete setting. The restrictions of the form Eq. \eqref{eq::f_hat_restricted} are only used in the construction of the bootstrap samples.} The theoretical properties of sophisticated schemes for AR(1) and AR(p) cases are explicitly addressed in Section 4 of Supplementary Online Appendix.
\begin{remark}
\textnormal{The decomposition in Eq. \eqref{eq::x} can be utilized without the need to fully specify the correlation structure in $\+\nu_{i,t}$ using some forms of (panel) Autoregressive Wild Bootstrap; see e.g., \citet{Juodis2020}. However, such approaches require stronger assumptions on the true DGP than utilized by the Algorithm \ref{algo::naive}.}
\end{remark}
\section{Empirical illustration}
\setcounter{equation}{0}
\label{section::empirical}
\subsection{Setup}
We apply the proposed methodology to re-investigate the dynamic effects of temperature shocks on Gross Domestic Product (GDP) growth. The premise was initially investigated by \citet{Dell2012Temperature} to assess the role of temperature in economic development and the impact of global warming on the future. Further work by \citet{Dell2014Weather} extended the model to capture unobserved heterogeneity by using a factor-augmented approach.
The climate panel data set contains $125$ countries observed over 1961-2003. Following \citet{Dell2012Temperature}, we account for the different effects of temperature of GDP growth for poor and rich countries. After applying this split and removing countries that have missing observations in either panel, we are left with a balanced panel of $93$ countries for the first panel, between 1962 and 1982.\footnote{See \citet{Dell2014Weather} who advocate for the splitting of the panel in two sub-periods, due to weather intensification or adaptation in recent years. The classification is based on an initial above or below median PPP-adjusted per capita GDP.} For the second panel over the period 1983-2003, the cross-section increases slightly to $118$ countries.
Given that $T$ is small in the resulting panels, \citet{de2021bias} advocate the use of their fixed-$T$ bias-corrected estimator, as temperature variable - the main regressor of interest - is expected to be strictly exogenous within the given time span. As indicated by the Monte Carlo results in the Supplementary Online Appendix, our proposed bootstrap procedure should be competitive for such values of $N,T$ even if it relies on the large $N,T$ asymptotic approximation.
Similarly to \citet{de2021bias}, we consider the following model
\begin{align}\label{eq:empirical}
g_{i,t} &= \alpha g_{i,t-1} + c_i + \beta_1 T_{i,t} + \beta_2 T_{i,t-1}+ \gamma_if_{t} + \varepsilon_{i,t}.
\end{align}
Here $g_{i,t}$ denotes the real per capita GDP growth, and $T_{i,t}$ denotes temperature. We further interact both temperature variables with a dummy variable indicating whether a country is rich or poor. This allows for heterogeneous exposure to temperature between developed and undeveloped economies. To proxy for the factors, we take the cross-sectional averages of all the available regressors.
The inclusion of $g_{i,t-1}$ in Equation \ref{eq:empirical} allows for persistent output growth, whereas the lagged temperature helps to distinguish the dynamic nature of the effect. The initial effect of a $1^{\circ}C$ increase in temperature on GDP growth is measured through $\beta_1$. A transitory shock has no permanent effect on output if $\beta_1 + \beta_2 = 0$. The implied cumulative growth effects (GE) in Equation \ref{eq:empirical} can be shown to be $\frac{\beta_1 + \beta_2}{1-\alpha}$. The vector $\widehat{\+\mathbf{f}}_{t}$ used in all our estimators is $[7\times 1]$. It includes averages of all right-hand-side (5 in total and the fixed effect) as well as the left-hand-side variable.
\subsection{Implementation}
Given the exogeneity of the temperature variable, the bias-correction from \citet{de2021bias} is the most obvious benchmark estimator in this setting (``DVS''). Moreover, we report the original CCE estimator (``CCEP''), as well as the Half-Panel-Jackknife procedure of \citet{dhaene2015split} (``HPJ''), and the analytical correction to the CCEP estimator from \citet{HahnKursteiner2011BiasReduction} (``AN''). Next to the label of the estimator considered, we also report the type of inference procedure used either ``(CS)'', that refers to cross-sectional bootstrap based inference methods or ``(RDn)'' that refers to the ``naive'' implementation of the recursive design bootstrap procedure.
For ``(CS)'' we follow \citet{de2021bias} and report the bootstrap-based standard errors (in the corresponding ``SE'' row), as well as bootstrap-based (equal-tailed) reverse percentile confidence intervals (the corresponding ``CI'' row). As advocated in Section \ref{section:Studentization}, for the recursive design based procedures ``(RDn)'', we report the leave-one-out jackknife based standard errors and double bootstrap-based confidence intervals (with studentization). Finally, following \citet{Higgins2024Bootstrap}, all point-estimates of the ``(RDn)'' estimators are bias-corrected using the median of the bootstrap distribution.
\subsection{Results}
In Tables \ref{tab:1962-1982}-\ref{tab:1983-2003}, we report the results for the first and the second sub-panels, respectively.
For the first sub-panel, we find that (focusing on the temperature variables) the main conclusions derived from the DVS and AN as well as CCEP (bias corrected by RDn) are comparable. If a poor country experiences a temperature shock it is found that there exists a statistically significant negative effect on GDP growth. However, this loss in growth is mostly compensated in the following year. As such, roughly 90\% of the growth loss is temporary, and the remainder being permanent. HPJ-based estimators generally result in smaller (in absolute value) temperature effect that is not found to be statistically significant (irrespective of the implementation used). This might serve as an indication of further time-series instabilities in the data beyond the original split suggested by \citet{Dell2014Weather}, or the fact that the length of the time-series is too short for precise estimation in every half-panel.\footnote{Note that for $T=20$ every half-panel has approximately $10$ observations, resulting in only $10-7$ effective degrees of freedom.}
For the second sub-panel, it is now found that the temperature effect is no longer statistically significant. In this implementation, both the contemporaneous and the lagged temperature variable are not statistically significant for a given statistical significance level. As \citet{Dell2012Temperature} argues, this could be seen as evidence that countries are adapting to more frequent swings in temperature. It is possible that countries which were classified as poor at the beginning of the sample, have developed into industries that are less affected by temperature shocks. This would reduce the exposure of the GDP growth rate of a country with respect to temperature. However, we find again, that the HPJ-based results tend to deviate the most from other estimators.
\begin{table}[h!]
\begin{center}
\caption{Estimation results for 1962-1982.}\label{tab:1962-1982}
\begin{adjustbox}{width=1\textwidth}
\begin{tabular}{l|cccc|ccc}
\hline
\hline
& CCEP (CS) & DVS (CS) & HPJ (CS) & AN (CS) & CCEP (RDn) & HPJ (RDn) & AN (RDn) \\
\hline
$g_{i,t-1}$ & 0.15** & 0.24** & 0.21** & 0.23** & 0.21** & 0.24** & 0.26** \\
SE & (0.08) & (0.08) & (0.09) & (0.08) & (0.09) & (0.10) & (0.09) \\
CI & (0.03, 0.32) & (0.11, 0.38) & (0.04, 0.38) & (0.10, 0.41) & (0.04, 0.52) & (0.06, 0.47) & (0.06, 0.50) \\
$T_{i,t}^{rich}$ & 0.47 & 0.48 & 0.18 & 0.48 & 0.44 & 0.10 & 0.42 \\
SE & (0.60) & (0.52) & (1.05) & (0.60) & (0.61) & (0.96) & (0.60) \\
CI & (-0.74, 1.49) & (-0.53, 1.41) & (-1.62, 2.45) & (-0.63, 1.56) & (-1.08, 1.78) & (-2.16, 2.29) & (-1.14, 1.76) \\
$T_{i,t-1}^{rich}$ & -0.35 & -0.39 & -0.82 & -0.38 & -0.49 & -0.99 & -0.48 \\
SE & (0.56) & (0.53) & (0.85) & (0.53) & (0.64) & (0.97) & (0.62) \\
CI & (-1.77, 0.45) & (-1.70, 0.40) & (-2.63, 0.50) & (-1.79, 0.45) & (-1.96, 0.86) & (-3.09, 1.05) & (-1.94, 0.94) \\
$T_{i,t}^{poor}$ & -1.94** & -1.93** & -0.48 & -1.93** & -2.05** & -0.34 & -2.10 \\
SE & (0.90) & (0.91) & (1.51) & (0.84) & (0.92) & (1.53) & (0.93) \\
CI & (-4.09, -0.40) & (-4.00, -0.33) & (-2.20, 3.31) & (-3.75, -0.31) & (-4.38, -0.09) & (-4.29, 2.80) & (-4.77, 0.24) \\
$T_{i,t-1}^{poor}$ & 1.76** & 1.84** & 1.77 & 1.83** & 1.95** & 1.86 & 1.92 \\
SE & (0.95) & (0.92) & (1.39) & (0.89) & (0.99) & (1.70) & (1.01) \\
CI & (0.13, 3.84) & (0.07, 3.73) & (-0.31, 5.26) & (0.22, 3.62) & (0.06, 4.35) & (-1.58, 6.94) & (-0.40, 4.39) \\
\hline
GE Rich Countries & 0.14 & 0.12 & -0.80 & 0.12 & -0.07 & -1.17 & -0.08 \\
SE & (1.11) & (1.08) & (1.92) & (1.14) & (1.36) & (2.21) & (1.39) \\
CI & (-2.31, 1.88) & (-2.24, 1.95) & (-3.99, 3.51) & (-2.04, 2.15) & (-3.71, 2.07) & (-4.74, 4.69) & (-3.49, 2.65) \\
GE Poor Countries & -0.21 & -0.12 & 1.62 & -0.13 & -0.12 & 2.00 & -0.24 \\
SE & (1.18) & (1.35) & (2.40) & (1.36) & (1.40) & (2.96) & (1.44) \\
CI & (-2.33, 2.11) & (-3.16, 2.63) & (-0.80, 8.48) & (-2.60, 2.60) & (-2.59, 2.98) & (-6.46, 6.19) & (-3.22, 2.81) \\
\hline
\hline
\end{tabular}
\end{adjustbox}
\end{center}
\footnotesize
\textbf{Note:} Bootstrapped standard deviations (SE) are shown in brackets. The 95\% confidence interval (CI) are shown as the tuple in brackets. ** denotes significance at level 5\%, using the confidence interval. \\ From left to right, the shown estimators are (CCEP) \citet{Pesaran2006}; (DVS) \citet{de2021bias}; (HPJ) \citet{dhaene2015split}; (AN) \citet{HahnKursteiner2011BiasReduction}. Estimators with the suffix (RDn) are supplied with the proposed recursive design bootstrap. In the remaining cases, (CS) is used to denote the cross-sectional bootstrap.
\end{table}
\begin{table}[h!]
\begin{center}
\caption{Estimation results for 1983-2003.}\label{tab:1983-2003}
\begin{adjustbox}{width=\textwidth}
\begin{tabular}{l|cccc|ccc}
\hline
\hline
& CCEP (CS) & DVS (CS) & HPJ (CS) & AN (CS) & CCEP (RDn) & HPJ (RDn) & AN (RDn) \\
\hline
$g_{i,t-1}$ & 0.07 & 0.22** & 0.27** & 0.16** & 0.16** & 0.23 & 0.21** \\
SE & (0.06) & (0.07) & (0.13) & (0.06) & (0.09) & (0.22) & (0.09) \\
CI & (-0.08, 0.17) & (0.07, 0.35) & (0.03, 0.53) & (0.05, 0.27) & (0.06, 0.48) & (-0.44, 0.72) & (0.08, 0.45) \\
$T_{i,t}^{rich}$ & 0.47 & 0.44 & 0.81 & 0.45 & 0.47 & 0.82 & 0.53 \\
SE & (0.39) & (0.38) & (0.65) & (0.38) & (0.45) & (0.83) & (0.43) \\
CI & (-0.29, 1.27) & (-0.26, 1.23) & (-0.34, 2.12) & (-0.34, 1.15) & (-0.60, 1.36) & (-1.23, 2.77) & (-0.52, 1.55) \\
$T_{i,t-1}^{rich}$ & 0.09 & 0.08 & 0.41 & 0.08 & 0.14 & 0.44 & 0.08 \\
SE & (0.33) & (0.35) & (0.60) & (0.33) & (0.39) & (0.80) & (0.37) \\
CI & (-0.56, 0.75) & (-0.64, 0.73) & (-0.70, 1.72) & (-0.55, 0.75) & (-0.85, 0.92) & (-1.09, 2.30) & (-0.64, 0.81) \\
$T_{i,t}^{poor}$ & -1.11 & -1.24 & -0.56 & -1.19** & -1.09 & -0.49 & -1.07 \\
SE & (0.67) & (0.66) & (1.00) & (0.62) & (0.75) & (1.22) & (0.77) \\
CI & (-2.31, 0.47) & (-2.47, 0.13) & (-2.39, 1.60) & (-2.39, -0.03) & (-2.56, 0.80) & (-2.89, 2.30) & (-2.70, 0.71) \\
$T_{i,t-1}^{poor}$ & 0.30 & 0.57 & 0.46 & 0.47 & 0.42 & 0.48 & 0.50 \\
SE & (0.72) & (0.70) & (1.23) & (0.72) & (0.81) & (1.74) & (0.83) \\
CI & (-1.26, 1.43) & (-0.78, 1.86) & (-2.08, 2.49) & (-0.91, 1.92) & (-1.40, 2.41) & (-3.44, 4.86) & (-1.37, 2.64) \\
\hline
GE Rich Countries & 0.60 & 0.66 & 1.67 & 0.64 & 0.72 & 1.63 & 0.78 \\
SE & (0.59) & (0.69) & (1.43) & (0.63) & (0.79) & (1.71) & (0.80) \\
CI & (-0.58, 1.70) & (-0.67, 1.98) & (-1.07, 4.67) & (-0.56, 1.88) & (-2.00, 1.11) & (-6.53, 1.76) & (-1.94, 1.24) \\
GE Poor Countries & -0.87 & -0.87 & -0.13 & -0.87 & -0.79 & -0.02 & -0.73 \\
SE & (0.85) & (0.87) & (2.08) & (0.86) & (0.88) & (2.60) & (0.90) \\
CI & (-2.49, 0.81) & (-2.49, 0.76) & (-4.38, 3.61) & (-2.49, 0.89) & (-1.13, 2.31) & (-6.05, 8.80) & (-1.16, 2.20) \\
\hline
\hline
\end{tabular}
\end{adjustbox}
\end{center}
\footnotesize
\textbf{Note:} Bootstrapped standard deviations (SE) are shown in brackets. The 95\% confidence interval (CI) are shown as the tuple in brackets. ** denotes significance at level 5\%, using the confidence interval. \\ From left to right, the shown estimators are (CCEP) \citet{Pesaran2006}, (DVS) \citet{de2021bias}, (HPJ) \citet{dhaene2015split}, and (AN) \citet{HahnKursteiner2011BiasReduction}. Estimators with the suffix (RDn) are supplied with the proposed recursive design bootstrap. In the remaining cases, (CS) is used to denote the cross-sectional bootstrap.
\end{table}
\section{Concluding remarks}
\setcounter{equation}{0}
\label{section::conclusions}
This paper studies the validity of the recursive-design wild bootstrap (RD-WB) for inference in linear dynamic panel-data models with common factors (interactive fixed effects). Our analysis is based on the CCEP estimator of \citet{Pesaran2006}, which approximates the factor structure using cross-sectional averages of observed variables. We establish the asymptotic validity of the bootstrap in the panel AR(1) setting. The Monte Carlo results show that the proposed procedure yields coverage rates that are at least comparable to, and often improve upon, those of commonly used alternatives across a range of designs.
We then examine this bootstrap algorithm in settings with additional weakly exogenous regressors. In the bootstrap world, we adopt an agnostic approach by holding the regressors fixed. This induces discrepancies in both the bias and variance between the sampling distribution and its bootstrap counterpart. For empirical applications, we propose a statistical toolkit that combines the RD-WB procedure with commonly used bias-correction methods and appropriate studentization. Although the asymptotic results for the ARX(1) setting are more nuanced, Monte Carlo evidence shows that RD-WB achieves coverage rates close to their nominal levels over a broad range of designs.
We restrict attention to the CCEP estimator of \citet{Pesaran2006} motivated by the fixed-$T$ consistency of factor estimates based on cross-sectional averages. This choice is primarily responsible for the negative aspects identified in Theorem \ref{theorem::arx1}. Alternatively, we could use the principal-components (PC) approach of \citet{GreenawayMcGrevy201248}; see also \citet{Westerlund2015} and \citet{juodis2026factoraugmentedpanelregressionsvarianceweighted}. We conjecture that the RD-WB is valid for this class of estimators, at least when the regressors are strictly exogenous. A formal analysis of this extension is beyond the scope of this study.
\section*{Acknowledgments}
\thanks{We thank Otilia Boldea, S\'{i}lvia Gon\c{c}alves, Ayden Higgens, and the participants of the NESG 2026 (Tilburg), Conference in Honour of James Mackinnon (Aarhus), IPDC 2026 (Exeter). Financial support from the Dutch Research Council (NWO) under research grant VI.Vidi.231E.030 is gratefully acknowledged by all authors. This paper benefited from the use of generative AI tools to assist with language editing and \LaTeX formatting. All output was carefully reviewed by the authors. All substantive content, results, and any remaining errors are the authors’ responsibility.}
\bibliographystyle{chicago}
\bibliography{biblio_ectj_new}
\clearpage