A Robust Residual-Based Test for Structural Changes in Factor Models

A Robust Residual-Based Test for Structural Changes in Factor Models

{8pt} {8pt}

{7pt} {7pt}

Introduction

As we are embracing the era of data rich environment, factor models have been very popular in the past three decades or so. See Jolliffe and FLL2021 for comprehensive reviews on a variety of applications of factor models. One of the most attractive features is possibly the low rank representation, which usually imposes that a set of high dimensional observations can be driven by a series of low dimensional vectors. More often than not, the low dimensional vectors are the foundation of many factor augmented analyses (see, e.g., STOCK2016415, and references therein).

To obtain these low dimensional vectors properly, extensive efforts have been devoted to estimating the number of factors, viz., the length of aforementioned low dimensional vectors. See, e.g., BN2002, onatski2010, NW2011, LM2012, ahn2013eigenvalue, among others. Notably, all of these works are based on the conventional factor models where the loadings are assumed to remain constant over time. As economists realize that the relationships between economic and financial variables may suffer from structural changes over time, the factor loadings, which capture the relationships between these variables and the latent common factors, could be varying over time or some state variables. For this reason, various methods have been proposed to test for the structural changes in factor models or to estimate the factor models with breaks. For example, breitung2011testing, chen2014detecting, han2015tests, yamamoto_tanaka2015, su2017time, su2020testing, BKW2021, fu2023testing, and BDH2024 consider various methods to detect structural changes in large dimensional factor models. Except su2017time, su2020testing and fu2023testing who allow for smooth changes, almost all the other papers in the literature consider one or multiple abrupt changes. We provide a detailed summary of these works in Section (ref) for the purpose of comparison. In terms of estimation, cheng2016shrinkage study the estimation of a factor models with a one-time break; ma_su2018 and BKW2021 consider the estimation of factor models with multiple abrupt breaks; su2017time and pelger_xiong2022 consider the kernel estimation of factor models with loadings changing over time and stated variables, respectively; fu_su_wang2024 study time-varying factor-augmented regression models.

In this paper, we are also interested in testing whether the loadings change over time or not. We aim to propose a new test statistic that has a superb power in detecting smaller local alternatives than almost all existing tests by utilizing richer information along both cross-sectional and time dimensions, and remains robust against both cross-sectional and serial dependence. One popular approach of detecting an abrupt change is via the use of a Lagrange multiplier (LM) type of test statistic. For example, chen2014detecting, BKW2021 and BDH2024 all consider LM-type of tests based on the information of the unknown factors. When doing so, the test statistic can only detect local alternatives that converge to the null at a slow rate (e.g., $T^{-1/2}$) associated with the time dimension $\left( T\right) $ alone, and it is generally not robust to misspecified factor numbers (e.g., see p. 359 in BKW2021). Another popular approach aims at detecting smooth structural changes and remains powerful against abrupt breaks too. The tests of su2017time,su2020testing and fu2023testing are capable of examining the parametric loadings against the nonparametric generalization via a smooth function. In general, such tests are able to utilize the information along both dimensions of the residuals, thereby detecting local alternatives converging to the null at a rate faster than $T^{-1/2}.$ In particular, su2017time,su2020testing consider a kernel-based smoothing nonparametric test for structure changes that can detect local alternatives converging to the null at the rate $ N^{-1/4}T^{-1/2}h^{-1/4},$ where $N$ and $T$ denote cross-sectional and time dimensions, respectively, and $h$ is a bandwidth parameter; fu2023testing consider a nonsmoothing test\footnote{See Chapter 13 in li2007nonparametric for the discussions on smoothing and nonsmoothing tests.} for the structural changes in factor models based on the discrete Fourier transform and their test can detect local alternatives at a rate $(NT)^{-1/2}$. But to work out the asymptotic distributional theory, early research typically assumes either cross sectional independence or serial uncorrelation for the error terms $\left\{ \varepsilon _{it}\right\} $ in the factor model. The reason is that one often has to deal with the following terms:

where $\left[ L\right] \coloneqq\left\{ 1,\ldots ,L\right\} $ for a given positive integer $L$. When aggregating the information from both the cross-sectional and time dimensions, early literature needs to balance the rate of convergence yielded from the aggregation of the terms in (ref) and the independence conditions imposed on $\left\{ \varepsilon _{it}\right\} $. To the best of our knowledge, no attempt has been made to construct a test statistic using (ref) and to study its limit null distribution while allowing for dependence along both dimensions of $\left\{ \varepsilon _{it}\right\} $. Despite the technical challenge in dosing so, the issue is a realistic one. As CP2015 remark, the weak dependence along both dimensions of $\left\{ \varepsilon _{it}\right\} $ is likely to be the rule rather than the exception practically.

In this paper, we focus on a residual-based test to examine whether one should go beyond a traditional type of factor model practically while accounting for dependence over both dimensions of $\left\{ \varepsilon _{it}\right\} $. We shall derive the asymptotic properties of the newly proposed test statistic, and conduct extensive numerical studies to confirm our theoretical findings. To proceed, we briefly review the literature on residual-based model specification tests. The basic idea is to use the estimated residuals (typically from a parametric or semiparametric model under the null) to construct a test statistic, in which

Under the null, the estimated residuals are typically sufficiently close to the true error terms so that the limiting distribution of a suitably constructed residual-based test statistic can be established. Under the alternative, the estimated residuals contain a term that reflects the deviation from the null hypothesis and thus contribute to the power of the test statistic. Along this line of research, early works focus on examining different types of time series models (e.g., DK1979, PP1988). Gradually, the idea is extended to justify different model specifications for a variety of data types such as cross-sectional data ( Li1994, JOHNXUZHENG1996263), time series data (GKLT2009), and panel data ( su2015specification; WSX2022).

To sum up, our main contributions are as follows. First, assuming the number of factors is correctly specified, we establish the asymptotic distribution of a residual-based test statistic under the null by allowing for dependence along both dimensions of $\left\{ \varepsilon _{it}\right\} $. As mentioned above, this is a highly technically challenging issue, and thus a new central limit theorem is given for the quadratic forms of panel data with dependence over both dimensions, thereby filling a gap in the literature. This result is of independent interest and can be used in various kernel-based specification tests (e.g., su2015specification,chen2018nonparametric). Second, we explore the local power of the proposed test and show that it can detect local alternatives converging to the null at a rate $(NT)^{-1/2}h^{-1/4},$ faster than almost all existing tests in the literature but fu2023testing who consider a nonsmoothing test. It is well known a nonsmoothing test (like that of fu2023testing) tends to be more powerful than a smoothing test (like ours) under low frequency local alternatives whereas\ a smoothing test can be more powerful than a nonsmoothing test under high frequency local alternatives as studied in rosenblatt1975quadratic. See the discussion on p. 400 of li2007nonparametric. Third, we show the robustness of our test by proving the limiting distributions under the null and local alternatives continue to hold true when the number of factors is over-specified. Fourth, we also study the global power behavior of the proposed test and our investigation helps to explain where the residual-based tests gain their power. Finally, we conduct numerical studies to examine the finite sample performance of the test and the results corroborate the theoretical findings.

The rest of the paper is organized as follows. Section (ref) presents the hypotheses and the test statistic. Section (ref) investigates the test statistic under the null, a sequence of Pitman local alternatives, and a global alternative. Section (ref) conducts numerical studies using simulated data and real data examples. Section (ref) concludes. The proofs of the main results are relegated to the appendix. The online supplement contains the proofs of the technical lemmas used in the proof of the main result, the verification of Assumption (ref), and some additional simulation results.

Notation. For a real-valued matrix $\mathbf{A=}\{a_{ij}\}$, let $ \Vert \mathbf{A}\Vert $ and $\Vert \mathbf{A}\Vert _{2}$ denote the Frobenius norm and the spectral norm, respectively. Let $\Vert \mathbf{A} \Vert _{1}=\max_{j}\sum_{i}|a_{ij}|$ and $\Vert \mathbf{A}\Vert _{\infty }=\max_{i}\sum_{j}|a_{ij}|$. When $\mathbf{A}$ has full column rank, let $ \mathbf{P}_{\mathbf{A}}\coloneqq \mathbf{A}(\mathbf{A}^{\top }\mathbf{A} )^{-1}\mathbf{A}^{\top }$ and $\mathbf{M}_{\mathbf{A}}\coloneqq \mathbf{I}- \mathbf{P}_{\mathbf{A}},$ where $^{\top }$ denotes transpose. When $\mathbf{A }$ is symmetric, we use $\lambda _{\min }(\mathbf{A})$ and $\lambda _{\max }( \mathbf{A})$ to denote its minimum and maximum eigenvalues, respectively. For a random vector $\mathbf{v}$, let $|\mathbf{v}|_{q}\coloneqq \left( E\Vert \mathbf{v}\Vert ^{q}\right) ^{1/q}$. For two constants $a$ and $b$, $ a\wedge b\coloneqq\min (a,b)$, $a\vee b\coloneqq\max (a,b)$, and $a\simeq b$ stands for $a=O(b)$ and $b=O(a)$. For two random variables $c$ and $d$, $ c\asymp d$ stands for $c=O_{P}(d)$ and $d=O_{P}(c)$. The operators $ \rightarrow _{P}$ and $\rightarrow _{D}$ stand for convergence in probability and convergence in distribution, respectively. $\left(N,T\right) \rightarrow \infty $ signifies that $N$ and $T$ diverge to infinity jointly.

The Null Hypothesis and the Test Statistic

In this section we introduce the null hypothesis and the test statistic.

Hypotheses

Consider a series of high dimensional vectors $\{\mathbf{x}_{t}:t\in \left[ T \right] \}$, and assume that $\{\mathbf{x}_{t}\}$ is generated by the following factor model that is potentially time-varying:

where $\pmb{\varepsilon}_{t}=(\varepsilon _{1t},\ldots ,\varepsilon _{Nt})^{\top }$ is a vector of zero-mean error terms, $\pmb{\Lambda}_{t}=( \pmb{\lambda}_{1t},\ldots ,\pmb{\lambda}_{Nt})^{\top }$ is an $N\times r$ time-varying loading matrix, and $\mathbf{f}_{t}$ is an $r\times 1$ latent factor. We focus on the case where $r$ is fixed, and $\left( N,T\right) \rightarrow \infty $. To examine whether $\pmb{\Lambda}_{t}$'s evolve over time, the null hypothesis is specified as follows:

where $\pmb{\Lambda}$ is a conformable constant matrix. The alternative hypothesis $\mathbb{H}_{1}$ is the negation of $\mathbb{H}_{0}.$ Obviously, under $\mathbb{H}_{0},$ the factor model in ((ref)) degenerates to the conventional time-invariant factor model that has been widely studied in the literature. Under $\mathbb{H}_{1},$ the loadings $\pmb{\lambda}_{it}$ must change over time for some $i\in \left[ N \right]$.

The Test Statistic

Under $\mathbb{H}_{0}$, we have

where $\mathbf{X=(\mathbf{x}}_{1},\ldots,\mathbf{\mathbf{x}}_{T}\mathbf{)} ^{\top }\mathbf{,}$ $\mathbf{F}=(\mathbf{f}_{1},\ldots ,\mathbf{f} _{T})^{\top },$ and $\mathcal{E}$ is defined accordingly. Following the lead of BN2002, we impose the identification condition $\frac{1}{T}\mathbf{ F}^{\top }\mathbf{F}=\mathbf{I}_{r}$ to conduct the standard principal component analysis (PCA):

where $\widehat{\mathbf{V}}$ is an $r\times r$ diagonal matrix including the first $r$ largest eigenvalues of $\mathbf{X}\mathbf{X}^{\top }$ on the main diagonal, and $\frac{1}{T}\widehat{\mathbf{F}}^{\top }\widehat{\mathbf{F}}= \mathbf{I}_{r}$. Accordingly, we have the residual matrix $\widehat{\mathcal{ E}}=\mathbf{X}-\widehat{\mathbf{F}}\widehat{\pmb{\Lambda}}^{\top }=\mathbf{M} _{\widehat{\mathbf{F}}}\mathbf{X},$ where $\widehat{\pmb{\Lambda}}=\frac{1}{T }\mathbf{X}^{\top }\widehat{\mathbf{F}}$. Having $\widehat{\mathcal{E}}=( \widehat{\mathcal{E}}_{1},\ldots ,\widehat{\mathcal{E}}_{N})=\{\widehat{ \varepsilon }_{it}\}_{T\times N}$ in hand, we construct the residual-based test statistic as follows:

where $\mathbf{K}_{h}=\{\frac{1}{h}K(\frac{t-s}{Th})\}_{T\times T}$ with $t,s\in [ T]$ and $h\coloneqq h_{T}$ is a bandwidth parameter. The test statistic in (ref) has a typical form of the residual-based test, and the expression further concurs with (ref) that one has to account for dependence yielded by $\varepsilon _{it}\varepsilon _{js}$ appearing in the leading term in the expansion of $L_{NT}.$

The Asymptotic Properties of the Test Statistic

In this section, we establish the limiting distribution of $L_{NT}$ under the null, and study its asymptotic local and global power properties. To facilitate the presentation, we first provide the basic assumptions used for our asymptotic analyses in Section (ref). In Section (ref) we establish the asymptotic distribution of the test statistic under the null, and then present the result with the number of factors over-specified. In Sections (ref) and (ref), we study the power properties of our test under different scenarios, of which one is also helpful to explain where the residual-based tests gain their global power. Finally, we point out that the newly proposed test can also be used to select the number of factors practically.

Basic Assumptions

Note that we can write $\pmb{\Lambda}=(\pmb{\lambda}_{1},\ldots , \pmb{\lambda}_{N})^{\top }$ under $\mathbb{H}_{0}.$ As in su2017time, we will assume that the loadings are nonrandom without loss of generality. To facilitate the asymptotic analyses, we make the following assumptions.

Assumption (ref) imposes some standard conditions on the kernel function $K\left( \cdot \right) $ and the bandwidth parameter $h.$ In addition, we allow for $N$ and $T$ diverging to infinity at different or same rates. Without loss of generality, Assumption (ref) assumes that the loadings under $\mathbb{H}_{0}$ are non-random; see, e.g., su2017time and fu2023testing and the references therein. It is easy to see that the asymptotic analyses remain valid when the loadings are random and independent of the factors and error terms, and satisfy some moment conditions as in BN2002.

Assumption (ref) imposes weak dependence conditions on $\{ \mathbf{f}_{t}\}$, allowing for serial correlations and conditional heteroscedasticity of unknown forms. The strength of temporal dependence is defined by the functional dependence measure of wu2005nonlinear, which nests many stationary processes as special cases (e.g., tong1990non; wu2005nonlinear). Here, we only require an algebraic decay rate of the temporal dependence. Intuitively, the physical dependence coefficients $\{\lambda_4^f(m)\}$ quantify the dependence of outputs $\{\mathbf{f}_{m}\}$ on $\mathbf{e}_0$.

Assumption (ref) regulates the decay rate of temporal dependence after taking weighted average cross-sectionally. By doing so, $\varepsilon _{it}$'s are allowed to exhibit cross-sectional dependence of various unknown forms. The constants \textquotedblleft $\alpha $\textquotedblright\ and \textquotedblleft $\beta $\textquotedblright\ play a role similar to the $\alpha $-mixing or $\beta $-mixing coefficients in the time series literature and we only require an algebraic decay rate of the temporal dependence. Assumption (ref)(c) requires the underlying data generating process to satisfy the following moment restriction (see Lemma (ref)(c) for details)

As a high level assumption, the above condition has been widely adopted in the literature on factor models (see, e.g., Assumption F.1 in bai2003inferential, Assumption A.1 (vii) in su2017time and Assumption A.7 in fu2023testing). Intuitively, Assumption (ref)(c) requires that $\mathbf{f}_{t}\frac{1}{\sqrt{N}} \sum_{i=1}^{N}\varepsilon _{it}\varepsilon _{is}$ and $\mathbf{e}_{t-m}$ are asymptotically independent as $m\rightarrow \infty $ for any fixed $s.$ In Section (ref) of the online supplement, we provide an example to show that Assumption (ref) is fulfilled by a well known data generating process.

Asymptotic Null Distribution

Define

where $\overline{\varepsilon }_{t,a}\coloneqq \frac{1}{\sqrt{N}} \sum_{i=1}^{N}(1-a_{i})\varepsilon _{it}$, $a_{i}\coloneqq \pmb{\lambda}_{i}^{\top }( \frac{1}{N}\pmb{\Lambda}^{\top }\pmb{\Lambda})^{-1}\overline{\pmb{\lambda}}$ , and $\overline{\pmb{\lambda}}\coloneqq\frac{1}{N}\sum_{i=1}^{N}\pmb{\lambda}_{i}$ . Clearly, $\sigma _{\varepsilon ,a}^{2}$ denotes a kernel-weighted version of the long-run variance of ${\varepsilon }_{it}.$ We are now ready to present the first main result of this paper.

As indicated in Theorem (ref), in order to achieve an asymptotically nondegenerate normal distribution, the long-run variance $ \sigma _{\varepsilon ,a}^{2}$ cannot be 0. This immediately rules out the following special case:

where $f_{t}$ is a common trend function with loadings given identically by 1. In this case, we have $1-a_{i}=0$ for all $i\in \lbrack N]$. Notably, the term $1-a_{i}$ is equivalent to $B_{i}$ of fu2023testing, which also suffers the issue raised here. The reason of having a structure as in $1-a_{i}$ is due to the use of projection matrix involved in the PCA. Our test has the usual asymptotic normal distribution when one has idiosyncratic trends as follows

or equivalently in scalar notations: $x_{it}=\lambda _{i}f_{t}+\varepsilon _{it},$ where $\lambda _{i}$ is not a constant across $i.$ Interestingly, the common trend issue has been addressed in another literature. Various tests have been proposed to test for common trend in the literature; see e.g., zhang2012testing and WSX2022. As discussed in zhang2012testing, high-dimensional time series rarely share a common trend in the real world and idiosyncratic trending behavior is a norm rather than the exception. For the factor model, if one wishes, one can first eliminate the possibility of having a common time trend in the dataset, say, by using the test of WSX2022. But this is subject to the usual pretesting issue. For this reason we do not recommend testing the presence of common trends before using our test.

Building on Theorem (ref), we investigate the case when the number of factors is over-specified. That said, suppose that when conducting PCA, $ \widetilde{r}\ (\geq r)$ is specified as the number of factors. Then the eigenvalue problem in ((ref)) becomes

where $\widetilde{\mathbf{V}}$ is an $\widetilde{r}\times \widetilde{r}$ diagonal matrix including the first $\widetilde{r}$ largest eigenvalues of $ \mathbf{X}\mathbf{X}^{\top }$ on the main diagonal, and $\frac{1}{T} \widetilde{\mathbf{F}}^{\top }\widetilde{\mathbf{F}}=\mathbf{I}_{\widetilde{r }}$. Accordingly, we have

where $\widetilde{\pmb{\Lambda}}=\frac{1}{T}\mathbf{X}^{\top }\widetilde{ \mathbf{F}}$ and $\widetilde{\mathcal{E}}=(\widetilde{\mathcal{E}} _{1},\ldots ,\widetilde{\mathcal{E}}_{N})=\{\widetilde{\varepsilon } _{it}\}_{T\times N}$. With a bit abuse of notation, we continue to use $ L_{NT}$ to denote the test statistic:

We summarize the asymptotic property of $L_{NT}$ under the null in the following theorem.

Theorem (ref) infers that when the number of factors is over-specified, we are still able to achieve the same limit null distribution as in Theorem (ref).

To implement our test, one needs to estimate the long-run variance $\sigma _{\varepsilon ,a}^{2}$. Define a panel heteroskedasticity and autocorrelation consistent (HAC) variance estimator by

where $\widehat{\sigma }_{\varepsilon ,a,k}^{2} \coloneqq \frac{1}{T} \sum_{t=1}^{T-|k|} \widehat{\overline{\varepsilon }}_{t}\widehat{\overline{ \varepsilon }} _{t+|k|}$, $\widehat{\overline{\varepsilon }}_{t}\coloneqq \frac{1}{ \sqrt{N}} \sum_{i=1}^{N}\widetilde{\varepsilon }_{it}$, $a(\cdot )$ is a kernel function, and $l\coloneqq l_{T}$ is a bandwidth parameter for the HAC estimation. We further impose the following assumption.

Assumption (ref) includes a set of typical conditions for the kernel function $a(\cdot)$ and the bandwidth $l$, and we omit the discussion and refer interested readers to andrews1991heteroskedasticity for details. Then the following result holds.

Proposition (ref) infers that $\widehat{L}_{T}$ is asymptotically pivotal under the null. By the first part of Proposition (ref), $ l\simeq T^{1/(2q+1)}$ corresponds to the optimal rate of bandwidth in terms of minimizing the asymptotic mean squared error of $\widehat{\sigma } _{\varepsilon ,a}^{2}.$

Asymptotic Local Power

To study the asymptotic local power property of our test, we consider a sequence of local alternatives:

where $a_{TN}\rightarrow 0$, $\mathbf{g}_{it}=\mathbf{g}_{i}(t/T)$, and $ \int_{0}^{1}\mathbf{g}_{i}(\tau )\mathrm{d}\tau =\mathbf{0}_{r}$ for the purpose of normalization. The term $a_{TN}\mathbf{g}_{it}$ characterizes the departure of $\pmb{\lambda}_{it}$ from the constant $\pmb{\lambda}_{i}$, and our setup allows for both abrupt and smooth structural changes. For example, if we have only one break at $T_{1}$ with $\tau _{1}=T_{1}/T:$

then we can define

with $\mathbf{g}_{i\left( 1\right) } =-\frac{\mathbf{\lambda } _{i\left( 2\right) }-\mathbf{\lambda }_{i\left( 1\right) }}{a_{TN}}\tau _{2}$, $\mathbf{g}_{i\left( 2\right) }=\frac{\mathbf{\lambda }_{i\left( 2\right) }-\mathbf{\lambda }_{i\left( 1\right) }}{a_{TN}}\tau _{1}$ and $\tau _{2}=1-\tau _{1}$ to ensure $\int_{0}^{1}\mathbf{g}_{i}(\tau )\mathrm{d}\tau =\tau _{1}\mathbf{ g}_{i\left( 1\right) }+\tau _{2}\mathbf{g}_{i\left( 2\right) }=0.$ Here, $ a_{TN}$ controls the speed at which $\mathbf{\lambda }_{i\left( 2\right) }- \mathbf{\lambda }_{i\left( 1\right) }$ shrinks to 0.

The following theorem indicates that our test statistic can detect a class of Pitman local alternatives at a departure rate of $(TN)^{-1/2}h^{-1/4}$. In addition, our test is robust to the over-specified number of factors under the local alternatives.

Theorem (ref) indicates that $L_{NT}$ has asymptotically non-negligible power in detecting local alternatives converging to the null at the rate $(TN)^{-1/2}h^{-1/4}$ provided $\mu _{1}>0.$ This rate is much faster than almost all existing tests in the literature by noticing that $ Nh^{1/2}\rightarrow \infty $ as $\left( N,T\right) \rightarrow \infty $. As mentioned in the introduction, existing time-series-based tests, such as breitung2011testing, chen2014detecting, han2015tests, yamamoto_tanaka2015, BKW2021, and BDH2024, have power in detecting local alternatives at rate $T^{-1/2}.$ Similarly, the kernel-based smoothing tests of su2017time and su2020testing can detect local alternatives at the rate $N^{-1/4}T^{-1/2}h^{-1/4}.$ The only test that can detect local alternatives converging to the null faster than ours is the non-smoothing test of fu2023testing, which is based on a discrete Fourier transform and can detect local alternatives at rate $ (NT)^{-1/2}.$ Again, neither the smoothing kernel test of su2017time and su2020testing nor the nonsmoothing test of fu2023testing allows for serial dependence in the error terms.

It is worth mentioning that the result in Theorem (ref) holds for any fixed $\widetilde{r}\geq r.$ That is, the over-specification of the number of factors does not affect the local power property of our test. Intuitively, this is due to the fact that the local alternatives converge to the null at a very fast rate (viz., $(TN)^{-1/2}h^{-1/4}$), and the over-specified $\widetilde{r}-r$ factors are essentially very weak factors so that they cannot be identified by the PCA and their effect is still present through the residuals in the PCA.

Note that the requirement that $\mu _{1}>0$ rules out the case where $E( \mathbf{f}_{t})=0$. We believe that this requirement is not restrictive for practical purpose. For example, in the case of asset pricing (e.g., lettau2020factors), the risk factors, say, in the Fama-French three-factor model (and also the estimated latent pricing factors) usually do not have zero mean due to the risk premium required by investors. The inclusion of an intercept may not address the issue of non-zero mean factors. As a matter of factor, when an idiosyncratic intercept term $\alpha _{i}$ is included in the regression, it can be absorbed into the vector of factor loadings which has the associated factor $\mathbf{1}_{T}$. In this case, the population mean of $\mathbf{f}_{t}$ is surely nonzero.

Below we show that when $E(\mathbf{f}_{t})=0$, our test statistic is able to detect a class of local alternatives at a departure rate of $ N^{-1/2}h^{1/4}. $ This rate converges to zero slower than $ (TN)^{-1/2}h^{-1/4}$, indicating the price we have to pay when all the factors in $\mathbf{f}_{t}$ have zero mean.

Global Power

In this section, we explain the source of power of our test when data include global structural breaks. Additionally, we point out that the newly proposed test can be used to select the number of factors practically.

To accommodate some typical settings of structural breaks such as those in breitung2011testing, we write the loadings as $\mathbb{A}=[ \pmb{\Lambda}_{1},\ldots ,\pmb{\Lambda}_{T}],$ where $\mathbb{A}$ is a $ N\times Tr$ matrix. Let $R\coloneqq(Tr)\wedge N$ and consider the singular value decomposition (SVD): $\mathbb{A}=\mathbf{U}\mathbf{S}\mathbf{V}^{\top } $, where $\mathbf{U}$ and $\mathbf{V}$ are $N\times R$ and $Tr\times R$ unitary matrices, and $\mathbf{S}$ denotes an $R\times R$ diagonal matrix with singular values $s_{TN,1},\ldots ,s_{TN,R}$ along the diagonal in descending order. Then we can rewrite (ref) in a compact form:

where $\mathbb{F}=\mathrm{diag}(\mathbf{f}_{1}^{\top },\ldots ,\mathbf{f} _{T}^{\top })$ is a $T\times Tr$ block diagonal matrix.

To facilitate the development, we add the following assumption.

This assumption encompasses many existing works as special cases. For example, the model (ref) with a single structural break as in breitung2011testing has $r<\mathrm{rank}(\mathbb{A})=J\leq2r.$ Assumption (ref)(a) regulates the factor loadings such that $\Vert \mathbb{A }\Vert /\sqrt{TN}=O(1)$. Thus, we have $\mathrm{tr}(\mathbb{A}^{\top } \mathbb{A}/(TN))=(TN)^{-1} \sum_{j=1}^{R}s_{TN,j}^{2}=O(1)$. Assumption (ref)(b) suggests that the first $J$ singular values diverge to infinity at the rate of $\sqrt{ TN}$.

Now, we rewrite the model in (ref) as

where $\mathbb{A}^{(J)}=\mathbf{U}^{(J)}\mathbf{S}^{(J)}(\mathbf{V}^{(J)} \mathbf{)}^{\top }$, $\mathbb{A}^{(-J)}=\mathbf{U}^{(-J)}\mathbf{S}^{(-J)}( \mathbf{V}^{(-J)})^{\top }$, $\mathcal{E}^{\dagger }=\mathcal{E}+\mathbb{F(} \mathbb{A}^{(-J)})^{\top }$, $\mathbf{U}^{(J)}$ (resp. $\mathbf{V}^{(J)}$) and $\mathbf{U}^{(-J)}$ (resp. $\mathbf{V}^{(-J)}$) contain the first $J$ and the last $R-J$ columns of $\mathbf{U}$ (resp. $\mathbf{V}$) respectively, and $\mathbf{S}^{(J)}$ (resp. $\mathbf{S}^{(-J)}$) contains the first $J$ (resp. the remaining $R-J$) singular values of $\mathbf{S}$. Therefore, (ref) has a representation involving the time-invariant factor loadings only:

where $\mathcal{F}=\sqrt{T}\mathbb{F}\mathbf{V}^{(J)}$ with $\Vert \mathcal{ F }^{\top }\mathcal{F}/T\Vert =O_{P}(1)$, and $\pmb{\Theta}=\mathbf{U}^{(J)} \mathbf{S}^{(J)}/\sqrt{T}$. Hence, the above analysis also explains where the power of LM-type tests come from under a general setup since these tests (e.g., chen2014detecting) try to test if the covariance matrix of the larger finite-dimensional set of estimated factors is stable.

Define the $J\times J$ symmetric positive definite matrix:

Let $\mathbb{V}$ be an $\widetilde{r}\times \widetilde{r}$ diagonal matrix containing the $\widetilde{r}$ largest eigenvalues of $\pmb{\Delta}$ in descending order and $\pmb{\Upsilon}$ be a $J\times \widetilde{r}$ eigenvector matrix of $\pmb{\Delta}$ corresponding to the $\widetilde{r}$ eigenvalues in $\mathbb{V}$. The following theorem studies the asymptotic behavior of $L_{NT}$ under the global alternative.

We make three comments on Theorem (ref). First, Theorem (ref) implies that $\widehat{L}_{T}$ defined in ((ref)) generally diverges to infinity at a rate $NTh^{1/2}$ if one selects $ \widetilde{r}<J$ under the alternative. Second, Theorem (ref) explains where the power of residual-based tests comes from. In the presence of $k$ breaks, $J\coloneqq$rank$\left( \mathbb{A}\right) \leq k\cdot r.$ Therefore as long as one specifies a factor model with $\widetilde{r}<J$ (so that the true rank of the represented model is under-specified), the residual-based tests would have power in detecting the global alternative. Third, in connection with the results in Section (ref), Theorem (ref) also indicates that we can propose sequential procedure to test the composite null hypotheses

to identify the number of factors practically if a low rank representation is believed to be true. The estimated value of factor number, $\widehat{r}$, will be the value of $j$ when we fail to reject the null at the first time. Of course, in the presence of abrupt structural changes, $\widehat{r}$ would over-estimate the underlying number of factors, $r$, for the original factor model.

To proceed, we emphasize that as explained in fu2023testing, the information criterion of su2017time still offers a consistent estimation of $r$ rather than $J$. We therefore suggest that one always turns to the information criterion of su2017time when estimating $r$ to avoid any power loss.

Simulating the Critical Value

Note that $\widehat{L}_{T}$ is asymptotically pivotal under $\mathbb{H}_{0}$ defined in ((ref)). It is well known that kernel-based nonparametric tests may not perform well in finite samples if one uses the normal critical values and they can be sensitive to the choice of bandwidth. To improve the finite sample performance of the test, we propose a simulation-based scheme for the selection of the critical value. A similar procedure has also been adopted by zhang2012inference for the same purpose in the context of testing univariate time-varying regression models.

The algorithm of a simulation-assisted testing procedure is as follows.

By Theorems (ref) and (ref), $\widehat{L}_{T}^{b}$ and $ \widehat{L}_{T}$ have the same asymptotic distribution under the null. This suggests that instead of using the normal critical value, we can obtain the simulated critical value $\widehat{q}_{1-\alpha }$ for $\widehat{L}_{T}$. Under the global alternatives, $\widehat{L}_{T}^{b}$ remains to be asymptotically normal whereas $\widehat{L}_{T}$ diverges to infinity in probability.

Note that in the above simulation procedure, we do not need to use the estimated factors and loadings under the null to generate $x_{it}^{\ast }.$ This is similar in spirit to the fixed-regressor wild bootstrap of hansen2000testing, where one does not need to mimic the exact dependence structure in the bootstrap world and may still obtain a valid bootstrap test that shares the same limit null distribution as the original test statistic. For our test, one may theoretically provide the approximation rate for the simulation procedure by establishing a deep Gaussian approximation theory for a quadratic form of panel data with complex dependence structure, based on which one can judge whether the simulation can provide a better approximation to the finite sample distribution than the asymptotic normal. But this certainly goes beyond the scope of the current paper, and has to be left for future research.

Numerical Studies

In this section, we conduct Monte Carlo simulations and consider a real data example to examine the theoretical findings.

Simulation

In this section, we examine the finite sample performance of the proposed testing procedure via extensive simulation experiments. We also compare our test with the parametric tests of chen2014detecting and han2015tests (referred to as CDG and HI, respectively) designed for a single structural break with an unknown break date in the factor loadings and the nonparametric tests of su2020testing and fu2023testing (referred to as SW and FHW, respectively) that allow for both single or multiple abrupt breaks and smooth changes under the alternative. We refer interested readers to these studies for detailed implementation, which is omitted here to save space.

Throughout the numerical studies, the Bartlett kernel is adopted. When estimating the long-run covariance, we follow stock2020introduction and set the bandwidth $\ell $ to be $\lceil 0.75T^{1/3}\rceil $, where $ \lceil \cdot \rceil $ denotes a ceiling function. For the bandwidth $h$, we use the rule of thumb bandwidth $h=(TN)^{-1/5}$. In Section (ref) of the online supplement, we also examine the performance of our test for different choices of bandwidth sequences, and the results reveal that our test is not sensitive to the choice of bandwidth. In addition, we consider the case with $\widetilde{r}=2,\ldots ,5$ to verify the theoretical results in Sections (ref)--(ref). For comparison, we also examine the performance of the above tests with the number of factors estimated by using su2017time's information criteria $IC_{h1}$.

Size Performance

We first study the size performance of the proposed test. We generate the data by using the following high-dimensional factor model with two common factors

where $\mathbf{f}_{t}=0.5+0.3\mathbf{f}_{t-1}+N(\mathbf{0},\mathbf{I}_{2})$. To examine the size performance, we consider the following designs for the factor loading $\pmb{\lambda}_{it}$ and the error term $\varepsilon _{it}$:

where each element of $\pmb{\lambda}_{i0}$ is from $N(1,1)$. Here, DGP.S1--DGP.S3 satisfy the null hypothesis of time-invariant factor loading and is used to examine the size performance of the test statistics. Specifically, DGP.S3 examines the performance under both time series autocorrelation (TSA) and cross-sectional dependence (CSD), while DGP.S1 assumes the absence of TSA and DGP.S2 does not include CSD. For each data generating process (DGP), we simulate 1000 data sets with $T=50,100,200$ and $N=50,100$ respectively, and conduct the test as aforementioned. For the simulation-assisted testing procedure, we set $B=1000$.

We first examine the size performance of our test (referred to as PSY). Table (ref) reports the empirical rejection rates based on 1000 replications at the 5% nominal level. We summarize some important findings from these tables. First, PSY tends to be slightly oversized (undersized) with the presence (absence) of time series correlation. But as the sample size $T$ increases, the empirical rejection rates of PSY get closer to the nominal levels. Second, when $\widetilde{r}$ varies from $2$ to $5$, the size of PSY remains stable if the sample size is not so small. This verifies the theoretical findings in Theorems (ref)--(ref). Third, CDG and HI are quite undersized, especially when the factor number is over specified and the sample size is small, which is consistent with the statements in BKW2021. Fourth, SW and FHW performs reasonably well in DGP.S1 in the absence of serial dependence in the error terms as required by their theories, but the sizes of SW and FHW lose control in the presence of serial dependence. Fifth, SW (resp. FHW) tends to be oversized (resp. undersized) when the factor number is over specified and the sample size is relatively small.

Local Power Performance

To examine the local power performance, we first consider the following sequence of local alternatives with smooth structural changes

where $a_{TN}=(TN)^{-1/2}h^{-1/4}$ and $\mathcal{G}(y;\varsigma ,\pmb{\beta} )=[1+\exp (-\varsigma \prod_{l=1}^{p}(y-\beta _{l}))]^{-1}$ is the logistic function with the scale parameter $\varsigma $ and location parameter vector $\pmb{\beta}=[\beta _{1},\ldots ,\beta _{p}]^{\top }$. The above setting has non-monotonic smooth structural changes, which is the same as DGP.P3 in fu2023testing. We also consider the following designs for the error term $ \varepsilon _{it}$:

In each case, $\mathbf{f}_{t}$ and $\pmb{\lambda}_{i0}$ are generated in the same way as in DGP.S1--DGP.S3.

Table (ref) reports the empirical rejection rates for the 5% test for DGP.L1--DGP.L3, which are smooth time-varying factor models with the break size shrinking to zero. We summarize some findings from these tables. First, as $\widetilde{r}$ increases from $2$ to $5$, the local power of PSY generally decreases for fixed $\left( N,T\right) $, which indicates the side effect of overspecifying the umber of factors. In theory (see Theorem (ref), Proposition (ref) and Theorem (ref)), the over-specification of the number of factors does not yield any difference in the asymptotic normal distribution under the null hypothesis and under the local alternatives. However, in finite sample study, especially when the sample size is small, it is expected that the empirical powers are affected. Although the powers are slightly different, the numerical behavior is indeed consistent with the theoretical results established in the paper when the sample size is relatively large. Specifically, for DGP.L1, when $N=100$ and $ T=200$, Table (ref) shows that the empirical rejection rates of PSY test are 0.312, 0.279, 0.301, 0.264 for $\widetilde{r}=2,\ldots ,5$, in which the differences are quite small. Second, CDG and HI seem to have no local power against smooth structural changes, which is not surprising since these two tests are designed to detect the presence of big structural breaks. Third, for DGP.L1, PSY has a larger (resp. smaller) power than SW (resp. FHW), which is consistent with the theory that our test is able to detect a class of local alternatives at the rate $(TN)^{-1/2}h^{-1/4}$, converging to 0 faster (resp. slower) than the rate $ T^{-1/2}N^{-1/4}h^{-1/4} $ (resp. $T^{-1/2}N^{-1/2}$) in SW (resp. FHW). Fourth, when the sample size is small ($T=50$), su2017time's information criteria $IC_{h1}$ tends to overestimate the number of common factors (indeed the effective sample size is $Th$, which is extremely small), in which case PSY test has better local power performance than FHW. Fifth, the local power of SW is very close to one if there is TSA in error terms, which is consistent with the simulation results for DGP.S2--DGP.S3.

To examine the local power performance of the proposed test under abrupt structural breaks, we also consider the following design for $\pmb{\lambda} _{it}$

where $a_{TN}=(TN)^{-1/2}h^{-1/4}$ and $\mathbf{b}=[2,2]^{\top }$. This setup is the same as DGP.P1 in fu2023testing. In addition, similar to DGP.L1--DGP.L3, we also consider the three designs for the error term $ \varepsilon _{it}$ denoted by DGP.L4--DGP.L6 respectively, while $\mathbf{f} _{t}$ and $\pmb{\lambda}_{i0}$ are generated in the same way as in DGP.S1--DGP.S3. Table (ref) reports the empirical rejection rates for the 5% test for DGP.L4--DGP.L6, from which we find similar results as those for DGP.L1--DGP.L3.

Global Power Performance

To examine the global power, we consider the following specification for the time-varying loadings:

where $\mathbf{b}=[0.25,0.25]^{\top }$. We also consider the following designs for the error term $\varepsilon _{it}$:

In each case, $\pmb{\lambda}_{i0}$ and $\mathbf{f}_{t}$ are generated in the same way as DGP.S1--DGP.S3. Clearly, we have a one-time big break in each of the above three DGPs. In this case, each DGP can be written as a time-invariant factor model with three factors since the break in the factor loadings is not idiosyncratic across $i$'s (see han2015tests for more discussions on this type of abrupt structural break). Hence, our test would have power to detect the breaks as long as $\widetilde{r}<3$.

Table (ref) reports the empirical rejection rates for the 5% test for DGP.G1--DGP.G3. We summarize some findings here. First, our test is powerful to reject the null hypothesis when $\widetilde{r}<J=3$, while the empirical rejection rates of both PSY and FHW converge to the 5% nominal level when $\widetilde{r}\geq J$, which is consistent with our theory in Section (ref). In addition, when the sample size is small ($T=50$), su2017time's information criteria $IC_{h1}$ tends to overestimate the number of common factors, in which case PSY and FHW tend to lose power. Second, CDG and HI have desirable power performance if the factor number is specified to be $\widetilde{r}=J=3$ and the sample size is not so small, and their tests lose power when $\widetilde{r}>3$ since their tests are designed to detect the structural breaks in the variance of common factors (c.f., p.39 in chen2014detecting). Third, SW is generally more powerful than PSY and FHW for detecting big structural breaks when the factor number is over-specified, although these empirical rejection rates decrease with the increase of sample size when the factor number is over-specified.

An Empirical Study

In this section, we revisit the U.S. macroeconomic dataset assembled by MN2016. Specifically, we examine the time period from October 2003 to September 2023. Therefore, in total, we have $T=240$ along the time dimension. After removing variables with missing values, $N=127$ macro variables are available for our investigation.

As in the simulation, we implement PSY, CDG, HI, SW and FHW tests. We let $ \widetilde{r}=1,\ldots ,8$ in what follows. To examine the necessity of accounting for CSD and TSA, for any given $\widetilde{r}$, we conduct the PCA and obtain $\{\widetilde{\varepsilon }_{it}\}$ as in Section (ref) . Subsequently, we consider the following two tests:

The results for the CD and LBQ tests are reported in the last two columns in Table (ref). It is obvious that both cross-sectional dependence and serial dependence are present in the residuals.

Now, we focus on the structural change tests in Table (ref) , where $\surd $ and $\times $ denote failing to detect break and detecting breaks, respectively at the 5% nominal level. First, most of the test results suggest the presence of structural changes in the factor models for $ \widetilde{r}\in \{1,\ldots,7\}$. Second, as in the simulations, PSY and SW are largely consistent in terms of rejecting the null with only one exception occurring at $\widetilde{r}=7$. Third, PSY suggests that there are eight common factors in this U.S. macroeconomic dataset if one wants to use a conventional time-invariant factor model to characterize the U.S. macro economy, which is consistent with the argument of MN2016. Fourth, FHW tends to reject the null hypothesis for any $\widetilde{r} \in \{1,\ldots,8\} $, possibly due to its higher power to detect smaller structural changes or its tendency to over-reject in the presence of serial correlations in the error terms. To further confirm the above findings, we use su2017time's information criteria $IC_{h1}$ to determine the number of factors as documented in the end of Section (ref), which indicates two common factors. Putting these numerical evidences together, we conclude that the dataset features a factor model with time-varying loadings.

Conclusion

In this paper, we propose an easy-to-implement residual-based specification testing procedure to detect structural changes in factor models, which is powerful against both smooth and abrupt structural changes with unknown break dates. The proposed test is robust against the over-specified number of factors under the null and local alternatives, and serially and cross-sectionally correlated error processes. A new central limit theorem is derived for the quadratic forms of panel data with dependence over both dimensions, thereby filling a gap in the literature. We establish the asymptotic properties of the proposed test statistic under the null and a sequence of Pitman local alternative, and accordingly develop a simulation-based scheme to select critical value in order to improve finite sample performance. Through extensive simulations and a real-world application, we confirm our theoretical results and demonstrate that the proposed test exhibits desirable size and power in practice.

Some extensions are possible. First, despite the superb asymptotic power properties of the proposed test among all existing kernel-based nonparametric smooth test, its power decreases in the case when all factors have zero mean in the model. It is worthwhile to consider a remedy to avoid such a power loss. Second, we only consider the case of pervasive factors in the paper. Given the recent burgeoning literature on weak factors ( bai2023approximate), it is interesting to extend the current work to allow for the presence of weak factors. We leave these topics for future research.

\setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{section}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{remark}{0} \setcounter{corollary}{0} \setcounter{assumption}{0}

{ {\bf Appendix A} }

{ This appendix includes three sections. Section (ref) summaries the details of the relevant tests in the literature. Section (ref) states some technical lemmas that are used in the proofs of the main results. Section (ref) contains the proofs of the main results in Section (ref). }

Relevant Tests

{ In what follows, we summarize the details of the relevant tests in the literature, and use the following notations. }

{

}

Some Technical Lemmas

Proofs of the Main Results in Section (ref)

{ In this section, we prove the main results in Section (ref). }

{ {3.pt plus 0ex} }

\setcounter{page}{1}

\setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{section}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{remark}{0} \setcounter{corollary}{0} \setcounter{assumption}{0}

{ This appendix includes three sections. Section (ref) contains the proofs of the technical lemmas of Appendix (ref). Section (ref) verifies Assumption (ref) in the paper. Section (ref) contains some additional simulation results. }

Proofs of the Technical Lemmas

{

}

{

}

{

}

{

}

{

}

{

}

{

}

{

}

{

}

{

}

Verification of Assumption (ref)

{ In this section we verify Assumption (ref) via a given example. Consider a high-dimensional MA($\infty $) process $\pmb{\varepsilon} _{t}=\sum_{j=0}^{\infty }\mathbf{B}_{j}\mathbf{e}_{t-j},$ where $\mathbf{B} _{j}$'s are $N\times N$ matrices. We make the following primitive conditions on $\left\{ \mathbf{B}_{j}\right\} ,$ $\{\mathbf{e}_{t-j}\}$ and $\{\mathbf{f }_{t}\}.$ }

{ Below we verify that $\{\pmb{\varepsilon}_{t}\}$ fulfills Assumption (ref)(a)--(b) under Assumption (ref)(a)--(b). Similarly, Assumption (ref)(c) is satisfied under Assumptions (ref)--(ref). }

Verification of Assumption (ref)(a)--(b)

{ By Assumption (ref)(a)-(b) and the independence between $\{\mathbf{e}_{t-j}\}$ and $\{\mathbf{f}_{t}\}$, $E(\boldsymbol{\varepsilon } _{t}|\mathbf{f}_{t})=E(\boldsymbol{\varepsilon }_{t})=0$ and $E[\boldsymbol{ \varepsilon }_{t}\boldsymbol{\varepsilon }_{t}^{\top }]=$$\sum_{j=0}^{\infty }\mathbf{B}_{j}\mathbf{B}_{j}^{\prime }\eqqcolon \Sigma _{\varepsilon }$. By the triangle inequality and the submultiplicative property of the spectral norm, we have

} { This verifies Assumption (ref)(a).}

{ Let $\mathbf{B}_{j}=\{B_{j,kl}\}_{k,l\in \lbrack N]}$ and $\mathbf{v} ^{\top }\mathbf{B}_{j}=(B_{j,\,\centerdot 1},\ldots ,B_{j,\,\centerdot N})$ for any $\Vert \mathbf{v}\Vert =1$. As $\{e_{it}\}$ are independent over $i$ , we can write

where the first inequality follows from some direct calculations, and the second inequality follows from the fact that $\sum_{l=1}^{N-1} \sum_{k=l+1}^{N}B_{j,\,\centerdot l}^{2}B_{j,\,\centerdot k}^{2}\leq (\sum_{l=1}^{N}B_{j,\,\centerdot l}^{2})^{2}$. Based on the above development, we have

This verifies Assumption (ref)(b).}

{Verification of Assumption (ref)(c)}

{ Let $\mathbf{f}_{t,t-m}$ and $\varepsilon _{it,t-m}$ be the coupled version of $\mathbf{f}_{t}$ and $\varepsilon _{it}$ with $\mathbf{e}_{t-m}$ replaced by $\mathbf{e}_{t-m}^{\prime }$. Since $\{\mathbf{f}_{t}\}$ is independent of $\{\pmb{\varepsilon}_{t}\}$, by Jensen inequality, we have

Noting that $|E(\mathbf{f}_{t}-\mathbf{f}_{t,t-m}\,|\, \mathscr{E} _{t-m})|_{2}\leq |\mathbf{f}_{m}-\mathbf{f}_{m}^{\ast }|_{2}=O(m^{-\alpha })$ and $|\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\varepsilon _{it}\varepsilon _{is}|_{2}=O(1)$, it is sufficient to show the order of

}

{ Let $\mathbf{B}_{j,i\cdot}$ denote the $i^{th}$ column of $\mathbf{B} _{j}^\top$. Hence, we have $\varepsilon_{it} = \sum_{j=0}^{\infty}\mathbf{B} _{j,i\cdot}^\top\mathbf{e}_{t-j}$ and $\varepsilon_{it,t-m} = \sum_{j=0,\neq m}^{\infty}\mathbf{B}_{j,i\cdot}^\top\mathbf{e}_{t-j} + \mathbf{B} _{m,i\cdot}^\top\mathbf{e}_{t-m}^\prime$. Then write

}

{ For $E_1$, we have

For $E_{11}$, by using Cauchy-Schwarz inequality and $\frac{1}{N} \sum_{l,i=1}^{N}B_{m,il}^2 \leq \|\mathbf{B}_{m} \|_2^2$, we have

Next, for $E_{12}$, we have

Similarly, we can show that $E_2$ is also $O(\|\mathbf{B}_{m} \|_2)=O(m^{-\alpha})$ for some $\alpha > 3$. }

{ Finally, we consider the third term $E_{3}$. By using arguments as used in the study of $E_{1}$, we have

This completes the verification of Assumption (ref)(c). }

Additional Simulation Results

{ This section reports some additional simulation results.}

Moving Block Bootstrap for the SW and FHW Tests

{ Following the suggestion of a referee, we examine the applicability of the moving block bootstrap (MBB) procedure for SW and FHW tests via simulation studies. The details of the moving block bootstrap (MBB) procedure can be found in Appendix C.6 of fu2023testing, in which they suggest a MBB procedure in order to deal with serially correlated errors. }

{ We generate the data under DGP.S3 where the error terms exhibit both serial and cross-sectional correlations. The simulation results for the 5% test are reported in Table (ref). Overall, we can see that the MBB procedure is not appropriate for obtaining the critical values of the SW and FHW tests. Specifically, the SW and FHW tests display severe size distortions and their sizes converge to zero with the increase of $T$. This result is similar to those in Appendix C.6 of fu2023testing, in which they find that SW and FHW tests tend to over-reject but their sizes decrease with the increasing sample size. }

{

}

Results for Zero-Mean Common Factors

In this section, we examine the performance of our test with zero-mean common factors. As implied by Proposition (ref), our test statistic is of order $Nh^{-1/2}$ under the global alternatives if the common factors have zero mean. Here, we focus on DGP.G1--DGP.G3 except that the common factors are generated by a vector autoregressive model with no intercept.

Table (ref) reports the empirical rejections rates of various tests under DGP.G1--DGP.G3 with zero-mean common factors where the nomial level is 5%. To conclude, although the empirical rejection rates of our test converge to $1$ at a slower rate than SW and FHW tests, our test is still powerful against DGP.G1--DGP.G3 with zero-mean common factors, especially when $r=2$.

{

}

Bandwidth Sensitivity

{ In this section, we investigate the sensitivity of our test to the choice of bandwidth. As before, we adopt the rule of thumb to set the benchmark bandwidth: $h_{C}=(TN)^{-1/5}.$ We further consider $ h_{L}=0.8h_{C} $ and $h_{R}=1.2h_{C}$ respectively in order to examine the sensitivity. Here, we focus on DGP.S3, DGP.L3, DGP.L6 and DGP.G3, all of which allows both TSA and CSD. This allows us examine the performance of the test under both TSA and CSD. }

{ Table (ref) presents the simulation results for $ T=100,200,300$ and $N=100,200$. As we can see clearly from these tables, the size and power performance of our test remains largely consistent with that in the main text. This suggests that our test results are not sensitive to the choice of bandwidth if the sample size is not so small. }

{

}