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.
110,535 characters
Bias Correction in Factor-Augmented Regression Models with Weak Factors
\maketitle
\begin{abstract}
In this paper, we study the asymptotic bias of the factor-augmented regression estimator and its reduction, which is augmented by the $r$ factors extracted from a large number of $N$ variables with $T$ observations. In particular, we consider general weak latent factor models with $r$ signal eigenvalues that may diverge at different rates, $N^{\alpha _{k}}$, $0<\alpha _{k}\leq 1$, $k=1,\dots,r$.
In the existing literature, the bias has been derived using an approximation for the estimated factors with a specific data-dependent rotation matrix $\hat{\mathbf{H}}$ for the model with $\alpha_{k}=1$ for all $k$, whereas we derive the bias for weak factor models.
In addition, we derive the bias using the approximation with a different rotation matrix $\hat{\mathbf{H}}_q$, which generally has a smaller bias than with $\hat{\mathbf{H}}$.
We also derive the bias using our preferred approximation with a purely signal-dependent rotation ${\mathbf{H}}$, which is unique and can be regarded as the population version of $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. Since this bias is parametrically inestimable, we propose a split-panel jackknife bias correction, and theory shows that it successfully reduces the bias.
The extensive finite-sample experiments suggest that the proposed bias correction works very well, and the empirical application illustrates its usefulness in practice.
\end{abstract}
\textbf{Keywords.} Factor model, Asymptotic bias, Jackknife, Cross-sectional dependence, 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}{\large\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_t$ is written as
\begin{align}\label{Augmodel}
y_{t+h}={\boldsymbol{\gamma}^*}'\mathbf{f}_t^*+\boldsymbol{\beta}' \mathbf{w}_t+\epsilon_{t+h}, \quad t=1,\dots,T,
\end{align}
where $\mathbf{f}_t^*$ is an $r \times 1$ vector of latent predictive factors, $\mathbf{w}_t$ is a $p \times 1$ vector of observable predictors, $(\boldsymbol{\gamma}^{*\prime},\boldsymbol{\beta}')'$ is an $(r+p)\times 1$ vector of their coefficients, and $\epsilon_{t+h}$ is an error term.
The latent $r$ factors drive a large number of $N$ predictors:
\begin{align}\label{factormodel}
x_{t,i}=\mathbf{b}_{i}^{*\prime} \mathbf{f}_t^* + e_{t,i},
\quad t=1,\dots,T,\quad i=1,\dots,N,
\end{align}
where $\mathbf{b}_{i}^{*}$ is an $r\times 1$ vector of unknown factor loadings and $e_{t,i}$ is an error term.
Since $\mathbf{f}_t^*$ is unobserved, it is typically replaced by the principal component (PC) estimator, $\hat{\mathbf{f}}_t$ such that $T^{-1}\sum_{t=1}^T\hat{\mathbf{f}}_t\hat{\mathbf{f}}_t^{\prime}=\mathbf{I}_r$, obtained as $\sqrt{T}$ times the $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})$.
Let $(\hat{\boldsymbol{\gamma}}',\hat{\boldsymbol{\beta}}')'$ be the least squares (LS) estimators of the regression of $y_{t+h}$ on $(\hat{\mathbf{f}}_t^{\prime}, \mathbf{w}_t^{\prime})^{\prime}$.
\cite{StockWatson2002JASA}, \cite{BaiNg2006} and \cite{GoncalvesPerron2014,gonccalves2020bootstrapping} employ the asymptotic approximation of the PC factor by rotated latent factors with a data dependent (but infeasible) rotation matrix $\hat{\mathbf{H}}$:
\begin{align}\label{"consis"}
\hat{\mathbf{f}}_t=\hat{\mathbf{H}}'\mathbf{f}_t^* +o_p(1),
\end{align}
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},\dots,\hat{\lambda}_{r})}$. Using \eqref{"consis"}, we approximate the first term on the right-hand side of \eqref{Augmodel} as $ \boldsymbol{\gamma}^{*\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 what $\hat{\boldsymbol{\gamma}}$ estimates.
\cite{BaiNg2006} show that so 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., no asymptotic bias).
When $N$ is not sufficiently large relative to $T$ to satisfy $\sqrt{T}/N \to c\in(0,\infty)$,
\cite{Ludvigson2011} show that $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ has an asymptotic bias and provide an analytical bias correction. \cite{GoncalvesPerron2014} refine the asymptotic bias expression and propose an analytical bias correction as well as a wild bootstrap method to remove the asymptotic bias. \cite{gonccalves2020bootstrapping} extend \cite{GoncalvesPerron2014} to allow the bias corrections for errors $e_{t,i}$ to be cross-correlated using the large covariance matrix estimator proposed by \cite{BickelLevina2008}.
The existing literature derives the asymptotic bias by assuming that all the $r$ largest eigenvalues of the sample covariance matrix of $x_{t,i}$, $(\hat{\lambda}_{1},\dots,\hat{\lambda}_{r})$, diverge proportionally to $N$ for sufficiently large samples. This is known as the \textit{strong factor} (SF) model. To the best of our knowledge, this paper is the first to study the derivation of asymptotic bias for the more general \textit{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,\dots,r$ for sufficiently large samples.
A growing body of literature suggests that such WFs are widely observed in real data; see \citet{BaileyEtAl2016,BaileyEtAl2021}, \cite{DeMol2008}, \cite{Freyaldenhoven21JoE}, \cite{Onatski2010}, \cite{UY2019,UY2019inference}, \cite{WeiZhang2023}, among many others.
In particular, we show that if $\sqrt{T}/N^{(3\alpha_r - \alpha_1)/2} \to c_1\in(0,\infty)$, $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$ has an asymptotic bias.
Using the derived asymptotic bias expression, we propose an analytical bias-corrected estimator, called the ``plug-in estimator'' in \cite{GoncalvesPerron2014,gonccalves2020bootstrapping}, using the POET estimator proposed by \cite{FanEtAl2013} and extended for WF models by \cite{RunyuEtAl2024}, allowing for cross-sectional and serial correlations in $e_{t,i}$.
As introduced in \cite{BaiNg2023} and \cite{jiang2023revisiting}, there are variants of asymptotically equivalent data-dependent rotations 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 show that if $\sqrt{T}/N^{\alpha_r}\to c_2\in(0,\infty)$, $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q})$ will have an asymptotic bias, where $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}=\hat{\mathbf{H}}_q^{-1}\boldsymbol{\gamma}^{*}$. Note that the rate $\sqrt{T}/N^{\alpha_r}$ is not slower than $\sqrt{T}/N^{(3\alpha_r - \alpha_1)/2}$.
Furthermore, it is found that the bias of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q})$ is generally smaller in magnitude than that of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}})$. We propose an analytical bias correction that makes the bias virtually zero even with small samples. In addition, we show that the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q})$ becomes \textit{exactly} zero when $\mathbf{f}_t^*$ and $\mathbf{w}_t$ are uncorrelated. To exploit this property, we also suggest projecting out $\mathbf{w}_t$ from $x_{t,i}$ and then extracting the factors, say $\hat{\mathbf{f}}_w$ to augment the regression. In fact, this is often done in practice; see \cite{YamamotoHara2022}, for example.
We have discussed the asymptotic biases of $\hat{\boldsymbol{\gamma}}$ relative to $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$, and a natural question that arises here is: ``Relative to which rotation matrix should the asymptotic bias of $\hat{\boldsymbol{\gamma}}$ be evaluated?'' Actually, for inferential purposes in particular, our preferred choice is neither ${\hat{\mathbf{H}}}$ nor ${\hat{\mathbf{H}}_q}$; instead, we argue for choosing $\mathbf{H}$ composed only of the signals $\mathbf{f}_t^*$ and $\mathbf{b}_i^*$ such that $T^{-1}\sum_{t=1}^T{\mathbf{f}}_t^0 {\mathbf{f}}_t^{0\prime}=\mathbf{I}_r$,
where
\begin{align}\label{f0}
\mathbf{f}_t^0:=\mathbf{H}'\mathbf{f}_t^*.
\end{align}
\cite{jiang2023revisiting} have shown that up to sign the PC estimator, $\hat{\mathbf{f}}_t$, is \textit{consistent} to $\mathbf{f}_t^0$ and $\mathbf{H}$ \textit{always} exists and is unique. Therefore, $\mathbf{H}$ can be considered as the population version of the data-dependent rotation matrices, $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. Substituting \eqref{f0} into \eqref{Augmodel}, we can rewrite the model as
\begin{align}\label{consis}
y_{t+h}={\boldsymbol{\gamma}^0}'\mathbf{f}_t^0+\boldsymbol{\beta}' \mathbf{w}_t+\epsilon_{t+h},
\end{align}
where $\boldsymbol{\gamma}^0= {\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*}$. Thus, a regression of $y_{t+h}$ on $(\hat{\mathbf{f}}_t,\mathbf{w}_t)$ consistently estimates the parameters $(\boldsymbol{\gamma}^{0\prime},\boldsymbol{\beta}')'$. We prefer to consider the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^0)$. Note that $\boldsymbol{\gamma}^0$ is a function of the signal parameters $\{\boldsymbol{\gamma}^*,\mathbf{f}_{t}^*,\mathbf{b}_{i}^*\}$, whereas $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ are functions of data, $(x_{t,i})$.
We show that if $\sqrt{T}/N^{(3\alpha_r - \alpha_1)/2} \to c_1\in(0,\infty)$, $\sqrt{T}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^0)$ has an asymptotic bias,
but it seems difficult to obtain the explicit analytical expression.
Therefore, we propose the \textit{split-panel Jackknife bias-correction}, and show that it effectively reduces the bias. Our Jackknife bias-correction is generally less computationally expensive than bootstrapping, while allowing for cross-sectional and serial correlations in $e_{t,i}$.
The finite sample behavior in terms of bias, standard deviation and tail distribution of the estimators $(\hat{\boldsymbol{\gamma}},\hat{\boldsymbol{\beta}})$ and their bias-corrected versions relative to ($\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$, $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$, $\boldsymbol{\gamma}^0$) and $\boldsymbol{\beta}$, are investigated for SF and WF models with cross-sectional and serially correlated errors, $e_{t,i}$. Throughout the design, the analytically bias-corrected estimator with respect to $(\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q},\boldsymbol{\beta})$ has the least bias, the least size distortion of t-tests, with the smallest standard errors, while the performance of the estimators with respect to $(\boldsymbol{\gamma}_{\hat{\mathbf{H}}},\boldsymbol{\beta})$ is far worse than others. Our preferred jackknife bias-corrected estimator with respect to the parameter $(\boldsymbol{\gamma}^0,\boldsymbol{\beta})$ comes in second, closely following the performance of the bias-corrected estimator with respect to $(\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q},\boldsymbol{\beta})$. The proposed jackknife method successfully reduces the bias with a small increase in the standard deviation, leading to the correct size of the test even for very weak factor models for sufficiently large sample sizes.
In addition, the empirical power curve of the significance test shows that, on balance, the split-panel jackknife estimator appears to be the most reliable among the estimators compared.
We apply the bias-corrected estimator to the factor-augmented forecast regression of bond yields $y_{t+h}$ on the factors extracted from 131 monthly macroeconomic series and one observed predictor, the \cite{CochranePiazzesi2005} factor, over the period January 1982 to December 2002. In the application, we have regarded $(\boldsymbol{\gamma}^0,\boldsymbol{\beta})$ as the parameter to be estimated, and the results show that the jackknife appears to effectively correct the bias of the LS estimator, thus providing more reliable inference. In particular, the jackknife estimator has given significantly different results from the rest of the estimators for testing the equal explanatory power of the factors.
The rest of the paper is organized as follows. Section \ref{sec:2} introduces models and assumptions, and also clarifies specific objects of our analysis.
Section \ref{sec:bias} derives the asymptotic biases of the augmented regression estimator with two different data dependent rotations, then derives the asymptotic bias with respect to the signal parameter and proposes a jackknife bias reduction.
To reduce the derived bias, Section \ref{sec:Mw} proposes to orthogonalize the predictor variables to the observable factor in the augmented regression before extracting the factors.
Section \ref{sec: MC} summarizes finite sample experiments and Section \ref{sec:emp} illustrates the practical usefulness of our proposed methods.
Section \ref{sec:con} contains some concluding remarks.
Mathematical proofs and additional experimental results can be found in the Online Appendix.
\noindent \textbf{Notations}:
For any matrix $\mathbf{M}=(m_{t,i})\in\mathbb{R}^{T\times N}$, we define the Frobenius norm
as
$\|\mathbf{M}\|_{\F}=(\sum_{t,i}m_{t,i}^2)^{1/2}$ and $\|\mathbf{M}\|_2$
denotes the square root of the largest eigenvalue for a positive semi-definite matrix $\mathbf{M}'\mathbf{M}$.
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$. We use $\lesssim$ ($\gtrsim$) to represent $\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 do not specifically mention it. $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}{\large\bf}{Preliminaries}\label{sec:2}
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Factor-augmented regression}
The factor-augmented regression model \eqref{Augmodel} can be rewritten in a matrix form as:
\begin{align}\label{augmodel_mat}
\mathbf{y}={\mathbf{F}}^*\boldsymbol{\gamma}^* + \mathbf{W} \boldsymbol{\beta} + \boldsymbol{\epsilon}
=\mathbf{Z}^* \boldsymbol{\delta}^* + \boldsymbol{\epsilon},
\end{align}
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,\dots, \mathbf{w}_T)'$, $\mathbf{Z}^* = (\mathbf{F}^*,\mathbf{W})$, and $\boldsymbol{\delta}^* = (\boldsymbol{\gamma}^{*\prime},\boldsymbol{\beta}')'$.
In line with \eqref{factormodel}, the latent factor model for the $T \times N$ matrix of predictors is given by
\begin{align}\label{factormodel_mat}
\mathbf{X}={\mathbf{F}}^*\mathbf{B}^{* \prime} + \mathbf{E},
\end{align}
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)$
be the $r$ largest eigenvalues of the signal covariance matrix of model \eqref{factormodel_mat}, $T^{-1}\mathbf{F}^*\mathbf{B}^{\ast\prime}\mathbf{B}^{\ast}\mathbf{F}^{\ast\prime}$, and set $\boldsymbol{\Lambda}=\diag(\lambda_1,\dots,\lambda_r)$. We allow them to diverge at different rates; namely, we assume $\lambda_k \asymp N^{\alpha_k}$ with $0<\alpha_k\leq 1$ for $k=1,\dots,r$. We call the model in \eqref{factormodel_mat} with such eigenvalues a weak factor (WF) model and that with $\alpha_r=1$ a strong factor (SF) model as a special case.
\cite{StockWatson2002JASA} proposed to extract PC factors from the predictors $\mathbf{X}$, and then use them in the forecast regression. The PC estimator, $(\hat{\mathbf{F}}, \hat{\mathbf{B}})$, is defined as a solution to minimization of $\left\|\mathbf{X}-\mathbf{F B}^{\prime}\right\|_{\mathrm{F}}^2$ subject to the $r^2$ restrictions, $T^{-1} \mathbf{F}^{\prime} \mathbf{F}=\mathbf{I}_r$ and $\mathbf{B}^{\prime} \mathbf{B}$ diagonal with rank $r$. This reduces to the eigen-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 obtained by $\hat{\mathbf{B}}=T^{-1} \mathbf{X}^{\prime} \hat{\mathbf{F}}$. By the construction, we have $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)$.
Setting $\hat{\mathbf{Z}}=(\hat{\mathbf{F}},\mathbf{W})$, we are interested in the LS estimator,
\begin{align}\label{ols}
\hat{\boldsymbol{\delta}}=(\hat\mathbf{Z}'\hat\mathbf{Z})^{-1}\hat\mathbf{Z}'\mathbf{y}.
\end{align}
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Rotation matrices}
The problem lies in the impossibility of separately identifying $(\mathbf{F}^*,\boldsymbol{\gamma}^*)$ and $(\mathbf{F}^*,\mathbf{B}^*)$ due to the rotation indeterminacy when we consider estimation. More precisely, the products ${\mathbf{F}}^*\boldsymbol{\gamma}^*$ and $\mathbf{F}^* \mathbf{B}^{*\prime}$ are observationally equivalent to $({\mathbf{F}}^* \mathbf{R}) (\mathbf{R}^{-1} \boldsymbol{\gamma}^*)$ and $(\mathbf{F}^* \mathbf{R}) (\mathbf{R}^{-1} \mathbf{B}^{*\prime})$, respectively, for any $r\times r$ invertible matrix $\mathbf{R}$ while the LS and PC estimators, $\hat{\boldsymbol{\delta}}$ and $(\hat{\mathbf{F}},\hat{\mathbf{B}})$, are uniquely determined. Therefore, replacing $\mathbf{F}^*\mathbf{R}$ with $\hat{\mathbf{F}}$ in \eqref{augmodel_mat} raises the question; which rotation matrices $\mathbf{R}$ justify approximation $\hat{\mathbf{F}}\approx \mathbf{F}^*\mathbf{R}$ and $\hat{\boldsymbol{\gamma}}\approx \mathbf{R}^{-1}\boldsymbol{\gamma}^*$? In this paper, we consider three rotation matrices, $\hat{\mathbf{H}}$, $\hat{\mathbf{H}}_q$, and $\mathbf{H}$, for such $\mathbf{R}$, where the first two matrices depend on data, but the last one consists only of $(\mathbf{F}^*,\mathbf{B}^*)$; see \cite{BaiNg2023} for more information on data-dependent rotations and \cite{jiang2023revisiting} for $\mathbf{H}$.
\subsubsection{First data-dependent rotation matrix: $\hat{\mathbf{H}}$}
\cite{BaiNg2002,BaiNg2006} and \cite{StockWatson2002JASA} consider the approximation
\begin{align}\label{FHhat}
\hat{\mathbf{F}}=\mathbf{F}^* \hat{\mathbf{H}} + o_p(1)~~~\text{with}~~~
\hat{\mathbf{H}} = {\mathbf{B}^*}'\mathbf{B}^*(T^{-1}{\mathbf{F}^*}'\hat{\mathbf{F}}) \hat{\boldsymbol{\Lambda}}^{-1},
\end{align}
where and hereafter $o_p(1)$ is understood to apply row-wise.
Using this approximation, we rewrite the model in \eqref{augmodel_mat} as
\begin{align}\label{augmodel_Hhat}
\mathbf{y}={\mathbf{F}^*}\hat{\mathbf{H}}\boldsymbol{\gamma}_{\hat{\mathbf{H}}} + \mathbf{W} \boldsymbol{\beta} + \boldsymbol{\epsilon}
=\hat{\mathbf{Z}} \boldsymbol{\delta}_{\hat{\mathbf{H}}} + \mathbf{u}
\end{align}
with $\mathbf{u} = \boldsymbol{\epsilon} - ({\hat\mathbf{F} -\mathbf{F}^* \hat{\mathbf{H}}})\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^*$, where $\boldsymbol{\delta}_{\hat{\mathbf{H}}}=(\boldsymbol{\gamma}_{\hat{\mathbf{H}}}',\boldsymbol{\beta}')'$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}=\hat{\mathbf{H}}^{-1} \boldsymbol{\gamma}^*$.
Thus, the LS estimator $\hat{\boldsymbol{\delta}}=(\hat{\boldsymbol{\gamma}}',\hat{\boldsymbol{\beta}}')'$, obtained by the regression of $\mathbf{y}$ on $\hat{\mathbf{Z}}$, is interpreted as the ``estimator'' of the data-dependent parameter $\boldsymbol{\delta}_{\hat{\mathbf{H}}}$.
The error term $\mathbf{u}$ contains $({\hat\mathbf{F} -\mathbf{F}^* \hat{\mathbf{H}}})\hat{\mathbf{H}}^{-1}\boldsymbol{\gamma}^*$ due to the approximation by $\hat{\mathbf{F}}$, which is correlated with the regressor $\hat{\mathbf{Z}}$ and can cause the asymptotic bias in $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}})$.
\subsubsection{Second data-dependent rotation matrix: $\hat{\mathbf{H}}_q$}
As pointed out by \cite{BaiNg2023} and \cite{jiang2023revisiting}, there are many asymptotically equivalent rotation matrices to $\hat{\mathbf{H}}$ for the approximation of $\hat{\mathbf{F}}$.
In particular, \cite{BaiNg2023}, \cite{jiang2023revisiting} and \cite{WeiZhang2023} suggested the approximation
\begin{align}\label{FHhatq}
\hat{\mathbf{F}}=\mathbf{F}^* \hat{\mathbf{H}}_q + o_p(1) ~~~\text{with}~~~
\hat{\mathbf{H}}_q = (T^{-1}\hat{\mathbf{F}}'{\mathbf{F}^*})^{-1}.
\end{align}
A similar discussion above will lead to a study of the asymptotic bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$,
where $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}=(\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}',\boldsymbol{\beta}')'$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}=\hat{\mathbf{H}}_q^{-1} \boldsymbol{\gamma}^*$.
\subsubsection{Population rotation matrix: $\mathbf{H}$}
We consider yet another rotation matrix, which is regarded as a population version of any data-dependent rotation matrix asymptotically equivalent to $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$.
\cite{jiang2023revisiting} have shown that there always exists the unique (up to sign) rotation matrix $\mathbf{H}$, which is a pure function of $(\mathbf{F}^* , \mathbf{B}^*)$, such that $T^{-1}\mathbf{F}^{0 \prime}\mathbf{F}^0 = \mathbf{I}_r$ and $\mathbf{B}^{0 \prime}\mathbf{B}^0 = \boldsymbol{\Lambda}$ with $\mathbf{F}^0:=\mathbf{F}^* \mathbf{H}$ and $\mathbf{B}^0:=\mathbf{B}^* \mathbf{H}^{-1\prime}$. Such a unique rotation matrix is indeed represented as
\begin{align}
\mathbf{H} = \mathbf{P}\mathbf{V}^{-1/2},
\end{align}
where $\mathbf{P}$ is the eigenvector matrix of ${\mathbf{B}^*}'\mathbf{B}^*\left(T^{-1}{\mathbf{F}^*}'\mathbf{F}^*\right)$ corresponding to $(\lambda_1,\dots\lambda_r)$, and $\mathbf{V}=\mathbf{P}\left(T^{-1}{\mathbf{F}^*}'\mathbf{F}^*\right)\mathbf{P}'$.
This rotation matrix is straightforwardly rewritten as
\begin{align}\label{HalaHhat}
{\mathbf{H}} = {\mathbf{B}^*}'\mathbf{B}^*(T^{-1}{\mathbf{F}^*}'{\mathbf{F}^0}) {\boldsymbol{\Lambda}}^{-1}
=(T^{-1}{\mathbf{F}^0}'{\mathbf{F}^*})^{-1},
\end{align}
insisting that $\mathbf{H}$ is the population counterpart of $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$. \cite{jiang2023revisiting} show that
\begin{align}
\hat{\mathbf{F}} = \mathbf{F}^0 + o_p(1)
\end{align}
and $\hat{\mathbf{B}} = \mathbf{B}^0 + o_p(1)$ in the WF setting. As in the same way, this approximation rewrites the model in \eqref{augmodel_mat} as
\begin{align}\label{pseudomodel}
\mathbf{y}={\mathbf{F}^0}\boldsymbol{\gamma}^0 + \mathbf{W} \boldsymbol{\beta} + \boldsymbol{\epsilon}
=\hat{\mathbf{Z}} \boldsymbol{\delta}^0 + \mathbf{u}^0
\end{align}
with $\mathbf{u}^0 = \boldsymbol{\epsilon}- ({\hat\mathbf{F} -\mathbf{F}^* {\mathbf{H}}})\boldsymbol{\gamma}^0$ where $\boldsymbol{\delta}^0=(\boldsymbol{\gamma}^{0\prime},\boldsymbol{\beta}')'$ and $\boldsymbol{\gamma}^0={\mathbf{H}}^{-1} \boldsymbol{\gamma}^*$.
We also study the asymptotic bias of
$\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$. Unlike the others, $\boldsymbol{\delta}^0$ does not depend on the data.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Assumptions}
For the asymptotic analysis that follows, we make the following assumptions.
\begin{ass}\normalfont\label{ass:eigen}
The smallest eigenvalues of ${\mathbf{B}^*}'\mathbf{B}^*$ and $T^{-1}{\mathbf{F}^*}'\mathbf{F}^*$ are bounded away from zero.
\end{ass}
\begin{ass}[Idiosyncratic errors]\normalfont\label{ass:errors}${ }^{ }$\\
(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}\})$.
\end{ass}
\begin{ass}[Signal strength]\normalfont\label{ass:signal}
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$.
\end{ass}
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}$}. 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 \citet[Section 5]{BaiNg2023} and/or $T^{-1}\mathbf{F}^{*\prime}\mathbf{F}^*=\mathbf{I}_r$ in \cite{Freyaldenhoven21JoE}.
As discussed earlier, the PC estimators $(\hat{\mathbf{F}},\hat{\mathbf{B}})$ are estimating the pseudo-true parameters $({\mathbf{F}^0},{\mathbf{B}^0})$. We directly impose the following assumptions on them.
\begin{ass}[Factors and Loadings]\normalfont\label{ass:factor and loadings}
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) As $N, T \rightarrow \infty$, $T^{-1} \sum_{t=1}^{T} \left( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i}\right) \left( \mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i}\right)'\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}\left(\mathbf{N}^{-\frac{1}{2}}\sum_{i=1}^{N} \mathbf{b}_i^0e_{t,i} \right)$.
\end{ass}
The moment restrictions in Assumption 4 (iii), (iv)
are essentially similar to Assumptions D, F2
in \cite{Bai2003}, and Assumption 4 (ii)
is similar moment restriction related for $\mathbf{b}_i^0$.
Assumption (v) is similar to Assumption 3(e) in \cite{GoncalvesPerron2014}.
Now we impose assumptions on the pseudo-true augmented model \eqref{pseudomodel}:
\begin{ass}[Weak dependence between idiosyncratic errors and regression errors]\normalfont\label{ass:2errors}
The $r \times r$ matrix satisfies $\E\left\|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}\right\|_2^2 \leq M$.
\end{ass}
\begin{ass}[Moments, parameters and CLT]
\normalfont\label{ass:Aug_errors}${ }^{ }$\\
(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.
\end{ass}
Assumptions \ref{ass:2errors} and \ref{ass:Aug_errors} are similar to Assumption 4 in \cite{GoncalvesPerron2014} and Assumption E in \cite{BaiNg2006}, respectively.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Bias Analysis}\label{sec:bias}
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Bias analysis with $\hat{\mathbf{H}}$}\label{subsec:Hhat}
The asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ has only been investigated in the literature, assuming the SF model.
\cite{BaiNg2006} show that the limiting distribution of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ is centered at zero (i.e., no asymptotic bias) so long as $\sqrt{T}/N \to 0$.
Under a relaxed condition of $\sqrt{T}/N \to c\in[0,\infty)$,
\cite{Ludvigson2011} show that $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ has an asymptotic bias and provide an analytical correction. \cite{GoncalvesPerron2014} refine the bias expression, and propose an analytical correction as well as a wild bootstrap method to remove the bias. \cite{gonccalves2020bootstrapping} extend \cite{GoncalvesPerron2014} to allow errors $e_{t,i}$ to be cross-correlated using the large covariance matrices estimator proposed by \cite{BickelLevina2008}.
We are ready to analyze the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ in the WF setting.
\begin{thm}\label{thm:bias_Hhat}
Suppose Assumptions \ref{ass:eigen}--\ref{ass:Aug_errors} 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
\begin{align*}
\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)
\end{align*}
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}^*
~~~ \text{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$.
\end{thm}
Theorem \ref{thm:bias_Hhat} establishes the asymptotic normality of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$, with the bias appearing in the location. Because $\nu$ takes the values one (when $\alpha_1=\alpha_r$) or zero (when $\alpha_1>\alpha_r$) only, when $c_1=0$ the asymptotic bias becomes zero and the distribution is identical to the results in \citet[Theorem 4]{jiang2023revisiting}.
Essentially, asymptotic bias is due to the correlation between the regressor matrix $\hat{\mathbf{F}}$ and the estimation error shown in \eqref{augmodel_Hhat}. Therefore, when $\mathbf{W}$ and $\mathbf{F}^*$ are orthogonal, $\hat{\boldsymbol{\beta}}$ will not have an asymptotic bias since $\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0}=\mathbf{0}$ and $\boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}$ is block diagonal.
Note that the bias term $\boldsymbol{\kappa}_{\boldsymbol{\delta}^*}$ is estimable, and the associated bias-corrected estimator is given in \eqref{bcHhats} in Section \ref{sec: MC}.
In general, $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$ can be consistently estimated by $\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\delta}}=(T^{-1}\hat{\mathbf{Z}}^{\prime}\hat{\mathbf{Z}})^{-1}
\hat{\boldsymbol{\Sigma}}_{\mathbf{Z}^0 \boldsymbol{\epsilon}}
(T^{-1}\hat{\mathbf{Z}}^{\prime}\hat{\mathbf{Z}})^{-1}
$
so long as $\hat{\boldsymbol{\Sigma}}_{\mathbf{Z}^0 \boldsymbol{\epsilon}}$ is consistent to ${\boldsymbol{\Sigma}}_{\mathbf{Z}^0 \boldsymbol{\epsilon}}$.
A choice of $\hat{\boldsymbol{\Sigma}}_{\mathbf{Z}^0 \boldsymbol{\epsilon}}$ depends on the property of $\epsilon_t$. For example, when $\epsilon_t$ is heteroskedastic, choose $\hat{\boldsymbol{\Sigma}}_{\mathbf{Z}^0 \boldsymbol{\epsilon}}=
T^{-1}\sum_{t=1+h}^{T+h}\hat\mathbf{z}_t \hat{\epsilon}_t^2 \hat\mathbf{z}_t'$, where $\hat{\epsilon}_{t+h}=y_{t+h}-\hat{\boldsymbol{\delta}}'\hat{\mathbf{z}}_t$.
The expressions $c_1 \mathbf{G} = \lim_{N, T \rightarrow \infty} \sqrt{T}\mathbf{N}^{\frac{1}{2}} \boldsymbol{\Gamma} \mathbf{D}^{-2} \mathbf{N}^{-\frac{3}{2}} $ and
$\nu \mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$
suggest a complicated asymptotic bias structure, depending on the structure of $(\alpha_1,\dots,\alpha_r)$. When all the divergence rates are identical, in which $\alpha=\alpha_1=\cdots=\alpha_r$, $c_1 = \lim_{N, T \rightarrow \infty} \sqrt{T}/N^{\alpha}$ and $\nu=1$ so that $c_1\mathbf{G} = c_1\boldsymbol{\Gamma} \mathbf{D}^{-2}$ and $\nu \mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1} = \mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$.
When $\alpha_1>\alpha_r$, $\nu=0$ so that the term $\mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$ disappears, while $c_1 = \lim_{N, T \rightarrow \infty}\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r}$ and the structure of the matrix $\mathbf{G}$ depends on how many factors diverge at the rate of $N^{\alpha_1}$ and $N^{\alpha_r}$. To see this, for simplicity, suppose that all the factors diverge in different rates (i.e., $\alpha_1>\alpha_2>\cdots>\alpha_r$). In $c_1 \mathbf{G}$, $\sqrt{T}\boldsymbol{\Gamma}\mathbf{D}^{-2}$ is pre-multiplied by $\mathbf{N}^{\frac{1}{2}}=\diag(N^{\frac{1}{2}\alpha_1},...,N^{\frac{1}{2}\alpha_r})$ and post-multiplied by $\mathbf{N}^{-\frac{3}{2}}=\diag(N^{-\frac{3}{2}\alpha_1},...,N^{-\frac{3}{2}\alpha_r})$.
The combination of elements that disappear at the slowest rate $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r}$ is given by the product of the first element in $\mathbf{N}^{\frac{1}{2}}$ and the last element in $\mathbf{N}^{-\frac{3}{2}}$. Therefore, the only non-zero element in $c_1\mathbf{G}$ is the $(1,r)$-th element, as other elements tend to zero.
Now let the first $u$ factors and the last $s$ factors ($u+s\leq r)$ share the equal divergence rates (i.e. $\alpha_1=\cdots=\alpha_u>\cdots>\alpha_{r-s+1}=\cdots=\alpha_r$). Then, a similar discussion leads to the conclusion that
the $u\times s$ sub-matrix in the upper right corner of $c_1\mathbf{G}$ is non-zero.
The structure of $c_1\mathbf{G}$ for different sets of $\alpha_k$, $k=1,\dots,r$, is summarized in the following corollary.
\begin{cor}\label{cor1}
Suppose that the conditions for the results in Theorem \ref{thm:bias_Hhat} are satisfied and $c_1 \in (0,\infty)$.
Consider an $r\times 1$ vector $\boldsymbol{\alpha}=(\alpha_1,\alpha_2,\ldots,\alpha_r)'$. Let an $r\times 1$ vector of a binary variable be $\mathbf{e}_{\alpha_1}$, which replaces elements in $\boldsymbol{\alpha}$ with $1$ if they are $\alpha_1$ and $0$ otherwise. Similarly, define an $r\times 1$ vector of a binary variable $\mathbf{e}_{\alpha_r}$ for $\alpha_r$. By construction, $||\mathbf{e}_{\alpha_1}||_2^2+||\mathbf{e}_{\alpha_r}||_2^2\leq r$.
If $\alpha_1=\alpha_r$, then
$c_1\mathbf{G} = c_1 \boldsymbol{\Gamma} \mathbf{D}^{-2}$.
If $\alpha_1>\alpha_r$, then
$c_1\mathbf{G} = c_1 (\mathbf{e}_{\alpha_1} \mathbf{e}_{\alpha_r}')\odot\boldsymbol{\Gamma} \mathbf{D}^{-2}$.
\end{cor}
Theorem \ref{thm:bias_Hhat} and Corollary \ref{cor1} can be viewed as the generalized version of the results in \cite{GoncalvesPerron2014,gonccalves2020bootstrapping} for WF models, which tell us the following. If $\alpha_1 =\alpha_r = 1$, then the model becomes a SF model, and the condition $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r}\to c_1$ reduces to $\sqrt{T}/N\to c_1$ and our result reduces to the result in Theorem 2.1 in \cite{GoncalvesPerron2014}.
The asymptotic bias expression remains the same so long as $\alpha_1=\alpha_r(=\alpha)$, whilst the associated condition $\sqrt{T}/N^{\alpha}\to c_1$ implies that, other things being equal, the weaker the model, the larger the bias.
For $\alpha_1\neq\alpha_r$, the condition $\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r}\to c_1$ implies that, all else being equal, the larger the deviation between $\alpha_1$ and $\alpha_r$ or the weaker the model, the larger the bias. Furthermore, all else being equal, the more heterogeneous the divergence rates are, the more the asymptotic bias deviates from that with $\alpha_1=\alpha_r$.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Bias analysis with $\hat{\mathbf{H}}_q$}
We next show the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$. We reveal that the bias is generally smaller than that of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$.
\begin{thm}\label{thm:bias_H3}
Suppose Assumptions \ref{ass:eigen}--\ref{ass:Aug_errors} 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}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}) \stackrel{d}{\longrightarrow} N\left(c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}, \boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right)
\]
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}} $.
If $\alpha_1=\alpha_r$, then $c_1=c_2$ and $c_2\bar{\mathbf{G}}=c_2\mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$.
If $\alpha_1 > \alpha_r$, then
$c_2\bar{\mathbf{G}} = c_2(\mathbf{e}_{\alpha_r} \mathbf{e}_{\alpha_r}')
\odot\mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$ .
\end{thm}
Theorem \ref{thm:bias_H3} establishes the asymptotic normality of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ with the bias generally smaller than that in Theorem \ref{thm:bias_Hhat}.
To see this, we compare their difference, $c_1{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}-c_2\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$, under the comparable condition of $c_1=c_2$ (and $\nu=1$), or, equivalently, $\alpha_1 = \alpha_r$.
Under the simplest error assumption, $\mathbf{e}_{t}\sim \text{i.i.d.}(\mathbf{0},\sigma_e^2 \mathbf{I}_N)$, we have $\boldsymbol{\Gamma}=\sigma_e^2 \mathbf{D}$, and hence $\mathbf{D}^{-1}\boldsymbol{\Gamma}\mathbf{D}^{-1}=\boldsymbol{\Gamma}\mathbf{D}^{-2}=\mathbf{G} = \bar{\mathbf{G}}$. Therefore, in this case a direct calculation yields
\begin{align}
\|{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}-\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}\|_2
=\left\|\boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}^{-1} \binom{{\mathbf{G}}+\mathbf{D}^{-1}\boldsymbol{\Gamma}\mathbf{D}^{-1}}{\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0} ({\mathbf{G}}-\bar{\mathbf{G}}) } \mathbf{H}_0^{-1} \boldsymbol{\gamma}^*\right\|_2
=\left\|\boldsymbol{\Sigma}_{\mathbf{Z}^0 \mathbf{Z}^0}^{-1} \binom{2\sigma_e^2\mathbf{D}^{-1}}{\mathbf{0} } \mathbf{H}_0^{-1} \boldsymbol{\gamma}^*\right\|_2 \geq 0,
\end{align}
meaning that the magnitude of the bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ is not smaller than that of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$.
We have ${\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}=\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$ only when $\boldsymbol{\gamma}^* = \mathbf{0}$.
It should also be noted that, other things being equal, the bias difference increases as the scaled noise to signal ratio ($\sigma_e^2\mathbf{D}^{-1}$) and/or $||\boldsymbol{\gamma}^*| |_2$ increases.
Regarding the conditions, only difference is that $a_{NT}:=\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r} \to c_1$ in Theorem \ref{thm:bias_Hhat} is replaced by $b_{NT}:=\sqrt{T}N^{-\alpha_r} \to c_2$ in Theorem \ref{thm:bias_H3}. We find that $b_{NT}$ is not slower than $a_{NT}$ by taking the ratio: $a_{NT}/b_{NT}=N^{\frac{1}{2}(\alpha_1-\alpha_r)}\geq 1$, and $c_1=c_2(\neq0)$ only when $\alpha_1=\alpha_r$.
This is essentially because the approximation $T^{-1}\mathbf{W}'(\hat{\mathbf{F}}-\mathbf{F}^* \hat{\mathbf{H}}_{{q}})=o_p(1)$ has a much faster convergence rate than the approximation $T^{-1}\mathbf{W}'(\hat{\mathbf{F}}-\mathbf{F}^* \hat{\mathbf{H}})=o_p(1)$.
The condition $\sqrt{T}/N^{\alpha_r}\to c_2$ in Theorem \ref{thm:bias_H3} implies that, all else being equal, the weaker the model, the larger the bias.
The bias term $\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$ is estimable and the associated bias-corrected estimator is given in \eqref{bcHhats} in Section \ref{sec: MC}.
The expression $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 suggests a complicated asymptotic bias structure, depending on the structure of $(\alpha_1,\dots,\alpha_r)$. When all the divergence rates are identical, in which $\alpha=\alpha_1=\cdots=\alpha_r$, $c_2 = \lim_{N, T \rightarrow \infty} \sqrt{T}/N^{\alpha}$ and $c_2 \bar{\mathbf{G}} = c_2\mathbf{D}^{-1}\boldsymbol{\Gamma} \mathbf{D}^{-1}$.
When $\alpha_1>\alpha_r$, $c_2 = \lim_{N, T \rightarrow \infty}\sqrt{T}/N^{\alpha_r}$ and the structure of the asymptotic bias depends on how many factors diverge at the rate of $N^{\alpha_r}$. In $c_2 \bar{\mathbf{G}}$, $\sqrt{T}\boldsymbol{\Gamma}$ is pre- and post-multiplied by the diagonal matrix $\mathbf{D}^{-1}\mathbf{N}^{-\frac{1}{2}}$ with $\mathbf{N}^{-\frac{1}{2}}=\diag(N^{-\frac{1}{2}\alpha_1},\dots,N^{-\frac{1}{2}\alpha_r})$.
If $s\leq r$ elements in $\boldsymbol{\alpha}$ take the value $\alpha_r$ (i.e. $\alpha_1>\dots>\alpha_{r-s+1}=\dots=\alpha_r$), the elements that disappear at the slowest rate in $\mathbf{N}^{-\frac{1}{2}}$ are the last $s$ diagonal elements. Therefore, the non-zero elements in $c_2\bar{\mathbf{G}}$ are the $s\times s$ block at its bottom-right corner, which leads to the last sentence in Theorem \ref{thm:bias_H3}.
Note also that the bias $\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$ is completely eliminated (i.e. $\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}=\mathbf{0}$) if the observed factor $\mathbf{W}$ is uncorrelated with $\mathbf{F}^*$ (i.e. $\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}^0}=\mathbf{0}$).
To exploit this property, we consider extracting the factor from the predictors after projecting out the observable factor, $\mathbf{w}_t$, in Section \ref{sec:Mw}.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Bias analysis with $\mathbf{H}$}\label{sec:jackknife}
In this section we consider the asymptotic bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$. For any choice of a data dependent invertible rotation matrix, say $\hat{\mathbf{R}}$, we can decompose $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)=\sqrt{T}(\hat{\boldsymbol{\delta}} - \boldsymbol{\delta}_{\hat{\mathbf{R}}})+\sqrt{T}({\boldsymbol{\delta}}_{\hat{\mathbf{R}}} - \boldsymbol{\delta}^0)$. Thus, in principle, for any choice of $\hat{\mathbf{R}}$, the asymptotic bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$ has to be the sum of the biases due to $\sqrt{T}(\hat{\boldsymbol{\delta}} - \boldsymbol{\delta}_{\hat{\mathbf{R}}})$ and $\sqrt{T}({\boldsymbol{\delta}}_{\hat{\mathbf{R}}} - \boldsymbol{\delta}^0)$.
In this section, we employ the decomposition
\begin{align}
\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).
\end{align}
The first term has been investigated in Theorem \ref{thm:bias_H3} and the second term is evaluated as
\begin{align}
\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^0)=\binom{\sqrt{T}(\hat{\mathbf{H}}_q^{-1}-\mathbf{H}^{-1})\boldsymbol{\gamma}^*}{\mathbf{0}}=O_p(\sqrt{T}N^{\frac{1}{2}\alpha_1-\frac{3}{2}\alpha_r}),
\end{align}
but no explicit bias expression is provided. We derive the asymptotic bias by assuming that $\sqrt{T}(\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}-\boldsymbol{\delta}^0)$ tends to a constant vector in probability, 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)$ as $N, T \to \infty$.
For the purpose of deriving the asymptotic bias, this additional assumption is not overly restrictive, since it includes the case in which the bias tends to zero.
\begin{thm}\label{thm:bias_H}
Suppose that Assumptions \ref{ass:eigen}--\ref{ass:Aug_errors} 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)$ as $N, T \to \infty$. Then, further assuming that $\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 whose last $p$ rows are zero, we have
\begin{align*}
& \sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0) \stackrel{d}{\longrightarrow}
N\left(c_1 \mathbf{h}_{\boldsymbol{\gamma}^*} + c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*} ,\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right).
\end{align*}
\end{thm}
The theorem tells that the bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$ has an additional bias due to the difference between $\hat{\mathbf{H}}_q$ and ${\mathbf{H}}$ compared to the bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$. When $\alpha_r=\alpha_1$, we have $c_2=c_1$ and the asymptotic bias is $c_2( \mathbf{h}_{\boldsymbol{\gamma}^*} + \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*})$. If $\alpha_r < \alpha_1$ and $c_1 \in (0,\infty)$, then $c_2=0$, hence the asymptotic bias is $c_1\mathbf{h}_{\boldsymbol{\gamma}^*}$. The sign and the magnitude of $\mathbf{h}_{\boldsymbol{\gamma}^*}$ cannot be identified analytically, but the experimental results below seem to suggest that it has the same sign as the sign of $\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$, rather than they canceling out.
Unlike the asymptotic biases $\boldsymbol{\kappa}_{\boldsymbol{\delta}^*}$ and $\bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}$ in Theorems \ref{thm:bias_Hhat} and \ref{thm:bias_H3},
the asymptotic bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$ is not parametrically estimable as it contains $\mathbf{h}_{\boldsymbol{\gamma}^*}$.
Here, we propose a subsampling method to reduce the bias, called the split-panel jackknife; this has been originally proposed to correct the $O(T^{-1})$ bias of fixed-effects panel data estimators by \cite{DhaeneJochman2015}, and then extended to correct the $O(N^{-1})$ bias as well by \cite{FERNANDEZVAL2016291}. To the best of our knowledge, this is the first application to the estimation of factor-augmented models.
To define the split-panel jackknife,
consider a partition of $\{1, \ldots, N\}$ into two half-panels, $\mathcal{N}_1 :=$ $\{1, \ldots, \lfloor N/ 2 \rfloor\}$ and $\mathcal{N}_2 :=\{\lfloor N/ 2 \rfloor+1, \ldots, N\}$.
Let $\hat{\mathbf{F}}_j$ and $\hat{\boldsymbol{\delta}}_j$ be the PC estimator using the subsample $\mathcal{N}_j$ with $T$ observations and the associated augmented regression estimator for $j=1,2$. Then, the split-panel jackknife bias-corrected estimator is given by\footnote{For finite samples, we introduce randomization of the order of the cross-sectional units to avoid potentially biased information on factors in $\mathcal{N}_j$. See \eqref{bcjkest_R} in Section \ref{sec: MC} for the procedure.}
\begin{align*}\label{bcjkest_1}
\hat{\boldsymbol{\delta}}_{bcjk} = 2 \hat{\boldsymbol{\delta}}-\frac{1}{2}\left(\hat{\boldsymbol{\delta}}_{1}+\hat{\boldsymbol{\delta}}_{2}\right).
\end{align*}
Now we derive the asymptotic bias of the jackknife bias-corrected estimator relative to the parameter $\boldsymbol{\delta}^0$.
\begin{thm}\label{thm:bias_JC}
Suppose that the same assumptions hold as in Theorem \ref{thm:bias_H}. Then, we have
\begin{align*}
& \sqrt{T}(\hat\boldsymbol{\delta}_{bcjk} - \boldsymbol{\delta}^0) \stackrel{d}{\longrightarrow}
N\left((2-2^{\frac{1}{2}(3\alpha_r-\alpha_1)})c_1 \mathbf{h}_{\boldsymbol{\gamma}^*} + (2-2^{\alpha_r})c_2 \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*} ,\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right).
\end{align*}
\end{thm}
The theorem tells that the proposed jackknife bias correction \textit{always} reduces the bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$, but the effectiveness of the bias reduction depends on the weakness of the factor model. For the SF model with $\alpha_r=1$, the asymptotic bias is removed completely. When $\alpha_r=\alpha_1$, the asymptotic bias is $(2-2^{\alpha_r})c_2(\mathbf{h}_{\boldsymbol{\gamma}^*} + \bar{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*})$, which is always smaller than the asymptotic bias of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$, but the correction becomes less effective as $\alpha_r$ deviates from unity. If $\alpha_r < \alpha_1$ and $c_1 \in(0,\infty)$, then the asymptotic bias is $(2-2^{\frac{1}{2}(3\alpha_r-\alpha_1)})c_1 \mathbf{h}_{\boldsymbol{\gamma}^*}$, which is again always smaller than that of $\sqrt{T}(\hat\boldsymbol{\delta} - \boldsymbol{\delta}^0)$, but the correction becomes less effective as $\frac{1}{2}(3\alpha_r - \alpha_1)$ deviates from unity.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Bias Analysis After Model Transformation}\label{sec:Mw}
So far, we have analyzed the asymptotic biases of the augmented regression estimator $\hat{\boldsymbol{\delta}}$ relative to the rotated parameter $\boldsymbol{\delta}^*$ with the different rotation matrices, $\hat{\mathbf{H}}$, $\hat{\mathbf{H}}_q$, and $\mathbf{H}$. All of the results indicate that the bias will be present due to the replacement of rotated $\mathbf{F}^*$ with the PC estimator $\hat{\mathbf{F}}$.
In this section, we introduce an estimation procedure that forces the latent factor and the observed factor to be uncorrelated. Such a procedure will eliminate the bias. In particular, the bias of the estimator relative to the rotated parameter with $\hat{\mathbf{H}}_q$ becomes zero.
To begin the discussion, let $
\mathbf{M}_w=\mathbf{I}_T- \mathbf{P}_w$ with $\mathbf{P}_w = \mathbf{W}\left(\mathbf{W}^{\prime} \mathbf{W}\right)^{-1} \mathbf{W}^{\prime}$. In the first step, we consider the transformed model of $\mathbf{X}$, $\mathbf{X}_w = \mathbf{M}_w \mathbf{X}$, which is written as
\begin{align}
\mathbf{X}_w=\mathbf{F}_w^* \mathbf{B}^{* \prime}+\mathbf{E}_w=\mathbf{F}_w^0 \mathbf{B}_w^{0 \prime}+\mathbf{E}_w,
\end{align}
where $\mathbf{F}_w^* = \mathbf{M}_w \mathbf{F}^*$, $\mathbf{F}_w^0=\mathbf{M}_w\mathbf{F}^* \mathbf{H}_w$, $\mathbf{B}_w^0= \mathbf{B}^* \mathbf{H}_w^{'-1}$, $\mathbf{E}_w=\mathbf{M}_w \mathbf{E}$, and the rotation matrix $\mathbf{H}_w$ is the analogous counterpart of $\mathbf{H}$ in the previous section but for the model of $\mathbf{X}_w$.
Suppose Assumptions \ref{ass:eigen}--\ref{ass:Aug_errors} hold to this model, but defining the variables and parameters with the subscript $w$ in Assumptions \ref{ass:signal}--\ref{ass:Aug_errors}.
$\hat\mathbf{F}_w$ is $\sqrt{T}$ times the eigenvectors associated with the first $r$ largest eigenvalues of $T^{-1}{\mathbf{X}}_w\mathbf{X}_w'$.
$\hat{\mathbf{H}}_{w}$ and $\hat{\mathbf{H}}_{q,w}$ are the data dependent rotation matrices corresponding to $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$ in the previous sections but for the model of $\mathbf{X}_w$.
Now consider the augmented regression model, which can be re-written using $\mathbf{F}_w^*$ as
\begin{align}
\nonumber
\mathbf{y}
=\mathbf{F}^* \boldsymbol{\gamma}^*+\mathbf{W} \boldsymbol{\beta}+\boldsymbol{\epsilon}
=\mathbf{F}_w^* \boldsymbol{\gamma}^*+\mathbf{W} \boldsymbol{\beta}_w +\boldsymbol{\epsilon}
=\mathbf{F}_w^0 \boldsymbol{\gamma}_w^0+\mathbf{W} \boldsymbol{\beta}_w +\boldsymbol{\epsilon},
\end{align}
where $\boldsymbol{\beta}_w = \left(\mathbf{W}^{\prime} \mathbf{W}\right)^{-1} \mathbf{W}^{\prime}\mathbf{F}^*\boldsymbol{\gamma}^*+\boldsymbol{\beta}$ and $\boldsymbol{\gamma}_w^0=\mathbf{H}_w^{-1}\boldsymbol{\gamma}^*$.
The relevant feasible augmented regression estimator is
\begin{align}
\hat{\boldsymbol{\delta}}_w = (\hat{\mathbf{Z}}_w'\hat{\mathbf{Z}}_w)^{-1}\hat{\mathbf{Z}}_w'\mathbf{y},
\end{align}
where $\hat{\mathbf{Z}}_w=(\hat{\mathbf{F}}_w,\mathbf{W})$. As easily seen, the approximations $\hat{\mathbf{F}}_w = \mathbf{F}_w \hat{\mathbf{H}}_w + o_p(1)$,
$\hat{\mathbf{F}}_w = \mathbf{F}_w \hat{\mathbf{H}}_{q,w} + o_p(1)$ and
$\hat{\mathbf{F}}_w = \mathbf{F}_w^0 + o_p(1)$ lead us to study the asymptotic biases of
$\sqrt{T}(\hat{\boldsymbol{\delta}}_w - {\boldsymbol{\delta}}_{\hat{\mathbf{H}}_w})$,
$\sqrt{T}(\hat{\boldsymbol{\delta}}_w - {\boldsymbol{\delta}}_{\hat{\mathbf{H}}_{q,w}})$ and
$\sqrt{T}(\hat{\boldsymbol{\delta}}_w - {\boldsymbol{\delta}}_w^0)$,
where
${\boldsymbol{\delta}}_{\hat{\mathbf{H}}_w}=(\boldsymbol{\gamma}_{\hat{\mathbf{H}}_w}',\boldsymbol{\beta}_w')'$ with $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_w}=\hat{\mathbf{H}}_w^{-1}\boldsymbol{\gamma}^{*}$,
${\boldsymbol{\delta}}_{\hat{\mathbf{H}}_{q,w}}=(\boldsymbol{\gamma}_{\hat{\mathbf{H}}_{q,w}}',\boldsymbol{\beta}_w')'$ with $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_{q,w}}=\hat{\mathbf{H}}_{q,w}^{-1}\boldsymbol{\gamma}^{*}$, and
${\boldsymbol{\delta}}_w^0=(\boldsymbol{\gamma}_w^{0\prime},\boldsymbol{\beta}_w')'$.
It is straightforward to prove the ``$\mathbf{X}_w$'' versions of Theorems \ref{thm:bias_Hhat}-\ref{thm:bias_JC}, putting $\operatorname*{plim}_{N,T\to\infty} T^{-1}\mathbf{W}'\mathbf{F}_w^0 =\boldsymbol{\Sigma}_{\mathbf{W} \mathbf{F}_w^0}= \mathbf{0}$.
There are a few comments to make. After the transformation, the augmented regression coefficient on $\mathbf{Z}_w^* = (\mathbf{F}_w^* , \mathbf{W})$ changes to $\boldsymbol{\delta}_w^* = (\boldsymbol{\gamma}^{*\prime},\boldsymbol{\beta}_w')'$, which is different from $\boldsymbol{\delta}^*$
unless $\mathbf{W}'\mathbf{F}^* = \mathbf{0}$.
Also, because of the orthogonality $\mathbf{W}' \mathbf{F}_w^0 = \mathbf{0}$, the bias of $\hat{\boldsymbol{\beta}}_w-{\boldsymbol{\beta}}_w$ is always zero, and the bias corrections are only applied to $\hat{\boldsymbol{\gamma}}_w$ with respect to the relevant ``parameters''.
We conclude this section by providing the result for the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}_w - {\boldsymbol{\delta}}_{\hat{\mathbf{H}}_{q,w}})$, which is completely eliminated due to the transformation.
\begin{thm}\label{thm:bias_Hw3}
Suppose Assumptions \ref{ass:eigen}--\ref{ass:Aug_errors} for the versions of $\mathbf{X}_w$ hold. If $ \alpha_r>\frac{1}{2}$, $\frac{N^{1-\alpha_r}}{\sqrt{T}} \to 0$, $\sqrt{T}N^{-\alpha_r} \to c_2 \in (0,\infty)$, as $N, T \to \infty$, we have
\begin{align*}
& \sqrt{T}(\hat\boldsymbol{\delta}_w - \boldsymbol{\delta}_{\hat{\mathbf{H}}_{q,w}}) \stackrel{d}{\longrightarrow}
N\left(\mathbf{0},\boldsymbol{\Sigma}_{\boldsymbol{\delta}_w}\right).
\end{align*}
\end{thm}
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Monte Carlo Experiments}\label{sec: MC}
In this section, we examine the finite sample performance of the estimators of
the factor-augmented regressions. In particular, we focus on the bias, standard deviation and size of the t-test.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Design}\label{sec:design}
$\mathbf{X=F}^{0}\mathbf{B}^{0\prime}+\mathbf{E}$, $\mathbf{X=(}
x_{t,i}\mathbf{)}$, $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}$ is generated 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}})$. Form a $T\times
N$ matrix $\mathbf{A}$ whose elements are independent draws from $N(0,1)$ for
each replication. Obtain the singular value decomposition $\mathbf{A=USV}
^{\prime}$, and set $\mathbf{F}^{0}$ the first $r$ columns of $\mathbf{U}$
multiplied by $\sqrt{T}$ and $\mathbf{B}^{0}$ the first $r$ columns of
$\mathbf{V}$ post-multiplied by $\mathbf{D}^{1/2}\mathbf{N}^{1/2}$. Then, given an invertible square matrix of order $r$, $\mathbf{H}$, set
$\mathbf{F}^{\ast}=\mathbf{F}^{0}\mathbf{H}^{-1}$ and $\mathbf{B}^{\ast
}=\mathbf{B}^{0}\mathbf{H}^{\prime}$.
We set $\mathbf{H}=\left(
\begin{smallmatrix}
1 & 1/2\\
1/2 & 2
\end{smallmatrix}\right)$.
The $t^{th}$ rows of $\mathbf{E}$ are generated as $\mathbf{e}_{t}=\rho_{e}\mathbf{e}_{t-1}+(1-\rho_{e}^{2})^{1/2}\boldsymbol{\Sigma}_{e}^{1/2}\boldsymbol{\xi}_{t}$ for $t=2,\dots,T$ with $\mathbf{e}_{1}\sim N(\mathbf{0},\mathbf{I}_N)$, where
$\boldsymbol{\xi}_{t}\sim i.i.d.N(\mathbf{0},\mathbf{I}_{N})$.
It is set to $\rho_{e}=0.2$.
We consider $\boldsymbol{\Sigma}
_{e}=\sigma_{e}^{2}\mathbf{R}_{s}$,
where $\mathbf{R}_{s}$ is the correlation matrix of $\left( \mathbf{I}
_{N}-\theta\mathbf{S}_{s}\right) \left( \mathbf{I}_{N}-\theta\mathbf{S}
_{s}\right) ^{\prime}$, $\mathbf{S}_{s}$ is the row normalized $s^{th}$-order
rook contiguity spatial matrix. We set $s=2$, $\sigma_{e}=0.5$ and $\theta=0.5$.
The factor-augmented regression is generated 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 vector of $\mathbf{F}
^{0}$ and $\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}+\sqrt{1-\rho_{fw}^{2}}\zeta_{t,\ell}),
\]
where $\zeta_{t,\ell}\sim i.i.d.N(0,1\mathbf{)}$, $\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 discussed in 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 to $\sigma_{w}^{2}=1$ and $\sigma_{\epsilon}^{2}=0.5$. We choose $r=2$ and $p=2$,
$(\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}$ for the model with $\alpha_{1}=\alpha_{2}$ are
necessary to ensure the identification of the two largest eigenvalues of
$\mathbb{E(}\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}\mathbb{)}$, $\lambda_{1}$
and $\lambda_{2}$. The experiments are conducted for
$(T,N)=(50,50),(100,100),(200,200)$ with 1,000 replications.
Supposing $\{\mathbf{x}_{t},\mathbf{w}_{t},\mathbf{y}_{t}\}$ are observable in
practice, $\mathbf{F}^{0}$ is estimated by PC using $\mathbf{X}$, which is the
first $r$ largest eigenvectors of $\mathbf{XX}^{\prime}/T$ multiplied by
$\sqrt{T}$. The PC estimator is denoted by $\mathbf{\hat{F}}$. The factor
augmented model is estimated by regressing $y_{t+1}$ on $\mathbf{\hat{z}}
_{t}=(\mathbf{\hat{f}}_{t}^{\prime},\mathbf{w}_{t}^{\prime})^{\prime}$, which
gives the estimates $\boldsymbol{\hat{\delta}}=(\boldsymbol{\hat{\gamma}
}^{\prime},\boldsymbol{\hat{\beta}}^{\prime})^{\prime}$.
Using the different rotation matrices, we compute the average (i.e. bias), the standard deviation, and two-sided t-test at the 5\% level for $\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 bias-corrected estimators. Let the bias estimates of $\hat{\boldsymbol{\delta}}-{\boldsymbol{\delta}}_{\hat{\mathbf{H}}}$ and $\hat{\boldsymbol{\delta}}-{\boldsymbol{\delta}}_{\hat{\mathbf{H}}_q}$ be
\begin{align*}
\hat{\boldsymbol{\kappa}}_{\boldsymbol{\delta}^*}
=-(T^{-1}\mathbf{\hat{Z}}^{\prime}
\mathbf{\hat{Z}})^{-1}
\begin{pmatrix}
\hat{\mathbf{G}}+ \hat{{\bar{\mathbf{G}}}}\\
T^{-1}\mathbf{W}^{\prime}\mathbf{\hat{F}\hat{G}}
\end{pmatrix}
\boldsymbol{\hat{\gamma}},\quad
\hat{\bar{\boldsymbol{\kappa}}}_{\boldsymbol{\delta}^*}
=(T^{-1}\mathbf{\hat{Z}}^{\prime
}\mathbf{\hat{Z}})^{-1}\left(
\begin{array}
[c]{c}
\mathbf{0}\\
T^{-1}\mathbf{W}^{\prime}\mathbf{\hat{F}}\hat{{\bar{\mathbf{G}}}}
\end{array}
\right)\boldsymbol{\hat{\gamma}}
\end{align*}
with
\begin{align*}
\mathbf{\hat{G}=\mathbf{\hat{B}}^{\prime}}\boldsymbol{\hat{\Sigma}}_{e}\mathbf{\mathbf{\hat{B}}(\hat{B}}^{\prime}\mathbf{\hat{B})}^{-2},\quad
\hat{{\bar{\mathbf{G}}}}=(\mathbf{\mathbf{\hat{B}}^{\prime}\mathbf{\hat{B}}
})^{-1}\mathbf{\mathbf{\hat{B}}^{\prime}}\boldsymbol{\hat{\Sigma}}
_{e}\mathbf{\mathbf{\hat{B}}}(\mathbf{\hat{B}}^{\prime}\mathbf{\hat{B}})^{-1},
\end{align*}
where $\boldsymbol{\hat{\Sigma}}_{e}$ is the POET estimator of \cite{RunyuEtAl2024}, which extends \cite{FanEtAl2013} for the estimation in WF models. The associated bias corrected estimators are defined by
\begin{align}\label{bcHhats}
\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}}=\hat{\boldsymbol{\delta}}-\hat{{\boldsymbol{\kappa}}}_{\boldsymbol{\delta}^*}, \quad
\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_{q}}=\boldsymbol{\hat{\delta}}-\hat{\bar{\boldsymbol{\kappa}}}_{\boldsymbol{\delta}^*},
\end{align}
and the bias, standard deviation, and size of the t-test at the 5\% level of $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}} - \boldsymbol{\delta}_{\hat{\mathbf{H}}}$ and $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}_q} - \boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$ are reported.
The panel-split jackknife bias corrected estimator is computed as
\begin{align}\label{bcjkest_R}
\boldsymbol{\hat{\delta}}_{bcjk}=2\boldsymbol{\hat{\delta}}-\boldsymbol{\hat
{\delta}}_{jk}\ \text{with }\boldsymbol{\hat{\delta}}_{jk}=R^{-1}
\sum\nolimits_{s=1}^{R}(\boldsymbol{\hat{\delta}}_{\mathcal{N}_{1}^{(s)}
}+\boldsymbol{\hat{\delta}}_{\mathcal{N}_{2}^{(s)}})/2
\end{align}
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$, with $\mathcal{N}_{1}^{(s)}$ is the
first half and $\mathcal{N}_{2}^{(s)}$ is the second half of $\mathbf{X}
^{(s)}$, whose $N$ columns are randomly re-ordered over $s=1,\dots,R$.
This randomization is to avoid potentially biased information on factors in $\mathcal{N}_j$.
When $N$ is odd, $\mathcal{N}_{1}^{(s)}$ and $\mathcal{N}_{2}^{(s)}$ are chosen to
contain one common index. The order and the sign of the columns of $\mathbf{\hat
{F}}_{\mathcal{N}_{j}^{(s)}}$ are determined 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$. We have chosen
$R=100$. The bias, standard deviation and size of the t-test at the 5\% level of $\hat{\boldsymbol{\delta}}_{bcjk} - \boldsymbol{\delta}^0$ are reported.
Using the same sample, we investigate similar statistics for the model with
the extracted factors from $\mathbf{M}_{w}\mathbf{X}$. Specifically, we define
$\mathbf{F}_{w}^{\ast}=\mathbf{M}_{w}\mathbf{F}^{\ast}$\ where $\mathbf{M}
_{w}=\mathbf{I}_{T}-\mathbf{P}_{w}$ with $\mathbf{P}_{w}=\mathbf{W(W}^{\prime
}\mathbf{W)}^{-1}\mathbf{W}^{\prime}$, so that $\mathbf{F}_{w}^{0}
=\mathbf{F}_{w}^{\ast}\mathbf{H}_{w}$ and $\mathbf{B}_{w}^{0}=\mathbf{B}
^{\ast}\mathbf{H}_{w}^{\prime-1}$ with $\mathbf{H}_{w}=\mathbf{L}
_{w}\mathbf{V}_{w}^{-1/2}\boldsymbol{\Pi}$, where $\mathbf{V}_{w}
=\mathbf{L}_{w}^{\prime}(T^{-1}\mathbf{F}_{w}^{\ast\prime}\mathbf{F}_{w}
^{\ast})\mathbf{L}_{w}$, $\mathbf{L}_{w}$ is the $r\times r$ eigenvector
matrix of $(\mathbf{\mathbf{B}^{\ast\prime}\mathbf{B}}^{\ast})(T^{-1}
\mathbf{F}_{w}^{\ast\prime}\mathbf{F}_{w}^{\ast})$, $\mathbf{\Pi}$ is a
diagonal matrix with elements either $-1$ or $1$, which makes the diagonal of
$\mathbf{H}_{w}$ positive. The PC estimator $\mathbf{\hat{F}}_{w}$ is
the $\sqrt{T}$ times $r$ eigenvectors corresponding to the $r$ largest
eigenvalues of $\mathbf{M}_{w}\mathbf{XX}^{\prime}\mathbf{M}_{w}/T$ and
$\mathbf{\hat{B}}_{w}=\mathbf{X}^{\prime}\mathbf{\hat{F}}_{w}/T$. The
augmented model of the $T\times1$ vector is re-written as
$\mathbf{y} =\mathbf{F}^{\ast}\boldsymbol{\gamma}^{\ast}+\mathbf{W}
\boldsymbol{\beta}+\boldsymbol{\epsilon} =\mathbf{F}_{w}^{*}\boldsymbol{\gamma}^{*}+\mathbf{W}
\boldsymbol{\beta}_{w}+\boldsymbol{\epsilon}$,
where $\boldsymbol{\beta}_{w}=\mathbf{(W}^{\prime}\mathbf{W)}^{-1}
\mathbf{W}^{\prime}\mathbf{F}^{\ast}\boldsymbol{\gamma}^{\ast}
+\boldsymbol{\beta}$.
Regression of $y$ on $\mathbf{z}
_{w,t}=(\mathbf{\hat{f}}_{w}^{\prime},\mathbf{w}_{t}^{\prime})^{\prime}$ gives
$\boldsymbol{\hat{\delta}}_{w}=(\boldsymbol{\hat{\gamma}}_{w}^{\prime
},\boldsymbol{\hat{\beta}}_{w}^{\prime})^{\prime}$, which is the estimator of
$\boldsymbol{\delta}_{w}^{0}=(\boldsymbol{\gamma}_w^{0\prime}
,\boldsymbol{\beta}_{w}^{\prime})^{\prime}$. We investigate analogous
counterpart statistics, $\boldsymbol{\delta}_{\hat{\mathbf{H}}_w}$, $\boldsymbol{\delta
}_{\hat{\mathbf{H}}_{q,w}}$, $\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_w}$,
$\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_{q,w}}$, $\boldsymbol{\hat{\delta}
}_{bcjk,w}$ with respect to the relevant ``parameters''.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Results}
\subsubsection{``Parameters'' }
Before discussing the performance of the estimators, we would like to draw our attention to the ``parameters" we are estimating. In the literature, including \cite{BaiNg2006} and \citet{GoncalvesPerron2014,gonccalves2020bootstrapping}, the LS estimator $\hat{\boldsymbol{\delta}}$ is considered to estimate $\boldsymbol{\delta}_{\hat{\mathbf{H}}}=(\boldsymbol{\gamma}^{* \prime}\hat{\mathbf{H}}^{-1} , \boldsymbol{\beta}')'$, which is noise-dependent via $\hat{\mathbf{H}}$. On the other hand, $\boldsymbol{\delta}^0$ is a pure function of the signals. Normalizing $\boldsymbol{\delta}^0 :=\mathbf{1}$, mean and standard deviation of the second elements of $\boldsymbol{\gamma}^0 :=\boldsymbol{\gamma}^{* \prime}{\mathbf{H}}^{-1}$, $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q} :=\boldsymbol{\gamma}^{* \prime}\hat{\mathbf{H}}_q^{-1}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}} :=\boldsymbol{\gamma}^{* \prime}\hat{\mathbf{H}}^{-1}$ over the replications are plotted over $N=T=50,100,200$ and $(\alpha_1,\alpha_2)=(1.0,1.0),(1.0,0.8),(0.8,0.6)$ in Figures \ref{fig:mean_gamma_hats} and \ref{fig:sd_gamma_hats}, respectively.
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{$\alpha_2 = 1.0$}
\label{fig:bias_f2_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_bias_rhowf0_a0.8_f2.pdf}
\caption{$\alpha_2 = 0.8$}
\label{fig:bias_f2_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_bias_rhowf0_a0.6_f2.pdf}
\caption{$\alpha_2 = 0.6$}
\label{fig:bias_f2_06_06}
\end{subfigure}
\includegraphics[width=0.33\textwidth]{newimages/legend_Parameters_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{Mean of ``parameters" $\gamma_2^0$,
${\gamma}_{\hat{\mathbf{H}}_{q},2}$ and
${\gamma}_{\hat{\mathbf{H}},2}$ over the replications}
\label{fig:mean_gamma_hats}
\end{figure}
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_sd_rhowf0.6_a1_f2.pdf}
\caption{$\alpha_2 = 1.0$}
\label{fig:bias_f2_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_sd_rhowf0.6_a0.8_f2.pdf}
\caption{$\alpha_2 = 0.8$}
\label{fig:bias_f2_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/Parameters_SigE3_rhoe2_sd_rhowf0.6_a0.6_f2.pdf}
\caption{$\alpha_2 = 0.6$}
\label{fig:bias_f2_06_06}
\end{subfigure}
\includegraphics[width=0.33\textwidth]{newimages/legend_Parameters_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{Standard deviation of ``parameters" $\gamma_2^0$,
${\gamma}_{\hat{\mathbf{H}}_{q},2}$ and
${\gamma}_{\hat{\mathbf{H}},2}$ over the replications}
\label{fig:sd_gamma_hats}
\end{figure}
As can be seen, the smaller the sample size and the weaker the model, the more $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ deviates upwards and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ deviates downwards from $\boldsymbol{\gamma}^0$ and their variations tend to increase.
It is a visual confirmation that the values of the ``parameters'' $\boldsymbol{\gamma}^0$, $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ can be very different even for relatively large sample sizes despite their asymptotic equivalence.
Furthermore, even when all the elements in $\boldsymbol{\gamma}^0$ are equal, it is unlikely that a similar property holds for $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$. For example, for $N=T=50$ and the model with $(\alpha_1,\alpha_2)=(1,1)$, in the above experiment, the average values of the $r\times1$ vectors $\boldsymbol{\gamma}^0$, $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ are $(1.00,1.00)'$, $(1.27,1.07)'$ and $(0.94,0.99)'$, respectively. Although all considered rotation matrices share the common probability limit, it does not seem plausible to treat $\boldsymbol{\gamma}^0$, $\boldsymbol{\gamma}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\gamma}_{\hat{\mathbf{H}}_q}$ indifferently in practice.
Therefore, researchers may want to be clear about which ``parameter'' is estimated by $\hat{\boldsymbol{\gamma}}$ in their analysis. Our preferred parameter is $\boldsymbol{\gamma}^0$, because it is the coefficient on the latent factor $\mathbf{f}_t^0$ that the PC estimator $\hat{\mathbf{f}}_t$ consistently estimates.
In what follows we will discuss the bias of the LS estimator $\hat{\boldsymbol{\delta}}$ and its bias-corrected versions relative to $\boldsymbol{\delta}^0$, $\boldsymbol{\delta}_{\hat{\mathbf{H}}}$ and $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$.
\subsubsection{Coefficient on the second factor, $\gamma_2$}
In this section we examine the finite sample performance of the estimated coefficient on the second factor $\hat{f}_{2,t}$, $\hat{\gamma}_2$ and its bias-corrected versions, relative to $\gamma_2^0$, $\gamma_{\hat{\mathbf{H}}_q,2}$ and $\gamma_{\hat{\mathbf{H}},2}$. Specifically, we report the bias, standard deviation (s.d.), and a t-test at the 5\% level of $\hat{\gamma}_2 - \gamma_2^0$, $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}}_q,2}$, $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}},2}$ and their bias-corrected versions, $\hat{\gamma}_{bcjk,2} - \gamma_2^0$, $\hat{\gamma}_{bc\hat{\mathbf{H}}_q,2} - \gamma_{\hat{\mathbf{H}}_q,2}$, $\hat{\gamma}_{bc\hat{\mathbf{H}},2} - \gamma_{\hat{\mathbf{H}},2}$, respectively. The bias-corrected versions are shown as dashed lines in the figures.
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:bias_f2_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:bias_f2_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:bias_f2_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a1_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:bias_f2_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:bias_f2_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:bias_f2_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.50\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{Bias of $\hat{\gamma}_2$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:bias.f2}
\end{figure}
The average bias over the replications is reported in Figure \ref{fig:bias.f2}. The bias $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}},2}$ is always negative and the largest of all in magnitude, and the bias becomes more significant for the smaller sample sizes and weaker factors. With the bias-corrected version, $\hat{\gamma}_{bc\hat{\mathbf{H}},2} - \gamma_{\hat{\mathbf{H}},2}$, always has a smaller bias than without bias-correction, but still other estimators have smaller biases in magnitude.
In contrast, $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}}_q,2}$ has virtually no bias when $w_t$ and $\mathbf{f}_t^*$ are uncorrelated (i.e. $\rho_{wf} = 0$). It has very small bias when $\rho_{wf} = 0.6$, but the bias-corrected version successfully reduces the bias; see $\hat{\gamma}_{bc\hat{\mathbf{H}}_q,2} - \gamma_{\hat{\mathbf{H}}_q,2}$ in the figure.
The most reasonable parameter that is considered to be estimated by $\hat{\gamma}_2$ is the parameter of the pure signals, $\gamma_2^0$, as we have argued. The bias $\hat{\gamma}_2 - \gamma_2^0$ is moderately negatively biased which gets worse for smaller sample sizes and weaker factors. In contrast, the proposed jackknife bias-correction reduces the bias very successfully; see $\hat{\gamma}_{bcjk,2} - \gamma_2^0$ in the figure.
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a1_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:sd_f2_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:sd_f2_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:sd_f2_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a1_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:sd_f2_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:sd_f2_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:sd_f2_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.50\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{Standard deviation of $\hat{\gamma}_2$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:sd.f2}
\end{figure}
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a1_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:test_f2_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:test_f2_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:test_f2_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a1_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:test_f2_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a0.8_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:test_f2_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a0.6_f2.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:test_f2_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.50\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_f2.pdf}
\caption{Size of the t-tests at the 5\% level using $\hat{\gamma}_2$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:test.f2}
\end{figure}
According to Figure \ref{fig:sd.f2}, the standard deviations of $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}}_q,2}$ and its bias-corrected one are the smallest, closely followed by that of $\hat{\gamma}_2 - \gamma_2^0$. For the strong factor model, the jackknife bias-correction does not increase the variation much, but it moderately does for the very weak factor model with the small sample size, which quickly goes down as the sample size increases. In most cases, the variation of $\hat{\gamma}_{bc\hat{\mathbf{H}},2} - \gamma_{\hat{\mathbf{H}},2}$ is the largest, followed by $\hat{\gamma}_{2} - \gamma_{\hat{\mathbf{H}},2}$.
Reflecting the bias and the variation inflation, the size of the t-test, which is reported in Figure \ref{fig:test.f2}, is affected. The size of the test based on $\hat{\gamma}_{bcjk,2} - \gamma_2^0$ and its bias-corrected versions are always around the nominal level. The size of the test based on $\hat{\gamma}_2 - \gamma_2^0$ is correct unless the model is very weak. The test based on the jackknife corrected estimator is correct except for the very weak factor and the very small sample size.
The tests based on $\hat{\gamma}_2 - \gamma_{\hat{\mathbf{H}},2}$ and its bias-corrected are the most unreliable, suffering from enormous size distortion.
\subsubsection{Coefficient on the observed factor, $\beta$}
Let us turn our attention to the performance of the estimated coefficient on the observed factor $w_{t}$, $\beta$ and its bias-corrected versions. Note that the elements of the ``parameters'' $\boldsymbol{\delta}^0$, $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$ and $\boldsymbol{\delta}_{\hat{\mathbf{H}}}$ corresponding to $w_t$ are all $\beta$, and their uncorrected estimators are identical. However, the corresponding elements of the bias-corrected estimators may have different values if the latent factors and the observable factors are correlated. Specifically, we report the bias, standard deviation (s.d.), and a t-test at the 5\% level of $\hat{\beta} - \beta$ and the bias-corrected estimators relative to $\beta$, $\hat{\beta}_{bcjk} - \beta$, $\hat{\beta}_{bc\hat{\mathbf{H}}_q} - \beta$, $\hat{\beta}_{bc\hat{\mathbf{H}}} - \beta$. The bias-corrected versions are shown as dashed lines in the figures.
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a1_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:bias_w_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a0.8_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:bias_w_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0_a0.6_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:bias_w_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a1_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:bias_w_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a0.8_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:bias_w_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_bias_rhowf0.6_a0.6_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:bias_w_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.45\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_w.pdf}
\caption{Bias of $\hat{\beta}$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:bias.w}
\end{figure}
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a1_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:sd_w_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a0.8_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:sd_w_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0_a0.6_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:sd_w_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a1_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:sd_w_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a0.8_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:sd_w_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_sd_rhowf0.6_a0.6_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:sd_w_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.45\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_w.pdf}
\caption{Standard deviation of $\hat{\beta}$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:sd.w}
\end{figure}
Figure \ref{fig:bias.w} shows the average bias over the replications. The bias of all the estimators is zero when $\mathbf{f}_t^*$ and $w_t$ are uncorrelated (i.e. $\rho_{wf}=0$) as expected, because the bias is caused by the estimation effect of $\hat{\mathbf{f}}_t$ and it is not transmitted to the estimator of $\beta$.
The picture changes dramatically when $\mathbf{f}_t^*$ and $w_t$ are correlated (i.e. $\rho_{wf}=0.6$). The LS estimator $\hat{\beta}$ is biased, and the magnitude of the bias tends to increase for weaker models and for smaller sample sizes. The bias-corrections $\hat{\beta}_{bc\hat{\mathbf{H}}_q}$ and $\hat{\beta}_{bc\hat{\mathbf{H}}}$ successfully reduce the bias, always by almost the same amount. The clear winner in terms of bias reduction is the jackknife estimator, $\hat{\beta}_{bcjk}$. The standard deviations of all the estimators shown in Figure \ref{fig:sd.w} are almost identical when $\mathbf{f}_t^*$ and $w_t$ are uncorrelated (i.e. $\rho_{wf}=0$), while the standard deviations of all the bias-corrected estimators are similar but very slightly larger than the non-corrected estimator when $\rho_{wf}=0.6$.
\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a1_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 1.0$}
\label{fig:test_w_00_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a0.8_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.8$}
\label{fig:test_w_00_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0_a0.6_w.pdf}
\caption{$\rho_{wf}=0.0,\alpha_2 = 0.6$}
\label{fig:test_w_00_06}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a1_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 1.0$}
\label{fig:test_w_06_10}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a0.8_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.8$}
\label{fig:test_w_06_08}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/SigE3_rhoe2_test_rhowf0.6_a0.6_w.pdf}
\caption{$\rho_{wf}=0.6,\alpha_2 = 0.6$}
\label{fig:test_w_06_06}
\end{subfigure}
\centering
\includegraphics[width=0.50\textwidth]{newimages/legend_SigE3_rhoe2_bias_rhowf0_a1_w.pdf}
\caption{Size of the t-tests using $\hat{\beta}$ and its bias corrected versions for cross and serially correlated $e_{t,i}$}
\label{fig:test.w}
\end{figure}
The size of the tests is summarized in Figure \ref{fig:test.w}. The size of the tests based on all the estimators is correct when $\rho_{wf}=0.0$, while the size of the tests based on the uncorrected estimators tends to deviate from the nominal level, which is successfully corrected by all the bias-correction methods.
Finally, the experimental results for the models augmented with factors extracted from the prediction variables orthogonalized to the observed factor, $\mathbf{M}_w \mathbf{X}$ are presented in Figures \ref{fig:bias.f2w}-\ref{fig:test.ww} in the online appendix. These results confirm that the $\hat{\boldsymbol{\delta}}_w - \boldsymbol{\delta}_{\hat{\mathbf{H}}_{q,w}}$ has little bias for all the designs, including the case with $\rho_{wf}=0.6$, for the small sample size and for the weakest factor model, as predicted by our theory. However, if researchers are interested in estimating the parameter $\boldsymbol{\delta}_w^0$ rather than the random vector $\boldsymbol{\delta}_{\hat{\mathbf{H}}_{q,w}}$, the jackknife estimator $\hat\boldsymbol{\delta}_{bcjk,w}$ may be preferred.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Size-adjusted power curve of the test of significance}
To see how the asymptotic bias can affect the power curve of the t-test of significance, $H_0 : \delta_{\mathbf{R},k}=0$ against $H_1 : \delta_{\mathbf{R},k}\neq 0$ for $\mathbf{R} = \mathbf{H}, \hat{\mathbf{H}}_q, \hat{\mathbf{H}}$, $k=1,\dots,r+p$, we conduct the following experiments. The data generating process (DGP) is identical to that described in Section \ref{sec:design}, except that we change the value of an element in $\boldsymbol{\delta}^0$ under the test between $-0.4$ and $0.4$ by $0.025$, keeping other elements at unity. We have chosen the case for $(\alpha_1, \alpha_2)=(0.8,0.6)$, $N=T=100$ and $\rho_{fw}=0.6$.
The t-ratios for $\gamma_{\mathbf{R},2}$ and $\beta$ based on $\hat{\boldsymbol{\delta}}$, $\hat{\boldsymbol{\delta}}_{bcjk}$, $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}_q}$ and $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}}$ are calculated.
Table \ref{table:sizequantile} reports the size of the significance tests for the second factor and $w_t$, using the critical value from the standard normal distribution.
There are moderate size distortions in the significance tests for the second factor based on $\hat{\gamma}_{bcjk,2}$, and higher size distortion based on $\hat{\gamma}_{bc\hat{\mathbf{H}},2}$.
In view of the size distortion, the two-sided tests are then implemented to compare the size-adjusted power curves using the 95\% quantiles of the absolute values of the t-ratios under the null over the replications as the critical values, which are also reported in the table.
\begin{table}[!htb]
\centering
\caption{The size of the significance tests for the second factor and $w_t$ and the 95\% quantile of the absolute values of the associated t-ratios based on different estimators}
\label{table:sizequantile}
\begin{tabular}{lcccclccc}
\hline
& & Size of the Test & 95\% quantile & & & & Size of the Test & 95\% quantile \bigstrut\\
\cline{1-4}\cline{6-9}
$\hat{\gamma}_2$ & & 6.1\% & 2.06 & & $\hat{\beta}$ & & 7.3\% & 2.14 \bigstrut[t]\\
$\hat{\gamma}_{bcjk,2}$ & & 7.9\% & 2.15 & & $\hat{\beta}_{bcjk}$ & & 5.1\% & 1.97 \\
$\hat{\gamma}_{bc\hat{\mathbf{H}}_q,2}$ & & 6.4\% & 2.04 & & $\hat{\beta}_{bc\hat{\mathbf{H}}_q}$ & & 5.2\% & 1.97 \\
$\hat{\gamma}_{bc\hat{\mathbf{H}},2}$ & & 10.0\% & 2.28 & & $\hat{\beta}_{bc\hat{\mathbf{H}}}$& & 5.2\% & 1.97 \bigstrut[b]\\
\hline
\end{tabular}
\end{table}
The estimated size-adjusted power curve for $\gamma_{\mathbf{R},2}$ and $\beta$ is shown in Figures \ref{fig:pwcrv}.
From Figure \ref{fig:pwcrv.g} we can see that the power curves of the significance tests for the second factor based on the three bias-corrected estimators are virtually identical, while the power curve of the test based on the least squares estimator is asymmetric most likely due to the bias in the parameter estimates.
The similarity of the power curves for the bias-corrected estimators breaks down for testing the significance of the observed factor $w_t$, which is shown in Figure \ref{fig:pwcrv.b}.
The power curve for the significance test based on the least squares estimator is biased, with the curve shifted substantially to the left.
The bias correction towards the `parameters' with the data dependent rotations $\hat{\mathbf{H}}_q$ and $\hat{\mathbf{H}}$ mitigates the bias of the power curve, but the location shift still remains. In contrast, the jackknife bias correction successfully corrects the bias and restores the symmetry of the power curve at $\beta=0$.
To summarize the power curve analysis, it is generally recommended to use the bias-corrected estimators for the significance test, and the split-panel jackknife estimator seems to be the most reliable among the bias-corrected estimators considered.
The recommendation does not apply in the special case of joint significance tests for all latent factors, since there is no bias in the least squares estimator under the null hypothesis; to see this, substitute $\boldsymbol{\gamma}^* = \mathbf{0}$ into the results in Theorems \ref{thm:bias_Hhat}-\ref{thm:bias_Hw3} and all asymptotic biases disappear.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/powercurve_gamma2.pdf}
\caption{Second factor}
\label{fig:pwcrv.g}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{newimages/powercurve_beta.pdf}
\caption{Observed factor $w_t$}
\label{fig:pwcrv.b}
\end{subfigure}
\caption{Size-adjusted power curve for testing the significance of the second factor and $w_t$ based on $\hat{\boldsymbol{\delta}}$, $\hat{\boldsymbol{\delta}}_{bcjk}$, $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}_q}$ and $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}}$}
\label{fig:pwcrv}
\end{figure}
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Empirical application}\label{sec:emp}
It is important to test which factors are significant in factor-augmented regressions. This is because the extracted PC factors are ordered by their importance in the covariation of the predictors, which may not correspond to their predictive power for the particular series of interest; see further discussion in \citet{BaiNg2008,BaiNg2009} and \cite{ChengHansen2015}.
We consider the factor-augmented forecast regression of bond yields $y_{t+h}$ on
the extracted factors from a large number of predictor variables $x_{t,i}$ and the
observed predictor, $w_{t}$. We use the dataset used in \cite{LudvigsonNg2009}, which is provided by Sydney Ludvigson's website. The data consists of
the continuously compounded (log) annual excess returns on an $2$-year
discount bond at month $t$, $y_{t+12}$, and a balanced panel of
$i=1,\dots,131$ monthly macroeconomic series at month $t$, $x_{t,i}$, spanning
the period $t=$January 1982,..., December 2002, which are standardized. For
$x_{t,i}$, the Edge Distribution (ED) estimator of \cite{Onatski2010} with $r_{\max}=9$ gives $\hat{r}=3$. Before running the regression, three PC factors, $\mathbf{\hat{f}}_{t}=(\hat{f}_{t1},\hat{f}_{t2},\hat{f}_{t3})^{\prime}$, are extracted from $x_{t,i}$ and the \cite{CochranePiazzesi2005} (observed) factor, $CP_{t}$, which is a linear combination of five forward Treasury yield spreads, is obtained. The correlations between $CP_{t}$ and $\mathbf{\hat{f}}_{t}^{\prime}$
were \{-0.11, 0.14, -0.11\}, respectively, which are significant or
insignificant on the borderline at the 10\% level test. We run a regression of
$y_{t+12}$ on $\{\mathbf{\hat{f}}_{t}^{\prime},CP_{t},1\}$. We also run a
similar regression but on $\{\mathbf{\hat{f}}_{wt}^{\prime},CP_{t},1\}$, where
$\mathbf{\hat{f}}_{wt}$ is extracted from the data variable from which $CP$ is projected out.
\begin{table}[!htb]
\centering
\caption{Prediction regression results of the yield of the two-year maturity bond}
\label{table:emp1}
\begin{tabular}
[c]{lrccccc}\hline
$y_{t+12}$ & & $\hat{f}_{t1}$ & $\hat{f}_{t2}$ & $\hat{f}_{t3}$ & $CP_{t}$ &
$R^{2}$\\\hline
$\boldsymbol{\hat{\delta}}$ & & 0.47*** & 0.13 & -0.01 & 0.36*** & 0.14\\
\multicolumn{1}{r}{} & & (3.24) & (0.92) & (-0.21) & (3.15) & \\
$\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}}$ & & 0.50*** & 0.14 & -0.01 &
0.36*** & 0.14\\
\multicolumn{1}{r}{} & & (3.43) & (1.00) & (-0.20) & (3.15) & \\
$\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_{q}}$ & & 0.47*** & 0.13 & -0.01 &
0.36*** & 0.14\\
\multicolumn{1}{r}{} & & (3.24) & (0.92) & (-0.21) & (3.15) & \\
$\boldsymbol{\hat{\delta}}_{bcjk}$ & & 0.49*** & 0.19 & 0.01 & 0.36*** &
0.14\\
\multicolumn{1}{r}{} & & (3.33) & (1.29) & (0.15) & (3.18) & \\\hline
\multicolumn{1}{r}{} & & $\hat{f}_{w,t1}$ & $\hat{f}_{w,t2}$ & $\hat
{f}_{w,t3}$ & $CP_{t}$ & $R^{2}$\\\hline
$\boldsymbol{\hat{\delta}}_{w}$ & & 0.47*** & 0.13 & -0.02 & 0.35*** & 0.14\\
\multicolumn{1}{r}{} & & (3.22) & (0.92) & (-0.28) & (3.09) & \\
$\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_w}$ & & 0.50*** & 0.14 & -0.02 &
0.35*** & 0.14\\
\multicolumn{1}{r}{} & & (3.41) & (1.00) & (-0.28) & (3.09) & \\
$\boldsymbol{\hat{\delta}}_{bc\hat{\mathbf{H}}_{q,w}}$ & & 0.47*** & 0.13 & -0.02 &
0.35*** & 0.14\\
\multicolumn{1}{r}{} & & (3.22) & (0.92) & (-0.28) & (3.09) & \\
$\boldsymbol{\hat{\delta}}_{bcjk,w}$ & & 0.49*** & 0.18 & 0.00 & 0.34*** &
0.14\\
\multicolumn{1}{r}{} & & (3.32) & (1.27) & (0.03) & (3.03) & \\\hline
\end{tabular}
\begin{minipage}{10cm}
\vspace{0.1cm}
\vspace{0.1cm}
\small Notes: Values in parentheses are t-ratio using HAC s.e. *, **, *** indicate significant at 10, 5 and 1\% level respectively.
\end{minipage}
\end{table}
The estimation results are summarized in Table \ref{table:emp1}.
The value in the parentheses is the t-ratio based on Newey-West heteroskedasticity and autocorrelation consistent (HAC) standard errors with the threshold being the integer part of $T^{1/4}$.
As can be seen, the augmented regressions with $\hat{\mathbf{F}}$ and $\hat{\mathbf{F}}_w$ give very similar results. The first factor and CP are significant at the 1\% level, while the test fails to reject $H_0 : \gamma_2 = 0$, but the magnitude of the t-statistic is larger with the jackknife estimator due to the positive upward correction.
It may be of practical interest to find out whether the orthonormal latent factors have the same explanatory power of $y_{t+12}$. In this application, it seems
reasonable to consider the pair $(f_{t1},f_{t2})$ for such a comparison.
Accordingly, we have computed the value of $\hat{\gamma}_{1} -\hat{\gamma}_{2}$ for different estimates and the corresponding t-test statistics for $H_{0}:\gamma_{1}=\gamma_{2}$ versus $H_{1}:\gamma_{1}\neq\gamma_{2}$. Very similar estimates are obtained for the transformed model, $\mathbf{M}_w \mathbf{X}$, which seems reasonable given that the correlations between $\hat{\mathbf{F}}$ and $CP$ are not strong.
The results are summarized in Table \ref{table:emp2}.
As can be seen, the test of the equal explanatory power of the first and the second factors is rejected with the estimators $\hat{\boldsymbol{\delta}}$, $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}}$ and $\hat{\boldsymbol{\delta}}_{bc\hat{\mathbf{H}}_q}$, but not rejected with $\hat{\boldsymbol{\delta}}_{bcjk}$. This is due to the jackknife bias-correction which makes the estimate of the difference $\gamma_1 - \gamma_2$ smaller.
Very similar comments apply to the results for the transformed model, $\mathbf{M}_w \mathbf{X}$.
\begin{table}[!htb]
\centering
\caption{The estimated difference $\gamma_1 - \gamma_2$ and the t-test statistic for $H_0 : \gamma_1 = \gamma_2$}
\label{table:emp2}
\begin{tabular}[c]{cccccc}
\hline
& \multicolumn{1}{c}{} & $\boldsymbol{\hat{\delta}}$ & $\boldsymbol{\hat
{\delta}}_{bc\mathbf{\hat{\mathbf{H}}}}$ & $\boldsymbol{\hat{\delta}}_{bc\mathbf{\hat
{\mathbf{H}}}_{q}}$ & $\boldsymbol{\hat{\delta}}_{bcjk}$\\\hline
$\gamma_{1}-\gamma_{2}$ & & 0.34* & 0.36* & 0.34* & 0.30\\
& & (1.75) & (1.84) & (1.75) & (1.55)\\\hline
& & $\boldsymbol{\hat{\delta}}_{w}$ & $\boldsymbol{\hat{\delta}
}_{bc\mathbf{\hat{\mathbf{H}}}_w}$ & $\boldsymbol{\hat{\delta}}_{bc\mathbf{\hat{\mathbf{H}}}
_{q,w}}$ & $\boldsymbol{\hat{\delta}}_{bcjk,w}$\\\hline
$\gamma_{w1}-\gamma_{w2}$ & & 0.35* & 0.37* & 0.35* & 0.31\\
& & (1.78) & (1.88) & (1.78) & (1.61)\\\hline
\end{tabular}
\begin{minipage}{10cm}
\vspace{0.1cm}
\vspace{0.1cm}
\small Notes: Values in parentheses are t-ratio using HAC s.e. *, **, *** indicate significant at 10, 5 and 1\% level respectively.
\end{minipage}
\end{table}
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Conclusion}\label{sec:con}
In this paper we have studied the asymptotic bias of the least squares (LS) estimator for the factor-augmented model, $y_{t+h}={\boldsymbol{\gamma}^*}'\mathbf{f}_t^*+\boldsymbol{\beta}' \mathbf{w}_t+\epsilon_{t+h}, t=1,\dots,T$,
replacing the latent factor $\mathbf{f}_t^*$ with the principal component (PC) estimator $\hat{\mathbf{f}}_t$ extracted from a large set of predictors, $\{x_{t,i}\}_{i=1}^N$. Unlike the existing literature, we allow for the predictors $x_{t,i}$ follow more general weak factor (WF) models, in which the $r$ largest eigenvalues of the sample covariance matrix of $x_{t,i}$ may diverge at different rates, $N^{\alpha _{k}}$, $0<\alpha _{k}\leq 1$, $k=1,\dots,r$. The literature typically assumes the strong factor (SF) model, in which $\alpha_1=\dots=\alpha_r=1$.
As discussed in \cite{BaiNg2023} and \cite{jiang2023revisiting}, there are choices of rotation matrices $\mathbf{R}$ for approximations $\hat{\mathbf{f}}_t = {\mathbf{R}}'\mathbf{f}_t^* + o_p(1)$. Accordingly, the first term of the augmented model is approximated by ${\boldsymbol{\gamma}^*}'\mathbf{f}_t^*=\boldsymbol{\gamma}_{{\mathbf{R}}}'\hat{\mathbf{f}}_t + o_p(1)$, where $\boldsymbol{\gamma}_{{\mathbf{R}}}={\mathbf{R}}^{-1}\boldsymbol{\gamma}^*$, hence resulting in different ``parameters'' $\boldsymbol{\delta}_{\mathbf{R}}=(\boldsymbol{\gamma}_{{\mathbf{R}}}',\boldsymbol{\beta}')'$ estimated by the least squares estimator $\hat{\boldsymbol{\delta}}=(\hat{\boldsymbol{\gamma}}',\hat{\boldsymbol{\beta}}')'$, obtained by regressing $y_{t+h}$ on $(\hat{\mathbf{f}}_t',\mathbf{w}_t')$.
Replacing the regressor $\mathbf{f}_t^*$ with $\hat{\mathbf{f}}_t$ generally results in non-zero correlation between the regressor $(\hat{\mathbf{f}}_t',\mathbf{w}_t')$ and the replacement error, which can lead to an asymptotic bias of $\hat{\boldsymbol{\delta}}$.
We have studied the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\mathbf{R}})$
for three different choices of $\mathbf{R}$ for the approximation: $\hat{\mathbf{H}}$ which is data ($x_{t,i}$) dependent and commonly used for the approximation in the literature including \cite{BaiNg2006}, \cite{GoncalvesPerron2014,gonccalves2020bootstrapping}; another data dependent matrix $\hat{\mathbf{H}}_q$, whose estimation error appears to be orthogonal to $\hat{\mathbf{f}}_t$; and $\mathbf{H}$, which is the population matrix $\mathbf{H}$ of $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$.
Our study has shown that if $\sqrt{T}/N^{(3\alpha_r - \alpha_1)/2} \to c_1\in[0,\infty)$ as $N,T\to\infty$, $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ has an asymptotic bias for WF models, which generalizes the results in \cite{GoncalvesPerron2014,gonccalves2020bootstrapping} for SF models. It is shown that the asymptotic bias expression for WF models is more complicated than for SF models, because the former depends on how many exponents $(\alpha_1,\dots,\alpha_r)$ are the same as $\alpha_1$ and $\alpha_r$.
We have also shown that if $\sqrt{T}/N^{\alpha_r} \to c_2\in[0,\infty)$ as $N,T\to\infty$, $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ has an asymptotic bias that is generally smaller in magnitude than that of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ and also the convergence rate is generally not slower than that of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$, and even the former is faster when $\alpha_1\neq\alpha_r$.
The structure of the asymptotic bias depends on how many exponents $(\alpha_1,\dots,\alpha_r)$ are the same as $\alpha_r$.
Importantly, the asymptotic biases of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ and $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ are parametrically estimable; thus, analytical bias corrections are feasible in practice.
Moreover, it turns out that the asymptotic bias of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ disappears completely when $\mathbf{w}_t$ and $\mathbf{f}_t^*$ are uncorrelated. To exploit this property, we propose to extract the factor from the predictors ($x_{t,i}$) after projecting out the observable factor, $\mathbf{w}_t$.
We have also studied the asymptotic bias with the population rotation matrix, $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}^0)$ with $\boldsymbol{\delta}^0:=\boldsymbol{\delta}_{{\mathbf{H}}}$, where $\boldsymbol{\delta}^0=(\boldsymbol{\gamma}^{0\prime},\boldsymbol{\beta}')'$, $\boldsymbol{\gamma}^0:= {\mathbf{H}}^{-1}\boldsymbol{\gamma}^{*}$. It is shown that if $\sqrt{T}/N^{(3\alpha_r - \alpha_1)/2} \to c_1\in[0,\infty)$ as $N,T\to\infty$, $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}^0)$ has an asymptotic bias which, unlike the biases with $\hat{\mathbf{H}}$ and $\hat{\mathbf{H}}_q$, cannot be estimated parametrically.
In view of this, we have proposed to use a subsampling method, called a split-panel jackknife bias correction, which is generally less computationally expensive than bootstrapping, while allowing for more general cross and serial correlations in $e_{t,i}$.
The finite sample evidence has shown that $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}_q})$ has the least bias, the least size distortion of t-tests, with the smallest standard errors, while the performance of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\hat{\mathbf{H}}})$ is far worse than others throughout the design. Our preferred jackknife bias-corrected estimator with respect to the parameter $\boldsymbol{\delta}^0$ comes in second, closely following the performance of the bias-corrected estimator with respect to $\boldsymbol{\delta}_{\hat{\mathbf{H}}_q}$.
In addition, the empirical power curve of the significance test shows that, on balance, the split-panel jackknife estimator appears to be the most reliable among the estimators compared.
We apply the bias-corrected estimator to the factor-augmented forecast regression of bond yields $y_{t+h}$ on the factors extracted from 131 monthly macroeconomic series and one observed predictor, the \cite{CochranePiazzesi2005} factor, over the period January 1982 to December 2002. The results show that the jackknife appears to effectively correct the bias of the LS estimator, thus providing more reliable inference.
A couple of implications follow from the results of this paper. First, the different approximation for $\hat{\mathbf{f}}_t$ with different rotation matrices that are asymptotically equivalent implies that the LS estimator $\hat{\boldsymbol{\delta}}$ estimates different ``parameters'' and they may have different asymptotic biases. Therefore, the researchers should clarify which ``parameter'' they are estimating with $\hat{\boldsymbol{\delta}}$. This is very important because in finite samples the ``parameters'' for different rotation matrices can take very different values, as shown in Figure \ref{fig:mean_gamma_hats}. Our recommendation is to primarily consider $\boldsymbol{\delta}^0$ as the parameter estimated by $\hat{\boldsymbol{\delta}}$.
Second, the bias in $\hat{\boldsymbol{\delta}}$ can be significant and should not be ignored in practice. The empirical results in Table \ref{table:emp1} illustrate the effectiveness of the jackknife bias correction of $\hat{\boldsymbol{\delta}}$ relative to $\boldsymbol{\delta}^0$. As shown in Table \ref{table:emp2}, the test results for parameter restrictions with jackknife bias correction are significantly different from those without bias correction.
Note that, as our theory tells that the bias correction is not necessary
in the special case of joint significance testing of all factors, since there is no bias in the LS estimator under the null hypothesis of $\boldsymbol{\gamma}^*=\mathbf{0}$ (all asymptotic biases in Theorems \ref{thm:bias_Hhat}-\ref{thm:bias_Hw3} disappear).
Finally, while we have focused on the analysis of the asymptotic bias of $\hat{\boldsymbol{\delta}}$ in this paper, it is also of great interest to extend our analysis to bootstrapping $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\mathbf{R}})$ for different choices of $\mathbf{R}$. In particular, the advantage of the bootstrapping is that it can provide a higher-order approximation not only for the asymptotic bias but also for the distribution of $\sqrt{T}(\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}_{\mathbf{R}})$. This line of research is pursued in a companion paper, \cite{jiang2024bootstrap}, to which interested readers may wish to refer.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}*{Acknowledgment}
We are grateful to Jia Chen, Naoko Hara, Yohei Yamamoto and Yang Zu for helpful discussions and useful comments.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}*{Funding}
This work was supported by JSPS KAKENHI (grant numbers 21H00700, 21H04397, 23K25501 and 24K16343).
\bibliographystyle{chicago}
\bibliography{references_wfr}
\newpage