EconBase
← Back to paper

Rolling Window Selection in FAR Models with Structural Instabilities

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

45,453 characters

Rolling Window Selection in FAR Models with Structural Instabilities




\def\spacingset#1{\renewcommand{1}
{#1}\small\normalsize} \spacingset{1}



\if00
{
  \title{\bf{Rolling Window Selection in FAR Models with Structural  Instabilities}\thanks{I am grateful for comments from  Silvia Goncalves, Jiaying Gu, Peter Hansen, Benoit Perron, Yuanyuan Yan, and seminar participants at the 2026 Africa Meeting of the Econometric Society, University of Toronto,   95th Southern Economic Association, Conference Celebrating James MacKinnon at Queen's University, 2025 RESA Conference, 2025 CIREQ Econometrics Conference,
  7th International Conference on Econometrics and Statistics, and 17th International Conference on Computational and Financial Econometrics. The author gratefully acknowledges
financial support from the Social Sciences and Humanities Research
Council (SSHRC) of Canada.}}
  \author{Antoine A. Djogbenou \hspace{.2cm}\\
    Department of Economics, York University\\
    }
  \maketitle
} \fi

\if10
{
  \bigskip
  \bigskip
  \bigskip
  \begin{center}
    {\LARGE\bf Window Selection in FAR Models with Structural  Instabilities}
\end{center}
  \medskip
} \fi


\begin{abstract}
The paper develops a theory for selecting the rolling window when generating out-of-sample forecasts with factor-augmented regression (FAR) models in the presence of structural instabilities. It shows how to select a rolling window by minimizing the conditional mean squared forecast error (MSFE) while accounting for uncertainty in factor estimation. Because the conditional MSFE is unobserved and the factors are latent, this paper proposes a feasible version of the criterion and derives conditions under which the new method is asymptotically loss-efficient. A simulation experiment documents the procedure's performance.
\end{abstract}


\noindent
\ \ \ \ \ \ \ \ {\it Keywords:}  Forecasting, factor models, structural instabilities, rolling window.

\noindent
\ \ \ \ \ \ \ \ {\it JEL Codes:}  C18, C38, C53, C55.
\vfill

\newpage
\spacingset{1.45}




\section{Introduction}
 Large datasets of financial and macroeconomic variables are currently easily accessible. In many policy-relevant works, empirical researchers have found it useful to extract common factors from a large set of economic variables. Originally used in economics by \citet*{sargent1977business} and \citet*{geweke1977dynamic}, factor models have been employed to model the dynamics of these common factors, which involve only a few time-series indexes that summarize information in a high-dimensional set of economic variables.
 \citet*{stock2002forecasting} popularized the use of factor-augmented regressions (FARs) or diffusion indices to model a large set of macroeconomic variables and produce forecasts. In practice, a common approach is to rely on the principal component method (PCM) to extract factors and use the estimated factors as predictors in a forecasting equation.
 For a complete discussion of the usefulness or further applications in applied macro-econometrics of factor-augmented models for forecasting, see, for example, \citet*{BB2005}, \citet*{DD2013}, \citet*{ABT2015}, \citet*{Djogbenou2020}, \citet*{Djogbenou2024}, and \citet*{DS2024}.
 PCM's popularity as a factor extraction procedure stems from its computational simplicity and well-known properties established in \citet*{BN2002}, \citet*{Bai2003}, and \citet*{BN2006}.
 \citet*{BN2002}  proposed criteria for selecting the optimal number of factors.  \citet*{Bai2003}  provided a theory for inference of the estimated factors and their loadings.  \citet*{BN2006} derived the conditions for the validity of the inference and prediction intervals for a regression model with estimated factors as predictors. \citet*{BN2013} discussed the conditions for identifying the factors. \citet*{GP2014} derived a formula for the asymptotic bias in the distribution of regression coefficients when using estimated rather than observed factors, and proposed bootstrap approaches that account for uncertainty in the estimated factors.  \citet*{DGP2015}  extends their results to account for time dependence in $h$-periods-ahead forecasting factor augmented regression models.

 As is well known, very good in-sample predictive ability does not by itself ensure good out-of-sample predictive ability.  \citet*{Rossi2005} pointed out that although common econometric methods of identifying potentially useful predictors rely on in-sample significance tests, these approaches provide little assurance that the identified predictors are reliable and stable over time. The author discusses optimal model selection tests for regression models with observed regressors in the presence of structural instability. \citet*{CM2003} also discussed the harmful effect of structural change on forecast comparisons in the presence of structural change in existing tests. A recent work by \citet*{GMP2017} showed that the encompassing test proposed by \citet*{CM2001} and the Granger causality test by \citet*{Mccracken2007} are asymptotically valid when factors are estimated recursively. However, when a rolling-window version of these tests is employed, the effect of structural instability in the factor loadings and regression coefficients is unclear. It is also unclear how much past information to use for forecasting.

Moreover, in the presence of structural instability, it is natural to consider the end of the sample as a solution to the economy's structural change. However, the proximity of the crises may not provide sufficient data, and existing tests of predictive ability rely on asymptotic arguments.  \citet*{HI2015} argued that factor dynamics change when structural breaks occur because factors' sensitivities change in large data sets. Thus, neglecting these structural changes may reduce efficiency and worsen forecast performance.  They proposed the Lagrange multiplier and Wald-type tests for structural breaks.
 \citet*{BE2011} also suggested an alternative Lagrange multiplier test based on the idea that the number of factors increases when a structural change occurs in the large panel. In this paper, we focus on determining the optimal forecasting window size in the presence of estimated factors and structural instability. Several papers have developed procedures for selecting the window size $R$ in linear regression models with observed regressors. For example, \citet{IJR2017} propose minimizing an approximation to the conditional MSFE, since it is unknown in practice.  However, none of the existing procedures discuss selecting $R$ in the presence of estimated factors. In FAR, the latent factors are unknown and can be identified only up to rotation.  \citet{BN2009} and \citet{Djogbenou2021} showed that the MSFE is inflated in finite samples by a factor-estimation uncertainty term. The implications of these aspects for selecting the window $R$ remain unstudied. This research fills this gap.


This paper develops a method for selecting the window size when forecasting in the presence of time-varying parameters and latent factors. It proposes a procedure that minimizes an estimator of the forecaster's
loss function, that is, the conditional MSFE
up to a constant. Because the forecaster's loss function depends on rotated parameters and rotated latent factors, we propose a feasible version that relies on local linear estimators of the regression coefficients and factor loadings.  The paper studies the asymptotic loss efficiency of the new approach and derives conditions for its validity. Simulation experiments document the approach's performance and show a substantial reduction in mean squared forecast error under the new procedure.

The remainder of the paper is organized as follows.  In \Cref{Settings}, we discuss the framework. In \Cref{Main}, we present the assumptions and the paper's main theoretical results. \Cref{Simulation} contains the simulation results, while
\Cref{Conclusion} concludes. The proofs are relegated to Appendix A. We denote
$\left\Vert \bm M\right\Vert =\left(
\mathrm{Trace}\left( \bm M^{\top} \bm M\right)  \right)  ^{1/2}$, the Euclidean
norm, $\bm Q$ $>0$, the positive definiteness for any square matrix $\bm Q$, and
$C$, a generic finite constant.

\section{Settings and Framework}
\label{Settings}
Consider the following time-varying coefficient factor-augmented regression model
\begin{equation}
y_{t+h}=(\bm\alpha(t/T))^\top \bm f_t  +(\bm\beta(t/T))^\top\bm w_t +\varepsilon_{t+h}, t=1,\ldots, T-h,
\label{FAR}
\end{equation}
where $t$ indicates the time period, $T$ is the time dimension,  $h$ is the horizon of forecast, and $\varepsilon_{t+h}$ is the forecast error.
The model assumes the dependent variable $y_{t+h}$ is explained by $q$ observed regressors, $\bm w_t$, and $r$ latent factors, $\bm f_t$, that summarize the comovement in a large set $N$ of variables, $x_{it},i=1,\ldots,N$, through the approximate factor model
\begin{equation}
x_{it}=(\bm\lambda_i(t/T))^\top \bm f_t +e_{it}, i=1,\ldots,N, t=1,\ldots, T,
\label{panel}
\end{equation}
where $e_{it}$ is the idiosyncratic error and the factor loading $\lambda_i(t/T)$ is an $r$-vector of unknown deterministic functions of $t$. The considered FAR model assumes the regression coefficients and the factor loadings smoothly change over time.

Our interest is to predict $y_{T+h}$ using the past information and the most recent relevant information that are available, i.e., $\left(y_t, w_t, x_{1t}, \ldots, x_{Nt}\right), t=T-R+1,\ldots, T,$ where $R$ needs to be determined. We contribute to the literature by developing a theory for selecting the window $R$ that is robust to factor-estimation uncertainty and structural changes in factor loadings and regression coefficients.

For any window $R$, the rolling window approach estimates parameters (factor loadings and regression coefficients) at the terminal date, where $t/T=1$,
and uses $R$ recent observations to generate a forecast.
Because the factors are latent, they must be estimated.
The rolling window-based sample analog of the latent factor for the last $R$ dates and the associated factor loadings are therefore obtained by minimizing the sum of squared idiosyncratic error, $\sum_{i=1}^{N}\sum_{t=T-R+1}^{T}\left(x_{it}-\bm\lambda_i(1)^\top\bm f_t\right)^2$, under the restriction $\frac{1}{T}\sum_{t=T-R+1}^{T}\bm f_t\bm f_t^\top=\bm I_r$. The solution of the optimization problem, known as the principal component estimator, is given by $\tilde{\bm F}_R=[\tilde{\bm f}_{T-R+1}, \ldots,\tilde{\bm f}_T]^\top: R\times r$ corresponding to the $r$ largest eigenvectors of $\frac{1}{R}\bm X_R^\top \bm X_R,$ where $\bm X_R=(x_{it})_{i=1,\ldots, N, t=T-R+1,\ldots, T}$. The estimated factor loadings are obtained by regressing the vector corresponding to each row of $\bm X_R$ on the estimated factors. Given the restriction on the estimated factors, the estimated $N\times r$-matrix of factor loadings is given by $\tilde{\bm \Lambda}_R(1)=[\tilde{\bm \lambda}_{1, R}(1), \ldots,\tilde{\bm \lambda}_{N, R}(1)]^\top=\frac{1}{R}\bm X_R \tilde{\bm F}_R$.\footnote{The estimation approach could be extended to allow for kernel smoothing. In that case, we would need to choose a kernel and address boundary issues. See, for example, \citet*{SW2017}, who developed local least-squares estimation for factor and factor loadings over time, for more details. We focus on our objective of estimating the factor at the end of the sample using a constant weight over a window that must be selected appropriately.}
 \citet{Bai2003} show that estimated factors converge to a rotation of the latent factor, where the rotation matrix is given by
 $\tilde{\bm V}_T^{-1}\frac{1}{T}\tilde{\bm F}_T\bm F_T\frac{1}{N}\bm \Lambda(1)^\top\bm \Lambda(1)$, when the factor loadings in the underlying factor model are time-invariant, where $\tilde{\bm V}_{R}$ is an $r\times r$-diagonal matrix, with the largest eigenvalues on its diagonal in decreasing order using a window of size $R$, and $\bm\Lambda(t/T)$ is the $N\times r$-matrix of factor loadings at time $t$.
 We extend that result to the time-varying setting. As we will see later $\tilde{\bm f}_{T,R}=\bm H_R\bm f_T+o_P(1)$,
where $\bm H_R=\tilde{\bm V}_R^{-1}\frac{1}{R}\sum_{t=T-R+1}^T\tilde{\bm f}_t\bm f_t^\top\frac{1}{N}\bm \Lambda(t/T)\bm \Lambda(1)$
is the new rotation matrix, which is unknown to the econometrician. We also show in the appendix that $\frac{1}{R}\sum_{t=T-R+1}^T\Vert \tilde{\bm f}_{t,R}-\bm H_{R}\bm f_t\Vert^2
=O_P(\delta_{NRT}^{-2})$,  which extends the finding in \citet{BN2002} that the order is
 $
O_P(\delta_{NT}^{-2})$ to the time-varying setting, where $\delta_{NRT}^2=\min[N,R,\frac{T^2}{R^2}]$ and $\delta_{NT}^2=\min[N,T]$.


Let $\bm\delta_R(1)=(\bm\alpha^\top(1)\bm H_R^{-1}, \bm\beta^\top(1))^\top$, $\hat{\bm z}_{t,R}=(\tilde{\bm f}_{t,R}^\top, \ \bm w_t^\top)^\top$.
The rolling window based ordinary least-square estimator of $\bm\delta_R(1)$, i.e.,  $\bm\delta(t/T)$ when $t=T$, is obtained by regressing $y_{t+h}$ on $\hat{\bm z}_{t,R}$, for $t=T-R+1,\ldots,T,$ and given by
\begin{equation}
	\hat{\bm\delta}_R(1)=\bigg(\sum_{t=T-R+1}^{T-h}\hat{\bm z}_{t,R}\hat{\bm z}_{t,R}^\top\bigg)^{-1}\sum_{t=T-R+1}^{T-h}\hat{\bm z}_{t,R} y_{t+h}.
	\label{OLS}
\end{equation}
The forecast $\hat{y}_{T+h,R}$ of $y_{T+h}$ depends on the choice of the window $R$ and is given by
$	\hat{y}_{T+h,R}=\hat{\bm\delta}_R(1)^\top\hat{\bm z}_{T,R}.$
We can write the forecast error $y_{T+h}-\hat y_{T+h,R}$ as
\begin{equation}
	\varepsilon_{T+h}-(\hat{\bm\delta}^\top_R(1)\hat{\bm z}_T-\bm\delta^\top(1)\bm z_T)=\varepsilon_{T+h}-(\hat{\bm\delta}_R(1)-\bm\delta_R(1))^\top \hat{\bm z}_T  -\bm\alpha^\top(1)\bm H_R^{-1}(\tilde{\bm f}_T-\bm H_R \bm f_T),
	\label{FAR_V2}
\end{equation}
where $\bm\delta(1)=(\bm\alpha(1)^\top, \bm\beta(1)^\top)^\top$ and $\bm z_T=(\bm f_T^\top, \bm w_T^\top)^\top$.
The first term is the $ h$-step-ahead error term. The second component captures parameter-estimation uncertainty. Its order in probability is driven by the order of $\hat{\bm\delta}_R(1)-\bm\delta_R(1)=O_P\big(\frac{1}{\sqrt R}+\frac{R}{T}\big)$ as shown in the next section.
The third term is specific to the factor-augmented regression model and arises because the factors are not observed but estimated. It is the factor estimation uncertainty component in the forecast error.
As we show below $\tilde{\bm f}_{R,T}-\bm H_R\bm f_T=O_P\big(\frac{1}{\sqrt N}+\frac{1}{\sqrt R\delta_{NRT}}\big)$. This third quantity can dominate the second one, and highlights the need to consider the impact of factor estimation when selecting the optimal rolling window.

To find an estimate $\hat R$ for $R$, we aim to use a procedure which minimizes the conditional MSFE, $\mathrm{E}((y_{T+h}-\hat y_{T+h,R})^2|\mathcal{F}_T)$, where $\mathcal{F}_T$ is the information set up to time $T$.
We write
\begin{align}
    \mathrm{E}((y_{T+h}-\hat y_{T+h,R})^2|\mathcal{F}_T)=&\mathrm{E}\left(\varepsilon_{T+h}^2|\mathcal{F}_T\right)+
    \big(\hat{\bm\delta}^\top_R(1)\hat{\bm z}_{T,R}-\bm\delta^\top(1)\bm z_T\big)^2.
	\label{MSFE}
\end{align}
Because the variance of the $h$-period ahead forecast error $\mathrm{E}(\varepsilon_{T+h}^2|\mathcal{F}_T)=\sigma_T^2$ is constant and does not depend on the choice of $R$, it can be ignored when searching for $R$. Therefore, we focus on the infeasible loss function that changes with $R$,
\begin{align}
	C_{NT}^0(R)=& \big(\hat{\bm\delta}^\top_R(1)\hat{\bm z}_{T,R}-\bm\delta^\top(1)\bm z_T\big)^2.
	\label{MSFE}
\end{align}
To obtain the sample counterpart for this unobserved quantity, we extend the approach of \citet*{IJR2017}, who only considers structural change in regression coefficients and observed regressors, to the estimated factor setting.\footnote{It is possible that when a large break occurs, the number of selected factors, using existing factor selection procedures, changes before the terminal date $T$. Our framework could be easily extended to accommodate such a situation by selecting the number of factors associated with each candidate window. However, we assume the number of factors to be fixed to keep the theoretical development in the appendix tractable.}



Let $R_0$ denote a given pilot window size such that $R_0 \gg R$. We rewrite $C_{NT}^0(R)=(\hat{\bm z}_{T,R}^\top(1)\hat{\bm\delta}_R(1)-\bm z_{T}^\top\bm\delta(1))^2$ as
\begin{align}
&\big(\hat{\bm z}_{T,R}^\top(1)\hat{\bm\delta}_R(1)-\big[(\bm H^{0}\bm f_{T})^\top(((\bm H^{0})^{-1})^\top\bm\alpha(1))+\bm w_{T}^\top\bm\beta(1)\big]\big)^2,\nonumber
\end{align}
where $\bm H^{0}$ is a rotation-associated matrix for pilot-window-based estimation.
To get a feasible version $C_{NT}(R)$ of $C_{NT}^0(R)$, local linear estimators
$\ddot{\bm \alpha}(1), \ddot{\bm \beta}(1)\text{   and    }\ddot{\bm f}_{T},$
are used to replace
$((\bm H^{0})^{-1})^\top\bm\alpha(1), \bm \beta(1)\text{   and   }\bm H^{0}\bm f_T.$  Denote $\tilde{\bm f}_{t}^0, t=T-R_0+1,\ldots,T$, the PC estimated factors based on the last $R_0$ observations.
 The proposed procedure minimizes a feasible version of $C_{NT}^0$, where sample estimators $\ddot{\bm \delta}(1)$ and  $\ddot{\bm f}_{T}$ replace a rotated version
of the unknown parameters
and the latent factors
at time $T$.
The criterion
 is therefore given by
\begin{align}
	C_{NT}(R)=
    &\big(\hat{\bm z}_{T,R}^\top(1)\hat{\bm\delta}_R(1)-\ddot{\bm z}_{T}^\top(1)\ddot{\bm\delta}(1)\big)^2,
	\nonumber
\end{align}
where $\ddot{\bm z}_{T}=(\ddot{\bm f}_{T}^\top,\bm w_T^\top)^\top $.

To find $\ddot{\bm f}_{T}$, we rely on the following local expansion for $\lambda_i(t/T)$ around one,
\begin{equation}
\bm\lambda_i\big(\frac{t}{T}\big)=\bm\lambda_i(1)+\bm\lambda_i^\prime(1)\big(\frac{t}{T}-1\big)+\bm\lambda_i^{\prime\prime}(c_i)\big(\frac{t}{T}-1\big)^2,
\label{loading_expansion}
\end{equation}
where $c_i\in (t/T,1)$ for $t\leq T$, and $\bm\lambda_i^\prime(\cdot)$ and $\bm\lambda_i^{\prime\prime}(\cdot)$ are first and second derivatives of $\bm\lambda_i(\cdot)$. Replacing $\bm\lambda_i\big(\frac{t}{T}\big)$ in the panel factor model and using the identity $\bm H^0\bm f_t=\tilde{\bm f}_{t}^0-(\tilde{\bm f}_{t}^0-\bm H^0\bm f_t)$, we get for $t=T-R_0+1,\ldots,T,$ that
\begin{equation}
x_{it}=(((\bm H^0)^{\top})^{-1}(1)\bm\lambda_i(1))^\top\tilde{\bm f}_{t}^0+((\bm H^0)^{\top -1}\bm\lambda_i^\prime(1))^\top\big(\frac{t}{T}-1\big)\tilde{\bm f}_{t}^0+u_{it},
\label{panel_expansion}
\end{equation}
where $u_{it}$ is the new error term.

Now, we can obtain an alternative estimator $\ddot{\bm \lambda}_{i}(1)$ for $\bm \lambda_{i}(1)$, by regressing $x_{it}$
on $\tilde{\bm f}_{t}^0 \text{ and } \big(t/T-1\big)\tilde{\bm f}_{t}^0, t=T-R_0+1,\ldots,T.$ Thus, the computed $\ddot{\bm \lambda}_{i}$ corresponds to the estimator of the factor loadings associated with $\ddot{\bm f}_T$. To find the alternative estimator $\ddot{\bm f}_T$ of $\bm H^{0}\bm f_T$,
we note that $\bm H^{0}\bm f_T$
equals
$$(\tilde{\bm V}^{0})^{-1}(R^{-1}\sum_{t=T-R+1}^T\tilde{\bm f}_{t}^0(\bm H^{0}\bm f_t)^\top )(N^{-1}\sum_{i=1}^N((\bm H^{0})^{-1 \top}\bm \lambda_i(t/T)) ((\bm H^{0})^{-1 \top}\bm \lambda_i)^\top)(\bm H^{0}\bm f_T),$$
where $\tilde{\bm V}^{0}$ is the matrix of eigenvalues from the PCM in decreasing order.
 To get $\ddot{\bm f}_{T}$, the alternative estimator of $\bm H^0\bm f_{T}$, we use the formula
$$
(\tilde{\bm V}^{0})^{-1}\frac{1}{R_0}\sum_{t=T-R_0+1}^T(\tilde{\bm f}_{t}^0) (\tilde{\bm f}_{t}^0)^{\top}\bigg(
\frac{1}{N}\sum_{i=1}^N[\ddot{\bm\lambda}_i(1)+\frac{(t-T)}{T}\ddot{\bm\lambda}_i^\prime(1)]\ddot{\bm\lambda}_i^\top(1)\bigg)
\tilde{\bm f}_{T}^0,$$
replacing $\bm H^{0}\bm f_t$ with $\tilde{\bm f}_{t}^0$ and $((\bm H^{0})^{-1})^{\top}\bm \lambda_i(t/T)$ with $\ddot{\bm \lambda}_{i}(1)+\frac{t-T}{T}\ddot{\bm \lambda}_{i}^\prime(1)$.


To obtain $\ddot{\bm \delta}(1)$, we follow \citet*{IJR2017} for the case without latent factors, and use a local linear expansion of $\bm\delta(\cdot)$, which is given by
\begin{equation}
\bm\delta\bigg(\frac{t}{T}\bigg)=\bm\delta(1)+\bm\delta^\prime(1)\bigg(\frac{t}{T}-1\bigg)+\bm\delta^{\prime\prime}(d)\bigg(\frac{t}{T}-1\bigg)^2,
\label{expansion}
\end{equation}
where $d\in (t/T,1)$ for $t\leq T$, and $\bm\delta^\prime(\cdot)$ and $\bm\delta^{\prime\prime}(\cdot)$ are first and second derivatives of $\bm\delta(\cdot)$.

Define $\hat{\bm z}_{t}^0=\begin{pmatrix} (\tilde{\bm f}_{t}^0)^\top, \bm w_t^\top \end{pmatrix}^\top, $ as well as the rotated parameters
$$  \bm\delta^0\big(\frac{t}{T}\big)=
 \begin{pmatrix}
  ((\bm H^0)^{ -1})^\top\bm\alpha\big(\frac{t}{T}\big) \\
  \bm\beta\big(\frac{t}{T}\big) \end{pmatrix}, \ \
  (\bm\delta^0)^\prime\big(\frac{t}{T}\big)=
 \begin{pmatrix}
  ((\bm H^0)^{ -1})^\top\bm\alpha^\prime\big(\frac{t}{T}\big) \\
  \bm\beta^\prime\big(\frac{t}{T}\big) \end{pmatrix}.$$
Using the expansion \Cref{expansion},
we rewrite \Cref{FAR} when $t$ is close to $T$ as
\begin{align}
y_{t+h}=(\bm\delta(1))^\top \hat{\bm z}_{t}^0+(\bm\delta^{\prime}(1))^\top \big(\frac{t}{T}-1\big) \hat{\bm z}_{t}^0+v_{t+h},
\label{expansion_y}
\end{align}
where $v_{t+h}$ is the corresponding error term.
By regressing $y_{t+h}$ on $\hat{\bm z}_t^0$ and $\frac{t-T}{T}\hat{\bm z}_t^0, t=T-R_0+1,\ldots, T-h$, we get $\ddot{\bm\delta}(1)$, which consists of the estimated coefficients for $\hat{\bm z}_t^0$. Hence, the first $r$ elements of $\ddot{\bm\delta}(1)$ provide $\ddot{\bm\alpha}(1)$.
Next, we present the convergence rate of the estimators in our new framework as well as the optimality result.

\section{Main Results}
\label{Main}

Before presenting the main results, we state the assumptions for any rolling window $R$.

\begin{assumption}
\label{Ass1}(Factor model and idiosyncratic errors)
\end{assumption}

\begin{description}
\item[(a)] $\mathrm{E}\left\Vert \bm{f}_{t}\right\Vert ^{8}\leq C$ and
$\bm{\Sigma}_{\bm{F}}=\mathrm{E}(\bm{f}_{t}\bm{f}_{t}^\top)=\bm I_r.$

\item[(b)] The factor loadings $\left\{  \bm{\lambda}_{i}(t/T)\right\}
_{i=1,\ldots,N, \ t=1,\ldots,T}$ are non-random and bounded. In addition,
$\frac{1}{N}
\sum_{i=1}^N\bm{\lambda}_{i}(1)\bm{\lambda}_{i}^{\top}(1)  =\bm{\Sigma
}_{\bm{\Lambda}}+o(1)$, where $\bm{\Sigma }_{\bm{\Lambda}}$ has eigenvalues that are finite and bounded away from zero.

\item[(c)] The eigenvalues of the $r\times r$ matrix
$\bm{\Sigma }_{\bm{\Lambda}}  $
are distinct.

\item[(d)] $\mathrm{E}\left(  e_{it}\right)  =0,$ $\mathrm{E}\left\vert
e_{it}\right\vert ^{8}\leq C.$

\item[(e)] $\mathrm{E}\left(  e_{it}e_{js}\right)  =\sigma_{ij,ts}
$,\ $\left\vert \sigma_{ij,ts}\right\vert \leq\overline{\sigma}_{ij}$ for all
$\left(  t,s\right)  $ and $\left\vert \sigma_{ij,ts}\right\vert \leq\tau
_{st}$ for all $\left(  i,j\right)  $, with $\frac{1}{N}\sum_{i,j=1}
^{N}\overline{\sigma}_{ij}\leq C,$ $\frac{1}{R}\sum_{t,s=T-R+1}^{T}\tau_{st}\leq
C$ and $\frac{1}{NR}\sum_{i,j=1}^N\sum_{t,s=T-R+1}^T\left\vert \sigma_{ij,ts}\right\vert \leq
C.$

\item[(f)] $\mathrm{E}\left\vert \frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(
e_{it}e_{is}-E\left(  e_{it}e_{is}\right)  \right)  \right\vert ^{4}\leq C$
for all $\left(  t,s\right)  .$

\item[(g)] $\sum_{s=T-R+1}^T \big\vert \frac{1}{N}\sum_{i=1}^N \sigma_{ii,st} \big\vert\leq C$, for any $t$.
\end{description}

\begin{assumption}
\label{Ass2}(Moment conditions and weak dependence for $\{\bm{z}_{t}\}$
and
$\{e_{it}\}$)
\end{assumption}

\begin{description}
\item[(a)] $
\mathrm{E}\left(  \frac{1}{N}\sum_{i=1}^{N}\left\Vert \frac{T^k}{R^k}\frac
{1}{\sqrt{R}}\sum_{t=T-R+1}^{T}\bm{f}_{t}e_{it}\frac{(T-t)^k}{T^k}\right\Vert
^{2}\right)  \leq C$,  $k=0,1,2$, where $\mathrm{E}\left(  \bm{f}_{t}
e_{it}\right)  =\bm 0$.

\item[(b)] For each $t,$ $
\mathrm{E}\left\Vert \frac{1}{\sqrt{RN}}\sum
_{s=T-R+1}^{T-h}\sum_{i=1}^{N}\bm{z}_{s}\left(  e_{it}e_{is}-\mathrm{E}\left(
e_{it}e_{is}\right)  \right)  \right\Vert ^{2}\leq C$, where $\bm z_s=(\bm f_s^\top \ \bm w_s^\top)^\top$.

\item[(c)] $
\frac{1}{R^2}\sum_{t=T-R+1}^{T}\sum
_{s=T-R+}^{T-h} \mathrm E\big(\big\Vert\frac{1}{\sqrt N}\sum_{i=1}^{N}\bm f_s^\top\bm \lambda_i\big(\frac{s}{T}\big)\bm e_{it}\big\Vert^2\big)\leq C$, for all $s$.

\item[(d)] $
\mathrm{E}\left\Vert \frac{T^k}{R^k} \frac{1}{\sqrt{RN}}\sum_{t=T-R+1}^{T-h}\sum
_{i=1}^{N}\bm{z}_{t}\bm{\lambda}_{i}^{\top}(\frac{s}{T})e_{it} \frac{(T-t)^k}{T^k}
\right\Vert
^{2}\leq C$, $k=0,1,2$, for all $s$, where $\mathrm{E}(\bm{z}_{t}e_{it})  =\bm0$.

\item[(e)] $
\mathrm{E}\left(
\left\Vert \frac
{1}{\sqrt{N}}\sum_{i=1}^{N}\bm{\lambda}_{i}(s/T)e_{it}\right\Vert
^{4}\right)  \leq C$, for all $(t,s)$.

\item[(f)] $
\mathrm{E}\left(
\left\Vert \frac
{1}{\sqrt{N}}\sum_{i=1}^{N}\bm{\lambda}_{i}^\prime(1)e_{it}\right\Vert
^{2}\right)  \leq C$
and
 $
\mathrm{E}\left(
\left\Vert \frac
{1}{\sqrt{N}}\sum_{i=1}^{N}\bm{\lambda}_{i}^{\prime\prime}(c_i)e_{it}\right\Vert
^{2}\right)  \leq C$
for all $(t,s)$.
\end{description}


\begin{assumption}
	\label{Ass3}(Moment conditions and weak dependence for $\{\varepsilon_{t+h}\}$
    and
	$\{e_{it}\}$)
\end{assumption}

\begin{description}

	\item[(a)] For each $t,$ $
    \mathrm{E}\left\Vert \frac{1}{\sqrt{RN}}\sum
	_{s=T-R+1}^{T-h}\sum_{i=1}^{N}\varepsilon_{s+h}\left(  e_{it}e_{is}-\mathrm{E}\left(
	e_{it}e_{is}\right)  \right)  \right\Vert ^{2}\leq C$.

	\item[(b)] $
    \mathrm{E}\left\Vert \frac{T^k}{R^k}\frac{1}{\sqrt{RN}}\sum_{t=T-R+1}^{T-h}\sum
	_{i=1}^{N}\varepsilon_{t+h}\bm{\lambda}_{i}^{\top}(s/T)e_{it}\frac{(T-t)^k}{T^k}\right\Vert
	^{2}\leq C$, $k=0,1,2$, for all $s$, where $\mathrm{E}\left(  \varepsilon_{t+h}
    e_{it}\right)  =\bm0$.
\end{description}

\begin{assumption}
	\label{Ass4}(Moments and dependence conditions for $\{v_t\}=\{(\bm z_{t-h}^\top, \varepsilon_t)^\top\}$)
\end{assumption}

\begin{description}
\item[(a)] $\{\varepsilon_{t+h} \}_{t=1}^T$ is a sequence such that $\mathrm{E}\left(\varepsilon_{t+h} |\mathcal{F}_t\right)=\bm 0,$ $\mathrm{E}\left(\varepsilon_{t+h}^2 |\mathcal{F}_t\right)=\sigma_t^2,$ is well-defined almost surely, and  $\mathcal{F}_t=\sigma(\ldots,x_{1(t-1)},x_{1t} ,\ldots, x_{N(t-1)},x_{Nt},\ldots, \bm z_{t-1}^\top, \bm z_{t}^\top, \ldots, y_{t-1}, y_{t}),$ an information set up to time $t$.

\item[(b)] For $m>2$, the sequence $\{v_t\}$,
$t=h+1,\ldots,T+h$, is $L_{4m/(m-1)}$-NED of size $-2$, with positive constants $d_t=O(\Vert \bm v_t -\mathrm{E}(v_t)\Vert_{4m/(m-1)})$, on the $\alpha$-mixing sequence  $\{w_t \}_{-\infty}^\infty$ of size $-2m/(m-2)$, and  $\Vert \mathrm{Vec}(\bm v_t \bm v_t^\top) \Vert_m\leq C$.

\item[(c)]  $
\frac{1}{R}\sum_{t=T-R+1}^{T-h}
\mathrm(\bm{z}_{t}\bm{z}_{t}^{\top})\overset{P}{\longrightarrow}
\bm{\Sigma }_{\bm{Z}}>0$, where $\bm{\Sigma }_{\bm{Z}}=\mathrm{E}( \bm z_t  \bm z_t^\top)$ is non-random, $\mathrm{E}(\Vert \bm z_t\Vert^8 )\leq C$ and $\mathrm{E}(\vert \varepsilon_{t+h}\vert^8 )\leq C$.

\item[(d)]  All eigenvalues of $\mathrm{E}(\bm z_t\bm z_t^\top)$ and $\mathrm{E}(\bm z_t\bm z_t^\top\varepsilon_{t+h}^2)$ are finite and bounded away from zero uniformly in $t$.
\end{description}

\begin{assumption}
	\label{Ass5}(Smoothness conditions)
\end{assumption}

\begin{description}
\item[(a)]
$\bm \delta(\cdot)
$ is twice continuously differentiable over $\mathbb R$, and $\Vert \bm \delta^{\prime}(t/T) \Vert$ and  $\Vert \bm \delta^{\prime\prime}(t/T) \Vert$  are bounded uniformly in $t$.

\item[(b)]  $\bm \lambda(\cdot)$ is twice continuously differentiable over $\mathbb R$, and $\Vert\bm \lambda^{\prime}(t/T)\Vert$ and  $\Vert\bm \lambda^{\prime\prime}(t/T)\Vert$  are bounded uniformly in $t$.

\end{description}

\Cref{Ass1} (a)-(f) adapts Assumptions 1 and 2 on the factor model and the idiosyncratic error of \citet{DGP2015} to the time-varying setting. The key difference in \Cref{Ass1} (a)-(c) comes from the normalization of the latent factor moment to an identity matrix similar to Assumption 1 (a) in \citet{GMP2017}. \Cref{Ass1} (g) is the same as in the same as Assumption E (1) on weak dependence in \citet{Bai2003}.  \Cref{Ass2,Ass3} require weak dependence among $\bm z_t, e_{it}$, and $\varepsilon_{t+h}$.
They modify the standard assumptions in \citet{GP2014} to accommodate local expansion around the terminal date in our derivation of the asymptotic results. \Cref{Ass5} contains conditions that resemble Assumptions 1-3 used by \citet{IJR2017} for the error term and the regressors in the time-varying forecasting equation.
\Cref{Ass5} states that the factor loadings and the regression coefficients are twice differentiable and the derivatives are bounded. The next theorem presents our first key results.


\begin{theorem}[Order from minimizing $C_{NT}^0(R)$]
If Assumptions \Cref{Ass1,Ass2,Ass3,Ass4,Ass5}  hold, $T,$ $
N\rightarrow \infty $, $\frac{\sqrt{N}}{R}\rightarrow 0$, and $\frac{\sqrt{R}}{N}\rightarrow c<\infty,$ then
\begin{align}
&\hat{\bm{\delta}}_R(1)-\bm{\delta}_R(1)=O_P\bigg(\frac{1}{\sqrt{R}}+\frac{1}{N}+\frac{R}{T}\bigg)=O_P\bigg(\frac{1}{\sqrt{R}}+\frac{R}{T}\bigg),
	\label{Thm11}\\
&\tilde{\bm f}_{T,R}-\bm H_R\bm f_T=O_P\bigg(\frac{1}{\sqrt{N}}+\frac{1}{\sqrt{R}\delta_{NRT}}\bigg),
	\label{Thm12}\\
& \hat{\bm z}_T^\top\hat{\bm{\delta}}_R(1)-\bm z_T^\top\delta(1)=O_P\bigg(\frac{1}{\sqrt R}+\frac{1}{\sqrt{N}}+\frac{R}{T}\bigg). \label{Thm13}
\end{align}
In  addition, the value of $R$ that minimizes
\begin{equation}
C_{NR}^0(R)=O_P\bigg(\frac{1}{R}+\frac{1}{N}+\frac{R^2}{T^2}\bigg).
\label{Thm14}
\end{equation}
is of order $T^{\frac{2}{3}}$
if $c=0$, and is of order $T^{\frac{4}{5}}$ otherwise, in probability.
	\label{Thm1}
\end{theorem}

\Cref{Thm1} \hyperref[{Thm11}]{\textup{\tagform@{\ref*{Thm11}}}} and \hyperref[{Thm12}]{\textup{\tagform@{\ref*{Thm12}}}} extend the order in probabilities for the asymptotic results from \citet{BN2006} and \citet{GP2014} to the time-varying context. In \Cref{Thm1} \hyperref[{Thm11}]{\textup{\tagform@{\ref*{Thm11}}}}, we show that the convergence rate of the rolling window-based estimator of the regression coefficients in the presence of estimated factors matches that in \citet{IJR2017} when the factors are observed. When the window size relative to the sample size is small enough, the convergence rate becomes $\sqrt R$. \Cref{Thm1} \hyperref[{Thm12}]{\textup{\tagform@{\ref*{Thm12}}}} is new and extends the $\sqrt N$ convergence rate from \citet{Bai2003} for the estimated factor at the end of the sample. It shows that convergence can be slower because the factor loadings vary over time. Combining these findings, we show in \Cref{Thm1} \hyperref[{Thm13}]{\textup{\tagform@{\ref*{Thm13}}}} that the number of variables $N$ can also affect the accuracy of the conditional mean estimate given information up to time $T$. \Cref{Thm1} \hyperref[{Thm14}]{\textup{\tagform@{\ref*{Thm14}}}} shows that when the window size is relatively small compared to $N$, such that $\sqrt R/N\rightarrow 0$, its value that minimizes the infeasible criterion is of order $T^{2/3}$, which is the same as in \citet{IJR2017}. When  $\sqrt R/N\rightarrow c$, such that $0<c<\infty$, the window rate that minimizes the infeasible criterion becomes $T^{4/5}$. To obtain the results related to the local expansion, we make the following additional assumptions.

\begin{assumption}
	\label{Ass6}(Uniform order in probability)
\end{assumption}


\begin{description}
\item[(a)] $\sup_{t,s}
\left(  \frac{1}{R}\sum_{t=T-R+1}^{T-h}\left\Vert \frac
{1}{\sqrt{N}}\sum_{i=1}^{N}\bm{\lambda}_{i}(s/T)e_{it}\right\Vert
^{4}\right) =O_P(1)$.
\item[(b)] $
   \sup_{t,s}\left\vert \frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(
e_{it}e_{is}-E\left(  e_{it}e_{is}\right)  \right)  \right\vert
=O_P(1)$.
\item[(c)] $\sup_R\frac{1}{R}\sum_{t=T-R+1}^{T}\Vert\frac{1}{\sqrt{RN}}\sum_{s=T-R+1}^{T-h}(\bm e_s^\top\bm e_t-\mathrm{E}(\bm e_s^\top\bm e_t))\bm z_s\Vert^2=O_P(1)$.
\item[(d)] $\sup_R\frac{1}{R}\sum_{t=T-R+1}^{T}\Vert\frac{1}{\sqrt{RN}}\sum_{s=T-R+1}^{T-h}(\bm e_s^\top\bm e_t-\mathrm{E}(\bm e_s^\top\bm e_t))\varepsilon_{s+h}\Vert^2=O_P(1)$.
\end{description}

\begin{assumption}
	\label{Ass7}(Conditions on $R$ and $R_0$)
\end{assumption}

\begin{description}
	\item[(a)] When $T\rightarrow \infty$, $R_0\rightarrow \infty$, $R\rightarrow \infty$, $R/\min(R_0,N) \rightarrow 0$, and $\frac{NR_0^2}{T^2}\rightarrow 0$

	\item[(b)] The cardinality of $\Theta_R \subseteq [\underline{R}, \  \overline{R}]\subset \mathbb{Z}^+$ satisfies $\underline{R}^\rho,$ for some $\rho\in (0,1)$.
	The conditions on $R$ in \Cref{Ass7} (a) are satisfied by $\underline{R}$ and $\overline{R}$.
\end{description}

In \Cref{Ass6}, we strengthen \Cref{Ass2} and \Cref{Ass3} to establish the optimality of the new procedure. In \Cref{Ass7}, we outline the conditions on the rolling window and the pilot window needed for implementation. \Cref{Ass7}
(b) is the same as in \citet{IJR2017}, while  \Cref{Ass7} (a)  is different and assumes the window to be smaller than the number of variables.

\begin{theorem}[Order for local linear Estimations]
If \Cref{Ass1,Ass2,Ass3,Ass4,Ass5} hold, $T,$ $
N\rightarrow \infty $, $\frac{\sqrt{N}}{R}\rightarrow 0$, and $\frac{\sqrt{R_0}}{N}\rightarrow c<\infty,$ then
\begin{align}
&\ddot{\bm{\delta}}(1)-\bm{\delta}^0(1)=O_P\bigg(\frac{1}{\sqrt{R_0}}+\frac{1}{N}+\frac{R_0}{T}
\bigg)=O_P\bigg(\frac{1}{\sqrt{R_0}}+\frac{R_0}{T}
\bigg),\label{Thm21}\\
&\ddot{\bm f}_{T}-\bm H^0\bm f_T=O_P\bigg(\frac{1}{\sqrt{N}}+\frac{1}{R_0}+\frac{R_0}{T}
\bigg)=O_P\bigg(\frac{1}{\sqrt{N}}+\frac{R_0}{T}
\bigg),\label{Thm23}\\
&
\ddot{\bm z}_T^\top\ddot{\bm{\delta}}(1)-\bm z_T^\top\bm{\delta}(1)=O_P\bigg(\frac{1}{\sqrt{R_0}}+\frac{1}{\sqrt{N}}+\frac{R_0}{T}
\bigg)\label{Thm24}
\end{align}
\label{Thm2}
\end{theorem}

In \Cref{Thm2} \hyperref[{Thm21}]{\textup{\tagform@{\ref*{Thm21}}}}, we find that the regression coefficient estimator based on local linear estimation is $O_P\big(\frac{1}{\sqrt{R_0}}+\frac{R_0}{T}
\big)$,
\Cref{Thm2} \hyperref[{Thm21}]{\textup{\tagform@{\ref*{Thm21}}}} shows that the order in probability of the factor estimator is inflated by the pilot window relative to the full sample, compared to the standard asymptotic case. \Cref{Thm2} \hyperref[{Thm23}]{\textup{\tagform@{\ref*{Thm23}}}} provides the resulting order for the conditional mean estimation. Our key result appears below.

\begin{theorem}[Optimality]
If Assumptions \Cref{Ass1,Ass2,Ass3,Ass4,Ass5,Ass6,Ass7} hold, $T,$ $
N\rightarrow \infty $, $\frac{\sqrt{N}}{R}\rightarrow 0$, and $\frac{\sqrt{R_0}}{N}\rightarrow c<\infty,$ then
\begin{equation}
\frac{C_{NT}\left( \hat{R}\right) }{C_{NT}^0\left(\hat R^0\right) }\overset{P}{\longrightarrow} 1,
\label{Thm3}
\end{equation}
where
\begin{equation*}
\hat{R}=\arg \min_{R\in \Theta_R}C_{NT}\left( R\right)\text{ \  and \ } \hat{R}^0=\arg \min_{R\in \Theta_R}C_{NT}^0\left( R\right).
\end{equation*}
\label{Thm3}
\end{theorem}
\Cref{Thm3} states that the proposed method is asymptotically loss efficient. It means that minimizing the criterion $C_{NT}(R)$ is equivalent to minimizing the infeasible loss function  $C_{NT}^0(R)$.
We provide simulation evidence of the performance of the proposed procedure in the next section.


\section{Simulation Experiments}
\label{Simulation}
These simulation experiments compare the procedure for selecting the window size to generate forecasts when factor-estimation uncertainty is accounted for and when it is not. Second, we show how the proposed criterion improves window-size estimation when structural instabilities arise in the loadings of the latent factor model and in the slope parameters of the forecasting equation.

To perform the simulation exercise, we consider a factor-augmented regression
that allows time-varying regression coefficients and factor loadings,
and consider several data generating processes (DGPs). The factor-augmented regression model is given by
$y_{t+1}=\alpha(t/T)f_t+\beta(t/T)+\varepsilon_{t+1}, t=0,\ldots,T,$
where $f_{t+1}=0.9f_t+u_{t+1}$,
the initial value
 and the error terms
 $\varepsilon_{t+1}$ and $u_{t+1}$ are $\mathrm{IIDN}(0, 1)$.
 The latent factors $f_t, t=1,\ldots,T$ summarize the comovement in the factor panel model
 $x_{it}=\lambda_i(t/T)f_t+e_{it}, t=1,\ldots, T, i=1,\ldots,N,$
where $\lambda_i(t/T)$ and $e_{it}$ are defined as follows.

The first DGP, called DGP 1, assumes $\alpha(t/T)=1, \beta(t/T)=1,$ and the factor loadings $\lambda_i(t/T)\sim \mathrm{IIDU}[0,1]$,
where the mean is positive to avoid zero exposure of the variables $x_{it}$  to the latent factor, and we let
the idiosyncratic errors  $e_{it}=\sigma_{i}v_{it}$, with $\sigma_{i}\sim \mathrm{IIDU}[0,1]$ and $v_{it}\sim \mathrm{IIDN}(0,1)$, to allow for heroskedasticity in the idiosyncratic errors.
In the second specification called DGP 2, we keep $\alpha(t/T)=1, \beta(t/T)=1,$  $\lambda_i\sim\mathrm{IIDU}[0,1]$, but let $\bm e_t=(e_{1t},\ldots,e_{Nt})^\top \sim \mathrm{IIDN}(\bm 0_{N\times 1},\bm \Sigma_{\bm e})$, where $\bm \Sigma_{\bm e}$ is an $N\times N$ Toeplitz matrix with $(i,j)$-element given by $0.5^{|i-j|}$, similar to \citet{GP2014}, to introduce correlation among the cross-sectional units at each date. In DGPs 3 and 4, we modify DGPs 1 and 2 to allow for structural change in the regression coefficients and the factor loadings. In particular, we set $\lambda_{i}(t/T)=\lambda_{i0}+2\mathbb{I}(t\geq 0.5 T)$, where $\lambda_{i0}\sim\mathrm{IIDU}[0,1]$, $\alpha=1+2\mathbb{I}(t\geq 0.5T)$, and $\beta=1+4\mathbb{I}(t\geq 0.5T)$.
To allow for smooth structural changes in the loadings, we let $\lambda_i(t/T)=\lambda_{i0}+2t/T$, with $\lambda_{i0}\sim\mathrm{IIDU}[0,1]$,
 $\alpha(t/T)=1+2t/T$, and $\beta(t/T)=1+4t/T$,
in DGPs 1 and 2, to have DGPs 5 and 6, respectively. For the analysis, we set $T=100,200,300,$, $N=100,200,300$,
and simulate the data $M=1,000$ times. In each iteration of the simulation experiment, we construct the large dataset $\bm X=[x_{it}]: N\times T$.

Next, we select a pilot rolling window using the cross-validation approach described in \citet{PT2007} and considered by \citet{IJR2017}. We set the minimum rolling window, $R_{\min}$, to the closest integer to $\max(0.10 T, T^{4/5})$, and employ $0.25T$ periods for the evaluation step in the cross-validation.
To select the number of factors, we apply the $IC_{p2}$ criterion proposed by \citet{BN2002} on recent data based on the pilot window size. Moreover, we compute the root mean square forecast error (RMSFE) under different rolling-window choices, relative to the RMSFE of the forecast based on the full sample.

To do so, we first use $\bm X$ to estimate the latent factors
and the factor loadings using PCM based on the full sample. We then obtain the least-squares estimator and forecast $y_{T+1}$ using the estimated factors at time $T$.
In addition, we construct the matrix $\bm X_R=[x_{it}]_{i=1,\ldots,N, t=T-R+1,\ldots,T}$, where $R_{min}\leq R\leq T$.  Therefore, we estimate the factors and factor loadings from the matrix $\bm X_R$. For each criterion discussed previously, we select the window $R$ and generate the forecasts. Furthermore, we select the rolling window using the approach for observed regressors (method $j=1$) in \citet{IJR2017}, the new selection procedure (method $j=2$), and the selection based on the minimization of the infeasible loss function (method $ j=3$).\footnote{Note that one can also test for structural instabilities before applying the rolling window selection procedures as in \citet{IJR2017} to exclude situations where the rolling window selection is irrelevant.} We call the forecast based on the full sample $j=0$.

For each simulated data, $m=1,\ldots,M,$ let $\hat y_{T+1}^{(m,j)}$ denote the forecasted values by the method $j=0,1,2,3$,
and the simulated values of $y_{T+1}^{(m)}$. The RMSFE for method $j$ is given by $\sqrt{\frac{1}{M}\sum_{m=1}^{M}(\hat y_{T+1}^{(m,j)}-y_{T+1}^{(m)})^2}$. The results below report the relative RMSFE for methods $j=1,2$ and $3$,
which is the RMSFE for each of these methods divided by the RMSFE based on the full-sample forecast.
The results are presented in
\Cref{tab1,tab2,tab3} for each value of $N$.




\begin{table}[h!]
\begin{center}
\caption{Relative RMSFE when $N=100$}
\label{tab1}
\begin{tabular}{lcccccc}
\hline\hline
$T=100$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.01 & 1.00 & 0.59 & 0.56 & 0.69 & 0.70 \\
Method 2 & 1.02 & 1.02 & 0.46 & 0.47 & 0.59 & 0.58 \\
Method 3 & 0.98 & 0.97 & 0.44 & 0.45 & 0.58 & 0.58 \\

\hline
$T=200$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.01 & 1.01 & 0.53 & 0.51 & 0.63 & 0.64 \\
Method 2 & 1.01 & 1.01 & 0.48 & 0.46 & 0.54 & 0.57 \\
Method 3 & 0.98 & 0.98 & 0.45 & 0.45 & 0.54 & 0.56 \\

\hline
$T=300$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.00 & 1.01 & 0.51 & 0.52 & 0.66 & 0.66 \\
Method 2 & 1.00 & 1.02 & 0.45 & 0.46 & 0.54 & 0.54 \\
Method 3 & 0.99 & 0.98 & 0.43 & 0.44 & 0.54 & 0.54 \\

\hline\hline
\end{tabular}
\end{center}
\textit{Note:} This table reports the MSFEs based on minimizing $C_{NT}(R)$ (called Method 2), $C_{NT}^0(R)$ (called Method 3), and the feasible version that ignores factor estimation (called Method 1),
relative to the MSFE for the full sample forecast.
\end{table}



\begin{table}[h!]
\begin{center}
\caption{Relative RMSFE when $N=200$}
\label{tab2}
\begin{tabular}{lcccccc}
\hline\hline
$T=100$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.00 & 1.00 & 0.58 & 0.60 & 0.69 & 0.70 \\
Method 2 & 1.01 & 1.01 & 0.48 & 0.47 & 0.58 & 0.59 \\
Method 3 & 0.98 & 0.98 & 0.46 & 0.44 & 0.57 & 0.59 \\

\hline
$T=200$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.02 & 1.00 & 0.52 & 0.51 & 0.65 & 0.66 \\
Method 2 & 1.02 & 1.01 & 0.47 & 0.46 & 0.55 & 0.56 \\
Method 3 & 0.99 & 0.99 & 0.45 & 0.44 & 0.55 & 0.55 \\

\hline
$T=300$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.00 & 1.01 & 0.56 & 0.52 & 0.67 & 0.65 \\
Method 2 & 1.01 & 1.01 & 0.49 & 0.45 & 0.55 & 0.54 \\
Method 3 & 0.99 & 0.99 & 0.47 & 0.44 & 0.54 & 0.54 \\

\hline\hline
\end{tabular}
\newline
\textit{Note:} See note for \Cref{tab1}.
\end{center}
\end{table}




\begin{table}[h!]
\begin{center}
\caption{Relative RMSFE when $N=300$}
\label{tab3}
\begin{tabular}{lcccccc}
\hline\hline
$T=100$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.01 & 1.01 & 0.59 & 0.60 & 0.69 & 0.71 \\
Method 2 & 1.03 & 1.02 & 0.46 & 0.47 & 0.58 & 0.59 \\
Method 3 & 0.99 & 0.99 & 0.43 & 0.44 & 0.57 & 0.58 \\

\hline
$T=200$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.00 & 1.00 & 0.51 & 0.50 & 0.66 & 0.66 \\
Method 2 & 1.00 & 1.01 & 0.46 & 0.45 & 0.57 & 0.55 \\
Method 3 & 0.99 & 0.98 & 0.44 & 0.44 & 0.57 & 0.55 \\

\hline
$T=300$  & DGP 1 & DGP 2 & DGP 3 & DGP 4 & DGP 5 & DGP 6 \\ \hline
Method 1 & 1.01 & 1.00 & 0.52 & 0.51 & 0.65 & 0.66 \\
Method 2 & 1.01 & 1.01 & 0.45 & 0.45 & 0.54 & 0.54 \\
Method 3 & 0.99 & 0.99 & 0.43 & 0.44 & 0.54 & 0.54 \\
\hline\hline
\end{tabular}
\newline
\textit{Note:} See note for \Cref{tab1}.
\end{center}
\end{table}



\clearpage
 \Cref{tab1} shows the results for the six DGPs when the number of variables is $N=100$. The table contains three panels for the different values of $T$. Regardless of the method and the sample size, the results for DGPs 1 and 2 show that the RMSFEs are close to one when there is no structural instability. These findings suggest that selecting
the rolling window yields the same results as forecasting based on the full sample. When there is a structural break (DGPs 3 and 4) or a smooth structural instability (DGPs 5 and 6), we note that the mean squared forecast error can be reduced by at least $30\%$. The most important improvement occurs when the gain in accuracy comes from minimizing the infeasible loss function, $C_{NT}^0(R)$. However, the econometrician does not know it, and it cannot be used in practice. The closest in performance is the proposed rolling window selection procedure, which extends the one proposed by \citet{IJR2017} to accommodate factor estimation uncertainty.  The \citet{IJR2017} approach performs reasonably well. Nevertheless, accounting for the fact that the factors are latent and must be estimated yields reduced mean-square forecast errors closer to those from the infeasible criterion. These patterns do not change when we have more observations and $T=200,300$. The conclusions are similar for $N=200$ in \Cref{tab2} and $N=300$ in \Cref{tab3}.




\section{Conclusion}
\label{Conclusion}
This paper explores rolling-window selection for forecasting with linear regression models with unobserved factors. It focuses on the situation where the latent factors are estimated using the principal component method, an unsupervised machine learning technique. The proposed approach is optimal in the sense that it minimizes the mean squared forecast error conditional on past information. The paper uses the squared distance between the fitted value and the unknown conditional mean at the end of our sample, which is the MSFE up to the constant $h$-step-ahead error, as a loss function, and
has two key theoretical results. First, we show that the window size that minimizes the infeasible loss function depends on factor-estimation uncertainty.
Second, the estimated rolling window based on the proposed criterion minimizes in probability the feasible version of our loss function. The new approach is theoretically valid, easy to implement, and yields excellent numerical results.





\clearpage