Supervised Dynamic PCA: Linear Dynamic Forecasting with Many Predictors

onehalfspacing\begin{titlepage} \thispagestyle{empty} \begin{abstract} This paper proposes a novel dynamic forecasting method using a new supervised Principal Component Analysis (PCA) when a large number of predictors are available. The new supervised PCA provides an effective way to bridge the gap between predictors and the target variable of interest by scaling and combining the predictors and their lagged values, resulting in an effective dynamic forecasting. Unlike the traditional diffusion-index approach, which does not learn the relationships between the predictors and the target variable before conducting PCA, we first re-scale each predictor according to their significance in forecasting the targeted variable in a dynamic fashion, and a PCA is then applied to a re-scaled and additive panel, which establishes a connection between the predictability of the PCA factors and the target variable. Furthermore, we also propose to use penalized methods such as the LASSO approach to select the significant factors that have superior predictive power over the others. Theoretically, we show that our estimators are consistent and outperform the traditional methods in prediction under some mild conditions. We conduct extensive simulations to verify that the proposed method produces satisfactory forecasting results and outperforms most of the existing methods using the traditional PCA. A real example of predicting U.S. macroeconomic variables using a large number of predictors showcases that our method fares better than most of the existing ones in applications. The proposed method thus provides a comprehensive and effective approach for dynamic forecasting in high-dimensional data analysis. \noindentKeywords: Dynamic Forecasting, Factor Analysis, Supervised Principal Components, Large-Dimension, LASSO \noindentJEL classification: C22, C23, C38, C53 \end{abstract} \end{titlepage} \section{Introduction} Dynamic forecasting is the process of analyzing time series data using statistical and dynamic modeling techniques to make predictions and to inform strategic decision-making. Recent advances in information technologies make it possible to collect large amounts of data over time, which naturally form high-dimensional time series and characterize many contemporary application problems in business, economics, finance, environmental sciences, and other scientific fields. Extracting useful information from such high-dimensional dependent data to make accurate predictions, dimension reduction becomes a necessity. In the past decades, there have been many dimension-reduction methods developed in the literature, and one of the most successful ones is the Principal Component Analysis (PCA). The PCA method is a versatile tool for dimension reduction and is capable of extracting useful features by transforming the observed variables into a few new uncorrelated factors while retaining much information in the data. See anderson1958introduction and anderson1963asymptotic for its basic concept and properties. The PCA has many applications in various scientific areas, including business, economics, finance, and management. An incomplete list of recent works includes the asset pricing in gu2020empirical, giglio2021asset, and he2022shrinking, the econometric modeling and forecasting in mullainathan2017machine, gaotsay2022, and huang2022scaled, the business and management in kim2005customer, among others. Despite its importance, existing works in these fields mainly focus on employing a static decomposition, where the extracted principal components represent the cross-sectional structure of a contemporaneous panel and ignore the information embedded in the past lagged variables of the series. Arguably the most widely used forecasting method that relies on PCA is the “diffusion-index” approach, which is also known as factor-augmented model developed by Stock2002b,Stock2002a and has received much attention among researchers and practitioners. In a data-rich environment, it has become common that the number of predictors is large and may exceed the number of data points. The factor-based diffusion-index approach is widely applicable under such situations. Specifically, for a high-dimensional observed predictor vector ${\mathbf x}_t\in R^{N}$, the goal is to predict $y_{t+h}$ based on the available information $\{{\mathbf x}_j, j=1,...,t\}$. To avoid the curse of dimensionality, the factor-based approach assumes that the predictors and the variable of interest admit the following structure: \begin{equation} \begin{array}{c} {\mathbf x}_t={\mathbf B}{\mathbf f}_t+{\mathbf u}_t, \\ y_{t+h}=\boldsymbol{\beta}'{\mathbf f}_t+\boldmath$\varepsilon$_{t+h}, \end{array} \end{equation} where ${\mathbf x}_t$ and $y_t$ are both taken to have mean zero, ${\mathbf f}_t$ is an $r$-dimensional vector of latent factors, which may also contain the dynamic factors in Stock2002b,Stock2002a, ${\mathbf B}$ is the associated factor loading matrix, ${\mathbf u}_t$ is an $N$-dimensional idiosyncratic term, which is uncorrelated with the factors ${\mathbf f}_t$, and $\boldsymbol{\beta}$ is an $r$-dimensional slope parameter linking the factors ${\mathbf f}_t$ to $y_{t+h}$ with $h > 0$. Motivated by the success of the diffusion-index approach and to explore further available information in ${\mathbf x}_t$ for predicting $y_{t+h}$, we consider the following four perspectives to develop the proposed approach: \begin{enumerate} • Information perspective: Because of the serial dependence in time-series data, past lagged variables of ${\mathbf x}_t$ are often helpful in predicting $y_{t+h}$. For example, the lagged variables are used in a factor-augmented VAR model in Bernanke2004 for studying the effect of monetary policy on the economy. Yet there is no unified approach available under the existing methods to incorporating lagged information. Different users make use of the lagged variables in different ways. Consider the diffusion approach. One can augment lagged values of the predictors to form an extended predictor vector before applying PCA to extract common factors. The extracted factors are then linear combinations of ${\mathbf x}_t$ and its lagged variables. Another possibility is to use lagged values of the extracted factors obtained from the original predictors. This procedure assumes {\em a priori} that the linear combinations of predictors carry over lags. Consequently, the lagged information is exploited differently and there are no guidelines available to best use the lagged variables. In addition, there is no simple way to select the number of past lagged variables needed in an application. • Adaptive perspective: Empirical applications are often interested in both short-term and long-term predictions so that various choices of $h$ in Model ((ref)) are employed. Yet the same common factors are used for all choices of $h$ because an unsupervised PCA does not consider the target variable in extracting common factors. It is possible in applications that useful predictors for $y_{t+h}$ may not be good ones for $y_{t+h'}$ under the setting in Model ((ref)), where $h\neq h'$. • Supervision perspective: In machine learning, applying principal component analysis to extract the latent factors from ${\mathbf x}_t$ is an unsupervised procedure. It does not require any information of the target variable $y_{t+h}$. This might become a drawback in prediction as the target variable is of main interest. In addition, PCA is not scale-invariant, and the extracted factors can change easily when scales of the components of ${\mathbf x}_t$ change. This may lead to inferior prediction if care is not exercised in applying diffusion index models. For this reason, the components of ${\mathbf x}_t$ are often standardized before applying PCA. • Accuracy in factor extraction: The traditional PCA method may not provide accurate estimation of the common factors, especially when some common factors are strong and the remaining ones are weak; see, for example, the explanation in huang2022scaled and bai2021approximate. For time series data, common factors that have stronger serial dependence tend to dominate those with weaker serial dependence as strong serial dependence often results in higher variances, which in turn leads to larger eigenvalues in the covariance matrix of ${\mathbf x}_t$. \end{enumerate} In view of the above discussions, we propose a new diffusion-index model in this paper for dynamic forecasting by developing a supervised dynamic PCA (sdPCA) method. First, our proposed sdPCA takes explicitly the lagged variables of predictors into account. Each predictor is allowed to have its own number of past lagged variables. In fact, the number of lagged variables for each predictor may be selected by an information criterion function such as the Akaike Information Criterion. Second, our new diffusion-index model can also employ lagged variables of the extracted common factors. Therefore, the proposed approach explores the dynamic dependence in two ways. The extracted factors by the proposed sdPCA contain lagged variables of the predictors and the new diffusion-index model also employs lagged variables of the common factors. Consequently, our new approach is different from the settings in Stock2002b and in bai2006, where, as discussed earlier, the lagged dependence is not fully explored. Third, our sdPCA incorporates information of the target variable $y_{t+h}$ in constructing the predictor vector for PCA so that the scale of each component of the predictor vector is properly adjusted. Consequently, the proposed new approach does not depend on the scales of observed predictor ${\mathbf x}_t$. This property is similar to that of the sPCA of huang2022scaled. Fourth, the sdPCA is adaptive in extracting common factors so that the common factors used to predict $y_{t+h}$ and $y_{t+h'}$ may be different. This can enhance the predictive power of the proposed new diffusion-index model. There are two steps in the proposed sdPCA that produces useful factors to serve as predictors. First, sdPCA runs predictive linear regressions of the target variable on each observed predictor and its past lagged variables with the number of lags being selected by an information criterion. In this way, useful lagged variables of each predictor are properly addressed. The fitted values of the individual predictive linear regressions form a new high-dimensional predictor vector on which PCA is carried out to extract common factors. Since the predictive linear regression automatically adjusts the scale of each predictor, the scaling effect of observed predictors on PCA is also properly addressed. As a matter of fact, one can think of the scale of each predictor being adjusted based on its predictive power to the target variable. This adjustment may depend on $h$, the number of step ahead prediction of interest. Second, sdPCA allows users to add lagged variables of extracted factors to the contemporaneous ones in the linear prediction model. This enables us to explore farther the dynamic information of the data to improve the forecasts. A similar technique is discussed in Stock2002a under different settings. In short, these two steps of the proposed sdPCA not only capture the dynamic information of the data in an additive manner, but also adjust automatically the scaling effect of the data. More importantly, the steps enable users to explore fully the linear dependence of the data to improve the accuracy of prediction. We show that, under some general conditions, the proposed sdPCA can outperform the traditional diffusion-index approach using unsupervised PCA, both in theory and in simulation. Asymptotic properties of the proposed sdPCA are established under the modern setting that the number of predictors $N$ and the sample size $T$ diverge to infinity. We also compare the proposed method with some commonly used methods, and derive the conditions under which, the proposed method can outperform the existing ones theoretically. In addition, to embrace the modern development of machine learning techniques, we also apply the well-known least absolute shrinkage and selection operator (Lasso) approach to select the most relevant factors that have predictive power for the target variable. Theoretically, we also establish the consistency of the Lasso estimators under some identification conditions. We illustrate and assess the performance of the proposed sdPCA and the new diffusion-index forecasting method with an application to macroeconomic index forecasting. We forecast the U.S. industrial production (IP) growth, change in the unemployment rate (UNRATE), the consumer price index: all (CPI-All), the S&P 500 index volatility change (Volatility Change), and the S&P 500 index return using 123 macroeconomic variables from FRED-MD, as that in huang2022scaled, mccracken2016fred, Stock2002b,Stock2002a, among others. Similarly to those in huang2022scaled, the sdPCA loadings have re-assigned the weights to the predictors and a smaller subset of the macro variables tend to have more predictive power compared with the unsupervised PCA loadings. Furthermore, the proposed sdPCA together with the Lasso procedure produces comparable or even better forecasting results compared with some commonly used factor-based forecasting methods, such as the traditional PCA method, the sPCA in huang2022scaled, and the diffusion-index model in Stock2002b,Stock2002a. The scaled PCA (sPCA) of huang2022scaled is the closest method to ours. Both methods reconstruct a new high-dimensional prediction vector before extracting common factors using PCA. However, our method differs from the sPCA in three ways. First, our proposed diffusion-index forecasts may include lagged variables of common factors as predictors while sPCA only uses the contemporaneous one. Thus, the predictive model in huang2022scaled is similar to that in Stock2002b whereas ours generalizes the traditional method to a dynamic fashion. As a matter of fact, sPCA is a special case of the sdPCA if one does not include any lagged variables of each predictor in constructing the high-dimensional prediction vector, nor use any lagged variables of the extracted factors. Second, instead of running simple linear regression of the target variable on each predictor alone, we run a time series regression by including lagged variables of the predictors with the number of lags being selected by an information criterion. In other words, the proposed sdPCA explores the dynamic dependence of the target variable on each observed predictor in constructing the new high-dimensional predictor vector for PCA. Third, the proposed prediction model uses the Lasso method to select the factors that have significant predictive power in prediction. This provides a data-driven approach to identifying relevant predictors. There are other related works in the literature concerning supervised learning or building connections between predictors and the target variable in prediction; see the references in huang2022scaled. We briefly discuss some differences between the proposed sdPCA approach and some related methods. bai2008forecasting first applied a screening method to select a subset of predictors that are tested to have relatively more predictive power to the target under either soft- or hard-thresholding rules, and the PCA is conducted on the selected subset to extract common factors for use in the diffusion-index forecasting. Our proposed method, on the other hand, assigns different weights to each predictor and its lagged variables according to their predictive power without using any thresholding rule. An incomplete list of other related works that share similar insights to ours includes the partial least squares (PLS) regression in wold1966estimation and kelly2015three with financial applications in kelly2013market, huang2015investor, and light2017aggregation, among others. The comparison of the aforementioned methods and the sPCA has been extensively studied in huang2022scaled, and the sPCA is shown to have advantages in forecasting. Therefore, we only compare our proposed method with the sPCA and some commonly used factor-based linear forecasting methods in this paper. The rest of the paper is organized as follows. Section (ref) introduces the sdPCA method and the new diffusion-index model, and presents their asymptotic properties and some comparison results. Section (ref) studies the finite-sample performance of the proposed approach via simulation, and (ref) illustrates the proposed procedure with an empirical application. Section (ref) concludes. All the proofs and derivations for the asymptotic results are relegated to an online Appendix. {\bf Notation:} We use the following notation. For a $p\times 1$ vector ${\mathbf u}=(u_1,..., u_p)'$, $\|{\mathbf u}\|_1=\sum_{i=1}^p|u_i|$ is the $\ell_1$-norm and $\|{\mathbf u}\|_\infty=\max_{1\leq i\leq p}|u_i|$ is the $\ell_\infty$-norm. ${\mathbf I}_p$ denotes the $p\times p$ identity matrix. For a matrix ${\mathbf H}$, its Frobenius norm is $\|{\mathbf H}\|=[\mathrm{trace}({\mathbf H}'{\mathbf H})]^{1/2}$ and its operator norm is $\|{\mathbf H} \|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$, where $\lambda_{\max} (\cdot) $ denotes the largest eigenvalue of a matrix, and $\|{\mathbf H}\|_{\min}$ is the square root of the minimum non-zero eigenvalue of ${\mathbf H}{\mathbf H}'$. The superscript ${'}$ denotes the transpose of a vector or matrix. We also use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$. \section{Methodology} \subsection{Model Setup} Let ${\mathbf x}_t=(x_{1,t},...,x_{N,t})'$ be an $N$-dimensional observable time series, for $t=1,...,T$, and $y_{t+h}$ be the target variable of interest, where $h \geq 1$. The goal is to predict $y_{t+h}$ using ${\mathbf x}_t$ and its past information. Based on the discussion in the {\it Introduction} section, we assume the data are centered and consider the following model: \begin{equation} \begin{array}{c} {\mathbf x}_t={\mathbf B}{\mathbf f}_t+{\mathbf u}_t, \\ y_{t+h}=\boldsymbol{\beta}(L)'{\mathbf f}_t+{\varepsilon}_{t+h}, \end{array} \end{equation} where $\boldsymbol{\beta}(L)=\boldsymbol{\beta}_0+\boldsymbol{\beta}_1L+...+\boldsymbol{\beta}_{q-1}L^{q-1}$ with $L$ being the backshift (or lag) operator such that $L{\mathbf f}_t = {\mathbf f}_{t-1}$, and $q\geq 1$ is the number of lagged variables of factors used in predictions. When $q=1$, Model ((ref)) reduces to the diffusion-index forecasting equation in ((ref)). The model in ((ref)) is similar to that in Stock2002b and bai2006, but the mechanisms to produce predictions of the models are fundamentally different. In Stock2002b, a dynamic factor model is considered for a large panel of time series and some dynamic factors are used as predictors, which are assumed to explain most of the variability of the panel, measured by the covariance matrix of the panel. Their factors not necessarily contain information of the past lagged values of ${\mathbf x}_t$. In Model ((ref)), the observed time series ${\mathbf x}_t$ admits a factor structure with static factor processes, and these static factors and their past lagged variables are used as predictors in the prediction equation. For similar reasons, Model ((ref)) is also different from the setting in bai2006 and huang2022scaled. On the other hand, there are some similar insights between Model ((ref)) and Model (2.3) in Stock2002a because both models attempt to explore the dynamic dependence of the data in the prediction. But the proposed model includes past lagged variables of the common factors whereas the one in Stock2002a employs the past lagged variables of $y_t$ in the prediction. A recursive expansion of $y_t$ in Model (2.3) of Stock2002a also leads to using past lagged variables of common factors as predictors albeit with some constraints in the coefficient parameters. For simplicity, we only consider Model ((ref)) in this paper and investigate its predictive ability. Under some identification conditions, it is natural to apply PCA to extract common factors from ${\mathbf x}_t$ and use the extracted factors and their lagged variables as predictors in predicting $y_{t+h}$. Alternatively, one may perform PCA on a stacked vector $({\mathbf x}_t',{\mathbf x}_{t-1}',...,{\mathbf x}_{t-q+1}')'$ and apply the extracted factors directly in forecasting. Since the factors so obtained already include the lagged variables of ${\mathbf x}_t$ so that no lagged variables of the factors are used. However, there are some drawbacks of this approach. First, the application of PCA in the first step does not learn any information from the target variable $y_{t+h}$ and, therefore, the factors extracted directly from the stacked vector may not have the best predictability for $y_{t+h}$. Second, the consistency of the PCA procedure is usually shown under the assumption of strong factors in the literature; see, for example, bai2002determining and fan2013large. In practice, there is no guarantee that all the factors are strong because the noise effect can be prominent when adding more variables to ${\mathbf x}_t$. As a result, the factors extracted from all the components of ${\mathbf x}_t$ (or stacked vector) may not have better predictability than those from a subset of the panel. See boivin2006more for further information. To overcome these drawbacks, we introduce a new supervised dynamic PCA method to extract the common factors and to employ relatively informative factors in a three-stage procedure to explore the dynamic dependence in the data. Our goal is to extract the factors that have more predictive power by learning from the target variable. In other words, we need to estimate the factors ${\mathbf f}_t$ using a new approach for predictions. Details are discussed next. \subsection{Estimation Procedure} In this section, we propose a new PCA useful for dynamic forecasting. The procedure consists of the following three steps: \begin{enumerate} • For $i=1,...,N$, estimate the slope parameters by regressing the target variable on the $i$-th predictor and its past lagged variables: \begin{equation} y_{t+h}\approx \widehat\mu_i+\widehat\gamma_{i,0}x_{i,t}+\widehat\gamma_{i,1}x_{i,t-1}+...+\widehat\gamma_{i,q_i-1}x_{i,t-q_i+1},\,\,t=q_i,...,T-h, \end{equation} where $q_i$ is selected by an information criterion such as the Akaike Information Criterion (AIC). Let $q=\max_{1\leq i \leq N}\{q_i\}$. • For $t=q,q+1,...,T$, let $\widehat{\mathbf x}_t=(\widehat x_{1,t},...,\widehat x_{N,t})'$ with $\widehat x_{i,t}=\widehat\gamma_{i,0}x_{i,t}+\widehat\gamma_{i,1}x_{i,t-1}+...+\widehat\gamma_{i,q_i-1}x_{i,t-q_i+1}$. Apply PCA to $\widehat{\mathbf x}_t$ and obtain the estimated factors $\widehat{\mathbf g}_t$, which contains the information of the original factors ${\mathbf f}_t$ as well as its lagged variables relevant to $y_{t+h}$. One can think of $\widehat{\mathbf g}_t$ as an estimator for $({\mathbf f}_t',...,{\mathbf f}_{t-q+1}')'$ in Model ((ref)) under some proper conditions. • For $t=q,q+1,..., T-h$, apply the Ordinary Least-Squares method to the linear regression of the target variable on $\widehat {\mathbf g}_t$ and obtain the estimated coefficients: \begin{equation} y_{t+h}\approx \widehat\alpha+ \widehat\boldsymbol{\beta}'\widehat{\mathbf g}_t,\,\, t=q,...,T-h. \end{equation} Finally, the prediction of $y_{T+h}$ is given by $\widehat y_{T+h}=\widehat\alpha+\widehat\boldsymbol{\beta}'\widehat{\mathbf g}_T$. \end{enumerate} Note that the intercepts $\widehat\mu_i$ and $\widehat\alpha$ in Step 1 and Step 3 above will be removed if the data are assumed to be centered as that in Model ((ref)). Some remarks are as follows. First, the idea of learning from the target variable using linear regression is similar to that in huang2022scaled. However, we focus on dynamic forecasting with lagged variables in the regression, while the sPCA of the aforementioned paper only uses contemporaneous information. Second, unlike the traditional method, the number of lagged variables of each observed predictor used in prediction is selected during the first step. Once the additive panel is formed according to the regression results in Step 1, the extracted factors in Step 2 will automatically contain the dynamic information and can be used as predictors in Step 3, which captures the dynamic dependence in time series forecasting. Third, we only use $\widehat{\mathbf g}_t$ in Step 3, but we can also employ some of its lagged variables if necessary. Furthermore, as will be seen in the next section, we apply Lasso regularization to obtain $\widehat\boldsymbol{\beta}$ so that adding additional lagged variables of $\widehat{\mathbf g}_t$ does not cause any difficulties to the proposed procedure. Fourth, when the number of predictors $N$ and the sample size $T$ are large, one can use a sufficiently large value $q$ in lieu of $q_i$ in Step 1 to simplify the computation. \subsection{When Diffusion-Index Forecast Meets Sparsity} As discussed before, even though we have constructed an $N$-dimensional vector of predictors in Step 2, it is possible that only a subset of the extracted common factors has predictive power for the target variable. Similar insights are also mentioned in huang2022scaled, but a statistical approach to automatically select the significant factors is still not available therein. In this section, we propose to use penalized regression methods to select the factors that have predictive power for the target variable. For simplicity, we assume the intercepts are zero and only introduce the Lasso approach since other penalized methods can be similarly established. Let $\widehat\boldsymbol{\beta}_{lasso}$ be the Lasso solution that solves the following optimization problem: \begin{equation} \widehat\boldsymbol{\beta}_{lasso}=\arg\min_{\boldsymbol{\beta}\in R^{rq}}\left\{\frac{1}{T}\sum_{t=q}^{T-h}\|y_{t+h}-\boldsymbol{\beta}'\widehat{\mathbf g}_t\|^2+\lambda_T\|\boldsymbol{\beta}\|_1\right\}, \end{equation} where $\widehat{\mathbf g}_t$ is the extracted factor process from Step 2 of the proposed procedure in Section (ref), and $\lambda_T>0$ is a penalty parameter to be determined later. The above optimization is a convex one and can be solved by many existing algorithms. See hastie2009elements for detailed illustrations. The basic idea of using regularization estimation, such as the Lasso regression, is that possibly only a subset of $\widehat{\mathbf g}_t$ has predictive power for $y_{t+h}$. This is particularly so if $N$ is large. In this case, $\boldsymbol{\beta}$ is a sparse vector with only a few non-zero elements. The penalty parameter $\lambda_T$ is used to control the number of significant factors to be used in the predictions. It is natural to ask why the proposed sdPCA and the new diffusion-index forecasts can outperform the traditional PCA and the scaled PCA in huang2022scaled in prediction. There are several reasons that can answer this question. First, consider the information aspect. Our proposed sdPCA generalizes the scaled PCA by including more relevant lagged variables of observed predictors. As such, the scaled PCA is a special case of the proposed method with $q_i=1$ at Step 1, for $i=1, \ldots, N$. Consequently, the proposed sdPCA captures more dynamic dependence information in the data and, hence, can improve the forecasting performance. Second, the proposed Lasso regression can further screen out irrelevant predictors, which in turn can improve parameter estimation and avoid the difficulty of over-parameterization, especially when $N$ is large. \begin{remark} The Lasso procedure could be very helpful in practice because we do not know how many factors should be included as predictors in Step 3 of Section (ref). Note that $\widehat{\mathbf g}_t$ in Step 3 can be treated as a proxy for $({\mathbf f}_t',...,{\mathbf f}_{t-q+1}')'$, which contains the dynamic information learned by regressing the target variable on lagged predictors in Step-1 of the proposed procedure. In fact, we may also add the lagged variables of $\widehat{\mathbf g}_t$ and adopt Lasso to select the significant factors in prediction. This expanded approach can collect additional dynamic information that may be missed in Step 1. We do not explore any further this approach to save space. \end{remark} \begin{remark} It is well-known that the Lasso estimate $\widehat\boldsymbol{\beta}_{lasso}$ is a biased estimate of the true coefficient $\boldsymbol{\beta}$. To correct the bias, one can re-run the linear regression using only those predictors selected by the Lasso procedure. This is referred to as a post-selection inference in the literature; see, for instance, belloni2013. \end{remark} \subsection{Selection of the Lag Parameters} We use information criteria in Step 1 of the proposed procedure to select the number $q_i$ of lagged variables for the predictor $x_{it}$; see Equation ((ref)). In applying any information criterion, there is a need to select the maximum order allowed. This maximum order is unknown in practice, but there are statistical methods available to guide the choice. We mention two possibilities in this section. First, from a time series analysis point of view, the maximum order provides an approximation to the {\em true} lagged linear dependence between the target $y_{t+h}$ and $\{x_{it},x_{i,t-1},x_{i,t-2},\ldots\}$. When the sample size $T$ is large, one can improve the accuracy in approximation by increasing the maximum order allowed. Therefore, a common practice in the literature is to use $q_{max} \approx \log(T)$. In this way, the information criterion used selects $q_i \in [1,\ldots, q_{max}]$. This approach can also be applied to the selection of $q$ in the diffusion-index forecasting model of Step 3. A second method is cross-validation, which is particularly useful when $N$ and $T$ are both large. We partition the data into two sub-samples, say $\{{\mathbf x}_1,...,{\mathbf x}_{T_1-h},y_{h},...,y_{T_1}\}$ and $\{{\mathbf x}_{T_1-h+1},...,{\mathbf x}_{T-h},y_{T_1+1},...,y_{T}\}$ for some $T_1<T$. For a small integer $q_{max} \geq 1$ and each $1\leq q\leq q_{max}$, we perform the proposed estimation of Section (ref) to the first sub-sample with $q_i = q$, for all $i$, and predict $y_{T_1+h}$ to compute the associated forecasting error. We then move the data $\{{\mathbf x}_{T_1-h+1},y_{T_1+1}\}$ from the second sub-sample to the first one and repeat the above procedure to calculate another forecasting error. This estimation-forecasting exercise is repeated until we obtain the forecast error of $y_T$. Then $\widehat q$ is chosen as the one that produces the smallest out-of-sample forecasting errors, where we may adopt the mean-squared forecasting error (RMSE) defined in ((ref)) in Section (ref) below. \begin{remark} In the empirical study of macroeconomic forecasting, we adopt the second approach mentioned above because the number of predictors is large, and many macroeconomic indicators only depend on a small number of the most recent lagged variables. The AIC or BIC criterion can be treated as another option if one wants to further improve the forecasting accuracy of some variables of interest, though information criteria may not always outperform cross-validation in lag selection. \end{remark} \subsection{Selection of the Number of Factors} The above analysis depends on a known number of factors $r$, which is unknown in practice. As discussed in bai2021approximate, if we want to estimate the number of factors with $\nu>0$, where $\nu$ is a strength parameter defined in Assumption (ref) below, the criteria in bai2002determining remain useful. In addition, there are other estimation methods such as the criteria of onatski2010determining and ahn2013eigenvalue which separate the bounded eigenvalues from diverging ones of the covariance matrix. For time series factor models with weak factors, lam2012factor proposed a multi-step eigenvalue-ratio method to estimate the number of factors with different strengths. This method remains valid under our framework and it is especially useful when the factors have different strengths of weaknesses. Therefore, we may apply those existing methods to estimate the number of factors. Note that it might be helpful to include more factors as predictors and let the Lasso procedure select the factors that have more predictive power. Limited simulation results suggest that the Lasso method can accurately identify the number of significant factors. We mention that there is an extreme case that the loading matrix associated with ${\mathbf g}_t$ in Step 2 may not be of full rank when $\boldsymbol{\beta}_i=\boldsymbol{\beta}_j$ for some $0\leq i,j\leq q-1$. Then the number of factors identified by the aforementioned methods would be fewer than the true one. But this is not an issue because we can treat ${\mathbf f}_{i,t}+{\mathbf f}_{j,t}$ as a new factor with components sharing a common regression coefficient in the linear forecasting step via a regression method. Simulation results in Section (ref) suggest that the proposed method still works well in out-of-sample prediction. \begin{remark} It is natural that the information criterion in bai2002determining should be refined to cover the case of weak factors. However, a valid criterion depends on an accurate estimation of the strength parameter $\nu$ in Assumption (ref) below, which is difficult to obtain in practice. This difficulty can be addressed by the Lasso procedure introduced in Section (ref). Specifically, we adopt a data-driven procedure by including more factors as predictors in our empirical studies. The Lasso procedure can then select the factors that play important roles in forecasting the target variable. \end{remark} \subsection{Assumptions} In this section, we introduce the assumptions needed to derive theoretical results of the asymptotic forecasting performance of the proposed method. Most assumptions below are commonly used in the PCA or approximate-factor modeling literature, and the derivation of the consistency of the LASSO estimates needs some slightly stronger assumptions. We use $c$ or $C$ to denote a generic positive constant the value of which may change at different places. \begin{assumption} The process $\{{\mathbf f}_t\}$ is $\alpha$-mixing with the mixing coefficients satisfying the condition $\alpha_N(k)<\exp(-k)$, where $\alpha_N(k)$ is defined as \begin{equation} \alpha_N(k)=\sup_{i}\sup_{A\in\mathcal{F}_{-\infty}^i,B\in \mathcal{F}_{i+k}^\infty}|P(A\cap B)-P(A)P(B)|, \end{equation} where $\mathcal{F}_i^j$ is the $\sigma$-field generated by $\{{\mathbf f}_t:i\leq t\leq j\}$. \end{assumption} \begin{assumption} For ${\mathbf g}_t=({\mathbf f}_t',{\mathbf f}_{t-1}',...,{\mathbf f}_{t-q+1}')'$ with a fixed $q \geq 1$, $\sup_tE\|{\mathbf g}_t\|^4\leq C$ and $\frac{1}{T}\sum_{t=q}^{T-h}{\mathbf g}_t{\mathbf g}_t'\rightarrow_p \boldsymbol{\Sigma}_g$, which is a $qr\times qr$ positive-definite matrix. \end{assumption} Note that the $i$-th row of loading ${\mathbf b}_i$ can be either a random or a fixed constant vector. Either way, we may define ${\mathbf b}_i={\bf 0}$ symbolically if $x_{i,t}$ does not depend on the common factors. Thus, ${\mathbf b}_i\neq {\bf 0}$ implies that $x_{i,t}$ depends on the common factors with some positive probability. We define $\mathcal{I}_b=\{i:{\mathbf b}_i\neq {\bf 0}, i=1,\ldots,N\}$, which consists of the indexes for which the corresponding predictors depend on the common factors with positive probabilities. The following assumption is related to the strength of the factor loading. \begin{assumption} $\sup_{1\leq i\leq N} E\|{\mathbf b}_i\|^4\leq C$ holds. The cardinality of the set $\mathcal{I}_b$ satisfies $\text{Card}(\mathcal{I}_b)\asymp N^\nu$, for some $0 < \nu \leq 1$, and $\frac{1}{N^\nu} \sum_{i=1}^N{\mathbf b}_i{\mathbf b}_i'\rightarrow_p\boldsymbol{\Sigma}_B$, which is an $r\times r$ positive definite matrix. \end{assumption} If $\nu = 1$ in Assumption 3, all factors are strong as that in bai2002determining and fan2013large. Similar to Assumption A2 in bai2021approximate, we exclude the case of $\nu=0$ because the factors and the idiosyncratic terms are indistinguishable in such a situation. \begin{assumption} The idiosyncratic term $u_{i,t}=\sigma_i e_{i,t}$ for some $c\leq \sigma_i\leq C$, where $e_{i,t}$ is independent and identically distributed over $i$ and $t$ with the eighth moment bounded. \end{assumption} \begin{assumption} $\{{\mathbf b}_i,1\leq i\leq N\}$, $\{{\mathbf f}_t,1\leq t\leq T\}$, and $\{{\mathbf u}_{i,t},1\leq i\leq N,1\leq t\leq T\}$ are mutually independent with each other. \end{assumption} \begin{assumption} $\{\mbox{\boldmath$\varepsilon$}_t\}$ is independent with the three sets of variables in Assumption (ref) and it is a martingale-difference sequence such that $E(\mbox{\boldmath$\varepsilon$}_{t+h}|\mathcal{F}_t)=0$ for any integer $h>0$, where $\mathcal{F}_t$ is the $\sigma$-field generated by $\{{\mathbf g}_t,{\mathbf u}_t,{\mathbf g}_{t-1},{\mathbf u}_{t-1},...\}$. Furthermore, $\sup_tE(\mbox{\boldmath$\varepsilon$}_t^4)\leq C$. \end{assumption} Assumption (ref) is standard to characterize the dynamic dependence of the factor processes. See, for example, gao2019banded. It is used to control the magnitude of joint partial sums in the derivations as well as the consistency of the Lasso estimators in Section (ref). Assumptions (ref)-(ref) are similar to those in huang2022scaled, and they also imply that the results in Assumptions A1-A3 of bai2021approximate hold, except for the assumption of distinct eigenvalues in A2(iii) therein. In fact, they are adequate for proving the consistency of the estimators and deriving the asymptotic forecasting performance. The illustrations of all the assumptions are stated in huang2022scaled or bai2021approximate, and we omit the details to save space. The assumption of distinct eigenvalues in A2(iii) of bai2021approximate is only used to show the limiting distributions of $\widetilde{\mathbf F}'{\mathbf F}/T$, where $\widetilde{\mathbf F}$ is the PC estimator for ${\mathbf F}$ therein, and we will make similar assumptions below in order to show the consistency of the Lasso estimators in Section (ref). It is worth mentioning that the independence assumption in Assumption (ref) is only made to simplify the theoretical derivations. It can be relaxed to a weaker assumption such as those in bai2021approximate, and those bounds and inequalities in Assumptions A1 and A3 therein can be verified if we impose some mixing condition or weak dependence assumption on the idiosyncratic vector both cross-sectionally over space and dynamically over time. For the consistency of the Lasso estimators, we also need the following assumptions. Let ${\mathbf B}_\gamma'=(\boldsymbol{\gamma}_1\otimes{\mathbf b}_1,...,\boldsymbol{\gamma}_N\otimes{\mathbf b}_N)$, where $\boldsymbol{\gamma}_i=({\mathbf I}_q\otimes {\mathbf b}_i')\boldsymbol{\beta}$, for $1\leq i\leq N$. \begin{assumption} ${\mathbf G}'{\mathbf G}/T={\mathbf I}_{rq}$ and ${\mathbf B}_\gamma'{\mathbf B}_\gamma$ is a diagonal matrix with distinct eigenvalues, where ${\mathbf G}=({\mathbf g}_q,...,{\mathbf g}_{T-h})'$ consists of ${\mathbf g}_t$'s defined in Assumption (ref) as its row vectors. \end{assumption} \begin{assumption} For any $1\leq i\leq r$, and $1\leq t\leq T$, $P(|f_{i,t}|>x)\leq C_0\exp(-C_1x)$ and $P(|{\varepsilon}_{t}|>x)\leq C_0\exp(-C_1x)$, where $C_0>0$ and $C_1>0$ are constants. \end{assumption} Assumption (ref) is an identification condition that guarantees the uniqueness of the estimated factors. See the illustration in bai2013principal for details. Assumption (ref) controls the tails of the factors and the random errors in the forecasting model, and is essentially a sub-exponential assumption. This assumption is stronger than the moment conditions above, but they are adequate for establishing the consistency of the Lasso estimators with the theory developed in merlevede2011bernstein. \subsection{Asymptotic Forecasting Performance} In this subsection, we present some theoretical properties of the proposed estimators and compare the asymptotic forecasting performance of the proposed method with some existing ones. Letting $\widehat{\mathbf G}_{\text{sdPCA}}$ be the estimated factors using the proposed method, we have the following consistency result. \begin{theorem} Suppose that Assumptions (ref)-(ref) hold. If $N^{1-\nu}/T^2\rightarrow 0$, there exits an invertible rotation matrix ${\mathbf H}_{\text{sdPCA}}$ such that the estimated factors satisfy \[\frac{1}{\sqrt{T}}\|\widehat{\mathbf G}_{\text{sdPCA}}-{\mathbf G}{\mathbf H}_{\text{sdPCA}}'\|=O_p(N^{-\nu/2}+T^{-1}+\frac{N^{1-\nu}}{T^2}).\] \end{theorem} The following proposition provides the condition under which the traditional PCA estimators are also consistent and the condition under which the traditional PCA does not produce consistent factor estimates. \begin{proposition} Suppose that Assumptions (ref)-(ref) hold.\\ (i) If $\lim\inf N^{1-\nu}/T \geq C$, for some positive constant $C>0$, then, for any invertible matrix ${\mathbf H}_{PCA}$, the following result holds, \[\frac{1}{\sqrt{T}}\|\widehat{\mathbf F}_{\text{PCA}}-{\mathbf F}{\mathbf H}'_{PCA}\|\geq C>0,\] implying that $\widehat{\mathbf F}_{\text{PCA}}$ is not a consistent estimator under such conditions.\\ (ii) If $N^{1-\nu}/T\rightarrow 0$, there exists a rotation matrix ${\mathbf H}_{PCA}$ such that, the tradition PCA estimator, denoted by $\widehat{\mathbf F}_{PCA}$, satisfies \[\frac{1}{\sqrt{T}}\|\widehat{\mathbf F}_{\text{PCA}}-{\mathbf F}{\mathbf H}'_{PCA}\|=O_p(N^{-\nu/2}+\frac{N^{1-\nu}}{T}).\] \end{proposition} Some remarks on the results in Theorem (ref) and Proposition (ref) are in order. First, The conditions in Theorem (ref) and in Proposition (ref) indicate that the requirement for the consistency of the factor estimation using the proposed method is relatively weak. For example, if $N^{1-\nu}/T^2\rightarrow 0$ but $\lim\inf N^{1-\nu}/T>C$ for some constant $C>0$, the extracted factors using the proposed method are consistent while those by the traditional method are not. Second, even if $N^{1-\nu}/T^2\rightarrow 0$, under which both methods produce consistent factors, the factors obtained by the proposed method can still produce more accurate predictions as shown in Theorem (ref) below, where the smaller term in ((ref)) converges to zero faster than that in ((ref)). Finally, if all the factors are strong ones, that is, $\nu=1$, the convergence rates of the factors extracted by two methods are the same. Next, we compare the asymptotic forecasting performance of several existing linear forecasting methods using factors. Define the mean-squares-forecast error (MSFE) as \begin{equation} MSFE={\frac{1}{T}\sum_{t=q}^{T-h}(y_{t+h}-\widehat y_{t+h})^2}, \end{equation} where $\widehat y_{t+h}=\widehat\boldsymbol{\beta}'\widehat{\mathbf g}_t$. Let $w_{1,NT}=N^{-\nu/2}+T^{-1}+\frac{N^{1-\nu}}{T^2}$, $w_{2,NT}=N^{-\nu/2}+\frac{N^{1-\nu}}{T}$, ${\mathbf y}=(y_{h+1},...,y_{T})'$ and $\widehat{\mathbf y}=(\widehat y_{h+1},...,\widehat y_{T})'$. In the following error analysis, we assume the true forecasting model is the one in ((ref)), and denote the MSFE produced by the proposed method, the method of huang2015investor, and the one in Stock2002b by $\text{MSFE}_{sdPCA}$, $\text{MSFE}_{sPCA}$, and $\text{MSFE}_{SW}$, respectively. The MSFE of the forecasting method that stacks all extracted factors by the traditional PCA without rescaling is denoted by $\text{MSFE}_{PCA}$. Note that both $\text{MSFE}_{sdPCA}$ and $\text{MSFE}_{PCA}$ make use of lagged factors in prediction while the other two methods only use contemporaneous factors. We have the following theorem concerning the forecasting performance of different methods. \begin{theorem} Suppose that Assumptions (ref)-(ref) hold.\\ (i) If $N^{1-\nu}/T^2\rightarrow 0$, \begin{equation} MSFE_{sdPCA}=\frac{1}{{T}}\|{\mathbf y}-\widehat{\mathbf y}\|^2= \frac{1}{{T}}\|({\mathbf I}-\frac{1}{T}{\mathbf G}{\mathbf G}')\boldmath$\varepsilon$\|^2+O_p(w_{1,NT}^2). \end{equation} (ii) If $N^{1-\nu}/T\rightarrow 0$, \begin{equation} \text{MSFE}_{PCA}=\frac{1}{{T}}\|{\mathbf y}-\widehat{\mathbf y}\|^2= \frac{1}{{T}}\|({\mathbf I}-\frac{1}{T}{\mathbf G}{\mathbf G}')\mbox{\boldmath$\varepsilon$}\|^2+O_p(w_{2,NT}^2). \end{equation} (iii) If $\lim\inf N^{1-\nu}/T\geq c>0$ and $N^{1-\nu}/T^2=o(1)$, then \[\text{MSFE}_{PCA}-\text{MSFE}_{sdPCA}\geq C>0,\] implying that the proposed method outperforms the traditional one in theory.\\ (iv) If $N^{1-\nu}/T\rightarrow 0$, we have \[\text{MSFE}_{SW}-\text{MSFE}_{sdPCA}\geq C>0,\] with a strictly positive probability, where, as defined earlier, $\text{MSFE}_{SW}$ denotes the mean-squares of forecasting errors using the diffusion-index model in Stock2002b. \\ (v) If $N^{1-\nu}/T^2\rightarrow 0$, then \[\text{MSFE}_{sPCA}-\text{MSFE}_{sdPCA}\geq C>0,\] with a strictly positive probability. \end{theorem} Some remarks on the results in Theorem (ref) are given below. First, the results in Theorem (ref)(iii) show that the proposed forecasting method outperforms the one using the traditional PCA if $\lim\inf N^{1-\nu}/T\geq c>0$ and $N^{1-\nu}/T^2=o(1)$. Second, even if $N^{1-\nu}/T\rightarrow 0$, the smaller term in ((ref)) converges to zero faster than that in ((ref)), implying that the proposed forecasting method still has a certain probability to make more accurate predictions. Third, Theorem (ref)(iv)-(v) show that our forecasting method can produce smaller errors than the ones in Stock2002b and huang2022scaled, where the latter two methods only include the contemporaneous factors and overlook the relevant information in the lagged variables. Next, we present the asymptotic behavior of the Lasso estimator in the forecasting step. \begin{theorem} Suppose that Assumptions (ref)-(ref) hold. If the sparsity is $0<s^*<rq$ and the penalty parameter $\lambda_T\geq\max\{M\sqrt{\frac{rq}{T}},Mw_{1,NT}\}$ for some $M>0$, then, with probability tending to one, we have \[\|\widehat\boldsymbol{\beta}_{lasso}-\boldsymbol{\beta}\|\leq C\lambda_T,\] where $C$ may depend on $s^*$. The Lasso estimate is consistent under the conditions in Theorem (ref) if we choose $\lambda_T=C^*\max(\sqrt{\frac{rq}{T}},w_{1,NT})$ for some appropriate constant $C^*>0$. \end{theorem} The result in Theorem (ref) is a classical one in the Lasso literature; see, for example, buhlmann2011statistics. When the penalty parameter is properly chosen, Theorem (ref) indicates that we can correctly recover the non-zero elements in the linear regression asymptotically. Simulation results in Section (ref) suggest that the Lasso procedure works sufficiently well in finite samples. In addition, it can also improve the forecasting performance as shown in the empirical example in Section (ref). \section{Simulation Studies} In this section, we use Monte-Carlo experiments to compare the forecasting performance of the proposed method with some existing ones. The data-generating process (DGP) used is given below. Consider a two-factor model in the experiment, that is, the number of factors $r=2$ in Model ((ref)). The factor process ${\mathbf f}_t$ is independently generated from normal distributions with zero mean and identity covariance, that is, ${\mathbf f}_t\sim N({\bf 0},{\bf {\mathbf I}_2})$. The idiosyncratic terms $u_{i,t}$ are independently and normally distributed with zero mean and unit variance, i.e., $u_{i,t}\sim N(0,1)$. The elements of the loading matrix ${\mathbf B}$ are drawn independently from the uniform distribution $U(-2,2)$. The target variables are generated by $y_{t+1}=\boldsymbol{\beta}_0'{\mathbf f}_t+\boldsymbol{\beta}_1'{\mathbf f}_{t-1}+{\varepsilon}_{t+1}$ with independent and identical errors ${\varepsilon}_t\sim N(0,1)$, that is, $q=2$ in Model ((ref)). We use 100 replications for each configuration $(T,N)$, where $T$ and $N$ are the sample size and the number of predictors, respectively. To make the results below replicable, the seed is set to be \texttt{1234} in the \texttt{R} programming. \subsection{In-Sample Forecasting} We first examine the in-sample forecasting errors with different competing methods. The coefficients in the linear forecasting model are set as $\boldsymbol{\beta}_0=(1,-0.8)'$ and $\boldsymbol{\beta}_1=(-1,2)'$. To create weak factors, we randomly choose $n$ rows in ${\mathbf B}$ as nonzero ones with $n\ll N$, and set all the remaining rows to zero. We consider two configurations of $(T,N)$ with $(T,N)=(200,300)$ and $(T,N)=(200,500)$, respectively. In practice, the time span of the predictors is from 1 to $T$ and that of the target variable $y$ is from $3$ to $T+1$. We compare the forecasting performance of the factors extracted by the proposed method (denoted by sdPCA), of the factors and their lagged variables extracted by the traditional PCA approach (denoted by PCA), of the factors extracted by the method in huang2022scaled (denoted by sPCA), and of the ones extracted by the diffusion-index model of Stock2002b (denoted by SW). Table (ref) presents the forecasting performance of the proposed supervised dynamic PCA method (sdPCA) with the aforementioned competing methods such as the PCA, sPCA, and SW, where the factors are of different degrees of weakness. We focus on their in-sample MSFEs defined as ((ref)) in this experiment under different settings, that is, we first estimate the factors and the coefficients in the forecasting model, and then examine the in-sample sum-of-squared residuals and report the MSFEs. We report the mean and the median of the in-sample MSFEs for each configuration of $(T, N,n)$ and each method in Table (ref). \begin{table}[h] \caption{The in-sample one-step ahead root mean squared forecast errors of $y_t$. Four methods are used: the proposed method (sdPCA), the one using the factors and their lagged ones extracted by PCA (denoted by PCA), the scaled PCA of huang2022scaled (denoted sPCA), and the original diffusion-index one in Stock2002b using the factors extracted by PCA without lagged ones (denoted by SW). Two configurations of the sizes are considered, $(T,N)=(200,300)$ and $(200,500)$. The factors are weak ones in the sense that only $n=40,30,20,10$ out of $N=300$ or $500$ predictors have non-zero loading coefficients on the common factors. The data generating process is that in ((ref)) with $r=2$, where the factors, idiosyncratic terms, and the random errors in the linear forecasting model are independently and identically drawn from the standard normal distribution. The non-zero loading elements are drawn independently from the Uniform distribution $U(-2,2)$. The regression coefficients are $\boldsymbol{\beta}_0=(1,-0.8)'$ and $\boldsymbol{\beta}_1=(-1,2)'$. 100 replications are used throughout the experiments.} {\begin{center} \begin{tabular}{cccccccccccc} \toprule \multicolumn{12}{c}{$(T,N)=(200,300)$}\\ \hline &\multicolumn{2}{c}{sdPCA}&&\multicolumn{2}{c}{PCA}&&\multicolumn{2}{c}{sPCA}&&\multicolumn{2}{c}{SW}\\ \cline{2-3}\cline{5-6}\cline{8-9}\cline{11-12} $n$&mean&median&&mean&median&&mean&median&&mean&median\\ 40&1.064&1.065&&1.080&1.084&&2.395&2.399&&2.414&2.419\\ 30&1.077&1.082&&1.112&1.113&&2.394&2.389&&2.415&2.410\\ 20&1.123&1.123&&1.203&1.205&&2.396&2.393&&2.425&2.422\\ 10&1.258&1.248&&1.858&1.844&&2.422&2.432&&2.526&2.529\\ \midrule \multicolumn{12}{c}{$(T,N)=(200,500)$}\\ \hline 40&1.054&1.055&&1.136&1.137&&2.406&2.392&&2.440&2.431\\ 30&1.093&1.089&&1.266&1.273&&2.403&2.380&&2.454&2.436\\ 20&1.129&1.123&&1.525&1.525&&2.387&2.376&&2.484&2.472\\ 10&1.325&1.329&&2.571&2.567&&2.303&2.307&&2.697&2.676\\ \bottomrule \end{tabular} \end{center}} \end{table} From Table (ref), we see that the proposed supervised dynamic PCA outperforms all the other three methods and it fits the model better in terms of the mean and median of the in-sample MSFEs. Specifically, the proposed sdPCA and the PCA methods produce the most accurate predictions in terms of the mean or median of the MSFEs, which is understandable since they both use the correct number of factors while the other two methods do not produce accurate predictions because they only use the contemporaneous factors without including any lagged information in prediction. Our proposed sdPCA method performs slightly better than the traditional PCA because the supervised procedure can strengthen the components that have non-zero loadings and mitigate the effect of the components with no loadings on the factors. Furthermore, we note that as $n$ decreases, i.e., when there are more zero-loadings or equivalently, the strength of the factors becomes weaker, the in-sample MSFEs tend to become larger, which is in agreement with our theory in the sense that the asymptotic rates in the Theorems of Section 3 will be higher if $\nu$ becomes smaller. In addition, the MSFEs tend to become larger when the number of predictors $N$ increases and this is also in line with the asymptotic theory. Finally, we see that the MSFEs produced by sPCA are smaller than those by SW, which is also reasonable since the scaled PCA may recover the first factor process more accurately using a supervised procedure than the latter unsupervised method. \subsection{Out-of-Sample Forecasting} Next, we examine the out-of-sample performance of the sdPCA and the other three methods using simulated data. The settings of all the parameters are the same as those in Section 3.1. For each configuration of $(T,N,n)$ and each iteration, we split the data into two sub-samples, one with the first $T_1=[3T/5]$ observations for training and the rest with $T_2=T-T_1$ data points for out-of-sample testing. Specifically, we adopt a rolling-window framework as follows. For the proposed method and the training samples $\{{\mathbf x}_1,...,{\mathbf x}_{T_1}\}$ and $\{y_3,...,y_{T_1}\}$, the supervised learning procedure and the estimation of the regression of the linear forecasting model are based on the training samples $\{{\mathbf x}_1,...,{\mathbf x}_{T_1-1}\}$ and $\{y_3,...,y_{T_1}\}$, because $h=1$. Then, the predictors are formed as $\{\widehat{\mathbf x}_2,...,\widehat{\mathbf x}_{T_1}\}$, where ${\mathbf x}_{T_1}$ of the forecast origin is used to obtain $\widehat{\mathbf x}_{T_1}$. Next, the factors are extracted as $\{\widehat{\mathbf g}_2,...,\widehat{\mathbf g}_{T_1}\}$. The regression coefficients in the linear forecasting model, denoted by $\widehat\boldsymbol{\beta}$, are estimated using $\{\widehat{\mathbf g}_2,...,\widehat{\mathbf g}_{T_1-1}\}$ and $\{y_3,...,y_{T_1}\}$, and finally, the forecast of $y_{T_1+1}$ is \[\widehat y_{T_1+1}=\widehat\boldsymbol{\beta}'{\mathbf g}_{T_1}.\] We continue this estimation and forecasting procedure by moving the next available data point to the first sub-sample and repeating the above process, and the root-MSFE is defined as \begin{equation} \text{RMSFE}=\sqrt{\frac{1}{T_2}\sum_{\tau=1}^{T_2}(y_{T_1+\tau+h-1}-\widehat y_{T_1+\tau+h-1})^2}, \end{equation} where $h=1$ is considered in the experiments. Table (ref) reports the mean and median of the out-of-sample RMSFEs for each method and each configuration of $(T,N,n)$, where $(T,N,n)$ are the same as those in Section 3.1. \begin{table}[h] \caption{The out-sample one-step ahead root mean squared forecast errors of $y_{t+1}$. Four methods are used: the proposed method (sdPCA), the one using the factors and their lagged ones extracted by PCA (denoted by PCA), the scaled PCA by huang2022scaled (denoted sPCA), and the original diffusion-index one in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). Two configurations of the sizes are considered: $(T,N)=(200,300)$ and $(200,500)$. The factors are weak ones in the sense that only $n=40,30,20,10$ out of $N=300$ or $500$ predictors have non-zero loadings on the common factors. The data generating process is that in ((ref)) with $r=2$, where the factors, idiosyncratic terms, and the random errors in the linear forecasting are independently and identically drawn from the standard normal distribution. The non-zero loading elements are drawn independently from the Uniform distribution $U(-2,2)$. The regression coefficients are $\boldsymbol{\beta}_0=(1,-0.8)'$ and $\boldsymbol{\beta}_1=(-1,2)'$. 100 replications are used throughout the experiments.} {\begin{center} \begin{tabular}{cccccccccccc} \toprule \multicolumn{12}{c}{$(T,N)=(200,300)$}\\ \hline &\multicolumn{2}{c}{sdPCA}&&\multicolumn{2}{c}{PCA}&&\multicolumn{2}{c}{sPCA}&&\multicolumn{2}{c}{SW}\\ \cline{2-3}\cline{5-6}\cline{8-9}\cline{11-12} $n$&mean&median&&mean&median&&mean&median&&mean&median\\ 40&1.120&1.126&&1.112&1.117&&2.427&2.425&&2.428&2.409\\ 30&1.134&1.137&&1.523&1.157&&2.428&2.440&&2.431&2.447\\ 20&1.188&1.181&&1.262&1.258&&2.438&2.467&&2.449&2.480\\ 10&1.136&1.350&&2.030&2.017&&2.508&2.508&&2.596&2.584\\ \midrule \multicolumn{12}{c}{$(T,N)=(200,500)$}\\ \hline 40&1.115&1.116&&1.179&1.183&&2.471&2.468&&2.482&2.479\\ 30&1.176&1.180&&1.334&1.343&&2.475&2.479&&2.490&2.505\\ 20&1.241&1.241&&1.644&1.654&&2.474&2.490&&2.543&2.539\\ 10&1.678&1.655&&2.639&2.612&&2.707&2.691&&2.734&2.714\\ \bottomrule \end{tabular} \end{center}} \end{table} From Table (ref), we see that the proposed sdPCA tends to outperform all the other three methods except for the case when $(n,T,N) = (40,200,300)$, which might be due to the errors incurred in finite samples when the idiosyncratic terms and the estimation errors in the supervised procedure are not dominated by the factors. As the factors become weaker, i.e., when $n$ becomes smaller, our proposed sdPCA outperforms all the other methods. Other findings are similar to those in the in-sample case in Section (ref), and we do not describe them to save space. Furthermore, we conduct simulations to confirm that the singularity of ${\mathbf B}_{\gamma}$ in the Assumption (ref) does not affect our estimation and forecasting. This case can happen when $\boldsymbol{\beta}_0=\boldsymbol{\beta}_1$ as mentioned in Section (ref). We consider the case of $\boldsymbol{\beta}_0=\boldsymbol{\beta}_1=(1,1)'$ in this experiment. Table (ref) presents the mean and median of the RMSFE for different methods and degrees of weakness with $(T,N)=(200,500)$. The settings of other parameters of the DGP are the same as before. One difference is that we recover two factors in our proposed sdPCA in this scenario while there are 4 factors extracted in the examples in Sections (ref) and (ref). From Table (ref), we can see that the performance of various methods is the same as those in Tables (ref) and (ref), and the sdPCA continues to work well under the extreme case that ${\mathbf B}_\gamma$ is singular. \begin{table}[h] \caption{The out-sample one-step ahead root mean squared forecast errors of $y_{t+1}$. Four methods are used: the proposed method (sdPCA), the one using the factors and their lagged ones extracted by PCA (denoted by PCA), the scaled PCA of huang2022scaled (denoted sPCA), and the original diffusion-index one in Stock2002b using the factors extracted by PCA without lagged ones (denoted by SW). Two configurations of the sizes are considered: $(T,N)=(200,300)$ and $(200,500)$. The factors are weak ones in the sense that only $n=40,30,20,10$ out of $N=300$ or $500$ predictors have non-zero loadings on the common factors. The data generating process is that in ((ref)) with $r=2$, where the factors, idiosyncratic terms, and the random errors in the linear forecasting are independently and identically drawn from the standard normal distribution. The non-zero loading elements are drawn independently from the Uniform distribution $U(-2,2)$. The regression coefficients are $\boldsymbol{\beta}_0=(1,1)'$ and $\boldsymbol{\beta}_1=(1,1)'$. 100 replications are used throughout the experiments.} {\begin{center} \begin{tabular}{cccccccccccc} \toprule \multicolumn{12}{c}{$(T,N)=(200,500)$,$\boldsymbol{\beta}=(1,1,1,1)'$}\\ \hline &\multicolumn{2}{c}{sdPCA}&&\multicolumn{2}{c}{PCA}&&\multicolumn{2}{c}{sPCA}&&\multicolumn{2}{c}{SW}\\ \cline{2-3}\cline{5-6}\cline{8-9}\cline{11-12} $n$&mean&median&&mean&median&&mean&median&&mean&median\\ 40&1.052&1.061&&1.079&1.083&&1.753&1.774&&1.759&1.780\\ 30&1.061&1.063&&1.106&1.110&&1.760&1.747&&1.773&1.762\\ 20&1.087&1.085&&1.193&1.203&&1.777&1.786&&1.810&1.808\\ 10&1.130&1.126&&1.465&1.466&&1.790&1.797&&1.903&1.913\\ \bottomrule \end{tabular} \end{center}} \end{table} \subsection{The LASSO Estimation} In this section, we conduct a simulation study to verify the efficacy of the proposed Lasso procedure. We set $\boldsymbol{\beta}_0=(6,3)'$ and $\boldsymbol{\beta}_1=(-5,0)'$ where the last coefficient in the linear forecasting model is zero. All the other settings of parameters and data-generating processes are the same as before. Due to the identification issues, there is a subtle change in the final regression step. Specifically, once we have obtained $\widehat\boldsymbol{\gamma}_i$ in the first step and hence a new loading matrix ${\mathbf B}_{\widehat\boldsymbol{\gamma}}$, we perform a singular-value decomposition on ${\mathbf B}_{\widehat\boldsymbol{\gamma}}$ and obtain the right singular matrix ${\mathbf V}$, then the final step is to perform a linear regression of $y_{t+1}$ on ${\mathbf V}\widehat{\mathbf g}_t$. We report the frequencies of the correct recoveries of the nonzero coefficients using Lasso through 100 replications in Table (ref). We see that, overall, the Lasso approach tends to recover the correct number of non-zero coefficients for moderately large $N$ and $T$ with weak factors. There are two findings from Table (ref). First, when the sample size increases with a fixed dimension $N$, the accuracy in recoveries may be improved and the frequencies may slightly decrease when the dimension $N$ increases for a fixed $T$, which is understandable as the Lasso approach can be more accurate when the sample size increases but it may create more errors when the dimension $N$ increases, which is in line with our asymptotic analysis in Theorem (ref). \begin{table} \caption{Frequencies of the correct recoveries of the nonzero coefficients using Lasso. Four configurations of the sizes are considered, $(T,N)=(100,250)$, $(200,250)$, $(200,400)$, and $(300,500)$. The factors are weak ones in the sense that only $n=60,50,40,30$ out of $N$ predictors have non-zero loadings on the common factors. The data generating process is that in ((ref)) with $r=2$, where the factors, idiosyncratic terms, and the random errors in the linear forecasting are independently and identically drawn from standard normal distributions. The non-zero loading elements are drawn independently from the Uniform distribution $U(-2,2)$. The regression coefficients are $\boldsymbol{\beta}_0=(6,3)'$ and $\boldsymbol{\beta}_1=(-5,0)'$, where only one position is zero. 100 replications are used throughout the experiments.} {\begin{center} \begin{tabular}{ccccc} \toprule & \multicolumn{4}{c}{$(T,N)$} \\ \cline{2-5} $n$ &$(100,250)$&$(200,250)$&$(200,400)$&$(300,500)$\\ \hline 60&0.94&0.94&0.86&0.66\\ 50&0.90&0.99&0.80&0.82\\ 40&0.87&0.91&0.79&0.85\\ 30&0.92&0.95&0.89&0.87\\ \bottomrule \end{tabular} \end{center}} \end{table} \section{Empirical Studies} In this section, we apply the proposed sdPCA and linear dynamic forecasting method to macroeconomic forecasting with the widely used U.S. monthly macroeconomic variables. To highlight the forecasting power of the proposed method, we compare the performance of sdPCA with some factor-based forecasting methods which are commonly used in the literature. Since the comparisons between the sPCA of huang2022scaled and the target PCA, PLS, and regularized methods have been studied in huang2022scaled, showing that the sPCA tends to dominate other methods in terms of the forecasting accuracy, we only compare our method with the sPCA, the factors and their lagged ones extracted by the traditional PCA, and the diffusion-index method in Stock2002b without including the lagged ones. \subsection{Data and In-Sample Results} We consider the macroeconomic variables studied by Stock2002a, mccracken2016fred, and huang2022scaled, among many others. The data are obtained from the FRED-MD data base which are maintained by St. Louis Fed\footnote{\url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}}. As described in mccracken2016fred, this data set extends the widely used Stock2002a set and covers broad economic categories including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). The groups of these variables are the same as those in mccracken2016fred while huang2022scaled re-grouped them into six ones. The detailed variables and transformation codes to ensure the stationarity of each macro variable are provided in the online data appendix. There are 127 variables in the online data set, but 4 of them are removed due to missing values therein. The remaining 123 macro variables are slightly different from those used in huang2022scaled because we only focus on the variables contained in the data file without adding new variables or replacing old ones with new variables. We consider the 123 macro variables spanning from July 1962 to December 2019 as all the series have no missing values during this period. Therefore, we have $N=123$ and $T=690$. We apply the proposed sdPCA to these 123 macro variables to forecast the 1-month ahead U.S. industrial production (IP) growth, change in the unemployment rate (UNRATE), change in the consumer price index: all (CPI-All), growth of real manufacturing and trade industries sales (M&T Sales), the S&P 500 index volatility change (Volatility Change), and the S&P 500 index return (Return), where the S&P 500 index volatility and the S&P 500 index return are obtained from the online data appendix of welch2008comprehensive\footnote{\url{https://sites.google.com/view/agoyal145}}. The S&P 500 index return are obtained from the \texttt{CRSP_SPvw} column, and the S&P 500 index volatility is the squared root of the \texttt{svar} column in the online data file, which is slightly different from those calculated in ludvigson2007empirical since the former does not subtract the risk-free rate when calculating the monthly volatility. To begin, we study the predictive power of each individual predictor and its lagged values on the variables of interest. In Figures (ref) and (ref), we plot the in-sample $R^2$s of predicting the 1-month ahead IP growth, change in UNRATE, change in M&T Sales, change in S&P 500 index volatility, and the S&P 500 index return by each of the 123 macro variables, respectively, where $q=2$ is used in Figure (ref) and $q=3$ in Figure (ref). Panels A and B indicate that the Labor Market (LM) conditions have the highest predictive power for future IP growth and unemployment rate, which is similarly found in huang2022scaled. In addition, the Output and income, Consumption, orders, and inventories, and the Interest and exchange rates also have higher predictive power than the remaining groups of variables. From Panel C, we see that the Prices have the highest predictive power for CPI-all, which is reasonable as the prices are directly related to inflation. Panel D indicates that the Out, LM, HS, COI, and IER have comparable predictive power for M&T sales, while the remaining variables do not have significant predictive power. For the Stock Market predictions in Panes E and F, we see that the Stock market prices have more predictive power for the volatility change, and the HS and LM conditions also have higher predictive power for the S&P 500 index returns. Furthermore, the plots in Figure (ref) suggest that the in-sample $R^2$ may be increased when we add more lagged variables in the linear forecasting. This is a common phenomenon in the autoregression context provided that the number of lagged values used is not too large. Overall, Figures (ref) and (ref) indicate that each predictor has different forecasting ability and their weights should be carefully assigned when extracting factors. \begin{figure} \begin{center} \caption{Bar charts of in-sample $R^2$s in predicting 1-month ahead industrial production growth (Panel A), unemployment rate (Panel B), consumer price index: all (Panel C), Real Manufacturing and Trade Industries Sales (Panel D), S&P 500 index Volatility Change (Panel E), and the S&P 500 index return (Panel F) by each of the 123 macro variables from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). We set $q=2$ in each regression. Macro variables are collected at a monthly frequency and the sample period is 1962:07--2019:12.} \end{center} \end{figure} \begin{figure} \begin{center} \caption{Bar charts of in-sample $R^2$s in predicting 1-month ahead industrial production growth (Panel A), unemployment rate (Panel B), consumer price index: all (Panel C), Real Manufacturing and Trade Industries Sales (Panel D), S&P 500 index Volatility Change (Panel E), and the S&P 500 index return (Panel F) by each of the 123 macro variables from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). We set $q=3$ in each regression. Macro variables are collected at a monthly frequency and the sample period is 1962:07--2019:12.} \end{center} \end{figure} Furthermore, we also use AIC to select the linear models used in Step 1 and plot the in-sample $R^2$s of predicting the 1-month ahead IP growth, change in UNRATE, change in M&T Sales, change in S&P 500 index volatility, and the S&P 500 index return by each of the 123 macro variables, respectively, in Figure (ref), where the maximal order is set to be $q=5$. From Figure (ref), we see that the overall pattern of the $R^2$'s in each plot is similar to its counterparts in Figures (ref) and (ref). There are some minor differences between the $R^2$ explained by the models with fixed lags and those selected by AIC. For example, the Consumption, orders, and inventories (COI) related variables have more predictive power for IP than those in Figures (ref)--(ref), and the Prices (PR) related variables produce slightly higher predictive power for the consumer price index: all (CPI-All). Nevertheless, the findings in Figures (ref)--(ref) remain valid in terms of variables that are related to the target variables according to their predictive power. \begin{figure} \begin{center} \caption{Bar charts of in-sample $R^2$s in predicting 1-month ahead industrial production growth (Panel A), unemployment rate (Panel B), consumer price index: all (Panel C), Real Manufacturing and Trade Industries Sales (Panel D), S&P 500 index Volatility Change (Panel E), and the S&P 500 index return (Panel F) by each of the 123 macro variables from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). We set $q=5$ in each regression and use AIC to select the model in Step 1. Macro variables are collected at a monthly frequency and the sample period is 1962:07--2019:12.} \end{center} \end{figure} Next, we consider the in-sample data analysis. First, we standardize each macroeconomic variable and calculate the eigenvalues of the resulting covariance of the 123 variables. That is, we perform eigen-value decomposition of the sample correlation matrix of the predictors. Based on the resulting eigenvalues, the first PCA factor explains about 18% of the total variation. When we apply the proposed sdPCA method to the data, the first sdPCA factor explains 21% to 52% of the total variation depending on the target variable of interest, which is higher than that explained by the first PCA factor. This suggests that the supervised PCA may improve the predictive power of the available predictors. We plot the loadings of the first to the sixth factors using the traditional PCA method in Figure (ref), where, for ease of reading, each loading vector is obtained by multiplying the corresponding eigenvector by 10. From the plot, we see that the first and the second PCA factors are more related to the real economic conditions and they have heavier loads on output and income, labor, and housing variables, followed by the interest and exchange rates. The third PCA factor depends mainly on price-related variables. The fourth factor has heavier loads on the interest rates and stock market conditions. The fifth factor has larger loadings on the interest rates while the sixth factor shows similar loading grouping as those of the first two factors. For comparison purposes, we also show the first six loadings of the proposed sdPCA in predicting the 1-month ahead IP growth, change of UNRATE, CPI-All, M&T Sales, Volatility Change, and the Return in Figures (ref)$-$(ref). From these plots, we see that the loadings are rather different from those in Figure (ref). For example, from Figure (ref), we see that the first factor in predicting IP growth has loads mainly on OUT, LM, and COI conditions while the effect of the HS seen in the first unsupervised PCA factor has decreased. Similar results are also found for UNRATE where the first sdPCA factor mainly depends on the OUT, LM, and COI. From Figure (ref), the first factor in predicting CPI-All is related mainly to the price variables resulting in certain adjustments to those in the first unsupervised PCA factor. The first factor in predicting M&T sales depends relatively heavier on the OUT, LM, and HS than on the others. For the stock market data predictions, we find from Figure (ref) that the first sdPCA factor in predicting the S&P volatility change has heavier loads on the SM variables, and this is understandable as they are more directly related. On the other hand, we also find that the first sdPCA factor in predicting the Return is more closely related to the HS conditions, implying that the stock returns and the Housing conditions are related, which, in turn, shows that investors may switch investment between the stock market and the real estate market. \begin{figure} \begin{center} \caption{Bar charts of loadings on the 123 macro variables of the first six PCA factors. The macro variables are collected at a monthly frequency from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). The sample period is 1962:07--2019:12.} \end{center} \end{figure} \begin{figure} \begin{center} \caption{Bar charts of loadings on the 123 macro variables of the first six factors of the proposed sdPCA in predicting the 1-month ahead industrial production growth (OP) and the unemployment rate (UNRATE). The macro variables are collected at a monthly frequency from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). The sample period is 1962:07--2019:12.} \end{center} \end{figure} \begin{figure} \begin{center} \caption{Bar charts of loadings on the 123 macro variables of the first six factors of the proposed sdPCA in predicting the 1-month ahead consumer price index: all (CPI-All) and Real Manufacturing and Trade Industries Sales (M&T Sales). The macro variables are collected at a monthly frequency from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). The sample period is 1962:07--2019:12.} \end{center} \end{figure} \begin{figure} \begin{center} \caption{Bar charts of the loadings on the 123 macro variables of the first six sdPCA factors in predicting the 1-month ahead S&P 500 index Volatility Change (Volatility Change), and the S&P 500 index return (Return). The macro variables are collected at a monthly frequency from the FRED-MD data set of mccracken2016fred, consisting of eight groups including the output and income (OUT), Labor market(LM), Housing (HS), Consumption, orders, and inventories (COI), Money and credit (MC), Interest and exchange rates (IER), Prices (PR), and Stock market (SM). The sample period is 1962:07--2019:12.} \end{center} \end{figure} \subsection{Out-of-Sample Forecasting} In this subsection, we assess the performance of the proposed method using out-of-sample forecasting experiments. For each target variable of interest, we split the sample into two subsamples, where the first one consists of the first 80% of the data for modeling and the second subsample of the remaining 20% of the data for out-of-sample prediction. We also adopt the rolling-window scheme as that in the simulation studies, that is, we train the factors and the forecasting coefficients using the first subsample to predict the next target data point. Then we repeat the above procedure after moving the next available observation of predictors and target variable from the second subsample to the first one to obtain the next prediction. This rolling-window scheme is terminated when there is no more observation to compute forecasting error. Similarly to the experiments in Section (ref), we compare the proposed method (sdPCA) with the traditional PCA (denoted by PCA), the sPCA in huang2022scaled (denoted by sPCA), and the diffusion-index forecast in Stock2002a (denoted by SW), where the lagged variables of the PCA factors are also included as predictors in the PCA method, and only the contemporaneous factors are used in the linear forecasting of the SW method. The forecasting performance is measured by the root MSFE defined in ((ref)). We use the autoregressive (AR) model with order 1 or 2 as benchmark methods in the comparison. For the forecasting of IP, UNRATE, CPI-All, and M&T Sales, we consider $h=1,2,3,4,$ and $5$-steps ahead predictions, and we only consider $h=1$ step ahead prediction for the financial data of volatility change and the stock market return, which are of major interest in most financial market predictions. Tables (ref)--(ref) report the results of $h=1,2,3,4$, and $5$-step ahead predictions of the IP, UNRATE, CPI-All, and M&T Sales, respectively. In the comparison, please note that Model ((ref)) implies that the number of factors used in sdPCA is $qs$ if the number of contemporaneous factors used in PCA, sPCA, and the SW is $s$ for an integer $s>0$, where $q$ is the number of lagged variables used in Steps 1 to 3 of the proposed procedure. For the proposed sdPCA, $q=2$ and $q=3$ are employed in the empirical studies, indicating the number of lagged variables used in the forecasting. The number of factors used ranges from 1 to 3 for the methods of PCA, sPCA, and the SW, and therefore, the number of factors used in sdPCA ranges from 1 to 6 when $q=2$ and from 1 to 9 if $q=3$ according to the above discussion. For each $h$, the smallest value is marked in boldface. The Lasso procedure is considered with the corresponding errors given in the parentheses if sdPCA does not beat other methods. We only report the results of the PCA method when $q=2$ because those with $q=3$ do not show any clear improvement in most cases due to the possibility of overfitting. On the other hand, our analysis suggests that the forecasts with more than 3 factors in PCA, sPCA, and the SW do not necessarily improve the prediction accuracy, and therefore, we only compare the results when the number of factors used ranges from 1 to 3 for the methods of PCA, sPCA, and the SW, and the corresponding number of factors used in sdPCA ranges from 1 to 6 when $q=2$ and from 1 to 9 if $q=3$. From Table (ref), we see that the smallest prediction error is achieved by the proposed sdPCA for each $h$, and the prediction using more lagged variables ($q=3$) can improve the forecasting ability ($q=2$) in most cases. Furthermore, the PCA method tends to produce the second smallest errors since it also includes the lagged variables as predictors, which implies that the lagged variables can be used to improve the forecasting performance of IP. In addition, all the methods outperform the benchmark AR methods. For the predictions of the UNRATE in Table (ref), we note that the sdPCA outperforms other competing methods for short-term predictions, such as the cases when $h=1,2$ and $3$, without using the Lasso approach. For $h=4$ and $5$, the errors accumulated by sdPCA are increasing, but the Lasso approach can significantly reduce the forecasting errors and a simple Lasso procedure can produce even smaller errors than all the other methods, implying that the penalized method is an effective way in selecting the factors that have more predictive power. Similar results are also found for the predictions of CPI in Table (ref), where the proposed sdPCA as well as the PCA with lagged factors can outperform other methods in short-term forecasting ($h=1$ and $2$), and the forecasting performance for long-term ahead predictions are not more accurate than those by simple AR approaches. In addition, the Lasso approach can also improve the forecasting performance as shown in the case of using sdPCA with $q=3$. Similar findings are also obtained in predicting the M&T sales in Table (ref), and we omit the details. From Tables (ref)--(ref), we see that the proposed sdPCA might produce increased errors for $h$ = 4 and 5, especially when the number of factors used increases. There are two possible explanations. First, for stationary time series, such as those considered in our example, the serial dependence decays exponentially so that as $h$ increases the information of the target variable embedded in the predictors decreases. The forecast errors may increase as $h$ increases with uncertainty approaching the unconditional variance of the target variable. Second, increasing the number of factors used also increases the possibility of overfitting, which may lead to inferior prediction. On the other hand, the Lasso procedure can significantly improve the forecasts, which confirms the overfitting issue when more factors are used without variable selection. This highlights the importance of using the Lasso procedure to select the relevant factors if we do not know how many factors to include in linear forecasting. Finally, the one-step-ahead predictions of the stock return and the volatility change are shown in Table (ref). From Table (ref), we see that the performance of the sdPCA is comparable with that of SW in forecasting the Stock return, and the sPCA method cannot beat the benchmark methods overall. For the predictions of the volatility change, we find that the proposed sdPCA outperforms all the other competing methods. The results in Table (ref) suggest that the proposed sdPCA could be helpful in forecasting financial data. \begin{table}[htbp] \caption{Out-of-sample root mean squared forecast errors in predicting IP. Four comparison methods are used: the proposed method (sdPCA), the one using the factors and their lagged variables extracted by traditional PCA (denoted by PCA), the scaled PCA by huang2022scaled (denoted sPCA), and the original diffusion-index method in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). We also use AR(1) and AR(2) models as benchmark methods in out-of-sample forecasting. For the proposed sdPCA, $q=2$ and $q=3$ are considered. For each $h$, the smallest MSFE is marked in boldface. } {\begin{center} \begin{tabular}{ccccccccccc} \toprule &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&0.080&0.078&0.079&0.080&0.079&0.076&&0.085&0.078&0.078\\ 2&0.081&0.079&0.080&\boxed{\bf 0.075}&0.076&0.079&&0.086&0.082&0.081\\ 3&0.085&0.082&0.084&0.082&0.083&0.083&&0.086&0.084&0.083\\ 4&0.088&\boxed{\bf 0.085}&0.088&0.086&0.093&0.127&&0.086&0.086&0.088\\ 5&0.090&0.087&0.089&0.087&0.089&0.228&&0.088&0.088&0.089\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&0.083&0.080&0.081&&0.085&0.079&0.080&&0.085&0.084\\ 2&0.086&0.082&0.084&&0.086&0.082&0.082&&0.089&0.085\\ 3&0.088&0.084&0.086&&0.086&0.084&0.083&&0.090&0.085\\ 4&0.089&0.086&0.088&&0.086&0.086&0.086&&0.091&0.088\\ 5&0.091&0.094&0.091&&0.087&0.089&0.089&&0.091&0.090\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts& 8fts&9fts&\\ 1&0.078&0.077&0.077&0.074&0.075&\boxed{\bf 0.073}&0.074&0.074&0.074\\ 2&0.080&0.079&0.079&0.076&0.077&0.077&0.079&0.078&0.078\\ 3&0.084&0.082&0.084&\boxed{\bf 0.080}&0.082&0.083&0.083&0.081&0.081\\ 4&0.086&\boxed{\bf 0.085}&0.087&\boxed{\bf 0.085}&0.089&0.232&0.181&0.182&0.134\\ 5&0.089&0.088&0.090&\boxed{\bf 0.086}&0.089&0.096&0.110&0.109&0.120\\ \bottomrule \end{tabular} \end{center}} \end{table} \begin{table}[htbp] \caption{Out-of-sample root mean squared forecast errors in predicting UNRATE. Four comparison methods are used: our proposed method (sdPCA), the one using the factors and their lagged variables extracted by PCA (denoted by PCA), the scaled PCA by huang2022scaled (denoted sPCA), and the original diffusion-index model in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). We also use AR(1) and AR(2) models as benchmark methods in out-of-sample forecasting. For the proposed sdPCA, $q=2$ and $q=3$ are considered. For each $h$, the smallest MSFE is marked in boldface. The Lasso procedure is applied and the corresponding errors are placed in parentheses if sdPCA does not outperform the other methods.} {\begin{center} \begin{tabular}{ccccccccccc} \toprule &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&1.713&0.694&1.703&1.714&1.682&1.627&&1.976&1.773&1.767\\ 2&1.772&1.764&1.795&1.781&1.757&\boxed{\bf 1.680}&&2.004&1.831&1.832\\ 3&1.821&1.897&1.900&1.834&1.867&1.852&&2.009&1.850&1.893\\ 4&8.769&3.082&3.657&4.013&4.887&4.783&&2.024&1.890&1.939\\ &&&&&&(1.914)\\ 5&2.115&4.370&4.730&4.850&4.903&4.774&&2.029&1.931&1.951\\ &&&&&&(\boxed{\bf 1.924})\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&1.736&1.734&1.707&&1.973&1.784&1.791&&2.044&1.941\\ 2&1.815&1.817&1.833&&1.996&1.859&1.852&&2.080&1.960\\ 3&1.859&1.919&1.912&&2.006&1.877&1.897&&2.090&2.027\\ 4&1.882&1.951&1.956&&2.014&1.906&1.937&&2.092&2.035\\ 5&2.272&4.782&4.752&&2.028&1.939&1.955&&2.093&2.069\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts&8 fts&9 fts&\\ 1&1.696&1.683&1.698&1.674&{1.600}&\boxed{\bf 1.594}&1.594&1.595&1.594\\ 2&1.755&1.804&1.808&1.756&1.733&1.696&1.708&1.711&1.704\\ 3&\boxed{\bf 1.803}&2.380&2.108&2.515&2.579&2.590&2.574&2.572&2.629\\ 4&10.011&2.604&3.936&4.207&4.942&4.986&4.798&4.849&4.868\\ &&&&&&(\boxed{\bf 1.881})\\ 5&8.353&4.835&5.299&5.286&5.082&4.976&4.991&5.032&5.129\\ \bottomrule \end{tabular} \end{center}} \end{table} \begin{table}[htbp] \caption{The out-sample root mean squared forecast errors for CPI-ALL. Four comparison methods are used: our proposed method (sdPCA), the one using the factors and their lagged variables extracted by PCA (denoted by PCA), the scales PCA by huang2022scaled (denoted sPCA), and the original diffusion-index method in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). We also use AR(1) and AR(2) models as benchmark methods to do out-of-sample forecasting. For the proposed sdPCA, $q=2$ and $q=3$ are considered. For each $h$, the smallest value is marked in boldface. The Lasso procedure is considered with the corresponding errors given in parentheses if sdPCA does not beat other methods.} {\begin{center} \begin{tabular}{ccccccccccc} \toprule &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&0.034&0.032&0.032&0.032&\boxed{\bf 0.031}&0.032&&0.034&0.034&0.034\\ 2&\boxed{\bf 0.032}&0.033&0.033&0.033&0.033&0.034&&0.034&0.034&\boxed{\bf 0.032}\\ 3&\boxed{\bf 0.034}&0.036&0.036&0.036&0.036&0.036&&\boxed{\bf 0.034}&0.034&0.034\\ 4&0.035&0.035&0.035&0.035&0.035&0.035&&0.034&0.034&0.034\\ &&&&&&(0.034)\\ 5&0.052&0.062&0.064&0.055&0.045&0.044&&\boxed{\bf 0.034}&0.034&0.035\\ &&&&&&(\boxed{\bf 0.034})\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&0.035&0.033&0.033&&0.034&0.034&0.035&&0.035&0.034\\ 2&0.033&0.033&0.033&&0.034&0.034&0.033&&0.035&0.034\\ 3&\boxed{\bf 0.034}&0.035&0.036&&\boxed{\bf 0.034}&0.034&0.034&&\boxed{\bf 0.034}&0.035\\ 4&0.034&0.034&0.034&&0.034&0.034&0.034&&0.034&0.034\\ 5&0.052&0.052&0.049&&\boxed{\bf 0.034}&0.034&0.034&&\boxed{\bf 0.034}&0.034\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts&8fts&9 fts\\ 1&0.033&0.032&0.031&0.032&0.032&0.032&0.033&0.032&0.032\\ 2&0.032&0.033&0.033&0.033&0.034&0.034&0.035&0.035&0.035\\ 3&0.035&0.034&0.036&0.036&0.036&0.034&0.038&0.039&0.039\\ 4&0.035&0.034&0.036&0.036&0.036&0.036&0.037&0.037&0.037\\ &&&&&&(\boxed{\bf 0.032})\\ 5&0.060&0.072&0.074&0.058&0.056&0.055&0.052&0.051&0.051\\ \bottomrule \end{tabular} \end{center}} \end{table} \begin{table}[htbp] \caption{The out-sample root mean squared forecast errors for Real Manu. and Trade Industries Sales. Four comparison methods are used: our proposed method (sdPCA), the one using the factors and their lagged variables extracted by PCA (denoted by PCA), the scales PCA by huang2022scaled (denoted sPCA), and the original diffusion-index method in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). We also use AR(1) and AR(2) models as benchmark methods to do out-of-sample forecasting. For the proposed sdPCA, $q=2$ and $q=3$ are considered. For each $h$, the smallest value is marked in boldface. The Lasso procedure is considered with the corresponding errors given in parentheses if sdPCA does not beat other methods.} {\begin{center} \begin{tabular}{ccccccccccc} \toprule &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&0.081&\boxed{\bf 0.078}&0.081&0.082&0.081&0.082&&0.081&\boxed{\bf 0.078}&0.080\\ 2&0.081&\boxed{\bf 0.079}&0.082&0.080&0.080&0.080&&0.081&0.081&0.083\\ 3&0.084&\boxed{\bf 0.081}&0.085&0.083&0.083&0.084&&0.082&0.082&0.083\\ 4&0.086&0.277&0.263&0.291&0.300&0.280&&\boxed{\bf 0.084}&0.084&0.084\\ &&&&&&(\boxed{\bf 0.084})\\ 5&0.087&0.286&0.255&0.289&0.284&0.295&&\boxed{\bf 0.084}&0.085&0.085\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&0.084&0.081&0.084&&0.082&0.081&0.081&&0.090&0.090\\ 2&0.084&0.080&0.083&&0.081&0.081&0.082&&0.086&0.086\\ 3&0.085&0.083&0.084&&0.082&0.083&0.084&&0.087&0.087\\ 4&0.087&0.090&0.088&&0.084&0.085&0.085&&0.087&0.087\\ 5&0.087&0.212&0.235&&\boxed{\bf 0.084}&0.085&0.086&&0.087&0.087\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts&8 fts&9 fts\\ 1&0.080&0.078&0.080&0.080&0.080&0.081&0.085&0.088&0.086\\ 2&0.081&0.080&0.083&0.080&0.080&0.080&0.083&0.083&0.083\\ 3&0.083&0.250&0.372&0.281&0.268&0.265&0.277&0.257&0.258\\ 4&0.085&0.285&0.286&0.334&0.337&0.316&0.328&0.332&0.334\\ 5&0.086&0.356&0.304&0.322&0.313&0.317&0.309&0.314&0.311\\ &&&&&&(\boxed{\bf 0.084})\\ \bottomrule \end{tabular} \end{center}} \end{table} \begin{table}[htbp] \caption{The out-sample root mean squared forecast errors for S&P 500 index returns and change of S&P 500 index volatility. Four comparison methods are used: our proposed method (sdPCA), the one using the factors and their lagged variables extracted by PCA (denoted by PCA), the scales PCA by huang2022scaled (denoted sPCA), and the original diffusion-index method in Stock2002b using the factors extracted by PCA without lagged variables (denoted by SW). We also use AR(1) and AR(2) models as benchmark methods to do out-of-sample forecasting. For the proposed sdPCA, $q=2$ and $q=3$ are considered. For each $h$, the smallest value is marked in boldface. The Lasso procedure is considered with the corresponding errors are placed in parentheses if sdPCA does not beat other methods.} {\begin{center} \begin{tabular}{ccccccccccc} \toprule \multicolumn{11}{c}{Stock return with dividends}\\ \hline &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&0.495&0.505&0.514&0.506&0.534&0.517&&0.495&0.503&0.510\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&0.498&0.506&0.511&&\boxed{\bf 0.494}&0.504&0.506&&0.498&0.499\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts&8 fts&9 fts\\ 1&\boxed{\bf 0.494}&0.506&0.508&0.510&0.518&0.525&0.510&0.515&0.522\\ \bottomrule \toprule \multicolumn{11}{c}{Change of Volatility}\\ \hline &\multicolumn{6}{c}{sdPCA ($q=2$)}&&\multicolumn{3}{c}{PCA ($q=2$)}\\ \cline{2-7}\cline{9-11} $h$&1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&&1 ft&2 fts&3 fts\\ 1&\boxed{\bf 0.270}&0.272&0.270&0.279&0.278&0.271&&0.279&0.276&0.278\\ \toprule &\multicolumn{3}{c}{sPCA}&&\multicolumn{3}{c}{SW}&&\multicolumn{2}{c}{AR}\\ \cline{2-4}\cline{6-8}\cline{10-11} &1 ft&2 fts& 3 fts&&1 ft& 2fts& 3fts&&AR(1)&AR(2)\\ 1&0.273&0.271&0.278&&0.280&0.280&0.280&&0.284&0.286\\ \toprule &\multicolumn{9}{c}{sdPCA ($q=3$)}&\\ \cline{2-10} &1 ft&2 fts &3 fts&4 fts&5 fts&6 fts&7 fts&8 fts&9 fts\\ 1&0.270&0.280&0.283&0.277&0.275&0.274&0.274&0.277&0.280\\ \bottomrule \end{tabular} \end{center}} \end{table} \section{Conclusion} This paper introduced a new supervised dynamic PCA method for linear forecasting with many predictors, which is commonly seen in the big-data environment. The new supervised PCA provides an effective way to bridge the gap between predictors and the targeted variables of interest by scaling and combining information of the predictors and their lagged variables, which is in line with dynamic forecasting. Furthermore, we also proposed to use penalized methods, such as the LASSO approach, to select the significant factors that have more predictive power than the others in the linear forecasting equation. To highlight the prediction power achieved by the proposed method, we showed that our estimators are consistent and outperform theoretically the traditional methods in prediction under some commonly used conditions. We conducted extensive simulations to verify that the proposed method produces satisfactory forecasting results and outperforms most of the existing methods using the traditional PCA. A real data example on predicting U.S. monthly macroeconomic variables using a large number of predictors shows that our method performs better than most of the existing ones in forecasting the U.S. industrial production (IP) growth, change in the unemployment rate (UNRATE), the consumer price index: all (CPI-All), the S&P 500 index volatility change (Volatility Change), and the S&P 500 index return with 123 macro variables from FRED-MD. Finally, the proposed sdPCA together with the Lasso procedure produces even more satisfactory results in many cases. Overall, the proposed procedure provides a comprehensive and effective method for dynamic forecasting when the data dimension is large.

titlepage\begin{center} { Internet Appendix for\\Supervised Dynamic PCA: Linear Dynamic Forecasting with Many Predictors} \end{center} \thispagestyle{empty} \begin{abstract} The first part of the Internet Appendix provides detailed proofs for all the theoretical statements in the main text, and the second part presents descriptions of the real data used in the main article. \noindentKeywords: Dynamic Forecasting, Factor Analysis, Supervised Principal Components, Large-Dimension, LASSO \noindentJEL classification: C22, C23, C38, C53 \end{abstract}

lemmaLet Assumptions (ref)$-$(ref) hold. For each $t$, (i) $E\|N^{-\nu/2}\sum_i {\mathbf b}_i u_{i,t}\|^2\le C$, (ii) for each $i$, $\frac{1}{NT} {\mathbf u}_t'{\mathbf U}'{\mathbf F}=O_p(\delta_{NT}^2)$, where $\delta_{NT}=\min\{\sqrt{N},\sqrt{T}\}$, (iii) $E\|T^{-1/2} \sum_t {\mathbf f}_t u_{i,t}||^2\le C$, (iv) $\frac{1}{N^\nu T}{\mathbf u}_i'{\mathbf U}{\mathbf B}=O_p(\frac 1 {N^\nu}) +O_p(\frac 1 {\sqrt{T N^\nu}})$; (v) ${\mathbf B}' {\mathbf U}' {\mathbf F} = \sum_{i=1}^N\sum_{t=1}^T {\mathbf b}_i {\mathbf f}_t' u_{i,t} =O_p( \sqrt{N^\nu T})$.

lemmaLet ${\mathbf Z}_i=(z_{i,q},...,z_{i,T-h})'$ and $\widetilde {\mathbf V}$ be the diagonal matrix consisting of the top $qr$ eigenvalues of $\sum_{i=1}^N{\mathbf Z}_i{\mathbf Z}_i'$ as its diagonal elements. Under Assumptions 1$-$6, if $N^{1-\nu}/T^2\rightarrow 0$, with probability tending to one, we have \[\widetilde{\mathbf V}\asymp N^{\nu}T.\]

lemmaUnder Assumptions 1$-$6, if $N^{1-\nu}/T^2\rightarrow 0$ and the identification condition ${\mathbf G}'{\mathbf G}/T={\mathbf I}_{rq}$ hold, then we have \[{\mathbf H}{\mathbf H}'-{\mathbf I}_{rq}=O_p(N^{-\nu/2}+T^{-1}+N^{1-\nu}/T^2)\,\,\text{and}\,\,{\mathbf H}'{\mathbf H}-{\mathbf I}_{rq}=O_p(N^{-\nu/2}+T^{-1}+N^{1-\nu}/T^2).\]

lemmaLet Assumptions 1$-$6 hold. If $N^{1-\nu}/T^2\rightarrow 0$, we have \[\frac{1}{T}(\widehat{\mathbf F}-{\mathbf F}{\mathbf H}')'{\mathbf F}=O_p(\frac{1}{N^\nu}+\frac{1}{T}),\] and \[\frac{1}{T}(\widehat{\mathbf F}-{\mathbf F}{\mathbf H}')'\widehat{\mathbf F}=O_p(\frac{1}{N^\nu}+\frac{1}{T}).\] As a result, Under the identification conditions that ${\mathbf G}'{\mathbf G}/T={\mathbf I}_{rq}$ and ${\mathbf B}_\gamma'{\mathbf B}_\gamma$ is a diagonal matrix with distinct diagonal elements, we also have \[{\mathbf H}={\mathbf I}_{rq}+O_p(\frac{1}{N^\nu}+\frac{1}{T}).\]