EconBase
← Back to paper

An alternative bootstrap procedure for factor-augmented regression models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

68,700 characters · 3 sections · 61 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

An alternative bootstrap procedure for factor-augmented regression models

abstractIn this paper, we propose a novel bootstrap algorithm that is more efficient than existing methods for approximating the distribution of the factor-augmented regression estimator for a rotated parameter vector. The regression is augmented by $r$ factors extracted from a large panel of $N$ variables observed over $T$ time periods. We consider general weak factor (WF) models with $r$ signal eigenvalues that may diverge at different rates, $N^{\alpha _{k}}$, where $0<\alpha _{k}\leq 1$ for $k=1,2,...,r$. We establish the asymptotic validity of our bootstrap method using not only the conventional data-dependent rotation matrix $\hat{\mathbf{H}}$, but also an alternative data-dependent rotation matrix, $\hat{\mathbf{H}}_q$, which typically exhibits smaller asymptotic bias and achieves a faster convergence rate. Furthermore, we demonstrate the asymptotic validity of the bootstrap under a purely signal-dependent rotation matrix ${\mathbf{H}}$, which is unique and can be regarded as the population analogue of both $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. Experimental results provide compelling evidence that the proposed bootstrap procedure achieves superior performance relative to the existing procedure.

Keywords. Factor model, Asymptotic bias, Bootstrap, Weak factors {\let\thefootnote\relax\footnote{$^*$Corresponding author. Email: [email removed]; Address: Faculty of Economics and Business Administration, Tokyo Metropolitan University, 1-1 Minami-Osawa, Hachioji-shi, Tokyo, Japan 192-0397}}

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Introduction}

Factor-augmented regressions are widely used in financial and economic research. They are often used to forecast macroeconomic and financial time series. The forecast regression is augmented with a few common factors extracted from a large set of predictors. Specifically, the $h$-ahead forecast regression of $y$ is written as

align[align omitted — 144 chars of source]

where $\mathbf{f}_t^*$ is an $r \times 1$ vector of latent predictive factors and $\mathbf{w}_t$ is a $p \times 1$ vector of observable predictors. Since $\mathbf{f}_t^*$ is unobserved, it is typically replaced by the principal component (PC) estimator, $\hat{\mathbf{f}}_t$, which satisfies $T^{-1}\sum_{t=1}^T\hat{\mathbf{f}}_t\hat{\mathbf{f}}_t^{\prime}=\mathbf{I}_r$, and is constructed from the $\sqrt{T}$ times $r$ eigenvectors corresponding to the $r$ largest eigenvalues ($\hat{\lambda}_{1}>\dots>\hat{\lambda}_{r}$) of the $T \times T$ sample covariance matrix of $N$ predictors, $\{x_{t,i}\}_{i=1}^N$, which are assumed to follow a latent factor structure:

align[align omitted — 124 chars of source]

Note that most of the existing literature assumes that the $r$ largest eigenvalues of the sample covariance matrix of $x_{t,i}$, $(\hat{\lambda}_{1},\dots,\hat{\lambda}_{r})$, diverge proportionally with $N$. This is known as the strong factor (SF) model. In contrast, we present results for more general, so-called weak factor (WF) models, in which each $\hat{\lambda}_{k}$ can diverge at a different rate $N^{\alpha_k}$, with $\alpha_{1}\geq\dots \geq{\alpha}_{r}$, $\alpha_k \in (0,1]$, $k=1,2,...,r$. A growing body of literature suggests that such weak factors are prevalent in real-world data. See, for example, BaileyEtAl2016,BaileyEtAl2021, DeMol2008, Freyaldenhoven21JoE, Onatski2010, UY2019,UY2019inference, WeiZhang2023, among many others.

Let $(\hat{\boldsymbol{\gamma}}',\hat{\boldsymbol{\beta}}')'$ be the least squares estimators of the regression of $y_{t+h}$ on $(\hat{\mathbf{f}}_t^{\prime}, \mathbf{w}_t^{\prime})^{\prime}$. For SF models, StockWatson2002JASA, BaiNg2006 and GoncalvesPerron2014,gonccalves2020bootstrapping employ an asymptotic approximation in which the PC factor approximates a rotated version of the latent factor, using a data-dependent (but infeasible) rotation matrix:

align[align omitted — 89 chars of source]

where $\hat{\mathbf{H}}= \sum_{i=1}^N{\mathbf{b}}_{i}^{*}{\mathbf{b}}_i^{*\prime}T^{-1}\sum_{t=1}^T{\mathbf{f}}_t^{*}\hat{\mathbf{f}}_t^{\prime} \hat{\boldsymbol{\Lambda}}^{-1}$ with $\hat{\boldsymbol{\Lambda}}=\diag{(\hat{\lambda}_{1}\cdots\hat{\lambda}_{r})}$. Note that $\hat{\mathbf{H}}$ is data dependent but not estimable as it depends on the unobserved $(\mathbf{f}_t^* ,\mathbf{b}_i^* )$. Using the rotation matrix $\hat{\mathbf{H}}$, the first term on the right-hand side of the forecast regression (ref) can be written as $\boldsymbol{\gamma}^{*\prime}\mathbf{f}_t^* = \boldsymbol{\gamma}^{*\prime}\hat{\mathbf{H}}^{-1\prime}\hat{\mathbf{H}}^{\prime}{\mathbf{f}}_t^* = \boldsymbol{\gamma}_{\hat{\mathbf{H}}}^{\prime}\hat{\mathbf{f}}_t + o_p(1)$, where $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}=\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*}$ is effectively what $\hat{\boldsymbol{\gamma}}$ estimates. BaiNg2006 show that as long as $\sqrt{T}/N \to 0$, the limiting distribution of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ is centered at zero (i.e., there is no asymptotic bias). Under a relaxed condition of $\sqrt{T}/N \to c\in(0,\infty)$, Ludvigson2011 show that $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ exhibits an asymptotic bias and provide an analytical bias correction for SF models. GoncalvesPerron2014 refine the asymptotic bias expression and propose an analytical bias correction. gonccalves2020bootstrapping extend the results of GoncalvesPerron2014 to allow for bias corrections when errors $e_{t,i}$ are cross-correlated, using the method for estimating large covariance matrices proposed by BickelLevina2008.

GoncalvesPerron2014,gonccalves2020bootstrapping propose a bootstrap procedure to correct the asymptotic bias. Noting that, in the bootstrap world, we can “observe” the population -- including $\hat{\mathbf{H}}$ -- and recalling that $\hat{\boldsymbol{\gamma}}$ can be viewed as an estimator of $\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*}$, it becomes possible to construct an estimator of $\boldsymbol{\gamma}^*$, namely $\hat{\mathbf{H}}\hat{\boldsymbol{\gamma}}$, in the bootstrap world . GoncalvesPerron2014,gonccalves2020bootstrapping essentially propose to obtain the empirical distribution of $\sqrt{T}(\hat{\mathbf{H}}\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^*)=\sqrt{T}\hat{\mathbf{H}}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ via bootstrap, to approximate the limiting distribution of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ given that $\hat{\mathbf{H}} \stackrel{p}{\longrightarrow} \mathbf{I}_r$ in the bootstrap world.

In this paper, we propose a simple and alternative bootstrap procedure, in which the bootstrap distribution of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ is directly constructed as is. We establish the asymptotic validity of this bootstrap procedure, and finite-sample experiments suggest that our method generally provides a more accurate distributional approximation.

Equipped with this new bootstrap procedure, we further consider bootstrapping the distribution of $\hat{\boldsymbol{\gamma}}$ relative to two alternative rotation matrices. As introduced by BaiNg2023 and jiang2023revisiting, there exist variants of asymptotically equivalent, data-dependent rotation matrices other than $\hat{\mathbf{H}}$. Among these, we consider $\hat{\mathbf{H}}_q=(T^{-1}\sum_{t=1}^T \hat{\mathbf{f}}_t{\mathbf{f}}_t^{*\prime})^{-1}$, and propose bootstrapping $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q})$, where $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}=\hat{\mathbf{H}}_{q}^{-1}\boldsymbol{\gamma}^{*}$. In addition, jiang2023revisiting show that a unique (up to sign) rotation matrix $\mathbf{H}$ always exists, which is a function of the signals $(\mathbf{f}_t^*,\mathbf{b}_i^*)$ for $t=1,...,T$ and $i=1,...,N$ only, such that

align[align omitted — 65 chars of source]

where $T^{-1}\sum_{t=1}^T{\mathbf{f}}_t^0 {\mathbf{f}}_t^{0\prime}=\mathbf{I}_r$. Indeed, $\mathbf{H}$ can be seen as the population version of $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. By substituting (ref) into (ref), the forecasting model can be equivalently expressed as

align[align omitted — 123 chars of source]

where $\boldsymbol{\gamma}^0= {\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*}$. Since $\hat{\mathbf{f}}_t$ is consistent to $\mathbf{f}_t^0$ (up to sign) as shown by jiang2023revisiting, regressing $y_{t+h}$ on $(\hat{\mathbf{f}}_t,\mathbf{w}_t)$ consistently estimates the parameter vector $(\boldsymbol{\gamma}^{0\prime},\boldsymbol{\beta}')'$. In this paper, we propose bootstrapping $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^0)$, which we recommend especially for inference on linear restrictions on $\boldsymbol{\gamma}^0$.

Some readers may wonder whether, to obtain a data-independent rotation matrix, it would be sufficient to consider the probability limit ${\mathbf{H}}_0$ of $\hat{\mathbf{H}}$ as $(N,T)\rightarrow\infty$. However, even if ${\mathbf{H}}_0$ is well-defined, an approach based on it requires an additional layer of approximation. To approximate $(\hat{\mathbf{f}}_t,\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^{})$ by $({\mathbf{H}}_0'{\mathbf{f}}_t^{},{\mathbf{H}}_0^{-1}\boldsymbol{\gamma}^{*})$, one must first invoke (ref), and then proceed with the approximation using the probability limit ${\mathbf{H}}_0$ as $(N,T)\rightarrow\infty$. In contrast, our approach directly approximates $(\hat{\mathbf{f}}_t,\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^{})$ by $({\mathbf{H}}'{\mathbf{f}}_t^{},{\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*})$, where $\mathbf{H}$ is given at finite $\{N,T\}$.

The finite sample performance of the proposed bootstrap bias correction is compared with the methods of GoncalvesPerron2014,gonccalves2020bootstrapping under both strong and weak factor models. The results confirm that our bootstrap procedure generally provides a more accurate approximation, leading to further bias reduction.

The rest of the paper is organized as follows. Section (ref) introduces models and estimators relative to the latent parameter vector rotated by $\hat{\mathbf{H}}$. Section (ref) proposes a new bootstrap procedure and introduces two alternative rotation matrices. Section (ref) states assumptions and presents theoretical results. Section (ref) discusses finite-sample experiments, and Section (ref) concludes. Mathematical proofs are provided in the Online Appendix.

Notations: Denote by $\lambda_k[\mathbf{A}]$ the $k$th largest eigenvalue of a square matrix $\mathbf{A}$. For any matrix $\mathbf{M}=(m_{t,i})\in\mathbb{R}^{T\times N}$, we define the Frobenius norm and $\ell_2$-induced (spectral) norm as $\|\mathbf{M}\|_{\F}=(\sum_{t,i}m_{t,i}^2)^{1/2}$ and $\|\mathbf{M}\|_2=\lambda_{1}^{1/2}(\mathbf{M}'\mathbf{M})$, respectively. We denote the identity matrix of order $s$ by $\mathbf{I}_s$ and $s\times 1$ vectors of ones and zeros by $\mathbf{1}_s$ and $\mathbf{0}_s$, respectively. $\lesssim$ ($\gtrsim$) represents $\leq$ ($\geq$) up to a positive constant factor. $\odot$ denotes the Hadamard product of matrices. For any positive sequences $a_n$ and $b_n$, we write $a_n \asymp b_n$ if $a_n \lesssim b_n$ and $a_n \gtrsim b_n$. All asymptotic results are for cases where $N,T\to\infty$, and we omit explicit mention of this unless necessary. $M$ denotes a positive constant which does not depend on $N$ and $T$.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Factor-Augmented Regression}

The factor-augmented regression model (ref) can be rewritten in matrix form as:

align[align omitted — 198 chars of source]

where $\mathbf{y}=(y_{1+h},\dots,y_{T+h})'$, $\boldsymbol{\epsilon}=(\epsilon_{1+h},\dots,\epsilon_{T+h})'$, $\mathbf{F}^{*} = (\mathbf{f}_1^{*},\dots,\mathbf{f}_T^{*})'$, $\mathbf{W}=(\mathbf{w}_1,\cdots, \mathbf{w}_T)'$, $\mathbf{Z}^* = (\mathbf{F}^*,\mathbf{W})$ and $\boldsymbol{\delta}^* = (\boldsymbol{\gamma}^{*\prime},\boldsymbol{\beta}')'$. In line with (ref), the latent factor model for the $T \times N$ matrix of predictors is given by

align[align omitted — 97 chars of source]

where $\mathbf{X} = (x_{t,i})$, $\mathbf{B}^* = (\mathbf{b}_1^* ,\dots,\mathbf{b}_N^*)'$ and $\mathbf{E} = (e_{t,i})$. Let $(\lambda_1>\cdots>\lambda_r)$ denote the $r$ largest eigenvalues of the signal component of the model (ref), namely $T^{-1}\mathbf{F}^*\mathbf{B}^{\ast\prime}\mathbf{B}^{\ast}\mathbf{F}^{\ast\prime}$, and define $\boldsymbol{\Lambda}=\diag(\lambda_1,\dots,\lambda_r)$. We allow the $r$ signal eigenvalues to diverge at different rates, specifically $\lambda_k \asymp N^{\alpha_k}$ with $0<\alpha_k\leq 1$ for $k=1,\dots,r$. We refer to model (ref) with $\alpha_r=1$ as a strong factor (SF) model, and the more general model without this restriction as a weak factor (WF) model.

jiang2023revisiting show that there always exists a unique (up to sign) rotation matrix

align[align omitted — 66 chars of source]

where $\mathbf{P}$ is the eigenvector matrix of ${\mathbf{B}^*}'\mathbf{B}^*(T^{-1}{\mathbf{F}^*}'\mathbf{F}^*)$ corresponding to $(\lambda_1,\dots,\lambda_r)$ and $\mathbf{V}=\mathbf{P}(T^{-1}{\mathbf{F}^*}'\mathbf{F}^*)\mathbf{P}'$, such that

align[align omitted — 113 chars of source]

which by construction satisfy the $r^2$ restrictions $T^{-1}\mathbf{F}^{0\prime}\mathbf{F}^0 = \mathbf{I}_r$ and $\mathbf{B}^{0\prime}\mathbf{B}^0 = \boldsymbol{\Lambda}$.

Therefore, $\mathbf{H}$ is a pure function of signals $(\mathbf{F}^* , \mathbf{B}^*)$. It can also be straightforwardly shown that

align[align omitted — 185 chars of source]

With this rotation, we can equivalently express models (ref) and (ref) in terms of $\mathbf{F}^0$, $\boldsymbol{\gamma}^0 = \mathbf{H}^{-1}\boldsymbol{\gamma}^*$, and $\mathbf{B}^0$, which define the pseudo-true models:

align[align omitted — 292 chars of source]

where $\mathbf{Z}^0 = (\mathbf{F}^0,\mathbf{W})$ and $\boldsymbol{\delta}^0=\boldsymbol{\Phi}_{\mathbf{H}}^{-1}\boldsymbol{\delta}^* with \boldsymbol{\Phi}_{{\mathbf{H}}}= \big(

smallmatrix{\mathbf{H}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_p

\big)$. \cite{StockWatson2002JASA} propose extracting principal component (PC) factors from the predictor matrix $\mathbf{X}$ and using them in the forecast regression. The PC estimator, $(\hat{\mathbf{F}}, \hat{\mathbf{B}})$, is defined as the solution to the minimization problem $\left\|\mathbf{X}-\mathbf{F B}^{\prime}\right\|_{\mathrm{F}}^2$ subject to the $r^2$ constraints: $T^{-1} \mathbf{F}^{\prime} \mathbf{F}=\mathbf{I}_r$ and $\mathbf{B}^{\prime} \mathbf{B}$ being a diagonal matrix with rank $r$. The constrained minimization reduces to the eigenvalue problem of $T^{-1} \mathbf{X X}^{\prime}$. The factor estimator $\hat{\mathbf{F}} \in \mathbb{R}^{T \times r}$ is obtained as $\sqrt{T}$ times the $r$ eigenvectors associated with the $r$ largest eigenvalues of $T^{-1} \mathbf{X} \mathbf{X}^{\prime}$ $(\hat{\lambda}_1>\cdots>\hat{\lambda}_r)$, and the loading estimator $\hat{\mathbf{B}} \in \mathbb{R}^{N \times r}$ is computed as $\hat{\mathbf{B}}=T^{-1} \mathbf{X}^{\prime} \hat{\mathbf{F}}$. By construction, $T^{-1} \hat{\mathbf{F}}^{\prime} \hat{\mathbf{F}}=\mathbf{I}_r$ and $\hat{\mathbf{B}}^{\prime} \hat{\mathbf{B}}=\hat{\boldsymbol{\Lambda}}=\operatorname{diag}(\hat{\lambda}_1, \ldots, \hat{\lambda}_r)$.

Then, regressing $\mathbf{y}$ on $\hat{\mathbf{Z}}=(\hat{\mathbf{F}},\mathbf{W})$ yields the least squares estimator

align[align omitted — 113 chars of source]

Hence, the PC estimator $\hat{\mathbf{F}}$ can naturally be viewed as an estimator of $\mathbf{F}^0$, and $\hat{\boldsymbol{\delta}}$ as an estimator of the parameter vector $\boldsymbol{\delta}^0$ in the pseudo-true models (ref) and (ref).

BaiNg2002,BaiNg2006, StockWatson2002JASA, consider the approximation

align[align omitted — 83 chars of source]

where $\hat{\mathbf{H}} = {\mathbf{B}^*}'\mathbf{B}^*(T^{-1}{\mathbf{F}^*}'\hat{\mathbf{F}}) \hat{\boldsymbol{\Lambda}}^{-1}$. Comparing this to (ref), we see that $\mathbf{H}$ is the population analogue of $\hat{\mathbf{H}}$. With respect to $\mathbf{F}^0$ in the pseudo true models (ref) and (ref), we can establish the following key identity:

align[align omitted — 95 chars of source]

where $\tilde{\mathbf{H}} :=\mathbf{H}^{-1}\hat{\mathbf{H}} = {\mathbf{B}^0}'\mathbf{B}^0(T^{-1}{\mathbf{F}^0}'\hat{\mathbf{F}}) \hat{\boldsymbol{\Lambda}}^{-1}$. In the same way that $\hat{\mathbf{H}}$ is considered an estimator of $\mathbf{H}$, $\tilde{\mathbf{H}}$ can be viewed as an estimator of the identity matrix $\mathbf{I}_r$. Using (ref), the first term on the right-hand side of the augmented model (ref) can be written as

align[align omitted — 216 chars of source]

where $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}=\hat{\mathbf{H}}^{-1} \boldsymbol{\gamma}^*=\tilde{\mathbf{H}}^{-1} \boldsymbol{\gamma}^0$.

Based on the approximation in (ref) and the identities (ref) and (ref), $\hat{\boldsymbol{\delta}}$ can be regarded as an estimator of $\boldsymbol{\delta}_{\hat{\mathbf{H}}}:=(\boldsymbol{\gamma}_{\hat{\mathbf{H}}}',\boldsymbol{\beta}')'$. jiang2024Mw derive the asymptotic distribution of $ \sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}}) $ together with its asymptotic bias, where

align[align omitted — 200 chars of source]

with $\boldsymbol{\Phi}_{\hat{\mathbf{H}}}= \big(

smallmatrix\hat{\mathbf{H}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_p

\big)$ and $\boldsymbol{\Phi}_{\tilde{\mathbf{H}}}= \big(

smallmatrix\tilde{\mathbf{H}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_p

\big)$. \@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{New Bootstrap Procedure}

Now consider bootstrapping $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}})$. Following jiang2023revisiting, it is natural to regard the PC estimators as estimators of the signal parameters in the pseudo-true models (ref) and (ref). We adopt these pseudo-true models in the bootstrap resampling because the PC parameters $(\hat{\mathbf{F}},\hat{\mathbf{B}})$, which serve as the “true” parameters in the bootstrap world, satisfy the same $r^2$ restrictions as $(\mathbf{F}^0,\mathbf{B}^0)$. The novelty of our bootstrap procedure is the use of the key identity (ref) to generate the rotation-dependent parameter vector $\boldsymbol{\delta}_{\hat{\mathbf{H}}}$ for the pseudo-true models.

We describe the bootstrap procedure for approximating the distribution of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}})$ as follows. Variables generated under the bootstrap law are denoted by the superscript “$\dag$”. The superscript $(b)$ refers to the $b$-th bootstrap sample, for $b=1,\dots,B$.

enumerate• Generate the bootstrap data $\mathbf{X}^{\dag(b)}$ using a resampled error matrix $\mathbf{E}^{\dag(b)}$: \begin{align} \mathbf{X}^{\dag(b)}=\mathbf{F}^{0\dag}\mathbf{B}^{0\dag\prime} + \mathbf{E}^{\dag(b)} \end{align} where $\mathbf{F}^{0\dag}:=\hat{\mathbf{F}}$ and $\mathbf{B}^{0\dag}:=\hat{\mathbf{B}}$. Using $\mathbf{X}^{\dag(b)}$, obtain the bootstrap PC estimators $\hat{\mathbf{F}} ^{\dag(b)}$ and $\hat{\mathbf{B}}^{\dag(b)}$, correcting their signs if necessary so that all sample correlations $cor(\hat{f}_{k,t}^{\dag(b)},f_{k,t}^{0\dag})$, $k=1,2,\dots,r$, are positive. Construct the bootstrap rotation matrix \begin{align} \tilde{\mathbf{H}}^{\dag(b)} =\mathbf{B}^{0\dag\prime}\mathbf{B}^{0\dag}(T^{-1}\mathbf{F}^{0\dag}{\hat{\mathbf{F}}}^{\dag(b)})(\hat{\mathbf{B}}^{\dag(b)\prime}\hat{\mathbf{B}}^{\dag(b)})^{-1}. \end{align} • Generate the bootstrap data $\mathbf{y}^{\dag(b)}$ using a resampled error vector $\boldsymbol{\epsilon}^{\dag(b)}$: \begin{align} \mathbf{y}^{\dag(b)}=\mathbf{F}^{0\dag}{\boldsymbol{\gamma}^{0\dag}} + \mathbf{W}\boldsymbol{\beta}^{\dag} +\boldsymbol{\epsilon}^{\dag(b)} =\mathbf{Z}^{0\dag}\boldsymbol{\delta}^{0\dag} + \boldsymbol{\epsilon}^{\dag(b)} \end{align} where ${\boldsymbol{\gamma}^{0\dag}}:=\hat{\boldsymbol{\gamma}}$, $\boldsymbol{\beta}^{\dag}:=\hat{\boldsymbol{\beta}}$, $\mathbf{Z}^{0\dag} = (\mathbf{F}^{0\dag},\mathbf{W})$ and $\boldsymbol{\delta}^{0\dag}=(\boldsymbol{\gamma}^{0\dag\prime},\boldsymbol{\beta}^{\dag\prime})'$. Using $\hat{\mathbf{Z}}^{\dag(b)} = (\hat{\mathbf{F}}^{\dag(b)},\mathbf{W})$, compute the bootstrap estimator $\hat{\boldsymbol{\delta}}^{\dag(b)}=(\hat{\mathbf{Z}}^{\dag(b)\prime}\hat{\mathbf{Z}}^{\dag(b)})^{-1}\hat{\mathbf{Z}}^{\dag(b)\prime}\mathbf{y}^{\dag(b)}$ and the corresponding bootstrap “parameter vector” $\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}}=(\boldsymbol{\gamma}^{0\dag\prime}\tilde{\mathbf{H}}^{\dag(b)\prime -1},\boldsymbol{\beta}^{\dag\prime})'=\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag(b)}}^{-1} \boldsymbol{\delta}^{0\dag}$ with $\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag(b)}}= \big( \begin{smallmatrix} \tilde{\mathbf{H}}^{\dag(b)} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_p \end{smallmatrix} \big)$ to obtain \begin{align} \sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag(b)} - \boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}}). \end{align} • Repeat Steps 1-2 for $b=1,2,\dots,B$ to construct the bootstrap distribution of (ref).

Note that the asymptotic justification for using the bootstrap statistic (ref) to mimic $\sqrt{T}(\hat{\boldsymbol{\delta}} - \boldsymbol{\delta}_{\hat{\mathbf{H}}})$ relies on two facts, which are overlooked in the literature: (i) the pseudo-true model is the unique (up to sign) transformation of the latent model that the PC estimators recover, and (ii) the identity ${\tilde{\mathbf{H}}}^{-1}\boldsymbol{\gamma}^0={\hat{\mathbf{H}}}^{-1}\gamma^*$ holds.

remGiven the bootstrap errors $(\mathbf{E}^{\dag},\boldsymbol{\epsilon}^{\dag})$, our bootstrap procedure differs from that of GoncalvesPerron2014,gonccalves2020bootstrapping in the way the objective statistics of interest are computed. Instead of using (ref), GoncalvesPerron2014 suggest computing $\sqrt{T}(\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}\hat{\boldsymbol{\delta}}^{\dag} - \boldsymbol{\delta}^{0\dag})=\sqrt{T}\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}(\hat{\boldsymbol{\delta}}^{\dag} - \boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}})$. Although exempt from the sign indeterminacy, even if $\tilde{\boldsymbol{\Phi}}^{\dag}\stackrel{p^{\dag}}{\longrightarrow}\mathbf{I}_{r+p}$, their bootstrap introduces additional randomness into $\sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag} - \boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}})$ through pre-multiplication by $\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}$. Therefore, our bootstrap procedure is expected to mimic the distribution of the objective statistic more efficiently in finite samples.

We can consider various bootstrap resampling methods for the elements of $\mathbf{E}^{\dag}$ and $\boldsymbol{\epsilon}^{\dag}$. To account for heteroskedastic errors, we employ the wild bootstrap, defined as $\mathbf{E}^{\dag}=(s_{t,i}^{\dag}\hat{e}_{t,i})$ and $\boldsymbol{\epsilon}^{\dag}=(\omega_{t}^{\dag}\hat{\epsilon}_{t})$, where $s_{t,i}^{\dag}$ and $\omega_{t}^{\dag}$ are i.i.d. random variables satisfying $\E^{\dag}[{s_{t,i}^{\dag}}]=0$, $\E^{\dag}[{s_{t,i}^{\dag2}}]=1$, $\E^{\dag}[{\omega_{t}^{\dag}}]=0$ and $\E^{\dag}[{\omega_{t}^{\dag 2}}]=1$. For bootstrap procedures designed to handle cross-correlated errors, or errors that are both cross- and serially correlated, see gonccalves2020bootstrapping and LiShenZhou2024.

The proposed procedure can be applied in various contexts, including asymptotic bias approximation, confidence interval construction, and hypothesis testing, under different choices of rotation matrices, as described next.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Bootstrapping for Different Rotation Matrices}

As shown by BaiNg2023 and jiang2023revisiting, there exist several rotation matrices other than $\hat{\mathbf{H}}$. In particular, jiang2024Mw consider the approximation and the identity

align[align omitted — 163 chars of source]

respectively, where $\hat{\mathbf{H}}_q=(T^{-1}\hat{\mathbf{F}}'\mathbf{F}^*)^{-1}$ and $\tilde{\mathbf{H}}_q=(T^{-1}\hat{\mathbf{F}}'\mathbf{F}^0)^{-1}$. Comparing this to (ref), we see that $\mathbf{H}$ is also the population analogue of $\hat{\mathbf{H}}_q$. jiang2024Mw further show that $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ is asymptotically normal, with an asymptotic bias generally different from that of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}})$. The approximate distribution of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ can be obtained using the same bootstrap procedure as before, but replacing $\tilde{\mathbf{H}}^{\dag(b)}$ in (ref) with $\tilde{\mathbf{H}}_q^{\dag(b)}=(T^{-1}\hat{\mathbf{F}}^{\dag(b)}\mathbf{F}^{0\dag})^{-1}$, and replacing $\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}}$ in (ref) with $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag(b)}}=(\boldsymbol{\gamma}^{0\dag\prime}\tilde{\mathbf{H}}_q^{\dag(b)\prime -1},\boldsymbol{\beta}^{\dag\prime})'$.

As argued in jiang2024Mw, it is natural to regard $(\boldsymbol{\delta}^0,{\mathbf{F}}^0)$ as the parameters estimated by $(\hat{\boldsymbol{\delta}},\hat{\mathbf{F}})$. In this context, the distribution of interest is $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$. To bootstrap this distribution, the same procedure described above can be used, with $\tilde{\mathbf{H}}^{\dag(b)}$ in (ref) replaced by $\mathbf{I}_r$ and $\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}}$ in (ref) replaced by $\boldsymbol{\delta}^{0\dag} (:= \hat{\boldsymbol{\delta}})$.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Theory} In this section, we establish the asymptotic validity of the proposed bootstrap procedure in approximating the distribution of the estimator $\hat{\boldsymbol{\delta}}$ relative to the rotated parameter vectors under different rotation matrices. \@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Assumptions} We begin with the assumptions underlying the non-bootstrap results, followed by the additional assumptions required for the bootstrap analysis. Assumptions (ref)--(ref) pertain to the non-bootstrap results and are identical to those in jiang2024Mw.

ass\normalfont The smallest eigenvalues of ${\mathbf{B}^*}'\mathbf{B}^*$ and $T^{-1}{\mathbf{F}^*}'\mathbf{F}^*$ are bounded away from zero.
ass[Signal strength]\normalfont There exist random or non-random variables $d_1,\dots,d_r>0$ and constants $0<\alpha_r\leq \dots \leq \alpha_1 \leq 1$ such that {$\lambda_k=d_kN^{\alpha_k}$} for $k=1,\dots, r$ with ordered $0<\lambda_r< \dots < \lambda_1$ for large $N$. If $d_k$'s are random, we have {$\E[d_k^2] \le M$} for all $k$.

Denote $\mathbf{N}=\diag(N^{\alpha_1}, \dots,N^{\alpha_r})$ and $\mathbf{D}=\diag(d_1,\dots,d_r)$, so that we can write {$\boldsymbol{\Lambda}=\mathbf{D} \mathbf{N}$}. Note that we do not require any specific structure in $(\mathbf{F}^{*}, \mathbf{B}^{*})$, such as diagonality of $\mathbf{N}^{-\frac{1}{2}}\mathbf{B}^{*\prime}\mathbf{B}^*\mathbf{N}^{-\frac{1}{2}}$ in BaiNg2023 and/or $T^{-1}\mathbf{F}^{*\prime}\mathbf{F}^*=\mathbf{I}_r$ in Freyaldenhoven21JoE.

ass[Idiosyncratic errors]\normalfont ${ }^{ }$\\ (i) $\E[e_{t,i}]=0$ and $\E[e_{t,i}^4]\le M$ for all $i$ and $t$;\\ (ii) For all $i$, $\left|\E[e_{s,i} e_{t,i}]\right| \leq\left|\gamma_{s, t}\right|$ for some $\gamma_{s, t}$ such that $\sum_{t=1}^T\left|\gamma_{ s, t}\right| \leq M$;\\ (iii) For all $t$, $\left|\E[e_{t,i} e_{t,j}]\right| \leq\left|\tau_{i, j}\right|$ for some $\tau_{i, j}$ such that $\sum_{j=1}^N\left|\tau_{ i,j}\right| \leq M$;\\ (iv) $\left\|\mathbf{E} \right\|_2^2=O_p(\max{\{N,T}\})$.

As discussed earlier, the PC estimators $(\hat{\mathbf{F}},\hat{\mathbf{B}})$ are viewed as estimators of the pseudo-true parameters $({\mathbf{F}^0},{\mathbf{B}^0})$. Accordingly, we impose the following assumptions directly on them.

ass[Factors and Loadings]\normalfont Denote $\mathbf{z}^0_t=(\mathbf{f}^{0\prime}_t, \mathbf{w}_t')'$. \\ (i) $\E\|{\mathbf{z}_t^0}\|_2^4 \le M$ and $\E\|\mathbf{b}_i^0\|_2^4 \le M$;\\ (ii) $\E\| \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}^0_ie_{t,i} \|_2^2 \le M$ for each $t$; \\ (iii) $\E\| T^{-\frac{1}{2}}\sum_{t=1}^{T}{\mathbf{z}_t^0}e_{t,i}\|_2^2 \le M$ for each $i$; \\ (iv) The $r \times r$ matrix satisfies $\E\| T^{-\frac{1}{2}}\mathbf{N}^{-\frac{1}{2}} \sum_{t=1}^{T}\sum_{i=1}^{N}\mathbf{b}_i^{0}e_{t,i}{\mathbf{z}_t^{0\prime}} \|_{{2}}^2 \le M$;\\ (v) $T^{-1} \sum_{t=1}^{T} ( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i}) ( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i})'\stackrel{p}{\longrightarrow} \boldsymbol{\Gamma} $, where $\boldsymbol{\Gamma} = \lim_{N, T \rightarrow \infty} T^{-1} \sum_{t=1}^T \boldsymbol{\Gamma}_t>0$, and $\boldsymbol{\Gamma}_t = \operatorname{Var}(\mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i})$.

The moment restrictions in Assumption 4 (iii), (iv) are essentially similar to Assumptions D, F2 in Bai2003, and Assumption 4 (ii) is similar moment restriction related for $\mathbf{b}_i^0$. Assumption (v) is similar to Assumption 3(e) in GoncalvesPerron2014.

Now we impose assumptions on the pseudo-true augmented model (ref):

ass[Weak dependence between idiosyncratic errors and regression errors]\normalfont The $r \times r$ matrix satisfies $\E||T^{-\frac{1}{2}}\mathbf{N}^{-\frac{1}{2}} \sum_{t=1}^{T} \sum_{i=1}^N \mathbf{b}_i^0 e_{t,i} \varepsilon_{t+h}||_2^2 \leq M$.

{

ass[Moments, parameters and CLT] \normalfont${ }^{ }$\\ (i) $\E[\epsilon_{t+h}] =0$ and $\E|\epsilon_{t+h}|^2 < M $;\\ (ii) $||\boldsymbol{\delta}^0||_{2}\leq M$ and $\mathbf{H}\stackrel{p}{\longrightarrow} \mathbf{H}_0$ which is fixed and invertible;\\ (iii) $\E||\mathbf{z}_t^0 ||^4 \leq M$, $T^{-1/2} {\mathbf{Z}}^{0\prime}\boldsymbol{\epsilon} \stackrel{d}{\longrightarrow} N(\mathbf{0},\boldsymbol{\Sigma}_{\mathbf{Z}^0 \boldsymbol{\epsilon}})$, $T^{-1} {\mathbf{Z}}^{0\prime}\mathbf{Z}^0 \stackrel{p}{\longrightarrow} \boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}$, where $\boldsymbol{\Sigma}_{\mathbf{Z}^0\boldsymbol{\epsilon}}$ and $\boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}$ are fixed, positive definite and bounded.

} Assumptions (ref) and (ref) are similar to Assumption 4 in GoncalvesPerron2014 and Assumption E in BaiNg2006, respectively. Under Assumption (ref), $\mathbf{H}$ is bounded in probability. Assumption (ref)(ii) further guarantees that its probability limit exists and coincides with that of other four data-dependent rotation matrices considered in jiang2023revisiting.

We now state the assumptions required for the bootstrap analysis, denoted by superscripts “$\dag$” in the assumption numbers.

remThroughout the paper, $\Pr^{\dag}$, $\E^{\dag}$ and $\operatorname{Var^{\dag}}$ denote probability, expectation and variance, conditional on the realization of the original sample, respectively. We use the symbols $o_{p^{\dag}}$ and $O_{p^{\dag}}$ for bootstrap sample asymptotics, which correspond to $o_p$ and $O_p$ for the original sample asymptotics.

\setcounter{assb}{2}

assb[Idiosyncratic errors] \normalfont${ }^{ }$\\ {(i) $\E^{\dag}[e_{t,i}^{\dag}]=0$ and $\E^{\dag}[e_{t,i}^{\dag4}]=O_p(1)$ for all $i$ and $t$;}\\ (ii) For all $i$, $|\E^{\dag}[e_{s,i}^{\dag} e_{t,i}^{\dag}]| \leq|\gamma_{s, t}^{\dag}|$ for some $\gamma_{s, t}^{\dag}$ such that $\sum_{t=1}^T|\gamma_{ s, t}^{\dag}| =O_p(1)$;\\ (iii) For all $t$, $|\E^{\dag}[e_{t,i}^{\dag} e_{t,j}^{\dag}]| \leq|\tau_{i, j}^{\dag}|$ for some $\tau_{i, j}^{\dag}$ such that $\sum_{j=1}^N|\tau_{ i,j}^{\dag}| =O_p(1)$;\\ (iv) $||\mathbf{E}^{\dag} ||_2^2=O_{p^{\dag}}(\max{\{N,T}\})$, in probability.
assb[Factors and Loadings]\normalfont ${ }^{ }$\\ (i) $\E^{\dag}\| \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \hat{\mathbf{b}}_ie_{t,i}^{\dag} \|_2^2 =O_p(1)$ for each $t$; \\ (ii) $\E^{\dag}\| T^{-\frac{1}{2}}\sum_{t=1}^{T}\hat{\mathbf{z}}_t e_{t,i}^{\dag}\|_2^2 =O_p(1)$ for each $i$; \\ (iii) The $r \times r$ matrix satisfies $\E^{\dag}\| T^{-\frac{1}{2}}\mathbf{N}^{-\frac{1}{2}} \sum_{t=1}^{T}\sum_{i=1}^{N} \hat{\mathbf{b}}_i e_{t,i}^{\dag}\hat{\mathbf{z}}_t^{\prime} \|_{2}^2 =O_p(1)$; \\ (iv) $T^{-1} \sum_{t=1}^{T} ( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \hat{\mathbf{b}}_i e_{t,i}^{\dag}) ( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \hat{\mathbf{b}}_i e_{t,i}^{\dag})'- \boldsymbol{\Gamma}^{\dag} =o_{p^{\dag}}(1) $, in probability, where $\boldsymbol{\Gamma}^{\dag} = T^{-1} \sum_{t=1}^T \operatorname{Var^{\dag}}(\mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \hat{\mathbf{b}}_i e_{t,i}^{\dag} )$ is positive definite almost surely.
assb[Weak dependence between idiosyncratic errors and regression errors] The $r \times r$ matrix satisfies $\E^{\dag} ||T^{-\frac{1}{2}}\mathbf{N}^{-\frac{1}{2}} \sum_{t=1}^{T} \sum_{i=1}^N \hat{\mathbf{b}}_i e_{t,i}^{\dag} \varepsilon_{t+h}^{\dag}||_2^2 =O_p(1)$.
assb[Moments and CLT for the Score Vector]\normalfont ${ }^{ }$\\ (i) $\E[\epsilon_{t+h}^\dag]=0$ and $T^{-1}\sum_{t=1}^T\E|\epsilon_{t+h}^{\dag}|^2=O_p(1)$;\\ (ii) $\boldsymbol{\Sigma}_{\hat{\mathbf{Z}} \boldsymbol{\epsilon}^{\dag}}^{-1/2} T^{-1/2} \hat{\mathbf{Z}}^{\prime}\boldsymbol{\epsilon}^{\dag} \stackrel{d^{\dag}}{\longrightarrow} N(\mathbf{0}, \mathbf{I}_{(r+p)})$, in probability, where $ \E^{\dag} || T^{-\frac{1}{2}}\sum_{t=1}^{T} \hat{\mathbf{z}}_t \epsilon_{t+h}^{\dag} ||_2^2 =O_p(1)$, and $\boldsymbol{\Sigma}_{\hat{\mathbf{Z}}\boldsymbol{\epsilon}^{\dag}}=\operatorname{Var^{\dag}}(T^{-\frac{1}{2}}\sum_{t=1}^{T} \hat{\mathbf{z}}_t \epsilon_{t+h}^{\dag} )$ is positive definite almost surely.
assb\normalfont $\operatorname*{plim} \boldsymbol{\Sigma}_{\hat{\mathbf{Z}} \boldsymbol{\epsilon}^{\dag}} = \boldsymbol{\Sigma}_{\mathbf{Z}^0 \boldsymbol{\epsilon} }$ and $\operatorname*{plim} \boldsymbol{\Gamma}^{\dag} =\boldsymbol{\Gamma}$.

Assumptions (ref)--(ref) are the bootstrap analogues of Assumption (ref)--(ref). Assumption (ref) is similar to Conditions E* and F* in GoncalvesPerron2014, which guarantees the consistency of the bootstrap, so that the relevant bootstrap and original statistics converge in probability to the same quantities. Since $\hat{\mathbf{z}}_t$ estimates $ \mathbf{z}_t^0$, $\boldsymbol{\Sigma}_{\hat{\mathbf{Z}} \boldsymbol{\epsilon}^{\dag}}$ is the sample analogue of $ \boldsymbol{\Sigma}_{\mathbf{Z}^0 \boldsymbol{\epsilon} } $ provided that $\epsilon_{t+h}^{\dag}$ is constructed to mimic the time series dependence of $\epsilon_{t+h}$. By Assumption (ref)(iv), $\boldsymbol{\Gamma}^{\dag} = T^{-1} \sum_{t=1}^T \operatorname{Var}^{\dag}(\mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \hat{\mathbf{b}}_i e_{t,i}^{\dag})$. Since $\hat{\mathbf{b}}_i$ estimates $ \mathbf{b}_i^0$, $\boldsymbol{\Gamma}^{\dag}$ is the sample analogue of $ \boldsymbol{\Gamma} $ if $e_{t,i}^{\dag}$ is constructed to mimic the cross-sectional dependence of $e_{t,i}$.

Given these assumptions, we now present our main theoretical results. \@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Main Results}

The Case of $\hat{\mathbf{H}}$

thmSuppose Assumptions (ref)–(ref) and (ref)--(ref) hold. If $\alpha_r>\frac{1}{2}$, $\frac{N^{1-\alpha_r}}{\sqrt{T}} \to 0$, and $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r} \to c_1 \in [0,\infty)$, as $N, T \to \infty$, we have \[ \sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}}) \stackrel{d^{\dag}}{\longrightarrow} N\left(-c_1 \boldsymbol{\kappa}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right) , \] in probability, with \begin{align*} \boldsymbol{\kappa}_{\boldsymbol{\delta}^*}= \boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}^{-1} \binom{\mathbf{G} + \nu\mathbf{D}^{-1}\boldsymbol{\Gamma}\mathbf{D}^{-1}}{\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0} \, {\mathbf{G}}} \, \mathbf{H}_0^{-1} \, \boldsymbol{\gamma}^* and \boldsymbol{\Sigma}_{\boldsymbol{\delta}}= \boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}^{-1} \boldsymbol{\Sigma}_{\mathbf{Z}^0 \boldsymbol{\epsilon}} \boldsymbol{\Sigma}_{\mathbf{Z}^0\mathbf{Z}^0}^{-1}, \end{align*} where $c_1 \mathbf{G} = \lim_{N,T\to\infty} \sqrt{T}\mathbf{N}^{\frac{1}{2}} \boldsymbol{\Gamma} \mathbf{D}^{-2} \mathbf{N}^{-\frac{3}{2}} $, $\nu = \lim_{N\to\infty} N^{-\frac{1}{2}(\alpha_1-\alpha_r)}$ and\\ $\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0} = \operatorname*{plim}_{N,T \rightarrow \infty} T^{-1}\mathbf{W}'\mathbf{F}^0$.
remTogether with the non-bootstrap asymptotic normality results established in jiang2024Mw, $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})\stackrel{d}{\longrightarrow} N\left(-c_1 \boldsymbol{\kappa}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right)$ under the same conditions, it is straightforward to show the bootstrap validity: \begin{align} \sup_{\mathbf{x}\in \mathbb{R}^{r+p}}{|{\Pr}^{\dag}[\sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}})<\mathbf{x}]-\Pr[\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})<\mathbf{x}]|}=o_p(1). \end{align} As analyzed in jiang2024Mw, the expression of the asymptotic bias suggests a complicated asymptotic bias structure, depending on the structure of the divergence rates, $(\alpha_1,\ldots,\alpha_r)$. Nonetheless, the result shows that the bootstrap procedure can mimic the non-central distribution of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ without requiring knowledge of the divergence rates, provided the conditions are satisfied.
remAs discussed in Remark (ref), GoncalvesPerron2014,gonccalves2020bootstrapping use the bootstrap statistic $\sqrt{T}(\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}^{0\dag}) = \sqrt{T}\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}})$ in our notation to mimic the distribution of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$; see Corollary 3.1 in GoncalvesPerron2014. Since $\tilde{\mathbf{H}}^{\dag}-\mathbf{I}_r=o_{p^{\dag}}(1)$, the asymptotic bootstrap validity of their procedure, $\sup_{\mathbf{x}\in \mathbb{R}^{r+p}} \big|{\Pr}^{\dag}[\sqrt{T}\boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag}}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}})<\mathbf{x}] -\Pr[\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})<\mathbf{x}]\big| = o_p(1)$, follows immediately. This is consistent with their result, though our framework clarifies that their procedure requires an additional estimation of $\mathbf{I}_r$ via $\tilde{\mathbf{H}}^{\dag}$.
remThere are three main differences between Theorem (ref) and the result in GoncalvesPerron2014 for SF models, which, ignoring sign indeterminacy, essentially states that $\sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}^{0\dag}) \stackrel{d^{\dag}}{\longrightarrow} N\!\left(-c_1 \boldsymbol{\kappa}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right)$, with $\boldsymbol{\delta}^{0\dag}:=\hat{\boldsymbol{\delta}}$ in our notation. First, their framework essentially treats $\mathbf{H}_0^{-1}{\boldsymbol{\gamma}}^*$, where $\mathbf{H}_0 = \operatorname*{plim}_{N,T\rightarrow\infty} \hat{\mathbf{H}}$, as the parameter of interest. By contrast, we consider $\boldsymbol{\gamma}^0:=\mathbf{H}^{-1}{\boldsymbol{\gamma}}^*$ for finite $\{N,T\}$. We therefore do not construct a bootstrap analogue of $\mathbf{H}_0$, as it has no counterpart in the original sample. Second, consider a decomposition $\sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}^{0\dag}) = \sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}}) + \sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}}-\boldsymbol{\delta}^{0\dag}).$ For their result to hold, it is necessary that $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}}-\boldsymbol{\delta}^{0\dag}) = o_{p^\dag}(1)$, but this appears to hold only under restrictive conditions.\footnote{For SF models, it can be shown that $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag}}-\boldsymbol{\delta}^{0\dag})=\sqrt{T}O_{p^{\dag}}(1/\min\{N,T\}))$, and this term does not necessarily vanish in probability when $\sqrt{T}/N \to c \in [0,\infty)$, which is precisely the condition under which a nonzero asymptotic bias may arise.} Third, a similar decomposition can be obtained via one of the alternative rotation matrices, such as $\hat{\mathbf{H}}_q^{\dag}$. However, as shown in Theorem (ref) below, the corresponding asymptotic bias differs from $-c_1 \boldsymbol{\kappa}_{\boldsymbol{\delta}^*}$, which introduces ambiguity in interpreting their result. We address these three points in Section (ref), where we derive the asymptotic distribution of $\sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}^{0\dag})$.

The Case of $\hat{\mathbf{H}}_q$

We next present the result for the rotated parameter vector under the alternative data-dependent rotation matrix $\hat{\mathbf{H}}_q$.

thmSuppose Assumptions (ref)–(ref) and (ref)--(ref) hold. If $\alpha_r>\frac{1}{2}$, $\frac{N^{1-\alpha_r}}{\sqrt{T}} \to 0$, and $\sqrt{T}N^{-\alpha_r} \to c_2 \in [0,\infty)$, as $N, T \to \infty$, we have \[ \sqrt{T}(\hat{\boldsymbol{\delta}}^{\dag}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag}}) \stackrel{d^{\dag}}{\longrightarrow} N\left(c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right) \] in probability, with \begin{align*} \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*} = \boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}^{-1} \binom{\mathbf{0}}{\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0} \bar{\mathbf{G}} } \mathbf{H}_0^{-1} \boldsymbol{\gamma}^*, \end{align*} where $c_2 \bar{\mathbf{G}} = \lim_{N, T \rightarrow \infty}\sqrt{T}\mathbf{N}^{-\frac{1}{2}} \mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}\mathbf{N}^{-\frac{1}{2}} $.

Again, together with the non-bootstrap asymptotic normality results for $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})\stackrel{d}{\longrightarrow} N\left(c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right)$, established in jiang2024Mw under the same conditions, it is straightforward to establish the bootstrap validity.

Theorem (ref) suggests that, in general, both $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ and its bootstrap counterpart converge to their limiting distribution faster and exhibit smaller asymptotic bias than $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ and its bootstrap analogue. Moreover, they are asymptotically unbiased when $\mathbf{F}^0$ and $\mathbf{W}$ are (asymptotically) uncorrelated. Therefore, for bootstrap asymptotic analysis in factor-augmented regressions, it is preferable to adopt the approximation $\hat{\mathbf{F}} = \mathbf{F}^* \hat{\mathbf{H}}_q + o_p(1)$ rather than $\hat{\mathbf{F}} = \mathbf{F}^* \hat{\mathbf{H}} + o_p(1)$, in both the bootstrap and original samples.

The Case of ${\mathbf{H}}$

Now let us investigate the bootstrap analogue of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$. Since $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ typically exhibits more favorable asymptotic properties than $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$, and is in fact the most favorable among the asymptotically equivalent rotation matrices considered in BaiNg2023, it is natural to consider the decomposition $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0) = \sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q}) + \sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^0)$. The first term has already been analyzed. For the second term, we obtain $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^0)= \big(

smallmatrix\sqrt{T}(\tilde{\mathbf{H}}_q^{-1}-\mathbf{I}_r)\mathbf{H}^{-1}\boldsymbol{\gamma}^*\\ \mathbf{0}_p

\big) =O_p(\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r})$, but no explicit bias expression is available. Following \cite{jiang2024Mw}, we assume that $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^0)$ converges in probability to a bounded constant vector, say $c_1\mathbf{h}_{\boldsymbol{\gamma}^*}$, when $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r} \to c_1 \in [0,\infty)$. For the bootstrap counterpart, we have $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag}}-\boldsymbol{\delta}^{0\dag})= \big(

smallmatrix\sqrt{T}(\tilde{\mathbf{H}}_q^{{\dag}-1}-\mathbf{I}_r)\hat{\boldsymbol{\gamma}}\\ \mathbf{0}_p

\big) =O_{p^\dag}(\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r})$, in probability. Since $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag}}-\boldsymbol{\delta}^{0\dag}$ is the bootstrap analogue of $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q} -\boldsymbol{\delta}^{0}$, we impose the analogous assumption that $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag}}-\boldsymbol{\delta}^{0\dag}) \stackrel{p^{\dag}}{\longrightarrow} c_1 {\mathbf{h}_{\boldsymbol{\gamma}}^\dag}$ in probability, where $\mathbf{h}_{\boldsymbol{\gamma}}^\dag$ depends on $\hat{\mathbf{F}}$ and $\hat{\boldsymbol{\gamma}}$. To ensure the bootstrap asymptotic validity, analogously to Assumption $7^\dag$, we further assume that $\operatorname*{plim}\mathbf{h}_{\boldsymbol{\gamma}}^\dag =\mathbf{h}_{\boldsymbol{\gamma}^*}$. This discussion leads to the following assumption.

assb$\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^{0}) \stackrel{p}{\longrightarrow} c_1 {\mathbf{h}_{\boldsymbol{\gamma}^*}}$, where $\mathbf{h}_{\boldsymbol{\gamma}^*}$ is a constant vector with its last $p$ entries equal to zero and $||\mathbf{h}_{\boldsymbol{\gamma}^*}||_2 \leq M$. In addition, $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q^{\dag}}-\boldsymbol{\delta}^{0\dag}) \stackrel{p^{\dag}}{\longrightarrow} c_1 {\mathbf{h}_{\boldsymbol{\gamma}}^\dag}$ in probability, where $\mathbf{h}_{\boldsymbol{\gamma}}^\dag=O_p(1)$ and $\operatorname*{plim}\mathbf{h}_{\boldsymbol{\gamma}}^\dag=\mathbf{h}_{\boldsymbol{\gamma}^*}$.

We are now ready to present the results on the bootstrap analogue of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}^0)$.

thmSuppose that Assumptions (ref)–(ref) and (ref)--(ref) hold, and that $\alpha_r>\frac{1}{2}$, $\frac{N^{1-\alpha_r}}{\sqrt{T}} \to 0$, $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r} \to c_1 \in [0,\infty)$, $\sqrt{T}N^{-\alpha_r} \to c_2 \in [0,\infty)$. Then, we have \begin{align*} & \sqrt{T}(\hat\boldsymbol{\delta}^{\dag} - \boldsymbol{\delta}^{0^{\dag}}) \stackrel{d^{\dag}}{\longrightarrow} N\left(c_1 {\mathbf{h}_{\boldsymbol{\gamma}^*}} + c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*} ,\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right), \; in probability. \end{align*}

Building on the non-bootstrap asymptotic normality result $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}^0)\stackrel{d}{\longrightarrow} N(c_1 \mathbf{h}_{\boldsymbol{\gamma}^*} +c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*},\boldsymbol{\Sigma}_{\boldsymbol{\delta}})$ established in jiang2024Mw under the same conditions, the bootstrap validity follows immediately.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Monte Carlo Experiments}

In this section, we examine the finite sample performance of the estimators of the factor-augmented regressions.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Design} We generate data according to $\mathbf{X=F}^{0}\mathbf{B}^{0\prime}+\mathbf{E}$, $\mathbf{X}=(x_{t,i})$, $i=1,\dots,N$, $t=1,\dots,T$, $\mathbf{F}^{0}\in\mathbb{R} ^{T\times r}$, $\mathbf{B}^{0}\in\mathbb{R}^{N\times r}$ are constructed as follows. Define a positive definite matrix $\mathbf{D}=\diag(d_{1},\dots,d_{r})$ and $\mathbf{N}=\diag(N^{\alpha_{1}},\dots,N^{\alpha_{r}})$. Construct a $T\times N$ matrix $\mathbf{A}$ with elements drawn independently from $N(0,1)$ in each replication. Let $\mathbf{A=USV}^{\prime}$ be its singular value decomposition. Set $\mathbf{F}^{0}$ to the first $r$ columns of $\mathbf{U}$ multiplied by $\sqrt{T}$, and $\mathbf{B}^{0}$ to the first $r$ columns of $\mathbf{V}$ post-multiplied by $\mathbf{D}^{1/2}\mathbf{N}^{1/2}$. Given an invertible $r \times r$ $\mathbf{H}$, define $\mathbf{F}^{\ast}=\mathbf{F}^{0}\mathbf{H}^{-1}$ and $\mathbf{B}^{\ast}=\mathbf{B}^{0}\mathbf{H}^{\prime}$. For experiments we set $\mathbf{H}=\left(

smallmatrix1 & 1/2\\ 1/2 & 2

\right)$. The error matrix $\mathbf{E}$ is cross-sectionally heteroskedastic but independent over $t$. Specifically, the $t^{th}$ row is generated as $\mathbf{e}_{t}=\boldsymbol{\Sigma}_{e}^{1/2}\boldsymbol{\xi}_{t}$ where $\boldsymbol{\xi}_{t}\sim i.i.d.N(\mathbf{0},\mathbf{I}_{N})$, and $\boldsymbol{\Sigma} _{e}=diag(\sigma_{e1}^{2},...,\sigma_{eN}^{2})$ with $\sigma_{ei}^{2}\sim i.i.d.U[0.5,1.5]$, $i=1,2,...,N$. The factor-augmented regression is specified as \[ y_{t+1}=\mathbf{f}_{t}^{0\prime} \boldsymbol{\gamma}^0+\mathbf{w}_{t} ^{\prime}\boldsymbol{\beta}+\epsilon_{t+1}\text{, }t=1,\dots,T, \] where $\mathbf{f}_{t}^{0\prime}$ is the $t^{th}$ row of $\mathbf{F}^{0}$, $\mathbf{w}_{t}=(w_{t,1},\dots,w_{t,p})^{\prime}$ with $w_{t,p}=1$, and $w_{t,\ell}=\sigma_{w}[\rho_{fw}\mathbf{f}_{t}^{0\prime}\mathbf{1}_{r} r^{-1/2}+(1-\rho_{fw}^{2})^{1/2}\zeta_{t,\ell}]$, with $\zeta_{t,\ell}\sim i.i.d.N(0,1\mathbf{)}$ for $\ell=1,\dots,p-1$, and $\epsilon_{t+1}\sim i.i.d.N(0,\sigma_{\epsilon}^{2})$. We set $\boldsymbol{\gamma}^0=\mathbf{1}_{r}$ and $\boldsymbol{\beta}=\mathbf{1}_{p}$, so that $\boldsymbol{\gamma}^{\ast}=\mathbf{H}\boldsymbol{\gamma}^0$.

As implied by the theory, the correlation between $\mathbf{w}_{t}$ and $\mathbf{f}_{t}$ affects the asymptotic bias of the estimator. We consider $\rho_{fw}=\{0,0.6\}$ while setting $\sigma_{w}^{2}=1$ and $\sigma_{\epsilon}^{2}=0.5$. We choose $r=2$ and $p=2$, and consider three factor models with different strengths: $(\alpha_{1},\alpha_{2})=(1,1)$, $(1,0.8),$ $(0.8,0.6)$, with $(d_{1} ,d_{2})=(0.05,0.2)$, $(0.2,0.2)$ and $(0.2,0.2)$, respectively. Different values of $d_{1}$ and $d_{2}$ are required in the case of $\alpha_{1}=\alpha_{2}$ to ensure identification of the two largest eigenvalues of $\E[\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}]$, denoted $\lambda_{1}$ and $\lambda_{2}$.

Suppose $\{\mathbf{x}_{t},\mathbf{w}_{t},\mathbf{y}_{t}\}$ are observable in practice. Then $\mathbf{F}^{0}$ is estimated by principal components, taken as the $r$ eigenvectors of $T^{-1}\mathbf{XX}^{\prime}$ corresponding to its $r$ largest eigenvalues, multiplied by $\sqrt{T}$. The resulting PC estimator of $\mathbf{F}^{0}$ is denoted by $\mathbf{\hat{F}}$. If necessary, the column signs of $\hat{\mathbf{F}}$ are adjusted so that all sample correlations $cor(\hat{f}_{k,t}, f_{k,t}^0)$, $k=1,2,\dots,r$, are positive.

The factor-augmented model is then estimated by regressing $y_{t+1}$ on $\mathbf{\hat{z}} _{t}=(\mathbf{\hat{f}}_{t}^{\prime},\mathbf{w}_{t}^{\prime})^{\prime}$, yielding $\boldsymbol{\hat{\delta}}=(\boldsymbol{\hat{\gamma} }^{\prime},\boldsymbol{\hat{\beta}}^{\prime})^{\prime}$. Using different rotation matrices, we evaluate the biases of the least squares estimators relative to alternative `parameter vectors'. Specifically, we compute the averages across replications of $\hat{\boldsymbol{\delta}}-{\boldsymbol{\delta}}_{\hat{\mathbf{H}}}$, $\hat{\boldsymbol{\delta}}-{\boldsymbol{\delta}}_{\hat{\mathbf{H}}_q}$ and $\hat{\boldsymbol{\delta}}-{\boldsymbol{\delta}}^0$.

In addition, we consider associated bootstrap bias-corrected estimators, defined as

align[align omitted — 372 chars of source]

where $\hat{\boldsymbol{b}}_{\hat{\mathbf{H}}}=B^{-1}\sum_{b=1}^{B} (\hat{\boldsymbol{\delta}}^{\dag(b)} - \boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}})$, $\hat{\boldsymbol{b}}_{{\hat{\mathbf{H}}}_q}=B^{-1}\sum_{b=1}^{B} (\hat{\boldsymbol{\delta}}^{\dag(b)} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_{q}^{\dag(b)}})$ and $\hat{\boldsymbol{b}}_{{\mathbf{H}}}=B^{-1}\sum_{b=1}^{B} (\hat{\boldsymbol{\delta}}^{\dag(b)} - \hat{\boldsymbol{\delta}})$. We report the biases of these estimators, namely $\hat{\boldsymbol{\delta}}_{bcb\hat{\mathbf{H}}} - \boldsymbol{\delta}_{\hat{\mathbf{H}}}$, $\hat{\boldsymbol{\delta}}_{bcb\hat{\mathbf{H}}_q} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$ and $\hat{\boldsymbol{\delta}}_{bcb} - \boldsymbol{\delta}^0$.

We also compare the performance of the proposed bootstrap algorithm with that of GoncalvesPerron2014, gonccalves2020bootstrapping. Their bias-corrected estimators are defined as

align[align omitted — 280 chars of source]

where $\hat{\boldsymbol{b}}_{GP\hat{\mathbf{H}}}=B^{-1}\sum_{b=1}^{B} \boldsymbol{\Phi}_{\tilde{\mathbf{H}}^{\dag (b)}}(\hat{\boldsymbol{\delta}}^{\dag(b)} - \boldsymbol{\delta}_{\hat{\mathbf{H}}^{\dag(b)}})$ and $\hat{\boldsymbol{b}}_{GP\hat{\mathbf{H}}_q}=B^{-1}\sum_{b=1}^{B} \boldsymbol{\Phi}_{\tilde{\mathbf{H}}_{q}^{\dag (b)}}(\hat{\boldsymbol{\delta}}^{\dag(b)} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_{q}^{\dag(b)}})$. The biases $\hat{\boldsymbol{\delta}}_{bcbGP\hat{\mathbf{H}}} - \boldsymbol{\delta}_{\hat{\mathbf{H}}}$ and $\hat{\boldsymbol{\delta}}_{bcbGP\hat{\mathbf{H}}_q} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$ are also reported.

In jiang2024Mw, instead of bootstrapping, the panel-split jackknife bias-corrected estimator for $\boldsymbol{\delta}^0$ is proposed. It is therefore of interest to compare its performance with that of the corresponding bootstrap procedure developed here. The panel-split jackknife estimator is defined as

align[align omitted — 301 chars of source]

where $\boldsymbol{\hat{\delta}}_{\mathcal{N}_{j}^{(s)}}$ is obtained regressing $\mathbf{y}$ on ($\mathbf{\hat{F}}_{\mathcal{N}_{j}^{(s)} },\mathbf{W}$), where $\mathbf{\hat{F}}_{\mathcal{N}_{j}^{(s)}}$ is the PC factor extracted from $\mathbf{X}_{\mathcal{N}_{j}^{(s)}}$, where $\mathbf{X}_{\mathcal{N}_{j}^{(s)}}=\{\mathbf{x}_{i\in\mathcal{N} _{j}^{(s)}}\}$, for $j=1,2$. Here, $\mathcal{N}_{1}^{(s)}$ and $\mathcal{N}_{2}^{(s)}$ denote the two halves of the $N$ columns of $\mathbf{X}^{(s)}$, which are randomly re-ordered in each replication $s=1,\dots,S$. Randomization helps avoid potentially biased information on the factors in $\mathcal{N}_j$. When $N$ is odd, $\mathcal{N}_{1}^{(s)}$ and $\mathcal{N}_{2}^{(s)}$ share one common index. The order and the sign of the columns of $\mathbf{\hat{F}}_{\mathcal{N}_{j}^{(s)}}$ are adjusted in line with those of $\mathbf{\hat{F}}$, based on the correlation between the pair $(\mathbf{\hat{F}}_{\mathcal{N}_{j}^{(s)}},\mathbf{\hat{F}})$, for each of $j=1,2$. The bias of the estimator, $\boldsymbol{\hat{\delta}}_{bcjk}-\boldsymbol{\delta}^0$, is reported.

The experiments are conducted for $(T,N)=(50,50),(100,100),(200,200)$ with 1,000 replications, $B=100$ and $S=100$.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Results} Figure (ref) plots the biases of the coefficient estimators for the second factor relative to their corresponding parameters. Biases are shown relative to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ (\textcolor{red}{red}), $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ (\textcolor{blue}{blue}), and $\boldsymbol{\gamma}^0$ (black). Thick solid lines denote uncorrected estimators; dashed lines indicate the jackknife bias-corrected estimator for $\boldsymbol{\gamma}^0$ (black) and the existing bootstrap bias-corrected estimator of GoncalvesPerron2014,gonccalves2020bootstrapping for $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ (\textcolor{red}{red}) and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ (\textcolor{blue}{blue}); short dotted lines correspond to the bootstrap bias-corrected estimators proposed here.

figure[figure omitted — 1,566 chars of source]

We begin with panels (a)–(c), where $\mathbf{f}_t$ and $w_t$ are uncorrelated (i.e., $\rho_{fw} = 0$). First, consider the red lines, which show biases relative to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$. The uncorrected estimator exhibits the largest bias, which worsens as the factor model weakens. Both bootstrap methods reduce this bias, with our proposed algorithm consistently outperforming that of GoncalvesPerron2014,gonccalves2020bootstrapping. Second, the blue lines correspond to biases relative to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$. Here, we observe a single flat line at zero, indicating that both $\hat{\boldsymbol{\gamma}}$ and the bias-corrected estimators are essentially unbiased with respect to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$. Third, the black lines plot biases relative to the true parameter $\boldsymbol{\gamma}^0$. For strong factor models, $\hat{\boldsymbol{\gamma}}$ shows negligible bias, but as the model weakens, a small bias emerges. Both the jackknife and the proposed bootstrap effectively correct for this.

Turning to panels (d)–(f), where $\mathbf{f}_t$ and $w_t$ are correlated (i.e., $\rho_{fw} = 0.8$), the biases of the uncorrected estimator $\hat{\boldsymbol{\gamma}}$ (solid lines) are consistently larger than in the uncorrelated case. Here, $\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ is no longer centered at zero, regardless of factor strength. Both bootstrap corrections reduce the bias, with our method again achieving greater reduction than GP’s. Bias properties relative to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ remain similar to the uncorrelated case. As before, the jackknife correction performs comparably to our bootstrap method.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}{Conclusion}

In this paper, we have proposed a novel bootstrap procedure that improves upon existing methods for replicating the asymptotic distribution of the factor-augmented regression estimator for a rotated parameter vector. The regression is augmented by $r$ factors extracted by the principal component (PC) method from a large panel of $N$ variables observed over $T$ time periods. We consider general weak factor (WF) models with $r$ signal eigenvalues that may diverge at different rates, $N^{\alpha _{k}}$, where $0<\alpha _{k}\leq 1$ for $k=1,2,...,r$.

We have established the asymptotic validity of our bootstrap method not only under the conventional data-dependent rotation matrix $\hat{\mathbf{H}}$, but also under an alternative data-dependent rotation matrix, $\hat{\mathbf{H}}_q$, which generally yields smaller asymptotic bias and achieves faster convergence. Moreover, we have shown bootstrap validity under a purely signal-dependent rotation matrix ${\mathbf{H}}$, which is unique and can be interpreted as the population analogue of both $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. This enables interpretation of the estimator's distribution relative to a parameter vector defined via ${\mathbf{H}}$, which is of practical importance.

While the asymptotic bias in WF models depends intricately on the structure of the divergence rates $(\alpha_1, \dots, \alpha_r)$, our results have shown that the bootstrap procedure can mimic the non-central distribution without requiring knowledge of these divergence rates.

Our theoretical contribution has also resolved a couple of ambiguities in the literature. First, existing approaches often define the parameter of interest via the probability limit of a data-dependent rotation matrix, $\mathbf{H}_0:=\operatorname*{plim}_{N,T\rightarrow\infty}\hat{\mathbf{H}}$, which is not observable or directly computable for bootstrapping. In contrast, we have proposed using a unique rotation matrix defined directly from the latent signal components at finite $\{N, T\}$ and constructible in bootstrap samples. Second, we have clarified the theoretical implications of using different data-dependent rotation matrices such as $\hat{\mathbf{H}}_q$, and highlighted the importance of properly accounting for the limiting behavior of $\sqrt{T}(\hat{\mathbf{H}} - \mathbf{H}_0)$ in establishing bootstrap validity.

One natural extension is to apply this bootstrap method to out-of-sample forecasting, where confidence intervals are often sensitive to normality assumptions. As emphasized in GodfreyOrme2000 and gonccalves2020bootstrapping, bootstrapping provides a practical alternative to such restrictive assumptions. Extending our approach to forecast evaluation remains an important direction for future research.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}*{Acknowledgment} We are grateful to Yoshimasa Uematsu for helpful discussions and useful comments.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\bf}*{Funding} This work was supported by JSPS KAKENHI (grant numbers 23H00804, 24K16343, 25K00625, 25K05036 and 25H00544).

{\setstretch{0.9} }